  Modern methods for Bayesian regression beyond the Gaussian response setting
are often computationally impractical or inaccurate in high dimensions. In
fact, as discussed in recent literature, bypassing such a trade-off is still an
open problem even in routine binary regression models, and there is limited
theory on the quality of variational approximations in high-dimensional
settings. To address this gap, we study the approximation accuracy of
routinely-used mean-field variational Bayes solutions in high-dimensional
probit regression with Gaussian priors, obtaining novel and practically
relevant results on the pathological behavior of such strategies in uncertainty
quantification, point estimation and prediction. Motivated by these results, we
further develop a new partially-factorized variational approximation for the
posterior of the probit coefficients which 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 approximation belongs
to a tractable class of unified skew-normal densities that crucially
incorporates skewness and, unlike for state-of-the-art mean-field solutions,
converges to the exact posterior density as p goes to infinity. To solve the
variational optimization problem, we derive a tractable CAVI algorithm that
easily scales to p in the tens of thousands, and provably requires a number of
iterations converging to 1 as p goes to infinity. Such findings are also
illustrated in extensive empirical studies where our novel solution is shown to
improve the approximation accuracy of mean-field variational Bayes for any n
and p, with the magnitude of these gains being remarkable in those
high-dimensional p>n settings where state-of-the-art methods are
computationally impractical.
