# Prior-preconditioned conjugate gradient method for accelerated Gibbs sampling in "large n & large p" Bayesian sparse regression

arXiv: https://arxiv.org/abs/1810.12437 · HTML: https://mlsanity.com/paper/1810.12437/

## Key facts

- Cited by 6 papers in the map (data through 14 Oct 2026).
- First posted on arXiv 29 Oct 2018.

## Built on by (most cited first, 6 of 6)

- [Query lower bounds for log-concave sampling](https://mlsanity.com/paper/2304.02599/) (2023)
- [Iterative Methods for Vecchia-Laplace Approximations for Latent Gaussian Process Models](https://mlsanity.com/paper/2310.12000/) (2023)
- [Gradient GA: Gradient Genetic Algorithm for Drug Molecular Design](https://mlsanity.com/paper/2502.09860/) (2025)
- [Subset selection for linear mixed models](https://mlsanity.com/paper/2107.12890/) (2021)
- [Conjugate gradient methods for high-dimensional GLMMs](https://mlsanity.com/paper/2411.04729/) (2024)
- [Scalable Krylov Subspace Methods for Generalized Mixed-Effects Models with Crossed Random Effects](https://mlsanity.com/paper/2505.09552/) (2025)

## Abstract (arXiv)

In a modern observational study based on healthcare databases, the number of observations and of predictors typically range in the order of 105 ~ 106 and of 104 ~ 105. Despite the large sample size, data rarely provide sufficient information to reliably estimate such a large number of parameters. Sparse regression techniques provide potential solutions, one notable approach being the Bayesian methods based on shrinkage priors. In the "large n & large p" setting, however, posterior computation encounters a major bottleneck at repeated sampling from a high-dimensional Gaussian distribution, whose precision matrix Phi is expensive to compute and factorize. In this article, we present a novel algorithm to speed up this bottleneck based on the following observation: we can cheaply generate a random vector b such that the solution to the linear system Phi β = b has the desired Gaussian distribution. We can then solve the linear system by the conjugate gradient (CG) algorithm through matrix-vector multiplications by Phi; this involves no explicit factorization or calculation of Phi itself. Rapid convergence of CG in this context is guaranteed by the theory of prior-preconditioning we develop. We apply our algorithm to a clinically relevant large-scale observational study with n = 72,489 patients and p = 22,175 clinical covariates, designed to assess the relative risk of adverse events from two alternative blood anti-coagulants. Our algorithm demonstrates an order of magnitude speed-up in posterior inference, in our case cutting the computation time from two weeks to less than a day.