Note
Go to the end to download the full example code.
Group Graphical Lasso experiment
We investigate the recovery performance of Group Graphical Lasso on Powerlaw networks, compared to estimating the precision matrices independently with SGL.
We generate a precision matrix with block-wise powerlaw networks. In each instance, one of the blocks is randomly set to zero. Hence, the true underlying precision matrices are group sparse.
import numpy as np
from matplotlib import pyplot as plt
import seaborn as sns
from sklearn.covariance import GraphicalLasso
from gglasso.solver.admm_solver import ADMM_MGL
from gglasso.helper.data_generation import group_power_network, sample_covariance_matrix
from gglasso.helper.experiment_helper import lambda_parametrizer, discovery_rate, error
from gglasso.helper.utils import get_K_identity
from gglasso.helper.experiment_helper import draw_group_heatmap, plot_fpr_tpr, plot_diff_fpr_tpr, surface_plot
from gglasso.helper.model_selection import aic, ebic
p = 100
K = 5
N = 80
M = 10
reg = 'GGL'
Sigma, Theta = group_power_network(p, K, M, scale = False, seed = 2345)
S, sample = sample_covariance_matrix(Sigma, N)
Parameter selection (GGL)
We do a grid search over \(\lambda_1\) and \(\lambda_2\) values. On each grid point we evaluate True/False Discovery Rate (TPR/FPR), True/False Discovery of Differential edges and AIC and eBIC.
Note: the package contains functions for doing this grid search, but here we also want to evaluate True and False positive rates on each grid points.
l2 = np.logspace(start = -0.5, stop = -2.5, num = 10, base = 10)
W2 = np.linspace(0.2, 0.5, 3)
l2grid, w2grid = np.meshgrid(l2, W2)
L1 = lambda_parametrizer(l2grid, w2grid)
L2 = l2grid.copy()
grid1 = L1.shape[0]; grid2 = L2.shape[1]
ERR = np.zeros((grid1, grid2))
FPR = np.zeros((grid1, grid2))
TPR = np.zeros((grid1, grid2))
DFPR = np.zeros((grid1, grid2))
DTPR = np.zeros((grid1, grid2))
AIC = np.zeros((grid1, grid2))
BIC = np.zeros((grid1, grid2))
Omega_0 = get_K_identity(K,p)
Theta_0 = get_K_identity(K,p)
X_0 = np.zeros((K,p,p))
for g2 in np.arange(grid2):
for g1 in np.arange(grid1):
lambda1 = L1[g1,g2]
lambda2 = L2[g1,g2]
sol, info = ADMM_MGL(S, lambda1, lambda2, reg , Omega_0, Theta_0 = Theta_0, X_0 = X_0, tol = 1e-8, rtol = 1e-8, verbose = False, measure = False)
Theta_sol = sol['Theta']
Omega_sol = sol['Omega']
X_sol = sol['X']
# warm start
Omega_0 = Omega_sol.copy()
Theta_0 = Theta_sol.copy()
X_0 = X_sol.copy()
dr = discovery_rate(Theta_sol, Theta)
TPR[g1,g2] = dr['TPR']
FPR[g1,g2] = dr['FPR']
DTPR[g1,g2] = dr['TPR_DIFF']
DFPR[g1,g2] = dr['FPR_DIFF']
ERR[g1,g2] = error(Theta_sol, Theta)
AIC[g1,g2] = aic(S, Theta_sol, N)
BIC[g1,g2] = ebic(S, Theta_sol, N, gamma = 0.1)
# get optimal lambda
ix= np.unravel_index(np.nanargmin(BIC), BIC.shape)
ix2= np.unravel_index(np.nanargmin(AIC), AIC.shape)
l1opt = L1[ix]
l2opt = L2[ix]
print("Optimal lambda values: (l1,l2) = ", (l1opt,l2opt))
ADMM terminated after 72 iterations with status: optimal.
ADMM terminated after 21 iterations with status: optimal.
ADMM terminated after 33 iterations with status: optimal.
ADMM terminated after 27 iterations with status: optimal.
ADMM terminated after 62 iterations with status: optimal.
ADMM terminated after 61 iterations with status: optimal.
ADMM terminated after 64 iterations with status: optimal.
ADMM terminated after 68 iterations with status: optimal.
ADMM terminated after 50 iterations with status: optimal.
ADMM terminated after 42 iterations with status: optimal.
ADMM terminated after 45 iterations with status: optimal.
ADMM terminated after 34 iterations with status: optimal.
ADMM terminated after 42 iterations with status: optimal.
ADMM terminated after 40 iterations with status: optimal.
ADMM terminated after 53 iterations with status: optimal.
ADMM terminated after 41 iterations with status: optimal.
ADMM terminated after 43 iterations with status: optimal.
ADMM terminated after 53 iterations with status: optimal.
ADMM terminated after 66 iterations with status: optimal.
ADMM terminated after 66 iterations with status: optimal.
ADMM terminated after 101 iterations with status: optimal.
ADMM terminated after 68 iterations with status: optimal.
ADMM terminated after 79 iterations with status: optimal.
ADMM terminated after 112 iterations with status: optimal.
ADMM terminated after 131 iterations with status: optimal.
ADMM terminated after 140 iterations with status: optimal.
ADMM terminated after 141 iterations with status: optimal.
ADMM terminated after 147 iterations with status: optimal.
ADMM terminated after 179 iterations with status: optimal.
ADMM terminated after 282 iterations with status: optimal.
Optimal lambda values: (l1,l2) = (0.22360679774997896, 0.31622776601683794)
Solving group sparse problems with SGL
We now solve K independent SGL problems and find the best \(\lambda_1\) parameter.
ALPHA = np.logspace(start = 0, stop = -1.5, num = 15, base = 10)
FPR_GL = np.zeros(len(ALPHA))
TPR_GL = np.zeros(len(ALPHA))
DFPR_GL = np.zeros(len(ALPHA))
DTPR_GL = np.zeros(len(ALPHA))
for a in np.arange(len(ALPHA)):
singleGL = GraphicalLasso(alpha = ALPHA[a], tol = 1e-4, max_iter = 50, verbose = False)
singleGL_sol = np.zeros((K,p,p))
for k in np.arange(K):
model = singleGL.fit(sample[k,:,:].T)
singleGL_sol[k,:,:] = model.precision_
dr = discovery_rate(singleGL_sol, Theta)
TPR_GL[a] = dr['TPR']
FPR_GL[a] = dr['FPR']
DTPR_GL[a] = dr['TPR_DIFF']
DFPR_GL[a] = dr['FPR_DIFF']
Solving
To demonstrate again how to call the ADMM solver, we solve to high accuracy again for the best values of \(\lambda_1\) and \(\lambda_2\).
Omega_0 = get_K_identity(K,p)
solA, infoA = ADMM_MGL(S, l1opt, l2opt, reg , Omega_0, tol = 1e-10, rtol = 1e-10, verbose = True, measure = True)
------------ADMM Algorithm for Multiple Graphical Lasso----------------
iter r_t s_t eps_pri eps_dual
0 10.67 11.06 2.528e-06 2.526e-06
1 7.213 4.359 2.528e-06 2.527e-06
2 4.717 3.05 2.527e-06 2.527e-06
3 3.163 2.026 2.527e-06 2.527e-06
4 2.199 1.274 2.527e-06 2.528e-06
5 1.587 0.8095 2.527e-06 2.528e-06
6 1.186 0.5182 2.527e-06 2.528e-06
7 0.9111 0.3452 2.527e-06 2.528e-06
8 0.7171 0.2376 2.527e-06 2.528e-06
9 0.576 0.169 2.527e-06 2.528e-06
10 0.4721 0.1264 2.527e-06 2.528e-06
11 0.391 0.09299 2.527e-06 2.528e-06
12 0.3283 0.07373 2.527e-06 2.528e-06
13 0.2784 0.0611 2.527e-06 2.528e-06
14 0.2384 0.04548 2.527e-06 2.528e-06
15 0.2037 0.03885 2.527e-06 2.528e-06
16 0.1745 0.03365 2.527e-06 2.528e-06
17 0.1515 0.02953 2.527e-06 2.528e-06
18 0.133 0.02199 2.527e-06 2.528e-06
19 0.1162 0.01751 2.527e-06 2.528e-06
20 0.1025 0.0163 2.527e-06 2.528e-06
21 0.0911 0.01324 2.527e-06 2.528e-06
22 0.08135 0.01065 2.527e-06 2.528e-06
23 0.0729 0.009075 2.527e-06 2.528e-06
24 0.06552 0.007872 2.527e-06 2.528e-06
25 0.05887 0.00688 2.527e-06 2.528e-06
26 0.05279 0.006114 2.527e-06 2.528e-06
27 0.04762 0.006792 2.527e-06 2.528e-06
28 0.0432 0.004981 2.527e-06 2.528e-06
29 0.0386 0.004196 2.527e-06 2.528e-06
30 0.03482 0.004869 2.527e-06 2.528e-06
31 0.03162 0.004036 2.527e-06 2.528e-06
32 0.02883 0.003097 2.527e-06 2.528e-06
33 0.02636 0.002668 2.527e-06 2.528e-06
34 0.02414 0.002366 2.527e-06 2.528e-06
35 0.02048 0.007638 2.527e-06 2.528e-06
36 0.01746 0.00634 2.527e-06 2.528e-06
37 0.0149 0.005258 2.527e-06 2.528e-06
38 0.01281 0.004802 2.527e-06 2.528e-06
39 0.01107 0.003738 2.527e-06 2.528e-06
40 0.00961 0.003076 2.527e-06 2.528e-06
41 0.008372 0.00259 2.527e-06 2.528e-06
42 0.007318 0.0022 2.527e-06 2.528e-06
43 0.006415 0.001879 2.527e-06 2.528e-06
44 0.005638 0.001613 2.527e-06 2.528e-06
45 0.004967 0.00139 2.527e-06 2.528e-06
46 0.004385 0.001203 2.527e-06 2.528e-06
47 0.003879 0.001044 2.527e-06 2.528e-06
48 0.003437 0.0009096 2.527e-06 2.528e-06
49 0.003049 0.0007948 2.527e-06 2.528e-06
50 0.002709 0.0006963 2.527e-06 2.528e-06
51 0.00241 0.0006115 2.527e-06 2.528e-06
52 0.002146 0.0005383 2.527e-06 2.528e-06
53 0.001913 0.0004748 2.527e-06 2.528e-06
54 0.001707 0.0004196 2.527e-06 2.528e-06
55 0.001524 0.0003714 2.527e-06 2.528e-06
56 0.001362 0.0003293 2.527e-06 2.528e-06
57 0.001217 0.0002923 2.527e-06 2.528e-06
58 0.001089 0.0002598 2.527e-06 2.528e-06
59 0.0009747 0.0002312 2.527e-06 2.528e-06
60 0.0008727 0.0002059 2.527e-06 2.528e-06
61 0.0007817 0.0001836 2.527e-06 2.528e-06
62 0.0007005 0.0001638 2.527e-06 2.528e-06
63 0.0006279 0.0001462 2.527e-06 2.528e-06
64 0.0005631 0.0001306 2.527e-06 2.528e-06
65 0.000505 0.0001168 2.527e-06 2.528e-06
66 0.0004531 0.0001045 2.527e-06 2.528e-06
67 0.0004066 9.348e-05 2.527e-06 2.528e-06
68 0.0003649 8.369e-05 2.527e-06 2.528e-06
69 0.0003276 7.496e-05 2.527e-06 2.528e-06
70 0.0002942 6.716e-05 2.527e-06 2.528e-06
71 0.0002642 6.02e-05 2.527e-06 2.528e-06
72 0.0002373 5.397e-05 2.527e-06 2.528e-06
73 0.0002132 4.841e-05 2.527e-06 2.528e-06
74 0.0001915 4.342e-05 2.527e-06 2.528e-06
75 0.0001721 3.896e-05 2.527e-06 2.528e-06
76 0.0001547 3.497e-05 2.527e-06 2.528e-06
77 0.000139 3.139e-05 2.527e-06 2.528e-06
78 0.000125 2.818e-05 2.527e-06 2.528e-06
79 0.0001123 2.531e-05 2.527e-06 2.528e-06
80 0.000101 2.273e-05 2.527e-06 2.528e-06
81 9.079e-05 2.042e-05 2.527e-06 2.528e-06
82 8.164e-05 1.834e-05 2.527e-06 2.528e-06
83 7.341e-05 1.648e-05 2.527e-06 2.528e-06
84 6.602e-05 1.481e-05 2.527e-06 2.528e-06
85 5.937e-05 1.331e-05 2.527e-06 2.528e-06
86 5.34e-05 1.196e-05 2.527e-06 2.528e-06
87 4.803e-05 1.075e-05 2.527e-06 2.528e-06
88 4.32e-05 9.666e-06 2.527e-06 2.528e-06
89 3.886e-05 8.69e-06 2.527e-06 2.528e-06
90 3.496e-05 7.813e-06 2.527e-06 2.528e-06
91 3.145e-05 7.025e-06 2.527e-06 2.528e-06
92 2.829e-05 6.317e-06 2.527e-06 2.528e-06
93 2.546e-05 5.681e-06 2.527e-06 2.528e-06
94 2.29e-05 5.109e-06 2.527e-06 2.528e-06
95 2.061e-05 4.595e-06 2.527e-06 2.528e-06
96 1.854e-05 4.133e-06 2.527e-06 2.528e-06
97 1.668e-05 3.718e-06 2.527e-06 2.528e-06
98 1.501e-05 3.344e-06 2.527e-06 2.528e-06
99 1.351e-05 3.009e-06 2.527e-06 2.528e-06
100 1.216e-05 2.707e-06 2.527e-06 2.528e-06
101 1.094e-05 2.435e-06 2.527e-06 2.528e-06
102 9.846e-06 2.191e-06 2.527e-06 2.528e-06
103 8.861e-06 1.971e-06 2.527e-06 2.528e-06
104 7.975e-06 1.773e-06 2.527e-06 2.528e-06
105 7.178e-06 1.596e-06 2.527e-06 2.528e-06
106 6.46e-06 1.436e-06 2.527e-06 2.528e-06
107 5.814e-06 1.292e-06 2.527e-06 2.528e-06
108 5.233e-06 1.163e-06 2.527e-06 2.528e-06
109 4.71e-06 1.046e-06 2.527e-06 2.528e-06
110 4.24e-06 9.416e-07 2.527e-06 2.528e-06
111 3.816e-06 8.473e-07 2.527e-06 2.528e-06
112 3.435e-06 7.626e-07 2.527e-06 2.528e-06
113 3.092e-06 6.863e-07 2.527e-06 2.528e-06
114 2.783e-06 6.177e-07 2.527e-06 2.528e-06
115 2.505e-06 5.559e-07 2.527e-06 2.528e-06
ADMM terminated after 116 iterations with status: optimal.
Plotting: TPR, FPR, differential edges
We plot FPR and TPR for all grid points for GGL and SGL. The circle-shape marker marks the point which would have been selected by eBIC. The diamond-shaped marker marks the point selected by AIC.
Differential edges are edges which are present in at least one but not all of the K precision matrices.
fig,ax = plot_fpr_tpr(FPR, TPR, ix, ix2, FPR_GL, TPR_GL, W2)
ax.set_xlim(-0.001, 0.2)
ax.set_ylim(0.1,1)
fig,ax = plot_diff_fpr_tpr(DFPR, DTPR, ix, ix2, DFPR_GL, DTPR_GL, W2)
Total running time of the script: (1 minutes 32.632 seconds)