Alternating Direction Methods for Latent Variable Gaussian Graphical Model Selection

Introduction

In this paper, we consider alternating direction methods with the theoretical guarantee of global convergence for computing the latent-variable graphical model selection . Graphical model selection is closely related to the inverse covariance matrix estimation problem, which is of fundamental importance in multivariate statistical inference. In particular, when data $X=(X_{1},\cdots,X_{p})^{\prime}$ follow a $p$-dimensional joint normal distribution with some unknown variance matrix $\Sigma$, the precision matrix $\Theta=\Sigma^{-1}$ can be directly translated into a Gaussian graphical model. The zero entries in the precision matrix \Theta=\bigl{(}\theta_{ij}\bigr{)}_{1\leq i,j\leq p} precisely capture the desired conditional independencies in the Gaussian graphical model , i.e. $\theta_{ij}=0$ if and only if X_{i}\mathchoice{\mathrel{\hbox to0.0pt{\displaystyle\perp\hss}\mkern 2.0mu{\displaystyle\perp}}}{\mathrel{\hbox to0.0pt{\textstyle\perp\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptstyle\perp\hss}\mkern 2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptscriptstyle\perp\hss}\mkern 2.0mu{\scriptscriptstyle\perp}}}X_{j}|~{}X_{-(i,j)}. The Gaussian graphical model has been successfully used to explore complex systems consisting of Gaussian random variables in many research fields, including gene expression genomics , image processing , macroeconomics determinants study , and social study .

The aforementioned Gaussian graphical model selection methods were proposed under the ideal setting without missing variables. The recent paper by considered a more realistic scenario where the full data consist of both observed variables and missing (hidden) variables. Let $X_{p\times 1}$ be the observed variables. Suppose that there are some hidden variables $Y_{r\times 1}$ ($r\ll p$) such that $(X,Y)$ jointly follow a multivariate normal distribution. Denote the covariance matrix by $\Sigma_{(X,Y)}$ and the precision matrix by $\Theta_{(X,Y)}$. Then we can write $\Sigma_{(X,Y)}=[\Sigma_{X},\Sigma_{XY};\Sigma_{YX},\Sigma_{Y}]$ and $\Theta_{(X,Y)}=[\Theta_{X},\Theta_{XY};\Theta_{YX},\Theta_{Y}]$. Given the hidden variables $Y$, the conditional concentration matrix of observed variables, $\Theta_{X}$, is sparse for a sparse graphical model. However, the marginal concentration matrix of observed variables, $\Sigma_{X}^{-1}=\Theta_{X}-\Theta_{XY}\Theta_{Y}^{-1}\Theta_{YX}$, might not be a sparse matrix but a difference between the sparse term $\Theta_{X}$ and the low-rank term $\Theta_{XY}\Theta_{Y}^{-1}\Theta_{YX}$. The problem of interest is to recover the sparse conditional matrix $\Theta_{X}$ based on observed variables $X$. accomplished this goal by solving a convex optimization problem under the assumption that $\Sigma_{X}^{-1}=S-L$ for some sparse matrix $S$ and low-rank matrix $L$. The low rank assumption on $L$ holds naturally since $r$ is much less than $p$. Motivated by the success of the convex relaxation for rank-minimization problem, introduced a regularized maximum normal likelihood decomposition framework called the latent variable graphical model selection (LVGLASSO) as follows.

where $\hat{\Sigma}_{X}$ is the sample covariance matrix of $X$ and $\mathbf{Tr}(L)$ denotes the trace of matrix $L$. In the high-dimensional setting, established the consistency theory for (1.2) concerning its recovery of the support and sign pattern of $S$ and the rank of $L$.

Solving the convex optimization problem (1.2) is very challenging, especially for large problems. considered (1.2) as a log-determinant semidefinite programming (SDP) problem, and used a Newton-CG based proximal point algorithm (LogdetPPA) proposed by to solve it. However, LogdetPPA does not take advantage of the special structure of the problem, and we argue that it is inefficient for solving large-scale problems. To illustrate our point, let us consider the special case of (1.2) with $L=0$, and then the latent variable graphical model selection (1.2) exactly reduces to the Gaussian graphical model selection (1.1). Note that (1.1) can be rewritten as

