Sparse Inverse Covariance Selection via Alternating Linearization Methods

Introduction

Note that (1) can be rewritten as $\min_{X\in S^{n}_{++}}\max_{\|U\|_{\infty}\leq\rho}-\log\det X+\langle\hat{\Sigma}+U,X\rangle,$ where $\|U\|_{\infty}$ is the largest absolute value of the entries of $U$. By exchanging the order of max and min, we obtain the dual problem $\max_{\|U\|_{\infty}\leq\rho}\min_{X\in S^{n}_{++}}-\log\det X+\langle\hat{\Sigma}+U,X\rangle,$ which is equivalent to

Both the primal and dual problems have strictly convex objectives; hence, their optimal solutions are unique. Given a dual solution $W$, $X=W^{-1}$ is primal feasible resulting in the duality gap

Although developed independently, our method is closely related to Yuan’s method . Both methods exploit the special form of the primal problem (1) by alternatingly minimizing one of the terms of the objective function plus an approximation to the other term. The main difference between the two methods is in the construction of these approximations. As we will show, our method has a theoretically justified interpretation and is based on an algorithmic framework with complexity bounds, while no complexity bound is available for Yuan’s method. Also our method has an intuitive interpretation from a learning perspective. Extensive numerical test results on both synthetic data and real problems have shown that our ALM algorithm significantly outperforms other existing algorithms, such as the PSM algorithm proposed by Duchi et al. and the VSM algorithm proposed by Lu . Note that it is shown in and that PSM and VSM outperform the BCD method in and $glasso$ in .

Organization of the paper. In Section 2 we briefly review alternating linearization methods for minimizing the sum of two convex functions and establish convergence and iteration complexity results. We show how to use ALM to solve SICS problems and give intuition from a learning perspective in Section 3. Finally, we present some numerical results on both synthetic and real data in Section 4 and compare ALM with PSM algorithm and VSM algorithm .

Alternating Linearization Methods

We consider here the alternating linearization method (ALM) for solving the following problem:

where $f$ and $g$ are both convex functions. An effective way to solve (4) is to “split” $f$ and $g$ by introducing a new variable, i.e., to rewrite (4) as

and apply an alternating direction augmented Lagrangian method to it. Given a penalty parameter $1/\mu$, at the $k$-th iteration, the augmented Lagrangian method minimizes the augmented Lagrangian function

with respect to $x$ and $y$, i.e., it solves the subproblem

and updates the Lagrange multiplier $\lambda$ via:

Since minimizing $\mathcal{L}(x,y;\lambda)$ with respect to $x$ and $y$ jointly is usually difficult, while doing so with respect to $x$ and $y$ alternatingly can often be done efficiently, the following alternating direction version of the augmented Lagrangian method (ADAL) is often advocated (see, e.g., ):

If we also update $\lambda$ after we solve the subproblem with respect to $x$, we get the following symmetric version of the ADAL method.

Algorithm (16) has certain theoretical advantages when $f$ and $g$ are smooth. In this case, from the first-order optimality conditions for the two subproblems in (16), we have that:

Substituting these relations into (16), we obtain the following equivalent algorithm for solving (4), which we refer to as the alternating linearization minimization (ALM) algorithm.

Algorithm 1 can be viewed in the following way: at each iteration we construct a quadratic approximation of the function $g(x)$ at the current iterate $y^{k}$ and minimize the sum of this approximation and $f(x)$. The approximation is based on linearizing $g(x)$ (hence the name ALM) and adding a “prox” term $\frac{1}{2\mu}\|x-y^{k}\|_{2}^{2}$. When $\mu$ is small enough ($\mu\leq 1/L(g)$, where $L(g)$ is the Lipschitz constant for $\nabla g$) this quadratic function, $g(y^{k})+\left\langle\nabla g(y^{k}),x-y^{k}\right\rangle+\frac{1}{2\mu}\|x-y^{k}\|_{2}^{2}$ is an upper approximation to $g(x)$, which means that the reduction in the value of $F(x)$ achieved by minimizing $Q_{g}(x,y^{k})$ in Step 1 is not smaller than the reduction achieved in the value of $Q_{g}(x,y^{k})$ itself. Similarly, in Step 2 we build an upper approximation to $f(x)$ at $x^{k+1}$, $f(x^{k+1})+\left\langle\nabla f(x^{k+1}),y-x^{k+1}\right\rangle+\frac{1}{2\mu}\|y-x^{k+1}\|_{2}^{2},$ and minimize the sum $Q_{f}(x^{k+1},y)$ of it and $g(y)$.

