The Graphical Lasso: New Insights and Alternatives

The graphical lasso \citep{FHT2007a} is an algorithm for learning the structure in an undirected Gaussian graphical model, using $\ell_1$ regularization to control the number of zeros in the precision matrix ${\B\Theta}={\B\Sigma}^{-1}$ \citep{BGA2008,yuan_lin_07}. The {\texttt R} package \GL\ \citep{FHT2007a} is popular, fast, and allows one to efficiently build a path of models for different values of the tuning parameter. Convergence of \GL\ can be tricky; the converged precision matrix might not be the inverse of the estimated covariance, and occasionally it fails to converge with warm starts. In this paper we explain this behavior, and propose new algorithms that appear to outperform \GL. By studying the"normal equations"we see that, \GL\ is solving the {\em dual} of the graphical lasso penalized likelihood, by block coordinate ascent; a result which can also be found in \cite{BGA2008}. In this dual, the target of estimation is $\B\Sigma$, the covariance matrix, rather than the precision matrix $\B\Theta$. We propose similar primal algorithms \PGL\ and \DPGL, that also operate by block-coordinate descent, where $\B\Theta$ is the optimization target. We study all of these algorithms, and in particular different approaches to solving their coordinate sub-problems. We conclude that \DPGL\ is superior from several points of view.


Introduction
Consider a data matrix X n×p , a sample of n realizations from a p-dimensional Gaussian distribution with zero mean and positive definite covariance matrix Σ.The task is to estimate the unknown Σ based on the n samples -a challenging problem especially when n ≪ p, when the ordinary maximum likelihood estimate does not exist.Even if it does exist (for p ≤ n), the MLE is often poorly behaved, and regularization is called for.The Graphical Lasso [Friedman et al., 2007] is a regularization framework for estimating the covariance matrix Σ, under the assumption that its inverse Θ = Σ −1 is sparse [Banerjee et al., 2008, Yuan and Lin, 2007, Meinshausen and Bühlmann, 2006].Θ is called the precision matrix; if an element θ jk = 0, this implies that the corresponding variables X j and X k are conditionally independent, given the rest.Our algorithms focus either on the restricted version of Θ or its inverse W = Θ −1 .The graphical lasso problem minimizes a ℓ 1 -regularized negative log-likelihood: minimize Θ≻0 f (Θ) := − log det(Θ) + tr(SΘ) + λ Θ 1 . (1) Here S is the sample covariance matrix, Θ 1 denotes the sum of the absolute values of Θ, and λ is a tuning parameter controlling the amount of ℓ 1 shrinkage.This is a semidefinite programming problem (SDP) in the variable Θ [Boyd and Vandenberghe, 2004].
In this paper we revisit the glasso algorithm proposed by Friedman et al. [2007] for solving (1); we analyze its properties, expose problems and issues, and propose alternative algorithms more suitable for the task.Some of the results and conclusions of this paper can be found in Banerjee et al. [2008], both explicitly and implicitly.We re-derive some of the results and derive new results, insights and algorithms, using a unified and more elementary framework.
Notation We denote the entries of a matrix A n×n by a ij .A 1 denotes the sum of its absolute values, A ∞ the maximum absolute value of its entries, A F is its Frobenius norm, and abs(A) is the matrix with elements |a ij |.For a vector u ∈ ℜ q , u 1 denotes the ℓ 1 norm, and so on.
From now on, unless otherwise specified, we will assume that λ > 0.
2 Review of the glasso algorithm.
We use the frame-work of "normal equations" as in Hastie et al. [2009], Friedman et al. [2007].Using sub-gradient notation, we can write the optimality conditions (aka "normal equations") for a solution to (1) as where Γ is a matrix of component-wise signs of Θ: (we use the notation γ jk ∈ Sign(θ jk )).Since the global stationary conditions of (2) require θ jj to be positive, this implies that where W = Θ −1 .glasso uses a block-coordinate method for solving (2).Consider a partitioning of Θ and Γ: where Θ 11 is (p − 1) × (p − 1), θ 12 is (p − 1) × 1 and θ 22 is scalar.W and S are partitioned the same way.Using properties of inverses of block-partitioned matrices, observe that W = Θ −1 can be written in two equivalent forms: glasso solves for a row/column of (2) at a time, holding the rest fixed.Considering the pth column of (2), we get − w 12 + s 12 + λγ 12 = 0.
The glasso algorithm solves (10) for β = θ 12 /θ 22 , that is where γ 12 ∈ Sign(β), since θ 22 > 0. ( 13) is the stationarity equation for the following ℓ 1 regularized quadratic program: minimize where W 11 ≻ 0 is assumed to be fixed.This is analogous to a lasso regression problem of the last variable on the rest, except the cross-product matrix S 11 is replaced by its current estimate W 11 .This problem itself can be solved efficiently using elementwise coordinate descent, exploiting the sparsity in β.From β, it is easy to obtain ŵ12 from (9).Using the lower-right element of ( 6), θ22 is obtained by Finally, θ12 can now be recovered from β and θ22 .Notice, however, that having solved for β and updated w 12 , glasso can move onto the next block; disentangling θ 12 and θ 22 can be done at the end, when the algorithm over all blocks has converged.The glasso algorithm is outlined in Algorithm 1.We show in Lemma 3 in Section 8 that the successive updates in glasso keep W positive definite.
2. Cycle around the columns repeatedly, performing the following steps till convergence: (a) Rearrange the rows/columns so that the target column is last (implicitly).
(b) Solve the lasso problem (14), using as warm starts the solution from the previous round for this column.
(d) Save β for this column in the matrix B.
3. Finally, for every row/column, compute the diagonal entries θjj using (15), and convert the B matrix to Θ.  (k) ) for the sequence of solutions produced by glasso on an example.Surprisingly, the curve is not monotone decreasing, as confirmed by the middle plot.If glasso were solving (1) by block coordinate-descent, we would not anticipate this behavior.
A closer look at steps ( 9) and (10) of the glasso algorithm leads to the following observations: (a) We wish to solve (8) for θ 12 .However θ 12 is entangled in W 11 , which is (incorrectly) treated as a constant.
(b) After updating θ 12 , we see from (7) that the entire (working) covariance matrix W changes. glasso however updates only w 12 and w 21 .These two observations explain the non-monotone behavior of glasso in minimizing f (Θ).Section 3 shows a corrected block-coordinate descent algorithm for Θ, and Section 4 shows that the glasso algorithm is actually optimizing the dual of problem (1), with the optimization variable being W.

