graphical_lasso.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "adult-source",
   "metadata": {},
   "source": [
    "# Graphical Lasso (GLasso)\n",
    "Reference: Friedman, Jerome, Trevor Hastie, and Robert Tibshirani. \"Sparse inverse covariance estimation with the graphical lasso.\" Biostatistics 9.3 (2008): 432-441."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fiscal-allocation",
   "metadata": {},
   "source": [
    "The objective function is $$max_\\Theta lndet\\Theta - tr(S\\Theta)-\\rho\\|\\Theta\\|_1,$$\n",
    "\n",
    "where $\\Theta$ is the precision (inverse covariance) matrix, $S=\\frac1nX^TX$ is empirical sample covariance, $\\|\\Theta\\|_1$ is the sum of absolute values of elements in $\\Theta$.\n",
    "\n",
    "The subgradient equation for objective function is $$\\Sigma-S-\\rho\\Gamma=0,$$\n",
    "\n",
    "where $\\Sigma$ is the covariance matrix, $\\Gamma_{ij}=sign(\\Theta_{ij})$ if $\\Theta_{ij}\\neq0$, $\\Gamma_{ij}\\in[-1,1]$ if $\\Theta_{ij}=0$.\n",
    "\n",
    "For diagonal elements, $\\Theta_{ii}>0$ then $\\Gamma_{ii}=1$, so $\\Sigma_{ii}=S_{ii}+\\rho$.\n",
    "\n",
    "Consider 1 versus (p-1) partition over $p$ variables, let\n",
    "$$\n",
    "\\left(\\begin{array}{cc} \\Sigma_{aa} & \\Sigma_{ab} \\\\ \\Sigma_{ba} & \\Sigma_{bb} \\end{array}\\right)\n",
    "\\left(\\begin{array}{cc} \\Theta_{aa} & \\Theta_{ab} \\\\ \\Theta_{ba} & \\Theta_{bb} \\end{array}\\right) = \n",
    "\\left(\\begin{array}{cc} I_{p-1} & \\textbf{1} \\\\ \\textbf{1}^T & 1 \\end{array}\\right)\n",
    "$$\n",
    "For the upper right block in the subgradient equation, $\\Sigma_{ab}-S_{ab}-\\rho\\Gamma_{ab}=0$, whose solution is equivalent with another optimization problem $$min_{\\beta} \\frac12\\|\\Sigma_{aa}^{1/2} \\beta - \\Sigma_{aa}^{-1} S_{ab} \\|^2+\\rho\\|\\beta\\|_1,$$\n",
    "whose subgradient equation is $\\Sigma_{aa}\\beta - S_{ab}+\\rho sign(\\beta)=0$, where $\\beta=\\Sigma_{aa}^{-1}\\Sigma_{ab}, sign(\\beta)=-\\Gamma_{ab}$.\n",
    "\n",
    "###### Note the connection and difference between column-wise solving 1 versus (p-1) conditional regression by lasso and the glasso over whole precision matrix is that $p$ seperate conditional regression is lasso($S_{aa},S_{ab},\\rho$), while glasso solves $p$ coupled lasso problems with the same $Sigma$ shared as lasso($\\Sigma_{ab}, S_{ab}, \\rho$)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fundamental-depression",
   "metadata": {},
   "source": [
    "  ###### Graphical Lasso algorithm\n",
    "Let $W=\\hat\\Sigma$ be the estimated covariance matrix,  is as follows\n",
    "1. $W=S+\\rho I_p$, then the diagonal remains unchanged\n",
    "2. Repeat $j=1,2,...,p,1,2,...,p,...$, for each column $j$ in $W$, solve $\\hat\\beta=$lasso($\\Sigma_{ab}, S_{ab}, \\rho$), update $W_{ab}=W_{aa}\\hat\\beta$, until convergence of $W$\n",
    "\n",
    "Details of step 2 in algorithm:\n",
    "1) circle through all the variables multiple times, for each $j$-th column in the original $W$ is switched to the last column then it is updated as the upper right block $W_{ab}$ is updated  \n",
    "2) lasso($\\Sigma_{ab}, S_{ab}, \\rho$) is solved by coordinate descent, where each entry in $\\beta$ is updated iteratively by $$\\hat\\beta_i \\leftarrow ST\\left( (S_{ab})_{i} - \\Sigma_{k\\neq i} (W_{aa})_{ki}\\hat\\beta_k, \\rho \\right) / { (W_{aa})_{ii} },$$ where ST is the soft-threshold operator $ST(x, \\rho)=sign(x)(|x|-t)_+$  \n",
    "3) Convergence of $W$ is verified by its average absolute change less than the multiplication between a fixed threshold $t=0.001$ and the average of absolute value of off-diagonal elements in empirical covariance matrix $S$  \n",
    "4) After convergence of $W$, the estimated precision matrix can be obtained by $\\hat\\beta$ s, according to $\\hat\\Theta_{bb}=1/(W_{bb}-W_{ab}^T\\hat\\beta)$ and $\\hat\\Theta_{ab}=-\\hat\\beta\\Theta_{bb}$.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "different-aruba",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table class=\"dataframe\">\n",
       "<caption>A matrix: 10 × 10 of type dbl</caption>\n",
       "<tbody>\n",
       "\t<tr><td> 1.27498518</td><td>-0.360013241</td><td>-0.13899649</td><td>-0.04053069</td><td>-0.108175498</td><td>-0.22897547</td><td> 0.43052975</td><td> 0.01709085</td><td>-0.39463630</td><td> 0.02215637</td></tr>\n",
       "\t<tr><td>-0.36001324</td><td> 3.134190094</td><td> 1.01737618</td><td>-0.77021533</td><td>-0.009104576</td><td> 0.37902948</td><td>-0.44378911</td><td> 0.56362750</td><td> 0.11966498</td><td>-0.46572278</td></tr>\n",
       "\t<tr><td>-0.13899649</td><td> 1.017376182</td><td> 1.38116548</td><td>-0.96414716</td><td>-0.030133155</td><td> 0.05454111</td><td>-0.25264663</td><td>-0.01056241</td><td> 0.05395003</td><td>-0.35386972</td></tr>\n",
       "\t<tr><td>-0.04053069</td><td>-0.770215326</td><td>-0.96414716</td><td> 3.52498851</td><td> 0.656934400</td><td>-0.39441460</td><td> 0.11595627</td><td>-0.58714774</td><td> 0.03967201</td><td>-0.27566292</td></tr>\n",
       "\t<tr><td>-0.10817550</td><td>-0.009104576</td><td>-0.03013316</td><td> 0.65693440</td><td> 1.878431591</td><td> 0.14829337</td><td>-0.20895989</td><td>-0.30118765</td><td> 0.03356446</td><td>-0.09224890</td></tr>\n",
       "\t<tr><td>-0.22897547</td><td> 0.379029484</td><td> 0.05454111</td><td>-0.39441460</td><td> 0.148293373</td><td> 1.30772526</td><td>-0.04191247</td><td> 0.11365076</td><td>-0.06444851</td><td>-0.14225305</td></tr>\n",
       "\t<tr><td> 0.43052975</td><td>-0.443789113</td><td>-0.25264663</td><td> 0.11595627</td><td>-0.208959889</td><td>-0.04191247</td><td> 1.32480422</td><td>-0.07923855</td><td>-0.22972662</td><td>-0.20156463</td></tr>\n",
       "\t<tr><td> 0.01709085</td><td> 0.563627505</td><td>-0.01056241</td><td>-0.58714774</td><td>-0.301187650</td><td> 0.11365076</td><td>-0.07923855</td><td> 1.38270476</td><td>-0.05168962</td><td> 0.31653789</td></tr>\n",
       "\t<tr><td>-0.39463630</td><td> 0.119664977</td><td> 0.05395003</td><td> 0.03967201</td><td> 0.033564464</td><td>-0.06444851</td><td>-0.22972662</td><td>-0.05168962</td><td> 0.93065700</td><td>-0.04880206</td></tr>\n",
       "\t<tr><td> 0.02215637</td><td>-0.465722778</td><td>-0.35386972</td><td>-0.27566292</td><td>-0.092248902</td><td>-0.14225305</td><td>-0.20156463</td><td> 0.31653789</td><td>-0.04880206</td><td> 1.77034700</td></tr>\n",
       "</tbody>\n",
       "</table>\n"
      ],
      "text/latex": [
       "A matrix: 10 × 10 of type dbl\n",
       "\\begin{tabular}{llllllllll}\n",
       "\t  1.27498518 & -0.360013241 & -0.13899649 & -0.04053069 & -0.108175498 & -0.22897547 &  0.43052975 &  0.01709085 & -0.39463630 &  0.02215637\\\\\n",
       "\t -0.36001324 &  3.134190094 &  1.01737618 & -0.77021533 & -0.009104576 &  0.37902948 & -0.44378911 &  0.56362750 &  0.11966498 & -0.46572278\\\\\n",
       "\t -0.13899649 &  1.017376182 &  1.38116548 & -0.96414716 & -0.030133155 &  0.05454111 & -0.25264663 & -0.01056241 &  0.05395003 & -0.35386972\\\\\n",
       "\t -0.04053069 & -0.770215326 & -0.96414716 &  3.52498851 &  0.656934400 & -0.39441460 &  0.11595627 & -0.58714774 &  0.03967201 & -0.27566292\\\\\n",
       "\t -0.10817550 & -0.009104576 & -0.03013316 &  0.65693440 &  1.878431591 &  0.14829337 & -0.20895989 & -0.30118765 &  0.03356446 & -0.09224890\\\\\n",
       "\t -0.22897547 &  0.379029484 &  0.05454111 & -0.39441460 &  0.148293373 &  1.30772526 & -0.04191247 &  0.11365076 & -0.06444851 & -0.14225305\\\\\n",
       "\t  0.43052975 & -0.443789113 & -0.25264663 &  0.11595627 & -0.208959889 & -0.04191247 &  1.32480422 & -0.07923855 & -0.22972662 & -0.20156463\\\\\n",
       "\t  0.01709085 &  0.563627505 & -0.01056241 & -0.58714774 & -0.301187650 &  0.11365076 & -0.07923855 &  1.38270476 & -0.05168962 &  0.31653789\\\\\n",
       "\t -0.39463630 &  0.119664977 &  0.05395003 &  0.03967201 &  0.033564464 & -0.06444851 & -0.22972662 & -0.05168962 &  0.93065700 & -0.04880206\\\\\n",
       "\t  0.02215637 & -0.465722778 & -0.35386972 & -0.27566292 & -0.092248902 & -0.14225305 & -0.20156463 &  0.31653789 & -0.04880206 &  1.77034700\\\\\n",
       "\\end{tabular}\n"
      ],
      "text/markdown": [
       "\n",
       "A matrix: 10 × 10 of type dbl\n",
       "\n",
       "|  1.27498518 | -0.360013241 | -0.13899649 | -0.04053069 | -0.108175498 | -0.22897547 |  0.43052975 |  0.01709085 | -0.39463630 |  0.02215637 |\n",
       "| -0.36001324 |  3.134190094 |  1.01737618 | -0.77021533 | -0.009104576 |  0.37902948 | -0.44378911 |  0.56362750 |  0.11966498 | -0.46572278 |\n",
       "| -0.13899649 |  1.017376182 |  1.38116548 | -0.96414716 | -0.030133155 |  0.05454111 | -0.25264663 | -0.01056241 |  0.05395003 | -0.35386972 |\n",
       "| -0.04053069 | -0.770215326 | -0.96414716 |  3.52498851 |  0.656934400 | -0.39441460 |  0.11595627 | -0.58714774 |  0.03967201 | -0.27566292 |\n",
       "| -0.10817550 | -0.009104576 | -0.03013316 |  0.65693440 |  1.878431591 |  0.14829337 | -0.20895989 | -0.30118765 |  0.03356446 | -0.09224890 |\n",
       "| -0.22897547 |  0.379029484 |  0.05454111 | -0.39441460 |  0.148293373 |  1.30772526 | -0.04191247 |  0.11365076 | -0.06444851 | -0.14225305 |\n",
       "|  0.43052975 | -0.443789113 | -0.25264663 |  0.11595627 | -0.208959889 | -0.04191247 |  1.32480422 | -0.07923855 | -0.22972662 | -0.20156463 |\n",
       "|  0.01709085 |  0.563627505 | -0.01056241 | -0.58714774 | -0.301187650 |  0.11365076 | -0.07923855 |  1.38270476 | -0.05168962 |  0.31653789 |\n",
       "| -0.39463630 |  0.119664977 |  0.05395003 |  0.03967201 |  0.033564464 | -0.06444851 | -0.22972662 | -0.05168962 |  0.93065700 | -0.04880206 |\n",
       "|  0.02215637 | -0.465722778 | -0.35386972 | -0.27566292 | -0.092248902 | -0.14225305 | -0.20156463 |  0.31653789 | -0.04880206 |  1.77034700 |\n",
       "\n"
      ],
      "text/plain": [
       "      [,1]        [,2]         [,3]        [,4]        [,5]         [,6]       \n",
       " [1,]  1.27498518 -0.360013241 -0.13899649 -0.04053069 -0.108175498 -0.22897547\n",
       " [2,] -0.36001324  3.134190094  1.01737618 -0.77021533 -0.009104576  0.37902948\n",
       " [3,] -0.13899649  1.017376182  1.38116548 -0.96414716 -0.030133155  0.05454111\n",
       " [4,] -0.04053069 -0.770215326 -0.96414716  3.52498851  0.656934400 -0.39441460\n",
       " [5,] -0.10817550 -0.009104576 -0.03013316  0.65693440  1.878431591  0.14829337\n",
       " [6,] -0.22897547  0.379029484  0.05454111 -0.39441460  0.148293373  1.30772526\n",
       " [7,]  0.43052975 -0.443789113 -0.25264663  0.11595627 -0.208959889 -0.04191247\n",
       " [8,]  0.01709085  0.563627505 -0.01056241 -0.58714774 -0.301187650  0.11365076\n",
       " [9,] -0.39463630  0.119664977  0.05395003  0.03967201  0.033564464 -0.06444851\n",
       "[10,]  0.02215637 -0.465722778 -0.35386972 -0.27566292 -0.092248902 -0.14225305\n",
       "      [,7]        [,8]        [,9]        [,10]      \n",
       " [1,]  0.43052975  0.01709085 -0.39463630  0.02215637\n",
       " [2,] -0.44378911  0.56362750  0.11966498 -0.46572278\n",
       " [3,] -0.25264663 -0.01056241  0.05395003 -0.35386972\n",
       " [4,]  0.11595627 -0.58714774  0.03967201 -0.27566292\n",
       " [5,] -0.20895989 -0.30118765  0.03356446 -0.09224890\n",
       " [6,] -0.04191247  0.11365076 -0.06444851 -0.14225305\n",
       " [7,]  1.32480422 -0.07923855 -0.22972662 -0.20156463\n",
       " [8,] -0.07923855  1.38270476 -0.05168962  0.31653789\n",
       " [9,] -0.22972662 -0.05168962  0.93065700 -0.04880206\n",
       "[10,] -0.20156463  0.31653789 -0.04880206  1.77034700"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Set the ground truth of precision matrix (\\Theta)\n",
    "n_aux = 20\n",
    "p = 10\n",
    "set.seed(0)\n",
    "X_aux = matrix( rnorm(n_aux*p), ncol=p ) # Auxillary data for generating truth of covariance\n",
    "Sigma_true = cov(X_aux)\n",
    "# Sigma_true\n",
    "Theta_true = solve(Sigma_true)\n",
    "Theta_true"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "chemical-debate",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Generate data from true precision (true covariance)\n",
    "n = 9\n",
    "X = MASS::mvrnorm(n=n, mu=rep(0,p), Sigma=Sigma_true)\n",
    "# Settings of glasso algorithm\n",
    "rho = 0.1\n",
    "t = 1e-4 # a fixed threshold\n",
    "S = cov(X)*(n-1)/n # empirical"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "hindu-amplifier",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "iter 1 0.1000172 \n",
      "iter 2 0.03564175 \n",
      "iter 3 0.001287652 \n",
      "iter 4 5.914316e-05 \n",
      "iter 5 1.559448e-05 \n"
     ]
    }
   ],
   "source": [
    "ave_S = mean(S[lower.tri(S)])\n",
    "threshold = t * ave_S\n",
    "max_iter = 100\n",
    "\n",
    "# Step 1: diagonal of W remains unchanged\n",
    "W = S + rho*diag(p)\n",
    "# print(W)\n",
    "# B = matrix(0, (p-1), p) # store all betas\n",
    "\n",
    "# Step 2:\n",
    "iter = 0\n",
    "while(iter < max_iter){\n",
    "  W_old = W\n",
    "  for(i in 1:p){ # iterate column in covariance\n",
    "#     W11 = S[-i, -i] # use approximation as column-wise lasso for test\n",
    "    W11 = W[-i, -i]\n",
    "    \n",
    "    s12 = S[-i,i]\n",
    "\n",
    "    w12 = W[-i, i]\n",
    "    beta = solve(W11) %*% w12\n",
    "#     beta = B[, i]\n",
    "    \n",
    "    # coordinate decent\n",
    "    max_iter_lasso = 20\n",
    "    iter_lasso = 0\n",
    "    while(iter_lasso < max_iter_lasso){\n",
    "      beta_old = beta\n",
    "      for(j in 1:(p-1)){ # iterate entries in beta\n",
    "        uj = s12[j]\n",
    "        for(k in 1:(p-1)){\n",
    "          if(k==j) next\n",
    "          uj = uj - W11[k,j] * beta[k]\n",
    "\n",
    "        }\n",
    "#         cat(j, uj, '  ')\n",
    "        if(abs(uj) > rho){\n",
    "          beta[j] = (abs(uj) - rho) / W11[j,j]\n",
    "          if(uj < 0) beta[j] = -beta[j]\n",
    "        } \n",
    "        else{\n",
    "          beta[j] = 0\n",
    "        }\n",
    "      }\n",
    "#       cat('\\n')\n",
    "      iter_lasso = iter_lasso + 1\n",
    "      \n",
    "      change_lasso = max( abs(beta - beta_old) )\n",
    "      if(change_lasso<t) break\n",
    "    }\n",
    "#     cat('iter_lasso: ', iter_lasso, '\\n')\n",
    "    \n",
    "    w12 = W11 %*% beta\n",
    "#     cat(i, 'beta: ', beta, '\\n', '  w12:', w12,'\\n')\n",
    "    W[-i, i] = W[i, -i] = w12\n",
    "#     B[, i] = beta\n",
    "#     print(W)\n",
    "  }\n",
    "  \n",
    "  iter = iter + 1\n",
    "  change = max( abs(W-W_old) )\n",
    "  if(change < threshold) break\n",
    "  else{\n",
    "    cat('iter', iter, change, '\\n')\n",
    "  }\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "rational-court",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table class=\"dataframe\">\n",
       "<caption>A matrix: 10 × 10 of type dbl</caption>\n",
       "<tbody>\n",
       "\t<tr><td> 1.44131720</td><td>-0.05468838</td><td> 0.05354554</td><td>-0.00000086</td><td>-0.56346842</td><td> 0.17435579</td><td> 0.32708832</td><td>-0.71906876</td><td> 0.00000000</td><td> 0.00152736</td></tr>\n",
       "\t<tr><td>-0.05468838</td><td> 2.24621882</td><td> 0.52364729</td><td>-0.00000084</td><td>-0.00000017</td><td>-0.00000011</td><td> 0.22789042</td><td> 0.14648430</td><td> 0.21936794</td><td>-0.38582536</td></tr>\n",
       "\t<tr><td> 0.05354554</td><td> 0.52364729</td><td> 2.04851276</td><td>-0.92237944</td><td>-0.00000013</td><td>-0.07054451</td><td>-0.07351795</td><td>-0.39155108</td><td> 0.43813365</td><td>-0.24805204</td></tr>\n",
       "\t<tr><td>-0.00000086</td><td>-0.00000084</td><td>-0.92237944</td><td> 3.51590910</td><td>-0.32891514</td><td>-0.78028081</td><td> 0.00000000</td><td> 0.00000002</td><td> 0.00000002</td><td>-0.42883836</td></tr>\n",
       "\t<tr><td>-0.56346842</td><td>-0.00000017</td><td>-0.00000013</td><td>-0.32891514</td><td> 3.10161854</td><td>-0.15999947</td><td> 0.53073839</td><td>-0.27070848</td><td>-0.15020936</td><td>-0.00649441</td></tr>\n",
       "\t<tr><td> 0.17435579</td><td>-0.00000011</td><td>-0.07054451</td><td>-0.78028081</td><td>-0.15999947</td><td> 1.10494806</td><td>-0.11286726</td><td>-0.09716585</td><td>-0.32038924</td><td> 0.00000000</td></tr>\n",
       "\t<tr><td> 0.32708832</td><td> 0.22789042</td><td>-0.07351795</td><td> 0.00000000</td><td> 0.53073839</td><td>-0.11286726</td><td> 1.28585940</td><td> 0.00000000</td><td> 0.00000000</td><td> 0.00000000</td></tr>\n",
       "\t<tr><td>-0.71906876</td><td> 0.14648430</td><td>-0.39155108</td><td> 0.00000002</td><td>-0.27070848</td><td>-0.09716585</td><td> 0.00000000</td><td> 1.76994044</td><td>-0.03693921</td><td> 0.28791870</td></tr>\n",
       "\t<tr><td> 0.00000000</td><td> 0.21936794</td><td> 0.43813365</td><td> 0.00000002</td><td>-0.15020936</td><td>-0.32038924</td><td> 0.00000000</td><td>-0.03693921</td><td> 0.86196105</td><td> 0.05180987</td></tr>\n",
       "\t<tr><td> 0.00152736</td><td>-0.38582536</td><td>-0.24805204</td><td>-0.42883836</td><td>-0.00649441</td><td> 0.00000000</td><td> 0.00000000</td><td> 0.28791870</td><td> 0.05180987</td><td> 0.75370995</td></tr>\n",
       "</tbody>\n",
       "</table>\n"
      ],
      "text/latex": [
       "A matrix: 10 × 10 of type dbl\n",
       "\\begin{tabular}{llllllllll}\n",
       "\t  1.44131720 & -0.05468838 &  0.05354554 & -0.00000086 & -0.56346842 &  0.17435579 &  0.32708832 & -0.71906876 &  0.00000000 &  0.00152736\\\\\n",
       "\t -0.05468838 &  2.24621882 &  0.52364729 & -0.00000084 & -0.00000017 & -0.00000011 &  0.22789042 &  0.14648430 &  0.21936794 & -0.38582536\\\\\n",
       "\t  0.05354554 &  0.52364729 &  2.04851276 & -0.92237944 & -0.00000013 & -0.07054451 & -0.07351795 & -0.39155108 &  0.43813365 & -0.24805204\\\\\n",
       "\t -0.00000086 & -0.00000084 & -0.92237944 &  3.51590910 & -0.32891514 & -0.78028081 &  0.00000000 &  0.00000002 &  0.00000002 & -0.42883836\\\\\n",
       "\t -0.56346842 & -0.00000017 & -0.00000013 & -0.32891514 &  3.10161854 & -0.15999947 &  0.53073839 & -0.27070848 & -0.15020936 & -0.00649441\\\\\n",
       "\t  0.17435579 & -0.00000011 & -0.07054451 & -0.78028081 & -0.15999947 &  1.10494806 & -0.11286726 & -0.09716585 & -0.32038924 &  0.00000000\\\\\n",
       "\t  0.32708832 &  0.22789042 & -0.07351795 &  0.00000000 &  0.53073839 & -0.11286726 &  1.28585940 &  0.00000000 &  0.00000000 &  0.00000000\\\\\n",
       "\t -0.71906876 &  0.14648430 & -0.39155108 &  0.00000002 & -0.27070848 & -0.09716585 &  0.00000000 &  1.76994044 & -0.03693921 &  0.28791870\\\\\n",
       "\t  0.00000000 &  0.21936794 &  0.43813365 &  0.00000002 & -0.15020936 & -0.32038924 &  0.00000000 & -0.03693921 &  0.86196105 &  0.05180987\\\\\n",
       "\t  0.00152736 & -0.38582536 & -0.24805204 & -0.42883836 & -0.00649441 &  0.00000000 &  0.00000000 &  0.28791870 &  0.05180987 &  0.75370995\\\\\n",
       "\\end{tabular}\n"
      ],
      "text/markdown": [
       "\n",
       "A matrix: 10 × 10 of type dbl\n",
       "\n",
       "|  1.44131720 | -0.05468838 |  0.05354554 | -0.00000086 | -0.56346842 |  0.17435579 |  0.32708832 | -0.71906876 |  0.00000000 |  0.00152736 |\n",
       "| -0.05468838 |  2.24621882 |  0.52364729 | -0.00000084 | -0.00000017 | -0.00000011 |  0.22789042 |  0.14648430 |  0.21936794 | -0.38582536 |\n",
       "|  0.05354554 |  0.52364729 |  2.04851276 | -0.92237944 | -0.00000013 | -0.07054451 | -0.07351795 | -0.39155108 |  0.43813365 | -0.24805204 |\n",
       "| -0.00000086 | -0.00000084 | -0.92237944 |  3.51590910 | -0.32891514 | -0.78028081 |  0.00000000 |  0.00000002 |  0.00000002 | -0.42883836 |\n",
       "| -0.56346842 | -0.00000017 | -0.00000013 | -0.32891514 |  3.10161854 | -0.15999947 |  0.53073839 | -0.27070848 | -0.15020936 | -0.00649441 |\n",
       "|  0.17435579 | -0.00000011 | -0.07054451 | -0.78028081 | -0.15999947 |  1.10494806 | -0.11286726 | -0.09716585 | -0.32038924 |  0.00000000 |\n",
       "|  0.32708832 |  0.22789042 | -0.07351795 |  0.00000000 |  0.53073839 | -0.11286726 |  1.28585940 |  0.00000000 |  0.00000000 |  0.00000000 |\n",
       "| -0.71906876 |  0.14648430 | -0.39155108 |  0.00000002 | -0.27070848 | -0.09716585 |  0.00000000 |  1.76994044 | -0.03693921 |  0.28791870 |\n",
       "|  0.00000000 |  0.21936794 |  0.43813365 |  0.00000002 | -0.15020936 | -0.32038924 |  0.00000000 | -0.03693921 |  0.86196105 |  0.05180987 |\n",
       "|  0.00152736 | -0.38582536 | -0.24805204 | -0.42883836 | -0.00649441 |  0.00000000 |  0.00000000 |  0.28791870 |  0.05180987 |  0.75370995 |\n",
       "\n"
      ],
      "text/plain": [
       "      [,1]        [,2]        [,3]        [,4]        [,5]        [,6]       \n",
       " [1,]  1.44131720 -0.05468838  0.05354554 -0.00000086 -0.56346842  0.17435579\n",
       " [2,] -0.05468838  2.24621882  0.52364729 -0.00000084 -0.00000017 -0.00000011\n",
       " [3,]  0.05354554  0.52364729  2.04851276 -0.92237944 -0.00000013 -0.07054451\n",
       " [4,] -0.00000086 -0.00000084 -0.92237944  3.51590910 -0.32891514 -0.78028081\n",
       " [5,] -0.56346842 -0.00000017 -0.00000013 -0.32891514  3.10161854 -0.15999947\n",
       " [6,]  0.17435579 -0.00000011 -0.07054451 -0.78028081 -0.15999947  1.10494806\n",
       " [7,]  0.32708832  0.22789042 -0.07351795  0.00000000  0.53073839 -0.11286726\n",
       " [8,] -0.71906876  0.14648430 -0.39155108  0.00000002 -0.27070848 -0.09716585\n",
       " [9,]  0.00000000  0.21936794  0.43813365  0.00000002 -0.15020936 -0.32038924\n",
       "[10,]  0.00152736 -0.38582536 -0.24805204 -0.42883836 -0.00649441  0.00000000\n",
       "      [,7]        [,8]        [,9]        [,10]      \n",
       " [1,]  0.32708832 -0.71906876  0.00000000  0.00152736\n",
       " [2,]  0.22789042  0.14648430  0.21936794 -0.38582536\n",
       " [3,] -0.07351795 -0.39155108  0.43813365 -0.24805204\n",
       " [4,]  0.00000000  0.00000002  0.00000002 -0.42883836\n",
       " [5,]  0.53073839 -0.27070848 -0.15020936 -0.00649441\n",
       " [6,] -0.11286726 -0.09716585 -0.32038924  0.00000000\n",
       " [7,]  1.28585940  0.00000000  0.00000000  0.00000000\n",
       " [8,]  0.00000000  1.76994044 -0.03693921  0.28791870\n",
       " [9,]  0.00000000 -0.03693921  0.86196105  0.05180987\n",
       "[10,]  0.00000000  0.28791870  0.05180987  0.75370995"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "Theta_glasso = solve(W)\n",
    "round(Theta_glasso,8)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "technological-standing",
   "metadata": {},
   "source": [
    "### Compare with the outcome of original implementaion of R package glasso (by Fortran internally)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "restricted-campaign",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table class=\"dataframe\">\n",
       "<caption>A matrix: 10 × 10 of type dbl</caption>\n",
       "<tbody>\n",
       "\t<tr><td> 1.441318313</td><td>-0.05468993</td><td> 0.05354667</td><td> 0.0000000</td><td>-0.563458507</td><td> 0.17435728</td><td> 0.32708792</td><td>-0.71907706</td><td> 0.00000000</td><td> 0.001524034</td></tr>\n",
       "\t<tr><td>-0.054688765</td><td> 2.24621888</td><td> 0.52364641</td><td> 0.0000000</td><td> 0.000000000</td><td> 0.00000000</td><td> 0.22788994</td><td> 0.14648668</td><td> 0.21936792</td><td>-0.385825164</td></tr>\n",
       "\t<tr><td> 0.053546259</td><td> 0.52364592</td><td> 2.04851200</td><td>-0.9223762</td><td> 0.000000000</td><td>-0.07054298</td><td>-0.07351858</td><td>-0.39155091</td><td> 0.43813384</td><td>-0.248052595</td></tr>\n",
       "\t<tr><td> 0.000000000</td><td> 0.00000000</td><td>-0.92237507</td><td> 3.5159082</td><td>-0.328924874</td><td>-0.78028279</td><td> 0.00000000</td><td> 0.00000000</td><td> 0.00000000</td><td>-0.428838587</td></tr>\n",
       "\t<tr><td>-0.563469507</td><td> 0.00000000</td><td> 0.00000000</td><td>-0.3289192</td><td> 3.101608921</td><td>-0.15999533</td><td> 0.53073901</td><td>-0.27071061</td><td>-0.15021093</td><td>-0.006493087</td></tr>\n",
       "\t<tr><td> 0.174353664</td><td> 0.00000000</td><td>-0.07054566</td><td>-0.7802807</td><td>-0.159991881</td><td> 1.10494852</td><td>-0.11286548</td><td>-0.09717106</td><td>-0.32038930</td><td> 0.000000000</td></tr>\n",
       "\t<tr><td> 0.327090480</td><td> 0.22789096</td><td>-0.07351854</td><td> 0.0000000</td><td> 0.530737365</td><td>-0.11286751</td><td> 1.28585876</td><td> 0.00000000</td><td> 0.00000000</td><td> 0.000000000</td></tr>\n",
       "\t<tr><td>-0.719067311</td><td> 0.14648621</td><td>-0.39155321</td><td> 0.0000000</td><td>-0.270710191</td><td>-0.09716916</td><td> 0.00000000</td><td> 1.76994907</td><td>-0.03693925</td><td> 0.287920282</td></tr>\n",
       "\t<tr><td> 0.000000000</td><td> 0.21936760</td><td> 0.43813413</td><td> 0.0000000</td><td>-0.150208957</td><td>-0.32038916</td><td> 0.00000000</td><td>-0.03693588</td><td> 0.86196136</td><td> 0.051809602</td></tr>\n",
       "\t<tr><td> 0.001529211</td><td>-0.38582500</td><td>-0.24805293</td><td>-0.4288401</td><td>-0.006491063</td><td> 0.00000000</td><td> 0.00000000</td><td> 0.28791951</td><td> 0.05180980</td><td> 0.753710017</td></tr>\n",
       "</tbody>\n",
       "</table>\n"
      ],
      "text/latex": [
       "A matrix: 10 × 10 of type dbl\n",
       "\\begin{tabular}{llllllllll}\n",
       "\t  1.441318313 & -0.05468993 &  0.05354667 &  0.0000000 & -0.563458507 &  0.17435728 &  0.32708792 & -0.71907706 &  0.00000000 &  0.001524034\\\\\n",
       "\t -0.054688765 &  2.24621888 &  0.52364641 &  0.0000000 &  0.000000000 &  0.00000000 &  0.22788994 &  0.14648668 &  0.21936792 & -0.385825164\\\\\n",
       "\t  0.053546259 &  0.52364592 &  2.04851200 & -0.9223762 &  0.000000000 & -0.07054298 & -0.07351858 & -0.39155091 &  0.43813384 & -0.248052595\\\\\n",
       "\t  0.000000000 &  0.00000000 & -0.92237507 &  3.5159082 & -0.328924874 & -0.78028279 &  0.00000000 &  0.00000000 &  0.00000000 & -0.428838587\\\\\n",
       "\t -0.563469507 &  0.00000000 &  0.00000000 & -0.3289192 &  3.101608921 & -0.15999533 &  0.53073901 & -0.27071061 & -0.15021093 & -0.006493087\\\\\n",
       "\t  0.174353664 &  0.00000000 & -0.07054566 & -0.7802807 & -0.159991881 &  1.10494852 & -0.11286548 & -0.09717106 & -0.32038930 &  0.000000000\\\\\n",
       "\t  0.327090480 &  0.22789096 & -0.07351854 &  0.0000000 &  0.530737365 & -0.11286751 &  1.28585876 &  0.00000000 &  0.00000000 &  0.000000000\\\\\n",
       "\t -0.719067311 &  0.14648621 & -0.39155321 &  0.0000000 & -0.270710191 & -0.09716916 &  0.00000000 &  1.76994907 & -0.03693925 &  0.287920282\\\\\n",
       "\t  0.000000000 &  0.21936760 &  0.43813413 &  0.0000000 & -0.150208957 & -0.32038916 &  0.00000000 & -0.03693588 &  0.86196136 &  0.051809602\\\\\n",
       "\t  0.001529211 & -0.38582500 & -0.24805293 & -0.4288401 & -0.006491063 &  0.00000000 &  0.00000000 &  0.28791951 &  0.05180980 &  0.753710017\\\\\n",
       "\\end{tabular}\n"
      ],
      "text/markdown": [
       "\n",
       "A matrix: 10 × 10 of type dbl\n",
       "\n",
       "|  1.441318313 | -0.05468993 |  0.05354667 |  0.0000000 | -0.563458507 |  0.17435728 |  0.32708792 | -0.71907706 |  0.00000000 |  0.001524034 |\n",
       "| -0.054688765 |  2.24621888 |  0.52364641 |  0.0000000 |  0.000000000 |  0.00000000 |  0.22788994 |  0.14648668 |  0.21936792 | -0.385825164 |\n",
       "|  0.053546259 |  0.52364592 |  2.04851200 | -0.9223762 |  0.000000000 | -0.07054298 | -0.07351858 | -0.39155091 |  0.43813384 | -0.248052595 |\n",
       "|  0.000000000 |  0.00000000 | -0.92237507 |  3.5159082 | -0.328924874 | -0.78028279 |  0.00000000 |  0.00000000 |  0.00000000 | -0.428838587 |\n",
       "| -0.563469507 |  0.00000000 |  0.00000000 | -0.3289192 |  3.101608921 | -0.15999533 |  0.53073901 | -0.27071061 | -0.15021093 | -0.006493087 |\n",
       "|  0.174353664 |  0.00000000 | -0.07054566 | -0.7802807 | -0.159991881 |  1.10494852 | -0.11286548 | -0.09717106 | -0.32038930 |  0.000000000 |\n",
       "|  0.327090480 |  0.22789096 | -0.07351854 |  0.0000000 |  0.530737365 | -0.11286751 |  1.28585876 |  0.00000000 |  0.00000000 |  0.000000000 |\n",
       "| -0.719067311 |  0.14648621 | -0.39155321 |  0.0000000 | -0.270710191 | -0.09716916 |  0.00000000 |  1.76994907 | -0.03693925 |  0.287920282 |\n",
       "|  0.000000000 |  0.21936760 |  0.43813413 |  0.0000000 | -0.150208957 | -0.32038916 |  0.00000000 | -0.03693588 |  0.86196136 |  0.051809602 |\n",
       "|  0.001529211 | -0.38582500 | -0.24805293 | -0.4288401 | -0.006491063 |  0.00000000 |  0.00000000 |  0.28791951 |  0.05180980 |  0.753710017 |\n",
       "\n"
      ],
      "text/plain": [
       "      [,1]         [,2]        [,3]        [,4]       [,5]         [,6]       \n",
       " [1,]  1.441318313 -0.05468993  0.05354667  0.0000000 -0.563458507  0.17435728\n",
       " [2,] -0.054688765  2.24621888  0.52364641  0.0000000  0.000000000  0.00000000\n",
       " [3,]  0.053546259  0.52364592  2.04851200 -0.9223762  0.000000000 -0.07054298\n",
       " [4,]  0.000000000  0.00000000 -0.92237507  3.5159082 -0.328924874 -0.78028279\n",
       " [5,] -0.563469507  0.00000000  0.00000000 -0.3289192  3.101608921 -0.15999533\n",
       " [6,]  0.174353664  0.00000000 -0.07054566 -0.7802807 -0.159991881  1.10494852\n",
       " [7,]  0.327090480  0.22789096 -0.07351854  0.0000000  0.530737365 -0.11286751\n",
       " [8,] -0.719067311  0.14648621 -0.39155321  0.0000000 -0.270710191 -0.09716916\n",
       " [9,]  0.000000000  0.21936760  0.43813413  0.0000000 -0.150208957 -0.32038916\n",
       "[10,]  0.001529211 -0.38582500 -0.24805293 -0.4288401 -0.006491063  0.00000000\n",
       "      [,7]        [,8]        [,9]        [,10]       \n",
       " [1,]  0.32708792 -0.71907706  0.00000000  0.001524034\n",
       " [2,]  0.22788994  0.14648668  0.21936792 -0.385825164\n",
       " [3,] -0.07351858 -0.39155091  0.43813384 -0.248052595\n",
       " [4,]  0.00000000  0.00000000  0.00000000 -0.428838587\n",
       " [5,]  0.53073901 -0.27071061 -0.15021093 -0.006493087\n",
       " [6,] -0.11286548 -0.09717106 -0.32038930  0.000000000\n",
       " [7,]  1.28585876  0.00000000  0.00000000  0.000000000\n",
       " [8,]  0.00000000  1.76994907 -0.03693925  0.287920282\n",
       " [9,]  0.00000000 -0.03693588  0.86196136  0.051809602\n",
       "[10,]  0.00000000  0.28791951  0.05180980  0.753710017"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "glasso::glasso(S, rho=rho)$wi"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "amber-temple",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "R",
   "language": "R",
   "name": "ir"
  },
  "language_info": {
   "codemirror_mode": "r",
   "file_extension": ".r",
   "mimetype": "text/x-r-source",
   "name": "R",
   "pygments_lexer": "r",
   "version": "3.6.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}