1 Introduction
Estimating covariance matrices in the highdimensional setting arises in many applications and has drawn considerable interest recently. Because the sample covariance matrix is typically poorly behaved in the highdimensional regime, regularizing the sample covariance based on some assumptions of the underlying true covariance is often essential to gain robustness and stability of the estimation.
One form of regularization that has gained popularity recently is to require the the underlying inverse covariance matrix to be sparse [1, 2, 3, 4]
. If the data follow a multivariate Gaussian distribution with covariance matrix
, the entries of the inverse covariance matrix (also known as concentration matrix or precision matrix) encode the information of conditional dependencies between variables: if the variables and are conditionally independent given all others. Therefore, the sparsity regularization is equivalent to the assumption that most of the variable pairs in the highdimensional setting are conditionally independent.To make the estimation problem computational tractable, one often adopts a convex relaxation of the sparsity constraint and uses the norm to promote the sparsity of the concentration matrix [3, 4, 5, 6]. Denote the empirical covariance. Under the maximum likelihood framework, the covariance matrix estimation problem is then formulated as solving the following optimization problem:
(1) 
where tr denotes the trace, is a sparsity regularization parameter, and denotes that is positive semidefinite. Due to the penalty term and the explicit positive definite constraint on , the method leads to a sparse estimation of the concentration matrix that is guaranteed to be positive definite. The problem is convex and many algorithms have been proposed to solve the problem efficiently in high dimension [4, 7, 8, 9].
However, in many real applications only a subset of the variables are directly observed, and no additional information is provided on both the number of the latent variables and their relationship with the observed ones. For instance, in the area of functional genomics it is often the case that only mRNAs of the genes can be directly measured, but not the proteins, which are correlated but have no direct correspondence to the mRNAs because of the prominent role of the postranscriptional regulation. Another example is the movie recommender system where the preference of a movie can be strongly influenced by latent factors such as advertisements, social environment, etc. In these and other cases, the observed variables can be densely correlated because of the marginalization over the unobserved hidden variables. Therefore, the sparsity regularization alone may fail to achieve the desire results.
We consider the setting in which the hidden () and the observed variables () are jointly Gaussian with covariance matrix . The marginal statistics of the observed variable are given by the marginal covariance matrix , which is simply a submatrix of the full covariance matrix . Let the concentration matrix . The marginal concentration matrix corresponding to the observed variables is given by the Schur complement [10]:
(2) 
where , and are the corresponding submatrices of the full concentration matrix. Based on the Schur complement, it is clear that the marginal concentration matrix of the observed variables can be decomposed into two components: one is , which specifies the conditional dependencies of the observed variables given both the observed and latent variables, and the other is , which represents the effect of marginalization over the hidden variables. One can now impose assumptions to the two underlying components separately.
By assuming that the matrix is sparse and the number of latent variables is small, the maximum likelihood estimation of the covariance matrix of the observed variables at the presence of latent variables can then be formulated as
(3) 
where we decompose with denoting and denoting . Because the number of the hidden variables is small, is of low rank, whose convex relaxation is the trace norm. There are two regularization parameters in this model: regularizes the sparsity of , and regularizes the rank of . Under certain regularity conditions, Chandrasekaran et al. showed that this model can consistently estimate the underlying model structure in the highdimensional regime in which the number of observed/hidden variables grow with the number of samples of the observed variables [10].
The objective function in (1) is strictly convex, so a global optimal solution is guaranteed to exist and be unique. Finding the optimal solution in the highdimension setting is computationally challenging due to the term appeared in the likelihood function, the trace norm, the nondifferentiability of the penalty, and the positive semidefinite constraints. For largescale problems, the stateoftheart algorithm for solving (1) is to use the special purpose solver LogdetPPA [11] developed for logdeterminant semidefinite programs. However, the solver LogdetPPA is designed to solve smooth problems. In order to use LogdetPPA, one has to reformulate (1) to a smooth problem. As a result, no optimal sparse matrix
can be generated and additional heuristic steps involving thresholding have to be applied in order to produce sparsity. In addition, LogdetPPA is not especially designed for (
1). We believe much more efficient algorithms can be generated by considering the unique structures of the model specifically.The main contribution of this paper contains two aspects. First, we present a new algorithm for solving (1) and show that the algorithm is significantly faster than the stateoftheart method, especially for largescale problems. The algorithm is derived by reformulating the problem and adapting the split Bregman method [8, 9]. We derive closed form solutions for each subproblem involved in the split Bregman iterations. Second, we apply the method to analyze a largescale gene expression data, and find that the model with latent variables explain the data much better than the one without assuming latent variables. In addition, we find that most of the correlations between genes can be explained by only a few latent factors, which provides a new aspect for analyzing this type of data.
The rest of the paper is organized as follows. In Section 2, we derive a split Bregman method, called SBLVGG, to solve the latent variable graphical model selection problem (1). The convergence property of the algorithm is also given. SBLVGG consists of four update steps and each update has explicit formulas to calculate. In Section 3, we illustrate the utility of our algorithm and compare its performance to LogdetPPA using both simulated data and gene expression data.
2 Split Bregman method for latent variable graphical model selection
The split Bregman method was originally proposed by Osher and coauthors to solve total variation based image restoration problems [12]. It was later found to be either equivalent or closely related to a number of other existing optimization algorithms, including DouglasRachford splitting [13], the alternating direction method of multipliers (ADMM) [14, 15, 12] and the method of multipliers [16]. Because of its fast convergence and the easiness of implementation, it is increasingly becoming a method of choice for solving largescale sparsity recovery problems [17, 18]. Recently, it is also used to solve (1) and find it is very successful [8, 9].
In this section, we first reformulate the problem by introducing an auxiliary variable and then proceed to derive a split Bregman method to solve the reformulated problem. Here we would like to emphasize that, although split Bregman method has been introduced to solve graphical model problems [8, 9], we have our own contributions. Firstly, it is our first time to use split Bregman method to solve (1) and we introduce an auxiliary variable for a data fitting term instead of penalty term which has been adopted in [8, 9]. Secondly, Secondly, the update three hasn’t been appeared in [8, 9] and we provide an explicit formula for it as well. Thirdly, instead of using eig (or schur) decomposition as done in previous work [8, 9]
, we use the LAPACK routine dsyevd.f (based on a divideandconquer strategy) to compute the full eigenvalue decomposition of a symmetric matrix which is essential for updating the first and third subproblems.
2.1 Derivation of the split bregman method for latent variable graphical model selection
The loglikelihood term and the regularization terms in (1) are coupled, which makes the optimization problem difficult to solve. However, the three terms can be decoupled if we introduce an auxiliary variable to transfer the coupling from the objective function to the constraints. More specially, the problem (1) is equivalent to the following problem
(4)  
The introduction of the new variable of is a key step of our algorithm, which makes the problem amenable to a split Bregman procedure to be detailed below. Although the split Bregman method originated from Bregman iterations, it has been demonstrated to be equivalent to the alternating direction method of multipliers (ADMM) [14, 15, 19]. For simplicity of presentation, next we derive the split Bregman method using the augmented Lagrangian method [20, 16].
We first define an augmented Lagrangian function of (4)
(5)  
where is a dual variable matrix corresponding to the equality constraint , and is a parameter. Compared with the standard Lagrangian function, the augmented Lagrangian function has an extra term , which penalizes the violation of the linear constraint .
With the definition of the augmented Lagrangian function (5), the primal problem (4) is equivalent to
(6) 
Exchanging the order of and in (6) leads to the formulation of the dual problem
(7) 
Note that the gradient can be calculated by the following [21]
(8) 
where .
Applying gradient ascent on the dual problem (7) and using equation (8), we obtain the method of multipliers [16] to solve (4)
(9) 
Here we have used as the step size of the gradient ascent. It is easy to see that the efficiency of the iterative algorithm (9) largely hinges on whether the first equation of (9) can be solved efficiently. Note that the augmented Lagrangian function still contains and can not easily be solved directly. But we can solve the first equation of (9) through an iterative algorithm that alternates between the minimization of and . The method of multipliers requires that the alternative minimization of and are run multiple times until convergence to get the solution . However, because the first equation of (9) represents only one step of the overall iteration, it is actually not necessary to be solved completely. In fact, the split Bregman method (or the alternating direction method of multipliers [14]) uses only one alternative iteration to get a very rough solution of (9), which leads to the following iterative algorithm for solving (4) after some reformulations,
(10) 
2.1.1 Convergence
2.1.2 Explicit formulas to update and
We first focus on the computation of the first equation of (10). Taking the derivative of the objective function and setting it to be zero, we get
(11) 
It is a quadratic equation where the unknown is a matrix. The complexity for solving this equation is at least because of the inversion involved in (11). Note that and , if is symmetric, so is . It is easy to check that the explicit form for the solution of (11) under constraint , i.e., , is
(12) 
where and denotes the square root of a symmetric positive definite matrix . Recall that the square root of a symmetric positive definite matrix
is defined to be the matrix whose eigenvectors are the same as those of
and eigenvalues are the square root of those of . Therefore, to get the update of , can first compute the eigenvalues and eigenvectors of , and then get the eigenvalues of according to (12) by replacing the matrices by the corresponding eigenvalues. We adopt the LAPACK routine dsyevd.f (based on a divideandconquer strategy) to compute the full eigenvalue decomposition of . It is about times faster than eig (or schur) routine when is larger than .For the second equation of (10), we have made the data fitting term separable with respect to the entries of . Thus, it is very easy to get the solution and the computational complexity would be for is also separable. Let be a soft thresholding operator defined on matrix space and satisfying
where . Then the update of is
For the update of , it can use the following Theorem.
Theorem 2.
Given a symmetric matrix and . Denote
Then , where are the eigenvalues of with being the corresponding eigenvector matrix and .
Proof.
Note that where
is the identity matrix. Thus,
Compute eigenvalue decomposition on matrix and get , where and is the diagonal matrix. ThenTogether with the fact that , should satisfy for and otherwise. Therefore, .
Using the operator defined in Theorem 2, it is easy to see that
(13) 
Here we also use the LAPACK routine dsyevd.f (based on a divideandconquer strategy) to compute the full eigenvalue decomposition of . Summarizing all together, we get SBLVGG to solve the latent variable Gaussian Graphical Model (1) as shown in Algorithm 1. [htb]
3 Numerical experiments
Next we illustrate the efficiency of the split Bregman method (SBLVGG) for solving (1) using time trials on artificial data as well as gene expression data. All the algorithms were implemented in Matlab and run on a 64bit linux desktop with Intel i3  3.2GHz QuadCore CPU and 8GB memory. To evaluate the performance of SBLVGG, we compare it with logdetPPA [11] which is stateofart solver for (1) in largescale case. LogdetPPA was originally developed for logdeterminant semidefinite programs with smooth penalties. In order to solve (1) using LogdetPPA , we need to reformulate (1) as a smooth problem as done in [10], which results in the derived sparse matrix not strictly sparse with many entries close to but not exactly . We also demonstrate that latent variable Gaussian graphical selection model (1) is better than sparse Gaussian graphical model (1) in terms of generalization ability using gene expression data.
Note that the convergence of Algorithm LABEL:algorithm_SBLVGG is guaranteed no matter what values of is used as shown in Theorem 1. The speed of the algorithm can, however, be influenced by the choices of as it would affect the number of iterations involved. In our implementation, we choose in for artificial data and for gene expression data.
3.1 Artificial data
Let with being the total number of variables in the graph, the number of observed variables and the number of hidden variables. The synthetic data are generated in a similar way as the one in Section 6.1 of [11]. First, we generate an random sparse matrix
with nonzero entries drawn from normal distribution
. Then set()  Method  Obj. Value  Rank  Sparse Ratio 

