Source-linked AI summary
Fast sampling with Gaussian scale-mixture priors in high-dimensional regression
Anirban Bhattacharya, Antik Chakraborty, Bani K. Mallick
TL;DR
High-dimensional regression with global-local priors requires repeated sampling from structured Gaussian conditionals, while Cholesky-based methods become costly as dimension grows. The paper presents an exact data-augmentation sampler using matrix multiplications and linear system solutions, and shows computational complexity linear in p. The method is intended for Gaussian scale-mixture priors and is illustrated with a horseshoe regression model.
Problem
Repeated Cholesky factorization makes sampling Gaussian conditional posteriors a bottleneck for large p, especially when the factor changes at every iteration.
Method
The paper introduces an exact data-augmentation algorithm for structured Gaussian distributions arising from conditionally Gaussian priors.
Results
The proposed algorithm has computational complexity that scales linearly in p, unlike the O(p3) complexity of Cholesky factorization.
Takeaways & Limitations
The algorithm is broadly applicable to high-dimensional models using Gaussian scale-mixture priors and is illustrated with horseshoe regression.
Takeaways & Limitations
The authors state that their algorithm does not extend to approximate factor models with idiosyncratic components.
Abstract
from arXiv · showhide
We propose an efficient way to sample from a class of structured multivariate Gaussian distributions which routinely arise as conditional posteriors of model parameters that are assigned a conditionally Gaussian prior. The proposed algorithm only requires matrix operations in the form of matrix multiplications and linear system solutions. We exhibit that the computational complexity of the proposed algorithm grows linearly with the dimension unlike existing algorithms relying on Cholesky factorizations with cubic orders of complexity. The algorithm should be broadly applicable in settings where Gaussian scale mixture priors are used on high dimensional model parameters. We provide an illustration through posterior sampling in a high dimensional regression setting with a horseshoe prior on the vector of regression coefficients.
1 Introduction
Global-local priors support blockwise Gibbs updates in high-dimensional regression, but sampling their Gaussian conditional posterior becomes costly as p grows. The paper introduces an exact data-augmentation algorithm whose computational complexity scales linearly in p.
- Motivation: Global-local priors use scale-mixture representations that enable block updates through Gibbs sampling across regression and other models.Applications include variable selection, wavelet denoising, factor models, covariance estimation, and time series.
- Motivation: In high-dimensional regression, p can greatly exceed n, with local and global scales controlling shrinkage and effective sparsity.The prior choices influence concentration when most regression coefficients are zero or near zero.
- Existing sampler: The conditional posterior of β is Gaussian, while local scale parameters have conditionally independent posteriors that can be updated in a block.The β update involves A = X^T X + Λ^-1 and the covariance σ^2A^-1.
- Computational challenge: Cholesky-based sampling requires factoring a matrix that changes at every iteration, creating a bottleneck for large p.The cited factorization has complexity O(p3), and precomputation is unavailable because Λ changes between iterations.
- Contribution: The proposed exact sampler uses data augmentation for the Gaussian conditional distribution and has computational complexity that scales linearly in p.The introduction presents the algorithm as an alternative to repeated Cholesky factorization.
2 The algorithm
The paper introduces an exact algorithm for sampling structured Gaussian distributions using matrix multiplications and linear-system solutions instead of Cholesky factorization. It is especially advantageous when p ≫ n, including Gaussian scale-mixture prior settings.
- Computational motivation: O(p3) Cholesky factorization becomes prohibitive because Q = ΦTΦ + D−1 lacks special structure in the present setting.The comparison concerns the competing approach based on Cholesky decomposition and subsequent linear-system solves.
- Algorithm 1 Proposed algorithm: Algorithm 1 samples u and δ independently, forms v = Φu + δ, solves (ΦDΦT + In)w = α − v, and returns θ = u + DΦTw.The procedure uses only Gaussian sampling, matrix operations, and a linear-system solution.
- Correctness: Proposition 2.1 establishes that Algorithm 1 produces θ ∼ N(µ, Σ), with µ and Σ matching the target Gaussian parameters.The result follows using the Sherman–Morrison–Woodbury identity and the covariance relationship cov(u, v) = DΦT.
- Computational complexity: O(n2p) is the proposed complexity when D is diagonal, whereas non-sparse D yields a dominating O(np2) term for calculating DΦT.The diagonal case covers global-local priors; the non-sparse bound reflects the worst-case matrix-multiplication cost.
- Empirical comparison: When p ≫ n, the algorithm offers its largest computational gains, including a speed-up factor over 250 at p = 5000 in the horseshoe-prior comparison.For 6000 Gibbs-sampler iterations with n = 100, the method takes roughly the same time as Rue (2001) at p = 500 and is over 250 times faster at p = 5000.
3 Frequentist operating characteristics in high dimensions
The simulations compare horseshoe-based and penalized estimators across covariance structures and signal strengths, while the uncertainty analysis evaluates frequentist interval coverage in high dimensions. The horseshoe performs strongly for estimation and provides adaptive intervals with nominal coverage across signal and noise variables.
- Estimation performance: The horseshoe prior is highly competitive across all simulation settings, particularly when the signal strength is weaker.Figures report ℓ1, ℓ2, and prediction errors across the two signal-strength settings.
- Estimation performance: The point-wise posterior median showed superior performance even under an ℓ2 loss.The paper relates this observation to similar findings for point-mass mixture priors.
- Estimation performance: The horseshoe prior had the best average performance in 22 of 24 simulation settings.The simulations were also repeated with p = 2500 and yielded similar conclusions.
- Uncertainty quantification: Shrinkage-prior interval lengths adapt between signal and noise variables while maintaining nominal coverage.The compared frequentist procedures produced approximately equal-sized intervals for signals and noise variables.
- Uncertainty quantification: The horseshoe and other shrinkage priors are free of tuning parameters, whereas comparator coverage depended on tuning choices.The default van de Geer et al. tuning produced poorer signal coverage, and substantial tuning was required for the reported coverage probabilities.
4 Discussion
The discussion presents the proposed algorithm as broadly applicable to high-dimensional shrinkage-prior models, while identifying extensions beyond linear regression and the need for further investigation.
- The numerical results warrant additional numerical and theoretical investigations.
- The proposed algorithm can be used for essentially all shrinkage priors in the literature and supports exhaustive prior comparisons.
- Its scope extends beyond linear regression to logistic, probit, multivariate, and high-dimensional additive models.Logistic and probit extensions use standard data augmentation, while multivariate regression uses block updates of vectorized coefficients; additive models update basis coefficients in blocks.
- The algorithm can be adapted to approximate factor models with dependent idiosyncratic errors by block updating vectorized loadings.
Appendix
The appendix constructs the sampling method through data augmentation and matrix identities, transforming an augmented Gaussian draw into a sample from the target covariance distribution.
- The Sherman–Morrison–Woodbury formula provides the algebraic identity underlying the proposed sampling construction.
- Data augmentation defines an augmented Gaussian vector whose covariance has a block structure linked to the target covariance Σ.The augmented vector ζ consists of v and u, and Γ is block diagonal with a lower p × p block related to Σ.
- The last p entries of ζ* provide a sample from N(0, Σ), after which further algebra produces the algorithm.