Let us now assume that $f(x)$ is in the class $C^{1,1}$ with Lipschitz constant $L(f)$, while $g(x)$ is simply convex. Then from the first-order optimality conditions for the second minimization in (16), we have $-\lambda_{y}^{k+1}\in\partial g(y^{k+1})$, the subdifferential of $g(y)$ at $y=y^{k+1}$. Hence, replacing $\nabla g(y^{k})$ in the definition of $Q_{g}(x,y^{k})$ by $-\lambda_{y}^{k+1}$ in (16), we obtain the following modified version of (16).

Algorithm 2 is identical to the symmetric ADAL algorithm (16) as long as $F(x^{k+1})\leq Q(x^{k+1},y^{k})$ at each iteration (and to Algorithm 1 if $g(x)$ is in $C^{1,1}$ and $\mu\leq 1/\max\{L(f),L(g)\}$). If this condition fails, then the algorithm simply sets $x^{k+1}\leftarrow y^{k}$. Algorithm 2 has the following convergence property and iteration complexity bound. For a proof see the Appendix.

Assume $\nabla f$ is Lipschitz continuous with constant $L(f)$. For $\beta/L(f)\leq\mu\leq 1/L(f)$ where $0<\beta\leq 1$, Algorithm 2 satisfies

where $x^{*}$ is an optimal solution of (4) and $k_{n}$ is the number of iterations until the $k-th$ for which $F(x^{k+1})\leq Q(x^{k+1},y^{k})$. Thus Algorithm 2 produces a sequence which converges to the optimal solution in function value, and the number of iterations needed is $O(1/\epsilon)$ for an $\epsilon$-optimal solution.

If $g(x)$ is also a smooth function in the class $C^{1,1}$ with Lipschitz constant $L(g)\leq 1/\mu$, then Theorem 2.1 also applies to Algorithm 1 since in this case $k_{n}=k$ (i.e., no “skipping” occurs). Note that the iteration complexity bound in Theorem 2.1 can be improved. Nesterov proved that one can obtain an optimal iteration complexity bound of $O(1/\sqrt{\epsilon})$, using only first-order information. His acceleration technique is based on using a linear combination of previous iterates to obtain a point where the approximation is built. This technique has been exploited and extended by Tseng , Beck and Teboulle , Goldfarb et al. and many others. A similar technique can be adopted to derive a fast version of Algorithm 2 that has an improved complexity bound of $O(1/\sqrt{\epsilon})$, while keeping the computational effort in each iteration almost unchanged. However, we do not present this method here, since when applied to the SICS problem, it did not work as well as Algorithm 2.

ALM for SICS

where $f(X)=-\log\det(X)+\langle\hat{\Sigma},X\rangle$ and $g(X)=\rho\|X\|_{1}$, is of the same form as (4). However, in this case neither $f(X)$ nor $g(X)$ have Lipschitz continuous gradients. Moreover, $f(X)$ is only defined for positive definite matrices while $g(X)$ is defined everywhere. These properties of the objective function make the SICS problem especially challenging for optimization methods. Nevertheless, we can still apply (16) to solve the problem directly. Moreover, we can apply Algorithm 2 and obtain the complexity bound in Theorem 2.1 as follows.

The $\log\det(X)$ term in $f(X)$ implicitly requires that $X\in S^{n}_{++}$ and the gradient of $f(X)$, which is given by $-X^{-1}+\hat{\Sigma}$, is not Lipschitz continuous in $S^{n}_{++}$. Fortunately, as proved in Proposition 3.1 in , the optimal solution of (19) $X^{*}\succeq\alpha I$, where $\alpha=\frac{1}{\|\hat{\Sigma}\|+n\rho}.$ Therefore, if we define $\mathcal{C}:=\{X\in S^{n}:X\succeq\frac{\alpha}{2}I\}$, the SICS problem (19) can be formulated as:

We can include constraints $X\in\mathcal{C}$ in Step 1 and $Y\in\mathcal{C}$ in Step 3 of Algorithm 2. Theorem 2.1 can then be applied as discussed in . However, a difficulty now arises when performing the minimization in $Y$. Without the constraint $Y\in\mathcal{C}$, only a matrix shrinkage operation is needed, but with this additional constraint the problem becomes harder to solve. Minimization in $X$ with or without the constraint $X\in\mathcal{C}$ is accomplished by performing an SVD. Hence the constraint can be easily imposed.

