1 Introduction↩︎

Numerical computation of integrals, often representing expectations or normalization constants, is a pervasive practical problem throughout machine learning, statistics, scientific computing, and engineering. In the form studied in this survey, numerical integration consists in approximating the definite integral \[I_P(f) = \int_Df({\boldsymbol{x}}) \dif P({\boldsymbol{x}})\] of a real-valued integrand function \(f\) with respect to a probability measure \(P\) on a set \(D\subseteq \mathbb{R}^d\). The approximation typically takes the form of a quadrature rule \(\sum_{i=1}^N w_i f({\boldsymbol{x}}_i)\), which is a weighted sum of integrand evaluations at some nodes \({\boldsymbol{x}}_i \in D\). Numerical integration can be done on arbitrary non-Euclidean domains and can also incorporate information other than point evaluations, generalisations that we discuss in 2.5. Over the past two centuries a vast number of numerical integration techniques have been developed from different starting points and for different types of integration problems. The most popular of these are classical Gaussian quadratures that are constructed using polynomial exactness criteria and admit an elegant mathematical theory [1]; easy-to-use and widely applicable Monte Carlo integration based on random sampling that is, no doubt, the most popular approach to approximating an integral [2]; and quasi-Monte Carlo methods that use low-discrepancy sequences and can compute high-dimensional integrals effectively [3]. Other approaches include sparse grid methods, Clenshaw–Curtis quadrature, the trapezoidal rule, and the topic of this survey, Bayesian quadrature.

Although its history can be traced further back, Bayesian quadrature was not popularized until the late 1980s by Persi Diaconis and Anthony O’Hagan [4], [5]. Since then, there has been a steady stream of interest and contributions from machine learning and statistics communities [6]–[9] and, more recently, from the point of view of numerical analysis and approximation theory [10]–[12]. Bayesian quadrature falls within probabilistic numerics [13], its defining characteristic being the interpretation of numerical integration as a statistical inference problem to which the Bayesian paradigm and methods can be brought to bear.

In Bayesian quadrature, a stochastic process prior placed on the integrand is conditioned on the “data” consisting of integrand evaluations. The mean of the resulting posterior distribution provides a point estimate for the integral, while the spread of the posterior quantifies uncertainty (see Figure 1). That explicit prior information, such as

smoothness or structural properties, is easy to encode both conceptually and in practice by selection of an appropriate prior for the integrand distinguishes Bayesian quadrature from other numerical integration techniques. Priors that are “correct” or “good” result in fast rates of convergence and reliable quantification of uncertainty (see Figure 2). Non-Bayesian integration methods too encode various types of prior information, but typically in a non-systematic and implicit manner. For example, Gaussian quadratures implicitly assume that the integrand is well approximated by polynomials; the trapezoidal rule that a sum of trapezoids approximates the integral; and Monte Carlo encodes no assumptions whatsoever besides square-integrability. While quasi-Monte Carlo methods use node placement to encode certain assumptions, Bayesian quadrature can couple any nodes with any prior information, making it extremely versatile. That any nodes proposed in the vast literature on numerical integration can be used in Bayesian quadrature is one of its great advantages. Another distinguishing feature is the statistically principled uncertainty quantification that Bayesian quadrature provides: the spread of the posterior tells one how much to trust the integral estimate. Note that this uncertainty quantification arises from the prior and is wholly distinct in character from the sample variance of Monte Carlo.

Much methodology and theory has been developed for Bayesian quadrature over the years (worthy mentions include [14]; [15]; and [11]), and applications have ranged from computer graphics [16] and engineering [17], [18] to cardiac and tsunami modelling [19], [20]. Despite a wealth of such contributions, no comprehensive treatment of Bayesian quadrature as a whole has ever been published. Perhaps partly for this reason, confusion continues to plague the field. For example, “Bayesian quadrature” is often used to refer exclusively to a single method (or some collection of methods) to tackle integration problems of a very specific form. Although comparisons of Bayesian quadrature methods to other numerical integration techniques abound in the literature, no systematic empirical study in a controlled environment on the effect of the multitude of choices on modelling, inference, and sampling that go into designing a Bayesian quadrature method has been undertaken. This is a pity, for the richness of these choices is one of the main attractions (and often a great source of frustration in application) of Bayesian quadrature. The purpose of this survey is to fill this gap and provide a comprehensive survey of Bayesian quadrature that covers both the theory and practice of Bayesian quadrature.

1.1 Contents↩︎

There are six main elements that make up the contents of this survey.

1.1.0.0.1 I — 2: Foundations.

This chapter contains a foundational introduction to Bayesian quadrature. In particular, the chapter contains what appears to be the first attempt at properly defining Bayesian quadrature (see Definition 2.1) in a way that encompasses the modern usage of the term in machine learning literature. In our parlance a Bayesian quadrature refers to numerical integration based on Gaussian process modelling. That is, a Bayesian quadrature method is any numerical integration method that computes or approximates a posterior distribution for the integral \(I_P(\mathsf{f})\) of a potentially transformed Gaussian process \(\mathsf{f}\) (i.e., \(\mathsf{f}= \varphi \circ \mathsf{g}\) for a mapping \(\varphi\) and a Gaussian process \(\mathsf{g}\)). The definition gives rise to a natural division of Bayesian quadrature methods in (a) conjugate methods for which \(\mathsf{f}\) is a Gaussian process (i.e., \(\varphi\) is affine) and (b) non-conjugate methods for which \(\mathsf{f}\) is not a Gaussian process (i.e., \(\varphi\) is non-affine). This anticipates the development of more detailed taxonomy in 3. We also discuss how the theory of reproducing kernel Hilbert spaces can be used to interpret Bayesian quadrature, an interpretation which underlies most theoretical analysis, and provide high-level theory that is independent of the underlying Gaussian process model.

1.1.0.0.2 II — 3: Taxonomy.

This chapter contains an extensive taxonomy for Bayesian quadrature that provides a well-defined framework for classifying and categorizing concrete algorithms and implementations. Inspiration for the taxonomy comes from the literature on global optimization [21]. We classify Bayesian quadrature methods along three main axes based on (i) what kind of model they employ, (ii) how they approach inference, and (iii) how they select, or sample, the nodes \({\boldsymbol{x}}_i\) at which the integrand is evaluated. 8 contains an overview of the taxonomy and its most important labels. The taxonomy illustrates that Bayesian quadrature is an umbrella term and that there is a great variety of different Bayesian quadrature methods that differ significantly from one another.

1.1.0.0.3 III — 4: Practice.

Like all statistical and numerical methods, Bayesian quadrature methods are faced with practical problems that limit their applicability. For Bayesian quadrature, the dominant problems are those of computational cost and computation of kernel means, as one needs to solve linear systems of \(N\) equations and compute integrals that involve the covariance function of the Gaussian process model. Although every practitioner of Bayesian quadrature is inevitably confronted by these problems, they are rarely emphasised or treated in a systematic manner. 4 reviews what we consider the most important practical problems that a practitioner needs to solve in order to implement a Bayesian quadrature method and offer some solutions. One should recognise that Bayesian quadrature can never hope to be as flexible or general as, for example, Monte Carlo integration. This is understandable, for Bayesian quadrature is based on principled statistical modelling, which must always have a cost.

1.1.0.0.4 IV — 5: Theoretical guarantees.

Various theoretical convergence guarantees, given in terms of upper bounds on the integration error under certain assumptions, have appeared in the literature [9], [11], [22]. In 5, we provide somewhat more practical versions of these guarantees by (i) making the dependency of the bounds explicit on those parameters of the covariance function that are commonly estimated from the data and by (ii) including the effect regularisation via a nugget term has. The proofs, which appear in 7, are based on results in the scattered data approximation literature and do not materially differ from those of the existing results.

1.1.0.0.5 V — 6: Empirical illustrations.

Here we put the taxonomy and theory from [sec:taxonomy,sec:guarantees] in action by conducting extensive and controlled numerical experiments along the three main axes of the taxonomy. That is, we systematically study how the model, inference methodology, and sampling scheme affect the accuracy and calibration of a Bayesian quadrature method, as well as how the different choices interact. In the experiments we use a collection of test integrands whose properties, such as differentiability, are well-understood. The experiments validate some of the theoretical guarantees. While almost every article on Bayesian quadrature contains an empirical study, these typically focus on comparing the new method to other candidates, such as Monte Carlo or quasi-Monte Carlo, in realistic—and hence opaque—integration problems that do not easily lend themselves to untangling the reasons behind one method outperforming the other. The purpose of this chapter is not to demonstrate the superiority of Bayesian quadrature in general or any particular method to alternative numerical methods but to highlight the effect of different approaches and choices within Bayesian quadrature.

1.1.0.0.6 VI — Comprehensive bibliography.

We do not conduct a systematic literature review. However, as a survey article such as this is naturally equipped with a long bibliography, we could not resist taking this to extremes and including every reference, no matter how obscure, on Bayesian quadrature or equivalent approaches known to us in machine learning, statistics, engineering, and mathematics literature. Consequently, in most cases we make little to no effort at providing more than a modicum of context for the references.

1.2 Notational conventions↩︎

The set of non-negative integers is denoted \(\mathbb{N}\). We strive to follow the following type conventions: objects in italics are scalars or scalar-valued functions (\(x\) or \(f\)); sans serif is used to indicate random variables and stochastic processes (\(\mathsf{f}\)); boldface objects are vectors (if lowercase; \({\boldsymbol{x}}\)) or matrices (if uppercase; \({\boldsymbol{K}}\)); and boldface sans serif stands for a random vector (\({\boldsymbol{\mathsf{f}}}\)). In 7.3.2, we use \(\lesssim\) to denote asymptotic inequality. This notation is equivalent to big \(O\) notation.

2 Bayesian quadrature↩︎

This chapter is an introduction to Bayesian quadrature, which we define (see Definition 2.1) as any method for computing or approximating the integral of a potentially transformed Gaussian process \(\mathsf{f}\) conditioned on a number of evaluations of the integrand \(f\). We further divide Bayesian quadrature methods to

• conjugate methods when the posterior over the integral is normally distributed and

• non-conjugate methods when the posterior is non-normal.

A detailed classification of the existing Bayesian quadrature methods is given in 3. To precede the construction of Bayesian quadrature in 2.2 we include a short review of Gaussian processes in 2.1. Moreover, in 2.3 we review the equivalence between conjugate Bayesian quadrature and worst-case optimal quadrature rules in reproducing kernel Hilbert spaces. This equivalence is crucial for deriving the convergence guarantees,as well as many other theoretical results, that we present in 5. Throughout this chapter it is assumed that both the set of nodes \(X\) at which the integrand is to be evaluated and the transformed Gaussian process \(\mathsf{f}\) (in particular, its covariance function \(k\)) used to model the integrand are fixed; different ways to select them are reviewed in [sec:taxonomy,sec:practical-issues].

2.1 Gaussian processes↩︎

[23] provide a standard introduction to Gaussian processes for machine learning. For more theoretical expositions, see [24] and [25]. Let \(k\colon D\times D\to \mathbb{R}\) be a symmetric positive-semidefinite covariance function (or kernel), which is to say that the kernel Gram matrix \[{\boldsymbol{K}}_{XX} = (k({\boldsymbol{x}}_i, {\boldsymbol{x}}_j))_{i,j=1}^N \in \mathbb{R}^{N \times N}\] is positive-semidefinite for any \(N \in \mathbb{N}\) and any \(X= \{{\boldsymbol{x}}_1, \ldots, {\boldsymbol{x}}_N\} \subseteq D\). Throughout this survey we tacitly assume that \(X\) are such that the Gram matrix is positive-definite and hence invertible. Most covariance function that we consider are positive-definite, in that their Gram matrix is positive-definite whenever the points \(X\) are pairwise distinct.

Let \(m\colon D\to \mathbb{R}\) be a function. A Gaussian process \(\mathsf{g}\) with mean \(m\) and covariance function \(k\) is denoted \[\mathsf{g}\sim \mathrm{GP}(m, k).\] Suppose that one wishes to do regression and is given a dataset \[{\cal D}= \{ ({\boldsymbol{x}}_1, y_1), \ldots ({\boldsymbol{x}}_N, y_N) \}\] consisting of some noiseless observations \({\boldsymbol{y}}= (y_1, \ldots, y_N) \in \mathbb{R}^N\) at pairwise distinct nodes \(X = \{{\boldsymbol{x}}_1, \ldots, {\boldsymbol{x}}_N\} \subseteq D\). The mean and covariance of the conditional Gaussian process, \(\mathsf{g}\mid {\cal D}\), are \[\begin{align} \label{eq:gp-posterior-mean} m_{\cal D}({\boldsymbol{x}}) = \mathbb{E}[ \mathsf{g}({\boldsymbol{x}}) \mid {\cal D}] = m({\boldsymbol{x}}) + {\boldsymbol{k}}_X({\boldsymbol{x}})^{\intercal}{\boldsymbol{K}}_{XX}^{-1} ( {\boldsymbol{y}}- {\boldsymbol{m}}_X) \end{align}\tag{1}\] and \[\begin{align} \label{eq:gp-posterior-covariance} k_{\cal D}( {\boldsymbol{x}}, {\boldsymbol{x}}' ) &= \mathrm{Cov}[ \mathsf{g}({\boldsymbol{x}}), \mathsf{g}({\boldsymbol{x}}') \mid {\cal D}] \\ &= k({\boldsymbol{x}}, {\boldsymbol{x}}') - {\boldsymbol{k}}_X({\boldsymbol{x}})^{\intercal}{\boldsymbol{K}}_{XX}^{-1} {\boldsymbol{k}}_X({\boldsymbol{x}}'), \end{align}\tag{2}\] where \({\boldsymbol{m}}_X= (m({\boldsymbol{x}}_1), \ldots, m({\boldsymbol{x}}_N)) \in \mathbb{R}^N\) and \({\boldsymbol{k}}_X\colon D \to \mathbb{R}^N\) is a vector-valued function defined as \[{\boldsymbol{k}}_X({\boldsymbol{x}}) = ( k({\boldsymbol{x}}, {\boldsymbol{x}}_1), \ldots, k({\boldsymbol{x}}, {\boldsymbol{x}}_N) ) \in \mathbb{R}^N.\]