where $\|W\|_{\infty}$ is the largest absolute value of the entries of $U$. The dual problem of (1.1) can be obtained by exchanging the order of max and min, i.e.,

Both the primal and the dual graphical Lasso problems (1.1) and (1.3) can be viewed as semidefinite programming problems and can be solved via interior point methods (IPMs) in polynomial time . However, the per-iteration computational cost and memory requirements of an IPM are prohibitively high for (1.1) and (1.3), especially when the size of the matrix is large. Customized SDP based methods such as the ones studied in and require a reformulation of the problem that increases the size of the problem and thus makes them impractical for solving large-scale problems. Therefore, most of the methods developed for solving (1.1) and (1.3) are first-order methods. These methods include block coordinate descent type methods , projected gradient method and variants of Nesterov’s accelerated method . Recently, alternating direction methods have been applied to solve (1.1) and shown to be very effective .

In this paper, we propose two alternating direction type methods to solve the latent variable graphical model selection. The first method is to apply the alternating direction method of multipliers to solve this problem. This is due to the fact that the latent variable graphical model selection can be seen as a special case of the consensus problem discussed in . The second method we propose is an alternating direction method with proximal gradient steps. To apply the second method, we first group the variables into two blocks and then apply the alternating direction method with one of the subproblems being solved inexactly by taking a proximal gradient step. Our methods exploit and take advantage of the special structure of the problem and thus can solve large problems very efficiently. Although the convergence results of the proposed methods are not very different from the existing results for alternating direction type methods, we still include the convergence proof for the second method in the appendix for completeness. We apply the proposed methods to solving problems from both synthetic data and gene expression data and show that our method outperform the state-of-the-art Newton-CG proximal point algorithm LogdetPPA significantly on both accuracy and CPU times.

The rest of this paper is organized as follows. In Section 2, we give some preliminaries on alternating direction method of multipliers and proximal mappings. In Section 3, we propose solving LVGLASSO (1.2) as a consensus problem using the classical alternating direction method of multipliers. We propose the proximal gradient based alternating direction method for solving (1.2) in Section 4. In Section 5, we apply our alternating direction method to solving (1.2) using both synthetic data and gene expression data. We draw some conclusions in Section 6.

Preliminaries

Problem (1.2) can be rewritten in the following equivalent form by introducing a new variable $R$:

where the indicator function $\mathcal{I}(L\succeq 0)$ is defined as

Note that we have dropped the constraint $R\succ 0$ since it is already implicitly imposed by the $\log\det R$ function.

Now since the objective function involves three separable convex functions and the constraint is simply linear, Problem (2.2) is suitable for alternating direction method of multipliers (ADMM). ADMM is closely related to the Douglas-Rachford and Peaceman-Rachford operator-splitting methods for finding zero of the sum of two monotone operators that have been studied extensively in . ADMM has been revisited recently due to its success in the emerging applications of structured convex optimization problems arising from image processing, compressed sensing, machine learning, semidefinite programming and statistics etc. (see e.g., ).

Problem (2.2) is suitable for alternating direction methods because the three convex functions involved in the objective function, i.e.,

The proximal mapping of $f(R)$ defined in (2.4) is

The first-order optimality conditions of (2.8) are given by

satisfies (2.9) and thus gives the optimal solution of (2.8), where $U\mathop{\operator@font diag}\nolimits(\sigma)U^{\top}$ is the eigenvalue decomposition of matrix $\xi\hat{\Sigma}_{X}-Z$ and

Note that (2.11) guarantees that the solution of (2.8) given by (2.10) is a positive definite matrix. The proximal mapping of $g(S)$ defined in (2.5) is

The proximal mapping of $h(L)$ defined in (2.6) is

It is easy to verify that the solution of (2.14) is given by

where $Z=U\mathop{\operator@font diag}\nolimits(\sigma)U^{\top}$ is the eigenvalue decomposition of $Z$ and $\gamma$ is given by

