fix issues in glasso reg path example

Author: Perceptronium

examples/plot_graphical_lasso.py

@@ -1,43 +1,33 @@
+# Authors: Can Pouliquen
+#          Mathurin Massias
+"""
+=======================================================================
+Regularization paths for the Graphical Lasso and its Adaptive variation
+=======================================================================
+Highlight the importance of using non-convex regularization for improved performance,
+solved using the reweighting strategy.
+"""
+
 import numpy as np
 from numpy.linalg import norm
 import matplotlib.pyplot as plt
 from sklearn.metrics import f1_score
-from sklearn.datasets import make_sparse_spd_matrix
-from sklearn.utils import check_random_state
 
+from skglm.utils.data import generate_GraphicalLasso_data
 from skglm.estimators import GraphicalLasso
 from skglm.estimators import AdaptiveGraphicalLasso
 
-# Data
+
 p = 100
 n = 1000
-rng = check_random_state(0)
-Theta_true = make_sparse_spd_matrix(
-    p,
-    alpha=0.9,
-    random_state=rng)
-
-Theta_true += 0.1*np.eye(p)
-Sigma_true = np.linalg.pinv(Theta_true, hermitian=True)
-X = rng.multivariate_normal(
-    mean=np.zeros(p),
-    cov=Sigma_true,
-    size=n,
-)
-
-S = np.cov(X, bias=True, rowvar=False)
-S_cpy = np.copy(S)
-np.fill_diagonal(S_cpy, 0.)
-alpha_max = np.max(np.abs(S_cpy))
-
+S, Theta_true, alpha_max = generate_GraphicalLasso_data(n, p)
 alphas = alpha_max*np.geomspace(1, 1e-4, num=10)
 
-
 penalties = [
     "L1",
-    "R-L1 (log)",
-    "R-L1 (L0.5)",
-    "R-L1 (MCP)",
+    "Log",
+    "L0.5",
+    "MCP",
 ]
 n_reweights = 5
 models_tol = 1e-4
@@ -67,9 +57,8 @@
 
 
 for i, (penalty, model) in enumerate(zip(penalties, models)):
-    print(penalty)
     for alpha_idx, alpha in enumerate(alphas):
-        print(f"======= alpha {alpha_idx+1}/{len(alphas)} =======")
+        print(f"======= {penalty} penalty, alpha {alpha_idx+1}/{len(alphas)} =======")
         model.alpha = alpha
         model.fit(S)
         Theta = model.precision_
@@ -78,66 +67,54 @@
 
         my_f1_score = f1_score(Theta.flatten() != 0.,
                                Theta_true.flatten() != 0.)
-        print(f"NMSE: {my_nmse:.3f}")
-        print(f"F1  : {my_f1_score:.3f}")
 
         my_glasso_nmses[penalty].append(my_nmse)
         my_glasso_f1_scores[penalty].append(my_f1_score)
 
 
 plt.close('all')
-fig, ax = plt.subplots(2, 1, sharex=True, figsize=(
-    [6.11, 3.91]), layout="constrained")
+fig, axarr = plt.subplots(2, 1, sharex=True, figsize=([6.11, 3.91]),
+                          layout="constrained")
 cmap = plt.get_cmap("tab10")
 for i, penalty in enumerate(penalties):
 
-    ax[0].semilogx(alphas/alpha_max,
-                   my_glasso_nmses[penalty],
-                   color=cmap(i),
-                   linewidth=2.,
-                   label=penalty)
-    min_nmse = np.argmin(my_glasso_nmses[penalty])
-    ax[0].vlines(
-        x=alphas[min_nmse] / alphas[0],
-        ymin=0,
-        ymax=np.min(my_glasso_nmses[penalty]),
-        linestyle='--',
-        color=cmap(i))
-    line0 = ax[0].plot(
-        [alphas[min_nmse] / alphas[0]],
-        0,
-        clip_on=False,
-        marker='X',
-        color=cmap(i),
-        markersize=12)
-
-    ax[1].semilogx(alphas/alpha_max,
-                   my_glasso_f1_scores[penalty],
-                   linewidth=2.,
-                   color=cmap(i))
-    max_f1 = np.argmax(my_glasso_f1_scores[penalty])
-    ax[1].vlines(
-        x=alphas[max_f1] / alphas[0],
-        ymin=0,
-        ymax=np.max(my_glasso_f1_scores[penalty]),
-        linestyle='--',
-        color=cmap(i))
-    line1 = ax[1].plot(
-        [alphas[max_f1] / alphas[0]],
-        0,
-        clip_on=False,
-        marker='X',
-        markersize=12,
-        color=cmap(i))
-
-
-ax[0].set_title(f"{p=},{n=}", fontsize=18)
-ax[0].set_ylabel("NMSE", fontsize=18)
-ax[1].set_ylabel("F1 score", fontsize=18)
-ax[1].set_xlabel(f"$\lambda / \lambda_\mathrm{{max}}$",  fontsize=18)
-
-ax[0].legend(fontsize=14)
-ax[0].grid(which='both', alpha=0.9)
-ax[1].grid(which='both', alpha=0.9)
-# plt.savefig(f"./non_convex_p{p}_n{n}.pdf")
+    for j, ax in enumerate(axarr):
+
+        if j == 0:
+            metric = my_glasso_nmses
+            best_idx = np.argmin(metric[penalty])
+            ystop = np.min(metric[penalty])
+        else:
+            metric = my_glasso_f1_scores
+            best_idx = np.argmax(metric[penalty])
+            ystop = np.max(metric[penalty])
+
+        ax.semilogx(alphas/alpha_max,
+                    metric[penalty],
+                    color=cmap(i),
+                    linewidth=2.,
+                    label=penalty)
+
+        ax.vlines(
+            x=alphas[best_idx] / alphas[0],
+            ymin=0,
+            ymax=ystop,
+            linestyle='--',
+            color=cmap(i))
+        line = ax.plot(
+            [alphas[best_idx] / alphas[0]],
+            0,
+            clip_on=False,
+            marker='X',
+            color=cmap(i),
+            markersize=12)
+
+        ax.grid(which='both', alpha=0.9)
+
+axarr[0].legend(fontsize=14)
+axarr[0].set_title(f"{p=},{n=}", fontsize=18)
+axarr[0].set_ylabel("NMSE", fontsize=18)
+axarr[1].set_ylabel("F1 score", fontsize=18)
+axarr[1].set_xlabel(r"$\lambda / \lambda_\mathrm{{max}}$",  fontsize=18)
+
 plt.show(block=False)