Efficient Latent Variable Graphical Model Selection via Split Bregman Method

We consider the problem of covariance matrix estimation in the presence of latent variables. Under suitable conditions, it is possible to learn the marginal covariance matrix of the observed variables via a tractable convex program, where the concentration matrix of the observed variables is decomposed into a sparse matrix (representing the graphical structure of the observed variables) and a low rank matrix (representing the marginalization effect of latent variables). We present an efficient first-order method based on split Bregman to solve the convex problem. The algorithm is guaranteed to converge under mild conditions. We show that our algorithm is significantly faster than the state-of-the-art algorithm on both artificial and real-world data. Applying the algorithm to a gene expression data involving thousands of genes, we show that most of the correlation between observed variables can be explained by only a few dozen latent factors.

💡 Research Summary

The paper addresses the challenging problem of estimating the covariance matrix of observed variables when latent (unobserved) variables influence the data. In many high‑dimensional applications—such as genomics, recommender systems, or any setting where only a subset of variables is measured—the marginal covariance of the observed variables can be dense because the effect of the hidden variables is “integrated out.” Classical Gaussian graphical model estimation assumes the precision (inverse covariance) matrix is sparse and solves a convex ℓ₁‑regularized log‑determinant program. However, this sparsity‑only assumption fails when latent factors induce strong correlations.

The authors adopt the latent‑variable graphical model introduced by Chandrasekaran et al., which decomposes the marginal precision matrix K̂_O of the observed variables into a sparse component S (capturing conditional independencies among observed variables) and a low‑rank component L (capturing the marginalization effect of a few latent variables). The resulting convex optimization problem is

  min_{S,L} –log det(S–L) + tr(Σ̂_O (S–L)) + λ₁‖S‖₁ + λ₂ tr(L)
  subject to S–L ≽ 0, L ≽ 0.

Here Σ̂_O is the empirical covariance of the observed data, λ₁ controls sparsity, and λ₂ controls the rank via the trace norm. This formulation is strictly convex and has a unique solution, but solving it efficiently is non‑trivial because the log‑determinant term, the ℓ₁ penalty, and the semidefinite constraints are all coupled.

State‑of‑the‑art solvers such as LogdetPPA handle smooth penalties and require reformulating the problem, which leads to loss of exact sparsity and additional heuristic thresholding. Moreover, LogdetPPA is not tailored to the specific structure of the latent‑variable model.

The main contribution of the paper is a new first‑order algorithm, called SBL‑VGG (Split Bregman for Latent Variable Gaussian Graphical model), derived from the split Bregman (or ADMM) framework. The authors introduce an auxiliary variable A = S – L to decouple the log‑determinant term from the regularizers. The augmented Lagrangian is

 L(A,S,L,U) = –log det A + tr(Σ̂_O A) + λ₁‖S‖₁ + λ₂ tr(L)
      + ⟨U, A – S + L⟩ + (μ/2)‖A – S + L‖_F²,

with dual matrix U and penalty parameter μ > 0. Alternating minimization over (A,S,L) followed by a dual ascent yields four simple update steps:

1. A‑update: solve –A⁻¹ + Σ̂_O + U_k + μ(A – S_k + L_k) = 0.
   This matrix equation admits a closed‑form solution using eigenvalue decomposition:
     A_{k+1} = K_k + √(K_k² + 4μI) / (2μ), where K_k = μ(S_k – L_k) – Σ̂_O – U_k.
   The square‑root of a symmetric positive‑definite matrix is computed by diagonalizing K_k² + 4μI. The authors use LAPACK’s dsyevd.f routine (divide‑and‑conquer) which is about ten times faster than MATLAB’s eig for dimensions > 500.

2. S‑update: minimize λ₁‖S‖₁ + (μ/2)‖A_{k+1} – S + L_k + U_k‖F².
   This is separable element‑wise, yielding a soft‑thresholding operation:
     S{k+1} = T_{λ₁/μ}(A_{k+1} + L_k + μ⁻¹U_k).

3. L‑update: minimize λ₂ tr(L) + (μ/2)‖A_{k+1} – S_{k+1} + L + U_k‖F² subject to L ≽ 0.
   The solution is the proximal operator of the trace norm with a PSD constraint: compute the eigen‑decomposition of X = S{k+1} – A_{k+1} – μ⁻¹U_k, then set
     L_{k+1} = V diag((λ_i – λ₂/μ)+) Vᵀ,
   where λ_i are eigenvalues of X and (·)+ denotes the positive part.

4. Dual update: U_{k+1} = U_k + μ(A_{k+1} – S_{k+1} + L_{k+1}).

The authors prove convergence of the iteration by invoking standard ADMM theory; the algorithm converges to the unique minimizer for any μ > 0.

Computational aspects

• The dominant cost is the eigen‑decomposition required in the A‑ and L‑updates, each O(p³) in the worst case, but highly optimized LAPACK routines and the possibility of parallelization make the method practical up to several thousand variables.
• The S‑update is O(p²) and trivially parallelizable.
• Memory usage is modest: only a few p×p matrices need to be stored.

Experimental evaluation
The authors implement SBL‑VGG in MATLAB and compare it against LogdetPPA on synthetic data (p = 200, 500, 1000) with varying sample sizes and numbers of latent factors. Across all settings, SBL‑VGG achieves the same or better estimation accuracy (measured by support recovery of S and rank recovery of L) while being 15–30 times faster. Convergence is typically reached within a few dozen iterations.

A real‑world test uses a gene expression dataset containing roughly 2,000 genes. After cross‑validating λ₁ and λ₂, the algorithm recovers a low‑rank component of rank ≈ 35, explaining the majority of the observed correlations. The sparse component yields a biologically plausible network of gene–gene interactions, while the low‑rank part suggests that a handful of latent biological processes (e.g., transcription factors, signaling pathways) drive most of the co‑expression structure.

Significance and limitations
The paper demonstrates that a split‑Bregman/ADMM approach can exploit the specific structure of latent‑variable graphical models to obtain a fast, scalable algorithm that preserves exact sparsity and low‑rankness. The closed‑form updates avoid inner iterative solvers, leading to substantial speedups over generic log‑determinant SDP solvers. However, the eigen‑decomposition step still scales cubically, which may become a bottleneck for problems with tens of thousands of variables. Future work could explore randomized low‑rank approximations or block‑coordinate schemes to further reduce computational cost, as well as extensions to non‑Gaussian data, dynamic networks, or online settings.

In summary, the paper provides a solid methodological contribution—an efficient split‑Bregman algorithm for latent‑variable Gaussian graphical model selection—backed by rigorous convergence analysis and compelling empirical results on both synthetic and large‑scale biological data. This work is likely to be of interest to researchers in high‑dimensional statistics, machine learning, and computational biology who need to disentangle sparse conditional dependencies from latent factor effects.