.. 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 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) .. 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.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. .. 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:** (0 minutes 43.164 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 `_