(0.0025, 0.21)  SBLVGG  5642.6678  8  5.56% 
lodgetPPA  5642.6680  8  99.97%  
(0.0025, 0.22)  SBLVGG  5642.4894  3  5.58% 
lodgetPPA  5642.4895  3  99.97%  
(0.0027, 0.21)  SBLVGG  5619.2744  16  4.14% 
lodgetPPA  5619.2746  16  99.97%  
(0.0027, 0.22)  SBLVGG  5619.0194  6  4.17% 
lodgetPPA  5619.0196  6  99.97% 
Note that is marginal precision matrix of observed variables. We generate Gaussian random samples from , and calculates its sample covariance matrix . In our numerical experiments, we set sparse ratio of around , and . The stopping criteria for SBLVGG is specified as follows. Let . We stop our algorithm if and with .
Figure 1(a) shows CPU time curve of SBLVGG and LogdetPPA with respect to the number of variable for the artificial data. For each fixed , the CPU time is averaged over 4 runs with four different pairs. We can see SBLVGG consistently outperform LogdetPPA. For dimension of 2500 or less, it is 3.5 times faster on average. For dimension 3000, it is 4.5 times faster. This also shows SBLVGG scales better to problem size than LogdetPPA. In terms of accuracy, Table 1 summarize performance of two algorithms at , in three aspects: objective value, rank of , sparsity of (ratio of nonzero offdiagonal elements). We find in terms of objective value and rank, both algorithms generate almost identical results. However, SBLVGG outperform LogdetPPA due to its soft thresholding operator in Algorithm LABEL:algorithm_SBLVGG for , while LogdetPPA misses this kind of operator and result in many nonzero but close to zero entries due to numerical error. We would like to emphasize that the results in lower dimensions are very similar to , . We omit the details here due space limitation.
3.2 Gene expression data
The gene expression data [23] contains measurements of mRNA expression level of the genes of S. cerevisiae (yeast) under different experimental conditions. First we centralize the data and choose three subset of the data, , and
genes with highest variances. Figure
1(b) shows CPU time of SBLVGG and LogdetPPA with different . We can see that SBLVGG consistently perform better than LogdetPPA: in 1000()  Algorithm  Obj. Value  Rank  # Non0 Entries 

