Hierarchical Bayesian Quadrature

Numerical integration is essential for scientific computing and machine learning, often requiring the evaluation of complex, black-box functions. Standard Bayesian Quadrature (BQ) uses Gaussian processes (GPs) with stationary kernels to estimate integrals and quantify uncertainty. However, these models struggle with nonstationary integrands—functions where the local variability changes significantly across the domain—leading to model misspecification and inefficient use of the evaluation budget. The authors ask: can we adaptively partition the integration domain to allow local stationary models to capture nonstationary behavior while maintaining the computational benefits of BQ?

Hierarchical Bayesian Quadrature (HBQ) addresses this by recursively partitioning the integration domain into a tree of axis-aligned hyperrectangles. Each leaf node in the tree fits a local stationary GP, allowing the model to adapt its lengthscale and smoothness parameters to the local behavior of the integrand. To compute the global integral, the authors introduce a tree-based conditioning mechanism that recombines local estimates while reintroducing cross-subdomain correlations, ensuring the global uncertainty estimate remains well-calibrated. The tree structure is grown adaptively using a Bayesian Information Criterion (BIC) to decide when to split a domain and an integral variance reduction (IVR) criterion to determine where to allocate additional function evaluations.

HBQ successfully automates the detection of nonstationary behavior, concentrating function evaluations in regions of high complexity (such as sharp ridges in likelihood surfaces) while using fewer points in flatter regions. The method requires no MCMC sampling and remains computationally efficient because domain splitting reduces the size of the local Gram matrices. Experimental results on benchmark problems and an epidemiological SIR model demonstrate that HBQ significantly outperforms standard BQ on nonstationary integrands and matches its performance on stationary ones, effectively bridging the gap between flexibility and computational tractability.