A Corrected glasso block coordinate-descent algorithm
Recall that ( 12) is a variant of (10), where the dependence of the covariance sub-matrix W 11 on θ 12 is explicit.With α = θ 12 w 22 (with w 22 ≥ 0 fixed), Θ 11 ≻ 0, ( 12) is equivalent to the stationary condition for minimize If α is the minimizer of ( 16), then θ12 = α/w 22 .To complete the optimization for the entire row/column we need to update θ 22 .This follows simply from (7) with w 22 = s 22 + λ.
To solve (16) we need Θ −1 11 for each block update.We achieve this by maintaining W = Θ −1 as the iterations proceed.Then for each block • once θ 12 is updated, the entire working covariance matrix W is updated (in particular the portions W 11 and w 12 ), via the identities in (7), using the known Θ −1 11 .
Both these steps are simple rank-one updates with a total cost of O(p 2 ) operations.We refer to this as the primal graphical lasso or p-glasso, which we present in Algorithm 2. The p-glasso algorithm requires slightly more work than glasso, since an additional O(p 2 ) operations have to be performed before and after each block update.In return we have that after every row/column update, Θ and W are positive definite (for λ > 0) and ΘW = I p .
2. Cycle around the columns repeatedly, performing the following steps till convergence: (a) Rearrange the rows/columns so that the target column is last (implicitly).
4 What is glasso actually solving?
Building upon the framework developed in Section 2, we now proceed to establish that glasso solves the convex dual of problem (1), by block coordinate ascent.We reach this conclusion via elementary arguments, closely aligned with the framework we develop in Section 2. The approach we present here is intended for an audience without much of a familiarity with convex duality theory Boyd and Vandenberghe [2004].
Figure 1 illustrates that glasso is an ascent algorithm on the dual of the problem 1.The red curve in the left plot shows the dual objective rising monotonely, and the rightmost plot shows that the increments are indeed positive.There is an added twist though: in solving the block-coordinate update, glasso solves instead the dual of that subproblem.

