🤖 machine learning

Hierarchical Bayesian Quadrature

This paper introduces Hierarchical Bayesian Quadrature, an adaptive algorithm that improves numerical integration for nonstationary integrands by partitioning the domain into local stationary models and recombining their estimates through a hierarchical Gaussian process structure, achieving superior performance over standard methods without requiring MCMC.

Original authors: Tim Weiland, Toni Karvonen, Philipp Hennig

Published 2026-07-14

📖 6 min read🧠 Deep dive

CC BY 4.0

Original authors: Tim Weiland, Toni Karvonen, Philipp Hennig

Original paper licensed under CC BY 4.0 (http://creativecommons.org/licenses/by/4.0/). ✨ This is an AI-generated explanation of the paper below. It is not written or endorsed by the authors. For technical accuracy, refer to the original paper. Read full disclaimer

Imagine you are trying to guess the total amount of water in a giant, irregularly shaped swimming pool. You can't measure the whole thing at once, so you have to take samples. In the world of math and computer science, this is called numerical integration. You want to know the total "volume" of a function, but the function is a black box that only tells you its height at specific points you choose.

For a long time, the smartest way to do this was Bayesian Quadrature (BQ). Think of BQ as a super-smart artist who draws a single, smooth, continuous painting of the pool's surface based on a few samples. This artist assumes the water behaves the same way everywhere: if it's wavy in one corner, it's probably wavy everywhere else. They use a "stationary" rule, meaning their brushstrokes (mathematical assumptions) are uniform across the whole canvas.

The Problem: The "One-Size-Fits-All" Trap
The paper points out a major flaw in this approach: what if the pool isn't uniform? What if one corner has calm, flat water, while another corner has a chaotic, churning whirlpool? If your artist tries to paint the whole pool with one style, they have to compromise. They might make the flat water look too bumpy or the whirlpool look too smooth. They end up wasting their time (and your computer's money) trying to get the easy parts perfect while missing the messy parts entirely.

The authors argue that standard BQ methods, which rely on these single, uniform rules, often fail when the function they are trying to measure changes its behavior from place to place (a behavior called nonstationarity). They explicitly rule out the idea that a single, global model can efficiently handle these complex, changing landscapes without making mistakes.

The Solution: The "Tree-Pruning" Team
Enter Hierarchical Bayesian Quadrature (HBQ), the new method proposed by Tim Weiland, Toni Karvonen, and Philipp Hennig.

Imagine instead of one artist painting the whole pool, you have a team of specialists. HBQ works like a smart tree that grows branches to divide the pool into smaller, manageable rooms.

The Split: The algorithm looks at the pool. If it sees a calm area and a stormy area, it says, "Okay, let's split this room in half." It draws a line down the middle.
The Specialists: Now, instead of one artist, you have two. One specializes in the calm, flat water (using a smooth, gentle brush). The other specializes in the churning whirlpool (using a rough, detailed brush). Each specialist fits a local model just for their own room.
The Reunion: Here is the clever part. Usually, if you just add up the guesses of two separate rooms, you lose the connection between them. But HBQ uses a special "hierarchical" trick. It remembers that the two rooms were once one big room. It uses a mathematical "glue" (called tree conditioning) to recombine the local guesses into a single, global answer. This ensures that the final result still knows how the calm water and the stormy water relate to each other, keeping the uncertainty estimates accurate.

How It Decides Where to Split
The team doesn't just split randomly. They use a strict rulebook (based on something called the Bayesian Information Criterion, or BIC). Before they cut a room in half, they ask: "Is the current artist struggling so much that splitting the room will actually make the answer better?" If the answer is no, they don't split. They just add more samples to the current room. This prevents the tree from growing too big and wasting resources on simple areas.

What the Experiments Showed
The authors tested this idea in several ways, and the results were promising but specific to certain conditions:

The "Ridge" Test: They tried to calculate the "model evidence" for a disease model (an SIR model). The math here had a strange shape: a long, curved ridge where the answer mattered, surrounded by flat, boring ground. Standard BQ spread its samples evenly, like sprinkling salt on a whole pizza. HBQ, however, noticed the ridge and concentrated almost all its samples right on that curve. In these simulations, HBQ achieved a relative error of less than 1% with just 100 samples, while the other methods were still stuck above 1%.
The "Genz" Challenge: They used a standard set of test functions. On four out of five types of functions (which were mostly smooth and uniform), HBQ performed just as well as the standard method. It didn't slow things down or make things worse. But on the "Corner Peak" function—where the action was concentrated in one tiny, sharp corner—HBQ was up to 6 times more accurate than the standard method.
The "Reaction-Diffusion" Test: They looked at a physics problem involving heat and chemicals spreading through a material that had both smooth and spiky regions. Here, HBQ pulled away dramatically. At 512 samples, HBQ was roughly 35 times more accurate than the standard method.