Note that (2.16) guarantees that $L$ given in (2.15) is a positive semidefinite matrix.

The discussions above suggest the following natural ADMM for solving (2.2) be efficient.

where the augmented Lagrangian function is defined as

$\Lambda$ is the Lagrange multiplier and $\mu>0$ is the penalty parameter. Note that the three subproblems in (2.17) correspond to the proximal mappings of $f$, $g$ and $h$ defined in (2.4), (2.5) and (2.6), respectively. Thus they are all easy to solve. However, the global convergence of ADMM (2.17) with three blocks of variables was ambiguous. Only until very recently, was it shown that (2.17) globally converges under certain conditions (see ). It should be noted, however, that the error bound condition required in is strong and only a few classes of convex function are known that satisfy this condition.

ADMM for Solving (2.2) as a Consensus Problem

Problem (2.2) can be rewritten as a convex minimization problem with two blocks of variables and two separable functions as follows:

with $f,g$ and $h$ defined in (2.4), (2.5) and (2.6), respectively. The ADMM applied to solving (3.1) can be described as follows:

The first-order optimality conditions of (3.3) are given by

where $\Gamma$ is the Lagrange multiplier associated with (3.3). Thus we get,

Substituting them into the equality constraint in (3.3), we get

By substituting (3.5) into (3.4) we get the solution to (3.3).

The ADMM (3.2) solves Problem (3.1) with two blocks of variables. It can be seen as a special case of the consensus problem discussed in . The global convergence result of (3.2) has also been well studied in the literature (see e.g., ).

A Proximal Gradient based Alternating Direction Method

In this section, we propose another alternating direction type method to solve (2.2). In Section 3, we managed to reduce the original problem with three blocks of variables (2.2) to a new problem with two blocks of variables (3.1). As a result, we can use ADMM for solving problems with two blocks of variables, whose convergence has been well studied. Another way to reduce the problem (2.2) into a problem with two blocks of variables is to group two variables (say $S$ and $L$) as one variable. This leads to the new equivalent form of (2.2):

where $W=[S;L]$ and $\varphi(W)=g(S)+h(L)$. Now the ADMM for solving (3.1) can be described as

where $\Lambda$ is the Lagrange multiplier associated with the equality constraint and $\mu>0$ is a penalty parameter. The first subproblem in (4.2) is still easy and it corresponds to the proximal mapping of function $f$. However, the second subproblem in (4.2) is not easy, because the two parts of $W$ are coupled together in the quadratic penalty term. To overcome this difficulty, we solve the second subproblem in (4.2) inexactly by one step of a proximal gradient method. Note that the second subproblem in (4.2) can be reduced to

One step of proximal gradient method solves the following problem

Since the two parts of $W=[S;L]$ are separable in the quadratic part now, (4.4) reduces to two problems

where $G_{R}^{k}=R^{k+1}-S^{k}+L^{k}-\mu\Lambda^{k}$. Both (4.5) and (4.6) are easy to solve as they correspond to the proximal mappings of functions $g$ and $h$, respectively. Thus, our proximal gradient based alternating direction method (PGADM) can be summarized as

The idea of incorporating proximal step into the alternating direction method of multipliers has been suggested by and . This idea has then been generalized by to allow varying penalty and proximal parameters. Recently, this technique has been used for sparse and low-rank optimization problems (see and ). More recently, some convergence properties of alternating direction methods with proximal gradient steps have been studied by , , and . However, for the seek of completeness, we include a global convergence proof for Algorithm 1 in the Appendix.

In Algorithm 1, we grouped $S$ and $L$ as one block of variable. We also implemented the other two ways of grouping the variables, i.e., group $R$ and $S$ as one block, and group $R$ and $L$ as one block. We found from the numerical experiments that these two alternatives yielded similar practical performance as Algorithm 1.

Numerical experiments

In this section, we present numerical results on both synthetic and real data to demonstrate the efficiency of the proposed methods: ADMM (3.2) and PGADM (Algorithm 1). Our codes were written in MATLAB. All numerical experiments were run in MATLAB 7.12.0 on a laptop with Intel Core I5 2.5 GHz CPU and 4GB of RAM.