Dual of the ℓ 1 regularized log-likelihood
We present below the following lemma, the conclusion of which also appears in Banerjee et al. [2008], but we use the framework developed in Section 2.
Notice that for the dual, the optimization variable is Γ, with S + Γ = Θ −1 = W.In other words, the dual problem solves for W rather than Θ, a fact that is suggested by the glasso algorithm.
Remark 1.The equivalence of the solutions to problems ( 19) and (1) as described above can also be derived via convex duality theory [Boyd and Vandenberghe, 2004], which shows that ( 19) is a dual function of the ℓ 1 regularized negative log-likelihood (1).Strong duality holds, hence the optimal solutions of the two problems coincide Banerjee et al. [2008].
We now consider solving ( 22) for the last block γ12 (excluding diagonal), holding the rest of Γ fixed.The corresponding equations are The only non-trivial translation is the θ 12 in the first equation.We must express this in terms of the optimization variable γ12 .Since s 12 + γ12 = w 12 , using the identities in ( 6), we have W −1 11 (s 12 + γ12 ) = −θ 12 /θ 22 .Since θ 22 > 0, we can redefine p12 = p 12 /θ 22 , to get The following lemma shows that a block update of glasso solves (24) (and hence ( 23)), a block of stationary conditions for the dual of the graphical lasso problem.Curiously, glasso does this not directly, but by solving the dual of the QP corresponding to this block of equations.
Note that the QP ( 27) is a (partial) optimization over the variable w 12 only (since s 12 is fixed); the sub-matrix W 11 remains fixed in the QP.Exactly one row/column of W changes when the blockcoordinate algorithm of glasso moves to a new row/column, unlike an explicit full matrix update in W 11 , which is required if θ 12 is updated.This again emphasizes that glasso is operating on the covariance matrix instead of Θ.We thus arrive at the following conclusion: Theorem 1. glasso performs block-coordinate ascent on the box-constrained SDP ( 19), the Lagrange dual of the primal problem (1).Each of the block steps are themselves box-constrained QPs, which glasso optimizes via their Lagrange duals.
In our annotation perhaps glasso should be called dd-glasso, since it performs dual block updates for the dual of the graphical lasso problem.Banerjee et al. [2008], the paper that inspired the original glasso article [Friedman et al., 2007], also operates on the dual.They however solve the block-updates directly (which are box constrained QPs) using interior-point methods.

A New Algorithmdp-glasso
In Section 3, we described p-glasso, a primal coordinate-descent method.For every row/column we need to solve a lasso problem ( 16), which operates on a quadratic form corresponding to the square matrix Θ −1 11 .There are two problems with this approach: • the matrix Θ −1 11 needs to be constructed at every row/column update with complexity O(p 2 ); • Θ −1 11 is dense.
We now show how a simple modification of the ℓ 1 -regularized QP leads to a box-constrained QP with attractive computational properties.
The KKT optimality conditions for (16), following (12), can be written as: Along the same lines of the derivations used in Lemma 2, the condition above is equivalent to for some vector (with non-negative entries) q12 .( 32) are the KKT optimality conditions for the following box-constrained QP: The optimal solutions of ( 33) and ( 31) are related by a consequence of (31), with α = θ12 • w 22 and w 22 = s 22 + λ.The diagonal θ 22 of the precision matrix is updated via (7): Algorithm 3 dp-glasso algorithm 1. Initialize Θ = diag(S + λI) −1 .
2. Cycle around the columns repeatedly, performing the following steps till convergence: (a) Rearrange the rows/columns so that the target column is last (implicitly).By strong duality, the box-constrained QP (33) with its optimality conditions ( 32) is equivalent to the lasso problem ( 16).Now both the problems listed at the beginning of the section are removed.The problem matrix Θ 11 is sparse, and no O(p 2 ) updating is required after each block.
The solutions returned at step 2(b) for θ12 need not be exactly sparse, even though it purports to produce the solution to the primal block problem ( 16), which is sparse.One needs to use a tight convergence criterion when solving (33).In addition, one can threshold those elements of θ12 for which γ is away from the box boundary, since those values are known to be zero.
Note that dp-glasso does to the primal formulation ( 1) what glasso does to the dual.dpglasso operates on the precision matrix, whereas glasso operates on the covariance matrix.

