An Order of Magnitude Time Complexity Reduction for Gaussian Graphical Model Posterior Sampling Using a Reverse Telescoping Block Decomposition

1. Problem Statement

The paper addresses the computational bottleneck in fully Bayesian posterior estimation for undirected Gaussian Graphical Models (GGMs) in high-dimensional regimes where the number of variables $p$ far exceeds the sample size $n$ ($p \gg n$).

• Context: GGMs model conditional independence relationships via the sparsity of the precision matrix $\Theta = \Sigma^{-1}$.
• The Challenge: While conjugate priors (like the Wishart distribution) allow for efficient posterior sampling with $O(p^3)$ complexity using Bartlett's decomposition, they do not naturally induce sparsity in the off-diagonal elements of $\Theta$.
• The Bottleneck: To enforce sparsity, researchers use element-wise non-conjugate priors (e.g., Graphical Horseshoe, Graphical Lasso). The state-of-the-art method for sampling under these priors, developed by Wang (2012) and extended by others, relies on a specific reparameterization to maintain the positive definiteness constraint ($\Theta \in \mathcal{M}^+_p$). However, this approach requires inverting a $(p-1) \times (p-1)$ matrix at every step of a cyclical Gibbs sampler, resulting in a per-iteration complexity of $O(p^4)$. This becomes computationally prohibitive when $p$ is large, even if $n$ is small.

2. Methodology

The authors propose a novel Reverse Telescoping Block Decomposition sampler that targets the exact posterior distribution without using variational approximations or pseudo-likelihoods.

Core Innovation: Reversing the Telescoping Decomposition

The method builds upon a decomposition recently introduced by Bhadra et al. (2024) for marginal likelihood calculation but repurposes it for sampling.

1. Forward Decomposition (Bhadra et al., 2024): Transforms the precision matrix $\Theta$ into a new parameterization $\tilde{\Theta}$ using a sequence of Schur complements. This allows for tractable marginal likelihood calculations but was previously too slow ($O(p^5)$) for sampling.
2. Reverse Decomposition (The Authors' Contribution): The authors invert this mapping. Instead of sampling $\Theta$ directly, they sample the transformed parameters $\tilde{\Theta} = \{\tilde{\theta}_j\}$, where each $\tilde{\theta}_j$ corresponds to a column of the precision matrix adjusted by the Schur complements of subsequent columns.
   • Likelihood Factorization: The joint likelihood $f(y|\Theta)$ is factorized into a product of conditional densities: $f(y|\Theta) = \prod_{j=1}^p f(y_j | y_{1:j-1}, \tilde{\theta}_j)$.
   • Conditional Structure: Under this parameterization, the conditional distribution of $y_j$ given previous variables $y_{1:j-1}$ depends only on the specific block $\tilde{\theta}_j$. This structure resembles a sequence of partial regressions.

The Sampling Algorithm (Block Gibbs Sampler)

The proposed sampler iterates through $j = p, \dots, 1$:

1. Induced Priors: The element-wise priors on $\Theta$ (e.g., Graphical Horseshoe) induce specific, tractable priors on the transformed blocks $\tilde{\theta}_j$.
2. Block Sampling: For each $j$, the sampler draws $\tilde{\theta}_j$ conditional on the data $y$, the current values of $\tilde{\theta}_{j+1}, \dots, \tilde{\theta}_p$, and the hyperparameters.
   • The diagonal term $\tilde{\theta}_{jj}$ is sampled from a Gamma or Generalized Inverse Gaussian (GIG) distribution.
   • The off-diagonal terms $\tilde{\theta}_{\cdot j}$ are sampled from a Multivariate Normal distribution. Crucially, the covariance matrix of this normal distribution can be updated efficiently using the method of Bhattacharya et al. (2016), which avoids direct matrix inversion by exploiting the structure of the data matrix $y$.
3. Reconstruction: After sampling all $\tilde{\theta}_j$, the algorithm applies the reverse telescoping update (Algorithm 1b) to reconstruct the original precision matrix $\Theta$.

3. Key Contributions

• Complexity Reduction: The primary contribution is reducing the per-iteration computational complexity from $O(p^4)$ to $O(p^3)$ in the regime where $n = o(p^{1/2})$. This brings the cost of non-conjugate element-wise priors down to the same level as the conjugate Wishart family.
• Exact Posterior Targeting: Unlike many fast methods that rely on pseudo-likelihoods (e.g., Meinshausen-Bühlmann) or variational approximations, this method targets the true posterior distribution exactly.
• General Applicability: The framework is generic and applies to any element-wise prior of the form of a normal scale mixture (including Graphical Horseshoe, Graphical Horseshoe-like, Graphical R2-D2, etc.).
• Preservation of Constraints: The reparameterization naturally enforces the positive definiteness of $\Theta$ without requiring rejection sampling or complex constraints, as the mapping between $\Theta$ and $\tilde{\Theta}$ is one-to-one and preserves the positive definite property.

4. Results

The authors validated the method through extensive simulations and a real-world application.

• Simulation Setup: They compared the proposed Reverse Telescoping (RT) sampler against the standard Cyclical sampler (Wang, 2012) under Graphical Horseshoe (GHS) and Graphical Horseshoe-like (GHSL) priors. Dimensions ranged from $p=100$ to $p=800$ with $n \approx 24 \log p$.
• Estimation Accuracy: The posterior estimates (Frobenius norm of the precision matrix and deviation from the truth) produced by the RT sampler were numerically indistinguishable from the Cyclical sampler, confirming the correctness of the algorithm.
• Computational Performance:
  • Runtime: The RT sampler demonstrated a massive speedup. At $p=400$, the RT sampler was 8.2 times faster in raw runtime than the Cyclical sampler.
  • Scalability: As $p$ increased, the Cyclical sampler failed to complete runs within 12 hours for $p > 400$, whereas the RT sampler scaled efficiently.
  • Mixing: The RT sampler achieved 2–5 times higher Effective Sample Size per minute (ESS/min) compared to the Cyclical sampler.
• Real Data Application: Applied to TCGA breast cancer RNA-seq data ($p=139, n=90$). The RT sampler was significantly faster (144s vs 700s for GHS) while inferring nearly identical gene interaction networks.

5. Significance

This work resolves a critical barrier in high-dimensional Bayesian statistics. By achieving $O(p^3)$ complexity for non-conjugate sparse priors, it makes fully Bayesian inference feasible for large-scale problems (e.g., genomics, neuroimaging) where $p$ is in the hundreds or thousands.

• Parity with Conjugate Models: It allows researchers to use sophisticated, sparsity-inducing priors (which are statistically superior for structure learning) without paying the prohibitive computational penalty previously associated with them.
• Algorithmic Advancement: The "Reverse Telescoping" concept provides a new tool for manipulating precision matrices in GGMs, potentially applicable to other problems like marginal likelihood estimation or multi-graph inference.
• Practical Impact: The availability of the RTsampler R package allows immediate adoption of these faster methods by the statistical community.