Integrals of linearly constrained multivariate Gaussian densities are a
frequent problem in machine learning and statistics, arising in tasks like
generalized linear models and Bayesian optimization. Yet they are notoriously
hard to compute, and to further complicate matters, the numerical values of
such integrals may be very small. We present an efficient black-box algorithm
that exploits geometry for the estimation of integrals over a small, truncated
Gaussian volume, and to simulate therefrom. Our algorithm uses the
Holmes-Diaconis-Ross (HDR) method combined with an analytic version of
elliptical slice sampling (ESS). Adapted to the linear setting, ESS allows for
rejection-free sampling, because intersections of ellipses and domain
boundaries have closed-form solutions. The key idea of HDR is to decompose the
integral into easier-to-compute conditional probabilities by using a sequence
of nested domains. Remarkably, it allows for direct computation of the
logarithm of the integral value and thus enables the computation of extremely
small probability masses. We demonstrate the effectiveness of our tailored
combination of HDR and ESS on high-dimensional integrals and on entropy search
for Bayesian optimization