remarkRemark hypothesisHypothesis claimClaim preconditioned ISTA for Graphical
LASSOG. Shalom, E. Treister, and I. Yavneh

# pISTA: preconditioned Iterative Soft Thresholding Algorithm for Graphical
Lasso

Gal Shalom Faculty of Computer Science, Technion — Israel Institute of
Technology, Israel (, ). Supported in part by the Israel Science Foundation,
grant No. 1639/19<EMAIL_ADDRESS><EMAIL_ADDRESS>Eran
Treister Department of Computer Science, Ben-Gurion University of the Negev,
Israel (). This work was supported in part by the Israeli Council for Higher
Education (CHE) via the Data Science Research Center, Ben-Gurion University of
the Negev, Israel<EMAIL_ADDRESS>Irad Yavneh11footnotemark: 1

###### Abstract

We propose a novel quasi-Newton method for solving the sparse inverse
covariance estimation problem also known as the graphical least absolute
shrinkage and selection operator (GLASSO). This problem is often solved using
a second order quadratic approximation. However, in such algorithms the
Hessian term is complex and computationally expensive to handle. Therefore,
our method uses the inverse of the Hessian as a preconditioner to simplify and
approximate the quadratic element at the cost of a more complex $\ell_{1}$
element. The variables of the resulting preconditioned problem are coupled
only by the $\ell_{1}$ sub-derivative of each other, which can be guessed with
minimal cost using the gradient itself, allowing the algorithm to be
parallelized and implemented efficiently on GPU hardware accelerators.
Numerical results on synthetic and real data demonstrate that our method is
competitive with other state-of-the-art approaches.

###### keywords:

Graphical LASSO, Sparse Precision Matrix Estimation, Proximal Methods,
Preconditioning.

90C25, 65D18, 65K10, 65F08

##  1 Introduction

Inverse covariance estimation is a fundamental problem in modern statistics.
Specifically, the inverse covariance matrices of multivariate normal
distributions are used in numerous applications. One of the most common uses
of the inverse covariance, also known as the Precision Matrix, is to describe
statistical models. That is, the graph structure of Gaussian graphical models
can be inferred from the Precision Matrix [19]. This inferred graph can be
used for analyzing gene networks [8], financial assets and stocks dependencies
[10, 30], social networks [13] and other inter-dependent data [5, 18].

The most straightforward approach to estimating the Precision Matrix would be
to use the inverse of the empirical covariance matrix. Evidently, inverting
the empirical covariance matrix requires a number of samples that equals or
exceeds the dimension of the matrix. However, many problems of interest are
high dimensional [17, 8], to the extent that they are significantly larger
than the number of available samples, so the empirical covariance matrix is
rank-deficient. In such cases, some kind of regularization is essential. One
of the most common ways to estimate the Precision Matrix is to solve a
$\ell_{1}$-regularized maximum likelihood estimation problem, known as the
graphical least absolute shrinkage and selection operator (Graphical LASSO, or
GLASSO) problem. The use of $\ell_{1}$-regularization aims to achieve a sparse
estimation while keeping the problem convex. Sparsity of the Precision Matrix
may be interpreted as simplicity—a sparse Precision Matrix implies a simple
inferred graph structure [15]. Thus, when a genuine empirical Precision Matrix
cannot be computed, a sparse valid estimation is often the preferred choice
[15, 31, 7]. We refer the reader to the first chapter of [15] and to [7] for
additional motivation.

As the problem is well studied, many algorithms and methods for its solution
have been presented over the years [1, 2, 6, 17, 14, 27, 12, 16, 26, 24, 22,
9, 21]. Some of these methods are computationally intensive or difficult to
parallelize, rendering parallel computation accelerators such as GPUs
ineffective. As far as we know, none of the algorithms in the literature were
designed specifically with an efficient GPU deployment in mind. Nevertheless,
various existing algorithms are able to take advantage of GPU processing power
with a simple implementation, e.g., G-ISTA [14], VSM [21], PSM [9], ALM [24],
Newton-Lasso [22] and Orthant-Based Newton [22]. Each of these has its
strengths and weaknesses. Other algorithms, such as BCD-IC [26], can be
deployed on a GPU as well, however, the method and the authors’ implementation
are complex, and the extensive usage of scalar operations may hinder its
performance.

In this work we introduce a preconditioned Iterative Soft Thresholding
Algorithm (pISTA) for solving the graphical LASSO problem. pISTA is a quasi-
Newton algorithm designed to be highly parallel and suitable for using GPU
capabilities to its benefit by exploiting the problem structure effectively. A
traditional second order quadratic approximation includes a complex Hessian
term which is computationally expensive to handle. Moreover, all the elements
of the Hessian term are coupled, making efficient GPU deployment challenging.
To this end, pISTA uses the inverse of the Hessian as a preconditioner to
simplify and approximate the quadratic element, at the cost of a more complex
$\ell_{1}$ element. This is highly beneficial since the Hessian inverse is
easily obtained for the GLASSO problem. The variables of the resulting
preconditioned problem are coupled only by the $\ell_{1}$ sub-derivative of
each other, which can be guessed with minimal cost using the gradient itself.
The resulting pISTA algorithm uses only matrix operations, making it easier
for exploiting the GPU efficiently.

This paper is organized as follows. In Section 2 we refer to related work, and
in Section 3 we formulate the problem. The pISTA algorithm is introduced in
Section 4, and its convergence is proved in Section 5. Experimental results
are presented in Section 6, and conclusions follow in Section 7.

##  2 Related Work

Many numerical algorithms have been developed, which are specifically designed
for solving the sparse inverse covariance estimation problem. G-ISTA [14] uses
a proximal gradient method which estimates the inverse covariance iteratively.
GLASSO [12] splits the problem into smaller LASSO problems and updates the
estimation by solving them separately. QUIC [16] uses a proximal Newton method
on the objective function. BCD-IC [26] uses a block coordinate descent method.
ALM [24] uses an alternating linearization technique by splitting the
objective function into two linear approximations and minimizing them
alternatingly. PSM [9] employs a projected gradient method. VSM [21] uses
Nesterov’s smooth optimization technique. Newton-Lasso [22] and Orthant-Based
Newton [22] use FISTA [3] and Conjugate Gradient methods, respectively, to
solve reduced second order approximations. Although some algorithms are
designed to be run in parallel using the multi-core model of the CPU, none of
them take into specific consideration the parallelism model of the GPU and the
usage of GPUs as a computation accelerator.

Additionally, some algorithm were developed specifically for large scale
matrices where the matrix can only reside in memory in sparse format. Big&Quic
and SQUIC [17, 4] extend QUIC [16] to large scales. Big&Quic uses on-demand
computation of the Hessian’s columns and a special procedure of the linesearch
conditions using Schur complements. SQUIC utilizes sparse Cholesky factors of
the iterations matrix. ML-BCD [27] takes BCD-IC [26], which is already
suitable to large scale problems, and accelerates it using a generic
multilevel framework, which was originally suggested in [28] for LASSO. This
acceleration can in principle be used with our proposed method as well. Large
scale algorithms may benefit from the GPU processing power, however memory
constraints and usage of sparse format make this task non-trivial.

Graphical LASSO can also be used in a mixture model setup when multiple sparse
inverse covariances are utilized [11]. This allows a richer statistical model,
but requires the estimation of multiple sparse Precision Matrices. An
efficient GPU-based solver is highly beneficial is such scenarios as well.

##  3 Background

###  3.1 Sparse Inverse Covariance Estimation

Estimating the parameters of multivariate Gaussian distributions is a
fundamental problem in statistics. Given $m$ independent samples
$\\{y_{i}\\}_{i=1}^{m}\in\mathbb{R}^{n}$, where
$y_{i}\sim\mathcal{N}(\mu,\Sigma)$, one would like to estimate the mean
$\mu\in\mathbb{R}^{n}$ and either the covariance $\Sigma$, or its inverse
$\Sigma^{-1}$, which is also called the Precision Matrix. Both the mean $\mu$
and the covariance $\Sigma$ are often estimated by the maximum likelihood
estimator (MLE). The MLE is given by the parameters $\mu$ and $\Sigma$ that
maximize the probability to sample the observed data $\\{y_{i}\\}_{i=1}^{m}$:

(1)
$\begin{split}\hat{\mu},\hat{\Sigma}&=\operatorname*{arg\,max}_{\mu,\Sigma}\prod^{m}_{i=1}\mathbb{P}(y_{i}|\Sigma,\mu)\\\
&=\operatorname*{arg\,max}_{\mu,\Sigma}\prod^{m}_{i=1}\frac{1}{\sqrt{(2\pi)^{m}det(\Sigma)}}\text{exp}\left(-\frac{1}{2}(y_{i}-\mu)^{T}\Sigma^{-1}(y_{i}-\mu)\right)\,.\end{split}$

MLE has an analytical solution:

(2)
$\hat{\mu}=\frac{1}{m}\sum_{i=1}^{m}y_{i},\quad\quad\hat{\Sigma}=S\overset{\Delta}{=}\frac{1}{m}\sum_{i=1}^{m}(y_{i}-\hat{\mu})(y_{i}-\hat{\mu})^{T}\,,$

where $\hat{\mu}$ and $\hat{\Sigma}$ are called the empirical mean and
empirical covariance, respectively.

Usually, one considers the log-likelihood objective (negative $log$ of Eq. 1),
minimized over the inverse covariance matrix, yielding the inverse covariance
MLE:

(3) $\hat{\Sigma}^{-1}=\operatorname*{arg\,min}_{A\succ
0}f(A)\overset{\Delta}{=}-\log(\det(A))+Tr(SA)\,.$

The solution is indeed the inverse of Eq. 2. However, in cases where the
number of available samples is smaller than the dimension of $y_{i}$ ($m<n$),
the matrix $S$ is rank-deficient and thus non-invertible, whereas the true
$\Sigma$ is assumed to have full rank and is positive definite. In other
words, we cannot estimate $\Sigma^{-1}$ by inverting $S$, and further
assumptions should be considered. A common choice in this case is to assume
that $\Sigma^{-1}$ is sparse. $\Sigma^{-1}$ can be interpreted as a
conditional dependence matrix, that is, its off-diagonal entries indicate the
dependence between the row and column variables, given all the remaining
variables [23].

A common approach is to regularize the log-likelihood objective Eq. 3 with a
sparsity promoting $\ell_{1}$-penalty [1, 2, 6]:

(4) $\hat{\Sigma}^{-1}=\operatorname*{arg\,min}_{A\succ
0}F(A)\overset{\Delta}{=}\operatorname*{arg\,min}_{A\succ
0}f(A)+\alpha||A||_{1}\,,$

