Bayesian quadrature

Bayesian quadrature with a Gaussian process conditioned on ${\displaystyle n=0,3,\ {\text{and}}\ 8}$ evaluations of the integrand (shown in black). Shaded areas in the left column illustrate the marginal standard deviations. The right figure shows the prior (${\displaystyle n=0}$) and posterior (${\displaystyle n=3,8}$) Gaussian distribution over the value of the integral, as well as the true solution.

Bayesian quadrature ^{[1] [2] [3] [4] [5]} is a method for approximating intractable integration problems. It falls within the class of probabilistic numerical methods. Bayesian quadrature views numerical integration as a Bayesian inference task, where function evaluations are used to estimate the integral of that function. For this reason, it is sometimes also referred to as "Bayesian probabilistic numerical integration" or "Bayesian numerical integration". The name "Bayesian cubature" is also sometimes used when the integrand is multi-dimensional. A potential advantage of this approach is that it provides probabilistic uncertainty quantification for the value of the integral.

Bayesian quadrature

Numerical integration

Let ${\displaystyle f:{\mathcal {X}}\rightarrow \mathbb {R} }$ be a function defined on a domain ${\displaystyle {\mathcal {X}}}$ (where typically ${\displaystyle {\mathcal {X}}\subseteq \mathbb {R} ^{d}}$). In numerical integration, function evaluations ${\displaystyle f(x_{1}),\ldots ,f(x_{n})}$ at distinct locations ${\displaystyle x_{1},\ldots ,x_{n}}$ in ${\displaystyle {\mathcal {X}}}$ are used to estimate the integral of ${\displaystyle f}$ against a measure ${\displaystyle \nu }$: i.e. ${\displaystyle \textstyle \nu [f]:=\int _{\mathcal {X}}f(x)\nu (\mathrm {d} x).}$ Given weights ${\displaystyle w_{1},\ldots ,w_{n}\in \mathbb {R} }$, a quadrature rule is an estimator of ${\displaystyle \nu [f]}$ of the form $\textstyle {\hat {\nu }}[f]:=\sum _{i=1}^{n}w_{i}f(x_{i}).$

Bayesian quadrature consists of specifying a prior distribution over ${\displaystyle f}$, conditioning this prior on ${\displaystyle f(x_{1}),\ldots ,f(x_{n})}$ to obtain a posterior distribution ${\displaystyle f}$, then computing the implied posterior distribution on ${\displaystyle \nu [f]}$. The name "quadrature" comes from the fact that the posterior mean on ${\displaystyle \nu [f]}$ sometimes takes the form of a quadrature rule whose weights are determined by the choice of prior.

Bayesian quadrature with Gaussian processes

The most common choice of prior distribution for ${\displaystyle f}$ is a Gaussian process as this permits conjugate inference to obtain a closed-form posterior distribution on ${\displaystyle \nu [f]}$. Suppose we have a Gaussian process with prior mean function ${\displaystyle m:{\mathcal {X}}\rightarrow \mathbb {R} }$ and covariance function (or kernel function) ${\displaystyle k:{\mathcal {X}}\times {\mathcal {X}}\rightarrow \mathbb {R} }$. Then, the posterior distribution on ${\displaystyle f}$ is a Gaussian process with mean ${\displaystyle m_{n}:{\mathcal {X}}\rightarrow \mathbb {R} }$ and kernel ${\displaystyle k_{n}:{\mathcal {X}}\times {\mathcal {X}}\rightarrow \mathbb {R} }$ given by: $${\displaystyle m_{n}(x)=m(x)+k(x,X)k(X,X)^{-1}f(X)\qquad {\text{and}}\qquad k_{n}(x,y)=k(x,y)-k(x,X)k(X,X)^{-1}k(X,y).}$$ where ${\displaystyle (k(X,X))_{ij}=k(x_{i},x_{j})}$, ${\displaystyle (f(X))_{i}=f(x_{i})}$, ${\displaystyle (k(\cdot ,X))_{i}=k(\cdot ,x_{i})}$ and ${\displaystyle (k(X,\cdot ))_{i}=k(x_{i},\cdot )}$.