What They Are Sure Of (and What They Aren't)
The authors are very confident that this method works for low-to-moderate dimensional problems (like 2D or 3D spaces) where the function changes its behavior in different spots. They have proved mathematically that the error of their method is controlled by the "worst" local part of the function, meaning if the local models are good, the whole answer is good.

However, they are careful to note a few limits:

They haven't solved the problem for very high dimensions (like 100 dimensions) yet. The method relies on splitting the space into rectangular boxes, which gets very hard to manage in huge dimensions.
They admit their math proof covers the final result if the tree is already built, but they haven't fully proved the math for how the tree builds itself in every single scenario yet.
They don't claim this is a magic bullet for every integration problem. If the function is already smooth and simple everywhere, HBQ just acts like the old method, doing the same job without extra magic.

In short, the paper suggests that by letting the computer "prune" the problem into smaller, specialized pieces and then carefully stitching the answers back together, we can solve complex math puzzles much faster and more accurately than before—especially when the puzzle has both calm and chaotic parts. It's a win for adaptability, not a replacement for the basics.

Technical Summary: Hierarchical Bayesian Quadrature

Problem Statement
Numerical integration is fundamental to scientific computing and probabilistic machine learning, particularly for computing model evidence and normalizing constants. Bayesian Quadrature (BQ) addresses this by modeling the integrand with a Gaussian Process (GP) surrogate, yielding a posterior distribution over the integral that includes a point estimate and quantified uncertainty. However, standard BQ relies on stationary covariance functions (kernels). This assumption leads to model misspecification when the integrand exhibits nonstationary behavior—where variability depends on location, such as sharp ridges in likelihood surfaces or multi-scale features in engineering simulations. A single stationary kernel must compromise between resolving complex regions and efficiently covering flat regions, often resulting in inefficient allocation of the evaluation budget and poor convergence.

Methodology
The authors propose Hierarchical Bayesian Quadrature (HBQ), a method that adaptively partitions the integration domain to handle nonstationarity while retaining the computational benefits of stationary models.

Domain Decomposition: The integration domain $\Omega$ is recursively partitioned into axis-aligned hyperrectangles (subdomains) via a binary tree structure. Each leaf node of the tree is assigned a local stationary GP surrogate with its own lengthscale and smoothness parameters, allowing the model to adapt to local integrand complexity.
Integral Recombination via Tree Conditioning: To obtain a global integral estimate, the method recombines local estimates from leaf nodes. Unlike simple independent summation, HBQ utilizes a hierarchical conditioning mechanism. When a node is split, the GP fitted to the parent domain is retained as a prior. The integral estimates from the children are treated as noisy observations of the parent's integral functional. By conditioning the parent GP on these child estimates, the method reintroduces cross-subdomain correlations, ensuring the global posterior reflects the structure of the entire domain.
Adaptive Tree Construction: The algorithm grows the tree on-the-fly based on two criteria:
- Splitting Criterion: A split is considered if a single stationary GP is misspecified. This is determined using the Bayesian Information Criterion (BIC) to compare the marginal likelihood of a single model against a split model. To prevent pathological splits, a minimum number of evaluation points per child is enforced.
- Refinement Strategy: The algorithm selects the next action (splitting a leaf or adding evaluation points to a leaf) by maximizing the expected reduction in the global integral variance per evaluation. This extends the Integral Variance Reduction (IVR) criterion to the hierarchical setting.
Hyperparameter Optimization: Local hyperparameters (lengthscale and smoothness) are optimized via profile marginal likelihood. A log-normal prior on the child lengthscale, anchored to the parent's fitted lengthscale, regularizes estimation in regions with few data points.

Key Contributions

Adaptive Nonstationarity Handling: HBQ automatically detects nonstationary behavior and allocates computational resources accordingly, avoiding the "one-size-fits-all" limitation of standard BQ.
Closed-Form Recombination: The method introduces a tree-conditioning framework that allows for the recombination of local integral estimates while preserving cross-domain correlations, all within a closed-form framework that avoids Markov Chain Monte Carlo (MCMC).
Computational Efficiency: By reducing the size of local Gram matrices through domain splitting, the method can accelerate GP inference. It requires no MCMC and introduces minimal computational overhead.
Theoretical Guarantees: The paper provides a convergence analysis showing that the global error is controlled by the worst-case local smoothness of the integrand (Theorem 1). It establishes that HBQ matches the convergence rates of standard BQ on stationary integrands while achieving superior rates on nonstationary ones.

Experimental Results
The authors evaluated HBQ on several benchmarks:

Genz Test Functions: On the "Corner Peak" function (characterized by sharp concentration near the origin), HBQ achieved up to $6\times$ lower error than standard BQ by concentrating evaluations in the complex region. On four other families (Oscillatory, Product Peak, Gaussian Peak, Continuous), HBQ performed comparably to standard BQ, demonstrating that the adaptive overhead does not degrade performance on stationary or well-behaved integrands.
Epidemiological Model Evidence: In computing the model evidence for an SIR model with a curved likelihood ridge, HBQ achieved sub-1% relative error with 100 evaluations, significantly outperforming standard BQ, Monte Carlo (MC), and Quasi-Monte Carlo (QMC). HBQ successfully concentrated evaluations along the ridge while providing coarse coverage to flat regions.
Structure Recovery: On a piecewise GP integrand with known nonstationary segments, HBQ correctly identified the number of segments and kernel orders in nearly all trials, yielding $3.3\times$ tighter uncertainty estimates than standard BQ.
Reaction-Diffusion PDE: On integrands arising from a PDE with multi-scale features (smooth vs. spiky regions), HBQ achieved a median error roughly $35\times$ lower than standard BQ at $N=512$ evaluations.

Significance and Scope
The paper positions HBQ as a principled Bayesian counterpart to the Genz–Malik adaptive cubature rule. While Genz–Malik uses embedded polynomial rules, HBQ replaces them with local GP surrogates, gaining principled uncertainty quantification and data-driven split decisions.

The authors claim that HBQ is particularly significant for expensive, low-to-moderate dimensional integration where integrands exhibit spatially varying behavior, such as model evidences over peaked likelihoods or simulator responses with multi-scale features. The method is designed to reduce to standard BQ when a single stationary model suffices, ensuring it does not sacrifice performance on stationary problems. The authors acknowledge limitations, noting that the method is restricted to axis-aligned hyperrectangles and kernels with closed-form mean embeddings, and that the convergence analysis is currently conditional on the final tree structure rather than the construction process itself.

Technical Summary: Hierarchical Bayesian Quadrature

More like this