(0.01,0.05)  SBLVGG  9793.3451  88  34 
LodgetPPA  9793.3452  88  8997000  
(0.01,0.1)  SBLVGG  9607.8482  60  134 
LodgetPPA  9607.8483  60  8997000  
(0.02,0.05)  SBLVGG  8096.2115  79  0 
LodgetPPA  8096.2115  79  8996998  
(0.02,0.1)  SBLVGG  8000.9047  56  0 
LodgetPPA  8000.9045  56  8997000 
dimension case, SBLVGG is times faster, while in 2000 and 3000 dimension case, almost 3 times faster. Table 2 summarize the accuracy for dimension case in three aspects: objective value, rank of , sparsity of (Number of nonzero offdiagonal elements) for four fixed pair of . Similar to artificial data, SBLVGG and LogdetPPA generate identical results in terms of objective value and number of hidden units. However, logdetPPA suffers from the floating point problem of not being able to generate exact sparse matrix. On the other hand, SBLVGG is doing much better in this aspect.
Exp. Number  LVGG  SGG  

Rank of  Sparsity of  Sparsity of  
1  48  30  2191.3  24734  1728.8 
2  47  64  2322.7  28438  1994.1 
3  50  58  2669.9  35198  2526.3 
4  52  64  2534.6  30768  2282.5 
5  48  0  2924.0  29880  2841.4 
6  51  52  2707.1  28754  2642.6 
7  45  0  2873.3  30374  2801.4 
8  49  0  2765.5  31884  2536.7 
9  48  54  2352.0  29752  2087.2 
10  47  0  2922.9  29760  2843.5 
We also investigated generalization ability of latent variable Gaussian graphical selection model (1) versus sparse Gaussian graphical model (1) using this data set. A subset of the data, 1000 genes with highest variances, are used for this experiment. The 300 samples are randomly divided into 200 for training and 100 for testing. Denote the negative log likelihood (up to a constant difference)
where is the empirical covariance matrix using observed sample data and is the estimated covariance matrix based on model (1) or model (1). It easy to see that is equivalent to negative Loglikelihood function up to some scaling. Therefore, we use as a criteria for crossvalidation or prediction. Regularization parameters for model (1) and for model (1) are selected by 10fold cross validation on training set. Table 3 shows that latent variable Gaussian graphical selection model (1) consistently outperform sparse Gaussian graphical model (1) in terms of generalization ability using criteria . We also note that latent variable Gaussian graphical selection model (1) tend to use moderate number of hidden units, and very sparse conditional correlation to explain the data. For , it tend to predict about hidden units, and the number of direct interconnections between observed variables are tens, and sometimes even 0. This suggests that most of the correlations between genes observed in the mRNA measurement can be explained by only a small number of latent factors. Currently we only tested the generalization ability of latent variable Gaussian graphical selection model using . The initial result with gene expression data is encouraging. Further work (model selection and validation) will be done by incorporating other prior information or by comparing with some known gene interactions.
4 Discussion
Graphical model selection in highdimension arises in a wide range of applications. It is common that in many of these applications, only a subset of the variables are directly observable. Under this scenario, the marginal concentration matrix of the observed variables is generally not sparse due to the marginalization of latent variables. A computational attractive approach is to decompose the marginal concentration matrix into a sparse matrix and a lowrank matrix, which reveals the conditional graphical model structure in the observed variables as well as the number of and effect due to the hidden variables. Solving the regularized maximum likelihood problem is however nontrivial for largescale problems, because of the complexity of the loglikelihood term, the trace norm penalty and norm penalty. In this work, we propose a new approach based on the split Bregman method (SBLVGG) to solve it. We show that our algorithm is at least three times faster than the stateofart solver for largescale problems.
We applied the method to analyze the expression of genes in yeast in a dataset consisting of thousands of genes measured over 300 different experimental conditions. It is interesting to note that the model considering the latent variables consistently outperforms the one without considering latent variables in term of testing likelihood. We also note that most of the correlations observed between mRNAs can be explained by only a few dozen latent variables. The observation is consistent with the module network idea proposed in the genomics community. It also might suggest that the postranscriptional regulation might play more prominent role than previously appreciated.
References
 Dempster [1972] A.P. Dempster. Covariance selection. Biometrics, 28(1):157–175, 1972. ISSN 0006341X.
 Peng et al. [2009] J. Peng, P. Wang, N. Zhou, and J. Zhu. Partial correlation estimation by joint sparse regression models. Journal of the American Statistical Association, 104(486):735–746, 2009. ISSN 01621459.