Instead of imposing constraint $Y\in\mathcal{C}$ we can obtain feasible solutions by a line search on $\mu$. We know that the constraint $X\succeq\frac{\alpha}{2}I$ is not tight at the solution. Hence if we start the algorithm with $X\succeq\alpha I$ and restrict the step size $\mu$ to be sufficiently small then the iterates of the method will remain in $\mathcal{C}$.

Note however, that the bound on the Lipschitz constant of the gradient of $f(X)$ is $1/\alpha^{2}$ and hence can be very large. It is not practical to restrict $\mu$ in the algorithm to be smaller than $\alpha^{2}$, since $\mu$ determines the step size at each iteration. Hence, for a practical approach we can only claim that the theoretical convergence rate bound holds in only a small neighborhood of the optimal solution. We now present a practical version of our algorithm applied to the SICS problem.

We now show how to solve the two optimization problems in Algorithm 3. The first-order optimality conditions for Step 1 in Algorithm 3, ignoring the constraint $X\in\mathcal{C}$ are:

Consider $V\mbox{Diag}(d)V^{\top}$ - the spectral decomposition of $Y^{k}+\mu_{k+1}(\Lambda^{k}-\hat{\Sigma})$ and let

Since $\nabla f(X)=-X^{-1}+\hat{\Sigma}$, it is easy to verify that $X^{k+1}:=V\mbox{Diag}(\gamma)V^{\top}$ satisfies (21). When the constraint $X\in\mathcal{C}$ is imposed, the optimal solution changes to $X^{k+1}:=V\mbox{Diag}(\gamma)V^{\top}$ with $\gamma_{i}=\max\left\{\alpha/2,\left(d_{i}+\sqrt{d_{i}^{2}+4\mu_{k+1}}\right)/2\right\},i=1,\ldots,n.$ We observe that solving (21) requires approximately the same effort ($O(n^{3})$) as is required to compute $\nabla f(X^{k+1})$. Moreover, from the solution to (21), $\nabla f(X^{k+1})$ is obtained with only a negligible amount of additional effort, since $(X^{k+1})^{-1}:=V\mbox{Diag}(\gamma)^{-1}V^{\top}$.

The first-order optimality conditions for Step 2 in Algorithm 3 are:

Since $g(Y)=\rho\|Y\|_{1}$, it is well known that the solution to (23) is given by

The $O(n^{3})$ complexity of Step 1, which requires a spectral decomposition, dominates the $O(n^{2})$ complexity of Step 2 which requires a simple shrinkage. There is no closed-form solution for the subproblem corresponding to $Y$ when the constraint $Y\in\mathcal{C}$ is imposed. Hence, we neither impose this constraint explicitly nor do so by a line search on $\mu_{k}$, since in practice this degrades the performance of the algorithm substantially. Thus, the resulting iterates $Y^{k}$ may not be positive definite, while the iterates $X^{k}$ remain so. Eventually due to the convergence of $Y^{k}$ and $X^{k}$, the $Y^{k}$ iterates become positive definite and the constraint $Y\in\mathcal{C}$ is satisfied.

Let us now remark on the learning based intuition behind Algorithm 3. We recall that $-\Lambda^{k}\in\partial g(Y^{k})$. The two steps of the algorithm can be written as

The SICS problem is trying to optimize two conflicting objectives: on the one hand it tries to find a covariance matrix $X^{-1}$ that best fits the observed data, i.e., is as close to $\hat{\Sigma}$ as possible, and on the other hand it tries to obtain a sparse matrix $X$. The proposed algorithm address these two objectives in an alternating manner. Given an initial “guess” of the sparse matrix $Y^{k}$ we update this guess by a subgradient descent step of length $\mu_{k+1}$: $Y^{k}+\mu_{k+1}\Lambda^{k}$. Recall that $-\Lambda^{k}\in\partial g(Y^{k})$. Then problem (24) seeks a solution $X$ that optimizes the first objective (best fit of the data) while adding a regularization term which imposes a Gaussian prior on $X$ whose mean is the current guess for the sparse matrix: $Y^{k}+\mu_{k+1}\Lambda^{k}$. The solution to (24) gives us a guess for the inverse covariance $X^{k+1}$. We again update it by taking a gradient descent step: $X^{k+1}-\mu_{k+1}(\hat{\Sigma}-(X^{k+1})^{-1})$. Then problem (25) seeks a sparse solution $Y$ while also imposing a Gaussian prior on $Y$ whose mean is the guess for the inverse covariance matrix $X^{k+1}-\mu_{k+1}(\hat{\Sigma}-(X^{k+1})^{-1})$. Hence the sequence of $X^{k}$’s is a sequence of positive definite inverse covariance matrices that converge to a sparse matrix, while the sequence of $Y^{k}$’s is a sequence of sparse matrices that converges to a positive definite inverse covariance matrix.

