Research Radar — 2026-11-05

Methodology Review: Bayesian quadrature: put uncertainty on the computation you have not finished

English · About 10 minutes

Evidence: four primary papers inspected in full text. Numerical example independently derived; published experiments were not reproduced. Suggested omics transfers are speculative.

The intuition

Imagine an expensive simulator that can tell you the value of a response at any chosen input, but only after a long computation. You want its average over a specified range. Evaluating a few inputs leaves many possible curves between them, and those curves can enclose different areas. The missing information is not necessarily measurement noise: even a perfectly deterministic simulator can be incompletely evaluated. This is the problem Bayesian quadrature makes explicit.

The method places a probability model over possible functions, conditions it on the evaluations already obtained, and integrates that model. It returns a distribution over the integral, not just a fitted curve. A Gaussian process is convenient because integration is linear. The uncertainty describes the numerical calculation under a chosen prior; it is not automatically biological variability or a guarantee about repeated experiments. Hennig and colleagues explain why familiar quadrature rules can emerge as special cases of this perspective. [R2]

A small worked example

Here is a deliberately small example derived for this review. We want the area under an unknown function on the interval from zero to one. Assume its value at zero is known to be zero. Assign a zero-mean Gaussian-process prior with covariance \(k(s,t)=\min(s,t)\), where s and t are positions in that interval. This is the Brownian covariance: it fixes the left endpoint and allows increasing uncertainty away from it. It is a convenient teaching assumption, not a claim that a simulator physically follows Brownian motion.

Evaluate the right endpoint and obtain \(f(1)=1\). The conditional mean curve is the straight line joining the two known endpoints, so its area is one half. The integral has prior variance one third; the endpoint has variance one and covariance one half with the integral. Gaussian conditioning therefore gives posterior integral variance \(1/3-(1/2)^2=1/12\). Even though both endpoint values are known exactly, the interior can still bend. The nonzero variance represents those unevaluated bends.

Now evaluate the midpoint and obtain \(f(1/2)=1/4\). The conditional mean becomes two straight segments. Their areas are one sixteenth and five sixteenths, giving a total of three eighths. In the matrix notation introduced below, the covariance matrix is \(K=\begin{pmatrix}1/2&1/2\\1/2&1\end{pmatrix}\) and the integrated-covariance vector is \(q=(3/8,1/2)^{\top}\). Solving \(Kw=q\) gives weights \(w=(1/2,1/4)^{\top}\). Thus the new estimate is \(\widehat I=(1/2)(1/4)+(1/4)(1)=3/8\).

Its variance is \(1/3-q^{\top}w=1/48\), four times smaller than before. An independent check comes from the two conditional Brownian bridges: each interval of length one half contributes its length cubed divided by twelve, and their contributions add. This also exposes where the uncertainty lives. Another evaluation of an already known endpoint would add no information in this noiseless calculation, whereas evaluating the interior divides an uncertain interval.

Suppose the hidden function really was \(f(x)=x^2\). Its exact integral is one third. The first estimate errs by one sixth and the second by one twenty-fourth; neither error equals the posterior standard deviation. These are different quantities. The model variance describes a distribution of possible remaining errors, while this particular quadratic supplies one realized error. A single successful example cannot establish calibration, and the same endpoint observations are compatible with much more difficult functions.

Finally multiply the covariance kernel by a positive constant c. Every entry of K and q scales by c, so the weights and both area estimates stay unchanged, but the posterior variances multiply by c. You can therefore have exactly the same numerical answer with very different uncertainty statements. This elementary calculation is a useful warning against judging an uncertainty model solely by whether its mean prediction looks sensible.

The necessary mathematics

Let \(I=\int f(x)p(x)\,\mathrm{d}x\), where f is the integrand and p is a known normalized density. For n exact evaluations, define the column vector y by \(y_i=f(x_i)\) and the covariance matrix by \(K_{ij}=k(x_i,x_j)\). With a zero-mean Gaussian-process prior, set \(q_i=\int k(x_i,x)p(x)\,\mathrm{d}x\) and \(r=\iint k(x,x^{\prime})p(x)p(x^{\prime})\,\mathrm{d}x\,\mathrm{d}x^{\prime}\). The prime denotes a second integration variable; i and j index evaluation locations. [R1]

If these integrals exist and K is invertible, Gaussian conditioning yields:

Here N denotes a normal distribution, the superscript top denotes transpose, and the two arguments are its mean and variance. Solve a linear system rather than explicitly forming the inverse. For a fixed kernel and locations, the variance does not depend on y. Computing q and r must itself be tractable; otherwise another numerical approximation enters. [R1]

Assumptions and failure modes