Banerjee et al. [2008]
O. Banerjee, L. El Ghaoui, and A. d’Aspremont.
Model selection through sparse maximum likelihood estimation for
multivariate gaussian or binary data.
The Journal of Machine Learning Research
, 9:485–516, 2008. ISSN 15324435.  Friedman et al. [2008] J. Friedman, T. Hastie, and R. Tibshirani. Sparse inverse covariance estimation with the graphical lasso. Biostatistics, 9(3):432–441, 2008. ISSN 14654644.
 Meinshausen and Bühlmann [2006] N. Meinshausen and P. Bühlmann. Highdimensional graphs and variable selection with the lasso. The Annals of Statistics, 34(3):1436–1462, 2006. ISSN 00905364.
 Yuan and Lin [2007] M. Yuan and Y. Lin. Model selection and estimation in the gaussian graphical model. Biometrika, 94(1):19–35, 2007. ISSN 00063444.
 Lu [2008] Z. Lu. Smooth optimization approach for sparse covariance selection. SIAM J. Optim., 19(4):1807–1827, 2008. ISSN 10526234.
 Scheinberg et al. [2010] K. Scheinberg, S. Ma, and D. Goldfarb. Sparse inverse covariance selection via alternating linearization methods. Advances in Neural Information Processing Systems (NIPS), 2010.
 Yuan [2009] X. Yuan. Alternating direction methods for sparse covariance selection. preprint, 2009.
 Chandrasekaran et al. [2010] V. Chandrasekaran, P.A. Parrilo, and A.S. Willsky. Latent variable graphical model selection via convex optimization. In Communication, Control, and Computing (Allerton), 2010 48th Annual Allerton Conference on, pages 1610–1613. IEEE, 2010.
 Wang et al. [2010] C. Wang, D. Sun, and K.C. Toh. Solving logdeterminant optimization problems by a newtoncg primal proximal point algorithm. SIAM J. Optimization, pages 2994–3013, 2010.
 Goldstein and Osher [2009] T. Goldstein and S. Osher. The split Bregman method for regularized problems. SIAM J. Imaging Sci., 2(2):323–343, 2009. ISSN 19364954.
 Wu and Tai [2010] C. Wu and X.C. Tai. Augmented Lagrangian method, dual methods, and split Bregman iteration for ROF, vectorial TV, and high order models. SIAM Journal on Imaging Sciences, 3:300, 2010.
 Gabay and Mercier [1976] D. Gabay and B. Mercier. A dual algorithm for the solution of nonlinear variational problems via finite element approximation. Computers and Mathematics with Applications, 2(1):17–40, 1976.
 [15] R. Glowinski and A. Marroco. Sur l’approximation, par éléments finis d’ordre un, et la résolution, par penalisationdualité, d’une classe de problèmes de dirichlet non linéaires. Rev. Franc. Automat. Inform. Rech. Operat.
 Rockafellar [1973] R. T. Rockafellar. A dual approach to solving nonlinear programming problems by unconstrained optimization. Math. Programming, 5:354–373, 1973. ISSN 00255610.
 Cai et al. [2009] J. F. Cai, S. Osher, and Z. Shen. Split bregman methods and frame based image restoration. Multiscale Model. Simul., 8(2):337–369, 2009.
 Candes et al. [2009] E.J. Candes, X. Li, Y. Ma, and J. Wright. Robust principal component analysis? Arxiv preprint arXiv:0912.3599, 2009.

Setzer [2009]
S. Setzer.
Split bregman algorithm, douglasrachford splitting and frame
shrinkage.
Scale space and variational methods in computer vision
, pages 464–476, 2009.  Hestenes [1969] M. R. Hestenes. Multiplier and gradient methods. J. Optimization Theory Appl., 4:303–320, 1969. ISSN 00223239.
 Bertsekas [1982] D.P. Bertsekas. Constrained optimization and lagrange multiplier methods. 1982.
 Eckstein and Bertsekas [1992] J. Eckstein and D.P. Bertsekas. On the douglas rachford splitting method and the proximal point algorithm for maximal monotone operators. Mathematical Programming, 55(1):293–318, 1992. ISSN 00255610.
 Hughes et al. [2000] T.R. Hughes, M.J. Marton, A.R. Jones, C.J. Roberts, R. Stoughton, C.D. Armour, H.A. Bennett, E. Coffey, H. Dai, Y.D. He, et al. Functional discovery via a compendium of expression profiles. Cell, 102(1):109–126, 2000. ISSN 00928674.