An important question is how to pick $\mu_{k+1}$. Theory tells us that if we pick a small enough value, then we can obtain the complexity bounds. However, in practice this value is too small. We discuss the simple strategy that we use in the next section.

Numerical Experiments

In this section, we present numerical results on both synthetic and real data to demonstrate the efficiency of our SICS ALM algorithm. Our codes for ALM were written in MATLAB. All numerical experiments were run in MATLAB 7.3.0 on a Dell Precision 670 workstation with an Intel Xeon(TM) 3.4GHZ CPU and 6GB of RAM.

Since $-\Lambda^{k}\in\partial g(Y^{k})$, $\|\Lambda^{k}\|_{\infty}\leq\rho$; hence $\hat{\Sigma}-\Lambda^{k}$ is a feasible solution to the dual problem (2) as long as it is positive definite. Thus the duality gap at the $k$-th iteration is given by:

We define the relative duality gap as: $Rel.gap:=Dgap/(1+|pobj|+|dobj|),$ where $pobj$ and $dobj$ are respectively the objective function values of the primal problem (19) at point $X^{k}$, and the dual problem (2) at $\hat{\Sigma}-\Lambda^{k}$. Defining $d_{k}(\phi(x))\equiv\max\{1,\phi(x^{k}),\phi(x^{k-1})\}$, we measure the relative changes of objective function value $F(X)$ and the iterates $X$ and $Y$ as follows:

Note that in (26), computing $\log\det(X^{k})$ is easy since the spectral decomposition of $X^{k}$ is already available (see (21) and (22)), but computing $\log\det(\hat{\Sigma}-\Lambda^{k})$ requires another expensive spectral decomposition. Thus, in practice, we only check (27)(i) every $N_{gap}$ iterations. We check (27)(ii) at every iteration since this is inexpensive.

A continuation strategy for updating $\mu$ is also crucial to ALM. In our experiments, we adopted the following update rule. After every $N_{\mu}$ iterations, we set $\mu:=\max\{\mu\cdot\eta_{\mu},\bar{\mu}\}$; i.e., we simply reduce $\mu$ by a constant factor $\eta_{\mu}$ every $N_{\mu}$ iterations until a desired lower bound on $\mu$ is achieved.

We compare ALM (i.e., Algorithm 3 with the above stopping criteria and $\mu$ updates), with the projected subgradient method (PSM) proposed by Duchi et al. in and implemented by Mark Schmidt The MATLAB can be downloaded from http://www.cs.ubc.ca/$\sim$schmidtm/Software/PQN.html and the smoothing method (VSM) The MATLAB code can be downloaded from http://www.math.sfu.ca/$\sim$zhaosong proposed by Lu in , which are considered to be the state-of-the-art algorithms for solving SICS problems. The per-iteration complexity of all three algorithms is roughly the same; hence a comparison of the number of iterations is meaningful. The parameters used in PSM and VSM are set at their default values. We used the following parameter values in ALM: $\epsilon_{gap}=10^{-3},\epsilon_{rel}=10^{-8},N_{gap}=20,N_{\mu}=20,\bar{\mu}=\max\{\mu_{0}\eta_{\mu}^{8},10^{-6}\},\eta_{\mu}=1/3,$ where $\mu_{0}$ is the initial $\mu$ which is set according to $\rho$; specifically, in our experiments, $\mu_{0}=100/\rho,$ if $\rho<0.5$, $\mu_{0}=\rho$ if $0.5\leq\rho\leq 10$, and $\mu_{0}=\rho/100$ if $\rho>10$.

From Table 1 we see that on these randomly created SICS problems, ALM outperforms PSM and VSM in both accuracy and CPU time with the performance gap increasing as $\rho$ increases. For example, for $\rho=1.0$ and $n=2000$, ALM achieves $Dgap=9.58e-4$ in about 1 hour and 15 minutes, while PSM and VSM need about 3 hours and 25 minutes and 10 hours and 23 minutes, respectively, to achieve similar accuracy.

2 Experiments on real data