Furthermore, the posterior distribution on ${\displaystyle \nu [f]}$ is a univariate Gaussian distribution with mean ${\displaystyle \mathbb {E} [\nu [f]]}$ and variance ${\displaystyle \mathbb {V} [\nu [f]]}$ given by $${\displaystyle \mathbb {E} [\nu [f]]=\nu [m]+\nu [k(\cdot ,X)]k(X,X)^{-1}f(X)\qquad {\text{and}}\qquad \mathbb {V} [\nu [f]]=\nu \nu [k]-\nu [k(\cdot ,X)]k(X,X)^{-1}\nu [k(X,\cdot )].}$$ The function ${\displaystyle \textstyle \nu [k(\cdot ,x)]=\int _{\mathcal {X}}k(y,x)\nu (\mathrm {d} y)}$ is the kernel mean embedding of ${\displaystyle k}$ and ${\displaystyle \textstyle \nu \nu [k]=\int _{\mathcal {X}}k(x,y)\nu (dx)\nu (\mathrm {d} y)}$ denotes the integral of ${\displaystyle k}$ with respect to both inputs. In particular, note that the posterior mean is a quadrature rule with weights ${\displaystyle \textstyle w_{i}=(\nu [k(\cdot ,X)]k(X,X)^{-1})_{i}.}$ and the posterior variance provides a quantification of the user's uncertainty over the value of ${\displaystyle \nu [f]}$.

In more challenging integration problems, where the prior distribution cannot be relied upon as a meaningful representation of epistemic uncertainty, it is necessary to use the data ${\displaystyle f(x_{1}),\ldots ,f(x_{n})}$ to set the kernel hyperparameters using, for example, maximum likelihood estimation. The estimation of kernel hyperparameters introduces adaptivity into Bayesian quadrature. ^{[6] [7]}

Example

Illustration of Bayesian quadrature for estimating ${\displaystyle \textstyle \nu [f]=\int _{0}^{1}f(x)\,\mathrm {d} x}$ where ${\displaystyle \textstyle f(x)=(1+x^{2})\sin(5\pi x)+8/5}$. The posterior distribution (blue) concentrates on the true integral when more data (the red points) is obtained of the integrand ${\displaystyle f}$.

Consider estimation of the integral $${\displaystyle \nu [f]=\int _{0}^{1}f(x)\,\mathrm {d} x\approx 1.79\quad {\text{ of the function }}\quad f(x)=(1+x^{2})\sin(5\pi x)+{\frac {8}{5}}}$$ using a Bayesian quadrature rule based on a zero-mean Gaussian process prior with the Matérn covariance function of smoothness ${\displaystyle 3/2}$ and correlation length ${\displaystyle \rho =1/5}$. This covariance function is ${\displaystyle \textstyle k(x,y)=(1+{\sqrt {3}}\,|x-y|/\rho )\exp(\!-{\sqrt {3}}\,|x-y|/\rho ).}$ It is straightforward (though tedious) to compute that $${\displaystyle \nu [k(\cdot ,x)]=\int _{0}^{1}k(y,x)\,\mathrm {d} y={\frac {4\rho }{\sqrt {3}}}-{\frac {1}{3}}\exp {\bigg (}{\frac {{\sqrt {3}}(x-1)}{\rho }}{\bigg )}{\big (}3+2{\sqrt {3}}\,\rho -3x{\big )}-{\frac {1}{3}}\exp {\bigg (}-{\frac {{\sqrt {3}}\,x}{\rho }}{\bigg )}{\big (}3x+2{\sqrt {3}}\,\rho {\big )}}$$$${\displaystyle \nu \nu [k]=\int _{0}^{1}\int _{0}^{1}k(x,y)\,\mathrm {d} x\,\mathrm {d} y={\frac {2\rho }{3}}{\Bigg [}2{\sqrt {3}}-3\rho +\exp {\bigg (}\!-{\frac {\sqrt {3}}{\rho }}{\bigg )}{\big (}{\sqrt {3}}+3\rho {\big )}{\Bigg ]}.}$$ Convergence of the Bayesian quadrature point estimate ${\displaystyle \mathbb {E} [\nu [f]]}$ and concentration of the posterior mass, as quantified by ${\displaystyle \mathbb {V} [\nu [f]]}$, around the true integral ${\displaystyle \nu [f]}$ as ${\displaystyle f}$ is evaluated at more and more points is displayed in the accompanying animation.

Advantages and disadvantages

Since Bayesian quadrature is an example of probabilistic numerics, it inherits certain advantages compared with traditional numerical integration methods:

• It allows uncertainty to be quantified and propagated through all subsequent computations to explicitly model the impact of numerical error. ^[8]
• It provides a principled way to incorporate prior knowledge by using a judicious choice of prior distributions for ${\displaystyle f}$, which may be more sophisticated compared to the standard Gaussian process just described. ^[7]
• It permits more efficient use of information, e.g. jointly inferring multiple related quantities of interest ^[9] or using active learning to reduce the required number of points. ^[10]

Despite these merits, Bayesian quadrature methods possess the following limitations:

• Although the Bayesian paradigm allows a principled treatment of the quantification of uncertainty, posterior inference over ${\displaystyle \nu [f]}$ is not always tractable, thus requiring a second-level estimation. E.g. for Bayesian quadrature with Gaussian processes, the kernel mean embedding ${\displaystyle \nu [k(\cdot ,x)]}$ has no closed-form expression for a general kernel ${\displaystyle k}$ and measure ${\displaystyle \nu }$.

• The computational cost of Bayesian quadrature methods based on Gaussian processes is in general ${\displaystyle {\mathcal {O}}(n^{3})}$ due to the cost of inverting ${\displaystyle n\times n}$ matrices, which may defy their applications to large-scale problems.

Algorithmic design

Prior distributions

The most commonly used prior for ${\displaystyle f}$ is a Gaussian process prior. This is mainly due to the advantage provided by Gaussian conjugacy and the fact that Gaussian processes can encode a wide range of prior knowledge including smoothness, periodicity and sparsity through a careful choice of prior covariance. However, a number of other prior distributions have also been proposed. This includes multi-output Gaussian processes, ^[9] which are particularly useful when tackling multiple related numerical integration tasks simultaneously or sequentially, and tree-based priors such as Bayesian additive regression trees, ^[10] which are well suited for discontinuous ${\displaystyle f}$. Additionally, Dirichlet processes priors have also been proposed for the integration measure ${\displaystyle \nu }$. ^[11]

Point selection

The points ${\displaystyle x_{1},\ldots ,x_{n}}$ are either considered to be given, or can be selected so as to ensure the posterior on ${\displaystyle \nu [f]}$ concentrates at a faster rate. One approach consists of using point sets from other quadrature rules. For example, taking independent and identically distributed realisations from ${\displaystyle \nu }$ recovers a Bayesian approach to Monte Carlo, ^[3] whereas using certain deterministic point sets such as low-discrepancy sequences or lattices recovers a Bayesian alternative to quasi-Monte Carlo. ^{[4] [12]} It is of course also possible to use point sets specifically designed for Bayesian quadrature; see for example the work of ^[13] who exploited symmetries in point sets to obtain scalable Bayesian quadrature estimators. Alternatively, points can also be selected adaptively following principles from active learning and Bayesian experimental design so as to directly minimise posterior uncertainty, ^{[14] [15]} including for multi-output Gaussian processes. ^[16]

Kernel mean and initial error

One of the challenges when implementing Bayesian quadrature is the need to evaluate the function ${\displaystyle \nu [k(\cdot ,x)]}$ and the constant ${\displaystyle \nu \nu [k]}$. The former is commonly called the kernel mean, and is a quantity which is key to the computation of kernel-based distances such as the maximum mean discrepancy. The latter is commonly called the initial error since it provides an upper bound on the integration error before any function values are observed. Unfortunately, the kernel mean and initial error can only be computed for a small number of ${\displaystyle (k,\nu )}$ pairs; see for example Table 1 in. ^[4]

Theory

There have been a number of theoretical guarantees derived for Bayesian quadrature. These usually require Sobolev smoothness properties of the integrand, ^{[4] [17] [18]} although recent work also extends to integrands in the reproducing kernel Hilbert space of the Gaussian kernel. ^[19] Most of the results apply to the case of Monte Carlo or deterministic grid point sets, but some results also extend to adaptive designs. ^{[20] [21]}

Software

• ProbNum: Probabilistic numerical methods in Python, including a Bayesian quadrature implementation.
• Emukit: Emulation and decision making under uncertainty in Python.
• QMCPy: Bayesian quadrature with QMC point sets in Python.

References

