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 71 iterations with status: optimal.
ADMM terminated after 23 iterations with status: optimal.
ADMM terminated after 34 iterations with status: optimal.
ADMM terminated after 66 iterations with status: optimal.
ADMM terminated after 65 iterations with status: optimal.
ADMM terminated after 65 iterations with status: optimal.
ADMM terminated after 70 iterations with status: optimal.
ADMM terminated after 69 iterations with status: optimal.
ADMM terminated after 35 iterations with status: optimal.
ADMM terminated after 84 iterations with status: optimal.
ADMM terminated after 32 iterations with status: optimal.
ADMM terminated after 32 iterations with status: optimal.
ADMM terminated after 39 iterations with status: optimal.
ADMM terminated after 38 iterations with status: optimal.
ADMM terminated after 51 iterations with status: optimal.
ADMM terminated after 38 iterations with status: optimal.
ADMM terminated after 49 iterations with status: optimal.
ADMM terminated after 53 iterations with status: optimal.
ADMM terminated after 64 iterations with status: optimal.
ADMM terminated after 65 iterations with status: optimal.
ADMM terminated after 100 iterations with status: optimal.
ADMM terminated after 68 iterations with status: optimal.
ADMM terminated after 90 iterations with status: optimal.
ADMM terminated after 111 iterations with status: optimal.
ADMM terminated after 130 iterations with status: optimal.
ADMM terminated after 140 iterations with status: optimal.
ADMM terminated after 136 iterations with status: optimal.
ADMM terminated after 147 iterations with status: optimal.
ADMM terminated after 180 iterations with status: optimal.
ADMM terminated after 283 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.7 11.06 2.528e-06 2.526e-06
1 7.265 4.367 2.527e-06 2.527e-06
2 4.743 3.073 2.527e-06 2.527e-06
3 3.183 2.048 2.527e-06 2.528e-06
4 2.225 1.288 2.527e-06 2.528e-06
5 1.616 0.8068 2.527e-06 2.528e-06
6 1.217 0.5227 2.527e-06 2.528e-06
7 0.9448 0.3482 2.527e-06 2.528e-06
8 0.7547 0.2442 2.527e-06 2.528e-06
9 0.6187 0.172 2.527e-06 2.528e-06
10 0.5182 0.1258 2.527e-06 2.528e-06
11 0.4411 0.09441 2.527e-06 2.528e-06
12 0.3806 0.07319 2.527e-06 2.528e-06
13 0.3265 0.05804 2.527e-06 2.528e-06
14 0.2808 0.05274 2.527e-06 2.528e-06
15 0.2442 0.04918 2.527e-06 2.528e-06
16 0.2145 0.03724 2.527e-06 2.528e-06
17 0.1891 0.02967 2.527e-06 2.528e-06
18 0.1659 0.02732 2.527e-06 2.528e-06
19 0.1458 0.02327 2.527e-06 2.528e-06
20 0.1292 0.02192 2.527e-06 2.528e-06
21 0.1146 0.01734 2.527e-06 2.528e-06
22 0.1018 0.0152 2.527e-06 2.528e-06
23 0.09104 0.01335 2.527e-06 2.528e-06
24 0.08164 0.01077 2.527e-06 2.528e-06
25 0.07358 0.009678 2.527e-06 2.528e-06
26 0.06654 0.007684 2.527e-06 2.528e-06
27 0.06029 0.006625 2.527e-06 2.528e-06
28 0.05471 0.005856 2.527e-06 2.528e-06
29 0.04963 0.005218 2.527e-06 2.528e-06
30 0.04512 0.004819 2.527e-06 2.528e-06
31 0.04109 0.004372 2.527e-06 2.528e-06
32 0.03748 0.003793 2.527e-06 2.528e-06
33 0.03424 0.003387 2.527e-06 2.528e-06
34 0.02884 0.01101 2.527e-06 2.528e-06
35 0.02437 0.009212 2.527e-06 2.528e-06
36 0.02058 0.00767 2.527e-06 2.528e-06
37 0.01739 0.006671 2.527e-06 2.528e-06
38 0.01463 0.005894 2.527e-06 2.528e-06
39 0.01237 0.005006 2.527e-06 2.528e-06
40 0.01052 0.004328 2.527e-06 2.528e-06
41 0.008984 0.003291 2.527e-06 2.528e-06
42 0.007694 0.002681 2.527e-06 2.528e-06
43 0.006603 0.002248 2.527e-06 2.528e-06
44 0.005678 0.001902 2.527e-06 2.528e-06
45 0.004892 0.001617 2.527e-06 2.528e-06
46 0.004223 0.001378 2.527e-06 2.528e-06
47 0.003651 0.001177 2.527e-06 2.528e-06
48 0.003162 0.001007 2.527e-06 2.528e-06
49 0.002744 0.0008635 2.527e-06 2.528e-06
50 0.002384 0.0007414 2.527e-06 2.528e-06
51 0.002075 0.0006376 2.527e-06 2.528e-06
52 0.001809 0.0005491 2.527e-06 2.528e-06
53 0.00158 0.0004737 2.527e-06 2.528e-06
54 0.001382 0.0004092 2.527e-06 2.528e-06
55 0.00121 0.000354 2.527e-06 2.528e-06
56 0.001062 0.0003067 2.527e-06 2.528e-06
57 0.0009332 0.0002661 2.527e-06 2.528e-06
58 0.0008213 0.0002313 2.527e-06 2.528e-06
59 0.0007239 0.0002012 2.527e-06 2.528e-06
60 0.0006391 0.0001754 2.527e-06 2.528e-06
61 0.000565 0.000153 2.527e-06 2.528e-06
62 0.0005003 0.0001337 2.527e-06 2.528e-06
63 0.0004436 0.0001171 2.527e-06 2.528e-06
64 0.0003939 0.0001026 2.527e-06 2.528e-06
65 0.0003503 9.006e-05 2.527e-06 2.528e-06
66 0.0003119 7.916e-05 2.527e-06 2.528e-06
67 0.0002781 6.969e-05 2.527e-06 2.528e-06
68 0.0002483 6.143e-05 2.527e-06 2.528e-06
69 0.0002219 5.424e-05 2.527e-06 2.528e-06
70 0.0001986 4.795e-05 2.527e-06 2.528e-06
71 0.0001779 4.245e-05 2.527e-06 2.528e-06
72 0.0001596 3.764e-05 2.527e-06 2.528e-06
73 0.0001433 3.342e-05 2.527e-06 2.528e-06
74 0.0001288 2.971e-05 2.527e-06 2.528e-06
75 0.0001158 2.645e-05 2.527e-06 2.528e-06
76 0.0001043 2.357e-05 2.527e-06 2.528e-06
77 9.398e-05 2.104e-05 2.527e-06 2.528e-06
78 8.475e-05 1.88e-05 2.527e-06 2.528e-06
79 7.649e-05 1.681e-05 2.527e-06 2.528e-06
80 6.908e-05 1.506e-05 2.527e-06 2.528e-06
81 6.242e-05 1.35e-05 2.527e-06 2.528e-06
82 5.645e-05 1.212e-05 2.527e-06 2.528e-06
83 5.107e-05 1.088e-05 2.527e-06 2.528e-06
84 4.623e-05 9.787e-06 2.527e-06 2.528e-06
85 4.187e-05 8.808e-06 2.527e-06 2.528e-06
86 3.794e-05 7.934e-06 2.527e-06 2.528e-06
87 3.44e-05 7.152e-06 2.527e-06 2.528e-06
88 3.119e-05 6.452e-06 2.527e-06 2.528e-06
89 2.83e-05 5.825e-06 2.527e-06 2.528e-06
90 2.568e-05 5.262e-06 2.527e-06 2.528e-06
91 2.332e-05 4.756e-06 2.527e-06 2.528e-06
92 2.117e-05 4.302e-06 2.527e-06 2.528e-06
93 1.923e-05 3.893e-06 2.527e-06 2.528e-06
94 1.748e-05 3.525e-06 2.527e-06 2.528e-06
95 1.588e-05 3.193e-06 2.527e-06 2.528e-06
96 1.444e-05 2.894e-06 2.527e-06 2.528e-06
97 1.313e-05 2.624e-06 2.527e-06 2.528e-06
98 1.194e-05 2.38e-06 2.527e-06 2.528e-06
99 1.086e-05 2.159e-06 2.527e-06 2.528e-06
100 9.877e-06 1.96e-06 2.527e-06 2.528e-06
101 8.986e-06 1.779e-06 2.527e-06 2.528e-06
102 8.177e-06 1.616e-06 2.527e-06 2.528e-06
103 7.442e-06 1.468e-06 2.527e-06 2.528e-06
104 6.773e-06 1.334e-06 2.527e-06 2.528e-06
105 6.166e-06 1.212e-06 2.527e-06 2.528e-06
106 5.613e-06 1.102e-06 2.527e-06 2.528e-06
107 5.111e-06 1.002e-06 2.527e-06 2.528e-06
108 4.654e-06 9.111e-07 2.527e-06 2.528e-06
109 4.238e-06 8.287e-07 2.527e-06 2.528e-06
110 3.86e-06 7.539e-07 2.527e-06 2.528e-06
111 3.515e-06 6.859e-07 2.527e-06 2.528e-06
112 3.202e-06 6.241e-07 2.527e-06 2.528e-06
113 2.917e-06 5.68e-07 2.527e-06 2.528e-06
114 2.657e-06 5.17e-07 2.527e-06 2.528e-06
115 2.421e-06 4.706e-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: (0 minutes 43.164 seconds)