We tested ALM on real data from gene expression networks using the five data sets from provided to us by Kim-Chuan Toh: (1) Lymph node status; (2) Estrogen receptor; (3) Arabidopsis thaliana; (4) Leukemia; (5) Hereditary breast cancer. See and references therein for the descriptions of these data sets. Table 2 presents our test results. As suggested in , we set $\rho=0.5$. From Table 2 we see that ALM is much faster and provided more accurate solutions than PSM and VSM.

3 Solution Sparsity

In this section, we compare the sparsity patterns of the solutions produced by ALM, PSM and VSM. For ALM, the sparsity of the solution is given by the sparsity of $Y$. Since PSM and VSM solve the dual problem, the primal solution $X$, obtained by inverting the dual solution $W$, is never sparse due to floating point errors. Thus it is not fair to measure the sparsity of $X$ or a truncated version of $X$. Instead, we measure the sparsity of solutions produced by PSM and VSM by appealing to complementary slackness. Specifically, the $(i,j)$-th element of the inverse covariance matrix is deemed to be nonzero if and only if $|W_{ij}-\hat{\Sigma}_{ij}|=\rho$. We give results for a random problem ($n=500$) and the first real data set in Table 3. For each value of $\rho$, the first three rows show the number of nonzeros in the solution and the last three rows show the number of entries that are nonzero in the solution produced by one of the methods but are zero in the solution produced by the other method. The sparsity of the ground truth inverse covariance matrix of the synthetic data is 6.76%.

From Table 3 we can see that when $\rho$ is relatively large ($\rho\geq 0.5$), all three algorithms produce solutions with exactly the same sparsity patterns. Only when $\rho$ is very small, are there slight differences. We note that the ROC curves depicting the trade-off between the number of true positive elements recovered versus the number of false positive elements as a function of the regularization parameter $\rho$ are also almost identical for the three methods.

Acknowledgements

We would like to thank Professor Kim-Chuan Toh for providing the data set used in Section 4.2. The research reported here was supported in part by NSF Grants DMS 06-06712 and DMS 10-16571, ONR Grant N00014-08-1-1118 and DOE Grant DE-FG02-08ER25856.

References

Appendix

where $\gamma_{\psi}(v)$ is any subgradient in the subdifferential $\partial\psi(v)$ of $\psi(v)$ at the point $v$, and

Let $\Phi(\cdot)=\phi(\cdot)+\psi(\cdot)$. For any $v$, if

Now since $\phi$ and $\psi$ are convex we have

where $\gamma_{\phi}(\cdot)$ is a subgradient of $\phi(\cdot)$ and $\gamma_{\phi}(p_{\psi}(v))$ satisfies the first-order optimality conditions for (28), i.e.,

Therefore, from (33), (36) and (37) it follows that

Let $I$ be the set of all iteration indices until $k-1$-st for which no skipping occurs and let $I_{c}$ be its complement. Let $I=\{n_{i}\},\ i=0,\ldots,k_{n}-1$. It follows that for all $n\in I_{c}$ $x^{n+1}=y^{n}$.

For $n\in I$ we can apply Lemma 5.1 to obtain the following inequalities. In (30), by letting $\psi=f$, $\phi=g$, $u=x^{*}$ and $u=x^{n+1}$, we get $p_{\psi}(v)=y^{n+1}$, $\Phi=F$ and

Similarly, by letting $\psi=g$, $\phi=f$, $u=x^{*}$ and $v=y^{n}$ in (30) we get $p_{g}(v)=x^{n+1}$, $\Phi=F$ and

Taking the summation of (38) and (39) we get

due to the fact that $x^{n+1}=y^{n}$ in this case.

Summing (40) and (41) over $n=0,1,\ldots,k-1$ we get

For any $n$, since Lemma 5.1 holds for any $u$, letting $u=x^{n+1}$ instead of $x^{*}$ we get from (38) that

Similarly, for $n\in I$ by letting $u=y^{n}$ instead of $x^{*}$ we get from (39) that

On the other hand, for $n\in I_{c}$, (45) also holds because $x^{n+1}=y^{n}$, and hence holds for all $n$.

(46) and (47) show that the sequences $F(y^{n})$ and $F(x^{n})$ are non-increasing. Thus we have,

From (44) we know that $F(x^{k})\geq F(y^{k})$. Thus (49) implies that

Also, for any given $\epsilon>0$, as long as $k\geq\frac{L(f)\|x^{0}-x^{*}\|^{2}}{2\beta\epsilon}$, we have from (18) that $F(y^{k})-F(x^{*})\leq\epsilon$; i.e., the number of iterations needed is $O(1/\epsilon)$ for an $\epsilon$-optimal solution. ∎