Conjugate gradient-accelerated Gibbs sampler for "large n and large p" sparse Bayesian logistic regression

In a modern observational study based on healthcare databases, the number of observations is typically in the order of 10^5 ~ 10^6 and that of the predictors in the order of 10^4 ~ 10^5. Despite the large sample size, the data rarely provide enough information to reliably estimate such a large number of parameters. Sparse regression provides a potential solution to this problem. There is a rich literature on desirable theoretical properties of the Bayesian approach based on shrinkage prior. On the other hand, the development of scalable methods for the required posterior computation has largely been limited to the p >> n case. While shrinkage priors are designed to make the posterior amenable to Gibbs sampling, a major computational bottleneck still arises from the need to sample from a high-dimensional Gaussian distribution at each iteration. The closed form expression for the precision matrix ¿ is available, but computing and factorizing such a large matrix is computationally expensive nonetheless. 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 of a linear system ¿¿ = b has the desired Gaussian distribution. An accurate solution of the linear system can then be found by the conjugate gradient algorithm with only a small number of the matrix-vector multiplications by ¿, without ever explicitly inverting ¿. [To read more click on the event link]