We first compared ADMM (3.2) with PGADM (Algorithm 1) on some synthetic problems. We compared ADMM and PGADM using two different ways of choosing $\mu$. One set of comparisons used a fixed $\mu=10$, and the other set of comparisons used a continuation scheme to dynamically change $\mu$. The continuation scheme we used was to set the initial value of $\mu$ as the size of the matrix $p$, and then multiply $\mu$ by $1/4$ after every 10 iterations.

We then compared the performance of PGADM (with continuation on $\mu$) with LogdetPPA proposed by and used in for solving (2.2).

We observed form the numerical experiments that the step size $\tau$ of the proximal gradient step in PGADM (Algorithm 1) can be slightly larger than $1/2$ and the algorithm produced very good results. We thus chose the step size $\tau$ to be $0.6$ in our experiments.

We computed the relative infeasibility of the sequence $(R^{k},S^{k},L^{k})$ generated by inexact ADMM using

In the comparison of ADMM and PGADM, the size of all problems was chosen as $p=1000$. For fixed $\mu=10$, we first ran the ADMM for 100 iterations, and recorded the objective function value and $infeas$. We then ran PGADM until it achieves an objective function value within relative error $10^{-5}$ compared with the objective function value given by ADMM, or it achieves an $infeas$ within relative error $10^{-5}$ compared with the $infeas$ given by ADMM. The number of iterations, CPU times, $infeas$ and objective function values for both ADMM and PGADM were reported in Table 1. From Table 1, we see that for fixed $\mu=10$, ADMM was faster than PGADM when $\alpha$ and $\beta$ are both small, and PGADM was faster than ADMM when $\alpha$ and $\beta$ are both large.

We then further compare ADMM and PGADM on synthetic data with the continuation scheme for $\mu$ discussed above. We terminated both ADMM and PGADM when $infeas<10^{-5}$. We reported the results in Table 2.

From Table 2, we see that the continuation scheme used really helped to speed up the convergence and produced much better results. Also, using this continuation scheme, PGADM was faster than ADMM with comparable residuals and objective function values. However, we should remark that PGADM was faster than ADMM using the specific continuation scheme. If other continuation schemes were adopted, the results could be quite different. In the comparison with LogdetPPA in the following sections, we only compare LogdetPPA with PGADM with this continuation scheme.

2 Comparison of PGADM and LogdetPPA on Synthetic Data

In this section, we compare PGADM with LogdetPPA on synthetic data created the same way as in the last section. LogdetPPA, proposed by Wang et al.in , is a proximal point algorithm for solving semidefinite programming problems with $\log\det(\cdot)$ function. The specialized MATLAB codes of LogdetPPA for solving (2.2) were downloaded from http://ssg.mit.edu/$\sim$venkatc/latent-variable-code.html.

We compared PGADM (with continuation on $\mu$) with LogdetPPA with different $\alpha$ and $\beta$. We reported the comparison results on objective function value, CPU time, sparsity of $S$ and infeas in Table 3. The sparsity of $S$ is denoted as

i.e., the percentage of nonzero entries. Since matrix $S$ generated by LogdetPPA is always dense but with many small entries, we also measure its sparsity by truncating small entries that less than $10^{-4}$ to zeros, i.e.,

All CPU times reported are in seconds. We report the speed up of PGADM over LogdetPPA in Table 4.

3 Comparison of PGADM and LogdetPPA on Gene expression data

To further demonstrate the efficacy of PGADM, we applied PGADM to solving (2.2) with two gene expression data sets. One data set is the Rosetta Inpharmatics Compendium of gene expression data (denoted as Rosetta) profiles which contains $301$ samples with $6316$ variables (genes). The other data set is the Iconix microarray data set (denoted as Iconix) from drug treated rat livers which contains $255$ samples with $10455$ variables.

