Scalable and Accurate Variational Bayes for High-Dimensional Binary Regression Models

State-of-the-art methods for Bayesian inference on regression models with binary responses are either computationally impractical or inaccurate in high dimensions. To cover this gap we propose a novel variational approximation for the posterior distribution of the coefficients in high-dimensional probit regression with Gaussian priors. Our method leverages a representation with global and local variables but, unlike for classical mean-field assumptions, it avoids a fully factorized approximation, and instead assumes a factorization only for the local variables. We prove that the resulting variational approximation belongs to a tractable class of unified skew-normal distributions that preserves the skewness of the actual posterior and, unlike for state-of-the-art variational Bayes solutions, converges to the exact posterior as the number of predictors p increases. A scalable coordinate ascent variational algorithm is proposed to obtain the optimal parameters of the approximating densities. As shown in theoretical studies and with an application to Alzheimer's data, this routine requires a number of iterations converging to one as p diverges to infinity, and can easily scale to large p settings where expectation-propagation and state-of-the-art Markov chain Monte Carlo algorithms are computationally impractical.