where $\alpha>0$ is a scalar and $||A||_{1}=\sum_{i,j}|A_{i,j}|$. Equation 4
is known as the _Graphical Lasso_ problem.

As seen from Eq. 4, the objective function is convex and is composed of two
parts—a smooth convex part $f(A)$ and a non-smooth convex part
$\alpha||A||_{1}$. Although the objective is convex, the non-smooth term makes
traditional algorithms ineffective, and more specialized solvers are needed.

###  3.2 Proximal Methods for Sparse Inverse Covariance Estimation

A common approach to solving convex problems comprised of smooth and non-
smooth parts is known as proximal methods, where we approximate the objective
function at each iteration as follows: the smooth part is approximated by a
quadratic function while the non-smooth part is kept unchanged. This approach
is especially attractive when the non-smooth part is a separable function
(e.g., point-wise).

Specifically, applying proximal methods to the objective function $F(A)$ of
Eq. 4, the smooth term $f(A)$ is approximated by a quadratic function while
the term $\alpha||A||_{1}$ remains unchanged. This approximation yields a
linear LASSO [25] problem ($\ell_{1}$-regularized quadratic objective). Once
defined, we minimize it approximately using a LASSO solver. Denoting the
descent direction at the $k$-th iteration by $D^{(k)}$, then at each iteration
we solve:

(5) $\begin{split}D^{(k)}&=\operatorname*{arg\,min}_{D}\tilde{F}(A^{(k)}+D)\\\
&=\operatorname*{arg\,min}_{D}f(A^{(k)})+\left\langle
g^{(k)},D\right\rangle+\frac{1}{2}\left\langle
D,H^{(k)}D\right\rangle+\alpha||A^{(k)}+D||_{1}\,,\end{split}$

where $g^{(k)}=\nabla f(A^{(k)})$ is the gradient at the $k$th iteration, and
$H^{(k)}$ depends on the particular method. For $H^{(k)}=I$ we get the
“Proximal Gradient” method, which is called G-ISTA [14]. At the other extreme,
for $H^{(k)}=\nabla^{2}f(A^{(k)})$, the Hessian at the $k$th iteration, we get
the “Proximal Newton” method known as QUIC [16].

The gradient and the Hessian are given by:

(6) $\nabla f(A)=S-A^{-1}\quad,\quad\nabla^{2}f(A)=A^{-1}\otimes A^{-1}\,,$

where $\otimes$ is the Kronecker product. We note that once the gradient has
been computed, we can compute Hessian-vector products at relatively low cost
by using Kronecker product properties. However, solving Eq. 5 with
$H^{(k)}=\nabla^{2}f(A^{(k)})$ is highly complex and it is typically done by
coordinate descent iterations [16]. This makes it hard to use computation
accelerators (such as GPUs) efficiently.

##  4 The preconditioned Iterative Soft Thresholding Algorithm (pISTA)

In this section we develop the preconditioned Iterative Soft Thresholding
Algorithm (pISTA). In [16], the Proximal Newton method for our problem, each
coordinate computation depends on the values of all the other coordinates due
to the Hessian in the quadratic element. In pISTA, we aim for a simpler and
easier to solve quadratic element at the cost of a more complex non-smooth
$\ell_{1}$ part.

First, in Section 4.1, we restrict our problem to a smaller set of variables
as done in other algorithms [16, 27, 22]. In Section 4.2, we simplify and
relax the quadratic part using a preconditioner, at the cost of making the
non-smooth part $\alpha||A||_{1}$ more complex. Lastly, in Section 4.3, we
solve the resulting problem, completing the development of pISTA. In contrast
to second order methods like QUIC, for example, in pISTA each coordinate
depends on the other variables only through the sub-gradient of the non-smooth
part ($\ell_{1}$ norm), which is more complex. However, the sub-gradient of
the $\ell_{1}$ term can be approximated easily, and relatively well, by the
sign of the elements of the current iterate (if it is nonzero) or of its
gradient (where the elements of the iterate vanish), allowing our quasi-Newton
method to be efficient.

###  4.1 Restricting the Updates to the Free-set

As introduced in [16] and used in [22, 27], we restrict the descent direction
$D$ to the free-set of $A^{(k)}$. That is, an element in $D$ which is not in
the free-set of $A^{(k)}$ is set to zero. Denote by $\mathcal{S}_{A}^{(k)}$
the free-set at iteration $k$, defined as follows:

(7) $\displaystyle\mathcal{S}_{A}^{(k)}=\bigg{\\{}(i,j)\big{|}\
A^{(k)}_{i,j}\neq 0\bigg{\\}}\cup\bigg{\\{}(i,j)\big{|}\ |\nabla
f(A^{(k)}|_{i,j}>\alpha\bigg{\\}}\,.$

According to Lemma 7 of [16], solving Eq. 5 restricted to the variables that
are not in the free-set of $A^{(k)}$ will result in a zero value in every
element in the direction $D$. Thus, restricting $D$ to the free-set is
equivalent to solving Eq. 5 in an alternating two-block manner: the first step
restricted to the variables not in the free-set (which makes no change in our
approximation), and the second step restricted to the free-set. The free set
is recomputed before the next iteration.

Define the restriction mask at the $k$-th iteration:

(8)
$[\mathcal{M}_{A}^{(k)}]_{i,j}=\bigg{\\{}\begin{array}[]{lr}1&(i,j)\in\mathcal{S}_{A}^{(k)}\\\
0&(i,j)\notin\mathcal{S}_{A}^{(k)}\end{array}\,.$

Then, the restricted problem we solve at the $k$th iteration can be written
as:

(9)
$\begin{split}D^{(k)}&=\operatorname*{arg\,min}_{D}\tilde{F}(A^{(k)}+(\mathcal{M}_{A}^{(k)}\odot
D))\end{split}\,,$

where $\odot$ is the Hadamard product. Note that, in the solution of Eq. 9,
the elements of $D^{(k)}$ that are not in the free-set of $A^{(k)}$ may have
nonzero values. However, we restrict the updates to the indices in the free-
set of $A^{(k)}$ only, therefore, our descent direction is
$\mathcal{M}_{A}^{(k)}\odot D^{(k)}$, which is equivalent to setting the
values of $D^{(k)}$ that are not in the free-set to zero.

###  4.2 Preconditioning of the Descent Direction

Note the following property of the Hessian in Eq. 6:

(10) $\left(\nabla^{2}f(A)\right)^{-1}=A\otimes A\,,$

where $\otimes$ is the Kronecker product. This means that the Hessian inverse
can be obtained without any computational overhead, using the typically sparse
$A$ instead of the typically dense $A^{-1}$. Here we are interested in using
this appealing property to accelerate the solution of Eq. 5, where
$H^{(k)}=\nabla^{2}f(A^{(k)})$. To this end, we use it as a preconditioner to
$D$:

(11)
$\text{vec}(D)=\big{(}\nabla^{2}f(A^{(k)})\big{)}^{-1}\text{vec}(\Delta)\Rightarrow
D=A^{(k)}\Delta A^{(k)}\,,$

where vec denotes the column-stacking of a matrix into a vector.

Denote $W^{(k)}=(A^{(k)})^{-1}$. Using the preconditioned descent direction
Eq. 11 in the restricted problem Eq. 9, with
$H^{(k)}=\nabla^{2}f(A^{(k)})=W^{(k)}\otimes W^{(k)}$, results in the
following equation:

(12)
$\begin{split}\Delta^{(k)}&=\operatorname*{arg\,min}_{\Delta}\tilde{F}\left(A^{(k)}+\mathcal{M}_{A}^{(k)}\odot(A^{(k)}\Delta
A^{(k)})\right)\\\ &=\operatorname*{arg\,min}_{\Delta}f(A^{(k)})+\left\langle
g^{(k)},\mathcal{M}_{A}^{(k)}\odot(A^{(k)}\Delta A^{(k)})\right\rangle\\\
&\quad\quad+\frac{1}{2}\left\langle\mathcal{M}_{A}^{(k)}\odot(A^{(k)}\Delta
A^{(k)}),W^{(k)}\left(\mathcal{M}_{A}^{(k)}\odot(A^{(k)}\Delta
A^{(k)})\right)W^{(k)}\right\rangle\\\
&\quad\quad+\alpha\left|\left|A^{(k)}+\mathcal{M}_{A}^{(k)}\odot(A^{(k)}\Delta
A^{(k)})\right|\right|_{1}\,.\end{split}$

After applying the mask restriction $\mathcal{M}$ and the preconditioning to
Eq. 5, the quadratic part (the second element in Eq. 12) is still complex. We
relax it by removing the mask restriction $\mathcal{M}$ and replacing it with
a multiplicative scalar:

(13) $\frac{1}{2}\left\langle\mathcal{M}_{A}^{(k)}\odot(A^{(k)}\Delta
A^{(k)}),W^{(k)}\left(\mathcal{M}_{A}^{(k)}\odot(A^{(k)}\Delta
A^{(k)})\right)W^{(k)}\right\rangle\\\ \approx\frac{1}{2t}\left\langle
A^{(k)}\Delta A^{(k)},W^{(k)}(A^{(k)}\Delta
A^{(k)})W^{(k)}\right\rangle=\frac{1}{2t}\left\langle A^{(k)}\Delta
A^{(k)},\Delta\right\rangle\,,$

where $t$ is a scalar that is computed using linesearch. Since we do not
change the first order terms (the gradient and $\ell_{1}$), this can still
result in a monotonically convergent method. Denote the resulting quasi-Newton
approximation by $P(\Delta;t)$, then the $k$-th iterate is given by:

(14)
$\begin{split}\Delta^{(k)}&=\operatorname*{arg\,min}_{\Delta}P(\Delta;t)\\\
&=\operatorname*{arg\,min}_{\Delta}f(A^{(k)})+\left\langle
g^{(k)},\mathcal{M}_{A}^{(k)}\odot(A^{(k)}\Delta
A^{(k)})\right\rangle+\frac{1}{2t}\left\langle A^{(k)}\Delta
A^{(k)},\Delta\right\rangle\\\
&\quad\quad+\alpha\left|\left|A^{(k)}+\mathcal{M}_{A}^{(k)}\odot(A^{(k)}\Delta
A^{(k)})\right|\right|_{1}\,.\end{split}$

###  4.3 pISTA

Our pISTA algorithm solves the problem defined by Eq. 14 in each iteration,
with an appropriately chosen $t$. First, we find a solution $\Delta^{(k)}$ to
Eq. 14. To obtain this, we need to formulate the sub-differential of the non-
smooth $\ell_{1}$ term $||A+D||_{1}$. First, denote the sub-differential of
$||A+D||_{1}$ by:

(15) $\displaystyle G^{(k)}(D)$ $\displaystyle\overset{\Delta}{=}$
$\displaystyle\frac{\partial||A^{(k)}+D||_{1}}{\partial D}$ (18)
$\displaystyle=$
$\displaystyle\Bigg{\\{}\mathcal{G}^{(k)}\Bigg{|}\begin{array}[]{lr}\mathcal{G}^{(k)}_{i,j}=[sign(A^{(k)}+D)]_{i,j}&[A^{(k)}+D]_{i,j}\neq
0\\\
\mathcal{G}^{(k)}_{i,j}\in[-1,1]&[A^{(k)}+D]_{i,j}=0\end{array}\Bigg{\\}}\,.$

Next, using the chain rule and some known derivatives with Eq. 15, we
formulate the sub-differential of the quasi-Newton approximation:

(19) $\frac{\partial
P(\Delta;t)}{\partial\Delta}=\Bigg{\\{}\frac{1}{t}A^{(k)}\Delta
A^{(k)}+A^{(k)}\left(g^{(k)}\odot\mathcal{M}_{A}^{(k)}\right)A^{(k)}+\alpha
A^{(k)}\left(\mathcal{G}^{(k)}\odot\mathcal{M}_{A}^{(k)}\right)A^{(k)}\\\
\Bigg{|}\ \mathcal{G}^{(k)}\in
G^{(k)}\left(\mathcal{M}_{A}^{(k)}\odot(A^{(k)}\Delta
A^{(k)})\right)\Bigg{\\}}\,,$

where $\mathcal{G}^{(k)}$ represents a sub-gradient of the $\ell_{1}$ term.

The desired $\Delta^{(k)}$ that will be used to compute the descent direction
is one which includes the matrix 0 in its sub-gradients:

(20) $\Delta^{(k)}:\quad\quad 0\in\frac{\partial
P(\Delta^{(k)};t)}{\partial\Delta}\,.$

Notice that after finding $\Delta^{(k)}$, it will be used to compute the
descent direction $\mathcal{M}_{A}^{(k)}\odot
D^{(k)}=\mathcal{M}_{A}^{(k)}\odot(A^{(k)}\Delta^{(k)}A^{(k)})$. Thus, we can
equivalently find $D^{(k)}=A^{(k)}\Delta^{(k)}A^{(k)}$ instead of
$\Delta^{(k)}$. In other words, we shall find $D^{(k)}$ which has 0 as one of
its sub-gradients:

(21)
$D^{(k)}:0\in\Bigg{\\{}\frac{1}{t}D^{(k)}+A^{(k)}\left(g^{(k)}\odot\mathcal{M}_{A}^{(k)}\right)A^{(k)}+\alpha
A^{(k)}\left(\mathcal{G}^{(k)}\odot\mathcal{M}_{A}^{(k)}\right)A^{(k)}\\\
\Bigg{|}\ \mathcal{G}^{(k)}\in G^{(k)}\left(\mathcal{M}_{A}^{(k)}\odot
D^{(k)}\right)\Bigg{\\}}\,.$

Note that the change of variables has made the first (Hessian) term trivial.
To find $D^{(k)}$, we split the sub-gradient to $n^{2}$ equations where the
$(i,j)$-th equation is composed by the $(i,j)$ element in the sub-gradient.
For each equation $(i,j)$, we define the equation variable as $D^{(k)}_{i,j}$.
The approximate solution (developed in Section 7) is:

(22) $D^{(k)}_{i,j}=-A^{(k)}_{i,j}+\text{SoftThreshold}\left(A^{(k)}-t\cdot
B^{(k)}_{i,j}\ ,\ t\cdot C^{(k)}_{i,j}\right)\,,$

where

(23)
$C^{(k)}_{i,j}=\left\\{\begin{array}[]{lc}\alpha\cdot\left(A^{(k)}_{i,i}\cdot
A^{(k)}_{j,j}\right),&i=j\\\ \alpha\cdot\left(A^{(k)}_{i,i}\cdot
A^{(k)}_{j,j}+A^{(k)}_{i,j}\cdot A^{(k)}_{j,i}\right),&i\neq
j\end{array}\right.\,,$

and

(24)
$B^{(k)}=A^{(k)}\left(g^{(k)}\odot\mathcal{M}_{A}^{(k)}\right)A^{(k)}+\alpha
A^{(k)}\left(\mathcal{G}^{(k)}\odot\mathcal{M}_{A}^{(k)}\right)A^{(k)}-C^{(k)}\odot\left(\mathcal{G}^{(k)}\odot\mathcal{M}_{A}^{(k)}\right)\,.$

So far, we developed a second order approximation for the GLASSO using the
Hessian inverse and a quadratic relaxation. We found a closed form direction
$D^{(k)}$ for each $(i,j)$ entry, assuming that the rest of the entries in
$\mathcal{G}^{(k)}$ are given. In the following subsection we show how we
approximate those entries.

Algorithm 1 describes the pISTA method using Eq. 22. One part that is still
missing from the description so far is the computation of $\mathcal{G}^{(k)}$
and the linesearch procedure for finding an appropriate $t$, which we describe
next. It is important to note that all the steps and operations in Algorithm 1
can be computed using matrix operations. This property allows us to easily
benefit from GPU computation acceleration and easy GPU deployment.

Algorithm 1 pISTA($S,\alpha,A^{(0)}$)

1: Result: $A$

2: Initialization: $k=0$;

3: while stop criteria not met do

4: Compute $\mathcal{S}_{A}^{(k)}$ according to Eq. 7;

5: Compute $\mathcal{M}_{A}^{(k)}$ according to Eq. 8;

6: Compute $\mathcal{G}^{(k)}$ according to Eq. 25;

7: Find an appropriate $t$ and compute $D^{(K)}$ according to Eq. 22;

8: $A^{(k+1)}=A^{(K)}+\mathcal{M}_{A}^{(k)}\odot D^{(K)}$;

9: end while

####  4.3.1 The Approximation of $\mathcal{G}$

The algorithm depends on a good approximation of $\mathcal{G}^{(k)}$,
especially at the beginning of the solution process. Recall that
$\mathcal{G}^{(k)}_{i,j}$ is the sub-derivative of
$|{A}^{(k)}_{i,j}+{D}^{(k)}_{i,j}|$, and although the sub-derivative cannot be
computed, it can be approximated by $sign({A}^{(k)}_{i,j}+{D}^{(k)}_{i,j})$.
Eventually, at the late iterations
$\mathcal{G}^{(k)}\odot\mathcal{M}_{A}^{(k)}$ will converge to
$sign(A^{(k)})$, as all zero elements in $A^{(k)}$ will be out of the free-set
$\mathcal{S}^{(k)}_{A}$. Therefore, their elements in
$\mathcal{G}^{(k)}\odot\mathcal{M}_{A}^{(k)}$ will be zero, and all non-zero
elements in $A^{(k)}$ will have a constant $\mathcal{G}^{(k)}$, which is their
sign. Moreover, at those later iterations, where $sign(A^{(k)})$ is not
expected to change at all, the above formula Eq. 22 solves Eq. 21 completely
as all the terms except $D^{(k)}$ are constant. On the other hand, at the
initial iterations $\mathcal{G}^{(k)}\odot\mathcal{M}_{A}^{(k)}$ cannot be
easily computed because some elements in the free set may change sign or be
zero in the final solution. To this end, we approximate $\mathcal{G}^{(k)}$ by
predicting the “next” sign of the elements. That is, we take the sign of the
non-zero entries, and if an entry is zero, then we take the opposite sign of
its gradient. Thus, a good approximation of $\mathcal{G}^{(k)}$ is:

(25)
$\mathcal{G}^{(k)}_{i,j}=\left\\{\begin{array}[]{lr}sign(A^{(k)}_{i,j})&A^{(k)}_{i,j}\neq
0\\\ -sign(g^{(k)}_{i,j})&A^{(k)}_{i,j}=0\end{array}\right.\,,$

where $g^{(k)}=\nabla f(A^{(k)})$ as defined before. This means that we choose
the next signs according to a small step of proximal gradient descent. The
above $\mathcal{G}^{(k)}$ together with the free-set restriction mask
$\mathcal{M}^{(k)}$ define the minimum sub-gradient of $F(A^{(k)})$ with
respect to $\ell_{2}$-norm.

####  4.3.2 Performing Linesearch to find $t$

At each iteration we perform linesearch to obtain a parameter $t$ satisfying:

(26) $\displaystyle A^{(k)}+\mathcal{M}_{A}^{(k)}\odot D^{(k)}\succ 0\,,$ (27)
$\displaystyle F(A^{(k)}+\mathcal{M}_{A}^{(k)}\odot D^{(k)})<F(A^{(k)})\,.$

##  5 Convergence Analysis

Throughout this section, unless stated otherwise, all the matrices are assumed
to be symmetric and real. The definitions of $f(A)$, $F(A)$, $P(\Delta;t)$,
$g$, $\mathcal{M}$, $G$ and $D$ are as above. We will denote $D_{t}^{*(k)}$ as
the direction $D^{(k)}$ which satisfies equation Eq. 21 for a given $t$:

(28)
$0\in\Bigg{\\{}\frac{1}{t}D_{t}^{*(k)}+A^{(k)}\left(g^{(k)}\odot\mathcal{M}_{A}^{(k)}\right)A^{(k)}+\alpha
A^{(k)}\left(\mathcal{G}^{(k)}\odot\mathcal{M}_{A}^{(k)}\right)A^{(k)}\\\
\Bigg{|}\ \mathcal{G}^{(k)}\in G^{(k)}\left(\mathcal{M}_{A}^{(k)}\odot
D_{t}^{*(k)}\right)\Bigg{\\}}\,.$

Also, denote by $\mathcal{G}_{t}^{*(k)}$ the $\mathcal{G}^{(k)}$ for which
$D_{t}^{*(k)}$ attains 0:

(29)
$0=\frac{1}{t}D_{t}^{*(k)}+A^{(k)}\left(g^{(k)}\odot\mathcal{M}_{A}^{(k)}\right)A^{(k)}+\alpha
A^{(k)}\left(\mathcal{G}_{t}^{*(k)}\odot\mathcal{M}_{A}^{(k)}\right)A^{(k)}\,.$

In the proofs only, we will omit the iteration symbol when the meaning is
eveident.

This section is organized as follows. In Section 5.1 we prove that if our
pISTA iteration does not update $A^{(k)}$, then
$A^{(k)}=\operatorname*{arg\,min}_{A\succ 0}F(A)$. In Section 5.2 we prove
that for every iteration there exists a $t\geq t_{pd}>0$ which satisfies Eq.
26 and keeps the matrix positive definite. Lastly, in Section 5.3 we prove
that for every iteration there exists such a $t\geq t_{dc}>0$ that also
satisfies Eq. 27 and decreases our function value. Combining all the theorems
and lemmas we get the following:

######  Corollary 5.1.

By Lemma 5.2, Theorem 5.5, Corollary 5.7, Lemma 5.15, Lemma 5.19 and Corollary
5.23 \- pISTA converges to $\operatorname*{arg\,min}_{A\succ 0}F(A)$.

###  5.1 Fixed point iteration at minimum

The following Lemma 5.2 states that if a pISTA iteration ends with no update
to $A^{(k)}$, then $A^{(k)}=\operatorname*{arg\,min}_{A\succ 0}F(A)$.

######  Lemma 5.2.

$D_{t}^{*(k)}=0\iff A^{(k)}=\operatorname*{arg\,min}_{A\succ 0}F(A)$.

######  Proof 5.3.

First, we will prove that if $A$ is the minimum of $F(\cdot)$, then pISTA will
result in no update. Assume $A^{(k)}=\operatorname*{arg\,min}_{A\succ 0}F(A)$,
therefore 0 is one of the sub-gradients of $F(A)$ at $A^{(k)}$:

(30) $0\in\frac{\partial F(A^{(k)})}{\partial
A}=\Bigg{\\{}g+\alpha\mathcal{G}\Bigg{|}\begin{array}[]{lr}\mathcal{G}_{i,j}=sign(A_{i,j})&A_{i,j}\neq
0\\\
\mathcal{G}_{i,j}\in[-1,1]&A_{i,j}=0\end{array}\Bigg{\\}}=\Bigg{\\{}g+\alpha\mathcal{G}\
\Bigg{|}\ \mathcal{G}\in G(0)\Bigg{\\}}\,.$

Let $\mathcal{G}^{*}$ be the sub-gradient $\mathcal{G}$ for which
$g+\alpha\mathcal{G}^{*}=0$. Thus:

$\displaystyle g+\alpha\mathcal{G}^{*}=0$ $\displaystyle\Rightarrow$
$\displaystyle(g+\alpha\mathcal{G}^{*})\odot\mathcal{M}_{A}=0$
$\displaystyle\Rightarrow$ $\displaystyle
A\left((g+\alpha\mathcal{G}^{*})\odot\mathcal{M}_{A}\right)A=0$
$\displaystyle\Rightarrow$ $\displaystyle
A\left(g\odot\mathcal{M}_{A}\right)A+\alpha
A\left(\mathcal{G}^{*}\odot\mathcal{M}_{A}\right)A=0\,.$

If we write this in a sub-differential form, we get:

$0\in\Bigg{\\{}A\left(g\odot\mathcal{M}_{A}\right)A+\alpha
A\left(\mathcal{G}\odot\mathcal{M}_{A}\right)A\ \Bigg{|}\ \mathcal{G}\in
G(0)\Bigg{\\}}\,,$

which, according to Eq. 28, means that $D_{t}^{*}=0$ has a 0 in its sub-
gradients.

Now, we will prove that if a pISTA iteration ends with no update, then $A$ is
the $F(\cdot)$ minimum. Assume $D_{t}^{*}=0$, then, according to Eq. 29,
$\mathcal{G}_{t}^{*}$ satisfies:

$\displaystyle A\left(g\odot\mathcal{M}_{A}\right)A+\alpha
A\left(\mathcal{G}_{t}^{*}\odot\mathcal{M}_{A}\right)A=0$
$\displaystyle\Rightarrow$ $\displaystyle
A\left((g+\alpha\mathcal{G}_{t}^{*})\odot\mathcal{M}_{A}\right)A=0\quad/\cdot
A^{-1}$ $\displaystyle\Rightarrow$
$\displaystyle(g+\alpha\mathcal{G}_{t}^{*})\odot\mathcal{M}_{A}=0\,.$

Let us now define $\hat{\mathcal{G}}$ so that
$g+\alpha\hat{\mathcal{G}}\in\frac{\partial F}{\partial A}$:

$\displaystyle\hat{\mathcal{G}}_{i,j}=\Bigg{\\{}\begin{array}[]{lr}[\mathcal{G}^{*}_{t}]_{i,j}&[\mathcal{M}_{A}]_{i,j}\neq
0\\\ \left[g/\alpha\right]_{i,j}&[\mathcal{M}_{A}]_{i,j}=0\end{array}\,.$

It is clear that $g+\alpha\hat{\mathcal{G}}=0$, because if
$[\mathcal{M}_{A}]_{i,j}=0$, then by definition $A_{i,j}=0$ and
$|g|\leq\alpha$. This implies that $A^{(k)}=\operatorname*{arg\,min}_{A\succ
0}F(A)$.

In the proof of Lemma 5.2, no assumptions on $t$ were made, therefore we have
the following:

######  Corollary 5.4.

$\exists t>0:\ D_{t}^{*(k)}=0\ \iff\ \forall t>0:D_{t}^{*(k)}=0$.

###  5.2 Satisfying positive definiteness

The following Theorem 5.5 states that if
$A^{(k)}\neq\operatorname*{arg\,min}_{A\succ 0}F(A)$ then there exists a
$t_{0}$ such that every step-size $t_{0}>t>0$ keeps the matrix positive
definite. Thus, as stated in Corollary 5.7, the linesearch parameter can be
bounded away from zero.

######  Theorem 5.5.

For every iteration $k$, there exists a $t_{0}$ such that for every
$t_{0}>t>0$: $A^{(k)}+\mathcal{M}_{A}^{(k)}\odot D_{t}^{*(k)}\succ 0$.

######  Proof 5.6.

We first define a new constant $\eta$:

(31)
$\eta=\max_{\mathcal{G}_{i,j}\in[-1,1]}\left|\left|A\left(g\odot\mathcal{M}_{A}\right)A+\alpha
A\left(\mathcal{G}\odot\mathcal{M}_{A}\right)A\right|\right|_{F}^{2}\,.$

Using $\eta$ we can upper-bound the Frobenius norm of the descent step
$\mathcal{M}_{A}\odot D_{t}^{*}$:

$\displaystyle||\mathcal{M}_{A}\odot D_{t}^{*}||_{F}^{2}$ $\displaystyle\leq$
$\displaystyle||D_{t}^{*}||_{F}^{2}$ $\displaystyle=$ $\displaystyle
t^{2}\left|\left|A\left(g\odot\mathcal{M}_{A}\right)A+\alpha
A\left(\mathcal{G}_{t}^{*}\odot\mathcal{M}_{A}\right)A\right|\right|_{F}^{2}\leq
t^{2}\eta\,.$

Let $t_{0}=\frac{\lambda_{min}(A)}{\sqrt{\eta}}$, then, for $0<t<t_{0}$:

$\displaystyle|\lambda_{min}(\mathcal{M}_{A}\odot
D_{t}^{*})|\leq||\mathcal{M}_{A}\odot D_{t}^{*}||_{F}\leq
t\sqrt{\eta}<\lambda_{min}(A)\,,$

and by Weyl’s inequality:

$\displaystyle\lambda_{min}(A+\mathcal{M}_{A}\odot
D_{t}^{*})\geq\lambda_{min}(A)+\lambda_{min}(\mathcal{M}_{A}\odot
D_{t}^{*})>0$ $\displaystyle\Rightarrow A+\mathcal{M}_{A}\odot D_{t}^{*}\succ
0\,.$

######  Corollary 5.7.

The linesearch parameter $t$ can be chosen to be bounded away from zero while
keeping the matrix positive definite, i.e., $t\geq\frac{t_{0}}{2}>0$.

###  5.3 Decrease in function objective

In this section we show that under suitable conditions the sequence
$\\{F(A^{(k)})\\}$ converges to $\min_{A\succ 0}F(A)$ which means that the
sequence $\\{A^{(k)}\\}$ converges to $\operatorname*{arg\,min}_{A\succ
0}F(A)$. Then, we prove that our pISTA iteration satisfies those conditions,
which completes the convergence proof. In this section we denote the space of
symmetric positive definite matrices by $S_{++}$.

Denote:

$\begin{split}T(\Delta;t)=P(\mathcal{M}_{A}^{(k)}\odot\Delta;t)&=f(A^{(k)})\\\
&\quad\quad+\left\langle
g^{(k)},\mathcal{M}_{A}^{(k)}\odot\left(A^{(k)}(\mathcal{M}_{A}^{(k)}\odot\Delta)A^{(k)}\right)\right\rangle\\\
&\quad\quad+\frac{1}{2t}\left\langle
A^{(k)}(\mathcal{M}_{A}^{(k)}\odot\Delta)A^{(k)},(\mathcal{M}_{A}^{(k)}\odot\Delta)\right\rangle\\\
&\quad\quad+\alpha\left|\left|A^{(k)}+\mathcal{M}_{A}^{(k)}\odot\left(A^{(k)}(\mathcal{M}_{A}^{(k)}\odot\Delta)A^{(k)}\right)\right|\right|_{1}\,.\end{split}$

We state a relation between $\operatorname*{arg\,min}_{\Delta}T(\Delta;t)$ and
the descent direction $\mathcal{M}_{A}^{(k)}\odot D_{t}^{*(k)}$:

######  Lemma 5.8.

Denote
$\Delta_{t}^{*(k)}\overset{\Delta}{=}\operatorname*{arg\,min}_{\Delta}T(\Delta;t)$,
then

$A^{(k)}\left(\mathcal{M}_{A}^{(k)}\odot\Delta_{t}^{*(k)}\right)A^{(k)}=\mathcal{M}_{A}^{(k)}\odot
D_{t}^{*(k)}\,.$

######  Proof 5.9.

The sub-differential of $T(\Delta;t)$ is defined as:

$\frac{\partial
T(\Delta;t)}{\partial\Delta}=\Bigg{\\{}\mathcal{M}_{A}\odot\left(\frac{1}{t}A\left(\mathcal{M}_{A}\odot\Delta\right)A+A\left(g\odot\mathcal{M}_{A}\right)A+\alpha
A\left(\mathcal{G}\odot\mathcal{M}_{A}\right)A\right)\\\ \Bigg{|}\
\mathcal{G}\in
G\left(\mathcal{M}_{A}\odot\big{(}A(\mathcal{M}_{A}\odot\Delta)A\big{)}\right)\Bigg{\\}}\,.$

Denote $Z=A(\mathcal{M}_{A}\odot\Delta)A$, then we get:

$\frac{\partial
T(\Delta;t)}{\partial\Delta}=\Bigg{\\{}\mathcal{M}_{A}\odot\left(\frac{1}{t}Z+A\left(g\odot\mathcal{M}_{A}\right)A+\alpha
A\left(\mathcal{G}\odot\mathcal{M}_{A}\right)A\right)\\\ \Bigg{|}\
\mathcal{G}\in G\left(\mathcal{M}_{A}\odot Z\right)\Bigg{\\}}\,.$

We notice that, for $Z_{t}^{*}=\mathcal{M}_{A}\odot D_{t}^{*}$, 0 is in the
sub-differential:

$0\in\Bigg{\\{}\mathcal{M}_{A}\odot\left(\frac{1}{t}Z_{t}^{*}+A\left(g\odot\mathcal{M}_{A}\right)A+\alpha
A\left(\mathcal{G}\odot\mathcal{M}_{A}\right)A\right)\\\ \Bigg{|}\
\mathcal{G}\in G\left(\mathcal{M}_{A}\odot Z_{t}^{*(k)}\right)\Bigg{\\}}\,.$

Thus, if ${\Delta}_{t}^{*}=\operatorname*{arg\,min}_{\Delta}T(\Delta;t)$, then
$A\left(\mathcal{M}_{A}\odot\Delta_{t}^{*}\right)A=\mathcal{M}_{A}\odot
D_{t}^{*}$.

We next prove that $f(A)$ has a Lipschitz continuous gradient under suitable
conditions and state an important property for this kind of functions.

######  Lemma 5.10.

For any square matrices $A\succ 0$ and $B$:
$||AB||_{F}\leq||A||_{2}||B||_{F}$.

######  Proof 5.11.

Denote by $b_{i}$ the $i$-th column of matrix B, then:

$||AB||_{F}^{2}=\sum_{i}||Ab_{i}||_{2}^{2}\leq||A||_{2}^{2}\sum_{i}||b_{i}||_{2}^{2}=||A||_{2}^{2}||B||_{F}^{2}\,.$

######  Lemma 5.12.

Let $\beta$ be a positive constant, then the function $f(A)$ has a Lipschitz
continuous gradient over the domain $B=\\{\beta I\preceq X\preceq\gamma
I\big{|}X\in S_{++}\\}$ with a Lipschitz constant of $\frac{1}{\beta^{2}}$.

######  Proof 5.13.

To prove that the function $f(A)$ has a Lipschitz continuous gradient with a
Lipschitz constant of $\frac{1}{\beta^{2}}$, we need to prove that $||\nabla
f(X)-\nabla f(Y)||_{F}\leq\frac{1}{\beta^{2}}||X-Y||_{F}$ for every $X,Y\in
B$. The gradient is given by $\nabla f(A)=S-A^{-1}$. Thus, for every $X,Y\in
B$:

$\displaystyle||\nabla f(X)-\nabla f(Y)||_{F}$
$\displaystyle=||Y^{-1}-X^{-1}||_{F}$
$\displaystyle=||Y^{-1}(X-Y)X^{-1}||_{F}$
$\displaystyle\leq||Y^{-1}||_{2}||X-Y||_{F}||X^{-1}||_{2}$
$\displaystyle=\frac{1}{\lambda_{min}(Y)}\frac{1}{\lambda_{min}(X)}||X-Y||_{F}$
$\displaystyle\leq\frac{1}{\beta^{2}}||X-Y||_{F}\,.$

######  Corollary 5.14.

[Properties of a function with Lipschitz continuous gradient] Let $\beta$ be a
positive constant and $B=\\{\beta I\preceq X\preceq\gamma I\big{|}X\in
S_{++}\\}$, then for every $Y,X\in B$:

$f(Y)\leq f(X)+\left\langle\nabla
f(X),(Y-X)\right\rangle+\frac{1}{2\beta^{2}}\left\langle
Y-X,Y-X\right\rangle\,.$

The rest of the proofs are inspired by [27], [16] and [29], following the
proofs of [27] with slight modifications.

######  Lemma 5.15 (corresponds to lemma 4.1 in [27]).

Let $\beta$ be a positive constant and $B=\\{\beta I\preceq X\preceq\gamma
I\big{|}X\in S_{++}\\}$. Assume that for every iteration $k$, $A^{(k)}\in B$,
then:

(32) $F(A^{(k)})-F(A^{(k+1)})\geq L\cdot||A^{(k)}-A^{(k+1)}||^{2}_{F}\,.$

Where $L>0$. Furthermore, the linesearch parameter $t$ can be chosen to be
bounded away from zero, i.e., $t\geq t_{min}>0$.

######  Proof 5.16.

According to pISTA iteration, $A^{(k+1)}=A^{(k)}+\mathcal{M}_{A}^{(k)}\odot
D_{t}^{*(k)}$. Therefore,

$\begin{split}F(A^{(k)})-F(A^{(k+1)})&=F(A^{(k)})-F(A^{(k)}+\mathcal{M}_{A}^{(k)}\odot
D_{t}^{*(k)})\end{split}\,,$

and from now on, we shall omit the iteration symbol $k$:

$\begin{split}F(A^{(k)})-F(A^{(k+1)})&=F(A)-F(A+\mathcal{M}_{A}\odot
D_{t}^{*})\\\ &=f(A)+\alpha||A||_{1}-f(A+\mathcal{M}_{A}\odot
D_{t}^{*})-\alpha||A+\mathcal{M}_{A}\odot D_{t}^{*}||_{1}\,.\end{split}$

Using Corollary 5.14 on $f(A+\mathcal{M}_{A}\odot D_{t}^{*})$:

$\begin{split}F(A^{(k)})-F(A^{(k+1)})&\geq\alpha||A||_{1}-\Big{(}\left\langle
g,\mathcal{M}_{A}\odot D_{t}^{*}\right\rangle\\\
&\quad+\frac{1}{2\beta^{2}}\left\langle\mathcal{M}_{A}\odot
D_{t}^{*},\mathcal{M}_{A}\odot
D_{t}^{*}\right\rangle+\alpha||A+\mathcal{M}_{A}\odot
D_{t}^{*}||_{1}\Big{)}\,.\end{split}$

By the definition of $T(\Delta;t)$ and Lemma 5.8:

$\begin{split}\left\langle g,\mathcal{M}_{A}\odot
D_{t}^{*}\right\rangle=T(\Delta^{*}_{t};t)-\frac{1}{2t}\left\langle\mathcal{M}_{A}\odot
D_{t}^{*},A^{-1}\left(\mathcal{M}_{A}\odot
D_{t}^{*}\right)A^{-1}\right\rangle-f(A)\,,\end{split}$

therefore,

$\begin{split}F(A^{(k)})-F(A^{(k+1)})&\geq
f(A)+\alpha||A||_{1}-\frac{1}{2\beta^{2}}\left\langle\mathcal{M}_{A}\odot
D_{t}^{*},\mathcal{M}_{A}\odot D_{t}^{*}\right\rangle\\\
&\quad+\frac{1}{2t}\left\langle\mathcal{M}_{A}\odot
D_{t}^{*},A^{-1}\left(\mathcal{M}_{A}\odot
D_{t}^{*}\right)A^{-1}\right\rangle-T(\Delta_{t}^{*};t)\,.\end{split}$

Because ${\Delta}_{t}^{*}=\operatorname*{arg\,min}_{\Delta}T(\Delta;t)$, we
know that $T(\Delta_{t}^{*};t)\leq T(0;t)=f(A)+\alpha||A||_{1}$:

$\begin{split}F(A^{(k)})-F(A^{(k+1)})&\geq
f(A)+\alpha||A||_{1}-\frac{1}{2\beta^{2}}\left\langle\mathcal{M}_{A}\odot
D_{t}^{*},\mathcal{M}_{A}\odot D_{t}^{*}\right\rangle\\\
&\quad+\frac{1}{2t}\left\langle\mathcal{M}_{A}\odot
D_{t}^{*},A^{-1}\left(\mathcal{M}_{A}\odot
D_{t}^{*}\right)A^{-1}\right\rangle-T(0;t)\\\
&\geq-\frac{1}{2\beta^{2}}\left\langle\mathcal{M}_{A}\odot
D_{t}^{*},\mathcal{M}_{A}\odot D_{t}^{*}\right\rangle\\\
&\quad+\frac{1}{2t}\left\langle\mathcal{M}_{A}\odot
D_{t}^{*},A^{-1}\left(\mathcal{M}_{A}\odot
D_{t}^{*}\right)A^{-1}\right\rangle\\\
&\geq(\frac{1}{2t\gamma^{2}}-\frac{1}{2\beta^{2}})\left\langle\mathcal{M}_{A}\odot
D_{t}^{*},\mathcal{M}_{A}\odot D_{t}^{*}\right\rangle\,.\end{split}$

Denote $L=\frac{1}{2t\gamma^{2}}-\frac{1}{2\beta^{2}}$, then $L$ is positive
as long as $t<(\frac{\beta}{\gamma})^{2}$. Furthermore, $t$ can be chosen to
be bounded away from zero, i.e., $t\geq
t_{min}=\frac{1}{2}\cdot(\frac{\beta}{\gamma})^{2}>0$.

######  Lemma 5.17 (corresponds to lemma 4.2 on [27]).

Let $\beta$ be a positive constant and $B=\\{\beta I\preceq X\preceq\gamma
I\big{|}X\in S_{++}\\}$, and assume that for every iteration $k$, $A^{(k)}\in
B$. Let $\\{A^{(k_{j})}\\}$ be any infinite and converging sub-series of
$\\{A^{(k)}\\}$ and let $\bar{A}$ denote its limit. Then
$\bar{A}=\operatorname*{arg\,min}_{A\succ 0}F(A)$.

######  Proof 5.18.

Since the sub-series $\\{A^{(k_{j})}\\}$ convergences to $\bar{A}$, then
$\\{F(A^{(k_{j})})\\}$ convergences to $F(\bar{A})$. According to Lemma 5.15,
the full series $\\{F(A^{(k)})\\}$ is monotone and hence convergent because
$\\{F(A^{(k_{j})})\\}$ is convergent. Therefore,
$\\{F(A^{(k_{j}+1)})-F(A^{(k_{j})})\\}\rightarrow 0$, which implies following
Eq. 32 that $||A^{(k_{j})}-A^{(k_{j}+1)}||^{2}_{F}\rightarrow 0$ and
$\lim_{j\rightarrow\infty}A^{(k_{j}+1)}=\bar{A}$. According to pISTA,
$A^{(k+1)}=A^{(k)}+\mathcal{M}_{A}^{(k)}\odot D_{t}^{*(k)}$ and because
$\mathcal{M}_{A}^{(k)}\neq 0$,
$||A^{(k_{j})}-A^{(k_{j}+1)}||^{2}_{F}\rightarrow 0$ implies
$||D_{t}^{*(k)}||^{2}_{F}\rightarrow 0$. We know that $D_{t}^{*(k)}$ satisfies
Eq. 28, by taking the limit as $j\rightarrow\infty$ and using Lemma 5.2 we get
that $\bar{A}=\operatorname*{arg\,min}_{A\succ 0}F(A)$.

######  Lemma 5.19.

[corresponds to lemma 4.3 on [27]] Let $\beta$ be a positive constant and
$B=\\{\beta I\preceq X\preceq\gamma I\big{|}X\in S_{++}\\}$. Assume that for
every iteration $k$, $A^{(k)}\in B$. Then the series $\\{A^{(k)}\\}$ has a
limit and it is given by $\operatorname*{arg\,min}_{A\succ 0}F(A)$.

######  Proof 5.20.

The domain $B$ is compact, the rest of the proof follows the proof of lemma
4.3 in [27] exactly.

We have shown that pISTA converges under the condition of bounded eigenvalues.
It remains to show that pISTA satisfies it.

######  Lemma 5.21 (corresponds to lemma 2 on [16]).

Let $U=\\{A|F(A)\leq F(A_{0})\text{ and }A\in S_{++}\\}$. Then, $U\subseteq
B=\\{\beta I\preceq X\preceq\gamma I\big{|}X\in S_{++}\\}$ where $\beta>0$.

######  Proof 5.22.

See lemma 2 in [16].

######  Corollary 5.23.

The series $\\{A^{(k)}\\}$ created by the pISTA algorithm satisfies Eq. 27,
hence, $A$ is contained in $U$ and in $B$.

##  6 Numerical Results

We compare the performance of our pISTA algorithm to the performance of GISTA
[14], Orthant-Based Newton (OBN) [22] and ALM [24] (without skipping step111As
implemented by the authors: https://www.math.ucdavis.edu/~sqma/ALM-SICS.html),
all of which can be implemented efficiently and easily on GPU based on matrix
operations. We note that Newton-Lasso [22], VSM [21] and PSM [9] can be
implemented in the same manner on GPU too, but all of them have worse
performance than OBN [22] and ALM [24] according to the authors’ measurements.
Other algorithms, such as BIG&QUIC [17] and BCD-IC [26], can be implemented on
GPU as well. However, their implementations might be highly complex due to the
differences between the CPU parallelism model and the GPU parallelism model,
or their implementations might have worse results than the CPU counterparts
due to the numerous scalar operations. The algorithms are all implemented in
Python as recommended by their authors using $32bit$ floating point precision.
Each algorithm is implemented using the GPU as efficiently as possible with
the Cupy library.

We initialize all methods using GISTA [14] initialization where:

(33) $A^{(0)}_{i,j}=\bigg{\\{}\begin{array}[]{lr}(S_{i,i}+\alpha)^{-1}&i=j\\\
0&i\neq j\end{array}\,.$

For ALM only, we consider any number which is less than or equal to
$10\cdot\epsilon_{machine}$ as zero and set those numbers to zero at end of
each iteration. For OBN we use 10 inner iterations, and for ALM we use the
hyperparamters used by [24] except that we set $\mu_{0}=\frac{1}{\alpha}$ if
$\alpha<0.5$. For our pISTA algorithm, we limit our linesearch over $t$ by
stopping if $t$ is less than $10^{-4}$. In that case, we use a step of
$t\leq\left(\frac{0.9}{cond(A^{(k)})}\right)^{2}$ which keeps the matrix
positive definite. As proved in Section 5.3, we can use a
$t<(\frac{\beta}{\gamma})^{2}$ where $\beta$ and $\gamma$ are lower and upper
bounds of the eigenvalues, respectively. However, we cannot compute those
bounds, so we assume that $\left(\frac{0.9}{cond(A^{(k)})}\right)^{2}$ is
sufficient and it should satisfy the linesearch criteria.

As a stopping criterion for all the method, we follow [27, 17, 26] and use
$\min_{z}||\partial F(A^{(k)}||_{1}<10^{-2}||A^{(k)}||_{1}$, where
$\min_{z}||\partial F(A^{(k)}||$ is the minimum sub-gradient norm.

All experiments were run on a machine with Intel(R) Xeon(R) Gold 6230 2.10GHz
processor with 4 cores, 20GB RAM, GeForce RTX 2080 Ti GPU and Ubuntu 18.04.5
operating system. We also use Python3.9.2, CUDA11.0, Numpy1.20.1 and Cupy9.0.
Our full python framework and code can be found at

https://github.com/GalSha/GLASSO_Framework.

###  6.1 Synthetic Experiments

First, we evaluate the algorithms on synthetic data. We use three different
types of matrices as our ground truth:

  * •

Chain graphs: as described in [16]. The ground truth matrix $\Sigma^{-1}$ is
set to be $\Sigma^{-1}_{i,i}=1$ and
$\Sigma^{-1}_{i,i+1}=\Sigma^{-1}_{i,i-1}=-0.5$.

  * •

Graphs with random sparsity structures: as described in [20]. We generate a
sparse matrix $U$ with non-zero elements set to be $+1$ or $-1$. We define the
ground truth as $\Sigma^{-1}=U^{T}U$ where all off diagonal elements are in
$[-1,1]$. We control the number of non-zeros in $U$ and tune it such that
$\Sigma^{-1}$ has approximately $0.5\%$ non-zeros.

  * •

Planar graphs: as described in [27]. First we create our graph $G(V,E)$. We
generate $n$ random points on the unit square to be our vertices $V$. We use
Delaunay triangulation to generate our edges $E$. Given the graph
connectivity, we define $\Sigma^{-1}$ as its graph Laplacian, i.e.,
$\Sigma^{-1}_{i,j}=-1$ if $(i,j)\in E$ and $\Sigma^{-1}_{i,i}=d_{i}$ where
$d_{i}$ is the degree of vertex $i$.

To ensure that the matrices are positive definite, we add a predefined
diagonal term of $\max\\{-1.2\lambda_{min},10^{-1}\\}\cdot I$.

We do two sets of experiments, with $n=1,000$ and with $n=10,000$. In each set
we draw $3\%\cdot n$ of samples and run the algorithms with two different
values of $\alpha$. We repeat each experiment five times and show the average
results in Table 1 and Table 2 for the first set and second set, respectively.
Also, we compare the results of a CPU only implementation of our pISTA
algorithm for both sets of experiments.

In the tables, we show the average time and iterations it took for each
algorithm to reach the stopping criterion. We also show the number of non-
zeros in the output matrix, and the minimum sub-gradient Frobenius norm.

Problem Parameters | pISTA (GPU) | pISTA (CPU) | GISTA | OBN | ALM
---|---|---|---|---|---
n $m/n$ $\alpha$ |  $\Sigma^{-1}$ type $|\Sigma^{-1}|_{0}$ |  time (iter) $\min_{z}||\partial F(\tilde{\Sigma}^{-1})||_{F}$ $||\tilde{\Sigma}^{-1}||_{0}$
1,000 3% 0.6 |  $Chain$ 2998 |  0.04s (2.0) 0.0017 2959.2 |  0.35s (2.0) 0.0017 2959.2 |  0.05s (3.0) 0.0275 2970.4 |  0.05s (2.0) 0.0020 2959.2 |  0.55s (11.0) 0.0312 2966.0
1,000 3% 0.4 |  $Chain$ 2998 |  0.10s (6.6) 0.0142 25307.2 |  1.40s (6.6) 0.0142 25307.2 |  0.15s (10.6) 0.0336 25302.0 |  0.11s (4.0) 0.0003 25246.0 |  2.47s (53.8) 0.5966 25246.8
1,000 3% 0.6 |  $Random$ 5936 |  0.04s (2.2) 0.0265 2184.0 |  0.32s (2.2) 0.0265 2184.0 |  0.06s (3.6) 0.0462 2191.2 |  0.06s (2.2) 0.0824 2185.2 |  0.61s (14.0) 0.8561 2209.0
1,000 3% 0.4 |  $Random$ 5936 |  0.10s (6.4) 0.0666 26335.2 |  1.05s (6.4) 0.0666 26335.2 |  0.21s (15.2) 0.0286 26400.4 |  0.13s (4.6) 0.0005 26246.4 |  2.41s (53.4) 2.5377 26256.8
1,000 3% 0.6 |  $Planar$ 6958 |  0.03s (2.0) 0.0008 2995.6 |  0.31s (2.0) 0.0008 2995.6 |  0.04s (2.4) 0.0269 3011.6 |  0.05s (2.0) 0.0028 2997.6 |  0.47s (10.2) 0.3156 3039.6
1,000 3% 0.4 |  $Planar$ 6958 |  0.24s (15.4) 0.0244 28495.2 |  3.05s (15.6) 0.0244 28471.6 |  0.28s (20.2) 0.0241 28623.2 |  0.15s (5.6) 0.0130 28290.4 |  3.13s (69.6) 1.8171 28312.8

Table 1: Results for $1,000\times 1,000$ Precision Matrix
Bold text marks the best results

Problem Parameters | pISTA (GPU) | pISTA (CPU) | GISTA | OBN | ALM
---|---|---|---|---|---
n $m/n$ $\alpha$ |  $\Sigma^{-1}$ type $|\Sigma^{-1}|_{0}$ |  time (iter) $\min_{z}||\partial F(\tilde{\Sigma}^{-1})||_{F}$ $||\tilde{\Sigma}^{-1}||_{0}$
10,000 3% 0.4 |  $Chain$ 29,998 |  5.69s (3.0) 0.0101 34685.2 |  249.18s (3.0) 0.0101 34685.2 |  10.80s (6.0) 0.1478 35104.0 |  13.85s (3.0) 0.0060 34839.2 |  103.44s (28.6) 28.5176 34879.6
10,000 3% 0.2 |  $Chain$ 29,998 |  12.89s (7.0) 0.0951 97765.6 |  578.37s (7.0) 0.0951 97765.6 |  27.53s (16.6) 0.8418 98432.8 |  20.02s (4.0) 0.0129 97881.2 |  168.46s (48.0) 24.0052 98339.2
10,000 3% 0.4 |  $Random$ 508,787 |  1.56s (0.6) 2.1558 10823.6 |  68.74s (0.6) 2.1558 10823.6 |  1.85s (0.6) 3.0901 10823.6 |  1.60s (0.6) 2.1558 10823.6 |  767.79s (237.4) 44.2140 11072.0
10,000 3% 0.2 |  $Random$ 508,787 |  5.67s (3.0) 4.3416 140699.6 |  263.86s (3.0) 4.3416 140699.6 |  15.95s (9.0) 0.8136 140869.6 |  13.97s (3.0) 1.0047 140715.6 |  153.16s (44.2) 21.1711 141009.8
10,000 3% 0.4 |  $Planar$ 69,949 |  3.92s (2.0) 1.2172 63249.6 |  174.70s (2.0) 1.2172 63249.6 |  8.84s (5.0) 0.0477 63446.8 |  13.73s (3.0) 0.1113 63324.0 |  100.34s (28.2) 6.6973 63572.8
10,000 3% 0.2 |  $Planar$ 69,949 |  20.91s (10.2) 1.4050 194261.6 |  943.43s (10.2) 1.4049 194261.2 |  63.17s (35.8) 0.5368 195714.0 |  41.54s (7.6) 1.3284 192740.0 |  269.97s (79.4) 12.4953 193075.0

Table 2: Results for $10,000\times 10,000$ Precision Matrix
Bold text marks the best results

From Table 1, we see that pISTA outperforms every other algorithm with respect
to running time on all the low-dimension matrices except one. We notice that
pISTA requires less iterations than GISTA and more iterations than OBN. This
is expected as GISTA is a first order method, OBN is a second order method,
and the pISTA algorithm is a quasi-Newton method. Thus, pISTA convergence rate
should be between GISTA and OBN. For small $\alpha$, OBN achieves a better
minimum sub-gradient and better sparsity, however, the difference is
negligible. Also, the GPU implementation of pISTA is up to $\times 14$ times
faster than its CPU counterpart, proving the desirability of GPU centric
algorithms.

In Table 2, we see similar results. pISTA outperforms every other algorithm
with respect to running time on all matrices. We notice that pISTA also
achieves better sparsity patterns than the other algorithms. As expected,
pISTA requires less iterations than GISTA and more iterations than OBN. For
high-dimension matrices, the speedup of the GPU over CPU is up to $\times 44$,
making it extremely valuable in high dimensions.

$0$$5$$10$$15$$20$$25$$30$$35$$40$$45$$50$$10^{-4}$$10^{-3}$$10^{-2}$$10^{-1}$$10^{0}$$10^{1}$$10^{2}$$10^{3}$$10^{4}$Iteration$F(\tilde{\Sigma}^{-1})-{min}F(A)$pISTAGISTAOBNALM
(a) Chain Precision Matrix

$0$$5$$10$$15$$20$$25$$30$$35$$40$$45$$10^{-4}$$10^{-3}$$10^{-2}$$10^{-1}$$10^{0}$$10^{1}$$10^{2}$$10^{3}$Iteration$F(\tilde{\Sigma}^{-1})-{min}F(A)$pISTAGISTAOBNALM
(b) Random Precision Matrix

$0$$10$$20$$30$$40$$50$$60$$70$$80$$10^{-4}$$10^{-3}$$10^{-2}$$10^{-1}$$10^{0}$$10^{1}$$10^{2}$$10^{3}$$10^{4}$Iteration$F(\tilde{\Sigma}^{-1})-{min}F(A)$pISTAGISTAOBNALM
(c) Planar Precision Matrix

Figure 1: Semilog plot of $F(\tilde{\Sigma}^{-1})-{min}F(A)$ as function of
iteration333The iterations are not equal in complexity or time between the
algorithms for $10,000\times 10,000$ Precision Matrix with $\alpha=0.2$

In Fig. 1, we show a semi-log plot of $F(\tilde{\Sigma}^{-1})-{min}F(A)$ as a
function of the iteration, however, it is important to note that the
iterations are not equal in complexity or required time. We only present the
plots of the $10,000\times 10,000$ Precision Matrices for $\alpha=0.2$ because
of space consideration. We define ${min}F(A)$ as the minimum value achieved
for $F$ among all the algorithms and iterations in that experiment. The plots
show what we expect, pISTA, as a quasi-Newton method, achieves a convergence
rate between linear (GISTA) and quadratic (OBN). Moreover, we see that the
convergence rate is quadratic in the first few iterations.

###  6.2 Real World Data Experiments

For real world data we use gene expression data sets that are available at the
Gene Expression Omnibus http://www.ncbi.nlm.nih.gov/geo/. We preprocess the
data to have zero mean and unit variance for each variable, i.e., $diag(S)=I$.
Table 3 shows the results for data sets of various sizes including the name
codes of the data sets used.

Problem Parameters | pISTA | GISTA | OBN | ALM
---|---|---|---|---
Data set |  n $m/n$ $\alpha$ |  time (iter) $\min_{z}||\partial F(\tilde{\Sigma}^{-1})||_{F}$ $||\tilde{\Sigma}^{-1}||_{0}$
GSE-3016 |  1322 4.77% (63) 0.85 |  0.09s (3) 0.1915 3478 |  0.11s (3) 0.0378 3560 |  0.16s (3) 0.2185 3490 |  1.01s (18) 0.2921 3596
GSE-3016 |  1322 4.77% (63) 0.75 |  0.14s (5) 0.0092 10796 |  0.35s (14) 0.0049 10890 |  0.32s (6) 0.0025 10806 |  2.06s (39) 8.5051 10809
GSE-3016 |  1322 4.77% (63) 0.65 |  1.02s (33) 0.0134 20682 |  1.55s (66) 0.0158 20660 |  0.66s (11) 0.0016 20386 |  3.17s (59) 0.0011 20440
GSE-26242 |  1536 6.25% (96) 0.85 |  0.08s (2) 0.0118 4876 |  0.15s (4) 0.0186 4962 |  0.14s (2) 0.0103 4862 |  1.67s (26) 0.0081 4912
GSE-26242 |  1536 6.25% (96) 0.75 |  0.45s (14) 0.0003 13896 |  0.70s (24) 0.0056 13900 |  0.47s (7) 0.0001 13798 |  2.39s (40) 3.5156 13810
GSE-26242 |  1536 6.25% (96) 0.65 |  1.81s (32) 0.0096 28494 |  2.43s (83) 0.0148 28942 |  1.01s (13) 0.0087 28148 |  4.28s (66) 6.7166 28182
GSE-7039 |  6264 3.99% (250) 0.9 |  2.31s (4) 0.0397 34286 |  8.52s (15) 0.0798 34364 |  5.92s (4) 0.0014 34250 |  74.19s (77) 33.9381 34266
GSE-7039 |  6264 3.99% (250) 0.8 |  10.42s (14) 0.0566 49200 |  39.91s (69) 0.1218 50020 |  9.39s (6) 0.0449 49054 |  *s ($>$1000)
GSE-7039 |  6264 3.99% (250) 0.7 |  41.87s (41) 0.0206 70700 |  101.95s (166) 0.1302 70820 |  16.39s (10) 0.0012 70572 |  *s ($>$1000)
GSE-52076 |  11064 5.39% (596) 0.9 |  9.78s (4) 0.0027 32840 |  30.65s (13) 0.1908 32850 |  Out of GPU memory |  103.10s (23) 23.4233 32990
GSE-52076 |  11064 5.39% (596) 0.8 |  10.05s (4) 0.1162 42852 |  55.54s (23) 0.2827 43086 |  Out of GPU memory |  168.27s (39) 71.0796 43471
GSE-52076 |  11064 5.39% (596) 0.7 |  17.34s (7) 0.3950 49612 |  167.93s (65) 0.3866 49644 |  Out of GPU memory |  317.99s (75) 58.4818 49840

Table 3: Results for real world data set
Bold text marks the best results

In Table 3, we see that there is no silver bullet. There are some problems
where pISTA will outshine and there are some problems where OBN will outshine.
However, OBN requires much more memory than pISTA, making it hard to use in
high dimension on GPUs.

##  7 Conclusions

In this work we presented a quasi-Newton method for solving sparse inverse
covariance estimation problem. Our method creates simplified approximate
second order optimization where the different variables dependency can be
guessed in a simple manner. Moreover, the method is designed to be implemented
using matrix operations which can be done on a GPU, allowing us to solve the
problem efficiently. We showed the desirability of GPU implementations and
that our pISTA method has better results for various problems structures.

###  .1 General solution for sign derivative equation

Define $\mathcal{T}(x)$:

$\displaystyle x\neq 0:\quad\quad\mathcal{T}(x)=sign(x)$ $\displaystyle
x=0:\quad\quad\mathcal{T}(x)\in[-1,1]\,.$

######  Lemma .1.

The solution for:

(34) $x+b+c\mathcal{T}(x+a)=0,\quad c>0$

is given by:

(35) $x=\left\\{\begin{array}[]{lr}-b+c&c<b-a\\\ -a&c>|b-a|\\\
-b-c&-c>b-a\end{array}\right.=-a+\text{SoftThreshold}(a-b,c)\,.$

######  Proof .2.

First, consider the case where $x>-a\Rightarrow\mathcal{T}(x+a)=1$:

$\displaystyle\left\\{\begin{array}[]{lr}x+b+c=0\\\ x>-a\\\
\end{array}\right.\Rightarrow\left\\{\begin{array}[]{lr}x=-b-c\\\ x>-a\\\
\end{array}\right.\Rightarrow\left\\{\begin{array}[]{lr}x=-b-c\\\ -b-c>-a\\\
\end{array}\right.\Rightarrow\left\\{\begin{array}[]{lr}x=-b-c\\\ -c>b-a\\\
\end{array}\right.\,.$

In a similar way, for $x<-a\Rightarrow\mathcal{T}(x+a)=-1$ and we get:

$\displaystyle\left\\{\begin{array}[]{lr}x=-b+c\\\ c<b-a\\\
\end{array}\right.\,.$

Lastly, consider the case that $x=-a$, then since $c>0$ we get:

$0\in-a+b+c\cdot t,\;\;t\in[-1,1]\\\ $

which leads to

$\displaystyle a-b\in[-c,c]\Rightarrow c\leq|a-b|\Rightarrow c\leq|b-a|\,.$

###  .2 Approximate solution for the pISTA sub-gradient equation

Consider the following problem:

(36) $D:\quad\quad
0\in\Bigg{\\{}\frac{1}{t}D+A\left(g\odot\mathcal{M}_{A}\right)A+\alpha
A\left(\mathcal{G}\odot\mathcal{M}_{A}\right)A\\\
\Bigg{|}\begin{array}[]{lr}\mathcal{G}_{i,j}=[sign(A+\mathcal{M}_{A}\odot
D)]_{i,j}&[A+\mathcal{M}_{A}\odot D]_{i,j}\neq 0\\\
\mathcal{G}_{i,j}\in[-1,1]&[A+\mathcal{M}_{A}\odot
D]_{i,j}=0\end{array}\Bigg{\\}}\,.$

Let us split the equation into $n^{2}$ scalar equations where the $(i,j)$-th
equation is composed of the elements with indices $(i,j)$. Denote the equation
variable $D_{i,j}=x$. The $(i,j)$ equation is:

(37)
$\frac{1}{t}x+\left[A\left(g\odot\mathcal{M}_{A}\right)A\right]_{i,j}+\alpha\left[A\left(\mathcal{G}\odot\mathcal{M}_{A}\right)A\right]_{i,j}=0\,,$

where $\mathcal{G}_{i,j}=\mathcal{T}\left(A_{i,j}+[\mathcal{M}_{A}]_{i,j}\cdot
x\right)$. As we consider equations where $[\mathcal{M}_{A}]_{i,j}\neq 0$ (the
rest are zeros), we can use
$\mathcal{G}_{i,j}=\mathcal{T}\left(A_{i,j}+x\right)$. Define:

(38) $\mathcal{G}^{-i,-j}_{k,l}=\left\\{\begin{array}[]{lr}0&k=i,l=j\\\
0&k=j,l=i\\\ \mathcal{G}_{k,l}&otherwise\end{array}\right.\,.$

We can write Eq. 37 as:

(39)
$\frac{1}{t}x+\left[A\left(g\odot\mathcal{M}_{A}\right)A\right]_{i,j}+\alpha\left[A\left(\mathcal{G}^{-i,-j}\odot\mathcal{M}_{A}\right)A\right]_{i,j}+C_{i,j}\mathcal{T}\left(A_{i,j}+x\right)=0\,,$

where

(40) $C_{i,j}=\left\\{\begin{array}[]{lc}\alpha\cdot\left(A_{i,i}\cdot
A_{j,j}\right),&i=j\\\ \alpha\cdot\left(A_{i,i}\cdot A_{j,j}+A_{i,j}\cdot
A_{j,i}\right),&i\neq j\end{array}\right.\,$

is the diagonal entry of the Kronecker matrix $\alpha\cdot A\otimes A$
corresponding to the entry $(i,j)$. Note that $C_{i,j}>0$ since $\alpha>0$ and
$A^{(k)}$ is symmetric positive definite, and its diagonal is strictly
positive. According to Lemma .1, the solution to Eq. 39 is:

$\displaystyle x=-A_{i,j}+\text{SoftThreshold}\Big{(}A_{i,j}-$ $\displaystyle
t\cdot\Big{(}\left[A\left(g\odot\mathcal{M}_{A}\right)A\right]_{i,j}$
$\displaystyle+$
$\displaystyle\alpha\left[A\left(\mathcal{G}^{-i,-j}\odot\mathcal{M}_{A}\right)A\right]_{i,j}\Big{)},t\cdot
C_{i,j}\Big{)}\,.$

Notice that:

$\displaystyle\alpha\left[A\left(\mathcal{G}^{-i,-j}\odot\mathcal{M}_{A}\right)A\right]_{i,j}$
$\displaystyle=$
$\displaystyle\alpha\left[A\left(\mathcal{G}\odot\mathcal{M}_{A}\right)A\right]_{i,j}-C_{i,j}\cdot\mathcal{G}_{i,j}\cdot[\mathcal{M}_{A}]_{i,j}$
$\displaystyle=$ $\displaystyle\left[\alpha
A\left(\mathcal{G}\odot\mathcal{M}_{A}\right)A-C\odot\left(\mathcal{G}\odot\mathcal{M}_{A}\right)\right]_{i,j}\,,$

where $\odot$ is the Hadmard product. To write (.2) more compactly, first
define:

$\displaystyle B=A\left(g\odot\mathcal{M}_{A}\right)A+\alpha
A\left(\mathcal{G}\odot\mathcal{M}_{A}\right)A-C\odot\left(\mathcal{G}\odot\mathcal{M}_{A}\right)\,,$

then we get

$\displaystyle x=-A_{i,j}+\text{SoftThreshold}(A_{i,j}-t\cdot B_{i,j},t\cdot
C_{i,j})$ $\displaystyle\Rightarrow D=-A+\text{SoftThreshold}(A-t\cdot
B,t\cdot C)\,.$

## References

  * [1] O. Banerjee, L. El Ghaoui, and A. d’Aspremont, Model selection through sparse maximum likelihood estimation for multivariate gaussian or binary data, J. Mach. Learn. Res., (2008).
  * [2] O. Banerjee, L. E. Ghaoui, A. d’Aspremont, and G. Natsoulis, Convex optimization techniques for fitting sparse gaussian graphical models, in Proceedings of the 23rd International Conference on Machine Learning, 2006, p. 89–96.
  * [3] A. Beck and M. Teboulle, A fast iterative shrinkage-thresholding algorithm for linear inverse problems, SIAM Journal on Imaging Sciences, (2009), pp. 183–202.
  * [4] M. Bollhofer, A. Eftekhari, S. Scheidegger, and O. Schenk, Large-scale sparse inverse covariance matrix estimation, SIAM Journal on Scientific Computing, (2019), pp. A380–A401.
  * [5] X. Chen, Y. Liu, H. Liu, and J. Carbonell, Learning spatial-temporal varying graphs with applications to climate data analysis, AAAI Conference on Artificial Intelligence, (2010).
  * [6] A. d’Aspremont, O. Banerjee, and L. Ghaoui, First-order methods for sparse covariance selection, SIAM Journal on Matrix Analysis and Applications, (2006).
  * [7] A. P. Dempster, Covariance selection, Biometrics, (1972), pp. 157–175.
  * [8] A. Dobra, C. Hans, M. Jones, J. Nevins, G. Yao, and M. West, Sparse graphical models for exploring gene expression data, Journal of Multivariate Analysis, (2004), pp. 196–212.
  * [9] J. Duchi, S. Gould, and D. Koller, Projected subgradient methods for learning sparse gaussians, in Uncertainty in Artificial Intelligence, 2008, p. 153–160.
  * [10] J. Fan, J. Zhang, and K. Yu, Vast portfolio selection with gross-exposure constraints, Journal of the American Statistical Association, (2012), pp. 592–606.
  * [11] S. E. Finder, E. Treister, and O. Freifeld, Effective learning of a GMRF mixture model, IEEE Access, (2022).
  * [12] J. Friedman, T. Hastie, and R. Tibshirani, Sparse inverse covariance estimation with the graphical lasso, Biostatistics, (2007), pp. 432–441.
  * [13] A. Goldenberg and A. W. Moore, Bayes net graphs to understand co-authorship networks?, in Proceedings of the 3rd International Workshop on Link Discovery, 2005, p. 1–8.
  * [14] D. Guillot, B. Rajaratnam, B. Rolfs, A. Maleki, and I. Wong, Iterative thresholding algorithm for sparse inverse covariance estimation, Advances in Neural Information Processing Systems, (2012).
  * [15] T. Hastie, R. Tibshirani, and M. Wainwright, Statistical Learning with Sparsity: The Lasso and Generalizations, Chapman & Hall/CRC, 2015.
  * [16] C.-J. Hsieh, M. A. Sustik, I. S. Dhillon, and P. Ravikumar, Quic: Quadratic approximation for sparse inverse covariance estimation, Journal of Machine Learning Research, (2014).
  * [17] C.-J. Hsieh, M. A. Sustik, I. S. Dhillon, P. K. Ravikumar, and R. Poldrack, Big & quic: Sparse inverse covariance estimation for a million variables, in Neural Information Processing Systems, 2013, pp. 3165–3173.
  * [18] T. Idé, A. C. Lozano, N. Abe, and Y. Liu, Proximity-Based Anomaly Detection using Sparse Structure Learning, 2009, pp. 97–108.
  * [19] S. Lauritzen, Graphical models, Oxford Statistical Science Series, Clarendon Press, 1996.
  * [20] L. Li and K.-C. Toh, An inexact interior point method for l1-regularized sparse covariance selection, Mathematical Programming Computation, (2010), p. 291–315.
  * [21] Z. Lu, Smooth optimization approach for sparse covariance selection, SIAM Journal on Optimization, (2009).
  * [22] P. A. Olsen, F. Oztoprak, J. Nocedal, and S. J. Rennie, Newton-like methods for sparse inverse covariance estimation, in Neural Information Processing Systems, 2012, p. 755–763.
  * [23] H. Rue and L. Held, Gaussian Markov random fields: theory and applications, CRC press, 2005.
  * [24] K. Scheinberg, S. Ma, and D. Goldfarb, Sparse inverse covariance selection via alternating linearization methods, in Neural Information Processing Systems, 2010, p. 2101–2109.
  * [25] R. Tibshirani, Regression shrinkage and selection via the lasso, Journal of the Royal Statistical Society. Series B (Methodological), (1996), pp. 267–288.
  * [26] E. Treister and J. S. Turek, A block-coordinate descent approach for large-scale sparse inverse covariance estimation, in Neural Information Processing Systems (NIPS), 2014.
  * [27] E. Treister, J. S. Turek, and I. Yavneh, A multilevel framework for sparse optimization with application to inverse covariance estimation and logistic regression, SIAM Journal on Scientific Computing, 38 (2016), pp. S566–S592.
  * [28] E. Treister and I. Yavneh, A multilevel iterated-shrinkage approach to $l\\_${$1$} penalized least-squares minimization, IEEE transactions on signal processing, 60 (2012), pp. 6319–6329.
  * [29] S. J. Wright, R. D. Nowak, and M. A. T. Figueiredo, Sparse reconstruction by separable approximation, IEEE Transactions on Signal Processing, (2009), pp. 2479–2493.
  * [30] X. Xuan and K. Murphy, Modeling changing dependency structure in multivariate time series, in Proceedings of the 24th International Conference on Machine Learning, 2007, p. 1055–1062.
  * [31] M. Yuan and Y. Lin, Model selection and estimation in the gaussian graphical model, Biometrika, (2007), pp. 19–35.