For a given number of observed variables $p$, we created the sample covariance matrix $\Sigma_{X}$ by the following procedure. We first computed the variances of all of variables using all the sample data. We then selected the $p$ variables with the highest variances and computed the sample covariance matrix $\Sigma_{X}$ of these $p$ variables using all the sample data. We reported the comparison results of PGADM and LogdetPPA in Tables 5 and 6 for the Rosetta and Iconix data sets, respectively. Table 7 summarizes the speed up of PGADM over LogdetPPA.

From Table 5, we again see that PGADM always generates solutions with comparable objective function values in much less time. For example, for $p=2000$, LogdetPPA needs 1 hour 16 minutes to solve it while PGADM takes just 7 minutes. From Table 6, we see that the advantage of PGADM is more obvious. For $p=200,500,1000$ and $2000$, PGADM always generates solutions with much smaller objective function values and it is always much faster than LogdetPPA. For example, for $p=2000$, LogdetPPA takes 3 hours 40 minutes to solve it while PGADM just takes about 12 minutes.

Conclusion

In this paper, we proposed alternating direction methods for solving latent variable Gaussian graphical model selection. The global convergence results of our methods were established. We applied the proposed methods for solving large problems from both synthetic data and gene expression data. The numerical results indicated that our methods were five to thirty five times faster than a state-of-the-art Newton-CG proximal point algorithm.

Acknowledgement

The authors thank Professor Stephen Boyd for suggesting solving (2.2) as a consensus problem (3.1) using ADMM. Shiqian Ma’s research is supported by the National Science Foundation postdoctoral fellowship through the Institute for Mathematics and Its Applications at University of Minnesota. Hui Zou’s research is supported in part by grants from the National Science Foundation and the Office of Naval Research.

References

Appendix A Global Convergence Analysis of PGADM

In this section, we establish the global convergence result of PGADM (Algorithm 1). This convergence proof is not much different with the one given by for compressed sensing problems. We include the proof here just for completeness.

We introduce some notation first. We define $W=\begin{pmatrix}S\\ L\end{pmatrix}$. We define functions $F(\cdot)$ and $G(\cdot)$ as

and our PGADM (Algorithm 1) can be rewritten as

Before we prove the global convergence result, we need to prove the following lemma.

Assume that $(R^{*},W^{*})$ is an optimal solution of (A.1) and $\Lambda^{*}$ is the corresponding optimal dual variable associated with the equality constraint $R+AW=0$. Assume the step size $\tau$ of the proximal gradient step satisfies $0<\tau<1/2$. Then there exists $\eta>0$ such that the sequence $(R^{k},W^{k},\Lambda^{k})$ produced by (A.2) satisfies

where $U^{*}=\begin{pmatrix}W^{*}\\ \Lambda^{*}\end{pmatrix}$, $U^{k}=\begin{pmatrix}W^{k}\\ \Lambda^{k}\end{pmatrix}$ and $H=\begin{pmatrix}\frac{1}{\mu\tau}I_{p\times p}&0\\ 0&\mu I_{p\times p}\end{pmatrix}$, and the norm $\|\cdot\|_{H}^{2}$ is defined as $\|U\|_{H}^{2}=\langle U,HU\rangle$ and the corresponding inner product $\langle\cdot,\cdot\rangle_{H}$ is defined as $\langle U,V\rangle_{H}=\langle U,HV\rangle$.

Since $(R^{*},W^{*},\Lambda^{*})$ is optimal to (A.1), it follows from the KKT conditions that the followings hold:

Note that the first-order optimality conditions for the first subproblem (i.e., the subproblem with respect to $R$) in (A.2) are given by

By using the updating formula for $\Lambda^{k}$, i.e.,

Combining (A.4) and (A.9) and using the fact that $\partial F(\cdot)$ is a monotone operator, we get

The first-order optimality conditions for the second subproblem (i.e., the subproblem with respect to $W$) in (A.2) are given by

Combining (A.5) and (A.12) and using the fact that $\partial G(\cdot)$ is a monotone operator, we get

Summing (A.10) and (A.13), and using $R^{*}=-AW^{*}$ and $R^{k+1}=\mu(\Lambda^{k}-\Lambda^{k+1})-AW^{k+1}$, we obtain,

