Coordinate descent algorithm for covariance graphical lasso

Abstract

Bien and Tibshirani (Biometrika, 98(4):807–820, 2011) have proposed a covariance graphical lasso method that applies a lasso penalty on the elements of the covariance matrix. This method is definitely useful because it not only produces sparse and positive definite estimates of the covariance matrix but also discovers marginal independence structures by generating exact zeros in the estimated covariance matrix. However, the objective function is not convex, making the optimization challenging. Bien and Tibshirani (Biometrika, 98(4):807–820, 2011) described a majorize-minimize approach to optimize it. We develop a new optimization method based on coordinate descent. We discuss the convergence property of the algorithm. Through simulation experiments, we show that the new algorithm has a number of advantages over the majorize-minimize approach, including its simplicity, computing speed and numerical stability. Finally, we show that the cyclic version of the coordinate descent algorithm is more efficient than the greedy version.

Appendix: Details of the majorize-minimize MM algorithm in Sect. 4

Consider $\sqrt{\sigma_{ij}^{2}+\epsilon}$ as an approximation to |σ _ij | for a small ϵ>0. Consequently, the original objective function (1) can be approximated by

$$ \operatorname{log}(\det\boldsymbol{\Sigma})+\operatorname {tr}\bigl(\mathbf{S}\boldsymbol{\Sigma}^{-1}\bigr) + 2\rho \sum _{i<j}\sqrt{\sigma_{ij}^2+ \epsilon} + \rho\sum_{i} \sigma_{ii}. $$

(9)

Note the inequality

$$\sqrt{\sigma_{ij}^2+\epsilon}\leq\sqrt{ \bigl(\sigma_{ij}^{(k)}\bigr)^2+\epsilon} + {\sigma_{ij}^2- (\sigma_{ij}^{(k)})^2 \over2 \sqrt{(\sigma_{ij}^{(k)})^2+\epsilon} }, $$

for a fixed $\sigma_{ij}^{(k)}$ and all σ _ij . Then (9) is majorized by

(10)

The minimize-step in MM then minimizes (10) along each column (row) of Σ. Without loss of generality, consider the last column and row. Partition Σ and S as in (2) and consider the same transformation from (σ ₁₂,σ ₂₂) to $(\boldsymbol{\beta}= \boldsymbol{\sigma}_{12},\gamma=\sigma_{22} - \boldsymbol{\sigma}_{12}^{\prime}\boldsymbol{\Sigma}_{11}^{-1} \boldsymbol{\sigma}_{12})$. The four terms in (10) can be written as functions of (β,γ)

where c ₁, c ₂ and c ₃ are constants not involving (β,γ). Dropping off c ₁,c ₂ and c ₃ from (10), we only have to minimize

(11)

For γ, it is easy to derive from (11) that the conditional minimum point given β is the same as in (5). For β, (11) can be written as a function of β,

$$Q\bigl(\boldsymbol{\beta}\mid\gamma, \varSigma^{(k)}\bigr) = \boldsymbol{\beta }^\prime\bigl(\mathbf{V}+\rho\mathbf{D}^{-1} \bigr) \boldsymbol{\beta}-2 \mathbf{u}^\prime\boldsymbol{\beta}, $$

where V and u are defined in (6). This implies that the conditional minimum point of β is β=(V+ρ D ⁻¹)⁻¹ u. Cycling through every column always drives down the approximated objective function (9). In our implementation, the value of ϵ is chosen as follows. The approximation error of (9) to (1) is

Note that the algorithm stops when the change of the objective function is less than 10⁻³. We choose ϵ such that $\rho p(p-1) \sqrt{\epsilon}=0.001$, i.e., ϵ=(0.001/(ρ(p−1)p))², to ensure that the choice of ϵ has no more influence on the estimated Σ than the stopping rule.