Computational Costs in Solving the Block QPs
The ℓ 1 regularized QPs appearing in ( 14) and ( 16) are of the generic form minimize for A ≻ 0. In this paper, we choose to use cyclical coordinate descent for solving (36), as it is used in the glasso algorithm implementation of Friedman et al. [2007].Moreover, cyclical coordinate descent methods perform well with good warm-starts.These are available for both ( 14) and ( 16), since they both maintain working copies of the precision matrix, updated after every row/column update.There are other efficient ways for solving (36), capable of scaling to large problemsfor example first-order proximal methods [Beck andTeboulle, 2009, Nesterov, 2007], but we do not pursue them in this paper.
The box-constrained QPs appearing in ( 27) and ( 33) are of the generic form: for some Ã ≻ 0. As in the case above, we will use cyclical coordinate-descent for optimizing (37).
In general it is more efficient to solve (36) than (37) for larger values of λ.This is because a large value of λ in (36) results in sparse solutions û; the coordinate descent algorithm can easily detect when a zero stays zero, and no further work gets done for that coordinate on that pass.If the solution to (36) has κ non-zeros, then on average κ coordinates need to be updated.This leads to a cost of O(qκ), for one full sweep across all the q coordinates.
On the other hand, a large λ for (37) corresponds to a weakly-regularized solution.Cyclical coordinate procedures for this task are not as effective.Every coordinate update of v results in updating the gradient, which requires adding a scalar multiple of a column of Ã.If Ã is dense, this leads to a cost of O(q), and for one full cycle across all the coordinates this costs O(q 2 ), rather than the O(qκ) for (36).However, our experimental results show that dp-glasso is more efficient than glasso, so there are some other factors in play.When Ã is sparse, there are computational savings.If Ã has κq non-zeros, the cost per column reduces on average to O(κq) from O(q 2 ).For the formulation (33) Ã is Θ 11 , which is sparse for large λ.Hence for large λ, glasso and dp-glasso have similar costs.
For smaller values of λ, the box-constrained QP (37) is particularly attractive.Most of the coordinates in the optimal solution v will pile up at the boundary points {−λ, λ}, which means that the coordinates need not be updated frequently.For problem (33) this number is also κ, the number of non-zero coefficients in the corresponding column of the precision matrix.If κ of the coordinates pile up at the boundary, then one full sweep of cyclical coordinate descent across all the coordinates will require updating gradients corresponding to the remaining q − κ coordinates.Using similar calculations as before, this will cost O(q(q − κ)) operations per full cycle (since for small λ, Ã will be dense).For the ℓ 1 regularized problem (36), no such saving is achieved, and the cost is O(q 2 ) per cycle.
Note that to solve problem (1), we need to solve a QP of a particular type (36) or (37) for a certain number of outer cycles (ie full sweeps across rows/columns).For every row/column update, the associated QP requires varying number of iterations to converge.It is hard to characterize all these factors and come up with precise estimates of convergence rates of the overall algorithm.However, we have observed that with warm-starts, on a relatively dense grid of λs, the complexities given above are pretty much accurate for dp-glasso (with warmstarts) specially when one is interested in solutions with small / moderate accuracy.Our experimental results in Section 9.1 and Appendix Section B support our observation.
We will now have a more critical look at the updates of the glasso algorithm and study their properties.

glasso: Positive definiteness, Sparsity and Exact Inversion
As noted earlier, glasso operates on W -it does not explicitly compute the inverse W −1 .It does however keep track of the estimates for θ 12 after every row/column update.The copy of Θ retained by glasso along the row/column updates is not the exact inverse of the optimization variable W. The precision matrix produced after every row/column update need not be the exact inverse of the working covariance matrix -the squared Frobenius norm of the error is being plotted across iterations.[Right Panel] The estimated precision matrix Θ produced by glasso need not be positive definite along iterations; plot shows minimal eigen-value.
In many real-life problems one only needs an approximate solution to (1): • for computational reasons it might be impractical to obtain a solution of high accuracy; • from a statistical viewpoint it might be sufficient to obtain an approximate solution for Θ that is both sparse and positive definite It turns out that the glasso algorithm is not suited to this purpose.Since the glasso is a block coordinate procedure on the covariance matrix, it maintains a positive definite covariance matrix at every row/column update.However, since the estimated precision matrix is not the exact inverse of W, it need not be positive definite.Although it is relatively straightforward to maintain an exact inverse of W along the row/column updates (via simple rank-one updates as before), this inverse W −1 need not be sparse.Arbitrary thresholding rules may be used to set some of the entries to zero, but that might destroy the positive-definiteness of the matrix.Since a principal motivation of solving (1) is to obtain a sparse precision matrix (which is also positive definite), returning a dense W −1 to (1) is not desirable.
Figures 2 illustrates the above observations on a typical example.
The dp-glasso algorithm operates on the primal (1).Instead of optimizing the ℓ 1 regularized QP ( 16), which requires computing Θ −1 11 , dp-glasso optimizes (33).After every row/column update the precision matrix Θ is positive definite.The working covariance matrix maintained by dp-glasso via w 12 := s 12 + γ need not be the exact inverse of Θ. covariance matrix estimates, if required, can be obtained by tracking Θ −1 via simple rank-one updates, as described earlier.
Unlike glasso, dp-glasso (and p-glasso) return a sparse and positive definite precision matrix even if the row/column iterations are terminated prematurely.

