Bayesian quadrature

Bayesian quadrature

Numerical integration

Let ${\displaystyle f:\mathcal{X}\to \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),\dots ,f(x_n)}$ at distinct locations ${\displaystyle x_1,\dots ,x_n}$ in ${\displaystyle \mathcal{X}}$ are used to estimate the integral of ${\displaystyle f}$ against a measure ${\displaystyle \nu }$: i.e. ${\displaystyle \nu [f]:={\int }_{\mathcal{X}}f(x)\nu (\mathrm{d}x).}$ Given weights ${\displaystyle w_1,\dots ,w_n\in \mathbb{R}}$, a quadrature rule is an estimator of ${\displaystyle \nu [f]}$ of the form $\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),\dots ,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}\to \mathbb{R}}$ and covariance function (or kernel function) ${\displaystyle k:\mathcal{X}\times \mathcal{X}\to \mathbb{R}}$. Then, the posterior distribution on ${\displaystyle f}$ is a Gaussian process with mean ${\displaystyle m_n:\mathcal{X}\to \mathbb{R}}$ and kernel ${\displaystyle k_n:\mathcal{X}\times \mathcal{X}\to \mathbb{R}}$ given by: $${\displaystyle m_n(x)=m(x)+k(x,X)k(X,X{)}^{-1}f(X)\text{and}{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)\text{and}\mathbb{V}[\nu [f]]=\nu \nu [k]-\nu [k(\cdot ,X)]k(X,X{)}^{-1}\nu [k(X,\cdot )].}$$ The function ${\displaystyle \nu [k(\cdot ,x)]={\int }_{\mathcal{X}}k(y,x)\nu (\mathrm{d}y)}$ is the kernel mean embedding of ${\displaystyle k}$ and ${\displaystyle \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 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),\dots ,f(x_n)}$ to set the kernel hyperparameters using, for example, maximum likelihood estimation. The estimation of kernel hyperparameters introduces adaptivity into Bayesian quadrature.

Example

Consider estimation of the integral $${\displaystyle \nu [f]={\int }_0^{1}f(x)\mathrm{d}x\approx 1.79\text{ of the function }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 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 (\frac{\sqrt{3}(x-1)}{\rho })(3+2\sqrt{3}\rho -3x)-\frac{1}{3}\exp (-\frac{\sqrt{3}x}{\rho })(3x+2\sqrt{3}\rho )}$$ $${\displaystyle \nu \nu [k]={\int }_0^1{\int }_0^{1}k(x,y)\mathrm{d}x\mathrm{d}y=\frac{2\rho }{3}[2\sqrt{3}-3\rho +\exp (-\frac{\sqrt{3}}{\rho })(\sqrt{3}+3\rho )].}$$ 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.
• 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.
• It permits more efficient use of information, e.g. jointly inferring multiple related quantities of interest or using active learning to reduce the required number of points.

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.