On the positivity and magnitudes of Bayesian quadrature weights

TL;DR

This paper analyzes the properties of Bayesian quadrature weights, focusing on conditions for positivity and bounds on their magnitudes, to improve stability and robustness of quadrature rules.

Contribution

It provides theoretical conditions for positive weights in univariate cases and bounds on weight magnitudes based on point distribution and kernel properties.

Findings

Weights are positive if points minimize posterior variance and kernel is totally positive.

Weight magnitudes are bounded by fill distance and separation radius in Sobolev spaces.

Numerical examples suggest potential for further generalizations and improvements.

Abstract

This article reviews and studies the properties of Bayesian quadrature weights, which strongly affect stability and robustness of the quadrature rule. Specifically, we investigate conditions that are needed to guarantee that the weights are positive or to bound their magnitudes. First, it is shown that the weights are positive in the univariate case if the design points locally minimise the posterior integral variance and the covariance kernel is totally positive (e.g., Gaussian and Hardy kernels). This suggests that gradient-based optimisation of design points may be effective in constructing stable and robust Bayesian quadrature rules. Secondly, we show that magnitudes of the weights admit an upper bound in terms of the fill distance and separation radius if the RKHS of the kernel is a Sobolev space (e.g., Mat\'ern kernels), suggesting that quasi-uniform points should be used. A number of numerical examples demonstrate that significant generalisations and improvements appear to be possible, manifesting the need for further research.

Figures

Equations

$$I_{\nu }(f):={\int }_{\Omega }fd\nu ,$$

$$\mathbb{E}\big{[}I_{\nu}(f_{\text{{\tiny GP}}})\mid\mathcal{D}\big{]}=\sum_{i=1}^{n}w_{X,i}^{\text{{\tiny BQ}}}f(\boldsymbol{x}_{i})\approx\int_{\Omega}f\operatorname{d\!}\nu,$$

$$[\boldsymbol{K}_{X}]_{ij}\coloneqq\mathrm{Cov}\big{[}f_{\text{{\tiny GP}}}(\boldsymbol{x}_{j}),f_{\text{{\tiny GP}}}(\boldsymbol{x}_{i})\big{]}=k(\boldsymbol{x}_{j},\boldsymbol{x}_{i})$$

$$\begin{split}\mu_{X,f}(\boldsymbol{x})&\coloneqq\mathbb{E}\big{[}f_{\text{{\tiny GP}}}(\boldsymbol{x})\mid\mathcal{D}\big{]}\\ &=\boldsymbol{k}_{X}(\boldsymbol{x})^{\textsf{T}}\boldsymbol{K}_{X}^{-1}\boldsymbol{f}_{X}\end{split}$$

$$\begin{split}\sigma^{2}_{X}(\boldsymbol{x},\boldsymbol{x}^{\prime})&\coloneqq\mathrm{Cov}\big{[}f_{\text{{\tiny GP}}}(\boldsymbol{x}),f_{\text{{\tiny GP}}}(\boldsymbol{x}^{\prime})\mid\mathcal{D}\big{]}\\ &=k(\boldsymbol{x},\boldsymbol{x}^{\prime})-\boldsymbol{k}_{X}(\boldsymbol{x})^{\textsf{T}}\boldsymbol{K}_{X}^{-1}\boldsymbol{k}_{X}(\boldsymbol{x}^{\prime}),\end{split}$$

$$\displaystyle\begin{split}I_{X}^{\text{{\tiny BQ}}}(f)&\coloneqq\mathbb{E}\big{[}I_{\nu}(f_{\text{{\tiny GP}}})\mid\mathcal{D}\big{]}\\ &=\int_{\Omega}\mathbb{E}\big{[}f_{\text{{\tiny GP}}}(\boldsymbol{x})\mid\mathcal{D}\big{]}\operatorname{d\!}\nu(\boldsymbol{x})\\ &=\boldsymbol{k}_{\nu,X}^{\textsf{T}}\boldsymbol{K}^{-1}_{X}\boldsymbol{f}_{X},\end{split}$$

$$\displaystyle\begin{split}\mathbb{V}_{X}^{\text{{\tiny BQ}}}&\coloneqq\text{Var}[I(f_{\text{{\tiny GP}}})\mid\mathcal{D}]\\ &=\int_{\Omega}\int_{\Omega}\mathrm{Cov}\big{[}f_{\text{{\tiny GP}}}(\boldsymbol{x}),f_{\text{{\tiny GP}}}(\boldsymbol{x}^{\prime})\mid\mathcal{D}\big{]}\operatorname{d\!}\nu(\boldsymbol{x})\operatorname{d\!}\nu(\boldsymbol{x}^{\prime})\\ &=I_{\nu}(k_{\nu})-\boldsymbol{k}_{\nu,X}^{\textsf{T}}\boldsymbol{K}_{X}^{-1}\boldsymbol{k}_{\nu,X},\end{split}$$

$$I_{\nu }(k_{\nu })={\int }_{\Omega }{k}_{\nu }(\mathbf{x})d\nu (x)={\int }_{\Omega }{\int }_{\Omega }k(\mathbf{x},{\mathbf{x}}^{\prime })d\nu ({\mathbf{x}}^{\prime })d\nu (\mathbf{x}).$$

$$I_X^{\text{BQ}}(f)=({\mathbf{w}}_X^{\text{BQ}}{)}^{\text{T}}{\mathbf{f}}_X=\sum _{i=1}^n{w}_{X,i}^{\text{BQ}}f({\mathbf{x}}_i),$$

$${\mathbf{w}}_X^{\text{BQ}}:=(w_{X,i}^{\text{BQ}}{)}_{i=1}^n:={\mathbf{K}}_X^{-1}{\mathbf{k}}_{\nu ,X}\in {\mathbb{R}}^n.$$

$$I_X^{\text{BQ}}(k_{{\mathbf{x}}_i})=I_{\nu }(k_{{\mathbf{x}}_i})=k_{\nu }({\mathbf{x}}_i)\text{ for each }i=1,\dots ,n,$$

$$\sum _{i=1}^{n}k({\mathbf{x}}_j,{\mathbf{x}}_i)w_{X,i}^{\text{BQ}}=k_{\nu }({\mathbf{x}}_j).$$

$$\begin{array}{rl}{\mathbb{V}}_X^{\text{BQ}}& =I_{\nu }(k_{\nu })-{\mathbf{k}}_{\nu ,X}^{\text{T}}{\mathbf{K}}_X^{-1}{\mathbf{k}}_{\nu ,X}\\ & =I_{\nu }(k_{\nu })-({\mathbf{w}}_X^{\text{BQ}}{)}^{\text{T}}{\mathbf{k}}_{\nu ,X}\\ & =I_{\nu }(k_{\nu })-I_X^{\text{BQ}}(k_{\nu }).\end{array}$$

$${\mu }_{X,f}(\mathbf{x})={\mathbf{u}}_X(\mathbf{x}{)}^{\text{T}}{\mathbf{f}}_X=\sum _{i=1}^{n}f({\mathbf{x}}_i)u_{X,i}(\mathbf{x}).$$

$$u_{X,i}({\mathbf{x}}_j)={\delta }_{ij}:=\{\begin{array}{l}1\text{ if }i=j,\\ 0\text{ if }i\ne j.\end{array}$$

$$I_X^{\text{BQ}}(f)=\sum _{i=1}^{n}f({\mathbf{x}}_i)I_{\nu }(u_{X,i}).$$

$$\begin{split}e_{\mathcal{H}(k)}(X,\boldsymbol{w})^{2}&\coloneqq\sup_{\lVert f\rVert_{\mathcal{H}(k)}\leq 1}\Bigg{\lvert}\int_{\Omega}f\operatorname{d\!}\nu-\sum_{i=1}^{n}w_{i}f(\boldsymbol{x}_{i})\Bigg{\rvert}\\ &=I_{\nu}(k_{\nu})-2\boldsymbol{w}^{\textsf{T}}\boldsymbol{k}_{\nu,X}+\boldsymbol{w}^{\textsf{T}}\boldsymbol{K}_{X}\boldsymbol{w}.\end{split}$$

$${\mathbf{w}}_X^{\text{BQ}}=\underset{\mathbf{w}\in {\mathbb{R}}^n}{\mathrm{argmin}}{e}_{\mathcal{H}(k)}(X,\mathbf{w})$$

$${\mathbb{V}}_X^{\text{BQ}}=e_{\mathcal{H}(k)}(X,{\mathbf{w}}_X^{\text{BQ}}{)}^2.$$

$$e_{\mathcal{H}(k)}(X,{\mathbf{w}}_X^{\text{BQ}})={\Vert k_{\nu }-k_Q\Vert }_{\mathcal{H}(k)},k_Q:=\sum _{i=1}^n{w}_{X,i}^{\text{BQ}}{k}_{{\mathbf{x}}_i}.$$

$$\begin{split}\mathinner{\lvert I_{\nu}(f)-I_{X}^{\text{{\tiny BQ}}}(f)\rvert}&=\mathinner{\!\bigl{\lvert}\langle k_{\nu}-k_{Q},f\rangle_{\mathcal{H}(k)}\bigr{\rvert}}\\ &\leq\mathinner{\lVert f\rVert}_{\mathcal{H}(k)}\mathinner{\lVert k_{\nu}-k_{Q}\rVert}_{\mathcal{H}(k)}\\ &=\mathinner{\lVert f\rVert}_{\mathcal{H}(k)}e_{\mathcal{H}(k)}(X,\boldsymbol{w}_{X}^{\text{{\tiny BQ}}}).\end{split}$$

$$\sum _{i=1}^n{w}_{i}P(x_i)={\int }_a^{b}P(x)d\nu (x)$$

$$0<{\int }_a^b{G}_i(x)d\nu (x)=\sum _{i=1}^n{w}_j{G}_i(x_j)=w_i.$$

$$P^{(j_{i})}(x_{i})\coloneqq\frac{\mathrm{d}^{j_{i}}}{\operatorname{d\!}x^{j_{i}}}P(x)\Bigr{|}_{x=x_{i}}=0$$

$$\phi^{(j)}(x)=\mathrm{e}^{x}\big{[}P(x)+c_{1}P^{(1)}(x)+\cdots+c_{j}P^{(j)}(x)\big{]}$$

$$\sum _{i=1}^n{\alpha }_i{\varphi }_i(x_j)=0\text{ for every }j=1,\dots ,n$$

$$s_{X,f}(x_j)=\sum _{i=1}^n{\alpha }_i{\varphi }_i(x_j)=f(x_j)$$

$$s_{X,f}(x)=\mathbf{\varphi }(x{)}^{\text{T}}\mathbf{\alpha }=\mathbf{\varphi }(x{)}^{\text{T}}{\mathbf{V}}_X^{-\text{T}}{\mathbf{f}}_X$$

$$\{(x_i,f^{(j_i)}(x_i)):i=1,\dots ,m\text{ and }{j}_i=0,\dots ,q_i-1\}.$$

$$s_{X,\mathbf{q},f}^{(j_i)}(x_i)=f^{(j_i)}(x_i)$$

$$s_{X,\mathbf{q},f}^{(j_i)}(x_i)=\sum _{l=1}^n{\alpha }_l{\varphi }_l^{(j_i)}(x_i)=f^{(j_i)}(x_i).$$

Full text

1 Introduction

This article is concerned with Bayesian quadrature OHagan1991 ; RasmussenGhahramani2002 ; Briol2017 , a probabilistic approach to numerical integration and an example of a probabilistic numerical method Larkin1972 ; Hennig2015 ; Cockayne2017 . Let $\Omega$ be a subset of $\mathbb{R}^{d}$, $d\geq 1$, and $\nu$ a Borel probability measure on $\Omega$. Given an integrand $f\colon\Omega\to\mathbb{R}$, the task is to approximate the integral

[TABLE]

the solution of which is assumed not to be available in closed form. In Bayesian quadrature, a user specifies a prior distribution over the integrand as a Gaussian process $f_{\text{{\tiny GP}}}\sim\mathcal{GP}(0,k)$ by choosing a positive-definite covariance kernel $k\colon\Omega\times\Omega\to\mathbb{R}$, so as to faithfully represent their knowledge about the integrand, such as its smoothness. The user then evaluates the true integrand at chosen design points $X=\{\boldsymbol{x}_{1},\dots,\boldsymbol{x}_{n}\}\subset\Omega$. By regarding the pairs $\mathcal{D}\coloneqq\{(\boldsymbol{x}_{i},f(\boldsymbol{x}_{i}))\}_{i=1}^{n}$ thus obtained as “observed data”, the posterior distribution $I_{\nu}(f_{\text{{\tiny GP}}})\mid\mathcal{D}$ becomes a Gaussian random variable. This posterior distribution is useful for uncertainty quantification and decision making in subsequent tasks; this is one factor that makes Bayesian quadrature a promising approach in modern scientific computation, where quantification of discretisation errors is of great importance Briol2017 ; Oates2017 .

In Bayesian quadrature, the mean of the posterior over the integral is used as a quadrature estimate. The mean given is as a weighted average of function values:

[TABLE]

where $w_{X,1}^{\text{{\tiny BQ}}},\dots,w_{X,n}^{\text{{\tiny BQ}}}\in\mathbb{R}$ are the weights computed with the kernel $k$, design points $X$ and the measure $\nu$ (see Section 2.1 for details). This form is similar to (quasi) Monte Carlo methods, where $\boldsymbol{x}_{1},\dots,\boldsymbol{x}_{n}$ are (quasi) random points from a suitable proposal distribution and $w_{1},\dots,w_{n}$ are the associated importance weights, positive by definition. This similarity naturally leads to the following question: Are the weights $w_{X,1}^{\text{{\tiny BQ}}},\dots,w_{X,n}^{\text{{\tiny BQ}}}$ of Bayesian quadrature positive? These weights are derived with no explicit positivity constraint, so in general some of them can be negative, which is observed in (HuszarDuvenaud2012, , Section 3.1.1). Therefore, the question can be stated as: Under which conditions on the points and the kernel are the weights guaranteed to be positive?

This question is important both conceptually and practically. On the conceptual side, positive weights are more natural, given that the weighted sample $(w_{i},\boldsymbol{x}_{i})_{i=1}^{n}$ can be interpreted as an approximation of the positive probability measure $\nu$; in fact, the Bayesian quadrature weights provide the best approximation of the representer of $\nu$ in the reproducing kernel Hilbert space (RKHS) of the covariance kernel, provided that $\boldsymbol{x}_{1},\dots,\boldsymbol{x}_{n}$ are fixed (see Section 2.2). Thus, if the weights are positive, then each weight $w_{i}$ can be interpreted as representing the “importance” of the associated point $\boldsymbol{x}_{i}$ for approximating $\nu$. This interpretation may be more acceptable to users familiar with Monte Carlo methods, encouraging them to adopt Bayesian quadrature.

On the practical side, quadrature rules with positive weights enjoy the advantage of being numerically more stable against errors in integrand evaluations. In fact, besides Monte Carlo methods, many other practically successful or in some sense optimal rules have positive weights. Some important examples include Gaussian (Gautschi2004, , Section 1.4.2) and Clenshaw–Curtis quadrature ClenshawCurtis1960 and their tensor product extensions. Other domains besides subsets of $\mathbb{R}^{d}$ have also received their share of attention. For instance, positive-weight rules on the sphere are constructed in Mhaskar2001 and interesting results connecting fill-distance and positivity of the weights of quadrature rules on compact Riemannian manifolds appear in Breger2018 . It is also known that, in some typical function classes, such as Sobolev spaces, optimal rates of convergence can be achieved by considering only positive-weight quadrature rules; see for instance (Novak1999, , Section 1) and references therein. Therefore, if one can find conditions under which Bayesian quadrature weights are positive, then these conditions may be used as guidelines in construction of numerically stable Bayesian quadrature rules.

This article reviews existing, and derives new, results on properties of the Bayesian quadrature weights, focusing in particular on their positivity and magnitude. One of our principal aims is to stimulate new research on quadrature weights in the context of probabilistic numerics. While convergence rates of Bayesian quadrature rules have been studied extensively in recent years Briol2017 ; Kanagawa2016 ; Kanagawa2018 , analysis of the weights themselves has not attracted much attention. On the other hand, the earliest work Larkin1970 ; RichterDyn1971a ; BarrarLoeb1976 (see Oettershagen2017 for a recent review) done in the 1970s on kernel-based quadrature already revealed certain interesting properties of the Bayesian quadrature weights. These results seem not well-known in the statistics and machine learning community. Moreover, there are some useful results from the literature on scattered data approximation DeMarchiSchaback2010 , which can be used to analyse the properties of Bayesian quadrature weights. The basics of Bayesian quadrature are reviewed in Section 2 while the main contents, including simulation results, of the article are presented in Sections 3 and 4.

In Section 3, we present results concerning positivity of the Bayesian quadrature weights. We discuss results on the number of the weights that must be positive, focusing on the univariate case and totally positive kernels (Definition 2). Corollary 1, the main result of this section, states that all the weights are positive if the design points are locally optimal. A practically relevant consequence of this result is that it may imply that the weights are positive if the design points are obtained by gradient descent, which is guaranteed to provide locally optimal points (see e.g. LeeSimJorRec16 ).

Section 4 focuses on results on the magnitudes of the weights. More specifically, we discuss the behaviour of the sum of absolute weights, $\sum_{i=1}^{n}\mathinner{\lvert w_{X,i}^{\text{{\tiny BQ}}}\rvert}$, that strongly affects stability and robustness of Bayesian quadrature. If this quantity is small, the quadrature rule is robust against misspecification of the Gaussian process prior Kanagawa2018 and errors in integrand evaluations Forster1993 and kernel means (SommarivaVianello2006, , pp. 298–300). This quantity is also related to the numerical stability of the quadrature rule. Using a result on stability of kernel interpolants by De Marchi and Schaback DeMarchiSchaback2010 , we derive an upper bound on the the sum of absolute weights for some typical cases where the Gaussian process has finite degree of smoothness and the RKHS induced by the covariance kernel is norm-equivalent to a Sobolev space.

2 Bayesian quadrature

This section defines a Bayesian quadrature rule as the integral of the posterior of Gaussian process used to model the integrand. We also discuss the equivalent characterisation of this quadrature rule as the worst-case optimal integration rule in the RKHS $\mathcal{H}(k)$ induced by the covariance kernel $k$ of the Gaussian process.

2.1 Basics of Bayesian quadrature

In standard Bayesian quadrature OHagan1991 ; Minka2000 ; Briol2017 , the deterministic integrand $f\colon\Omega\to\mathbb{R}$ is modelled as a Gaussian process. The integrand is assigned a zero-mean Gaussian process prior $f_{\text{{\tiny GP}}}\sim\mathcal{GP}(0,k)$ with a positive-definite covariance kernel $k$. This is to say that for any $n\in\mathbb{N}$ and any distinct points $X=\{\boldsymbol{x}_{1},\ldots,\boldsymbol{x}_{n}\}\subset\Omega$ we have $(f_{\text{{\tiny GP}}}(\boldsymbol{x}_{1}),\ldots,f_{\text{{\tiny GP}}}(\boldsymbol{x}_{n}))\sim\mathrm{N}(\boldsymbol{0},\boldsymbol{K}_{X})$, with

[TABLE]

the $n\times n$ positive-definite (and hence invertible) kernel matrix. Conditioning on the data $\mathcal{D}=\{(\boldsymbol{x}_{i},f(\boldsymbol{x}_{i}))\}_{i=1}^{n}$, consisting of evaluations $\boldsymbol{f}_{X}\coloneqq(f(\boldsymbol{x}_{i}))_{i=1}^{n}\in\mathbb{R}^{n}$ of $f$ at points $X$, yields a Gaussian posterior process with the mean

[TABLE]

and covariance

[TABLE]

where the $n$-vector $\boldsymbol{k}_{X}(\boldsymbol{x})$ has the elements $[\boldsymbol{k}_{X}(\boldsymbol{x})]_{i}=k(\boldsymbol{x},\boldsymbol{x}_{i})$. Note that the posterior covariance only depends on the points, not on the integrand, and that the posterior mean interpolates the data (i.e., $\mu_{X,f}(\boldsymbol{x_{i}})=f(\boldsymbol{x}_{i})$ for $i=1,\ldots,n$). Accordingly, the posterior mean often goes by the name kernel interpolant or, if the kernel is isotropic, radial basis function interpolant.

Due to the linearity of the integration operator, the posterior of the integral becomes a Gaussian distribution $I(f_{\text{{\tiny GP}}})\mid\mathcal{D}\sim\mathcal{N}(I_{X}^{\text{{\tiny BQ}}}(f),\mathbb{V}_{X}^{\text{{\tiny BQ}}})$ with the mean and variance

[TABLE]

where $k_{\nu}(\boldsymbol{x})\coloneqq\int_{\Omega}k(\cdot,\boldsymbol{x})\operatorname{d\!}\nu(\boldsymbol{x})$ is the kernel mean Smola2007 , $\boldsymbol{k}_{\nu,X}\in\mathbb{R}^{n}$ with $[\boldsymbol{k}_{\nu,X}]_{i}=k_{\nu}(\boldsymbol{x}_{i})$ and

[TABLE]

The integral mean $I_{X}^{\text{{\tiny BQ}}}(f)$ is used to approximate the true intractable integral $I_{\nu}(f)$ while the variance $\mathbb{V}_{X}^{\text{{\tiny BQ}}}$ is supposed to quantify epistemic uncertainty, due to partial information being used (i.e., a finite number of function evaluations) inherent to this approximation.

The integral mean $I_{X}^{\text{{\tiny BQ}}}(f)$ indeed takes the form a quadrature rule, a weighted sum of function evaluations:

[TABLE]

where $w_{X,1}^{\text{{\tiny BQ}}},\dots,w_{X,n}^{\text{{\tiny BQ}}}$ are the Bayesian quadrature weights given by

[TABLE]

The purpose of this article is to analyse the properties of these weights.

A particular property of a Bayesian quadrature rule is that the $n$ kernel translates $k_{\boldsymbol{x}_{i}}\coloneqq k(\cdot,\boldsymbol{x}_{i})$ are integrated exactly:

[TABLE]

which is derived from the fact that the $j$th equation of the linear system $\boldsymbol{K}_{X}\boldsymbol{w}_{X}^{\text{{\tiny BQ}}}=\boldsymbol{k}_{\nu,X}$ defining the weights is

[TABLE]

The left-hand side is precisely $I_{X}^{\text{{\tiny BQ}}}(k_{\boldsymbol{x}_{j}})$ while on the right-hand side we have $k_{\nu}(\boldsymbol{x}_{j})=I_{\nu}(k_{\boldsymbol{x}_{j}})$. Note also that the integral variance is the integration error of the kernel mean:

[TABLE]

Occasionally, it is instructive to interpret the weights as integrals of the Lagrange cardinal functions ${\boldsymbol{u}_{X}=(u_{X,i})_{i=1}^{n}}$ (see e.g. (Wendland2005, , Chapter 11)). These functions are defined as $\boldsymbol{u}_{X}(\boldsymbol{x})=\boldsymbol{K}_{X}^{-1}\boldsymbol{k}_{X}(\boldsymbol{x})$, from which it follows that

[TABLE]

Consequently, the cardinality property

[TABLE]

is satisfied, as can be verified by considering the interpolant $\mu_{X,g_{i}}$ to any function $g_{i}$ such that $g_{i}(\boldsymbol{x}_{j})=\delta_{ij}$. Since the integral mean is merely the integral of $\mu_{X,f}$, we have from (5) that

[TABLE]

That is, the $i$th Bayesian quadrature weight is the integral of the $i$th Lagrange cardinal function: $w_{X,i}^{\text{{\tiny BQ}}}=I_{\nu}(u_{X,i})$.

2.2 Reproducing kernel Hilbert spaces

An alternative interpretation of Bayesian quadrature weights is that they are, for the given points, the worst-case optimal weights in the reproducing kernel Hilbert space $\mathcal{H}(k)$ induced by the covariance kernel $k$. The material of this section is contained in, for example, (Briol2017, , Section 2), (Oettershagen2017, , Section 3.2) and (Karvonen2018, , Section 2). For a comprehensive introduction to RKHSs, see the monograph of Berlinet and Thomas-Agnan Berlinet2011 .

The RKHS induced by $k$ is the unique Hilbert space of functions characterised by (i) the reproducing property $\langle k_{\boldsymbol{x}},f\rangle_{\mathcal{H}(k)}=f(\boldsymbol{x})$ for every $f\in\mathcal{H}(k)$ and $\boldsymbol{x}\in\Omega$ and (ii) the fact that $k_{\boldsymbol{x}}\in\mathcal{H}(k)$ for every $\boldsymbol{x}\in\Omega$. The worst-case error in $\mathcal{H}(k)$ of a quadrature rule with points $X$ and weights $\boldsymbol{w}\in\mathbb{R}^{n}$ is

[TABLE]

It can be then shown that the Bayesian quadrature weights $\boldsymbol{w}_{X}^{\text{{\tiny BQ}}}$ are the unique minimiser of the worst-case error among all possible weights for these points:

[TABLE]

and

[TABLE]

Furthermore, the worst-case error can be written as the RKHS error in approximating the integration representer $k_{\nu}$ that satisfies $I_{\nu}(f)=\langle k_{\nu},f\rangle_{\mathcal{H}(k)}$ for all $f\in\mathcal{H}(k)$:

[TABLE]

From this representation and the Cauchy–Schwarz inequality it follows that

[TABLE]

For analysis of convergence of Bayesian quadrature rules as $n\to\infty$ it is therefore sufficient to analyse how the worst-case error (i.e., integral variance) behaves—as long as the integrand indeed lives in $\mathcal{H}(k)$. Convergence will be discussed in Section 4.

3 Positivity

This section reviews existing results on the positivity of the weights of Bayesian quadrature that can be derived in one dimension when the covariance kernel is totally positive. This assumption, given in Definition 2, is stronger than positive-definiteness but is satisfied by, for example, the Gaussian kernel. For most of the section we assume that $d=1$ and $\Omega=[a,b]$ for $a<b$. Furthermore, the measure $\nu$ is typically assumed to admit a density function with respect to the Lebesgue measure,111This can be usually relaxed to $I_{\nu}$ being a positive linear functional: $I_{\nu}(f)>0$ whenever $f$ is almost everywhere positive. an assumption that implies $I_{\nu}(f)>0$ if $f(x)>0$ for almost every $x\in\Omega$.

Positivity of the weights was actively investigated during the 1970s Richter1970 ; RichterDyn1971a ; RichterDyn1971b ; BarrarLoebWerner1974 ; BarrarLoeb1976 , and these results have been recently refined and collected by Oettershagen (Oettershagen2017, , Section 4). To simplify presentation, some of the results in this section are given in a slightly less general form than possible. Two of the most important results are

• •

Theorem 3.1: At least one half of the weights of any Bayesian quadrature rule are positive (Theorem 3.1).

• •

Corollary 1: All the weights are positive when the points are selected so that the integral posterior variance in (2) is locally minimised in the sense that each of its $n$ partial derivatives with respect to the integration points vanishes (Definition 3).

The latter of these results is particularly interesting since (i) it implies that points selected using a gradient descent algorithm may have positive weights and (ii) the resulting Bayesian quadrature rule is a positive linear functional and hence potentially well-suited for integration of functions that are known to be positive—a problem for which a number of transformation-based methods have been developed recently Osborne2012 ; Gunter2014 ; Chai2018 .

As no multivariate extension of the theory used to prove the aforementioned results appears to have been developed, we do not provide any general theoretical results on the weights in higher dimensions. However, some special cases based on, for example, tensor products are discussed in Sections 3.7 and 3.9 and two numerical examples are used to provide some evidence for the conjectures that multivariate versions of Theorem 3.1 and Corollary 1 hold.

It will turn out that optimal Bayesian quadrature rules are analogous to classical Gaussian quadrature rules in the sense that, in addition to being exact for kernel interpolants (recall (4)), they also exactly integrate Hermite interpolants (see Section 3.2.2). We thus begin by reviewing the argument used to establish positivity of the Gaussian quadrature weights.

3.1 Gaussian quadrature

Under the assumption that $\nu$ admits a density222This can be generalised to the cumulative distribution function having infinitely many points of increase. there exist unique weights $w_{1},\ldots,w_{n}$ and points $x_{1},\ldots,x_{n}\in[a,b]$ such that

[TABLE]

for every polynomial $P$ of degree at most $2n-1$ (Gautschi2004, , Chapter 1). This quadrature rule is known as a Gaussian quadrature rule (for the measure $\nu$). One can show the positivity of the weights of a Gaussian rule as follows.

Proposition 1

Assume that $\nu$ admits a Lebesgue density. Then the weights $w_{1},\dots,w_{n}$ of the Gaussian quadrature (7) are positive.

Proof

For each $i=1,\ldots,n$ there exists a unique polynomial $L_{i}$ of degree $n-1$ such that $L_{i}(x_{j})=\delta_{ij}$. This property is shared by the function $G_{i}\coloneqq L_{i}^{2}\geq 0$ that, being of degree $2n-2$, is also integrated exactly by the Gaussian rule. Because $G_{i}$ is almost everywhere positive, it follows from the assumption that $\nu$ admits a density that

[TABLE]

The positivity of the weights is thus concluded. ∎

This proof may appear to be based on the closedness of the set of polynomials under exponentiation. Closer analysis reveals a structure that can be later generalised.

To describe this, recall that one of the basic properties of polynomials is that a polynomial $P$ of degree $n$ can have at most $n$ zeroes, when counting multiplicities (for some properties of polynomials and interpolation with them, see e.g. (Atkinson1989, , Chapter 3)). This is to say that, if for some points $x_{1},\ldots,x_{m}$ it holds that

[TABLE]

for $j_{i}=0,\ldots,q_{i}-1$, with $q_{i}$ being the multiplicity of the zero $x_{i}$ of $P$, then $\sum_{i=1}^{m}q_{i}\leq n$. This fact on zeroes of polynomials can be used to supply a proof of positivity of the Gaussian quadrature weights that does not explicitly use of the fact that square of a function is non-negative. By the chain rule, the derivative of $G_{i}$ vanishes at each $x_{j}$ such that $j\neq i$. That is, $G_{i}$ has a double zero at each of these $n-1$ points (i.e., $G_{i}(x_{j})=0$ and $G_{i}^{(1)}(x_{j})=0$), for the total of $2n-2$ zeroes. Being a polynomial of degree $2n-2$, $G_{i}$ cannot have any other zeroes besides these. Since all the zeroes of $G_{i}$ are double, it cannot hence have any sign changes. This is because, in general, a function $g$ that satisfies $g(x)=g^{(1)}(x)=0$ but $g^{(2)}(x)\neq 0$ at a point $x$ cannot change its sign at $x$, since its derivative changes sign at the point. From $G_{i}(x_{i})=1>0$ it then follows that $G_{i}$ is almost everywhere positive.

3.2 Chebyshev systems and generalised Gaussian quadrature

The argument presented above works almost as such when the polynomials are replaced with generalised polynomials and the Gaussian quadrature rule with a generalised Gaussian quadrature rule. Much of the following material is covered by the introductory chapters of the monograph by Karlin and Studden KarlinStudden1966 . In the following $C^{m}([a,b])$ stands for the set of functions that are $m$ times continuously differentiable on the open interval $(a,b)$.

Definition 1 (Chebyshev system)

A collection of functions $\{\phi_{i}\}_{i=1}^{m}\subset C^{m-1}([a,b])$ constitutes an (extended) Chebyshev system if any non-trivial linear combination of the functions, called a generalised polynomial, has at most $m-1$ zeroes, counting multiplicities.

Remark 1

Some of the results we later present, such as Proposition 3, are valid even when a less restrictive definition, that does not require differentiability of $\phi_{i}$, of a Chebyshev system is used. Of course, in this case the definition is not given in terms of multiple zeroes. The above definition is used here to simplify presentation. The simplest relaxation is to require that $\{\phi_{i}\}_{i=1}^{m}$ are merely continuous and that no linear combination can vanish at more than $m-1$ points.

By selecting $\phi_{i}(x)=x^{i-1}$ we see that polynomials are an example of a Chebyshev system. Perhaps the simplest example of a non-trivial Chebyshev system is given by the following example.

Example 1

Let $\phi_{i}(x)=\mathrm{e}^{x}x^{i-1}$ for $i=1,\dots,m$. Then $\{\phi_{i}\}_{i=1}^{m}$ constitute a Chebyshev system. To verify this, observe that any linear combination $\phi$ of $\phi_{1},\ldots,\phi_{m}$ is of the form $\phi(x)=\mathrm{e}^{x}P(x)$ for a polynomial $P$ of degree at most $m-1$ and that the $j$th derivative of this function takes the form

[TABLE]

for certain integer coefficients $c_{1},\ldots,c_{j}$. We observe that $\phi(x_{0})=0$ for a point $x_{0}$ if and only if $P(x_{0})=0$. If also $\phi^{(1)}(x_{0})=0$, then it follows from (8) that $P^{(1)}(x_{0})=0$, and, generally, that $\phi^{(i)}(x_{0})=0$ for $i=0,\ldots,j$ if and only if $P^{(i)}(x_{0})=0$. That is, the zeroes of $\phi$ are precisely those of $P$ and, consequently, the functions $\phi_{i}$ constitute a Chebyshev system.

3.2.1 Interpolation using a Chebyshev system

A crucial property of generalised polynomials is that unique interpolants can be constructed using them, as we next show. For any Chebyshev system $\{\phi_{i}\}_{i=1}^{n}$ and a set of distinct points $X=\{x_{1},\ldots,x_{n}\}\subset[a,b]$, we know that there cannot exist $\boldsymbol{\alpha}=(\alpha_{1},\ldots,\alpha_{n})\neq\boldsymbol{0}$ such that

[TABLE]

since $\alpha_{1}\phi_{1}+\cdots+\alpha_{n}\phi_{n}$ can have at most $n-1$ zeroes. Equivalently, the only solution $\boldsymbol{\beta}\in\mathbb{R}^{n}$ to the linear system $\boldsymbol{V}_{X}^{\textsf{T}}\boldsymbol{\beta}=\boldsymbol{0}$ defined by the $n\times n$ matrix $[\boldsymbol{V}_{X}]_{ij}=\phi_{i}(x_{j})$ is $\boldsymbol{\beta}=\boldsymbol{0}$. That is, $\boldsymbol{V}_{X}$ is invertible.

For any data $\{(x_{i},f(x_{i}))\}_{i=1}^{n}$, the above fact guarantees the existence and uniqueness of an interpolant $s_{X,f}$ such that (i) $s_{X,f}$ is in $\mathrm{span}\{\phi_{1},\ldots,\phi_{n}\}$ and (ii) $s_{X,f}(x_{j})=f(x_{j})$ for each $j=1,\ldots,n$. These two requirements imply that

[TABLE]

for and some $\boldsymbol{\alpha}\in\mathbb{R}^{n}$ and every $j=1,\dots,n$. In matrix form, these $n$ equations are equivalent to $\boldsymbol{V}_{X}^{\textsf{T}}\boldsymbol{\alpha}=\boldsymbol{f}_{X}$. Hence $\boldsymbol{\alpha}=\boldsymbol{V}_{X}^{-\textsf{T}}\boldsymbol{f}_{X}$ and the interpolant is

[TABLE]

for $[\boldsymbol{\phi}(x)]_{i}=\phi_{i}(x)$ an $n$-vector.

3.2.2 Hermite interpolants

A Hermite interpolant $s_{X,\boldsymbol{q},f}$ is based on data containing also derivative values (see (Atkinson1989, , Section 3.6) for polynomial and (Fasshauer2007, , Chapter 36) for kernel-based Hermite interpolation). In this setting, the point set $X$ contains $m$ points and $\boldsymbol{q}\in\mathbb{N}_{0}^{m}$ is a vector of multiplicities such that $\sum_{i=1}^{m}q_{i}=n$. The data to be interpolated is

[TABLE]

That is, the interpolant is to satisfy

[TABLE]

for each $i=1,\ldots,m$ and $j_{i}=0,\ldots,q_{i}-1$. Note that the interpolant $s_{X,f}$ is a Hermite interpolant with $m=n$ and $q_{1}=\cdots=q_{n}=1$. If the interpolant is to lie in $\mathrm{span}\{\phi_{1},\ldots,\phi_{n}\}$, we must have, for some $\alpha_{1},\ldots,\alpha_{n}$,

[TABLE]

Again, these $n$ equations define a linear system that is invertible because $\{\phi_{i}\}_{i=1}^{n}$ constitute a Chebyshev system. The Hermite interpolant can be written in the form (9) with $\boldsymbol{V}_{X}$ replaced with a version involving also derivatives of $\phi_{i}$ (see e.g. (Oettershagen2017, , Section 2.3.1)).

3.2.3 Generalised Gaussian quadrature

A generalised Gaussian quadrature rule is a quadrature rule that uses $n$ points to integrate exactly all functions in the span of $\{\phi_{i}\}_{i=1}^{2n}$ constituting a Chebyshev system:

[TABLE]

for every $\phi\in\mathrm{span}\{\phi_{1},\ldots,\phi_{2n}\}$. The existence and uniqueness of the points and weights is guaranteed under fairly general assumptions Barrow1978 . We prove positivity of the weights by constructing a function ${F_{i}\in\mathrm{span}\{\phi_{1},\ldots,\phi_{2n}\}}$ analogous to $G_{i}$ in Section 3.1.

Proposition 2

Assume that $\nu$ admits a Lebesgue density. Then the weights $w_{1},\dots,w_{n}$ of the generalised Gaussian quadrature rule (10) are positive.

Proof

Let $F_{i}$ be the Hermite interpolant to the data

[TABLE]

An example is depicted in Figure 1. As there are $2n$ data points, $F_{i}$ indeed exists since $\{\phi_{i}\}_{i=1}^{2n}$ are a Chebyshev system. Moreover, $F_{i}$ has $2n-1$ zeroes. Because all its zeroes occurring on $(a,b)$ are double, $F_{i}$ cannot have sign-changes. Since $F_{i}(x_{i})=1>0$, we conclude $F_{i}$ is almost everywhere positive. Consequently, $w_{i}=I_{\nu}(F_{i})>0$. ∎

Next we turn our attention to kernels whose translates and their derivatives constitute Chebyshev systems.

3.3 Totally positive kernels

We are now ready to begin considering kernels and Bayesian quadrature. A concept related to Chebyshev systems is that of totally positive kernels whose theory is covered by the monograph of Karlin Karlin1968 . For a sufficiently differentiable kernel, define the derivatives

[TABLE]

If the derivative

[TABLE]

exists and is continuous for every $j\leq m$, the kernel is said to be $m$ times continuously differentiable, which we denote by writing $k\in C^{m}([a,b]^{2})$. In this case, $f\in C^{m}([a,b])$ if $f\in\mathcal{H}(k)$ and that the kernel derivatives (11) act as representers for differentiation (i.e., $\langle f,k^{(j)}(\cdot,x)\rangle_{\mathcal{H}(k)}=f^{(j)}(x)$ for $f\in\mathcal{H}(k)$ and $j\leq m$); see (Steinwart2008, , Corollary 4.36) and its proof.

Definition 2 (Totally positive kernel)

A kernel $k\in C^{\infty}([a,b]^{2})$ is (extended) totally positive of order $q\in\mathbb{N}$ if the collection

[TABLE]

constitutes a Chebyshev system for any $m\in\mathbb{N}$, any distinct $x_{1},\ldots,x_{m}\in\Omega$ and any multiplicities $q_{1},\ldots,q_{m}\leq q$ of these points.

The class of totally positive kernels is smaller than that of positive-definite kernels. For the simplest case of $q=1$ and $m=n$ the total positivity condition is that the kernel translates $k_{x_{1}},\ldots,k_{x_{n}}$ constitute a Chebyshev system. This implies that the $n\times n$ matrix $[\boldsymbol{K}_{Y,X}]\coloneqq k(y_{j},x_{i})$, which is just the matrix $\boldsymbol{V}_{Y}$ considered in Section 3.2 for the Chebyshev system $\phi_{i}=k_{x_{i}}$, is invertible for any ${Y=\{y_{1},\ldots,y_{n}\}\subset[a,b]}$. Positive-definiteness of $k$ only guarantees that $\boldsymbol{K}_{Y,X}$ is invertible when $Y=X$.

Basic examples of totally positive kernels are the Gaussian kernel

[TABLE]

with length-scale $\ell>0$ and the Hardy kernel ${k(x,x^{\prime})=r^{2}/(r^{2}-xx^{\prime})}$ for $r>0$. Both of these kernels are totally positive of any order. There is also a convenient result that guarantees total positivity (Burbea1976, , Proposition 3): $k$ is totally positive if there are positive constants $a_{m}$ and a positive increasing function $v\in C^{\infty}([a,b])$ such that

[TABLE]

for all $x,x^{\prime}\in\Omega$. More examples are collected in Karlin1968 ; Burbea1976 .

3.4 General result on weights

The following special case of the theory developed in (KarlinStudden1966, , Chapter 2) appears in, for instance, (RichterDyn1971a, , Lemma 2). Its proof is a generalisation of the proof for the case $m=2n$ that was discussed in Section 3.2.

Proposition 3

Suppose that $\{\phi_{i}\}_{i=1}^{m}\subset C^{m-1}([a,b])$ constitute a Chebyshev system, that $\nu$ admits a Lebesgue density and that $Q(f)\coloneqq\sum_{i=1}^{n}w_{i}f(x_{i})$ for $x_{1},\ldots,x_{m}\in\Omega$ is a quadrature rule such that ${Q(\phi_{i})=I_{\nu}(\phi_{i})}$ for each $i=1,\ldots,m$. Then at least $\lfloor(m+1)/2\rfloor$ of the weights $w_{1},\ldots,w_{n}$ are positive.

An immediate consequence of this proposition is that a Bayesian quadrature rule based on a totally positive kernel has at least one half of its weights positive.

Theorem 3.1

Suppose that the kernel $k\in C^{\infty}([a,b]^{2})$ is totally positive of order 1. Then, for any points, at least $\lfloor(n+1)/2\rfloor$ of the Bayesian quadrature weights $w_{X,1}^{\text{{\tiny BQ}}},\ldots,w_{X,n}^{\text{{\tiny BQ}}}$ are positive.

Proof

Since the kernel is totally positive of order 1, the translates $\{k_{x_{i}}\}_{i=1}^{n}$ constitute a Chebyshev system. The exactness condition (4) holds for each of these functions. The claim follows by setting $m=n$ in Proposition 3. ∎

3.5 Weights for locally optimal points

Recall the definition of the Bayesian quadrature variance:

[TABLE]

The variance can be considered a function $X\mapsto\mathbb{V}_{X}^{\text{{\tiny BQ}}}$ defined on the simplex

[TABLE]

We introduce the following definition of locally optimal points. For this purpose, define the function

[TABLE]

and its partial derivatives

[TABLE]

Definition 3

Let $m\leq n$. A Bayesian quadrature rule with points $X\subset[a,b]$ is locally $m$-optimal if $X\in\mathcal{S}^{n}$ and there is an index set $\mathcal{I}_{m}^{*}\subset\{1,\ldots,n\}$ of $m$ indices such that

[TABLE]

A locally $n$-optimal rule is called locally optimal. The point set of a locally $m$-optimal Bayesian quadrature rule is also called locally $m$-optimal.

When the kernel is totally positive of any order, it has been shown that any local minimiser of $\mathbb{V}_{X}^{\text{{\tiny BQ}}}$ is locally optimal in the sense of above definition. That is, no point in a point set that locally minimises the variance can be located on the boundary of the integration interval nor can any two points in the set coalesce.333Coalescence is possible because $\mathbb{V}_{X}^{\text{{\tiny BQ}}}$ is in fact a continuous function of $X$ defined on the whole of $\Omega^{n}$, not merely on $\mathcal{S}^{n}$ (Oettershagen2017, , Proposition 5.5). Coalescence of some of the points would result in a quadrature rule that uses also evaluations of derivatives of the integrand. These results, the origins of which can be traced to the 1970s BarrarLoebWerner1974 ; BarrarLoeb1976 ; Bojanov1979 , have been recently collated by Oettershagen (Oettershagen2017, , Corollary 5.13).

A locally $m$-optimal Bayesian quadrature rule is, in addition to the kernel translates at $X$, exact for translate derivatives at $x_{j}$ with $j\in\mathcal{I}_{m}^{*}$ (it is worth noting that Bayesian quadrature rules with derivative evaluations have been recently considered in PruherSarkka2016 ; WuAoiPillow2018 ). When $m=n$, this is analogous to the interpretation of classical Gaussian quadrature rules as integrated Hermite interpolants RichterDyn1971b . This result first appeared in Larkin1970 . Its proof is typically based on considering the RKHS representation

[TABLE]

of the variance; see (RichterDyn1971a, , Section 3) or (Oettershagen2017, , Section 5.1.3). We present a mainly linear algebraic proof.

Proposition 4

Let $m\leq n$. Suppose that the $n$-point set $X\in\mathcal{S}^{n}$ is locally $m$-optimal. If the kernel $k$ is once continuously differentiable, then

[TABLE]

where $k_{x}^{(1)}$ is the kernel derivative defined in (11).

Proof

By definition of local $m$-optimality, the partial derivatives

[TABLE]

must vanish for each $j\in\mathcal{I}_{m}^{*}$. Let $\partial_{i}\boldsymbol{g}(X)\in\mathbb{R}^{n}$ stand for the $i$th partial derivative of a vector-valued function $\boldsymbol{g}\colon\mathbb{R}^{n}\to\mathbb{R}^{n}$ evaluated at $X$. From the explicit expression (2) for the variance we compute

[TABLE]

where the inverse matrix derivative formula

[TABLE]

and the weight expression $\boldsymbol{w}_{X}^{\text{{\tiny BQ}}}=\boldsymbol{K}_{X}^{-1}\boldsymbol{k}_{\nu,X}$ have been used. The two partial derivatives appearing in the equation for $E_{j}(X)$ can be explicitly computed. First, only the $j$th element of $\boldsymbol{k}_{\nu,X}$ depends on $x_{j}$. Thus,

[TABLE]

Secondly, only the $j$th row and column of $\boldsymbol{K}_{X}$ have dependency on $x_{j}$. For $l\neq j$ we have

[TABLE]

where the first equality is consequence of symmetry of the kernel. The diagonal element is a total derivative:

[TABLE]

Therefore $\partial_{j}\boldsymbol{K}_{X}$ is a zero matrix except for the $j$th row and column that are

[TABLE]

and its transpose, respectively. Hence

[TABLE]

If $w_{X,j}^{\text{{\tiny BQ}}}\neq 0$, then $E_{j}(X)=0$ so that the form of $E_{j}$ above implies that $I_{X}^{\text{{\tiny BQ}}}(k_{x_{j}}^{(1)})=I_{\nu}(k_{x_{j}}^{(1)})$. This concludes the proof. ∎

Remark 2

Proposition 4 admits an obvious multivariate extension (Gavrilov1998, , Теорема 2): when $d>1$, the $md$ partial derivative representers

[TABLE]

for $j=1,\ldots,d$ and $i\in\mathcal{I}_{m}^{*}$ are integrated exactly by a locally $m$-optimal Bayesian quadrature rule, defined by requiring a gradient version of (13). See also Gavrilov2007 . However, there appear to exist no generalisations of Chebyshev systems and Proposition 3 to higher dimensions.

Theorem 3.2

Let $k\in C^{\infty}([a,b]^{2})$ be a totally positive kernel of order 2 and $m\leq n$. Suppose that the point set $X\in\mathcal{S}^{n}$ is locally $m$-optimal with an index set $\mathcal{I}_{m}^{*}\subset\{1,\ldots,n\}$ and that the weights associated with $q\leq m$ indices in $\mathcal{I}_{m}^{*}$ are non-zero. Then at least $\lfloor(n+2m-q+1)/2\rfloor$ of the weights are non-negative, and $q$ must satisfy $2m-n\leq q$.

Proof

By (14), the Bayesian quadrature rule in the statement is exact for $n$ kernel translates and $q$ of their derivatives. By the total positivity of the kernel, the collection of these $n+q$ functions constitutes a Chebyshev system. By Proposition 3, at least $\lfloor(n+q+1)/2\rfloor$ of the weights are positive. Since the weights associated with $m-q$ indices in $\mathcal{I}_{m}^{*}$ are zero, it follows that at least $\lfloor(n+q+1)/2\rfloor+m-q=\lfloor(n+2m-q+1)/2\rfloor$ of the weights are non-negative. The lower-bound for $q$ follows because $\lfloor(n+2m-q+1)/2\rfloor\leq n$ implies that $n+2m-q+1\leq 2n+1$. ∎

The main result of this section follows by setting $m=n$ in the preceding theorem and observing that this implies $q=n$, which means that there can be no zero weights.

Corollary 1

If $k\in C^{\infty}([a,b]^{2})$ is totally positive of order 2 and $X\in\mathcal{S}^{n}$ is locally optimal, then all the Bayesian quadrature weights $w_{X,1}^{\text{{\tiny BQ}}},\ldots,w_{X,n}^{\text{{\tiny BQ}}}$ are positive.

Remark 3

A key consequence of Corollary 1 is the following: If $w_{X,1}^{\text{{\tiny BQ}}},\ldots,w_{X,n}^{\text{{\tiny BQ}}}$ contain negative values, then the design points $X$ are not locally optimal. In other words, in this case there is still room for improvement by optimising these points using, for example, gradient descent. In this way, the signs of the weights can provide information about the quality of the design point set.

A positive-weight quadrature rule is a positive linear functional (i.e., every positive function is mapped to a positive real). A locally optimal Bayesian quadrature rule may therefore be appropriate for numerical integration of functions that are a priori known to be positive, such as likelihood functions. Theoretical comparison to warped models Osborne2012 ; Gunter2014 ; Chai2018 that encode positivity of the integrand by placing the GP prior on, for example, square root of the integrand would be an interesting topic of research.

3.6 Greedily selected points

The optimal points discussed in the preceding section cannot be constructed efficiently (see (Oettershagen2017, , Section 5.2) for what appears to be the most advanced published algorithm). In practice, points selected by greedy minimisation of the integral variance are often used. This approach is known as sequential Bayesian quadrature CookClayton1998 ; HuszarDuvenaud2012 . Assuming for a moment that $d$ is arbitrary and an $n$-point set $X_{n}\subset\Omega$ has been already generated, sequential Bayesian quadrature proceeds by selecting a new point $\boldsymbol{x}_{n+1}\in\Omega$ by minimising the integral variance:

[TABLE]

In higher dimensions there is little that we are able to say about qualitative properties of the resulting quadrature rules. However, when $d=1$ we can invoke Theorem 3.2 since $X_{n}\cup\{x_{n+1}\}$ is locally $1$-optimal.

Proposition 5

Suppose that $k\in C^{\infty}([a,b]^{2})$ is totally positive of order 2. If $X_{n}\cup x_{n+1}\in\mathcal{S}^{n}$, then at least $\lfloor(n+3)/2\rfloor$ of the weights of a $n+1$ point sequential Bayesian quadrature rule are positive.

3.7 Other kernels and point sets

A number of combinations of kernels and point sets, that are not covered by the theory above, have been shown, either theoretically or experimentally, to yield positive Bayesian quadrature weights:

• •

The GP posterior mean for the Brownian motion kernel $k(x,x^{\prime})=\min(x,x^{\prime})$ on $[0,1]$ is a piecewise linear interpolant. As this implies that the Lagrange cardinal functions $u_{X,i}$ are non-negative, it follows from the identity $w_{X,i}^{\text{{\tiny BQ}}}=I_{\nu}(u_{X,i})$ that the weights are positive. See Diaconis1988 and (Ritter2000, , Lemma 8 in Section 3.2, Chapter 2) for more discussion.

• •

Suitably selected priors give rise to Bayesian quadrature rules whose posterior mean coincides with a classical rule, such a Gaussian quadrature Karvonen2017a ; Karvonen2018c . Analysis of the weights and their positivity naturally reduces to that of the reproduced classical rule.

• •

There is convincing numerical evidence that the weights are positive if the nodes for the Gaussian kernel and measure on $\mathbb{R}$ are selected by suitable scaling the classical Gauss–Hermite nodes Karvonen2018b .

• •

Uniform weighting (i.e., $w_{X,i}^{\text{{\tiny BQ}}}=1/n$) can be achieved when certain quasi-Monte Carlo point sets and shift invariant kernels are used Jagadeeswaran2018 .

3.8 Upper bound on the sum of weights

We summarise below a simple yet generic result that has an important consequence on the stability of Bayesian quadrature in Section 4.

Lemma 1

Let $\Omega\subset\mathbb{R}^{d}$. If the Bayesian quadrature weights $w_{X,1}^{\text{{\tiny BQ}}},\dots,w_{X,n}^{\text{{\tiny BQ}}}$ are non-negative, then we have

[TABLE]

Proof

The claim immediately follows from the property (4) that $\sum_{i=1}^{n}w_{X,i}^{\text{{\tiny BQ}}}k_{\boldsymbol{x}_{j}}(\boldsymbol{x}_{i})=I_{\nu}(k_{\boldsymbol{x}_{j}})$ for each $j=1,\ldots,n$. ∎

Combined with Corollary 1, we get a bound on the sum of absolute weights $\sum_{i=1}^{n}|w_{X_{n},i}^{\text{{\tiny BQ}}}|$, which is the main topic of discussion in Section 4.

Corollary 2

Let $\Omega=[a,b]\subset\mathbb{R}$. If $k\in C^{\infty}([a,b]^{2})$ is totally positive of order 2 and design points $X\in\mathcal{S}^{n}$ are locally optimal, then we have

[TABLE]

Most importantly, Corollary 2 is applicable to the Gaussian kernel, for which the upper-bound is finite. This result will be discussed in Section 4.4 in more detail. One may see a supporting evidence in Fig. 2, where the sum of weights seems to converge to a value around $1$.

3.9 Higher dimensions

As far as we are aware of, there are no extensions of the theory of Chebyshev systems to higher dimensions. Consequently, it is not possible to say much about positivity of the weights when $d>1$. Some simple cases can be analysed, however.

Let $\Omega_{1}=[a,b]$, $\nu_{1}$ be a measure on $\Omega_{1}$, $\Omega=\Omega_{1}^{d}\subset\mathbb{R}^{d}$ and $\nu=\nu_{1}^{d}$. That is, $\Omega=\Omega_{1}\times\cdots\times\Omega_{1}$ and $\operatorname{d\!}\nu(\boldsymbol{x})=\operatorname{d\!}\nu_{1}(x_{1})\times\cdots\times\operatorname{d\!}\nu_{1}(x_{d})$, where there are $d$ terms in the products. Suppose that

1. (i)

the point set $X$ is now a Cartesian product of one-dimensional sets $X_{1}=\{x_{1}^{1},\ldots,x_{n}^{1}\}\subset\Omega_{1}$: $X=X_{1}^{d}$; 2. (ii)

the kernel is of product form: $k(\boldsymbol{x},\boldsymbol{x}^{\prime})=\prod_{i=1}^{d}k_{1}(x_{i},x_{i}^{\prime})$ for some kernel $k_{1}$ on $\Omega_{1}$.

A quadrature rule using Cartesian product points is called a tensor product rule. For such points the Bayesian quadrature weights $\boldsymbol{w}_{X}^{\text{{\tiny BQ}}}$ are products of the one-dimensional weights $\boldsymbol{w}_{X_{1}}^{\text{{\tiny BQ}}}$: the weight for the point $(x_{i(1)},\ldots,x_{i(d)})\in X$ is $\prod_{j=1}^{d}w_{X_{1},i(j)}^{\text{{\tiny BQ}}}$ (Oettershagen2017, , Section 2.4). In particular, if $k_{1}$ is totally positive and $X_{1}$ is a locally optimal set of points, then all the $n^{d}$ weights $\boldsymbol{w}_{X}^{\text{{\tiny BQ}}}$ are positive.444Note that a tensor product rule based on an optimal one-dimensional point set need not be locally optimal for $\Omega$, $\nu$ and $k$. Analysis of more flexible sparse grid and symmetry based methods Karvonen2018 might yield more interesting results.

We conclude this section with two numerical examples. Both of them involve the standard Gaussian measure

[TABLE]

on $\Omega=\mathbb{R}^{d}$ and the Gaussian kernel (12).

Locally optimal points

First, we investigated positivity of weights for locally optimal points. We set $\ell=1$ and $d=2$ and used a gradient-based quasi-Newton optimisation method (MATLAB’s fminunc) to find points that locally minimise the integral variance for $n=2,\ldots,20$. Optimisation was initialised with a set of random points. The point set output by the optimiser was then randomly perturbated and optimisation repeated for 20 times, each time initialising with the point set giving the smallest Bayesian quadrature variance so far. The weights were always computed directly from (3). However, to improve numerical stability, the kernel matrix $\boldsymbol{K}_{X}$ was replaced by $\boldsymbol{K}_{X}+10^{-6}\boldsymbol{I}$, where $\boldsymbol{I}$ is the $n\times n$ identity matrix, during point optimisation. Some point sets generated using the same algorithm have appeared in (Sarkka2016, , Section IV) (for other examples of optimal points in dimension two, see OHagan1992 ; Minka2000 and, in particular, (Oettershagen2017, , Chapter 6)). The point sets we obtained appear sensible and all of them are associated with positive weights; four sets and their weights are depicted in Fig. 2. For $n=20$, the maximal value of a partial derivative of $\mathbb{V}_{X}^{\text{{\tiny BQ}}}$ at the computed points was $9\times 10^{-10}$.

Random points

Secondly, we investigated the validity of Theorem 3.1 in higher dimensions. We set $\ell=1.5$ and $d=4$ and counted the number of positive weights for $n=2,\ldots,1000$ when each $n$-point set is generated by drawing Monte Carlo samples from $\nu$. Random samples are often used in Bayesian quadrature RasmussenGhahramani2002 ; Briol2017 ; Briol2017a and they also function as a suitable test case where structurality of point sets has little role in constraining behaviour of some subsets of the weights as happens when product or symmetric point designs are used. Fig. 3 shows the proportion of positive weights; it appears that at least half of the weights for randomly drawn points are always positive. This supports the obvious conjectural extension to higher dimensions of Theorem 3.1.

4 Magnitudes of weights and the stability

This section studies the magnitudes of the weights in a Bayesian quadrature rule and discusses how they are related to stability and robustness of the quadrature rule. We are in particular interested in the following quantity, which we call the Bayesian quadrature stability constant:

[TABLE]

To make dependency on $n$ more explicit, the quadrature point set is denoted by $X_{n}$ instead of $X$ in this section. The terminology is motivated by the close connection of $\Lambda_{X_{n}}^{\text{{\tiny BQ}}}$ to the Lebesgue constant $\Lambda_{X_{n}}$, a quantity that characterises the stability of an interpolant. For kernel interpolants the Lebesgue constant is

[TABLE]

where $u_{X_{n},i}$ are Lagrange cardinal functions from Section 2.1. The connection to (15) arises from the fact that $w_{X_{n},i}^{\text{{\tiny BQ}}}=I_{\nu}(u_{X_{n},i})$ for $i=1,\dots,n$.

The importance of the stability constant (15) is illustrated by the following argument. Let $\mu_{f}^{*}$ be an optimal approximant to the integrand function $f\colon\Omega\to\mathbb{R}$ in the span of $\{k_{\boldsymbol{x}_{i}}\}_{i=1}^{n}$ in the sense that

[TABLE]

where $\mathinner{\lVert f-\mu_{f}\rVert}_{\infty}\coloneqq\sup_{\boldsymbol{x}\in\Omega}\mathinner{\lvert f(\boldsymbol{x})-\mu_{f}(\boldsymbol{x})\rvert}$ is the uniform norm. Note that $\mu_{f}^{*}$ does not in general interpolate $f$ at $X_{n}$ nor coincide with the Gaussian process posterior mean $\mu_{X,f}$. Then

[TABLE]

where we have used the fact that $I_{X_{n}}^{\text{{\tiny BQ}}}(g)=I_{\nu}(g)$ if $g\in\mathrm{span}\{k_{\boldsymbol{x}_{i}}\}_{i=1}^{n}$. That is, the approximation error by a Bayesian quadrature rule can be related to that by the best uniform approximant via the stability constant.

The stability constant also controls the error introduced by inaccurate funtion evaluations. Suppose that the function evaluations contain errors (which may be numerical or stochastic), denoted by $\epsilon_{i}$ and modelled as indendent zero-mean random variables with variance $\sigma^{2}$. Then the mean-square error (where the expectation is w.r.t. $\epsilon_{1},\dots,\epsilon_{n}$) of Bayesian quadrature is given by

[TABLE]

This implies a small stability constant (15) suppresses the additional error caused by the perturbations $\epsilon_{i}$. A third motivating example will be given in Section 4.2, after introducing necessary notation.

It is clear from Lemma 1 that if the weights are positive for every $n$, the stability constant remains uniformly bounded. However, the results on positivity in the preceding section are valid only when $d=1$ and the kernel is totally positive. This section uses a different technique to analyse the stability constant. The results are based on those in DeMarchiSchaback2010 , which are applicable to kernels that induce Sobolev-equivalent RKHSs (e.g., Matérn kernels). Accordingly, we mainly focus on such kernels in this section. We begin by reviewing basic properties of Sobolev spaces in Section 4.1 and convergence results for Bayesian quadrature in Section 4.2. The main results, Theorem 4.3 and Corollary 3, on the magnitudes of quadrature weights and the stability constant appear in Section 4.3. We discuss a relevant stability issue, known as the Runge phenomenon, for infinitely smooth kernels such as the Gaussian kernel in Section 4.4. Finally, simulation results in Section 4.5 demonstrate that the obtained upper bound is conservative; there is much room for improving the results.

Notation and basic definitions.

The Fourier transform $\hat{f}$ of a Lebesgue integrable $f\colon\mathbb{R}^{d}\to\mathbb{R}$ is defined by

[TABLE]

Two normed vector spaces $\mathcal{F}_{1}$ and $\mathcal{F}_{2}$ are norm-equivalent if $\mathcal{F}_{1}=\mathcal{F}_{2}$ as a set and there exist constants $C_{1},C_{2}>0$ such that

[TABLE]

4.1 Kernels inducing Sobolev-equivalent RKHSs

Let $\Phi\colon\mathbb{R}^{d}\to\mathbb{R}$ be a continuous and integrable positive-definite function with Fourier transform satisfying

[TABLE]

for $r>d/2$, some positive constants $c_{1}$ and $c_{2}$, and for all $\boldsymbol{\xi}\in\mathbb{R}^{d}$. In this section, we consider shift-invariant kernels on $\mathbb{R}^{d}$ of the form $k(\boldsymbol{x},\boldsymbol{x}^{\prime})=\Phi(\boldsymbol{x}-\boldsymbol{x}^{\prime})$. For instance, a Matérn kernel (Rasmussen2006, , Section 4.2.1)

[TABLE]

with smoothness parameter $\rho\mathrel{\mathop{\ordinarycolon}}=r-d/2$ and length-scale parameter $\ell>0$ satisfies (16).555Note that the smoothness parametrisation $\rho=r$ is often used. With this parametrisation $k_{\rho}$ would satisfy (16) with the exponent $-(r+d/2)$ and its RKHS would be norm-equivalent to $H^{r+d/2}(\mathbb{R}^{d})$. Here $\mathrm{K}_{\rho}$ is the modified Bessel function of the second kind of order $\rho$. Another notable class of kernels satisfying (16) are Wendland kernels (Wendland2005, , Theorem 10.35).

By (Wendland2005, , Corollary 10.13), the RKHS $\mathcal{H}(k)$ of any kernel $k$ satisfying (16) is norm-equivalent to the Sobolev space $H^{r}(\mathbb{R}^{d})$ of order $r>d/2$ on $\mathbb{R}^{d}$, which is a Hilbert space consisting of square-integrable and continuous functions $f\colon\mathbb{R}^{d}\to\mathbb{R}$ such that

[TABLE]

As can be seen from this expression, $r$ quantifies the smoothness of functions in $H^{r}(\mathbb{R}^{d})$: as $r$ increases, function in $H^{r}(\mathbb{R}^{d})$ become smoother.

The Sobolev space $H^{r}(\Omega)$ on a general measurable domain $\Omega\subset\mathbb{R}^{d}$ can be defined as the restriction of $H^{r}(\mathbb{R}^{d})$ onto $\Omega$. The kernel $k$ satisfying (16), when seen as a kernel on $\Omega$, then induces an RKHS that is norm-equivalent to $H^{r}(\Omega)$ (Wendland2005, , Theorems 10.12, 10.46 and 10.47).666The reader may ask whether $\Omega$ needs to have a Lipschitz boundary for this norm-equivalence, but this assumption is indeed not needed. The assumption that $\Omega$ has a Lipschitz boundary is required when using Stein’s extension theorem (Ste70, , p. 181) for Sobolev spaces defined using weak derivatives (see the proof of (Wendland2005, , Corollary 10.48)). On the other hand, we consider here a Sobolev space defined in terms of the Fourier transform, and the norm-equivalence follows from the extension and restriction theorems for a generic RKHS (Wendland2005, , Theorems 10.46 and 10.47) and the expression of the RKHS norm in terms of Fourier transforms (Wendland2005, , Theorem 10.12).

4.2 Convergence for Sobolev-equivalent kernels

Recall from Section 2.2 that integration error by a Bayesian quadrature rule for functions in $\mathcal{H}(k)$ satisfies

[TABLE]

so that in convergence analysis only the behaviour of the worst-case error needs to be considered. If the RKHS is norm-equivalent to a Sobolev space, rates of convergence for Bayesian quadrature can be established. These results follow from (Arcangeli2007, , Corollary 4.1). See (Wendland2005, , Corollary 11.33) or (WendlandRieger2005, , Proposition 3.6) for earlier and slightly more restricted results that require $\lfloor r\rfloor>d/2$ and (Kanagawa2018, , Proposition 4) for a version specifically for numerical integration. Some assumptions, satisfied by all domains of interest to us, are needed; see for instance (Kanagawa2018, , Section 3) for precise definitions.

Assumption 4.1

The set $\Omega\subset\mathbb{R}^{d}$ is a bounded open set that satisfies an interior cone condition and has a Lipschitz boundary.

This assumption essentially says that the boundary of $\Omega$ is sufficiently regular (Lipschitz boundary) and that there is no “pinch point” on the boundary of $\Omega$ (interior cone condition). Convergence results are expressed in terms of the fill-distance

[TABLE]

that quantifies the size of the largest “hole” in an $n$-point set $X_{n}$. We use $\lesssim$ to denote an inequality that is valid up to a constant independent of $n$, number of points, and $f$, the integrand. That is, for generic sequences of functionals $g_{n}$ and $h_{n}$, $g_{n}(f)\lesssim h_{n}(f)$ means that there is a constant $C>0$ such that $g_{n}(f)\leq Ch_{n}(f)$ for all $n\in\mathbb{N}$ and any $f$ in a specified class of functions.

Theorem 4.2

Suppose that (i) $\Omega$ satisfies Assumption 4.1 (ii) that the measure $\nu$ has a bounded (Lebesgue) density function and that (iii) the kernel $k$ satisfies (16) for a constant $r$ such that $r>d/2$. Then

[TABLE]

for any $f\in H^{r}(\Omega)$ when the fill-distance is sufficiently small.

The following simple result is an immediate consequence of this theorem.

Proposition 6

Suppose that the assumptions of Theorem 4.2 are satisfied. Then $\mathinner{\lvert 1-\sum_{i=1}^{n}w_{X_{n},i}^{\text{{\tiny BQ}}}\rvert}\lesssim h_{X_{n},\Omega}^{r}$ when the fill-distance is sufficiently small.

Proof

Under the assumptions, constant functions are in $H^{r}(\Omega)$. Setting $f\equiv 1$ in (17) verifies the claim. ∎

Note that the same argument can be used whenever a general rate of convergence for functions in an RKHS is known and constant functions are contained in the RKHS. However, this is not always the case; for example, the RKHS of the Gaussian kernel (12) does not contain polynomials (Minh2010, , Theorem 2).

Rates explicitly dependent on the number of points are achieved for point sets that are quasi-uniform, which is to say that

[TABLE]

for some constants $\tilde{c}_{1},\tilde{c}_{2}>0$ independent of $n$. Here

[TABLE]

is the separation distance. In dimension $d$ quasi-uniform sets satisfy $h_{X_{n},\Omega}=\mathcal{O}(n^{-1/d})$ as $n\to\infty$ (e.g., regular product grids). In Theorem 4.2 we thus obtain the rate

[TABLE]

for $f\in H^{r}(\Omega)$ when the point sets are quasi-uniform and $n$ is sufficiently large.

Of course, it is the stability constant $\Lambda_{X_{n}}^{\text{{\tiny BQ}}}=\sum_{i=1}^{n}\mathinner{\lvert w_{X_{n},i}^{\text{{\tiny BQ}}}\rvert}$ that we analyse next whose behaviour is typically more consequential. However, the above proposition may be occasionally interesting if one desires to interpret Bayesian quadrature as a weighted Dirac approximation $\nu_{\text{{\tiny BQ}}}\coloneqq\sum_{i=1}^{n}w_{X_{n},i}^{\text{{\tiny BQ}}}\delta_{\boldsymbol{x}_{i}}\approx\nu$ of a probability measure (i.e., $\nu_{\text{{\tiny BQ}}}(\Omega)\approx 1$). Note that there is also a simple way to ensure summing up to one of the weights by inclusion of a non-zero prior mean function for the Gaussian process prior; see OHagan1991 and (Karvonen2018c, , Section 2.3)

Finally, we provide a third example that highlights the importance of analysing the stability constant. Kanagawa et al. (Kanagawa2018, , Section 4.1) (see also Kanagawa2016 ) studied convergence rates of kernel-based quadrature rules in Sobolev spaces when the integrand is potentially rougher (i.e., $f\in H^{s}(\Omega)$ for some $s\leq r$) than assumed. If $s<r$, the integrand $f$ may not belong to the Sobolev space $H^{r}(\Omega)$ that is assumed by the user when constructing the quadrature rule; therefore this is a misspecifed setting. Under certain conditions, they showed (Kanagawa2018, , Corollary 7) that if $\Lambda_{X_{n}}^{\text{{\tiny BQ}}}\lesssim n^{c}$ for a constant $c\geq 0$, then

[TABLE]

when $X_{n}$ are quasi-uniform.

The condition $\Lambda_{X_{n}}^{\text{{\tiny BQ}}}\lesssim n^{c}$ means that the stability constant $\Lambda_{X_{n}}^{\text{{\tiny BQ}}}$ should not grow quickly as $n$ increases. The bound (18) shows that the error in the misspecified setting becomes small if $c$ is small. This implies that if the stability constant $\Lambda_{X_{n}}^{\text{{\tiny BQ}}}$ does not increase quickly, then the quadrature rule becomes robust against the misspecification of a prior. This provides a third motivation for understanding the behaviour of $\Lambda_{X_{n}}^{\text{{\tiny BQ}}}$.

4.3 Upper bounds for absolute weights

We now analyse magnitudes of individual weights and the stability constant (15). We first derive an upper bound on the magnitude of each weight $w_{X_{n},i}^{\text{{\tiny BQ}}}$. The proof of this result is based on an upper bound on the $L^{2}(\Omega)$ norm of Lagrange functions derived in DeMarchiSchaback2010 .

Theorem 4.3

Suppose that (i) $\Omega$ satisfies Assumption 4.1, that (ii) the measure $\nu$ has a bounded (Lebesgue) density function and that (iii) the kernel $k$ satisfies (16) for a constant $r$ such that $r>d/2$. Then

[TABLE]

for all $i=1,\dots,n$, provided that $h_{X_{n},\Omega}$ is sufficiently small. When $X_{n}$ are quasi-uniform, this becomes

[TABLE]

for $n$ large enough.

Proof

It is proved in (DeMarchiSchaback2010, , Theorem 1) that each of the Lagrange functions $u_{X_{n},i}$ admits the bound

[TABLE]

provided that $h_{X_{n},\Omega}$ is sufficiently small. Let $\mathinner{\lVert\nu\rVert}_{\infty}<\infty$ stand for the supremum of the density function of $\nu$. Then it follows from $w_{X_{n},i}^{\text{{\tiny BQ}}}=I_{\nu}(u_{X_{n},i})$ that

[TABLE]

Inequality (19) now follows from (21). When $X_{n}$ are quasi-uniform, the ratio $h_{X_{n},\Omega}/q_{X_{n}}$ remains bounded and $h_{X_{n},\Omega}$ behaves like $n^{-1/d}$. ∎

An important consequence of Theorem 4.3 is that the magnitudes of quadrature weights decrease uniformly to zero as $n$ increases if the design points are quasi-uniform and $\nu$ has a density. In other words, none of the design points will have a constant weight that does not decay. This is similar to importance sampling, where the weights decay uniformly at rate $1/n$. As a direct corollary of Theorem 4.3 we obtain bounds on the stability constant $\Lambda_{X_{n}}^{\text{{\tiny BQ}}}$.

Corollary 3

Under the assumptions of Theorem 4.3 and provided that $h_{X_{n},\Omega}$ is sufficiently small we have

[TABLE]

When $X_{n}$ are quasi-uniform and $n$ sufficiently large this becomes

[TABLE]

While the bounds of Corollary 3 are somewhat conservative (as will be demonstrated in Section 4.5), they are still useful in understanding the factors affecting stability and robustness of Bayesian quadrature. That is, inequality (22) shows that the stability constant can be made small if the ratio $h_{X_{n},\Omega}/q_{X_{n}}$ is kept small; this is possible if the point set is sufficiently uniform.

Another important observation concerns the exponent ${r-d/2}$ of the ratio $h_{X_{n},\Omega}/q_{X_{n}}$: if the smoothness $r$ of the kernel is large, then the stability constant may also become large if the points are not quasi-uniform. This is true because $h_{X_{n},\Omega}/q_{X_{n}}\geq 1$ for any configuration of $X_{n}$, as can be seen easily from the definitions of $q_{X_{n}}$ and $h_{X_{n},\Omega}$. This observation implies that the use of a smoother kernel may lead to higher numerical instability. Accordingly, we next discuss stability of infinitely smooth kernels and the Runge phenomenon that manifests itself in this setting.

4.4 On infinitely smooth kernels

While the theoretical results of this section only concern kernels of finite smoothness, we make a few remarks on the stability of Bayesian quadrature when using infinitely smooth kernels, such as the Gaussian kernel. When using such a kernel, Bayesian quadrature rules suffer from the famous Runge phenomenon: if equispaced points are used, then Lebesgue constants and the stability constants grow rapidly; see  (Oettershagen2017, , Section 4.3) and Platte2005 ; Platte2011 . This effect is demonstrated in Fig. 4, and can be seen also in (SommarivaVianello2006b, , Table 1).

A key point is that Runge phenomenon typically occurs when the design points are quasi-uniform (e.g., equispaced). This means that quasi-uniformity of the points does not ensure stability of Bayesian quadrature when the kernel is infinitely smooth. Care has to be taken if a numerically stable Bayesian quadrature rule is to be constructed with such a kernel. One possibility is to use locally optimal design points from Section 3.5. Corollary 2 then guarantees uniform boundedness of the stability constant, at least when $d=1$.

4.5 A numerical example

Numerical examples of the behaviour of kernel Lebesgue constants can be found in DeMarchiSchaback2008 , where it was observed that the theoretical bounds similar to (23) are conservative: the Lebesgue constant appears to remain uniformly bounded. Bayesian quadrature weights are no different. We experimented with the Matérn kernel

[TABLE]

with length-scale $\ell=0.5$ and the uniform measure on the interval $[0,1]$. When uniformly spaced points were used, all weights remained positive and their sum quickly converged to one when $n$ was increased. In contrast, Corollary 3 provides the, up to a constant, upper bound $\sqrt{n}$ that is in this case clearly very conservative. When points were drawn from the uniform distribution on $[0,1]$, more interesting behaviour was observed (Fig. 5). As expected, the magnitude of $\Lambda_{X_{n}}^{\text{{\tiny BQ}}}$ was closely related to the ratio $h_{X_{n},\Omega}/q_{X_{n}}$. Nevertheless, increase in $n$ did not generally correspond to increase in $\Lambda_{X_{n}}^{\text{{\tiny BQ}}}$.

Note that the results of Section 3 do not explain why the weights became positive in this experiment, because Matérn kernels do not appear to be totally positive even if the differentiability requirements were to be relaxed and only single zeroes counted (recall Remark 1). We have numerically observed that selecting $n>\rho+1/2$ and point sets such that $\max X<\min Y$ makes the matrix $\boldsymbol{K}_{Y,X}$ discussed in Section 3.3 singular for the Matérn kernel $k_{\rho}$. This implies that there is a non-trivial linear combination of the $n$ Matérn translates at $X$ that vanishes at more than $n-1$ points. When $\rho=1/2$ and $\ell=1$ (so that $k_{\rho}(x,x^{\prime})=\mathrm{e}^{-\mathinner{\lvert x-y\rvert}}$), an analytical counterexample can be constructed by setting $X=\{x_{1},x_{2}\}$ and $Y=\{x_{1}+h,x_{2}+h\}$ with $h>x_{2}-x_{1}$. Then

[TABLE]

which is not invertible because multiplying the first row by $\mathrm{e}^{-(x_{1}-x_{2})}$ yields the second row. Therefore the positivity of quadrature weights for Matérns and other kernels with finite smoothness requires a further research.

Acknowledgements.

TK was supported by the Aalto ELEC Doctoral School. MK acknowledges support by the European Research Council (StG Project PANAMA). SS was supported by the Academy of Finland project 313708. This material was developed, in part, at the Prob Num 2018 workshop hosted by the Lloyd’s Register Foundation programme on Data-Centric Engineering at the Alan Turing Institute, UK, and supported by the National Science Foundation, USA, under Grant DMS-1127914 to the Statistical and Applied Mathematical Sciences Institute. Any opinions, findings, conclusions or recommendations expressed in this material are those of the author(s) and do not necessarily reflect the views of the above-named funding bodies and research institutions.