Warm Starts and Path-seeking Strategies
Since we seldom know in advance a good value of λ, we often compute a sequence of solutions to (1) for a (typically) decreasing sequence of values λ 1 > λ 2 > . . .> λ K .Warm-start or continuation methods use the solution at λ i as an initial guess for the solution at λ i+1 , and often yield great efficiency.It turns out that for algorithms like glasso which operate on the dual problem, not all warm-starts necessarily lead to a convergent algorithm.We address this aspect in detail in this section.
The following lemma states the conditions under which the row/column updates of the glasso algorithm will maintain positive definiteness of the covariance matrix W. Lemma 3. Suppose Z is used as a warm-start for the glasso algorithm.If Z ≻ 0 and Z−S ∞ ≤ λ, then every row/column update of glasso maintains positive definiteness of the working covariance matrix W.
Proof.Recall that the glasso solves the dual (19).Assume Z is partitioned as in (5), and the pth row/column is being updated.Since Z ≻ 0, we have both Since Z 11 remains fixed, it suffices to show that after the row/column update, the expression ( ŵ22 − ŵ21 (Z 11 ) −1 ŵ12 ) remains positive.Recall that, via standard optimality conditions we have ŵ22 = Combining the above along with the fact that ŵ22 ≥ z 22 we see which implies that the new covariance estimate W ≻ 0.
Remark 3. If the condition Z − S ∞ ≤ λ appearing in Lemma 3 is violated, then the row/column update of glasso need not maintain PD of the covariance matrix W.
We have encountered many counter-examples that show this to be true, see the discussion below.
The R package implementation of glasso allows the user to specify a warm-start as a tuple (Θ 0 , W 0 ).This option is typically used in the construction of a path algorithm.
If ( Θ λ , W λ ) is provided as a warm-start for λ ′ < λ, then the glasso algorithm is not guaranteed to converge.It is easy to find numerical examples by choosing the gap λ − λ ′ to be large enough.Among the various examples we encountered, we briefly describe one here.Details of the experiment/data and other examples can be found in the online Appendix A.1.We generated a data-matrix X n×p , with n = 2, p = 5 with iid standard Gaussian entries.S is the sample covariance matrix.We solved problem (1) using glasso for λ = 0.9 × max i =j |s ij |.We took the estimated covariance and precision matrices: W λ and Θ λ as a warm-start for the glasso algorithm with λ ′ = λ × 0.01.The glasso algorithm failed to converge with this warm-start.We note that W λ − S ∞ = 0.0402 λ ′ (hence violating the sufficient condition in Lemma 4) and after updating the first row/column via the glasso algorithm we observed that "covariance matrix" W has negative eigen-values -leading to a non-convergent algorithm.The above phenomenon is not surprising and easy to explain and generalize.Since W λ solves the dual ( 19), it is necessarily of the form W λ = S + Γ, for Γ ∞ ≤ λ.
In the light of Lemma 3 and also Remark 3, the warm-start needs to be dual-feasible in order to guarantee that the iterates W remain PD and hence for the sub-problems to be well defined convex programs.Clearly W λ does not satisfy the box-constraint W λ − S ∞ ≤ λ ′ , for λ ′ < λ.However, in practice the glasso algorithm is usually seen to converge (numerically) when λ ′ is quite close to λ.
The following lemma establishes that any PD matrix can be taken as a warm-start for p-glasso or dp-glassoto ensure a convergent algorithm.
Lemma 4. Suppose Φ ≻ 0 is a used as a warm-start for the p-glasso (or dp-glasso) algorithm.Then every row/column update of p-glasso (or dp-glasso) maintains positive definiteness of the working precision matrix Θ.
Note that the block Φ 11 remains fixed; only the pth row/column of Θ changes.φ 21 gets updated to θ21 , as does θ12 .From (7) the updated diagonal entry θ22 satisfies: Thus the updated matrix Θ remains PD.The result for the dp-glasso algorithm follows, since both the versions p-glasso and dp-glasso solve the same block coordinate problem.
Remark 4. A simple consequence of Lemmas 3 and 4 is that the QPs arising in the process, namely the ℓ 1 regularized QPs ( 14), ( 16) and the box-constrained QPs ( 27) and ( 33) are all valid convex programs, since all the respective matrices W 11 , Θ −1 11 and W −1 11 , Θ 11 appearing in the quadratic forms are PD.
As exhibited in Lemma 4, both the algorithms dp-glasso and p-glasso are guaranteed to converge from any positive-definite warm start.This is due to the unconstrained formulation of the primal problem (1).
glasso really only requires an initialization for W, since it constructs Θ on the fly.Likewise dp-glasso only requires an initialization for Θ.Having the other half of the tuple assists in the block-updating algorithms.For example, glasso solves a series of lasso problems, where Θ play the role as parameters.By supplying Θ along with W, the block-wise lasso problems can be given starting values close to the solutions.The same applies to dp-glasso.In neither case do the pairs have to be inverses of each other to serve this purpose.
If we wish to start with inverse pairs, and maintain such a relationship, we have described earlier how O(p 2 ) updates after each block optimization can achieve this.One caveat for glasso is that starting with an inverse pair costs O(p 3 ) operations, since we typically start with W = S + λI.For dp-glasso, we typically start with a diagonal matrix, which is trivial to invert.

Experimental Results & Timing Comparisons
We compared the performances of algorithms glasso and dp-glasso (both with and without warmstarts) on different examples with varying (n, p) values.While most of the results are presented in this section, some are relegated to the online Appendix B. Section 9.1 describes some synthetic examples and Section 9.2 presents comparisons on a real-life micro-array data-set.

