Asymptotically Exact 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. 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 we show with both theoretical results and an application to Alzheimer's data, such a routine requires a number of iterations converging to one as p goes 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.