A covariance function is stationary if it only depends on the difference of the inputs and isotropic if it depends only on their distance. That is, a covariance \(k\) is stationary if there is a function \(\Phi_k \colon D\to \mathbb{R}\) such that \(k({\boldsymbol{x}}, {\boldsymbol{x}}') = \Phi_k({\boldsymbol{x}}- {\boldsymbol{x}}')\) for all \({\boldsymbol{x}}, {\boldsymbol{x}}' \in D\) and isotropic if \(\Phi_k\colon [0, \infty) \to \mathbb{R}\) and \(k({\boldsymbol{x}}, {\boldsymbol{x}}') = \Phi_k( \norm[0]{ {\boldsymbol{x}}- {\boldsymbol{x}}'} )\). Commonly used covariance functions on any \(D \subseteq \mathbb{R}^d\) include the square exponential \[\label{eq:se-kernel} k({\boldsymbol{x}}, {\boldsymbol{x}}') = \sigma^2 \exp\bigg( \! - \frac{\norm[0]{{\boldsymbol{x}}- {\boldsymbol{x}}'}^2}{2 \ell^2} \bigg)\tag{3}\] and covariance functions of the Matérn class \[\label{eq:matern-kernel} k({\boldsymbol{x}}, {\boldsymbol{x}}') = \sigma^2 \frac{2^{1-\nu}}{\Gamma(\nu)} \bigg( \frac{\sqrt{2\nu} \norm[0]{{\boldsymbol{x}}- {\boldsymbol{x}}'}}{\ell} \bigg)^\nu K_\nu \bigg( \frac{\sqrt{2\nu} \norm[0]{{\boldsymbol{x}}- {\boldsymbol{x}}'}}{\ell} \bigg),\tag{4}\] where \(\Gamma\) denotes the gamma function and \(K_\nu\) the modified Bessel function of the second kind of order \(\nu > 0\).

Both the square exponential and all Matérn covariance functions are isotropic. Product Matérn covariance functions, which are defined as dimension-wise products of one-dimensional Matérn kernels (and are therefore stationary), \[\label{eq:product-matern-kernel} k({\boldsymbol{x}}, {\boldsymbol{x}}') = \sigma^2 \prod_{i=1}^d \frac{2^{1-\nu}}{\Gamma(\nu)} \bigg( \frac{\sqrt{2\nu} \abs[0]{x_i - x_i'}}{\ell} \bigg)^\nu K_\nu \bigg( \frac{\sqrt{2\nu} \abs[0]{x_i - x_i'}}{\ell} \bigg),\tag{5}\] are also often used. Note that the square exponential kernel 3 is the product of one-dimensional square-exponential kernels because \(\norm[0]{{\boldsymbol{x}}}^2 = x_1^2 + \cdots + x_d^2\).¹ The Brownian motion covariance function \[\label{eq:brownian-kernel} k(x, x') = \min\{x, x'\},\tag{6}\] which is defined on \(D = [0, T] \subseteq \mathbb{R}\) for any \(T > 0\) is important due to its connections to spline interpolation and the trapezoidal rule. The Brownian motion kernel is not stationary. 1 shows three different covariance functions and GPs. Although most covariance functions used in Bayesian quadrature are stationary, there are several exceptions besides the Brownian motion kernel. For example, [27] use a piecewise constant covariance function.

2.2 Construction of Bayesian quadrature↩︎

Let \(\varphi\colon \mathbb{R}\to \mathbb{R}\) be a function. Bayesian quadrature methods proceed by taking a Gaussian process prior \(\mathsf{g}\sim \mathrm{GP}(m, k)\) with a prior mean \(m\) and covariance function \(k\) and defining a prior \(\mathsf{f}\) for the integrand \(f\) by mapping \(\mathsf{g}\) through \(\varphi\): \[\label{eq:fGP-transformation} \mathsf{f}= \varphi\circ \mathsf{g}.\tag{7}\] Let the dataset consist of noiseless evaluations of the integrand \(f\) at nodes \({\boldsymbol{x}}_1, \ldots, {\boldsymbol{x}}_N\): \[{\cal D}= \{ ({\boldsymbol{x}}_1, f({\boldsymbol{x}}_1)), \ldots, ({\boldsymbol{x}}_N, f({\boldsymbol{x}}_N)) \}.\] We call any method that computes or approximates the posterior distribution over \(I_P(\mathsf{f}) \mid {\cal D}\) a Bayesian quadrature (BQ) method.

The integral posterior mean \(\mathbb{E}[ I_P(\mathsf{f}) \mid {\cal D}]\) or its approximation is used as a point estimate of the true integral \(I_P(f)\). The integral posterior variance, \(\mathrm{Var}[ I_P(\mathsf{f}) \mid {\cal D}]\), is typically used to summarise the spread of \(I_P(\mathsf{f}) \mid {\cal D}\). We shall always tacitly assume that the process \(\mathsf{f}\), the domain \(D\), and the probability measure \(P\) are such that the posterior distribution in ?? is well-defined and has finite mean and variance. If the function \(\varphi\) is affine, it suffices that \(m\) and \(k(\cdot, {\boldsymbol{x}})\) be integrable with respect to \(P\) for every \({\boldsymbol{x}}\in D\) and . For example, the integral is finite if the covariance function is stationary and \(P\) has a bounded density on a bounded \(D\) since in that case \(k({\boldsymbol{x}}, {\boldsymbol{x}}) = \Phi_k(\boldsymbol{0})\) for all \({\boldsymbol{x}}\in D\).

The selection of \(\mathsf{g}\) and \(\varphi\) is supposed to be informed about any prior information one may have about the true integrand \(f\). The prior information that is typically encoded includes different structural properties, such as stationarity, non-negativity or periodicity, and smoothness, described by the number of continuous derivatives that one expects \(f\) to possess. For example, when \(\mathsf{f}= \mathsf{g}\) and the covariance function of \(\mathsf{g}\) is a Matérn in 4 with smoothness \(\nu\), the samples of \(\mathsf{f}\) essentially have smoothness \(\nu\) in a Sobolev sense; see [31]–[33], as well as [34]. However, model misspecification, to be understood here in a broad sense that \(f\) “looks” different than the samples from \(\mathsf{f}\), is inevitable in anything more complicated than toy examples. Practical selection of the prior is a general challenge in Bayesian inference and not limited to Bayesian quadrature or inference with Gaussian processes. Prior selection, when simplified to mean the selection of a limited number of kernel parameters, such as \(\sigma\) and \(\ell\) of the square exponential kernel in 3 , is briefly reviewed in 4.3.

2.2.1 Conjugate Bayesian quadrature (affine \(\varphi\))↩︎

The prior \(\mathsf{f}\) is a GP if the function \(\varphi\) is affine. In this case Gaussian conjugate inference is possible and we accordingly call the resulting methods conjugate Bayesian quadrature methods. This class of methods is sometimes called standard Bayesian quadrature [35] or vanilla Bayesian quadrature [36]. Because an affine \(\varphi\) merely shifts and scales the mean and covariance of \(\mathsf{g}\), we can set \(\varphi\) to be the identity function without loss of generality. This gives the GP prior \[\mathsf{f}= \mathrm{id} \circ \mathsf{g}= \mathsf{g}\sim \mathrm{GP}(m, k).\] Due to the linearity of integration, the posterior for \(I_P(\mathsf{f})\) is a normal random variable that is obtained by simply integrating the posterior process \(\mathsf{f}\mid {\cal D}= \mathsf{g}\mid {\cal D}\): \[I_P(\mathsf{f}) \mid {\cal D}= I_P(\mathsf{f}\mid {\cal D}) = \int_D[ \mathsf{f}({\boldsymbol{x}}) \mid {\cal D}] \dif P({\boldsymbol{x}}) \sim \mathrm{N}( \mu_{{\cal D}}, \Sigma_{{\cal D}} ).\] The integral mean \[\label{eq:bq-mean} \begin{align} \mu_{{\cal D}} = \mathbb{E}[ I_P(\mathsf{f}) \mid {\cal D}] &= \int_Dm_{\cal D}({\boldsymbol{x}}) \dif P({\boldsymbol{x}}) \\ &= m_P+ {\boldsymbol{k}}_{PX}^{\intercal}{\boldsymbol{K}}_{XX}^{-1} ( {\boldsymbol{f}}_X- {\boldsymbol{m}}_X), \end{align}\tag{8}\] which is the integral of the posterior mean in 1 , is expressed in terms of the integral \(m_P= I_P(m)\) of the prior mean function \(m\), the integrand evaluations \({\boldsymbol{f}}_X= (f({\boldsymbol{x}}_1), \ldots, f({\boldsymbol{x}}_N)) \in \mathbb{R}^N\), and the evaluations \({\boldsymbol{k}}_{PX} = I_P({\boldsymbol{k}}_X) = ( k_P({\boldsymbol{x}}_1), \ldots k_P({\boldsymbol{x}}_N)) \in \mathbb{R}^N\) of the kernel mean embedding \[\label{eq:kernel-mean} k_P({\boldsymbol{x}}) = \int_Dk({\boldsymbol{x}}', {\boldsymbol{x}}) \dif P({\boldsymbol{x}}').\tag{9}\] The integral variance \[\begin{align} \label{eq:bq-var} \Sigma_{{\cal D}} = \mathrm{Var}[ I_P(\mathsf{f}) \mid {\cal D}] &= \int_D\int_Dk_{\cal D}({\boldsymbol{x}}, {\boldsymbol{x}}') \dif P({\boldsymbol{x}}) \dif P({\boldsymbol{x}}') \\ &= k_{PP} - {\boldsymbol{k}}_{PX}^{\intercal}{\boldsymbol{K}}_{XX}^{-1} {\boldsymbol{k}}_{PX}, \end{align}\tag{10}\] where \[\label{eq:kernel-var} k_{PP} = \int_D\int_Dk({\boldsymbol{x}}, {\boldsymbol{x}}') \dif P({\boldsymbol{x}}) \dif P({\boldsymbol{x}}') = \int_Dk_P({\boldsymbol{x}}) \dif P({\boldsymbol{x}})\tag{11}\] is the initial variance, is the integral of the posterior covariance in 2 . We encapsulate conjugate Bayesian quadrature and its posterior equations in the following definition.

Both the integral mean and variance can be expressed in terms of the weight vector \({\boldsymbol{w}}= {\boldsymbol{K}}_{XX}^{-1} {\boldsymbol{k}}_{PX}\) as \[\label{eq:conj-bq-mean-var} \mu_{\cal D}= m_P+ {\boldsymbol{w}}^{\intercal}( {\boldsymbol{f}}_X- {\boldsymbol{m}}_X) \quad \text{ and } \quad \Sigma_{\cal D}= k_{PP} - {\boldsymbol{w}}^{\intercal}{\boldsymbol{k}}_{PX}.\tag{12}\] Some of the weights \({\boldsymbol{w}}= (w_1, \ldots, w_N)\) can be negative [37], [38]. For a zero-mean prior the mean is simply \(\mu_{\cal D}= {\boldsymbol{w}}^{\intercal}{\boldsymbol{f}}_X\), which means that it is a quadrature rule, an approximation of the integral in the form of a weighted sum of integrand evaluations. Although it is customary to use the term cubature rule to describe such an approximation when \(d > 1\) [39], this distinction is not consistently made in the literature on Bayesian quadrature. For articles where Bayesian cubature is used, see [10], [40]–[42].

2.2.2 Non-conjugate Bayesian quadrature (non-affine \(\varphi\))↩︎

When \(\varphi\) is not affine, we speak of non-conjugate Bayesian quadrature.

Non-conjugate Bayesian quadrature is usually motivated by a desire to model the non-negativity of the integrand, a property which cannot be encoded by a GP model [8], [45]–[47]. Functions used to produce a non-negative \(\mathsf{f}\) include \[\varphi(x) = \alpha + x^2 \quad \text{ and } \quad \varphi(x) = \exp(x),\] where \(\alpha > 0\), and a general framework for encoding range constraints is described by [15]. See [48] and [49] for applications to engineering and neural networks. However, when \(\varphi\) is non-affine (as the two functions above are), the prior \(\mathsf{f}= \varphi\circ \mathsf{g}\) is usually not a Gaussian process and the posterior \(I_P(\mathsf{f}) \mid {\cal D}\) is not in general available in closed form, which precludes exact inference. Intractability of the posterior and its moments introduces a significant complication over conjugate Bayesian quadrature, forcing one to approximate the posterior. Gaussian approximations based on linearisation or moment-matching are typically used, but one can employ any posterior approximation available in the literature.

2.2.3 Bayesian probabilistic numerical integration↩︎

In the definition of Bayesian quadrature (2.1) the integrand is modelled as a potentially non-linearly transformed Gaussian process \(\mathsf{f}= \varphi\circ \mathsf{g}\). Nothing stops one from using any suitable stochastic process \(\mathsf{f}\). The following definition is adapted from [50].

All Bayesian quadrature methods are BPNI methods. Due to computational tractability issues, few non-Gaussian BPNI methods exist. Student’s \(t\)-processes, which differ little from Gaussian processes in practice, have been used in numerical integration by [5] and [51]. Integration methods based on Bayesian additive regression trees introduced by [50] appear to be the only non-Gaussian BPNI methods that substantially differ from Bayesian quadrature. There is much room for the development of new BPNI methods, particularly when integrand evaluations constitute the computational bottleneck. One can consider even more general and difficult integration problems in which the distribution \(P\) is unknown and accessible only via sampling or black-box evaluation of the density function [19], [52].

2.3 Connections to approximation theory↩︎

Given the well known equivalences of Gaussian process regression to optimal approximation in a reproducing kernel Hilbert space (RKHS) and to kernel interpolation [34], [53], [54], it is no surprise that these equivalences extend to conjugate Bayesian quadrature. We briefly describe the equivalences because of their historical value and because much of the theory for Bayesian quadrature, summarised in Sections 2.4 and 5, exploits them. See [9] and [11] for earlier reviews on the topic.

2.3.1 Optimal integration in a reproducing kernel Hilbert space↩︎

Let \(\mathcal{H}(k)\) denote the RKHS of \(k\). This is a certain Hilbert space consisting of real-valued functions defined on \(D\) in which the kernel has the following reproducing property: \(\langle g, k(\cdot, {\boldsymbol{x}}) \rangle_{\mathcal{H}(k)} = g({\boldsymbol{x}})\) for every \(g \in \mathcal{H}(k)\) and \({\boldsymbol{x}}\in D\). For example, the RKHS of a Matérn kernel 4 of order \(\nu\) is, up to an equivalent norm, the Sobolev space \(H^{\nu + d/2}(\mathbb{R}^d)\), a fact that is discussed in more detail in Section 5.1.1. We refer to [25], [55], and [56] for the theory of RKHSs. Section 7.1 contains a very concise review.

The worst-case integration error of any weights \({\boldsymbol{v}}\in \mathbb{R}^N\) and points \(X= \{{\boldsymbol{x}}_1, \ldots, {\boldsymbol{x}}_N\} \subset D\) in \(\mathcal{H}(k)\) is defined as \[\label{eq:wce-definition} e_k({\boldsymbol{v}}, X) = \sup_{ \norm[0]{g}_{\mathcal{H}(k)} \leq 1} \, \Big \lvert \int_Dg({\boldsymbol{x}}) \dif P({\boldsymbol{x}}) - \sum_{i=1}^N v_i g({\boldsymbol{x}}_i) \Big \rvert.\tag{13}\] It follows from the reproducing property that the worst-case error can be expressed in terms of the kernel.

See, for example, [57] or [58] for generalisations of Lemma 2.1. The worst-case error measures the largest possible integration error by a quadrature rule of the form \(\sum_{i=1}^N v_i g({\boldsymbol{x}}_i)\) for functions \(g \in \mathcal{H}(k)\) normalised in the RKHS norm. If the function \(f\) the quadrature rule is used to integrate is not an element of \(\mathcal{H}(k)\), the worst-case error is still well-defined, as is evident from ?? . However, in this case the worst-case error needs to bear no relation to the integration error \(\abs[0]{\int_Df({\boldsymbol{x}}) \dif P({\boldsymbol{x}}) - \sum_{i=1}^N v_i f({\boldsymbol{x}}_i)}\). As squared worst-case error is quadratic in \({\boldsymbol{v}}\), it is easy to see that, for fixed points, Equation ?? is minimised by the weights \({\boldsymbol{v}}= {\boldsymbol{w}}= {\boldsymbol{K}}_{XX}^{-1} {\boldsymbol{k}}_{PX}\), which are precisely the weights that define the integral mean and variance of conjugate Bayesian quadrature in 12 . Moreover, for these weights the worst-case error simplifies to \(e_k({\boldsymbol{w}}, X) = \Sigma_{\cal D}^{1/2} = \sqrt{\smash[b]{k_{PP} - {\boldsymbol{w}}^{\intercal}{\boldsymbol{k}}_{PX}}}\). The posterior mean of a conjugate Bayesian quadrature method is thus given by weights which minimise the worst-case error in the RKHS of the kernel \(k\) and the posterior variance equals the corresponding minimal squared worst-case error. These facts are summarised in the following proposition.

If the prior mean \(m\) is zero and the integrand \(f\) is an element of \(\mathcal{H}(k)\), we have the convenient error estimate \[\label{eq:rkhs-error-bound} \abs[0]{ I_P(f) - \mu_{\cal D}} = \langle f , k_P- \textstyle\sum_{i=1}^N w_i k(\cdot, {\boldsymbol{x}}_i) \rangle_{\mathcal{H}(k)} \leq \norm[0]{f}_{\mathcal{H}(k)} \Sigma_{\cal D}^{1/2},\tag{15}\] which follows from \(\langle f, k_P\rangle_{\mathcal{H}(k)} = I_P(f)\) (see the proof of Lemma 2.1, Section 7.1, or [59] Lem. 3.1) and the Cauchy–Schwarz inequality. This bound allows one to use the magnitude of \(\Sigma_{\cal D}\) to assess the quality of \(\mu_{\cal D}\) as a point estimator of the integral; its use is the reason behind the assumption \(f \in \mathcal{H}(k)\) in several theorems in Section 5. If \(f \notin \mathcal{H}(k)\), the bound 15 is not valid because the norm \(\norm[0]{f}_{\mathcal{H}(k)}\) does not exist.

The construction of general worst-case optimal quadrature rules in RKHSs goes back at least to the work of Larkin in the 1970s [60]–[62], while the related average-case framework was developed by [63], [64].⁴ For early analyses, mostly in RKHSs of analytic functions, of these quadratures we refer to [68]–[70]. See also [71], [58], [72]–[74]. A different tradition may be found in the quasi-Monte Carlo literature [3], where spaces of periodic functions are often considered. The optimal weights in these spaces can be occasionally computed in closed form [75]. Certain quasi-Monte carlo rules can be interpreted as Bayesian quadratures [10], [41], [76].

Whenever making use of the RKHS framework one should keep in mind that the samples from a Gaussian process are not in the RKHS of its kernel. Specifically, under fairly non-restrictive assumptions one can show that if \(\mathsf{f}\sim \mathrm{GP}(m, k)\), then \(\mathbb{P}[ \mathsf{f}\in \mathcal{H}(k)] = 0\); see [77] and [34]. That is, the RKHS is too small to contain the samples. One cannot therefore simply substitute the process \(\mathsf{f}\) for the RKHS element \(f\) in results such as the error estimate 15 .

2.3.2 Kernel interpolation↩︎

A different tradition may be found in the literature on kernel interpolation or, when the kernel is isotropic, radial basis function interpolation [78]. Given a positive-definite kernel \(k\) and evaluations of \(f\) at a collection \(X= \{{\boldsymbol{x}}_1, \ldots, {\boldsymbol{x}}_N\}\) of pairwise distinct points, the kernel interpolant \(s_{f, X}\) is the unique function in the span of the kernel translates \(\{k(\cdot, {\boldsymbol{x}}_i)\}_{i=1}^N\) that interpolates \(f\) at \(X\), in the sense that \(s_{f, X}({\boldsymbol{x}}_i) = f({\boldsymbol{x}}_i)\) for every \(i = 1, \ldots, N\). From the requirement that the interpolant be in the span of the translates it is easy to deduce that \(s_{f, X}({\boldsymbol{x}}) = {\boldsymbol{k}}_X({\boldsymbol{x}})^{\intercal}{\boldsymbol{K}}_{XX}^{-1} {\boldsymbol{f}}_X\). We see that the kernel interpolant coincides with the GP conditional mean 1 under a zero-mean prior when \({\boldsymbol{y}}= {\boldsymbol{f}}_X\). The kernel interpolant may be alternatively defined as the minimum-norm interpolant in the RKHS of \(k\): \[\label{eq:minimum-norm-interpolation} s_{f, X} = \mathop{\mathrm{arg\,min}}_{s \in \mathcal{H}(k)} \Set[]{ \norm[0]{s}_{\mathcal{H}(k)} }{ s({\boldsymbol{x}}_i) = f({\boldsymbol{x}}_i) \text{ for } i = 1, \ldots, N}.\tag{16}\] See, for example, Theorem 3.5 in [34]. Due to the minimum-norm property kernel interpolants are often called spline interpolants ([71], Sec. 5.7; [79]). As the kernel interpolant \(s_{f, X}\) depends on the value of \(f\) at a finite point set, it is well-defined whether \(f\) is in \(\mathcal{H}(k)\) or not. However, the useful inequality , a consequence of the minimum-norm interpolation property, holds only if \(f\in \mathcal{H}(k)\).

Integrating the kernel interpolant yields the integral approximation \(I_P(f) \approx I_P(s_{f, X}) = {\boldsymbol{k}}_{PX}^{\intercal}{\boldsymbol{K}}_{XX}^{-1} {\boldsymbol{f}}_X\), which is equal to the conjugate BQ integral mean 8 in the mean-zero case. [80] appears to have been the first to suggest this approach. More recent work on integrated kernel interpolants includes [81]–[91]. [92] discusses the interpretation of spline-based quadrature as BQ. The connection to kernel interpolation is important because many estimates for the integration error \(\lvert I_P(f) - \mu_{\cal D}\rvert\) of Bayesian quadrature are derived from \(L^p\)-error estimates for kernel interpolation [11], [22].

2.3.3 Kernel quadrature↩︎

Quadrature rules with the weights 14 that minimise the worst-case error in an RKHS are often called kernel quadratures [93]–[96]. However, this nomenclature is more commonly (and somewhat ambiguously) applied to any quadrature rule that seeks to achieve small worst-case error in an RKHS and does this by making use of the kernel or objects derived from it [97], [98].⁵ For example, [99], [102]–[109] use Mercer expansions of \(k\) to construct quadrature rules. See also [110], [111].

2.3.4 Bayesian quadrature in disguise↩︎

Many classical quadrature rules constructed with non-probabilistic arguments can be cast as conjugate Bayesian quadrature if the covariance is selected appropriately [40], [105], [112], [113] or recovered by considering the length-scale limit \(\ell \to \infty\) [5], [92], [114], [115]. See also [35] and [13].

It is a well-known and established fact that the integral mean of conjugate Bayesian quadrature is equivalent to the classical trapezoidal rule if \(k(x, x') = \min\{x, x'\}\) is the covariance function of the Brownian motion [4], [116]. Recall that the trapezoidal rule is the integral approximation \[\int_0^1 f(x) \dif x \approx \sum_{i=1}^N \frac{f(x_i) + f(x_{i-1})}{2} (x_i - x_{i-1}),\] where \(0 = x_0 < x_1 < \cdots < x_{N-1} < x_N = 1\). The geometric interpretation is that the integral is approximated with \(N\) trapezoids with base lengths \(x_i - x_{i-1}\). [117] is an accessible source for computations that constitute the proof of the following proposition.

We refer to [118] for more explicit computations of this kind for the Brownian motion covariance.

Some other constructions have a hint of artificiality, relying as they do on finite-dimensional covariance functions [113] or limiting arguments that can be used to “force” Bayesian quadrature to take any form desired [40].

2.4 Fundamental theory for Bayesian quadrature↩︎

With the results and tools introduced above it is possible to establish some fundamental theory for conjugate Bayesian quadrature that holds regardless of the properties the covariance kernel has. The following proposition consists of the intuitively obvious fact that new data decrease the integral variance.

Note that both \(\Sigma_{\cal D}\) and \(\Sigma_{{\cal D}'}\) in Proposition [prop:monotone-var] are based on the same covariance function. In practice it can happen that the variance increases as more data are obtained if the covariance parameters are being learned from the data. For example, new data could indicate that the integrand is much rougher than previously thought. We discuss parameter estimation in Section 4.3. The next proposition describes how the covariance function affects the variance via its RKHS.

Inclusions between RKHSs of standard covariance functions are well-understood [34]. For example, if \(k_\nu\) is a Matérn covariance in 4 with regularity \(\nu\) and \(k_\mathrm{se}\) the square exponential in 3 , then \[\mathcal{H}(k_{\nu_1}) \subseteq \mathcal{H}(k_{\nu_2}) \subsetneq \mathcal{H}(k_\mathrm{se})\] if \(\nu_1 \leq \nu_2\) and the domain \(D\) has non-empty interior. The kernels in this chain of inclusions are allowed to have different length-scales. The final proposition of this section says that conjugate Bayesian quadrature integrates the kernel translates \(k(\cdot, {\boldsymbol{x}}_i)\) exactly.

2.5 Generalisations↩︎

We have elected to work in the setting described above to retain a degree of concreteness and to avoid introduction of overbearing notation. However, a number of conceptually straightforward (though not necessarily so in practice) generalisations are possible and appear often in the literature. These include the following:

• Domain. There is nothing in Bayesian quadrature or Gaussian processes that ties them to subsets of \(\mathbb{R}^d\). In general, one can let \(D\) be any set and \(P\) any measure on \(D\). This generality of Bayesian quadrature has been recognised since the very beginning [5], [61]. A recurring example in the literature is the computation of global illumination integrals on the unit sphere \(\mathbb{S}^2 = \Set{ {\boldsymbol{x}}\in \mathbb{R}^3}{ \norm[0]{{\boldsymbol{x}}} = 1}\), to which Bayesian quadrature has been applied in [9], [14], [16], [119]–[121]. Other domains may be found in [122] and [123]. Integration on a sphere has also been considered in the literature on radial basis function integration [81], [86], [87] whose connection to Bayesian quadrature we described in Section 2.3.2.

• Information. It is straightforward to incorporate any linear information (i.e., observations that have been obtained via application of a linear functional) in Gaussian process regression [124] and, consequently, Bayesian quadrature. We refer to [58] for a rigorous treatment of the topic in the context of approximation theory, the connection to Gaussian processes of which we reviewed in Section 2.3. Derivative observations, which occur frequently in Gaussian process literature [125], have been used in Bayesian quadrature by [126] and [127].

• Noise. Incorporating observations that are perturbed by Gaussian noise, such as when \(f\) is the output of a stochastic simulation, to Bayesian quadrature is straightforward. In the typical Gaussian conjugate setting this simply amounts to incrementing the diagonal of the covariance matrix with the noise variances (see the equations for the regularised mean and variance in 4.4); this is supported by most open source implementations of Bayesian quadrature. Even though the setting of noisy observations is dominant in the Gaussian process regression literature, Bayesian quadrature with noisy observations appears to have received little explicit attention. The only substantial analysis the noisy case that we are aware of is by [128]. [129] also incorporate noise in their method.

• Prior. In our definition of Bayesian quadrature the integrand needs to be modelled as a potentially non-linearly transformed Gaussian process \(\mathsf{f}= \varphi\circ \mathsf{g}\). As discussed in Section 2.2.3, one is in principle free to use any suitable stochastic process \(\mathsf{f}\) to model the integrand. Due to computational tractability issues, priors other than transformed Gaussian processes are rarely used. Student’s \(t\)-processes, which differ little from Gaussian processes in practice, have been used in numerical integration by [5] and [51]. What appears to be the most interesting example of non-Gaussian priors in numerical integration is due to [50], who use Bayesian additive regression trees.

3 Taxonomy of Bayesian quadrature↩︎

Definition 2.1 constructed Bayesian quadrature as a method for computing or approximating the distribution \(I_P(\mathsf{f})\mid{\cal D}\). This abstract formulation admits several concrete realizations, each corresponding to a particular Bayesian quadrature method or algorithm. In recent years, the literature has proposed a wide variety of approaches consistent with Definition 2.1, differing mainly in their methodological choices and practical implementation. This chapter offers a concise, high-level overview of these alternatives. It may be read as a taxonomy of Bayesian quadrature methods or used as a reference for terminology. The taxonomy is illustrated visually in Figure 8, and Table ¿tbl:tab:references? classifies selected algorithms from the literature.

To motivate the taxonomy, we inspect a high-level pseudocode for Bayesian quadrature shown in Algorithm 4. The pseudocode reveals three core components: a model for \(\mathsf{f}\) specified by \(\mathsf{g}\) and \(\varphi\) (see 7 ); a sampling procedure (line 1) for selecting the nodes \(X\) that constitute the dataset; and an inference procedure (line 4) that computes the return object \(res\) representing the distribution \(I_P(\mathsf{f}) \mid {\cal D}\).⁶

This naturally leads to three main axes along which Bayesian quadrature methods can be organized:

• a model axis, which characterizes the Gaussian process prior \(\mathsf{g}\) and the transformation \(\varphi\) (3.1);

• an inference axis, which specifies how the result \(res\) is computed (3.2);

• a sampling axis, which describes how the nodes \(X\) are selected (3.3).

Along each axis, we introduce multi-label subcategories intended to reflect the current landscape of methods and algorithms. Here, “multi-label” indicates that categories are not necessarily mutually exclusive (when they are, this will be stated explicitly); consequently, a method may be assigned multiple labels per axis. An illustrative example is provided in 3.5.

Not all Bayesian quadrature methods fit neatly into the pseudocode outlined above. This primarily concerns approaches that actively manage data acquisition and computational budgets (3.3). Inevitably, any taxonomy involves a degree of subjectivity; ours is guided by practical considerations motivated by applications in applied science.

We occasionally use the big-\({\mathcal{O}}\) notation when referring to computational complexity or memory requirements. Accessible implementations of Bayesian quadrature algorithms are unfortunately still scarce and not on par with the richness of the literature. Some are available in the open source libraries EmuKit (for Python; [130], [131]), GAIL (for MATLAB; [132]), and ProbNum (for Python; [133]).

3.1 Model↩︎

The first axis of the taxonomy is concerned with the prior model \(\mathsf{f}= \varphi\circ \mathsf{g}\), in particular the choice of the function \(\varphi\). We differentiate between two mutually distinct classes: “Conjugate” and “non-conjugate” that in essence categorize the posterior \(I_{P}(\mathsf{f}) \mid {\cal D}\) as Gaussian or non-Gaussian.

• Conjugate: Conjugate models are arguably the most well-known and the most straightforward Bayesian quadrature models. They are described in Section 2.2.1, where \(\varphi\) is affine and the Gaussian process \(\mathsf{g}\) implies a Gaussian process \(\mathsf{f}\) and a Gaussian posterior \(I_{P}(\mathsf{f}) \mid {\cal D}\). Conjugate Bayesian quadrature models inherit the conveniences of GP regression models, and \(I_{P}(\mathsf{f}) \mid {\cal D}\) can be computed exactly under certain well-understood circumstances (Section 3.2). For these reasons, conjugate Bayesian quadrature is sometimes also called standard Bayesian quadrature [35], [134] or vanilla Bayesian quadrature [36], [130], [135]. As \(\mathsf{f}\) and \(I_{P}(\mathsf{f}) \mid {\cal D}\) are confined to Gaussians, conjugate models are limited in their expressiveness.

• Non-conjugate: Non-conjugate models are often motivated by the desire to encode physically meaningful, non-negative integrands \(f\) and integral values \(I_{P}(\mathsf{f}) \mid {\cal D}\) which Gaussians cannot provide. Steering away from Gaussianity is possible with an appropriate function \(\varphi\). Non-conjugate models are described in Section 2.2.2 where \(\varphi\) is not affine, \(\mathsf{f}\) is a random process other than a Gaussian process and the implied posterior \(I_{P}(\mathsf{f}) \mid {\cal D}\) is non-Gaussian. For example, for \(\varphi(x)=\alpha + x^2\) as in [47], \(\mathsf{f}\) is a process whose marginals are noncentral \(\chi^2\)-distributions of degree 1. Inference, however, is usually less straightforward and may require approximations (3.2). Due to the presence of \(\varphi\), non-conjugate Bayesian quadrature is sometimes also referred to as warped⁷ Bayesian quadrature [47], [130]. While the space of possible \(\varphi\) is rich in principle, only a limited number have been explored in the context of BQ.

Further, a kernel \(k\) and mean function \(m\) identify the Gaussian process \(\mathsf{g}\). As mentioned in 2.1, common choices for \(k\) are the square exponential kernel and kernels of the Matérn family, either isotropic or of product form. Any model may additionally be characterized by hyper-parameters (such as kernel parameters) that are either given or set by an appropriate method (4.3).

3.2 Inference↩︎

Once the model \(\mathsf{f}\) has been defined and a dataset \({\cal D}\) has been obtained (node sampling forms the third axis of the taxonomy and is the content of 3.3), one can attempt to compute the solution \(I_P(\mathsf{f}) \mid {\cal D}\). The second axis of the taxonomy—inference—distinguishes the quality of the solution as either “exact” or “approximate”. A third label “scalable” emphasizes computational complexity of the inference algorithm with respect to the dataset size, \(N\).

• Exact: The inference method is called exact if \(I_P(\mathsf{f}) \mid {\cal D}\) can be computed without approximation error (\(res\) in Algorithm 4 is an exact representation). This is possible if \(I_P(\mathsf{f}) \mid {\cal D}\) has closed form and admits numerically exact computations. Exact inference can for example be realized under the conjugate BQ model if the measure \(P\) and kernel \(k\) pair yields closed form solutions of Eqs. 8 and 10 ; see Section 4.1. Since the expressions have closed form, exact inference in Bayesian quadrature is also called analytic. 4.4 discusses some common practical considerations concerning numerical precision.

• Approximate: The inference method is called approximate if the returned distribution (\(res\) in Algorithm 4) differs from \(I_P(\mathsf{f}) \mid {\cal D}\) due to methodological choices.⁸ Such approximations to \(I_P(\mathsf{f}) \mid {\cal D}\) are necessary if (i) the model \(\mathsf{f}\mid{\cal D}\) does not admit analytic integration, or (ii) if evaluating the analytic expressions to numerical precision is prohibitively expensive. In practice, inference methods for non-conjugate BQ models often entail approximations due to the former reason because the process \(\mathsf{f}\) can be arbitrarily complicated. Some recent approaches use Gaussian approximations to \(\mathsf{f}\) via linearisation or moment matching which are integrated instead of \(\mathsf{f}\). Conjugate Bayesian quadrature, where \(\mathsf{f}\) is a Gaussian process, may also require approximations when the kernel \(k\) and integration measure \(P\) pair has no closed form kernel mean embedding (Section 4.1). Approximations of kernel embeddings have been scarcely studied in the context of BQ, but in principle, any means of approximation is admissible. Another reason arises when the closed form solution is too expensive to compute for large dataset sizes \(N\) because the computational complexity scales unfavorably. Ways to remedy this are discussed next.

Algorithm 5 shows a pseudocode of the exact inference method for a conjugate Bayesian quadrature model. From the pseudocode it is apparent that exact inference for the conjugate model is feasible if the kernel mean embedding \(k_{P}\), the initial variance \(k_{PP}\) and the prior mean integral \(m_{P}\) are analytic. Algorithm 5 first computes \(k_{PP}\) and \(m_{P}\) (lines 2–3). Then, it evaluates the kernel \(k\), mean function \(m\) and kernel mean embedding \(k_{P}\) at the given nodes (lines 4–5) and constructs the kernel Gram matrix \({\boldsymbol{K}}_{XX}\) (line 6). The bulk of the computation is performed in line 8 where the algorithm solves a linear system of size \(N\) in order to retrieve the weights \({\boldsymbol{w}}= {\boldsymbol{K}}_{XX}^{-1} {\boldsymbol{k}}_{PX}\). Lines 9–10 then perform simple vector multiplications and summations that yield the final solution.

Algorithm 5 shows arguably the most well-known BQ setting which inherits its quadratic memory \({\mathcal{O}}(N^2)\) (line 6) and cubic compute cost \({\mathcal{O}}(N^3)\) (line 8) from standard Gaussian process regression. Thus, the final label of the inference axis describes the scalability of the BQ inference method with respect to \(N\) compared to the standard cubic cost. This label is not mutually exclusive to the previous two; in fact, both an “exact” and “approximate” inference method can be either scalable or not.

• Scalable: The Bayesian quadrature inference method is called scalable if its computational complexity is less than cubic in the number of datapoints \(N\) and thus scales favorably compared to standard Gaussian process regression. Interestingly, so far scalable methods have been proposed primarily for exact inference under the conjugate BQ model. For instance, [94] exploit structure in the posterior mean equation for certain kernels, measures and a custom node sampling scheme; their algorithm scales sub-cubically. [10] use fast Fourier transform to obtain computational cost \({\mathcal{O}}(N \log N)\). Scalable methods for approximate inference under the conjugate BQ model have not been explored systematically.⁹ [142] use a scalable approximate BQ method based on structured kernel interpolation [143] to estimate policy gradients and [144] use sparse GPs within variational BQ. We are not aware of scalable methods for non-conjugate BQ models. Scalable Bayesian quadrature is sometimes also referred to as fast Bayesian quadrature [10], [41]. This is not to be confused with the use of the term fast BQ as in [47] or [135] where the term highlights sample efficiency or end-to-end wall-clock time performance of the BQ method and data retrieval which our taxonomy does not address. Here we also mention a recent quantum algorithm for Bayesian quadrature by [145].

3.3 Sampling (node selection)↩︎

The first two axes described the BQ model and inference method when a dataset \({\cal D}\) is provided. The third and last axis of the taxonomy focuses on the dataset \({\cal D}\), particularly emphasizing the selection of the nodes \(X\) (as the corresponding function values \({\boldsymbol{f}}_{X}\) can be retrieved by evaluating \(f\) at \(X\)). We start by observing that traditionally, to solve an intractable integral, practitioners may refer to classic quadrature rules where nodes are predefined and cannot be chosen freely, or Monte Carlo estimators, where nodes must be random. Bayesian quadrature is positioned in between these methods as it generally can use any nodes \({\boldsymbol{x}}_1,\dots,{\boldsymbol{x}}_N\) as long they are pairwise distinct. This fortunate property of BQ leads to a rich variety of practical sampling options wherefore this section is divided into several subsections.

3.3.1 Mode of sampling↩︎

We begin with three high-level categories that are fundamental to the mode in which the nodes are obtained: “Deterministic”, “random”, and “fixed”. The former two are mutually distinct, while the third is simply undetermined.

• Deterministic: Deterministic node selection schemes are those schemes that ideally do not use random or pseudo-random numbers, except potentially for an initialization step such as random selection of the generating vector for lattice nodes [10], [75]. In practice, deterministic schemes often comprise a degree of randomness, for example when the nodes are selected by maximizing an acquisition function and the optimizer makes use of randomness.¹⁰

• Random: Random sampling schemes are those that explicitly use random or pseudo-random numbers to generate the nodes. In contrast to “deterministic” sampling as described above, “random” sampling has a degree of randomness intrinsic to its design and even under ideal practical conditions. An early use of random node sampling for BQ was proposed as part of the Bayesian Monte Carlo method by [7] who sampled from the integration measure \(P\). In the meantime, also other distributions have been used [93], [97].

• Fixed / unknown: Fixed or unknown sampling describes a set of nodes whose generating process is either unknown or beyond the control of the Bayesian quadrature practitioner. A fixed set of nodes could for example be provided by a (black-box) service.

The above three categories describe the degree of control that the practitioner has (or rather needs to have) over the nodes, strictly within the context of Bayesian quadrature: When the nodes are fixed, one has no control; when they are random, one has only distributional control; and when they are deterministic, one has full control over the nodes (at least in principle).¹¹ The remaining categories provide more granularity to node selection.

3.3.2 Contextual sampling↩︎

We now inspect line 1 of Algorithm 4 more closely and observe that the model \(\mathsf{f}\) is an optional argument to the sample method which means that some sampling procedures may depend on their context. Thus the first granular label of the sampling axis is called “model-dependent”.

• Model-dependent: It is easy to fathom sampling schemes that are independent of the model \(\mathsf{f}\), such as grids or draws from the integration measure \(P\). If the sampling method does not request the model \(\mathsf{f}\) (the sample method in line 1 of Algorithm 4 does not require \(\mathsf{f}\) as argument) we call it “model-independent” and if it does, “model-dependent”. As of now, this definition is purely algorithmic in the sense that the nodes \(X\) could be created either before or after \(\mathsf{f}\) is chosen. As we will see soon, the distinction is somewhat more intricate than it may initially appear since it disentangles (more so than other labels) the method from the algorithm. This is best explained with an example: It is a well-known fact [13] that for \(\mathsf{f}\) a Brownian motion and \(P\) the uniform measure on the interval \([0, 1]\), equidistant nodes are the maximally-informative design in the sense that \(\mathrm{Var}[ I_P(\mathsf{f}) \mid {\cal D}]\) is minimized. The creation of equidistant nodes on \([0, 1]\) of course does not require \(\mathsf{f}\) explicitly. However, since the practitioner promises maximally informative nodes, knowledge about \(\mathsf{f}\) (i.e., it being a Brownian motion) is required before the nodes are created. In other words, if the model \(\mathsf{f}\) changes, the sampling method will need to change as well. If the practitioner were to use the equidistant nodes no matter \(\mathsf{f}\) (no “link” between \(\mathsf{f}\) and the sampling method is claimed), the BQ method would be model-independent. To conclude the example, maximally informative nodes for processes \(\mathsf{f}\) other than a Brownian motion are not generally analytic and their computation may require \(\mathsf{f}\) explicitly in which case it being model-dependent is apparent. Model-dependent sampling often implies a notion of informativity of the nodes under the model. A recent example of model-dependent node selection is given in [135] who construct a sampling distribution that depends on \(\mathsf{f}\) as well as \(P\).¹² Further, active sampling (Section 3.3.4) is model-dependent by definition as the acquisition function depends on \(\mathsf{f}\). It is highly unclear when model-dependent or independent sampling for Bayesian quadrature should be preferred and it may be contingent on the application at hand. Although model-independent node selection has no knowledge about \(\mathsf{f}\) it often performs well, is generally fast as nodes can be pre-computed, and lends itself well to theoretical examination. Model-dependent node selection has the benefit of being tailored to the integration problem at hand.

Model-dependent sampling, as mentioned, usually implies that the nodes are chosen due to their informativity under the model \(\mathsf{f}\). Another type of dependence may occur if the sampler is interlaced with the model additionally or exclusively via the inference method. Such nodes may or may not be chosen due to their informative content, but most likely they are chosen for making the inference method “work”. An example are inference methods that are only scalable (have computational cost less than cubic in \(N\)) if the nodes follow a specific sampling strategy. The inference method can even be connected to the sampler in such a way that choosing any other sampler would not lead to an executable algorithm. To highlight the association we will call these sampling schemes “inference-dependent”. This is a good point to reiterate that the taxonomy labels should be attached to specific instances of BQ methods that comprise a model, an inference and a sampling method. Hence labels need to be interpreted with respect to the BQ method as a whole and in the context of the practitioner’s aim. A context-free taxonomy with merit is almost impossible as BQ methods only allow to plug-and-play under restrictions.¹³

• Inference-dependent: Some BQ methods require specific node sampling schemes that are motivated by the unique properties of their inference method. We call those samplers “inference-dependent”. As mentioned, examples are the sampling schemes of the exact, scalable BQ methods in [94] and [10]. The concept is similar in spirit (albeit not as strict) to how nodes in classical Gaussian quadrature rules are derived.

3.3.3 Sequential sampling↩︎

So far, we have worked with Algorithm 4 where the nodes \(X\) are sampled (line 1) as one large batch of size \(N\). We are now ready to study “sequential” sampling where nodes are acquired repeatedly in a loop in smaller sized batches (and often even only one at a time). Algorithm 6 shows the general layout of a sequential method.

We observe that Algorithm 6 is essentially a loop (lines 3–7). Before each step (one step equals executing lines 4–6) a stopping_condition checks if the iteration should be halted (line 3). The stopping condition may use information of the dataset \({\cal D}\) containing the data collected thus far, the model \(\mathsf{f}\), and the measure \(P\). During each step, a sampler produces a batch of nodes \(X_\mathrm{batch}\) of size \(n\) (line 4) at which the integrand is evaluated (line 5). Finally, the dataset \({\cal D}\) is augmented with the newly acquired nodes and observations \({\boldsymbol{f}}_{X_\mathrm{batch}}\) (line 6). Once the stopping condition is fulfilled, the loop terminates (the while-statement in line 3 is false and the loop breaks) and the solution \(res\) is computed and returned (lines 8–9).

Algorithm 6 not only has different structure than Algorithm 4 but also requires different inputs, most prominently a stopping condition, which alludes to one of the general motivations of sequential BQ methods: the sampling (and hence the evaluation of \(f\)) is paced in order to test for some user-defined stopping criterion. For example, the practitioner may want to terminate the BQ method when the variance \(\mathrm{Var}[ I_P(\mathsf{f}) \mid {\cal D}]\) falls below a tolerance value. The BQ method samples as many (small) batches as needed until the stopping condition is fulfilled and it is unknown a priori how many nodes the method eventually collects in total (this is elaborated further under dynamic budget in Section 3.4). A second, orthogonal motivation for sequential node acquisition emerges if the sampling method in line 4 is model-dependent and changes with the collected dataset up to that point (since \(\mathsf{f}\mid {\cal D}\) changes). A sequential approach is then able to gather informative nodes at each step of the loop depending on the current evidence and prior model. An example of the latter are some versions of ‘active’ sampling introduced in the next subsection.

• Sequential: Sequential sampling produces nodes in an iterative fashion in small batches of size \(n\) as sketched in Algorithm 6. Sequential sampling requires a stopping condition that indicates when the data acquisition should terminate. Often, nodes are sampled one at a time (\(n=1\)); if \(n>1\) then sequential sampling is also called batched. [148] were the first to focus on sequential BQ methods.

3.3.4 Active sampling↩︎

The last label of the sampling axis of the taxonomy is “active” sampling for BQ. The term “active” is used in various machine learning fields and generally describes that the model can query useful or informative data. The idea is that the model may perform better if it can actively select data rather than just passively accept data that is handed to it. Choosing the data thus constitutes a decision problem that the model needs to solve. In the context of BQ this means that the model \(\mathsf{f}\) influences the sampling method and “requests” \(f\) to be evaluated at certain nodes. For probabilistic models (the BQ model is a probabilistic model), an explicit formulation of this decision problem is possible via Bayesian decision theory. There, node selection is cast as the solution of the following optimization problem (we will use this formulation throughout for active sampling): \[\label{eq:active-sampling} \begin{align} {\boldsymbol{x}}_1,\dots, {\boldsymbol{x}}_N &=\mathop{\mathrm{arg\,max}}_{\tilde{{\boldsymbol{x}}}_1,\dots, \tilde{{\boldsymbol{x}}}_N} \, \mathbb{E}_{\mathsf{f}}[u(\tilde{{\boldsymbol{x}}}_1,\dots, \tilde{{\boldsymbol{x}}}_N\mid\mathsf{f}, P)]. \end{align}\tag{17}\] In the above formula, the nodes \({\boldsymbol{x}}_1,\dots, {\boldsymbol{x}}_N\) are the maximizer of an expected utility called an acquisition function \[a(\tilde{{\boldsymbol{x}}}_1,\dots, \tilde{{\boldsymbol{x}}}_N\mid P) \mathrel{\vcenter{:}}= \mathbb{E}_{\mathsf{f}}[u(\tilde{{\boldsymbol{x}}}_1,\dots, \tilde{{\boldsymbol{x}}}_N\mid\mathsf{f}, P)]\] that, together with the \(\mathop{\mathrm{arg\,max}}\) operator, defines the node selection policy and hence the sampling method of the active BQ method. The function \(u(\cdot\mid\mathsf{f}, P)\) is a utility function that quantifies the usefulness of evaluating the integrand \(f\) at certain nodes under the model \(\mathsf{f}\).¹⁴

While the utility function is arbitrary in the sense that it describes the practitioner’s conviction of usefulness, the active BQ method internally assumes the resulting sampling policy to provide the best course of action. Once again, this is best explained through an example. Let \(\tilde{{\cal D}} = \{ (\tilde{{\boldsymbol{x}}}_1, \mathsf{f}(\tilde{{\boldsymbol{x}}}_1)), \ldots (\tilde{{\boldsymbol{x}}}_N, \mathsf{f}(\tilde{{\boldsymbol{x}}}_N)) \}\) denote a hypothetical dataset, meaning that the targets are the hypothetical values \(\mathsf{f}(\tilde{{\boldsymbol{x}}}_i)\) under the model [instead of the integrand values \(f(\tilde{{\boldsymbol{x}}})\)], and let \(\mathsf{f}\mid\tilde{{\cal D}}\) be the model conditioned on the hypothetical dataset. If we choose \(u(\cdot\mid\mathsf{f}, P) = -(I_{P}(\mathsf{f}\mid \tilde{{\cal D}}) - \mu_{\tilde{{\cal D}}})^2\) to be the negative [Eq. 17 is a maximization problem] square error between the integral \(I_{P}(\mathsf{f}\mid \tilde{{\cal D}})\) of the model and its expectation \(\mu_{\tilde{{\cal D}}}\) (the square error is a natural and decent choice of utility), the resulting acquisition function \(a(\tilde{{\boldsymbol{x}}}_1,\dots, \tilde{{\boldsymbol{x}}}_N\mid P)=-\mathrm{Var}[ I_P(\mathsf{f}) \mid \tilde{{\cal D}}]\) equals the negative variance of the integral over the conditioned model and its maximization yields maximally informative nodes [13]. Likewise, other utility functions can be linked to well-known BQ acquisition functions, some of which will be discussed further in Section 4.2.3.

Eq. 17 assumes that the practitioner aims to sample a set of \(N\) nodes at once. Often, active node sampling is deployed in a sequential fashion due to the mentioned potential benefits of sequential node acquisition of model-dependent sampling schemes. In such cases, the utility function takes into account sampling decisions under the model in future steps of the iteration when making the sampling decision in the current step [150]. In practice, for tractability reasons often only an isolated decision is considered, meaning that the utility function encodes the usefulness of sampling the current step (or even just a single next node) as if it were the last step/node that the model can ever request. Such an approach is called one-step-lookahead or (maximally) myopic sequential active sampling and is usually the default over more elaborate schemes due to its practicality and tractability [150].¹⁵ Due to its importance, Algorithm 7 shows this approach.

Let us inspect the return object (line 5) of Algorithm 7 which is a single node \({\boldsymbol{x}}\) since the active sampler is maximally myopic with batch size \(n=1\).¹⁶ The sampled node is the return object of a maximum acquisition policy (line 4). The policy, which is highlighted in a separate function on lines 8–11, takes an acquisition function as input that is defined like a lambda function on lines 2 and 3 inside the sampling method based on \(\mathsf{f}\), \(P\) and \({\cal D}\). The dependence of \(u(\cdot \mid \mathsf{f}, P, {\cal D})\) on the current dataset \({\cal D}\) in Algorithm 7 is due to the sequential nature of the sampler and was not present in Eq. (17 ). The maximizer of the optimization problem on line 9 can either be implemented directly if an analytically solution is available, or be found by off-the-shelf optimizers, especially if gradients of \(\mathsf{f}\) are accessible.¹⁷

We are finally ready to define the “active” label.

• Active: Actively sampled nodes are the solution to an optimization problem as defined in Eq. 17 . They require the choice of a utility function \(u\) by the practitioner. Often active samplers are deployed in combination with an iterative, sequential sampling scheme. Similar concepts exist in experimental design [154] and Bayesian optimization [150], [155] where the utility functions are tailored to the respective quantity of interest. Section 4.2.3 discusses utility functions specific to Bayesian quadrature. If the resulting acquisition function depends on some previously collected dataset \({\cal D}\), we further distinguish implicitly and explicitly active sampling which describes the type of dependence.

  • Implicitly active: Active sampling is implicit if the acquisition function depends on the observations \({\boldsymbol{f}}_X\) only via a data-dependent model for \(f\). This is for example the case when the kernel-parameters of a conjugate BQ model are estimated with empirical Bayes after each step, while performing sequential active sampling (Section 4.3 will discuss kernel-parameter estimation). Implicitly active sampling is sometimes also called semi-active, non-adaptive, open-loop [13] or simply active [36].

  • Explicitly active: Active sampling is explicit if the acquisition function depends on \({\boldsymbol{f}}_X\) directly, beyond a data-dependent model. An example is uncertainty sampling in the WSABI method [47]. In contrast to implicitly active sampling, explicitly active sampling schemes cannot be precomputed if data-dependency of the model is dropped. Explicitly active sampling thus naturally combines with sequential sampling. Explicitly active sampling is also referred to as adaptive¹⁸ or closed-loop [13] sampling.

The taxonomy has illustrated that Bayesian quadrature should be understood as an umbrella term (under Definition 2.1) and that Bayesian quadrature methods differ significantly in model assumption, inference method, sampling procedure and even algorithmic requirements and structure. They indeed range from methods that accept fixed datasets to agent-like methods that interact with the data collection mechanism. Figure 8 illustrates the taxonomy and Table ¿tbl:tab:references? categorizes some existing BQ methods.

3.4 Other terms↩︎

This section lists terms that are occasionally used to describe some characteristics of Bayesian quadrature algorithms but are not part of our taxonomy. We list them for completeness.

• Structured: Bayesian quadrature is sometimes called structured¹⁹ when the model \(\mathsf{f}\) encodes a known structural property, such as an invariance or range constraint, of the integrand \(f\) via the transformation \(\varphi\) or otherwise. This is to highlight that the method is designed to address certain prior information that may be available about \(f\). For example, [134] encode input invariances of the integrand \(f\) in a conjugate BQ model. [15] and [47] use exponential and square transforms \(\varphi\), respectively, to encode non-negativity in a non-conjugate model (see Section 2.2.2).

• Fixed & dynamic budget: The computational budget of a Bayesian quadrature algorithm is fixed if the total number of integrand evaluations \(N\) is determined before the algorithm runs. This is the case if \(N\) is set by the practitioner by, for example, providing a fixed set of nodes \({\boldsymbol{x}}_1, \dots, {\boldsymbol{x}}_N\) or by providing \(N\) directly alongside a sampler. The budget is dynamic if the number of evaluations of \(f\) is not known before run time, and the algorithm determines an appropriate \(N\) automatically.²⁰ Methods with a dynamic budget are sometimes called automatic as for example in [10], [41], [158]. A dynamic budget does not imply that node selection is active. Whether a fixed or dynamic budget is preferable in practice depends on the requirements and constraints of the application.

• Multi-source: Multi-source Bayesian quadrature is an umbrella term that implies that not only one, but multiple, usually related integrals are being considered simultaneously. Several settings exist: Multi-fidelity BQ assumes a natural ordering of the integrals (e.g., from cheap to expensive) and usually aims to infer the integral with the highest fidelity (e.g., the most expensive one) from integrand evaluations of all fidelity levels. The ordering can either be encoded in the model directly [20] or implicitly via the sampling strategy [36]. Multi-source BQ [36] generalizes multi-fidelity BQ in that it aims to learn a primary integral but does not assume a natural ordering of the related integrals. [14] assume no primary integral and aim to learn all related integrals equally well. Active multi-source BQ [36] can, via its acquisition function, naturally be tailored towards aiming to learn a primary, some, or all related integrals. Other BQ methods that feature multiple integrals include Bayesian quadrature optimization methods which seek to optimise functions of the form \(g({\boldsymbol{y}}) = \int_Df({\boldsymbol{x}}, {\boldsymbol{y}}) \dif P({\boldsymbol{x}})\) [159]–[162] and the conditional BQ method of [163] that estimates integrals with respect to a parametric family of distributions \(\{ P_{\boldsymbol{\theta}}\}_{{\boldsymbol{\theta}}\in \Theta}\). While multi-source BQ does not strictly comply with Definition 2.1, resulting algorithms fit well into our taxonomy.

3.5 Example↩︎

Here we present an example of a BQ method and examine how modififying the method would affects its classification under the taxonomy. The method is only a conceptual example. We have not implemented it and there may be practical issues that render it impractical. Proper numerical experiments are found in 6.

Let \(D= [0, 1]^d\) and suppose that integration problem is to compute \(I_P(f) = \int_D f({\boldsymbol{x}}) \dif {\boldsymbol{x}}\) for some integrand function \(f\colon D \to \mathbb{R}\). Our illustrative BQ method is the following:

• Model: To model \(f\), we select the zero-mean GP model \(\mathsf{f}\sim \mathrm{GP}(0, k_{\boldsymbol{\theta}})\) whose covariance function is the square exponential kernel in 3 and \({\boldsymbol{\theta}}= (\sigma, \ell)\) are the kernel parameters. This choice means that the method is conjugate (3.1).

• Sampling: Let \(\{x_1, \ldots, x_n\}^d \subset [0, 1]^d\) be the \(d\)-ary Cartesian product of a set points \(\{x_1, \ldots, x_n \} \subset [0, 1]\). Select the number of iterations \(i_{\max} \geq 1\) and sample the nodes by iterating the following steps for \(i = 1, \ldots, i_{\max}\):

  1. Denote \(X(x) = \{x_1, \ldots, x_{i-1}, x\}^d \subset [0, 1]^d\) and solve the optimization problem \[\label{eq:acquisition-example} x_i = \mathop{\mathrm{arg\,max}}_{x \in [0, 1]} \rho^2(x)\tag{18}\] for \[\rho^2(x) = {\boldsymbol{k}}_{{\boldsymbol{\theta}},PX(x)}^{\intercal}{\boldsymbol{K}}_{{\boldsymbol{\theta}},X(x)X(x)}^{-1}{\boldsymbol{k}}_{{\boldsymbol{\theta}},PX(x)} / k_{{\boldsymbol{\theta}},PP}.\] That the quantities defined in 2.2.1 depend on \({\boldsymbol{\theta}}\) is made explicit here.

  2. Set \(X= X(x_i)\) and and note that \(X\) contains \(N_i = i^d\) points. In the notation of 6, and the batch size \(n = i^d - (i - 1)^d = N_i - N_{i-1}\) depends on the step.

  3. Use maximum likelihood to update the kernel parameters (see 4.3): \[{\boldsymbol{\theta}}= \mathop{\mathrm{arg\,max}}_{ \tilde{{\boldsymbol{\theta}}} \in \Theta} \bigg\{ -\frac{1}{2} {\boldsymbol{f}}_{X}^{\intercal}{\boldsymbol{K}}_{\tilde{{\boldsymbol{\theta}}},XX}^{-1} {\boldsymbol{f}}_{X} -\frac{1}{2} \log{\det {\boldsymbol{K}}_{\tilde{{\boldsymbol{\theta}}},XX}} \bigg\},\] where \(\Theta\) is a suitable subset of \([0, \infty) \times (0, \infty)\).

  The method is sequential and batched (3.3.3). It is also deterministic (3.3.1) because a perfect implementation would contain no randomness, though in practice there may be a degree of randomness in the optimization performed during Steps 1 and 3. Because \(x_i\) in Step 1 depends on the model \(\mathsf{f}\) via the covariance function \(k\), the method is model-dependent (3.3.2). In fact, as we shall see in 4.2.3, on each step the optimization problem 18 can be written in the form 17 , which makes the method active (and, in this case, myopic; see 3.3.3). Since there is no “direct” dependency on the observations \({\boldsymbol{f}}_X\) in Step 1, the method is not explicitly active. But because \(\rho^2\) depends on the kernel parameters, which in turn depend on \({\boldsymbol{f}}_X\) via Step 3, the method is implicitly active.

• Inference: After the nodes have been sampled, the method returns the posterior \(I_P(\mathsf{f}) \mid {\cal D}\sim \mathrm{N}(\mu_{\cal D}, \Sigma_{\cal D})\) by computing \(\mu_{\cal D}\) and \(\Sigma_{\cal D}\) from 8 and 10 . Because both \(k_{{\boldsymbol{\theta}},P}({\boldsymbol{x}})\) and \(k_{{\boldsymbol{\theta}},PP}\) are available in closed form for the square exponential kernel and the uniform measure on \([0, 1]^d\) (see 4.1), inference is exact (3.2) potential numerical problems are ignored. The method becomes scalable (3.2) but remains exact if one does not naively solve 8 and 10 but instead uses the well-known Kronecker decompositions of \({\boldsymbol{K}}_{{\boldsymbol{\theta}},XX}\) and \({\boldsymbol{k}}_{{\boldsymbol{\theta}},PX}\) for Cartesian product nodes and product kernels, such as the square exponential [5]. Sampling would no longer be model-dependent, active, or sequential should the method be modified so that \(X\) were simply some pre-specified product design, such as the \(d\)-ary product of equispaced points on \([0, 1]\). However, this sampling could be considered inference-dependent (3.3.2) as a product design is needed for scalable inference.

Next we turn to practical issues, many of them present in the method above, that one needs to negotiate in order to successfully design and implement a BQ method.

4 Practical issues↩︎

In this chapter we describe a number of practical issues that one encounters and how to solve them when implementing Bayesian quadrature algorithms. Some, such as numerical ill-conditioning and the need for kernel parameter estimation, need to be taken into account in all applications of Gaussian process regression and interpolation, while other, in particular the computation of kernel means and variances, are intrinsic to Bayesian quadrature. We conclude the section with an overview of what we view as the main limitations of Bayesian quadrature and how these have affected—and continue to affect—algorithmic development and the range of integration problems that Bayesian quadrature can or should solve.

4.1 Kernel means and variances↩︎

Conjugate Bayesian quadrature requires the computation of \(\mu_{{\cal D}}\) and \(\Sigma_{{\cal D}}\) as in Eq. 12 in order to characterize \(I_{P}(\mathsf{f})\mid{\cal D}\). Eq. 12 poses the challenge of computing the vector \({\boldsymbol{k}}_{PX}= ( k_P({\boldsymbol{x}}_1), \ldots k_P({\boldsymbol{x}}_N))\in\mathbb{R}^N\) of kernel means in Eq. 9 and the scalar kernel variance \(k_{PP}\) in Eq. 11 . A straightforward approach is to use a kernel \(k\) that admits an analytic solution of both \(k_{PP}\) and \(k_P({\boldsymbol{x}})\). A well-known example is the square exponential kernel in Eq. 3 integrated with respect to the standard Lebesgue measure on the domain \(D=[0, 1]^d\), which, for \(d=1\) yields \[\begin{align} k_{P}(x) &= \Big(\frac{\pi\ell^2}{2}\Big)^{\frac{1}{2}}\Big[ \erf\Big(\frac{1 - x}{\ell\sqrt{2}}\Big) + \erf\Big(\frac{x}{\ell\sqrt{2}}\Big) \Big], \\ k_{PP} & = 2\ell^2(e^{-\frac{1}{2\ell^2}}-1) + (2\pi\ell^2)^{\frac{1}{2}}\erf\Big(\frac{1}{\ell\sqrt{2}}\Big). \end{align}\] Comprehensive lists of closed form kernel embeddings can be found in [164] and [13]. The former is accompanied with a kernel mean embedding library in Python²¹. A variety of kernel embeddings are also implemented in the Python libraries EmuKit [130], [131] and ProbNum [133]. Approximations to non-analytic kernel embeddings for use in BQ have been less explored in the literature. Notable exceptions are [165] and [166], who use fast Fourier transforms to approximate kernel means of stationary kernels. See also [167] for a different approximation for isotropic Matérns.

Non-conjugate BQ often requires the computation of further quantities, such as \(\int_{D}\int_{D}k^2({\boldsymbol{x}},{\boldsymbol{x}}')\dif P({\boldsymbol{x}})\dif P({\boldsymbol{x}}')\), the integral of the kernel squared, and other combinations of kernel products [47]. Hence, methods based on non-conjugate BQ models are often specific to measure-kernel combinations which meet these requirements. Such integrals are also required in application of Bayesian quadrature to non-linear Kalman filtering [168]–[174].

4.2 Sampling and node selection↩︎

Section 3.3 described how the nodes may be sampled (i.e., selected). This section details practical choices and issues one is faced with when designing a sampling procedure for BQ.

4.2.1 Model-independent deterministic designs↩︎

If computational or other reasons prevent active sampling and better coverage of the domain than what is achievable via random sampling is desired, it is common to use deterministic designs that are model-independent²² (see Section 3.3.2). If the domain is a hypercube (say, \(D= [0, 1]^d\)), arguably the simplest deterministic design is the regular tensor grid \[\label{eq:regular-tensor-grid} \begin{align} X_{N_1}^d &= X_{N_1} \times \cdots \times X_{N_1} \\ &= \bigg\{ \frac{1}{N_1+1}, \ldots, \frac{N_1}{N_1+1} \bigg\} \times \cdots \times \bigg\{ \frac{1}{N_1+1}, \ldots, \frac{N_1}{N_1+1} \bigg\}, \end{align}\tag{19}\] which consists of \(N = N_1^d\) nodes. More advanced designs are obtained by replacing each one-dimensional node sets \(X_{N_1}\) with the nodes of an appropriate Gaussian quadrature rule, such as Gauss–Hermite if \(P\) is standard Gaussian on \(\mathbb{R}^d\) or Gauss–Legendre if \(P\) is uniform on \([0, 1]^d\). The resulting designs are sometimes called classical in the BQ literature. The computation of these nodes is efficiently implemented by most scientific computing libraries, including SciPy for Python. As the number of points depends exponentially on the dimension, tensor grids are inflexible and difficult to use in high dimensions. Sparse grids [94] as well as lattice nodes [10] and other designs from quasi-Monte Carlo literature may be used in high dimensions. Some of the above designs are depicted in Figure 9 for \(d \in \{2, 3\}\). It is a greate advantage of Bayesian quadrature that one may use any nodes that have been proposed in the literature.

In using non-nested nodes, many deterministic designs implicitly assume that the computational budget is fixed (see Section 3.4). For example, \(X_{N_1} \not\subset X_{N_1+1}\) if \(X_{N_1}\) are the equispaced nodes in 19 or any classical Gaussian quadrature nodes [176]. A non-nested design can be used as an initial design in sequential Bayesian quadrature where more nodes are sampled (e.g., actively) until some stopping criterion is met. Nested deterministic designs include the extensible integration lattices used in [10] and the Clenshaw–Curtis nodes [177] when \(N\) is a power of two (i.e., \(X_{2^n} \subset X_{2^{n+1}}\)).

4.2.2 Random nodes↩︎

An early instance of a BQ method—Bayesian Monte Carlo [7]—uses nodes that are random samples from the integration measure \(P\). As such it puts higher value on observing the integrand in areas where the density \(p\) of \(P\) is large than where it is small. This sampling method is model-uninformed. [135] recently proposed a sampling distribution for BQ that is inspired by importance sampling and also covers areas where \(f\) is supposedly large (even if the density \(p\) is small). This sampling method, in contrast, is model-dependent. Purposely sampling the nodes from a distribution different from \(P\) has also been proposed and studied [93], [97].

The space-filling Latin hypercube design [178] is a random, scalable equivalent to a grid. Being model-independent, it is used as initial design in sequential active BQ methods. Determinantal point processes, used in BQ by [102], provide model-dependent random designs constructed using the kernel function. Nodes sampled from a determinantal point process are repulsive and space-filling.

4.2.3 Actively sampled nodes↩︎

As promised in Section 3.3.4, this section reviews some commonly used utility functions for active BQ and their corresponding acquisition functions. For the conjugate model we give explicit expressions whenever possible. Recall from Section 3.3.4 that active nodes maximize an expected utility according to Eq. 17 . The equation is repeated here for convenience: \[\begin{align} {\boldsymbol{x}}_1,\dots, {\boldsymbol{x}}_N &\in \mathop{\mathrm{arg\,max}}_{\tilde{{\boldsymbol{x}}}_1,\dots, \tilde{{\boldsymbol{x}}}_N} \mathbb{E}_{\mathsf{f}}[u(\tilde{{\boldsymbol{x}}}_1,\dots, \tilde{{\boldsymbol{x}}}_N\mid\mathsf{f}, P)]. \end{align}\] Again, the function \(u\) is the utility function and \[a(\tilde{{\boldsymbol{x}}}_1,\dots, \tilde{{\boldsymbol{x}}}_N\mid P) \mathrel{\vcenter{:}}= \mathbb{E}_{\mathsf{f}}[u(\tilde{{\boldsymbol{x}}}_1,\dots, \tilde{{\boldsymbol{x}}}_N\mid\mathsf{f}, P)]\] the corresponding acquisition function. The above equation, at heart, considers what would happen under the model \(\mathsf{f}\) if the integrand \(f\) were hypothetically evaluated at some yet unseen nodes \(\tilde{{\boldsymbol{x}}}_1,\dots, \tilde{{\boldsymbol{x}}}_N\), expressed quantitatively as an average utility \(a\). The active sampler returns the maximizer of \(a\). In other words, the sampler produces nodes that yield the maximal average utility under the model \(\mathsf{f}\). A question of course arises: How to choose \(u\)?

From an information theoretic perspective, a meaningful acquisition function is the mutual information (MI) between the integral \(I_{P}(\mathsf{f})\) and yet unseen integrand values \({\boldsymbol{\mathsf{f}}}_{\tilde{X}} \mathrel{\vcenter{:}}= (\mathsf{f}(\tilde{{\boldsymbol{x}}}_1),\dots, \mathsf{f}(\tilde{{\boldsymbol{x}}}_N))\in\mathbb{R}^N\) at nodes \(\tilde{X} = \{\tilde{{\boldsymbol{x}}}_1,\dots, \tilde{{\boldsymbol{x}}}_N\}\). Let \(\tilde{{\cal D}} = \{ (\tilde{{\boldsymbol{x}}}_1, \mathsf{f}(\tilde{{\boldsymbol{x}}}_1)), \ldots (\tilde{{\boldsymbol{x}}}_N, \mathsf{f}(\tilde{{\boldsymbol{x}}}_N)) \}\) be the corresponding hypothetical dataset. The MI acquisition function is \[\begin{align} a(\tilde{{\boldsymbol{x}}}_1,\dots, \tilde{{\boldsymbol{x}}}_N\mid P) &= \mathrm{MI}(I_{P}(\mathsf{f}); {\boldsymbol{\mathsf{f}}}_{\tilde{X}}) \\ &\mathrel{\vcenter{:}}= \mathbb{E}_{\mathsf{f}}[ \mathrm{KL}\kldiv{ {\boldsymbol{\mathsf{f}}}_{\tilde{X}}\mid I_{P}(\mathsf{f}), \tilde{X} }{{\boldsymbol{\mathsf{f}}}_{\tilde{X}}\mid\tilde{X}} ] \\ &= \mathbb{E}_{{\boldsymbol{\mathsf{f}}}_{\tilde{X}}}[\mathrm{KL}\kldiv{I_{P}(\mathsf{f})\mid \tilde{{\cal D}}}{I_{P}(\mathsf{f})}], \end{align}\] where \(\mathrm{KL}\kldiv{Q}{Q'}\) is the Kullback–Leibler (KL) divergence between two distributions \(Q\) and \(Q'\). The MI acquisition function measures the expected change in distribution between the integral \(I_{P}(\mathsf{f})\) over the prior model and the integral \(I_{P}(\mathsf{f})\mid\tilde{{\cal D}}\) over the model that is condition on the hypothetical dataset \(\tilde{{\cal D}}\). Maximizing it yields maximally informative nodes (in the Shannon sense) about the integral value under the model \(\mathsf{f}\). For conjugate Bayesian quadrature, the mutual information is analytic and has the form \[\begin{align} \label{eq:mi-conj} \mathrm{MI}(I_{P}(\mathsf{f}); {\boldsymbol{\mathsf{f}}}_{\tilde{X}}) = -\frac{1}{2}\log( 1 - \rho^2(\tilde{{\boldsymbol{x}}}_1,\dots, \tilde{{\boldsymbol{x}}}_N) ), \end{align}\tag{20}\] where \[\begin{align} \label{eq:rho2} \rho^2(\tilde{{\boldsymbol{x}}}_1,\dots, \tilde{{\boldsymbol{x}}}_N) \mathrel{\vcenter{:}}= ({\boldsymbol{k}}_{P\tilde{X}}^{\intercal}{\boldsymbol{K}}_{\tilde{X}\tilde{X}}^{-1}{\boldsymbol{k}}_{P\tilde{X}})/k_{PP} \in [0, 1] \end{align}\tag{21}\] is the scalar squared correction. Here \(k_{PP}\) is the initial variance defined in Eq. 11 , \({\boldsymbol{k}}_{P\tilde{X}} = ( k_P(\tilde{{\boldsymbol{x}}}_1), \ldots k_P(\tilde{{\boldsymbol{x}}}_N)) \in \mathbb{R}^N\) is a vector, and \({\boldsymbol{K}}_{\tilde{X}\tilde{X}} = (k(\tilde{{\boldsymbol{x}}}_i, \tilde{{\boldsymbol{x}}}_j))_{i,j=1}^N \in \mathbb{R}^{N \times N}\) a matrix. For example [179] and [180] use information-theoretic acquisition functions.

Other frequently used BQ acquisition functions include the integral variance reduction (IVR) and the negative integral variance (NIV) mentioned in 3.3.4 that arise from the square error utility: \[\begin{align} \mathrm{IVR}(&I_{P}(\mathsf{f}); I_{P}(\mathsf{f})\mid\tilde{{\cal D}}) \\ &\mathrel{\vcenter{:}}= s^{-1}\mathbb{E}_{\mathsf{f}}\big[ ( I_P(\mathsf{f}) - \mathbb{E}_{\mathsf{f}}[ I_P(\mathsf{f})] )^2 - ( I_P(\mathsf{f}) \mid \tilde{{\cal D}} - \mathbb{E}_{\mathsf{f}}[ I_P(\mathsf{f}) \mid \tilde{{\cal D}} ] )^2 \big]\\ &= s^{-1}(s - \mathrm{Var}[ I_P(\mathsf{f}) \mid \tilde{{\cal D}} ]) \end{align}\] and \[\begin{align} \mathrm{NIV}(I_{P}(\mathsf{f})\mid\tilde{{\cal D}}) & \mathrel{\vcenter{:}}= -s^{-1}\mathbb{E}_{\mathsf{f}}\big[ ( I_P(\mathsf{f}) \mid \tilde{{\cal D}} - \mathbb{E}_{\mathsf{f}}[ I_P(\mathsf{f}) \mid \tilde{{\cal D}}] )^2 \big]\\ & = - s^{-1}\mathrm{Var}[ I_P(\mathsf{f}) \mid \tilde{{\cal D}} ]. \end{align}\] The inverse of the factor \(s = \mathrm{Var}[ I_P(\mathsf{f})]\) (not necessarily equal to \(k_{PP}\) for arbitrary \(\mathsf{f}\)) is a positive constant with respect to \(\tilde{X}\) and serves as normalization. The IVR is simply the (normalized) difference between the variance \(s\) of the integral over the prior model and the variance of the integral over the model conditioned on the hypothetical dataset \(\tilde{{\cal D}}\). Hence IVR selects nodes that shrink the variance of the integral maximally. Analogously, NIV chooses nodes that minimize the variance of the integral. The IVR and NIV yield the same maximizer and hence the same active nodes as both functions are strictly increasing transformations of each other. For the conjugate model, these two acquisition functions are analytic and simplify to \[\begin{align} {3} \tag{22} \mathrm{IVR}(I_{P}(\mathsf{f}); I_{P}(\mathsf{f})\mid\tilde{{\cal D}}) &= \rho^2(\tilde{{\boldsymbol{x}}}_1,\dots, \tilde{{\boldsymbol{x}}}_N)\quad &&\in [0, 1], \\ \tag{23} \mathrm{NIV}(I_{P}(\mathsf{f})\mid\tilde{{\cal D}}) &= \rho^2(\tilde{{\boldsymbol{x}}}_1,\dots, \tilde{{\boldsymbol{x}}}_N) - 1 \quad &&\in [-1, 0]. \end{align}\] By comparing Eqs. 20 , 22 , and 23 it is apparent that, in the conjugate case, all three acquisition functions (IVR, NIV, MI) can be expressed as functions of the squared correlation \(\rho^2\) and as strictly increasing transformations of each other. Hence all three yield the same nodes [36]. Figure 10 shows the acquisitions as functions of \(\rho^2\).²³

For sequential active sampling as introduced in Algorithms 6 and 7 in Sections 3.3.3 and 3.3.4, the formulas need only slight modification. Recall that a sequential sampler produces a small batch of nodes of size \(n\) at every step of the iteration which is then added to the current dataset. Hence, denote the current dataset as \({\cal D}\) and the hypothetical batch dataset as \(\tilde{{\cal D}} = \{ (\tilde{{\boldsymbol{x}}}_1, \mathsf{f}(\tilde{{\boldsymbol{x}}}_1)), \ldots (\tilde{{\boldsymbol{x}}}_n, \mathsf{f}(\tilde{{\boldsymbol{x}}}_n)) \}\), which overloads notation in this paragraph. Further, denote the current posterior as \(\mathsf{f}_{{\cal D}}=\mathsf{f}\mid{\cal D}\). The acquisition function can be defined autoregressively (meaning that \(\mathsf{f}_{{\cal D}}\) acts as the prior model in the current step) according to \[\begin{align} \mathrm{MI}(I_{P}(\mathsf{f}_{{\cal D}}); {\boldsymbol{\mathsf{f}}}_{\tilde{X}}) &= \mathbb{E}_{{\boldsymbol{\mathsf{f}}}_{\tilde{X}}}\big[\mathrm{KL}\kldiv{ I_{P}(\mathsf{f}_{{\cal D}})\mid \tilde{{\cal D}}}{ I_{P}(\mathsf{f}_{{\cal D}})}\big], \\ \mathrm{IVR}(I_{P}(\mathsf{f}_{{\cal D}}); I_{P}(\mathsf{f}_{{\cal D}})\mid\tilde{{\cal D}}) &= s^{-1}_{{\cal D}}(s_{{\cal D}} - \mathrm{Var}[ I_P(\mathsf{f}_{{\cal D}}) \mid \tilde{{\cal D}} ]), \\ \mathrm{NIV}(I_{P}(\mathsf{f}_{{\cal D}})\mid\tilde{{\cal D}}) & = - s^{-1}_{{\cal D}}\mathrm{Var}[ I_P(\mathsf{f}_{{\cal D}}) \mid \tilde{{\cal D}} ], \end{align}\] with \(s_{{\cal D}} = \mathrm{Var}[ I_P(\mathsf{f}_{{\cal D}})]\) again a constant with respect to \(\tilde{X}\) that only serves to normalize the acquisition functions. The acquisition functions are myopic by definition as they do not consider the impact of future steps. For the conjugate model the autoregressive squared correlation is \[\begin{align} \label{eq:rho2-seq} \rho^2(\tilde{{\boldsymbol{x}}}_1,\dots, \tilde{{\boldsymbol{x}}}_n) &= ({\boldsymbol{k}}_{{\cal D}; P\tilde{X}}^{\intercal}{\boldsymbol{K}}_{{\cal D};\tilde{X}\tilde{X}}^{-1}{\boldsymbol{k}}_{{\cal D}; P\tilde{X}})/\Sigma_{{\cal D}}, \end{align}\tag{24}\] where \({\boldsymbol{k}}_{{\cal D}; P\tilde{X}} = I_P({\boldsymbol{k}}_{{\cal D}; \tilde{X}}) = ( k_{{\cal D}; P}(\tilde{{\boldsymbol{x}}}_1), \ldots k_{{\cal D}; P}(\tilde{{\boldsymbol{x}}}_n))\) is the kernel mean embedding of the posterior kernel function as in Eq. 2 as a function of the \(n\) hypothetical nodes, \({\boldsymbol{K}}_{{\cal D};\tilde{X}\tilde{X}}\) is the corresponding Gram matrix, and \(\Sigma_{\cal D}\) is the integral variance in Eq. 10 . For the prevalent case \(n=1\), the sequential squared correlation simplifies to \[\rho^2(\tilde{{\boldsymbol{x}}}) = \frac{k_{{\cal D};P}(\tilde{{\boldsymbol{x}}})^2}{\Sigma_{{\cal D}}k_{{\cal D}}(\tilde{{\boldsymbol{x}}}, \tilde{{\boldsymbol{x}}})}.\]

We mention two acquisition functions that are only indirectly related to the integral. Uncertainty sampling (US), popular in experimental design, is defined for \(n=1\) only by \[\begin{align} \mathrm{US}(\mathsf{f}_{{\cal D}}) \mathrel{\vcenter{:}}= \mathbb{E}_{\mathsf{f}}\big[ ( \mathsf{f}_{{\cal D}}(\tilde{{\boldsymbol{x}}}) - \mathbb{E}_{\mathsf{f}_{{\cal D}}}[\mathsf{f}_{{\cal D}}(\tilde{{\boldsymbol{x}}}) ] )^2p^2(\tilde{{\boldsymbol{x}}}) \big] = \mathrm{Var}[\mathsf{f}_{{\cal D}}(\tilde{{\boldsymbol{x}}})]p^2(\tilde{{\boldsymbol{x}}}). \end{align}\] For the conjugate model we get \(\mathrm{US}(\mathsf{f}_{{\cal D}}) = k_{{\cal D}}(\tilde{{\boldsymbol{x}}}, \tilde{{\boldsymbol{x}}})p^2(\tilde{{\boldsymbol{x}}})\). Uncertainty sampling yields nodes at locations where the variance of \(\mathsf{f}_{{\cal D}}\) and the density of the integration measure \(p(\tilde{{\boldsymbol{x}}})\) are both large. By focusing on the variance of the model \(\mathsf{f}\) instead of its integral \(I_{P}(\mathsf{f})\), US in a sense disregards the integral as the quantity of interest and may hence be more tractable or accessible than many other acquisition functions. Uncertainty sampling is for example used in [47] and in a related (non-active) context in [135]. The posterior variance contribution (PVC) acquisition function [181] \[\begin{align} \mathrm{PVC}(\mathsf{f}_{{\cal D}}) \mathrel{\vcenter{:}}= k_{{\cal D}; P}(\tilde{{\boldsymbol{x}}}) p(\tilde{{\boldsymbol{x}}}) \end{align}\] is a heuristic defined only for the conjugate model. Like US, \(\mathrm{PVC}\) yields nodes likely to contribute much to the integral variance, but ignores spatial correlations between candidate nodes. One can also proceed to the opposite direction and use an acquisition function developed for Bayesian quadrature to select evaluation locations for a purpose other than integration [182]–[185].

Implementations of sequential active sampling loops and BQ acquisition functions for \(n=1\) can for example be found in the Python library EmuKit [130], [131] which we use in the experiments in Chapter 6.

Recall from 2.3 that BQ is equivalent to worst-case optimal integration in the RKHS \(\mathcal{H}(k)\) of the covariance function. This and related equivalences allow expressing all acquisition functions above as quantities involving different worst-case errors. For example, we saw that the MI, IVR and NIV acquisition functions yield nodes that maximize the squared correlation in Eq. 21 Equivalently, they minimize the integral variance \(\Sigma_{\tilde{{\cal D}}} = k_{PP} - {\boldsymbol{k}}_{P\tilde{X}}^{\intercal}{\boldsymbol{K}}_{\tilde{X}\tilde{X}}^{-1} {\boldsymbol{k}}_{P\tilde{X}}\). By Proposition [prop:wce], for given nodes the integral variance equals the smallest possible squared worst-case error in \(\mathcal{H}(k)\) obtainable by a quadrature rule using these nodes. Nodes \({\boldsymbol{x}}_1, \ldots, {\boldsymbol{x}}_N\) sampled by MI, IVR or NIV are thus worst-case optimal: \[\label{eq:obq-nodes} \begin{align} {\boldsymbol{x}}_1, \ldots, {\boldsymbol{x}}_N \in \mathop{\mathrm{arg\,min}}_{\tilde{{\boldsymbol{x}}}_1, \ldots, \tilde{{\boldsymbol{x}}}_N} \, \min_{ {\boldsymbol{v}}\in \mathbb{R}^N} e_k({\boldsymbol{v}}, \tilde{X}), \end{align}\tag{25}\] where \[e_k({\boldsymbol{v}}, \tilde{X}) = \sup_{ \norm[0]{g}_{\mathcal{H}(k)} \leq 1} \, \abs[3]{\int_Dg({\boldsymbol{x}}) \dif P({\boldsymbol{x}}) - \sum_{i=1}^N v_i g(\tilde{{\boldsymbol{x}}}_i)} .\] These nodes provide maximal information for integration in \(\mathcal{H}(k)\). Accordingly, optimal Bayesian quadrature is sometimes used to refer to BQ with the nodes in Eq. 25 . These nodes are called optimal also in numerical analysis and approximation theory, where the corresponding worst-case error is known as the \(N\)th minimal error; see Chapters 2 and 3 in [186] and Chapter 5 in [58].

4.3 Kernel parameter estimation↩︎

In Bayesian inference one is confronted with the challenging task of selecting an appropriate prior. While other options exist, in Gaussian process modelling the prevalent approach is to select a parametrised kernel, such as the square exponential in 3 or a Matérn kernel in 4 with a fixed smoothness parameter \(\nu\), and estimate its parameters from the data [23]. Most Bayesian quadrature methods use maximum likelihood estimation to set the kernel parameters. Marginalisation and cross-validation are also occasionally used [9], [10]. As our experiments in 6 use maximum likelihood estimation, we briefly review it in the context of conjugate Bayesian quadrature; that is, when \(\mathsf{f}= \mathsf{g}\sim \mathrm{GP}(m, k_{\boldsymbol{\theta}})\) for a kernel with parameters \({\boldsymbol{\theta}}\). As the conjugate Bayesian quadrature integral variance in 11 does not directly depend on the integrand evaluations \({\boldsymbol{f}}_X\), estimation of the kernel or its parameters from the data is essential if the variance is to measure integration uncertainty.

Let \({\boldsymbol{\theta}}\) be the kernel parameters and let the zeroth parameter \(\theta_0 = \sigma^2\) denote the kernel scale. Let \({\cal D}\) be the dataset consisting of the nodes \(X= \{ {\boldsymbol{x}}_1, \ldots, {\boldsymbol{x}}_N \}\) and corresponding observations \({\boldsymbol{f}}_{X}\). Maximum likelihood (ML) estimation is a popular way to fit kernel parameters. The ML estimate of \({\boldsymbol{\theta}}\) is \[\begin{align} \hat{{\boldsymbol{\theta}}}_\mathrm{ML} = \mathop{\mathrm{arg\,max}}_{{\boldsymbol{\theta}}} L({\boldsymbol{\theta}}; {\cal D}) &= \mathop{\mathrm{arg\,max}}_{{\boldsymbol{\theta}}} p({\boldsymbol{f}}_{X}\mid X, {\boldsymbol{\theta}}), \end{align}\] where \(L({\boldsymbol{\theta}}; {\cal D})\) is the likelihood of the parameters. From a probabilistic perspective, the ML estimator maximizes the probability \(p({\boldsymbol{f}}_{X}\mid X, {\boldsymbol{\theta}})\) that the given dataset is observed under the marginalized model \(\mathsf{f}\). Equivalently, the estimator \(\hat{{\boldsymbol{\theta}}}_{ML}\) maximizes the log-likelihood \(\log L({\boldsymbol{\theta}}; {\cal D})\). For the conjugate BQ model the logarithm admits the closed form \[\label{eq:hyper-likelihood} \begin{align} \log L({\boldsymbol{\theta}}; {\cal D}) ={}& -\frac{1}{2}({\boldsymbol{f}}_{X} - {\boldsymbol{m}}_{X})^{\intercal}{\boldsymbol{K}}_{XX}^{-1}({\boldsymbol{\theta}})({\boldsymbol{f}}_{X} - {\boldsymbol{m}}_{X}) \\ &-\frac{1}{2}\log{\det({\boldsymbol{K}}_{XX}({\boldsymbol{\theta}}))} -\frac{N}{2}\log{2\pi}, \end{align}\tag{26}\] where \(\det({\boldsymbol{K}}_{XX}({\boldsymbol{\theta}}))\) is the determinant of the Gram matrix \({\boldsymbol{K}}_{XX}({\boldsymbol{\theta}})\) and the dependence on \({\boldsymbol{\theta}}\) has been made explicit. Combining the two equations yields \[\label{eq:hyper-ml} \begin{align} \hat{{\boldsymbol{\theta}}}_\mathrm{ML} ={}& \mathop{\mathrm{arg\,max}}_{{\boldsymbol{\theta}}} \bigg\{ -\frac{1}{2}({\boldsymbol{f}}_{X} - {\boldsymbol{m}}_{X})^{\intercal}{\boldsymbol{K}}_{XX}^{-1}({\boldsymbol{\theta}})({\boldsymbol{f}}_{X} - {\boldsymbol{m}}_{X}) \\ &-\frac{1}{2}\log{\det({\boldsymbol{K}}_{XX}({\boldsymbol{\theta}}))} \bigg\}. \end{align}\tag{27}\] In general, an analytic solution of Eq. 27 only exists for the scale parameter \[\hat{\sigma}^2_\mathrm{ML} = \frac{1}{N}({\boldsymbol{f}}_{X} - {\boldsymbol{m}}_{X})^{\intercal}{\boldsymbol{K}}_{XX}^{-1}({\boldsymbol{\theta}}_0)({\boldsymbol{f}}_{X} - {\boldsymbol{m}}_{X}),\] where \({\boldsymbol{\theta}}_0 = (\theta_0 = 1, \theta_1, \ldots)\), so that \({\boldsymbol{K}}_{XX}({\boldsymbol{\theta}}_0)\) is the unscaled kernel Gram matrix. If all parameters \({\boldsymbol{\theta}}\) need to be estimated, Eq. 27 can be solved by an off-the-shelf optimizer as gradients of \(\log L({\boldsymbol{\theta}}; {\cal D})\) with respect to \({\boldsymbol{\theta}}\) are usually available.

The kernel parameters \({\boldsymbol{\theta}}\) are also called hyper-parameters to distinguish them from the (possibly infinite) parameters of the model \(\mathsf{f}\). Thus, maximizing \(L({\boldsymbol{\theta}}; {\cal D})\) is also called maximum likelihood type II (instead of the standard type I). Other common descriptors are maximum marginal likelihood [since \(p({\boldsymbol{f}}_{X}\mid X, {\boldsymbol{\theta}})\) describes a marginal probability], or evidence maximization [as \(p({\boldsymbol{f}}_{X}\mid X, {\boldsymbol{\theta}})\) is the evidence of the probabilistic model].

As \(\hat{{\boldsymbol{\theta}}}_\mathrm{ML}\) is an estimator based on a limited dataset \({\cal D}\) overfitting can occur. As a consequence, BQ methods with sequential active node selection as in [alg:bq-seq,alg:active-sampl] that compute the ML-estimate \(\hat{{\boldsymbol{\theta}}}_\mathrm{ML}\) at every step anew require initial designs of appropriate size that are independent of the model \(\mathsf{f}\). This is conceptually analogous to initial designs that are required by sequential active experimental design algorithms and Bayesian optimization.

4.4 Numerical ill-conditioning↩︎

In practically all variants of Bayesian quadrature the computation of the posterior \(I_P(\mathsf{f}) \mid {\cal D}\) and the fitting of the kernel parameters require one to solve—repeatedly in the case of parameter fitting—the linear system \({\boldsymbol{K}}_{XX} {\boldsymbol{z}}= {\boldsymbol{a}}\) for some vector \({\boldsymbol{a}}\in \mathbb{R}^N\). Unless the kernel or the nodes have properties that can be exploited to simplify this task, one has to use some standard linear solver, the accuracy of which depends on floating-point accuracy and the magnitude of the condition number of the kernel Gram matrix, \[\mathrm{cond}({\boldsymbol{K}}_{XX}) = \frac{\lambda_\mathrm{max}({\boldsymbol{K}}_{XX})}{\lambda_\mathrm{min}({\boldsymbol{K}}_{XX})} .\] Here \(\lambda_\mathrm{max}({\boldsymbol{A}})\) and \(\lambda_\mathrm{min}({\boldsymbol{A}})\) denote the maximal and minimal eigenvalues of a symmetric positive-definite matrix \({\boldsymbol{A}}\). When the condition number is large, solving the linear system is prone to numerical errors and hence numerically unstable. The condition numbers of commonly used kernels can grow fast when the nodes fill the domain. Let \(D= [0, 1]^d\) and \[h = \min_{i = 1, \ldots, N } \, \sup_{ {\boldsymbol{x}}\in [0, 1]^d } \, \min_{j \neq i} \norm[0]{{\boldsymbol{x}}- {\boldsymbol{x}}_j},\] which quantifies how well the nodes fill the domain, being the radius of the largest open ball in \([0, 1]^d\) that intersects no subset of \(X\) containing \(N-1\) nodes. This quantity is closely related to the fill-distance that we shall use in 5; see Eq. 31 . The following theorem that links \(h\) to \(\mathrm{cond}({\boldsymbol{K}}_{XX})\) may be found in [187]. See also [78], [188], and [189]. The general principle is that the smoother the kernel is, the faster the condition number grows as \(h\) tends to zero, which is to say that the nodes fill the domain.

The assumption \(D= [0, 1]^d\) can be relaxed as described in Remark [rmk:domain-generalisations]. Numerical ill-conditioning issues are likewise present in Gaussian process interpolation and kriging, where they are usually dealt with by introducing a parameter \(\lambda > 0\) called the nugget (alternatively, regularisation parameter or jitter) and replacing each occurrence of the Gram matrix \({\boldsymbol{K}}_{XX}\) with \({\boldsymbol{K}}_{XX} + \lambda {\boldsymbol{I}}_N\), where \({\boldsymbol{I}}_N\) is the \(N \times N\) identity matrix ([190]; [23, Sec. 3.4.3]; [191], [192]). [193] are strongly in favour of nuggets. If the kernel is stationary and \(k({\boldsymbol{x}}, {\boldsymbol{x}}) = 1\) for every \({\boldsymbol{x}}\in D\), it is straightforward to compute that \[\mathrm{cond}( {\boldsymbol{K}}_{XX} + \lambda {\boldsymbol{I}}_N) = \frac{\lambda_\mathrm{max}({\boldsymbol{K}}_{XX}) + \lambda}{\lambda_\mathrm{min}({\boldsymbol{K}}_{XX}) + \lambda} \leq \frac{N + \lambda}{\lambda},\] which means that the condition number may increase at most linearly in the number of nodes. In the context of conjugate Bayesian quadrature, the use of a nugget corresponds to replacing the mean \(\mu_{\cal D}\) and variance \(\Sigma_{\cal D}\) in 8 and 10 with their regularised versions \[\begin{align} \mu_{{\cal D}, \lambda} &= m_P+ I_P({\boldsymbol{k}}_X)^{\intercal}({\boldsymbol{K}}_{XX} + \lambda {\boldsymbol{I}}_N )^{-1} ( {\boldsymbol{f}}_X- {\boldsymbol{m}}_X), \label{eq:bq-mean-regularised} \\ \Sigma_{{\cal D}, \lambda} &= k_{PP} - {\boldsymbol{k}}_{PX}^{\intercal}({\boldsymbol{K}}_{XX} + \lambda {\boldsymbol{I}}_N)^{-1} {\boldsymbol{k}}_{PX}. \nonumber \end{align}\tag{28}\] The nugget can be seen as an additional hyper-parameter and fitted alongside other parameters as described in Section 4.3.

Theorem 4.1 shows that nuggets are necessary when \(N\) is sufficiently large. However, in 5 we shall see that the nugget has a limiting effect on the accuracy of conjugate Bayesian quadrature, which suggests that its use ought to be avoided when very high accuracy is being sought. That is, there is a trade-off between numerical stability and accuracy known as Schaback’s uncertainty principle [187].²⁴ Whether to use a nugget in Bayesian quadrature therefore depends on the modelling choices, the desired accuracy of the particular method, and the required level of code robustness. Most GP libraries hardcode a minimal default nugget that may depend on the type (at the time of writing GPy uses \(\lambda = 10^{-8}\) while GPtorch uses \(\lambda = 10^{-6}\) for floats and \(\lambda = 10^{-8}\) for doubles) plus possibly a dynamic nugget that depends on the Gram matrix \({\boldsymbol{K}}_{XX}\) and is increased successively if solving the linear system fails. 6.1.3 describes how nuggets are used in our experiments.

4.5 Prior mean estimation↩︎

The prevalent approach in the literature on BQ is to use a zero-mean prior \(\mathsf{g}\sim \mathrm{GP}(0, k)\). However, like the kernel or its parameters, also the prior mean \(m\) can be estimated from the data. Suppose that we are in the conjugate setting described in Section 2.2.1 where \(\mathsf{f}= \mathsf{g}\). The prior mean is typically estimated by postulating a parametric form \[m({\boldsymbol{x}}) = \sum_{i=1}^q \beta_i \, g_i({\boldsymbol{x}})\] for \(q \in \mathbb{N}\) fixed basis functions \(g_i \colon D\to \mathbb{R}\) and coefficients \({\boldsymbol{\beta}}= (\beta_1, \ldots, \beta_q) \in \mathbb{R}^q\) which are assigned some prior.²⁵ A common approach is to place a non-informative prior \(p({\boldsymbol{\beta}}, \sigma^2) \propto \sigma^{-2}\) on \({\boldsymbol{\beta}}\) and the kernel scale parameter \(\sigma^2\) and subsequently marginalize out these parameters [5]. Doing so yields an integral posterior which is Student’s \(t\) with \(N - q\) degrees of freedom. The mean \(\tilde{\mu}_{\cal D}\) and variance \(\tilde{\Sigma}_{\cal D}\) of the posterior have closed form expressions similar to those of \(\mu_{\cal D}\) and \(\Sigma_{\cal D}\) in Section 2.2.1, but they also involve evaluations of the basis functions \(g_i\) at the nodes and their integrals [5].

When only an unknown constant trend is postulated (i.e., \(m \equiv \beta_1\)) the posterior mean is \(\tilde{\mu}_{\cal D}= {\boldsymbol{w}}_{\mathbf{1}}^{\intercal}{\boldsymbol{f}}_X\), where \[\label{eq:normalised-bq-mean} {\boldsymbol{w}}_{\mathbf{1}} = \big( {\boldsymbol{I}}_N - \frac{{\boldsymbol{K}}_{XX}^{-1} \mathbf{1}_N \mathbf{1}_N^{\intercal}}{\mathbf{1}_N^{\intercal}{\boldsymbol{K}}_{XX}^{-1} \mathbf{1}_N } \Big) {\boldsymbol{K}}_{XX}^{-1} {\boldsymbol{k}}_{PX} + \frac{{\boldsymbol{K}}_{XX}^{-1} \mathbf{1}_N^{\intercal}}{\mathbf{1}_N^{\intercal}{\boldsymbol{K}}_{XX}^{-1} \mathbf{1}_N }\tag{29}\] and \(\mathbf{1}_N\) is the \(N\)-vector of ones, and one obtains normalized Bayesian quadrature. The terminology is motivated by the weights \({\boldsymbol{w}}_{\mathbf{1}}\) in 29 summing up to one, which does not generally happen in conjugate Bayesian quadrature with a fixed prior mean. Weight normalization is particularly useful as a precaution against poorly set length-scale in high dimensions. If the length-scale \(\ell\) of a stationary kernel, such as square exponential in 3 or Matérn in 4 , is very small relative to the typical distance \(\norm[0]{{\boldsymbol{x}}_i - {\boldsymbol{x}}_j}\) of the nodes, the covariance matrix is almost diagonal and \({\boldsymbol{w}}_{\mathbf{1}} \approx (\frac{1}{N}, \ldots, \frac{1}{N})\). In contrast, the unnormalized weights in 12 are very close to zero in this case, which means that the posterior is close to the prior.

Prior mean estimation can be alternatively achieved via the use of a flat prior limit [40] or maximum likelihood estimation [10]. When maximum likelihood estimation is used to set the covariance scale parameter \(\sigma^2\) (see Section 4.3), the resulting posteriors, except for being Gaussians, differ little from the Student’s \(t\) posterior described above; many relevant posteriors are collected in [197]. In addition to the works already mentioned, prior mean estimation or normalized Bayesian cubature has been used in [6], [126], [198]. We also refer to [80] and [199] for non-probabilistic versions, whose relation to Bayesian quadrature has been described in Section 2.3, of prior mean estimation. Finally, we note that the theoretical guarantees available in the literature usually assume that the prior mean is fixed.

4.6 Main limitations↩︎

Several factors that limit the applicability of Bayesian quadrature emerge from the practical issues reviewed in [sec:kern-mean-var,sec:sampling,sec:param-estimation,sec:num-ill-cond,sec:prior-mean-estimation]. The main limitations that have restricted Bayesian quadrature to a relatively narrow range of applications are (i) the need to have access to kernel means and variances and (ii) the cubic computational complexity in the number of nodes associated to Gaussian process methods. Mainly due to these two challenges a typical application of Bayesian quadrature features a small number of evaluations of a low-dimensional integrand and a simple integration measure.

4.6.1 Kernel means and variances↩︎

That the kernel mean \(k_P\) and variance \(k_{PP}\) are available in closed forms for a limited collection of kernel-measure pairs is a severe limitation. Most Bayesian quadrature applications involve integration with respect to the uniform or Gaussian measure, and sometimes Bayesian quadrature is interpreted in a narrow sense to mean an algorithm or a class of algorithms for integration with respect to Gaussian measures. In this respect Bayesian quadrature resembles Gaussian quadratures and other classical approaches that are designed for integration with respect to particularly tractable measures.

While [19] have used Dirichlet process mixtures to extend Bayesian quadrature for intractable integration distributions that can only be sampled from, their approach has not caught on and is bound to sacrifice many desirable properties of Bayesian quadrature, such as a straightforward theoretical guarantees. An approach that is in principle general and powerful is to use the measure to design a kernel for which \(k_P\) and \(k_{PP}\) are easy to compute. [200] and [201] use Stein’s method to design such kernels. However, this simply transfers the trouble from computing integrals to computing a kernel and results in a model that is not easy to interpret. We believe that no satisfactory solution can exist that would obviate the need to compute kernel means or other typically intractable integrals. To construct a quadrature rule that encodes some properties of the integrand one needs access to some functions that effect the encoding. When \(N\) is small (e.g., the integrand is a computationally expensive simulation), it may in some cases make sense to approximate \(k_P\) and \(k_{PP}\) with, for example, Monte Carlo, if doing so has negligible computational cost in comparison to the integrand evaluations.

4.6.2 Cubic computational complexity↩︎

Computational cost is the second major limitation. As discussed in 4.4, to implement Bayesian quadrature one has to solve linear systems defined by the typically dense \(N \times N\) kernel matrix \({\boldsymbol{K}}_{XX}\). The cubic computational complexity in the number of datapoints that this gives rise to is a well-known computational bottleneck in all Gaussian process methods. Over the past 20 years much effort has gone into the development of computationally scalable Gaussian process methods [141]. These methods are predominantly approximate in the sense of the definition given in 3.2 and developed under the assumption that the data are noisy, so that the linear systems that need to be solved are defined by the regularised kernel matrix \({\boldsymbol{K}}_{XX} + \lambda {\boldsymbol{I}}_N\) rather than \({\boldsymbol{K}}_{XX}\). Because the nugget \(\lambda\) is typically either zero or very small in Bayesian quadrature, many scalable Gaussian process methods cannot be imported to Bayesian quadrature as such. For example, methods based on low-rank approximations of the kernel matrix require a positive nugget. It is nevertheless surprising that no systematic study is available on the feasibility of applying existing approaches to scalable Gaussian process regression to Bayesian quadrature.

There exist a number of scalable approaches specifically to Bayesian quadrature (see the scalable label in 3.2). While some have been highly successful, all are tied to certain measure-kernel pairs and do not generalise easily or at all. Except for some special cases [10], Bayesian quadrature is presently not a serious competitor to Monte Carlo, quasi-Monte Carlo, and related approaches to numerical integration in high dimensions or in other situations in which a large number of integrand evaluations are necessary. We direct the reader to [13] for an assessment of computational challenges that is in alignment with ours.

Because they incorporate prior information about the object being modelled, Bayesian methods can be expected to perform better than other methods when the amount of data is limited, a setting in which solving linear systems is not a computational problem. This has resulted in a marked emphasis on applications in which integrand evaluations are expensive [9], [30]. For example, [9] “focus on numerical integrals where the cost of evaluating the integrand forms a computational bottleneck, and for which numerical error is hence more likely to be non-negligible.” Model selection naturally plays all the more important a role in the low-data regime in which non-conjugate Bayesian quadrature methods that would become intractable even for moderately large datasets may shine.

5 Theoretical guarantees↩︎

This chapter contains a number of mathematical results that provide theoretical guarantees on the accuracy of conjugate Bayesian quadrature under certain assumptions on the covariance function \(k\), the integrand \(f\) and the coverage of the domain \(D\) by the nodes. We shall work with a sequence \((X_N)_{N=1}^\infty\) of sets \(X_N = \{ {\boldsymbol{x}}_1, \ldots, {\boldsymbol{x}}_N\} \subset D\) consisting of \(N\) pairwise distinct nodes. In spite of our notation, the node sets need not be nested (i.e., despite both being denoted \({\boldsymbol{x}}_1\), the first nodes in \(X_N\) and \(X_{N+1}\) need not be equal). To further simplify notation, we often replace subscripted \(X_N\) and \({\cal D}\) with \(N\) to emphasise the dependency on the number of nodes. For example, under this convention \(\mu_N\) stands for \(\mu_{\cal D}\) in 8 with \(X = X_N\) and \(\mu_{N,\lambda}\) for its regularized version in 28 . Most theoretical guarantees we give state that \[\label{eq:generic-guarantee} \abs[0]{ I_P(f) - \mu_{N, \lambda} } \leq C( r_N + \sqrt{\lambda} \, )\tag{30}\] for a positive constant \(C\) and a decreasing positive sequence \((r_N)_{N=1}^\infty\) which tends to zero. The rate with which \(r_N\) tends to zero reflects the smoothness of \(k\) and \(f\). A large collection of theoretical guarantees for Bayesian quadrature may be found in references such as [9], [11], [35], [202]–[204]. Proofs are relegated to 7, which also contains a more detailed comparison of our results to those in the literature.

We prioritize simplicity and ease of presentation over generality. Accordingly, throughout this section we assume that the prior mean \(m\) is identically zero, an assumption that could be relaxed to \(m\) being an element of the RKHS of \(k\). Furthermore, in most results of this section we make use of the following assumption on the domain and the integration distribution.

5.1 Guarantees for generic nodes↩︎

We first provide guarantees for generic nodes, by which we mean that no reference is made to the process by which the nodes are obtained and upper bounds of the form 30 are given in terms of the computable fill-distance \(h_X\) and separation radius \(q_X\), which are defined as \[\label{eq:fill-distance} h_{X} = \sup_{ {\boldsymbol{x}}\in [0, 1]^d } \, \min_{i = 1,\ldots, N} \norm[0]{ {\boldsymbol{x}}- {\boldsymbol{x}}_i} \quad \text{ and } \quad q_X= \frac{1}{2} \min_{ i \neq j } \norm[0]{ {\boldsymbol{x}}_i - {\boldsymbol{x}}_j }\tag{31}\] for the domain \(D= [0, 1]^d\). Fill-distance equals the radius of the largest ball in \([0, 1]^d\) that does not contain any of the points in \(X\) while the separation radius is half the smallest distance between any two points in \(X\). Because the fill-distance and separation radius are, at least in principle, computable for any set of nodes, error bounds expressed in terms of these quantities apply to all relevant Bayesian quadrature methods regardless of how the nodes have been sampled (3.3). Choices made in the sampling axis of taxonomy affect the bounds through fill-distance and separation radius.

If, for a sequence of node sets \((X_N)_{N=1}^\infty\), the fill-distance and separation radius are asymptotically comparable as \(N \to \infty\), the point design is said to be quasi-uniform. This means that there exists a constant \(c_\mathrm{qu} > 0\) such that \[q_N \leq h_N \leq c_\mathrm{qu} \, q_N\] for every \(N \geq 1\). This implies that there are positive constants \(c_1\) and \(c_2\) such that \[\label{eq:quasi-uniformity-rate} c_1 N^{-1/d} \leq h_N \leq c_2 N^{-1/d} \quad \text{ and } \quad c_1 N^{-1/d} \leq q_N \leq c_2 N^{-1/d}\tag{32}\] for every \(N \geq 1\) [78]. Regular grids are quasi-uniform. For example, the equispaced nodes \(X_N = \{1/N, 2/N, \ldots, (N-1)/N, 1\}\) on \(D= [0, 1]\) have \(q_N = h_N = 1/N\) for every \(N \geq 1\). How the fill-distance and certain related quantities behave for random points is discussed in [205].

We begin with two generic propositions that guarantee the convergence of Bayesian quadrature under very general conditions. Recall from Section 2.3.1 that \(\mathcal{H}(k)\) is the RKHS of \(k\) on the domain \(D\). The RKHS consists of certain functions defined on \(D\).

Proposition [prop:generic-guarantee-1] guarantees that the integral mean of Bayesian quadrature tends to the true integral if the point design is space-filling in the sense that \(h_N \to 0\). However, the proposition gives us no inkling on the rate of this convergence. The next proposition ensures that at least the familiar square-root rate of Monte Carlo can be attained. That is, conjugate Bayesian quadrature can always beat Monte Carlo.

Propositions [prop:generic-guarantee-1] and [prop:generic-guarantee-2] provide slow rates because they make few assumptions about the integrand, covariance function, and the nodes. The rest of this section shows how smoothness assumptions on \(k\) and \(f\) improve the square-root rate in 33 .

5.1.1 Isotropic Matérn kernels↩︎

The Sobolev space \(H^\alpha([0, 1]^d)\) of order \(\alpha \in \mathbb{N}\) is a Hilbert space that consist of those functions \(g\) for which the weak derivative \(\mathrm{D}^{{\boldsymbol{n}}} g\) exists and is square-integrable for every non-negative multi-index such that \(\abs[0]{{\boldsymbol{n}}} = n_1 + \cdots + n_d \leq \alpha\). That is, \(g \in H^\alpha([0, 1]^d)\) if and only if \[\label{eq:sobolev-norm-weak-derivative} \bigg( \sum_{ \abs[0]{{\boldsymbol{n}}} \leq \alpha} \int_{[0, 1]^d} [ \mathrm{D}^{{\boldsymbol{n}}} g({\boldsymbol{x}}) ]^2 \dif {\boldsymbol{x}}\bigg)^{1/2} < \infty.\tag{34}\] We use the Fourier transform to define the norm \(\lVert \cdot \rVert_{H^\alpha([0,1]^d)}\) of \(H^\alpha([0,1]^d)\); see 41 in 7.2. The resulting norm is equivalent to the norm defined by 34 . Any function whose classical derivatives up to order \(\alpha\) exist and are square-integrable on \([0, 1]^d\) is an element of \(H^\alpha([0, 1]^d)\). Conversely, if \(\alpha > r + d/2\) for \(r \in \mathbb{N}\), then all functions in \(H^\alpha([0, 1]^d)\) are \(r\) times differentiable in the classical sense. In particular, all functions in \(H^\alpha([0, 1]^d)\) are continuous if \(\alpha > d/2\). Sobolev spaces \(H^\alpha([0, 1]^d)\) of fractional smoothness \(\alpha > 0\) may be easily defined via Fourier transforms; details are again given in Section 7.2. Theorems 5.1 and 5.2 provide theoretical guarantees for isotropic Matérn kernels based on the fill-distance. For the purposes of these theorems, define the mesh ratio \[\label{eq:mesh-ratio} \rho_X= \frac{h_X}{q_X} \geq 1.\tag{35}\]

Theorem 5.1 provides an error estimate for any given set of nodes \(X\) in terms of \(h_X\), \(q_X\), and \(\rho_X\). Using 32 to control the fill-distance and separation radius and making sure that constants that depend on kernel parameters do not blow up yields the following variant of Theorem 5.1, which applies to a sequence of node sets \((X_N)_{N=1}^\infty\) produced by a hypothetical Bayesian quadrature method that is run indefinitely without a stopping criterion.

Theorems 5.1 and 5.2 comprise two somewhat different cases. Firstly, if \(\beta = \alpha\), then \(f\) is an element of the RKHS of \(k\) (see Section 2.3.1). In this case the terms \(\smash{\rho_X^{\alpha-\beta}}\) and \(\smash{q_X^{-(\alpha-\beta)}}\) in ?? and \(\smash{N^{(\alpha-\beta)/d}}\) in ?? , which slow down convergence, vanish and the bounds depend on the nodes only via the fill-distance. The theorems do not yield improved rates of convergence if \(f \in H^\gamma([0, 1]^d)\) for \(\gamma > \alpha\) (i.e., when the model undersmooths). Since \(H^\gamma([0, 1]^d) \subset H^\alpha([0, 1]^d)\), in this case one simply applies the theorems with \(\beta = \alpha\).²⁶ Secondly, if \(\beta < \alpha\) and \(f \notin H^\alpha([0, 1]^d)\), then \(f\) is not in \(\mathcal{H}(k)\), so that the model oversmooths. The bound 15 links integration error to the BQ variance but holds only for functions in the RKHS. This gives rise to the distinction between the cases \(f\in \mathcal{H}(k)\) and \(f\notin \mathcal{H}(k)\).

When \(\lambda_N = 0\), Theorem 5.2 gives the rate of convergence \(N^{-\beta/d}\) that is worst-case optimal in the Sobolev space \(H^\beta([0,1]^d)\) [211]. Bayesian quadrature therefore adapts to misspecified smoothness, being capable of attaining the optimal rate of convergence even when \(\alpha > \beta\). This is convenient, as eliciting the smoothness of \(f\) tends to be difficult. As evidenced by the term \(N^{-1/2 + (\alpha-\beta)/d}\), for a positive nugget the rate given by Theorem 5.2 depends on the degree of smoothness misspecification, \(\alpha - \beta \geq 0\). In particular, if \(\alpha \geq \beta + d/2\), Theorem 5.2 does not guarantee convergence unless \(\lambda_N\) tends to zero sufficiently fast. This may be due to sub-optimality of some estimates the proof uses.

For a machine learning practitioner or a statistician, randomly sampled nodes are undoubtedly the most relevant ones. The following theorem, which has not appeared in the context of Bayesian quadrature before, guarantees that for Matérns the expected error decays with the worst-case optimal rate.

5.1.2 Product Matérn kernels↩︎

For a product Matérn kernel in 5 , theoretical guarantees can be derived from the results in [212] if the nodes form a sparse grid on \(D= [0, 1]^d\) and the integrand is in a mixed-order Sobolev space [213]–[215]. For \(\alpha \in \mathbb{N}\), the mixed-order Sobolev space \(H^{\alpha}_\mathrm{mix}([0, 1]^d)\) is the tensor product of univariate Sobolev spaces \(H^\alpha([0,1])\) and consists of \(g \in L^2([0, 1]^d)\) whose weak derivative \(\mathrm{D}^{{\boldsymbol{n}}} g\) exists and is in \(L^2([0,1]^d)\) for every multi-index \({\boldsymbol{n}}= (n_1, \ldots, n_d)\) such that \(n_i \leq \alpha\) for every \(i\). If \(k\) is a product Matérn kernel 5 with smoothness \(\nu > 0\) such that \(\alpha = \nu + 1/2\) is a positive integer \(f \in H^{\beta}_\mathrm{mix}([0, 1]^d)\) for an integer \(\beta\) such that \(1 \leq \beta \leq \alpha\), it follows from [212] that there is \(C > 0\), which does not depend on \(f\) or \(N\), such that \[\abs[0]{ I_P(f) - \mu_{N} } \leq C N^{-\beta} (\log N)^{(\beta + 1)(d-1)} \norm[0]{f}_{H^\beta_\mathrm{mix}([0,1]^d)}\] when the nodes form a certain sparse grid (see Figure 9 for an example). Due to the requirement that the nodes form a sparse grid, this result is much more restrictive than Theorems 5.1 to 5.3. We do not go into further details so as to avoid having to review the construction sparse grids. In general, error analysis for product Matérns is much less mature than for isotropic Matérns.

5.1.3 Square exponential kernel↩︎

Here we provide theoretical guarantees for conjugate Bayesian quadrature based on the square exponential kernel in 3 . Unlike those for Matérn kernels, these guarantees apply only to functions in the RKHS of the kernel, \(\mathcal{H}(k)\). One can show [216] that the RKHS on any \(D\subseteq \mathbb{R}^d\) with non-empty interior consists precisely of the functions \[g({\boldsymbol{x}}) = \exp\Big( \! - \frac{\norm[0]{{\boldsymbol{x}}}^2}{2 \ell^2} \Big) \sum_{ {\boldsymbol{n}}\in \mathbb{N}_0^d } c_{{\boldsymbol{n}}} {\boldsymbol{x}}^{\boldsymbol{n}}\quad \text{ such that } \quad \sum_{{\boldsymbol{n}}\in \mathbb{N}_0^d} \ell^{2\abs[0]{{\boldsymbol{n}}}} {\boldsymbol{n}}! c_{\boldsymbol{n}}^2 < \infty,\] where the sums are over all non-negative \(d\)-dimensional multi-indices \({\boldsymbol{n}}= (n_1, \ldots, n_d)\) and we denote \({\boldsymbol{x}}^{\boldsymbol{n}}= x_1^{n_1} \cdots x_d^{n_d}\) and \({\boldsymbol{n}}! = n_1 ! \cdots n_d!\). The RKHS is minuscule: it consists of infinitely differentiable functions but contains no polynomials except the zero function. The current state of the theory is inadequate to provide theoretical guarantees for \(f\notin \mathcal{H}(k)\).

For Matérns, the length-scale \(\ell\) only affects the norm of the RKHS, not which functions are contained therein; see 48 and 49 . However, for the square exponential kernel \(\ell\) does affect the size of the RKHS and thus plays a role similar to the smoothness parameter \(\nu\) in Matérn kernels. We therefore require that \(\ell\) be kept fixed in the results below.

Like Theorem 5.1, the above theorem applies to any given set of nodes. By instead considering a quasi-uniform sequence of node sets generated by some sampler we obtain an explicit rate of convergence given by the following theorem.

Note that the terms dependent on the nugget in Theorems 5.4 and 5.5 blow up very fast as \(h_X \to 0\) or \(N \to \infty\). These theorems are useful only when the nugget is close to zero (or tends to zero sufficiently fast). When \(\lambda_N = 0\), the convergence in Theorem 5.5 is super-exponential and hence much faster than in Theorems 5.2 to 5.3, which yield polynomial rates for Matérns. This is because the class of functions these theorems apply to, the RKHS of the square exponential kernel, is significantly smaller than any Sobolev space. Theorems 5.4 and 5.5 do not therefore demonstrate that the square exponential kernel is in any way superior to a Matérn (or any other) kernel unless one is certain that \(f\) is in \(\mathcal{H}(k)\). We note that on dimension one and for product designs in higher dimensions it is possible to obtain slightly improved rates with explicit constants \(c_1\) and \(c_2\) by using results in [217], [218].

5.2 Guarantees for sequential Bayesian quadrature↩︎

Consider a sequential setting where the nodes are obtained by sequential maximization of an acquisition function that is a strictly increasing transformation \(\psi \colon [0, \infty) \to \mathbb{R}\) of the scalar squared correction in 21 . That is, for each \(n \geq 1\), the \((n+1)\)th node \({\boldsymbol{x}}_{n+1}\) satisfies \[\label{eq:sequential-selection-for-guarantees} \begin{align} {\boldsymbol{x}}_{n+1} &\in \mathop{\mathrm{arg\,max}}_{\tilde{{\boldsymbol{x}}}} \psi(\rho^2(\tilde{{\boldsymbol{x}}})) \\ &= \mathop{\mathrm{arg\,max}}_{\tilde{{\boldsymbol{x}}}} \psi\big( {\boldsymbol{k}}_{PX(\tilde{{\boldsymbol{x}}})}^{\intercal}{\boldsymbol{K}}_{X(\tilde{{\boldsymbol{x}}})X(\tilde{{\boldsymbol{x}}})}^{-1}{\boldsymbol{k}}_{PX(\tilde{{\boldsymbol{x}}})} \big), \end{align}\tag{36}\] where \(X(\tilde{{\boldsymbol{x}}}) = X \cup \{ \tilde{{\boldsymbol{x}}} \} = \{{\boldsymbol{x}}_1, \ldots, {\boldsymbol{x}}_n, \tilde{{\boldsymbol{x}}}\}\). Examples include the MI, IVR, and NIV acquisition functions from 4.2.3. The following theorem demonstrates that sequential conjugate Bayesian quadrature converges at least as fast as Monte Carlo.

The rate ?? is likely conservative [12]. See [203], [219]–[221] for results on the convergence of Bayesian quadrature and related methods with acquisition functions other than 36 .

5.3 Other guarantees↩︎

Various other guarantees for Bayesian quadrature are available in the literature. These include Bayesian quadrature with (i) nodes having strong symmetry properties [94]; (ii) nodes obtained by optimization [203], [222], [223]; (iii) nodes sampled from distributions different from \(P\) ([93]; see also [97]); (iv) nodes sampled from a determinantal point process [102], [104], [105]; and (v) nodes selected using a certain recombination algorithm [135]. [224] have proved error bounds for BQ applied to the computation of normalizing constants of the form \(\int_D\exp(-\lambda f({\boldsymbol{x}})) \dif {\boldsymbol{x}}\) with \(\lambda > 0\). Convergence guarantees for multi-output and multilevel Bayesian quadrature are provided in [14] and [20], respectively. [225] consider subsampling kernel mean embeddings approximated with Monte Carlo. The equivalences reviewed in 2.3.1 allow one to apply the wealth of results collected in [58] to Bayesian quadrature.

6 Empirical illustrations↩︎

As mentioned in the introduction, the aim of our experiments is not to re-produce the performance of Bayesian quadrature methods, this having been done in the literature we have cited, nor to propose and evaluate new methods. Rather, the aim is to illustrate the taxonomy of 3 and the qualitative influence of its axes, sampling, model and inference. In particular, the experiments show ablation curves that disentangle the effect of the axes for several integrands \(f\) that differ in their degree of smoothness or (non-)stationarity. Where possible, we tie the qualitative or quantitative behavior of the BQ methods to the theoretical results from 5. We hope the illustrations provide additional intuition on BQ methods for practitioners.

The available open source code for BQ is unfortunately still scarce in comparison to the diversity of methods proposed in the literature. Under this constraint, we have opted to implement the experiments in the open source Python library EmuKit [130], [131] since, at the time of writing, it provided the highest diversity of BQ algorithms. It further enabled us to code in a single language. In particular, EmuKit contains a variety of kernel embeddings for the prevalent conjugate BQ model and a variety of node selection schemes. Another contender, the Python library ProbNum [133], did not contain adequate hyper-parameter tuning capability, nor state-of-the-art active sampling loops at the time of writing.²⁷ Unfortunately, scalable BQ methods are currently not available in EmuKit and are thus omitted here; some are implemented in the MATLAB library GAIL [132]. All experiments use exact inference. Section 6.1 describes the experimental setup and integrands (called test functions) in more detail. Section 6.2 illustrates and discusses results.

6.1 Experimental setup↩︎

This section describes the experimental setup.

6.1.1 Kernels (model axis)↩︎

As all the experiments use the conjugate model, the models only differ by their choice of the kernel \(k\). We select 5 different kernels to denote 5 different models. Three of the kernels are Matérn kernels in 4 with smoothness parameter \(\nu \in \{\frac{1}{2}, \frac{3}{2}, \frac{5}{2}\}\). As described in Sections 2.3.1 and 5.1.1, these kernels are linked to Sobolev spaces. Matérn kernels are stationary. Further, we use the non-stationary and non-smooth Brownian motion kernel as in 6 and the stationary and infinitely smooth square exponential kernel in 3 . All selected kernels have analytic kernel embeddings with respect to \(P\) that are available in EmuKit. Table ¿tbl:tab:exp-hypers? summarizes the kernels and their parameters.

6.1.2 Node selection (sampling axis)↩︎

Along the sampling axis we use 5 sampling schemes: two random ones and three deterministic ones to cover the mode of sampling. The two random schemes are random sampling from \(P\) and the space-filling Latin hypercube sampling (Section 4.2.2). Two of the three deterministic samplers produce model-independent designs: the space-filling Sobol sequence and the classically motivated Legendre roots (Section 4.2.1). The final deterministic node selection scheme uses a sequential, active sampler and the IVR acquisition function (Section 4.2.3). The samplers are summarized in Table ¿tbl:tab:exp-sampler?. Further details can be found in Section 6.3.

6.1.3 Implementation details↩︎

As the experiments use the conjugate model throughout, \(\mathsf{f}=\mathsf{g}\) is a Gaussian process. To represent \(\mathsf{f}\) we use the GPy Python library [226] and for the kernel means and variances \(k_{P}\) and \(k_{PP}\) we use the Python library EmuKit [130], [131] and its existing GPy quadrature wrappers. The kernel parameters \({\boldsymbol{\theta}}\) (Table ¿tbl:tab:exp-hypers?) are fitted with maximum marginal likelihood (Section 4.3), which GPy defaults to.

We use the default settings for the nugget of both EmuKit and GPy. The nugget \(\lambda\) (Section 4.4) is composed of three terms: \[\lambda = \lambda_{\mathrm{Emu}} + \lambda_{\mathrm{GPy}} + \lambda_{\mathrm{GPyDyn}} .\] The first two terms, \(\lambda_{\mathrm{Emu}}=10^{-10}\) and \(\lambda_{\mathrm{GPy}}=10^{-8}\), are the fixed default nuggets of EmuKit and GPy respectively. While \(\lambda_{\mathrm{Emu}}\) can be adjusted via the EmuKit interface, \(\lambda_{\mathrm{GPy}}\) is not directly accessible to the user. The dynamic nugget term \(\lambda_{\mathrm{GPyDyn}}\) is computed as follows. Let \({\boldsymbol{A}}= {\boldsymbol{K}}_{XX} + (\lambda_{\mathrm{Emu}} + \lambda_{\mathrm{GPy}} ){\boldsymbol{I}}_N\) be the kernel Gram matrix augmented with the fixed nugget terms and let \(\bar{a}=\frac{1}{N}\sum_{i=1}^NA_{ii}\) be the mean of its diagonal elements. The third summand, \(\lambda_{\mathrm{GPyDyn}}\), is a dynamic nugget that GPy adds automatically when the Cholesky decomposition of \({\boldsymbol{A}}\) (used to solve the linear system of the kernel Gram matrix) fails. \(\lambda_{\mathrm{GPyDyn}}\) is initialized to zero and then increased up to 5 times, taking the successive values \(0\), \(\bar{a}10^{-6}\) , \(\bar{a}10^{-5}\) , \(\bar{a}10^{-4}\), \(\bar{a}10^{-3}\), and \(\bar{a}10^{-2}\) until the Cholesky decomposition is successful. If the Cholesky decomposition fails with the largest dynamic nugget \(\bar{a}10^{-2}\), the method finally gives up and throws an error. In contrast to the fixed nugget terms, the dynamic nugget depends on \({\boldsymbol{K}}_{XX}\) and is reset to zero after each successful Cholesky decomposition. Hence every Gram matrix in the experiments incorporates a nugget that is at least \(10^{-8}\) but may vary from one matrix to the next.

All nodes, except the ones obtained actively via IVR, can be pre-computed as their sampling is model-independent. In the settings where we pre-compute nodes we (i) first sample the nodes, (ii) then hand them to the EmuKit BQ wrapper as a fixed dataset, (iii) fit the hyper-parameters \({\boldsymbol{\theta}}\), and (iv) finally compute \(I_{P}(\mathsf{f})\mid{\cal D}\) using EmuKit’s kernel embeddings. When the nodes are random, we use the same NumPy seed of 0 across all experiments. This means that random, model-independent nodes are identical across models and integrands \(f\) for a given \(N\), which may be reflected in some curves showing similar fluctuations. We elaborate on this further when describing the results below. For active node selection with IVR, we fit the hyper-parameters successively after every newly obtained node. IVR as an active node acquisition method further requires an initial design which we choose to be \(N_{\mathrm{init}}=\min\{N, 5\}\) random nodes from \(P\). This ensures that hyper-parameter optimization by maximum marginal likelihood does not fail. The initial design uses the same NumPy seed of 0.

6.1.4 Integrands (test functions)↩︎

An integration problem is composed of an integrand or, test function, \(f\) and a measure \(P\). For the experiments we choose \(P\) to be the standard Lebesgue (i.e., uniform) measure on the unit interval \(D=[0, 1]\). The test functions \(f\) cover a smooth and non-smooth as well as stationary and non-stationary²⁸ integrands. All test functions are continuous and do not exhibit “jumps”.

In order to create distinct test functions with similar properties, we choose parameterized families of functions and sample the parameters for an individual test function from some distribution. That is, \(f\) is deterministic for each fixed set of parameters. In total, for every parameterized family we produce \(T=50\) deterministic integrands for evaluation. The two families are:

• Univariate, smooth, stationary integrands: We use smooth test functions as in [227]. The integrand is defined as a finite Fourier series \[\label{eq:fourier-test-func} f(x) = \sqrt{2}\sum_{j=0}^J{[a_j\cos(2\pi j x/L) + u_j\sin(2\pi j x/L)]}\tag{37}\] with length-scale \(L=5\) and \(J=25\) terms and scalar coefficients \(a_j\) and \(u_j\) for \(j=0,\dots, J\). The chosen values for \(L\) and \(J\) ensure that the function is reasonably variable in the domain \(D\) yet not too easy to learn. All functions are infinitely smooth. Each of the \(T=50\) sets of coefficients \(\{a_j, u_j\}_{j=1}^J\) is produced according to \(a_j, u_j\sim\mathrm{N}(0, 2(J+1))\). We set \(a_0=u_0 = 0\) to define an offset.

• Univariate, non-smooth, continuous, non-stationary integrands: These are approximate Brownian motion paths represented by a truncated Karhunen–Loève expansion of the Brownian motion: \[\label{eq:BM-test-func} f(x) = \sqrt{2}\sum_{j=1}^J a_j \frac{\operatorname{sin}((j-\frac{1}{2})\pi x)}{(j-\frac{1}{2})\pi} .\tag{38}\] We choose \(J=500\) terms. Each of the \(T=50\) sets of sets of coefficients \(\{a_j\}_{j=1}^J\) is sampled according to \(a_j\sim\mathrm{N}(0, 1)\). The constructed functions resemble non-smooth, continuous functions quite well.

All of the above integrands can be represented as function handle \(f(\cdot)\) in code as opposed to some fixed vector of observations \({\boldsymbol{f}}_{X}\). This is achieved by constructing \(f(\cdot)\) as a finite sum of basis functions. In our experiments, an explicit function handle is only really needed by active IVR sampling; the BQ model and inference do not require a function handle per se.

The ground truth integrals are compute with SciPy’s integrate.quad method to high precision. Table ¿tbl:tab:exp-f? summarizes the tests functions. Illustrations of some of the \(T=50\) integrands for each family are shown in Figure 11.

6.1.5 Description of individual experiments↩︎

We pair each of the five conjugate BQ models of Section 6.1.1 with the five node selection schemes of Section 6.1.2. This yields a total of 25 BQ methods, which we identify by the kernel name and sampling method. Then we evaluate each of the methods on the \(T=50\) test integrands of each family of integrands of Section 6.1.4, which amounts to 1250 combinations of test integrands and BQ methods per family.

Furthermore, we evaluate each of those 1250 combinations on several fixed budgets of size \(N\) ranging from \(N=1\) to some maximal \(N_{f}\) in increments of 1. The maximal budget \(N_f\) depends on the integrand family. We choose \(N_f = 30\) and \(N_f = 50\) for the smooth and non-smooth univariate problems, respectively. These numbers are chosen to illustrate the performance change of the BQ methods across a relevant range of \(N\). In total, we execute \(100\,000\) BQ runs, which were run sequentially on a general purpose m1.xlarge instance with 4 vCPUs.

6.1.6 Scores↩︎

Recall from Eqs. 8 and 10 that \(\mu_{{\cal D}}\) and \(\Sigma_{{\cal D}}\) stand for the posterior mean and variance, respectively, of the conjugate BQ model. The mean \(\mu_{{\cal D}}\) serves as an estimator for the true integral value \(I_{P}(f)\) and \(\Sigma_{{\cal D}}\) describes the spread of the posterior distribution \(I_P(\mathsf{f}) \mid {\cal D}\).

In the experiments, we evaluate the behavior of the BQ methods with two empirical scores. The first score is the error score \[\frac{1}{T}\sum_{t=1}^T\left|\frac{I_{P}(f_t) - \mu_{{\cal D}, t}}{I_{P}(f_t)}\right|.\] The error score quantifies the performance of \(\mu_{{\cal D}}\) as an estimator for the ground truth \(I_{P}(f)\). More precisely, the error score is the average relative error between \(\mu_{{\cal D}}\) and \(I_{P}(f)\) over the test functions. The subscript \(t\) indicates the \(t\)^th test integrand; the average is computed over the \(T=50\) integrands for each family separately.²⁹

The second score we use is the calibration score \[\frac{1}{T}\sum_{t=1}^T\frac{|I_{P}(f_t) - \mu_{{\cal D}, t}|}{\Sigma_{{\cal D}, t}^{1/2}}.\] The calibration score quantifies the calibration of the BQ method, or how well the credible interval of the posterior describes the error between the estimator \(\mu_{{\cal D}}\) and the true integral value \(I_{P}(f)\). The calibration score equals the average absolute error in units of the corresponding credible intervals. For each test function, the 68% credible interval centered around \(I_{P}(\mathsf{f}) \mid {\cal D}= I_{P}(f)\) and 1 standard deviation wide to each side is , or equivalently for the shifted value \(I_{P}(\mathsf{f}) \mid {\cal D}- \mu_{{\cal D}} = I_{P}(f) - \mu_{{\cal D}}\). Taking the absolute value lets us compare the size of \(|I_{P}(f) - \mu_{{\cal D}}|\) (enumerator) to \(\Sigma_{{\cal D}}^{1/2}\) (denominator). For an individual test function \(f\), the score is smaller than 1 if the true value \(I_{P}(f)\) lies inside the credible interval and larger than 1 if it lies outside. Hence, a somewhat constant calibration score across different \(N\) (e.g., hovering around \(1\)) suggests that the BQ method is well calibrated. An increasing or decaying score suggests that the method is prone to over- or underconfidence, respectively [204], [228]. It is usually not as harmful for a method to be underconfident (too uncertain) than overconfident (too certain).³⁰

6.2 Results↩︎

This section describes the qualitative behavior of the BQ methods on the families of test functions.

6.2.1 Smooth test functions (sampling axis)↩︎

Figure [fig:sampling-config003] depicts the performance of the 5 models along the sampling axis for the family of univariate smooth test functions. Each row shows a single model with kernels of increasing smoothness from top to bottom. The first 4 rows show models of finite, and hence of lesser, smoothness than the test functions (Brownian motion, Matérn \(\frac{1}{2}\), Matérn \(\frac{3}{2}\), and Matérn \(\frac{5}{2}\)), while the last row shows a model with matching (infinite) smoothness (square exponential). The left and right columns show the error score and the calibration score, respectively. From the plots we can make several qualitative observations.

First, it is apparent that the sampling choice has a major effect on the error score (left column). In this particular example, Legendre roots seem to outperform all other sampling methods considered, followed by IVR. Random samples from \(P\) perform the worst across all models. The Sobol sequence and Latin hypercube perform similarly, and generally better than random samples from \(P\). While not covered by the theory in Section 5, the superiority of Legendre roots is likely related to the well-known phenomenon in numerical analysis that nodes which are denser near the boundary of the domain than in its interior are good for approximating and integrating infinitely smooth functions [177], [230]. In Figure [fig:sampling-config017] we observe that the performance of Legendre nodes differs little from other sampling choices when the integrands are non-smooth.

Second, from the right column we see that all models are well calibrated or underconfident for , as the calibration scores hover around a constant value or tend to zero. There is a “kink” in the calibration score for random nodes (red) at \(N \approx 18\) for nodes sampled from \(P\) that vanishes at \(N=20\) (coinciding with a drop of the corresponding error scores). This is an artifact of a “gap” in the randomly sampled nodes that gets filled at \(N = 20\). The other, space-filling designs do not exhibit this fluctuation. The “kink” is less pronounced for more well-specified (i.e., smoother) models. The “kinks” do not average out across the \(T=50\) integrands because the random nodes \(X\) are identical for a given \(N\) and across all experiments as the same random seed was used (see Section 6.1.3). When the model does not oversmooth a univariate stationary integrand, the sampling choice therefore appears to mostly affect the performance of the mean estimator \(\mu_{{\cal D}}\). Calibration is decent and largely unaffected by sampling. This makes the BQ methods usable in this particular setting.

The third observation we make is the flattening of the performance score in the bottom left plot of Figure [fig:sampling-config003] for the Legendre nodes at around \(10^{-7}\). This is likely due to the positive nugget used in the experiments (see Section 6.1.3). [thm:matern-generic,thm:matern-quasi-uniform,thm:se-generic,thm:se-quasi-uniform] suggest that a positive nugget may either greatly slow the rate of convergence or altogether prevent convergence after a certain accuracy has been reached.

6.2.2 Smooth test functions (model axis)↩︎

Figure [fig:kernel-config003] shows the same curves as Figure [fig:sampling-config003] but this time along the model axis. That is, rows are now grouped according to the sampling choice such that the effect of the kernel can be observed in each subplot. The nodes \(X\) are identical on each row and \(N\) across all test functions if sampling is model-independent (rows 1–4). For deterministic nodes this is fulfilled trivially; for random nodes this is because the same random seed was used across experiments (see Section 6.1.3). Hence, solely the effect of the model can be observed in each plot. The nodes on row 5 differ across \(N\), as well as models and test functions, as the IVR sampler is active and model-dependent. Again we make some qualitative observations.

First, the right column shows that all models are well calibrated or underconfident. The “kink” mentioned in Section 6.2.1 at \(N \approx 18\) for random nodes from \(P\) can be observed in the top right plot. There are no substantial differences between the models, except that the infinitely smooth square exponential model is generally the least underconfident. Because the test integrands are infinitely smooth, this is not surprising.

Second, the error score (left column) decreases as the kernel becomes smoother. In particular, the infinitely smooth square exponential kernel performs better than other kernels across all sampling choices. This is consistent with the theoretical guarantees of Section 5, which indicate that smoother kernels can attain faster rates of convergence provided that the integrand is sufficiently smooth. In conclusion, models that undersmooth the truth are decently calibrated but have larger integration errors than models that do not undersmooth.

6.2.3 Non-smooth test functions (sampling axis)↩︎

Figure [fig:sampling-config017] depicts the performance of the 5 models along the sampling axis for the family of approximate Brownian motion paths. Again, each row shows a single model with kernels of increasing smoothness from top to bottom. The first two rows show models that are continuous but non-differentiable and match the smoothness (or roughness) of the test functions. The models on rows 3–5 are increasingly smooth and thus oversmooth the truth. In particular, the square exponential kernel on row 5 is “infinitely smoother” than the test functions. Again, the left and right columns show the error score and the calibration score, respectively. From the plots we can make the following qualitative observations.

First, the sampling choice has no apparent effect on the error score for the models with matching smoothness (rows 1–2); solely random sampling from \(P\) performs slightly poorer. This is quite surprising and unlike the behavior we observed in Section 6.2.1 for the family of smooth test functions. As expected, calibration is good for all sampling schemes.

Second, models that oversmooth (rows 3–5) perform successively worse in error score. Interestingly, random sampling from \(P\) reinforces the adverse effect of oversmoothing. This may be explained by the presence of the mesh ratio term \(\rho_X^{\alpha - \beta}\) [mesh ratio is define in Eq. 35 ] in 5.1: when there is significant oversmoothing (\(\alpha - \beta\) is large), this term is large. Moreover, the mesh ratio is large if the nodes exhibit clustering, which explains the behavior of random sampling. The calibration score is less stable and indicates increasing overconfidence across sampling schemes for models that oversmooth.

Unlike undersmoothing, oversmoothing the truth can have a detrimental effect on error and calibration scores across all sampling choices. Oversmoothing the truth slightly can have an undesirable effect on the calibration score while the error score is mostly unaffected unless random sampling from \(P\) is used. At the same time, the choice of sampler is less important for models with matching smoothness. Random sampling seems to be the poorest choice for all models, a fact that may be explained by the observation above on the magnitude of the mesh ratio for random nodes.

6.2.4 Non-smooth test functions (model axis)↩︎

Figure [fig:kernel-config017] shows the same curves as Figure [fig:sampling-config017] but this time along the model axis. That is, rows are now grouped according to the sampling choice such that the effect of the kernel can be observed in each subplot. Again, the nodes \(X\) are identical on each row and \(N\) across all test functions if sampling is model-independent (rows 1–4). Hence, solely the effect of the model can be observed in each plot. Again, the nodes in the bottom row 5 differ across \(N\), as well as models and test functions, as the IVR sampler is active and model-dependent. We observe the following qualitative behavior.

Unsurprisingly, models that match the roughness of the integrand (Brownian and Matérn \(\frac{1}{2}\)) perform best on both error and calibration scores for all sampling schemes. Models that oversmooth slightly (Matérn \(\frac{3}{2}\) and Matérn \(\frac{5}{2}\)) have equally good error scores (exception described below) but appear to be somewhat overconfident. The square exponential kernel that oversmooths heavily performs poorly across all sampling schemes. Random sampling from \(P\) (row 1) is the only sampling scheme where even slight oversmoothing has an adverse effect on the the error score. This effect was given theoretical explanation in 6.2.3. Random sampling should therefore be avoided whenever there is a risk of oversmoothing. However, random sampling can be expected to perform well and be well-calibrated if smoothness is correctly specified.

6.3 Further experimental details↩︎

This section contains further details on the implementation of the experiments. The samplers of Table ¿tbl:tab:exp-sampler? were implemented as follows:

• Random from \(P\) (random): Nodes are random samples from the integration measure \(P\). We use the method \[\text{\texttt{uniform} \: of \: {NumPy}'s \texttt{random.default\_rng}}\] when \(P\) is the uniform measure on the unit interval \([0, 1]\).

• Latin hypercube (random): Nodes are samples of a Latin hypercube [178], a distribution of space-filling pseudo hyper-grids. We use the method \[\text{ \texttt{sample} \: of \: {SciPy}'s \texttt{stats.qmc.LatinHypercube} }\] with default arguments.

• Sobol sequence (deterministic): Nodes are the Sobol sequence in dimension \(d\). We use the method \[\text{\texttt{random\_base2} \: of \: {SciPy}'s \texttt{stats.qmc.Sobol} }\] with scramble=False. A Sobol sequence produces points for \(N\), an integer power of 2. For \(N\) which are not powers of 2, we use the next larger sequence and truncate it.

• Legendre roots (deterministic): Nodes are the roots of some Legendre polynomial. We use the method \[\text{\texttt{roots\_legendre} \: of \: {SciPy}'s \texttt{special}}\] module.

• Integral variance reduction (sequential, active, model-dependent): Nodes are collected sequentially and myopically with batch size \(n=1\). The acquisition function \(a(\tilde{{\boldsymbol{x}}}\mid P)\) is the integral variance reduction (IVR) acquisition function from Section 4.2.3. We use \[\begin{align} &{\tt IntegralVarianceReduction} \\ &\text{of \: {EmuKit}'s {\tt quadrature.acquisitions}} \end{align}\] module as well as its sequential acquisition loop.

7 Proofs for Chapter 5↩︎

This chapter contains proofs for the results in 5. We begin with some preliminaries.

7.1 Reproducing kernel Hilbert spaces↩︎

Some basic properties of RKHSs have appeared already in 2.3. A kernel \(k\colon D\times D\to \mathbb{R}\) is positive-semidefinite if the kernel Gram matrix \({\boldsymbol{K}}_{XX} = (k({\boldsymbol{x}}_i, {\boldsymbol{x}}_j))_{i,j=1}^N\) is positive-semidefinite for any \(N \in \mathbb{N}\) and any set \(X= \{{\boldsymbol{x}}_1, \ldots, {\boldsymbol{x}}_N \} \subseteq D\) of points. The classical Moore–Aronszajn theorem [231] states that for every positive-semidefinite kernel there is a unique Hilbert space \(\mathcal{H}(k)\), the RKHS of \(k\), which consists of real-valued functions defined on \(D\) and is equipped with an inner product \(\langle \cdot, \cdot \rangle_{\mathcal{H}(k)}\) and the corresponding norm \(\norm{\cdot}_{\mathcal{H}(k)}\) such that the kernel has the following reproducing property: \[\label{eq:reproducing-property} \langle g, k(\cdot , {\boldsymbol{x}}) \rangle_{\mathcal{H}(k)} = g({\boldsymbol{x}}) \quad \text{ for all } \quad g \in \mathcal{H}(k) \text{ and } {\boldsymbol{x}}\in D.\tag{39}\] Whenever it is necessary to be precise about the domain of the functions in the RKHS, we say that \(\mathcal{H}(D)\) is an RKHS on \(D\). From the Riesz representation theorem it follows that the reproducing property extends to any continuous linear functional on \(\mathcal{H}(k)\). Namely, if \(L \colon \mathcal{H}(k) \to \mathbb{R}\) is a continuous linear functional, then \[L(g) = \langle g, k_L \rangle_{\mathcal{H}(k)} \quad \text{ for all } \quad g \in \mathcal{H}(k),\] where \(k_L\), the representer of \(L\), is defined as \(k_L({\boldsymbol{x}}) = L(k(\cdot, x))\). Selecting \(L = I_P\) yields \[I_P(g) = \langle g, k_P\rangle_{\mathcal{H}(k)},\] where \(k_P({\boldsymbol{x}}) = \int_Dk({\boldsymbol{x}}', {\boldsymbol{x}}) \dif P({\boldsymbol{x}}')\) is the kernel mean embedding in 9 .

It is usually not straightforward to decide whether or not a given function is an element of the RKHS of a given kernel. There are two characterisations of the RKHS that are occasionally helpful. If there are functions \(\{\varphi_i\}_{i=1}^\infty\) such that \[k({\boldsymbol{x}}, {\boldsymbol{x}}') = \sum_{i=1}^\infty \varphi_i({\boldsymbol{x}}) \varphi_i({\boldsymbol{x}}') \quad \text{ for all } \quad {\boldsymbol{x}}, {\boldsymbol{x}}' \in D,\] then every \(g \in \mathcal{H}(k)\) may be written as \(g = \sum_{i=1}^\infty c_i \varphi_i\) (the convergence is both in \(\norm{\cdot}_{\mathcal{H}(k)}\) and pointwise) for some \(c_i \in \mathbb{R}\) such that \(\sum_{i=1}^\infty c_i^2 < \infty\) [56]. Conversely, if no such expansion for square-summable \(c_i\) exists, then \(g \notin \mathcal{H}(k)\).

If the kernel is stationary on \(\mathbb{R}^d\), in that there is a function \(\Phi_k\colon \mathbb{R}^d \to \mathbb{R}\) such that \(k({\boldsymbol{x}}, {\boldsymbol{x}}') = \Phi_k({\boldsymbol{x}}- {\boldsymbol{x}}')\) for all \({\boldsymbol{x}}, {\boldsymbol{x}}' \in \mathbb{R}^d\), then the RKHS of \(k\) on \(\mathbb{R}^d\) may be expresssed in terms of Fourier transforms ([53]; [78] Thm. 10.12). Namely, if \(\Phi_k\) is continuous and integrable, then \(\mathcal{H}(k)\) is \[\label{eq:RKHS-Fourier-Rd} \Set[\bigg]{ g \in C(\mathbb{R}^d) \cap L^2(\mathbb{R}^d) }{ \norm{g}_{\mathcal{H}(k)}^2 = \frac{1}{(2\pi)^{d/2}} \int_{\mathbb{R}^d} \frac{\abs[0]{\widehat{g}({\boldsymbol{\xi}})}^2}{\widehat{\Phi}_k({\boldsymbol{\xi}})} \dif {\boldsymbol{\xi}}< \infty},\tag{40}\] where \(\widehat{g}({\boldsymbol{\xi}}) = (2\pi)^{-d/2} \int_{\mathbb{R}^d} g({\boldsymbol{x}}) e^{-\mathrm{i} {\boldsymbol{\xi}}^{\intercal}{\boldsymbol{x}}} \dif {\boldsymbol{x}}\) denotes the Fourier transform. To use the Fourier characterisation on a subset of \(\mathbb{R}^d\) we use a standard result on restrictions of RKHSs [25]. Let \(\mathcal{H}(k, \mathbb{R}^d)\) denote the RKHS in 40 and let \(\mathcal{H}(k)\) be the RKHS on \(D\subset \mathbb{R}^d\). Then \(\mathcal{H}(k)\) consists of functions \(g \colon D\to \mathbb{R}\) for which there exists an extension \(g_e \in \mathcal{H}(k, \mathbb{R}^d)\) such that \(g_e({\boldsymbol{x}}) = g({\boldsymbol{x}})\) for all \({\boldsymbol{x}}\in D\) (i.e., \(g_e|_D= g\)). Moreover, \[\norm[0]{g}_{\mathcal{H}(k)} = \min\Set{ \norm[0]{g_e}_{\mathcal{H}(k, \mathbb{R}^d)} }{ g_e \in \mathcal{H}(k, \mathbb{R}^d) \text{ is s.t. } g_e|_D= g}.\]

7.2 Sobolev spaces↩︎

In Section 5.1.1, we defined the Sobolev space \(H^\alpha([0, 1]^d)\) of integer smoothness \(\alpha\) as the space of functions whose weak derivatives up to order \(\alpha\) exist and are in \(L^2([0, 1]^d)\). To allow fractional smoothness we may turn to Fourier transforms. For any \(\alpha > 0\), we let \(H^{\alpha}(\mathbb{R}^d)\) be the Hilbert space of functions \(g \in L^2(\mathbb{R}^d)\) such that \[\label{eq:sobolev-norm-fourier} \norm[0]{ g }_{H^{\alpha}(\mathbb{R}^d)}^2 = \int_{\mathbb{R}^d} (1 + \norm[0]{{\boldsymbol{\xi}}}^2)^\alpha \abs[0]{ \widehat{g}({\boldsymbol{\xi}}) }^2 \dif {\boldsymbol{\xi}}< \infty.\tag{41}\] To illustrate the relationship of this norm to condition in 34 expressed using weak derivatives, consider \(\alpha \in \mathbb{N}\) and \(d = 1\). We may then use the binomial theorem, the properties of the Fourier transform, and Parseval’s identity to compute \[\begin{align} \norm[0]{ g }_{H^{\alpha}(\mathbb{R})}^2 = \int_{\mathbb{R}} (1 + \xi^2)^\alpha \abs[0]{ \widehat{g}(\xi) }^2 \dif \xi &= \sum_{n=0}^\alpha \binom{\alpha}{n} \int_{\mathbb{R}} \abs[0]{ \xi^{n} \widehat{g}(\xi) }^2 \dif \xi \\ &= \sum_{n=0}^\alpha \binom{\alpha}{n} \int_{\mathbb{R}} [ \mathrm{D}^n g(x) ]^2 \dif x, \end{align}\] which is, up to binomial coefficient factors for each integral, equal to the left-hand side of 34 defined on \(\mathbb{R}\) instead of on \([0, 1]\). A fractional Sobolev space and its norm on \([0, 1]^d\) can be defined via restrictions in a way analogous to how we related \(\mathcal{H}(k)\) to \(\mathcal{H}(k, \mathbb{R}^d)\) in Section 7.1. One can then show that \(g\) is in the resulting space if and only if it satisfies 34 .

7.3 Proofs for Section 5.1↩︎

This section contains proofs for the results in Section 5.1.

7.3.1 Proofs for generic nodes↩︎

It is a basic result [34] that the regularised posterior mean \[m_{{\cal D}, \lambda}({\boldsymbol{x}}) = {\boldsymbol{k}}_X({\boldsymbol{x}})^{\intercal}({\boldsymbol{K}}_{XX} + \lambda {\boldsymbol{I}}_N )^{-1} {\boldsymbol{y}}\] is the unique solution to a regularised least-squares problem: \[\label{eq:mean-regularised-least-squeres} m_{{\cal D}, \lambda} = \mathop{\mathrm{arg\,min}}_{s \in \mathcal{H}(k)} \Big\{ \sum_{i=1}^N [f({\boldsymbol{x}}_i) - s({\boldsymbol{x}}_i)]^2 + \lambda \norm[0]{s}_{\mathcal{H}(k)}^2 \Big \}.\tag{42}\] Because \(m_{\cal D}\) is an interpolant to \(f\), which is to say that \(m_{\cal D}({\boldsymbol{x}}_i) = f({\boldsymbol{x}}_i)\) for every \(i = 1,\ldots,N\), and an element of \(\mathcal{H}(k)\), we have from 42 that \[\label{eq:regularisation-norm-bound} \sum_{i=1}^N [f({\boldsymbol{x}}_i) - m_{{\cal D},\lambda}({\boldsymbol{x}}_i)]^2 + \lambda \norm[0]{m_{{\cal D},\lambda}}_{\mathcal{H}(k)}^2 \leq \lambda \norm[0]{m_{{\cal D}}}_{\mathcal{H}(k)}^2.\tag{43}\] Consequently, \[\label{eq:regularised-nodeset-error} \Big(\sum_{i=1}^N [f({\boldsymbol{x}}_i) - m_{{\cal D},\lambda}({\boldsymbol{x}}_i)]^2 \Big)^{1/2} \leq \sqrt{\lambda} \norm[0]{m_{{\cal D}}}_{\mathcal{H}(k)}.\tag{44}\] We also need to recall the minimum-norm interpolation property 16 of the conditional mean, which states that \[\label{eq:minimum-norm-interpolation-proofs} m_{\cal D}= \mathop{\mathrm{arg\,min}}_{s \in \mathcal{H}(k)} \Set[]{ \norm[0]{s}_{\mathcal{H}(k)} }{ s({\boldsymbol{x}}_i) = f({\boldsymbol{x}}_i) \text{ for every } i = 1, \ldots, N},\tag{45}\] and its particular consequence \(\norm[0]{m_{\cal D}}_{\mathcal{H}(k)} \leq \norm[0]{f}_{\mathcal{H}(k)}\) if \(f\in \mathcal{H}(k)\).

Observe that 8 and its regularised variant, followed by Jensen’s inequality, yield \[\label{eq:bq-integral-bound-interpolant} \begin{align} \abs[0]{ I_P(f) - \mu_{{\cal D}, \lambda} } &= \Big\lvert \int_D[ f({\boldsymbol{x}}) - m_{{\cal D},\lambda}({\boldsymbol{x}}) ] \dif P({\boldsymbol{x}}) \Big\rvert \\ &\leq \int_D\abs[0]{ f({\boldsymbol{x}}) - m_{{\cal D},\lambda}({\boldsymbol{x}}) } \dif P({\boldsymbol{x}}). \end{align}\tag{46}\] Assumption 5.1 and a second application of Jensen’s inequality then produce \[\label{eq:error-L1L2} \begin{align} \abs[0]{ I_P(f) - \mu_{N, \lambda} } &\leq p_\mathrm{max} \int_{[0, 1]^d} \abs[0]{ f({\boldsymbol{x}}) - m_{{\cal D},\lambda}({\boldsymbol{x}}) } \dif {\boldsymbol{x}}\\ &\leq p_\mathrm{max} \norm[0]{f- m_{{\cal D},\lambda}}_{L^2([0,1]^d)}. \end{align}\tag{47}\] To bound the integration error it is therefore sufficient to bound the error of approximating \(f\) with \(m_{{\cal D},\lambda}\) in \(L^2([0,1]^d)\) norm. This is what is done in most of the proofs below.

Proposition [prop:generic-guarantee-2] is likely an old result. We follow the proof of Corollary 2.8 in [232].

7.3.2 Proofs for Section 5.1.1↩︎

In this section we use \(\lesssim\) to denote an inequality which holds up to a constant coefficient independent of \(f\), \(X\), the nugget \(\lambda\), the length-scale \(\ell\) and the scale \(\sigma\). Let \(\alpha > d/2\) and \(\nu = \alpha - d/2 > 0\). Here we assume that \(k\) is the isotropic Matérn kernel of order \(\nu\) given by \[\begin{align} k({\boldsymbol{x}}, {\boldsymbol{x}}') &= \Phi_k(\norm[0]{{\boldsymbol{x}}- {\boldsymbol{x}}'}) \\ &= \sigma^2 \frac{2^{1-\nu}}{\Gamma(\nu)} \Big( \frac{\sqrt{2\nu} \norm[0]{{\boldsymbol{x}}- {\boldsymbol{x}}'}}{\ell} \bigg)^\nu K_\nu \Big( \frac{\sqrt{2\nu} \norm[0]{{\boldsymbol{x}}- {\boldsymbol{x}}'}}{\ell} \Big), \end{align}\] where \(\Gamma\) denotes the gamma function and \(K_\nu\) the modified Bessel function of the second kind. The length-scale \(\ell > 0\) controls how “wiggly” the processes are, while the regularity parameter \(\nu > 0\) determines the differentiability of a process with Matérn covariance. Because the Fourier transform of a Matérn kernel is \[\widehat{\Phi}_k({\boldsymbol{\xi}}) = \sigma^2 \frac{2^{\:d/2} \Gamma(\nu + d/2)}{\Gamma(\nu)} \Big( \frac{2\nu}{\ell^2} \Big)^\nu \Big( \frac{2\nu}{\ell^2} + \norm[0]{{\boldsymbol{\xi}}}^2 \Big)^{-(\nu + d/2)},\] one may compute from the Fourier definition of the Sobolev norm that the RKHS of a Matérn kernel is norm-equivalent to a Sobolev space of order \(\nu\), which in this case means that \(\mathcal{H}(k) = H^\alpha([0, 1]^d)\) as sets of functions and \[\label{eq:norm-equivalence} C_k \norm[0]{g}_{H^\alpha([0,1]^d)} \leq \norm[0]{g}_{\mathcal{H}(k)} \leq C_k' \norm[0]{g}_{H^\alpha([0,1]^d)}.\tag{48}\] for the positive constants \[\label{eq:matern-norm-constants} C_k = c \frac{1}{\sigma \ell^{d/2}} \max\Big\{1, \frac{\sqrt{2\nu}}{\ell} \Big\}^{-1}, \quad C_k' = c' \frac{1}{\sigma \ell^{d/2}} \max\Big\{1, \frac{\ell}{\sqrt{2\nu}} \Big\},\tag{49}\] where \(c\) and \(c'\) are positive constants that depend only on \(d\) and \(\nu\). See, for example, Lemma 3.4 and its proof in [212] for the details of this derivation. Consequently, \[\label{eq:matern-norm-constants-ratio}C_k^{-1} C_k' \lesssim \max\Big\{1, \frac{\ell}{\sqrt{2\nu}} \Big\} \max\Big\{1, \frac{\sqrt{2\nu}}{\ell} \Big\} = \max\Big\{ \frac{\ell}{\sqrt{2\nu}}, \frac{\sqrt{2\nu}}{\ell} \Big\}.\tag{50}\]

The following result is Theorem 3.2 in [233] with \(n = d\), \(p = q = 2\), \(s = \tau\), \(r = \kappa\), \(\varkappa = 2\) and \(\Omega = [0, 1]^d\). The requirement in [233] that \(h_X\) be sufficiently small can be eliminated by inflating the constant \(C\).

The following result is a version of Theorem 4.2 in [234] that makes explicit the dependency on the norm-equivalence constants. See also Theorem 16 in [22]. A proof is contained in the proof of Theorem 4.2 in [204].

7.3.3 Proofs for Section 5.1.3↩︎

Theorems 5.4 and 5.5 have previously appeared as Theorem 2.20 in [35]. Theoretical guarantees for explicitly active Bayesian quadrature based on the square exponential kernel may be found in [203]. Further guarantees, including on the unbounded domain \(D= \mathbb{R}^d\), could be derived from the results in [218], [236]–[238].

7.4 Proofs for Section 5.2↩︎

Acknowledgements↩︎

MM’s work was mainly performed while at the University of Tübingen. MM gratefully acknowledges financial support by the European Research Council through ERC StG Action 757275 / PANAMA; the DFG Cluster of Excellence “Machine Learning - New Perspectives for Science”, EXC 2064/1, project number 390727645; the German Federal Ministry of Education and Research (BMBF) through the Tübingen AI Center (FKZ: 01IS18039A); and funds from the Ministry of Science, Research and Arts of the State of Baden-Württemberg. TK was generously supported by the Research Council of Finland projects 338567 (“Scalable, adaptive and reliable probabilistic integration”), 359183 (“Flagship of Advanced Mathematics for Sensing, Imaging and Modelling”), and 368086 (“Inference and approximation under misspecification”). TK acknowledges the research environment provided by ELLIS Institute Finland. We thank Vesa Kaarnioja for correcting some of our misapprehensions about Bayesian quadrature with periodic covariance.

References↩︎

---

1. In fact, the square exponential covariance function is the only isotropic covariance function that is the product of its one-dimensional versions [26].↩︎

2. The term Bayesian quadrature* was popularised by [5]. The earliest occurrence that we have found is in [28]. [4] indicates an origin no later than 1985. In the early literature, Bayes–Hermite quadrature is used to refer to conjugate BQ when the GP prior has a square exponential covariance and \(P\) is a Gaussian distribution [5], [29], [30].*↩︎

3. See [43] and [44] for some theorems on preservation of Gaussianity under transformations.↩︎

4. [65] provide a historical account on the pioneering role played by Suldin and Larkin in the development of probabilistic numerical computation. We note that earlier constructions of optimal quadrature rules that do not utilise RKHSs may be found in [66] and [67].↩︎

5. Kernel quadrature is occasionally used to refer any quadrature rule designed so as to minimize its worst-case error in an RKHS [99]. Such unconditional usage of the term is somewhat ambitious: Every Hilbert space of functions that one can construct is an RKHS because such a Hilbert space is an RKHS if and only if the point evaluation functionals \(f \mapsto f({\boldsymbol{x}})\) are bounded, the existence of unbounded functionals depends on the axiom of choice [100], [101]. For example, under this definition quasi-Monte Carlo integration would fall under the umbrella of kernel quadrature [3].↩︎

6. As \(I_P(\mathsf{f}) \mid {\cal D}\) is a distribution, any datatype that can represent a distribution is valid. For Gaussians, a common way is to return the sufficient statistics: mean \(\mathbb{E}[ I_P(\mathsf{f}) \mid {\cal D}]\) and variance \(\mathrm{Var}[ I_P(\mathsf{f}) \mid {\cal D}]\).↩︎

7. The name is likely inspired by [136] who learn the transformation \(\varphi\) from data in the context of GP regression and call their method “warped Gaussian processes”. In the context of Bayesian quadrature, the transformation \(\varphi\) is usually fixed, but could in principle also be learned in a similar way.↩︎

8. Bayesian quadrature is related to the field of probabilistic numerics [13], [137], [138] which frames integration, among other numerical algorithms, as inference. If we were to follow [138] who characterize a probabilistic numerical method as “Bayesian” or “non-Bayesian”, Bayesian quadrature as in Definition 2.1 admitting exact or approximate inference according to our taxonomy would naturally be called Bayesian Bayesian quadrature or non-Bayesian Bayesian quadrature, respectively.↩︎

9. There is a vast literature on scalable algorithms, such as random Fourier features, for approximate Gaussian process inference [139]–[141]. Such algorithms have seen little use in the context of Bayesian quadrature.↩︎

10. We are not aware of explicit use of stochastic optimizers in BQ, although stochastic optimizers have been used in the related field of Bayesian optimization [146]. However, even deterministic optimizers usually contain a degree of randomness. For example, the initial point of the iteration needs to be chosen which may affect which maximizer is found and hence which node is sampled.↩︎

11. The mode of sampling at this point simply describes characteristics of the node selection option. Unlike data-dependent estimators that may inherit the randomness property from the mode of sampling, Bayesian models (the BQ model is a Bayesian model) do not. However, estimators derived from the BQ model such as the mean or MAP estimator may be analyzed in this way.↩︎

12. See [147] for an application of this method to the selection of lithium-ion battery models.↩︎

13. The situation is slightly different in the related Bayesian optimization (BO) method where plug-and-play is often feasible. This is because the inference method in BO has a less prominent role and influences node sampling but not the return object, which is \(\mathop{\mathrm{arg\,min}}_{ {\boldsymbol{x}}_i \in X} f({\boldsymbol{x}}_i)\). For the same reason, BO methods always require model-dependent sampling in order to be meaningful.↩︎

14. The formulation of the active sampler is supported by the von Neumann–Morgenstern utility theorem [149] which shows, under some conditions, that the decision of a rational agent under uncertainty can be constructed as a policy that maximizes an expected utility function. For an in-depth derivation and motivation of active sampling in the related field of Bayesian optimization see, for example, Section 5.1 in [150]. The treatment in BQ is analogous and not repeated here. For Bayesian decision theory for GPs, see Section 2.4 in [23].↩︎

15. Non-myopic sequential, active sampling has so far received little attention in BQ. The only approach we are aware of is due to [151]. A simple example can illustrate the difference between myopic and non-myopic sampling: Consider a sampling horizon of \(N=2\) and let \(\mathsf{f}\) be a Brownian motion as in the example of Section 3.3.2. A non-myopic sampler would place the first node at \(x=4/5\), anticipating a future node on \(x=2/5\), while the initial node of the myopic sampler would be placed sub-optimally at \(x=2/3\) [13].↩︎

16. Note that non-myopic sampling and batched sampling are orthogonal concepts. In principle, there exist batched (\(n>1\)) and non-batched (\(n=1\)) methods that can be either myopic or not. Batched BQ methods that use maximally myopic acquisition functions are proposed in [152] and [153]. Further, batched and non-batched methods need of course not even be active.↩︎

17. It is not necessary to have access to gradients of the integrand \(f\) since the gradients of the acquisition function \(a\) depend only on gradients of the model \(\mathsf{f}\), which is under the control of the practitioner.↩︎

18. The term adaptive has several additional meanings in numerical integration literature: in adaptive quadrature the integral is approximated by a sum of quadratures rules on sub-domains that are adaptively refined until a desired error tolerance is reached [156], while dimension-adaptive quadrature attempts to find the dimensions which contribute most to the integral [157]. Although the modelling framework would be well-suited to the task, no dimension-adaptive Bayesian quadrature methods appear to have been developed so far.↩︎

19. Our use of the term has nothing in common with structured kernel interpolation that is related to scalable GP regression [143].↩︎

20. Apart from a safeguard condition that caps the number of \(f\)-evaluations at some large \(N_\mathrm{max}\).↩︎

21. https://github.com/mmahsereci/kernel_embedding_dictionary↩︎

22. However, it is always good to keep in mind some properties of the prior when selecting the nodes. For example, equispaced nodes are not well suited for approximation of infinitely differentiable functions [175]. One should avoid equispaced nodes if the kernel is, for example, the square exponential.↩︎

23. Active sampling typically assumes that the cost of evaluating \(f({\boldsymbol{x}})\) is the same for all \({\boldsymbol{x}}\) making maximizing the acquisition function reasonable. However, if the cost varies with \({\boldsymbol{x}}\), an alternative approach may be preferable. The only one we are aware of in the context of BQ is due to [36] who maximize acquisition rates \(r(\tilde{{\boldsymbol{x}}})= a(\tilde{{\boldsymbol{x}}})/c(\tilde{{\boldsymbol{x}}})\) based on some positive cost function \(c(\tilde{{\boldsymbol{x}}})\). For example, if \(a\) is MI and \(c\) measures the cost in time or monetary value, then \(r\) would yield nodes that maximize bits per time or bits per dollar. In this cost-sensitive scenario, even if \(\mathop{\mathrm{arg\,max}}a_1(\tilde{{\boldsymbol{x}}}) = \mathop{\mathrm{arg\,max}}a_2(\tilde{{\boldsymbol{x}}})\) for two different acquisition functions \(a_1\) and \(a_2\) (such as IVR and NIV), generally \(\mathop{\mathrm{arg\,max}}r_1(\tilde{{\boldsymbol{x}}}) \neq \mathop{\mathrm{arg\,max}}r_2(\tilde{{\boldsymbol{x}}})\), even if they share the same cost function \(c\). In fact, certain acquisition functions, like NIV, cannot be converted to meaningful acquisition rates at all. Cost-sensitive evaluations further play a role in multi-source and multi-fidelity active BQ where the acquisition strategy requires careful consideration.↩︎

24. This principle only applies when one uses the basis \(\{k(\cdot, {\boldsymbol{x}}_i)\}_{i=1}^N\) of kernel translates. Numerically more stable bases exist [194], [195].↩︎

25. The approach that we describe is closely related to ordinary and universal kriging [196].↩︎

26. However, bounds of order \(h_X^{2\alpha}\) or \(N^{-2\alpha}\) are available if \(f\) has smoothness \(2\alpha\) and satisfies certain boundary conditions [204], [206]–[210]. This case is known as superconvergence or improved rate of kernel-based interpolation.↩︎

27. As the library seems to be under active development, this may have changed.↩︎

28. When referring to functions we use the word “stationary” as in [23] who write “informally, stationarity means that the functions look similar at all \(x\) locations”.↩︎

29. Occasional outliers, mainly from failed runs, have been discarded when computing the averages. Outliers are defined as scores that lie more than \(1.5\times \Delta_{1-3}\) outside of the interquartile range \(\Delta_{1-3}\) of all scores. This is the same definition that SciPy uses for its default boxplots.↩︎

30. The calibration score is a simple proxy to get a general idea of the calibration and provides by no means a complete picture. We refer to [229] for a formal treatment of calibration (see in particular Remark 9).↩︎