Synthetic Experiments
In this section we present examples generated from two different covariance models -as characterized by the covariance matrix Σ or equivalently the precision matrix Θ.We create a data matrix X n×p by drawing n independent samples from a p dimensional normal distribution MVN(0, Σ).The sample covariance matrix is taken as the input S to problem (1).The two covariance models are described below: Type-1 The population concentration matrix Θ = Σ −1 has uniform sparsity with approximately 77 % of the entries zero.
We created the covariance matrix as follows.We generated a matrix B with iid standard Gaussian entries, symmetrized it via 1 2 (B + B ′ ) and set approximately 77% of the entries of this matrix to zero, to obtain B (say).We added a scalar multiple of the p dimensional identity matrix to B to get the precision matrix Θ = B + ηI p×p , with η chosen such that the minimum eigen value of Θ is one.Type-2 This example, taken from Yuan and Lin [2007], is an auto-regressive process of order two -the precision matrix being tri-diagonal: For each of the two set-ups Type-1 and Type-2 we consider twelve different combinations of (n, p): For every (n, p) we solved (1) on a grid of twenty λ values linearly spaced in the log-scale, with λ i = 0.8 i × {0.9λ max }, i = 1, . . ., 20, where λ max = max i =j |s ij |, is the off-diagonal entry of S with largest absolute value.λ max is the smallest value of λ for which the solution to (1) is a diagonal matrix.
Since this article focuses on the glasso algorithm, its properties and alternatives that stem from the main idea of block-coordinate optimization, we present here the performances of the following algorithms: Dual-Cold glasso with initialization W = S + λI p×p , as suggested in Friedman et al. [2007].

Dual-Warm
The path-wise version of glasso with warm-starts, as suggested in Friedman et al. [2007].Although this path-wise version need not converge in general, this was not a problem in our experiments, probably due to the fine-grid of λ values.
Primal-Warm The path-wise version of dp-glasso with warm-starts.
We did not include p-glasso in the comparisons above since p-glasso requires additional matrix rank-one updates after every row/column update, which makes it more expensive.None of the above listed algorithms require matrix inversions (via rank one updates).Furthermore, dp-glasso and p-glasso are quite similar as both are doing a block coordinate optimization on the dual.Hence we only included dp-glasso in our comparisons.We used our own implementation of the glasso and dp-glasso algorithm in R. The entire program is written in R, except the inner block-update solvers, which are the real work-horses: • For glasso we used the lasso code crossProdLasso written in FORTRAN by Friedman et al. [2007]; • For dp-glasso we wrote our own FORTRAN code to solve the box QP.
An R package implementing dp-glasso will be made available in CRAN.
In the figure and tables that follow below, for every algorithm, at a fixed λ we report the total time taken by all the QPs -the ℓ 1 regularized QP for glasso and the box constrained QP for dpglasso till convergence All computations were done on a Linux machine with model specs: Intel(R) Xeon(R) CPU 5160 @ 3.00GHz.
Convergence Criterion: Since dp-glasso operates on the the primal formulation and glasso operates on the dual -to make the convergence criteria comparable across examples we based it on the relative change in the primal objective values i.e. f (Θ) (1) across two successive iterations: where one iteration refers to a full sweep across p rows/columns of the precision matrix (for dp-glasso ) and covariance matrix (for glasso ); and TOL denotes the tolerance level or level of accuracy of the solution.To compute the primal objective value for the glasso the precision matrix is computed from W via direct inversion (the time taken for inversion and objective value computation is not included in the timing comparisons).
Computing the objective function is quite expensive relative to the computational cost of the iterations.In our experience convergence criteria based on a relative change in the precision matrix for dp-glasso and the covariance matrix for glasso seemed to be a practical choice for the examples we considered.However, for reasons we described above, we used criterion 40 in the experiments.
Observations: Figure 4 presents the times taken by the algorithms to converge to an accuracy of TOL = 10 −4 on a grid of λ values.
The figure shows eight different scenarios with p > n, corresponding to the two different covariance models Type-1 (left panel) and Type-2 (right panel).It is quite evident that dp-glasso with warmstarts (Primal-Warm) outperforms all the other algorithms across all the different examples.All the algorithms converge quickly for large values of λ (typically high sparsity) and become slower with decreasing λ.For large p and small λ, convergence is slow; however for p > n, the non-sparse end of the regularization path is really not that interesting from a statistical viewpoint.Warm-starts apparently do not always help in speeding up the convergence of glasso ; for example see Figure 4 with (n, p) = (500, 1000) (Type 1) and (n, p) = (500, 800) (Type 2).This probably further validates the fact that warm-starts in the case of glasso need to be carefully designed, in order for them to speed-up convergence.Note however, that glasso with the warm-starts prescribed is not even guaranteed to converge -we however did not come across any such instance among the experiments presented in this section.
Based on the suggestion of a referee we annotated the plots in Figure 4 with locations in the regularization path that are of interest.For each plot, two vertical dotted lines are drawn which correspond to the λs at which the distance of the estimated precision matrix Θ λ from the population precision matrix is minimized wrt to the • 1 norm (green) and • F norm (blue).The optimal λ corresponding to the • 1 metric chooses sparser models than those chosen by • F ; the performance gains achieved by dp-glasso seem to be more prominent for the latter λ.
Table 1 presents the timings for all the four algorithmic variants on the twelve different (n, p) combinations listed above for Type 1.For every example, we report the total time till convergence on a grid of twenty λ values for two different tolerance levels: TOL ∈ {10 −4 , 10 −5 }.Note that the dp-glasso returns positive definite and sparse precision matrices even if the algorithm is terminated at a relatively small/moderate accuracy level -this is not the case in glasso .The rightmost column presents the proportion of non-zeros averaged across the entire path of solutions Θ λ , where Θ λ is obtained by solving (1) to a high precision i.e. 10 −6 , by algorithms glasso and dp-glasso and averaging the results.
Again we see that in all the examples dp-glasso with warm-starts is the clear winner among its competitors.For a fixed p, the total time to trace out the path generally decreases with increasing n.There is no clear winner between glasso with warm-starts and glasso without warm-starts.It is often seen that dp-glasso without warm-starts converges faster than both the variants of glasso (with and without warm-starts).
Table 2 reports the timing comparisons for Type 2. Once again we see that in all the examples Primal-Warm turns out to be the clear winner.
For n ≤ p = 1000, we observe that Primal-Warm is generally faster for Type-2 than Type-1.This however, is reversed for smaller values of p ∈ {800, 500}.Primal-Cold is has a smaller overall computation time for Type-1 over Type-2.In some cases (for example n ≤ p = 1000), we see that Primal-Warm in Type-2 converges much faster than its competitors on a relative scale than in Type-1 -this difference is due to the variations in the structure of the covariance matrix.