1. Diaconis, P. (1988). "Bayesian Numerical Analysis". Statistical Decision Theory and Related Topics IV: 163–175. doi:10.1007/978-1-4613-8768-8_20 (inactive 1 November 2024). ISBN   978-1-4613-8770-1.
2. O’Hagan, A. (2002). "Bayes–Hermite quadrature". Journal of Statistical Planning and Inference (29): 245–260.
3. Rasmussen, C.; Ghahramani, Z. (2002). "Bayesian Monte Carlo". Neural Information Processing Systems: 489–496.
4. Briol, F.-X.; Oates, C. J.; Girolami, M.; Osborne, M. A.; Sejdinovic, D. (2019). "Probabilistic integration: A role in statistical computation? (with discussion and rejoinder)". Statistical Science. 34 (1): 1–22.
5. Hennig, P.; Osborne, M. A.; Kersting, H. P. (2022). Probabilistic Numerics (PDF). Cambridge University Press. pp. 63–122. ISBN   978-1107163447.
6. Jagadeeswaran, R.; Hickernell, Fred J. (2019-09-10). "Fast automatic Bayesian cubature using lattice sampling". Statistics and Computing. 29 (6): 1215–1229. arXiv: 1809.09803 . doi:10.1007/s11222-019-09895-9. ISSN   0960-3174. S2CID   119709309.
7. Fisher, Matthew; Oates, Chris; Powell, Catherine; Teckentrup, Aretha (2020-06-03). "A Locally Adaptive Bayesian Cubature Method". International Conference on Artificial Intelligence and Statistics. PMLR: 1265–1275. arXiv: 1910.02995 .
8. Cockayne, Jon; Oates, Chris; Sullivan, Tim; Girolami, Mark (2017). "Probabilistic numerical methods for PDE-constrained Bayesian inverse problems". Bayesian Inference and Maximum Entropy Methods in Science and Engineering. AIP Conference Proceedings. 1853 (1). Author(s): 060001. arXiv: 1701.04006 . Bibcode:2017AIPC.1853f0001C. doi:10.1063/1.4985359. hdl:10044/1/66123. S2CID   17349210.
9. Xi, X.; Briol, F.-X.; Girolami, M. (2018). "Bayesian quadrature for multiple related integrals". International Conference on Machine Learning: 8533–8564. arXiv: 1801.04153 .
10. Zhu, H.; Liu, X.; Kang, R.; Shen, Z.; Flaxman, S.; Briol, F.-X. (2020). "Bayesian probabilistic numerical integration with tree-based models". Neural Information Processing Systems: 5837–5849. arXiv: 2006.05371 .
11. Oates, C. J.; Niederer, S.; Lee, A.; Briol, F.-X.; Girolami, M. (2017). "Probabilistic models for integration error in the assessment of functional cardiac models". Neural Information Processing Systems: 110–118. arXiv: 1606.06841 .
12. Jagadeeswaran, R.; Hickernell, F. J. (2019). "Fast automatic Bayesian cubature using lattice sampling". Statistics and Computing. 29 (6): 1215–1229. arXiv: 1809.09803 . doi:10.1007/s11222-019-09895-9. S2CID   119709309.
13. Karvonen, T.; Särkkä, S. (2018). "Fully symmetric kernel quadrature". SIAM Journal on Scientific Computing. 40 (2): 697–720. arXiv: 1703.06359 . Bibcode:2018SJSC...40A.697K. doi:10.1137/17M1121779. S2CID   9707762.
14. Gunter, T.; Garnett, R.; Osborne, M. A.; Hennig, P.; Roberts, S. (2014). "Sampling for inference in probabilistic models with fast Bayesian quadrature". Neural Information Processing Systems: 2789–2797. arXiv: 1411.0439 .
15. Briol, F.-X.; Oates, C. J.; Girolami, M.; Osborne, M. A. (2015). "Frank-Wolfe Bayesian quadrature: Probabilistic integration with theoretical guarantees". Neural Information Processing Systems: 1162–1170. arXiv: 1506.02681 .
16. Gessner, A.; Gonzalez, J.; Mahsereci, M. (2019). "Active multi-information source Bayesian quadrature". Uncertainty in Artificial Intelligence. arXiv: 1903.11331 .
17. Kanagawa, M.; Sriperumbudur, B. K.; Fukumizu, K. (2020). "Convergence analysis of deterministic kernel-based quadrature rules in misspecified settings". Foundations of Computational Mathematics. 20: 155–194. arXiv: 1709.00147 . doi: 10.1007/s10208-018-09407-7 . S2CID   11717907.
18. Wynne, G.; Briol, F.-X.; Girolami, M. (2020). "Convergence guarantees for Gaussian process means with misspecified likelihoods and smoothness". Journal of Machine Learning Research. 22 (123): 1–40. arXiv: 2001.10818 .
19. Karvonen, T.; Oates, C. J.; Girolami, M. (2021). "Integration in reproducing kernel Hilbert spaces of Gaussian kernels". Mathematics of Computation. 90 (331): 2209–2233. arXiv: 2004.12654 . doi: 10.1090/mcom/3659 . S2CID   216552869.
20. Kanagawa, M.; Hennig, P. (2019). "Convergence guarantees for adaptive Bayesian quadrature methods". Neural Information Processing Systems: 6237–6248. arXiv: 1905.10271 .
21. Adachi, M.; Satoshi, H.; J{\o}rgensen, M.; Oberhauser, H.; Osborne, M. A. (2022). "Fast Bayesian inference with batch Bayesian quadrature via kernel recombination". Neural Information Processing Systems: 16533–16547. arXiv: 2206.04734 .