The prior is substantive. Smoothness, scale, boundary behavior and unobserved narrow peaks all affect whether its uncertainty is credible. Our example silently enforced a zero left endpoint. A function with a nonzero value there would contradict the model before any algebra began. Check known constraints first, then test uncertainty against functions or simulations with accurately known integrals. [R2]

The simple Gaussian formula is not a positivity constraint. A likelihood is nonnegative, but an unconstrained Gaussian model can assign probability to negative values. Gunter and colleagues address this with a square-root transformation and linearization or moment matching. Those are additional modeling and approximation choices, not the same exact Gaussian conditioning calculation. Their computational results do not establish universal superiority over Monte Carlo. [R4]

Keep three budgets separate: evaluations of the expensive function, fitting and updating the surrogate, and checking its reliability. A method that uses fewer simulator calls can still take longer overall. Record both counts and elapsed time. If evaluations themselves contain stochastic simulation error, do not label them exact; our noiseless example would need an observation model. These are practical checks proposed here, not extra guarantees supplied by the integral posterior.

Verified applications

Osborne and colleagues apply active Bayesian quadrature to model evidence for detecting a damped Lyman-alpha absorber in quasar spectra, alongside synthetic benchmarks. The paper evaluates numerical error and uncertainty calibration under a limited evaluation budget. This is a concrete scientific use beyond toy integration, but it is astronomy rather than omics and does not validate biological simulator uncertainty. [R3]

Gunter and colleagues evaluate warped quadrature on probabilistic inference tasks, including Gaussian-process regression and classification. The extension is useful reading because it connects integrand constraints to the algorithm rather than treating the kernel as a neutral default. Reported benchmark advantages remain specific to their tasks, implementations and baselines. [R4]

Possible transfers — speculative

A speculative omics use would average an expensive mechanistic model response over a deliberately specified distribution of a small number of uncertain parameters. First build a low-dimensional test whose integral can be computed accurately. Compare quadrature estimates and intervals under multiple kernel assumptions, then repeat with hidden sharp-response regions. This would test numerical reliability before interpreting the average as a biological prediction. It is a proposed experiment, not a verified application from the readings.

A second speculative connection to automated experimentation is budget allocation inside the analysis software: decide which additional simulator evaluation most reduces uncertainty in a predefined integrated score. That differs from choosing which wet-lab experiment has the highest expected biological payoff. The target functional, uncertainty model and evaluation cost must be stated separately. A confidently computed average can still be the wrong scientific objective, and numerical precision cannot repair an unsupported mechanistic model.

What to try and what to read next

Reproduce the one-endpoint and midpoint calculations using fractions before building an optimizer. Print the estimate, posterior variance and actual error separately. Change the kernel amplitude and confirm that the estimate stays fixed while uncertainty changes. Then replace the quadratic by a function with a narrow bump between sampled locations. The central lesson is that uncertainty can arise from work you have not computed, but quantifying it requires assumptions that deserve testing. Start with the main reading for the integral posterior, the background paper for interpretation, and the two application papers for adaptive evaluation and positivity constraints.

Reading list

• [R1] Primary reading: Probabilistic Integration: A Role in Statistical Computation?. François-Xavier Briol, Chris J. Oates, Mark Girolami, Michael A. Osborne and Dino Sejdinovic (2019). Main reading: sections 2.1 and 4 explain the posterior and calibration. Journal 2019; first preprint 2015-12-03, v6 2017-10-18. Earlier title: Probabilistic Integration: A Role for Statisticians in Numerical Analysis?

• [R2] Related reading: Probabilistic numerics and uncertainty in computations. Philipp Hennig, Michael A. Osborne and Mark Girolami (2015). Conceptual background: numerical uncertainty is uncertainty about an incompletely computed quantity. Journal publication 2015-07-08; first and only arXiv version 2015-06-03. Focus on sections 2.1 and 2.3.

• [R3] Related reading: Active Learning of Model Evidence Using Bayesian Quadrature. Michael A. Osborne, David Duvenaud, Roman Garnett, Carl E. Rasmussen, Stephen J. Roberts and Zoubin Ghahramani (2012). Application and adaptive evaluation: inspect section 6 for the astronomy experiment and its restricted evaluation budget. NIPS 2012 proceedings publication; no unverified precise publication day or separate preprint alias is asserted.

• [R4] Related reading: Sampling for Inference in Probabilistic Models with Fast Bayesian Quadrature. Tom Gunter, Michael A. Osborne, Roman Garnett, Philipp Hennig and Stephen J. Roberts (2014). Extension and caveat: WSABI handles nonnegative integrands through a square-root transformation, with approximation choices that matter. NIPS 2014; first and only preprint version 2014-11-03. Read sections 3 and 5 rather than assuming every reported uncertainty is exact.