Using the notation of $U^{k}$, $U^{*}$ and $H$, (A.14) can be rewritten as

Let $\xi:=\tau+1/2$, then we know that $2\tau<\xi<1$ since $0<\tau<1/2$. Let $\rho:=\mu\xi$. Then from Cauchy-Schwartz inequality we have

where the $\lambda_{\max}(A^{\top}A)$ denotes the largest eigenvalue of matrix $A^{\top}A$ and the equality is due to the fact that $\lambda_{\max}(A^{\top}A)=2$. Combining (A.17) and (A.18) we get

where $\eta:=\min\{1-\frac{2\mu\tau}{\rho},1-\frac{\rho}{\mu}\}=\min\{1-\frac{2\tau}{\xi},1-\xi\}>0$. This completes the proof. ∎

We are now ready to give the main convergence result of Algorithm (A.2).

The sequence $\{(R^{k},W^{k},\Lambda^{k})\}$ produced by Algorithm (A.2) from any starting point converges to an optimal solution to Problem (A.1).

(i) $\|U^{k}-U^{k+1}\|_{H}\rightarrow 0$;

(ii) $\{U^{k}\}$ lies in a compact region;

(iii) $\|U^{k}-U^{*}\|_{H}^{2}$ is monotonically non-increasing and thus converges.

It follows from (i) that $\Lambda^{k}-\Lambda^{k+1}\rightarrow 0$ and $W^{k}-W^{k+1}\rightarrow 0$. Then (A.8) implies that $R^{k}-R^{k+1}\rightarrow 0$ and $R^{k}+AW^{k}\rightarrow 0$. From (ii) we obtain that, $U^{k}$ has a subsequence $\{U^{k_{j}}\}$ that converges to $\hat{U}=(\hat{W},\hat{\Lambda})$, i.e., $\Lambda^{k_{j}}\rightarrow\hat{\Lambda}$ and $W^{k_{j}}\rightarrow\hat{W}$. From $R^{k}+AW^{k}\rightarrow 0$ we also get that $R^{k_{j}}\rightarrow\hat{R}:=-A\hat{W}$. Therefore, $(\hat{R},\hat{W},\hat{\Lambda})$ is a limit point of $\{(R^{k},W^{k},\Lambda^{k})\}$.

(A.20), (A.21) and $\hat{R}+A\hat{W}=0$ imply that $(\hat{R},\hat{W},\hat{\Lambda})$ satisfies the KKT conditions for (A.1) and thus is an optimal solution to (A.1). Therefore, we showed that any limit point of $\{(R^{k},W^{k},\Lambda^{k})\}$ is an optimal solution to (A.1).

To complete the proof, we need to show that the limit point is unique. Let $\{(\hat{R}_{1},\hat{W}_{1},\hat{\Lambda}_{1})\}$ and $\{(\hat{R}_{2},\hat{W}_{2},\hat{\Lambda}_{2})\}$ be any two limit points of $\{(R^{k},W^{k},\Lambda^{k})\}$. As we have shown, both $\{(\hat{R}_{1},\hat{W}_{1},\hat{\Lambda}_{1})\}$ and $\{(\hat{R}_{2},\hat{W}_{2},\hat{\Lambda}_{2})\}$ are optimal solutions to (A.1). Thus, $U^{*}$ in (A.19) can be replaced by $\hat{U}_{1}:=(\hat{R}_{1},\hat{W}_{1},\hat{\Lambda}_{1})$ and $\hat{U}_{2}:=(\hat{R}_{2},\hat{W}_{2},\hat{\Lambda}_{2})$. This results in

and we thus get the existence of the limits

Thus we must have $\|\hat{U}_{1}-\hat{U}_{2}\|_{H}^{2}=0$ and hence the limit point of $\{(R^{k},W^{k},\Lambda^{k})\}$ is unique. ∎

We now immediately have the global convergence result for Algorithm 1 for solving Problem (2.2).

The sequence $\{(R^{k},S^{k},L^{k},\Lambda^{k})\}$ produced by Algorithm 1 from any starting point converges to an optimal solution to Problem (2.2).