Micro-array Example
We consider the data-set introduced in Alon et al. [1999] and further studied in Rothman et al. [2008], Mazumder and Hastie [2012].In this experiment, tissue samples were analyzed using an Affymetrix Oligonucleotide array.The data was processed, filtered and reduced to a subset of 2000 gene expression values.The number of Colon Adenocarcinoma tissue samples is n = 62.For the purpose of the experiments presented in this section, we pre-screened the genes to a size of p = 725.We obtained this subset of genes using the idea of exact covariance thresholding introduced in our paper [Mazumder and Hastie, 2012].We thresholded the sample correlation matrix obtained from the 62 × 2000 microarray data-matrix into connected components with a threshold of 0.003641 the genes belonging to the largest connected component formed our pre-screened gene pool of size p = 725.This (subset) data-matrix of size (n, p) = (62, 725) is used for our experiments.
The results presented below in Table 3 show timing comparisons of the four different algorithms: Primal-Warm/Cold and Dual-Warm/Cold on a grid of fifteen λ values in the log-scale.Once again we see that the Primal-Warm outperforms the others in terms of speed and accuracy.Dual-Warm performs quite well in this example.

Conclusions
This paper explores some of the apparent mysteries in the behavior of the glasso algorithm introduced in Friedman et al. [2007].These have been explained by leveraging the fact that the glasso algorithm is solving the dual of the graphical lasso problem (1), by block coordinate ascent.Each block update, itself the solution to a convex program, is solved via its own dual, which is equivalent  tolerance levels (TOL).We took a grid of fifteen λ values, the average % of zeros along the whole path is 90.8.
to a lasso problem.The optimization variable is W, the covariance matrix, rather than the target precision matrix Θ.During the course of the iterations, a working version of Θ is maintained, but it may not be positive definite, and its inverse is not W. Tight convergence is therefore essential, for the solution Θ to be a proper inverse covariance.There are issues using warm starts with glasso, when computing a path of solutions.Unless the sequence of λs are sufficiently close, since the "warm start"s are not dual feasible, the algorithm can get into trouble.
We have also developed two primal algorithms p-glasso and dp-glasso.The former is more expensive, since it maintains the relationship W = Θ −1 at every step, an O(p 3 ) operation per sweep across all row/columns.dp-glasso is similar in flavor to glasso except its optimization variable is Θ.It also solves the dual problem when computing its block update, in this case a box-QP.This box-QP has attractive sparsity properties at both ends of the regularization path, as evidenced in some of our experiments.It maintains a positive definite Θ throughout its iterations, and can be started at any positive definite matrix.Our experiments show in addition that dp-glasso is faster than glasso.
An R package implementing dp-glasso will be made available in CRAN.

