.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "auto_examples/plot_powerlaw_ggl.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code. .. rst-class:: sphx-glr-example-title .. _sphx_glr_auto_examples_plot_powerlaw_ggl.py: 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. .. GENERATED FROM PYTHON SOURCE LINES 9-33 .. code-block:: Python 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) .. GENERATED FROM PYTHON SOURCE LINES 34-42 Parameter selection (GGL) ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ We do a grid search over :math:`\lambda_1` and :math:`\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. .. GENERATED FROM PYTHON SOURCE LINES 42-100 .. code-block:: Python 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)) .. rst-class:: sphx-glr-script-out .. code-block:: none 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) .. GENERATED FROM PYTHON SOURCE LINES 101-107 Solving group sparse problems with SGL ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ We now solve K independent SGL problems and find the best :math:`\lambda_1` parameter. .. GENERATED FROM PYTHON SOURCE LINES 107-129 .. code-block:: Python 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'] .. GENERATED FROM PYTHON SOURCE LINES 130-135 Solving ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ To demonstrate again how to call the ADMM solver, we solve to high accuracy again for the best values of :math:`\lambda_1` and :math:`\lambda_2`. .. GENERATED FROM PYTHON SOURCE LINES 135-140 .. code-block:: Python 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) .. rst-class:: sphx-glr-script-out .. code-block:: none ------------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. .. GENERATED FROM PYTHON SOURCE LINES 141-148 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. .. GENERATED FROM PYTHON SOURCE LINES 148-155 .. code-block:: Python 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) .. rst-class:: sphx-glr-horizontal * .. image-sg:: /auto_examples/images/sphx_glr_plot_powerlaw_ggl_001.png :alt: Discovery rate for different regularization strengths :srcset: /auto_examples/images/sphx_glr_plot_powerlaw_ggl_001.png :class: sphx-glr-multi-img * .. image-sg:: /auto_examples/images/sphx_glr_plot_powerlaw_ggl_002.png :alt: Discovery of differential edges :srcset: /auto_examples/images/sphx_glr_plot_powerlaw_ggl_002.png :class: sphx-glr-multi-img .. rst-class:: sphx-glr-timing **Total running time of the script:** (1 minutes 32.632 seconds) .. _sphx_glr_download_auto_examples_plot_powerlaw_ggl.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: plot_powerlaw_ggl.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: plot_powerlaw_ggl.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: plot_powerlaw_ggl.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_