A Online Appendix
This section complements the examples provided in the paper with further experiments and illustrations.
A With q denoting the maximum off-diagonal entry of S (in absolute value), we solved (1) using glasso at λ = 0.9 × q.The covariance matrix for this λ was taken as a warm-start for the glasso algorithm with λ ′ = λ × 0.01.The smallest eigen-value of the working covariance matrix W produced by the glasso algorithm, upon updating the first row/column was: −0.002896128, which is clearly undesirable for the convergence of the algorithm glasso .This is why the algorithm glasso breaks down.

Example 2:
The example is similar to above, with (n, p) = (10, 50), the seed of random number generator in R being set to set.seed(2008) and X n×p is the data-matrix with iid Gaussian entries.If the covariance matrix W λ which solves problem (1) with λ = 0.9 × max i =j |s ij | is taken as a warm-start to the glasso algorithm with λ ′ = λ × 0.1 -the algorithm fails to converge.Like the previous example, after the first row/column update, the working covariance matrix has negative eigen-values.

B Further Experiments and Numerical Studies
This section is a continuation to Section 9, in that it provides further examples comparing the performance of algorithms glasso and dp-glasso .The experimental data is generated as follows.For a fixed value of p, we generate a matrix A p×p with random Gaussian entries.The matrix is symmetrized by A ← (A + A ′ )/2.Approximately half of the off-diagonal entries of the matrix are set to zero, uniformly at random.All the eigen-values of the matrix A are lifted so that the smallest eigen-value is zero.The noiseless version of the precision matrix is given by Θ = A + τ I p×p .We generated the sample covariance matrix S by adding symmetric positive semi-definite random noise N to Θ −1 ; i.e. S = Θ −1 + N, where this noise is generated in the same manner as A. We considered four different values of p ∈ {300, 500, 800, 1000} and two different values of τ ∈ {1, 4}.
For every p, τ combination we considered a path of twenty λ values on the geometric scale.For every such case four experiments were performed: Primal-Cold, Primal-Warm, Dual-Cold and Dual-Warm (as described in Section 9).Each combination was run 5 times, and the results averaged, to avoid dependencies on machine loads.Figure 4 shows the results.Overall, dp-glasso with warm starts performs the best, especially at the extremes of the path.We gave some explanation for this in Section 6.For the largest problems (p = 1000) their performances are comparable in the central part of the path (though dp-glasso dominates), but at the extremes dp-glasso dominates by a large margin.

Figure 1 (
Figure1(left panel, black curve) plots the objective f (Θ(k) ) for the sequence of solutions produced by glasso on an example.Surprisingly, the curve is not monotone decreasing, as confirmed by the middle plot.If glasso were solving (1) by block coordinate-descent, we would not anticipate this behavior.A closer look at steps (9) and (10) of the glasso algorithm leads to the following observations:

Figure 1 :
Figure 1: [Left panel] The objective values of the primal criterion (1) and the dual criterion (19) corresponding to the covariance matrix W produced by glasso algorithm as a function of the iteration index (each column/row update).[Middle Panel] The successive differences of the primal objective values -the zero crossings indicate non-monotonicity.[Right Panel] The successive differences in the dual objective valuesthere are no zero crossings, indicating that glasso produces a monotone sequence of dual objective values.

Figure 2 Figure 2 :
Figure2: Figure illustrating some negative properties of glasso using a typical numerical example.[Left Panel] The precision matrix produced after every row/column update need not be the exact inverse of the working covariance matrix -the squared Frobenius norm of the error is being plotted across iterations.[Right Panel] The estimated precision matrix Θ produced by glasso need not be positive definite along iterations; plot shows minimal eigen-value.

Figure 3 :
Figure 3: The timings in seconds for the four different algorithmic versions: glasso (with and without warmstarts) and dp-glasso (with and without warm-starts) for a grid of λ values on the log-scale.[Left Panel] Covariance model for Type-1, [Right Panel] Covariance model for Type-2.The horizontal axis is indexed by the proportion of zeros in the solution.The vertical dashed lines correspond to the optimal λ values for which the estimated errors Θ λ − Θ 1 (green) and Θ λ − Θ F (blue) are minimum.

Table 2 :
Table showing comparative timings of the four algorithmic variants of glasso and dp-glasso for the covariance model in Type-2.This table is similar to Table1, displaying results for Type-1.dp-glasso with warm-starts consistently outperforms all its competitors.

Table 3 :
Comparisons among algorithms for a microarray dataset with n = 62 and p = 725, for different .1 Examples: Non-Convergence of glasso with warm-starts This section illustrates with examples that warm-starts for the glasso need not converge.This is a continuation of examples presented in Section 8. We took (n, p) = (2, 5) and setting the seed of the random number generator in R as set.seed(2008)we generated a data-matrix X n×p with iid standard Gaussian entries.The sample covariance matrix S is given below: