Gaussian kernel quadrature at scaled Gauss-Hermite nodes

TL;DR

This paper presents a precise, stable method for computing kernel quadrature weights using scaled Gauss-Hermite nodes for Gaussian kernels, with proven exponential convergence for functions in the associated RKHS.

Contribution

It introduces an explicit, numerically stable approximation to kernel quadrature weights on Gaussian measure using scaled Gauss-Hermite nodes and Mercer eigendecomposition truncation.

Findings

Kernel quadrature weights are positive.

The approximation achieves exponential convergence.

Numerical results confirm stability and accuracy.

Abstract

This article derives an accurate, explicit, and numerically stable approximation to the kernel quadrature weights in one dimension and on tensor product grids when the kernel and integration measure are Gaussian. The approximation is based on use of scaled Gauss-Hermite nodes and truncation of the Mercer eigendecomposition of the Gaussian kernel. Numerical evidence indicates that both the kernel quadrature and the approximate weights at these nodes are positive. An exponential rate of convergence for functions in the reproducing kernel Hilbert space induced by the Gaussian kernel is proved under an assumption on growth of the sum of absolute values of the approximate weights.

Equations

$$k(x,y)=\exp\bigg{(}\!-\frac{(x-y)^{2}}{2\ell^{2}}\bigg{)}$$

$$Q_k(f):=\sum _{n=1}^N{w}_{k,n}f(x_n)\approx \mu (f):=\frac{1}{\sqrt{2\pi }}{\int }_{\mathbb{R}}f(x)e^{-x^2/2}dx,$$

$$K{w}_k=k_{\mu },$$

$${\tilde{x}}_n:=\frac{1}{\sqrt{2}\alpha \beta }{x}_n^{\text{GH}}$$

$$\widetilde{w}_{k,n}\coloneqq\bigg{(}\frac{1}{1+2\delta^{2}}\bigg{)}^{1/2}w_{n}^{\text{\tiny{GH}}}\operatorname{e}^{\delta^{2}\tilde{x}_{n}^{2}}\sum_{m=0}^{\lfloor(N-1)/2\rfloor}\frac{1}{2^{m}m!}\bigg{(}\!\frac{2\alpha^{2}\beta^{2}}{1+2\delta^{2}}-1\bigg{)}^{m}\mathrm{H}_{2m}(x_{n}^{\text{\tiny{GH}}}),$$

$$Q_k^d=Q_{k,1}\otimes \cdots \otimes Q_{k,d},$$

$$k(x,y)=\sum _{n=0}^{\infty }{\lambda }_n{\phi }_n(x){\phi }_n(y)$$

$$d{\mu }_{\alpha }(x):=\frac{\alpha }{\sqrt{\pi }}{e}^{-{\alpha }^2{x}^2}dx,$$

$${\mathrm{H}}_n(x):=(-1{)}^n{e}^{x^2/2}\frac{d{}^n}{d{x}^n}{e}^{-x^2/2}$$

$$\varepsilon=\frac{1}{\sqrt{2}\ell},\hskip 14.22636pt\beta=\bigg{(}1+\bigg{(}\frac{2\varepsilon}{\alpha}\bigg{)}^{2}\bigg{)}^{1/4},\hskip 14.22636pt\text{and}\hskip 14.22636pt\delta^{2}=\frac{\alpha^{2}}{2}(\beta^{2}-1)$$

$$\lambda_{n}^{\alpha}\coloneqq\sqrt{\frac{\alpha^{2}}{\alpha^{2}+\delta^{2}+\varepsilon^{2}}}\bigg{(}\frac{\varepsilon^{2}}{\alpha^{2}+\delta^{2}+\varepsilon^{2}}\bigg{)}^{n}$$

$$\varphi_{n}^{\alpha}(x)\coloneqq\sqrt{\frac{\beta}{n!}}\operatorname{e}^{-\delta^{2}x^{2}}\mathrm{H}_{n}\big{(}\sqrt{2}\alpha\beta x\big{)}.$$

$$\mu(\varphi_{2m+1}^{\alpha})=0\hskip 14.22636pt\text{and}\hskip 14.22636pt\mu(\varphi_{2m}^{\alpha})=\bigg{(}\frac{\beta}{1+2\delta^{2}}\bigg{)}^{1/2}\frac{\sqrt{(2m)!}}{2^{m}m!}\bigg{(}\frac{2\alpha^{2}\beta^{2}}{1+2\delta^{2}}-1\bigg{)}^{m}$$

$$\mathrm{H}_{2m}(x)=\frac{(2m)!}{2^{m}}\sum_{p=0}^{m}\frac{(-1)^{m-p}}{(2p)!(m-p)!}\big{(}\sqrt{2}x\big{)}^{2p},$$

$${\int }_{\mathbb{R}}{x}^{2p}{e}^{-{\delta }^2{x}^2}d\mu (x)=\frac{1}{\sqrt{2\pi }}{\int }_{\mathbb{R}}{x}^{2p}{e}^{-({\delta }^2+1/2)x^2}dx=\frac{(2p)!}{2^{p}p!(1+2{\delta }^2{)}^{p+1/2}},$$

$$\begin{split}\mu(\varphi_{2m}^{\alpha})&=\frac{\sqrt{(2m)!\beta}}{2^{m}}\sum_{p=0}^{m}\frac{(-1)^{m-p}}{(2p)!(m-p)!}(2\alpha\beta)^{2p}\int_{\mathbb{R}}x^{2p}\operatorname{e}^{-\delta^{2}x^{2}}\operatorname{d\!}\mu(x)\\ &=\frac{(-1)^{m}\sqrt{(2m)!\beta}}{2^{m}\sqrt{1+2\delta^{2}}}\sum_{p=0}^{m}\frac{1}{p!(m-p)!}\bigg{(}-\frac{2\alpha^{2}\beta^{2}}{1+2\delta^{2}}\bigg{)}^{p}\\ &=\frac{(-1)^{m}\sqrt{(2m)!\beta}}{2^{m}m!\sqrt{1+2\delta^{2}}}\sum_{p=0}^{m}{m\choose p}\bigg{(}-\frac{2\alpha^{2}\beta^{2}}{1+2\delta^{2}}\bigg{)}^{p}\\ &=\bigg{(}\frac{\beta}{1+2\delta^{2}}\bigg{)}^{1/2}\frac{\sqrt{(2m)!}}{2^{m}m!}\bigg{(}\frac{2\alpha^{2}\beta^{2}}{1+2\delta^{2}}-1\bigg{)}^{m}.\end{split}$$

$$\widetilde{w}_{k}^{M}\coloneqq\big{(}\Phi\Lambda\Phi^{\mathsf{T}}\big{)}^{-1}\Phi\Lambda\varphi_{\mu}.$$

$$\Phi =QR:=Q\left[\begin{array}{cc}{R}_1& R_2\end{array}\right]$$

$$\widetilde{w}_{k}^{M}=\big{(}QR\Lambda R^{\mathsf{T}}Q^{\mathsf{T}}\big{)}^{-1}QR\Lambda\varphi_{\mu}=Q\big{(}R\Lambda R^{\mathsf{T}}\big{)}^{-1}R\Lambda\varphi_{\mu}.$$

$$\Lambda =\left[\begin{array}{cc}{\Lambda }_1& 0\\ 0& {\Lambda }_2\end{array}\right]$$

$$R\Lambda R^{\mathsf{T}}=R_{1}\Lambda_{1}\big{(}R_{1}^{\mathsf{T}}+\Lambda_{1}^{-1}R_{1}^{-1}R_{2}\Lambda_{2}R_{2}^{\mathsf{T}}\big{)}.$$

$$\widetilde{w}_{k}^{M}=Q\big{(}R_{1}^{\mathsf{T}}+\Lambda_{1}^{-1}R_{1}^{-1}R_{2}\Lambda_{2}R_{2}^{\mathsf{T}}\big{)}^{-1}\begin{bmatrix}I_{N}&\Lambda_{1}^{-1}R_{1}^{-1}R_{2}\Lambda_{2}\end{bmatrix}\varphi_{\mu},$$

$$\sum _{n=1}^N{w}_n^{\text{GH}}p(x_n^{\text{GH}})=\mu (p)$$

$$w_{k}\approx\widetilde{w}_{k}\coloneqq\widetilde{w}_{k}^{N}=\big{(}\Phi\Lambda\Phi^{\mathsf{T}}\big{)}^{-1}\Phi\Lambda\varphi_{\mu}=\Phi^{-\mathsf{T}}\varphi_{\mu}.$$

$$[V]_{ij}\coloneqq\frac{1}{\sqrt{(j-1)!}}\mathrm{H}_{j-1}\big{(}\sqrt{2}\alpha\beta x_{i}\big{)}$$

$${\tilde{w}}_k=\frac{1}{\sqrt{\beta }}E{V}^{-\mathsf{T}}{\phi }_{\mu }.$$

$${\tilde{x}}_n:=\frac{1}{\sqrt{2}\alpha \beta }{x}_n^{\text{GH}}.$$

$${\tilde{w}}_k=\frac{1}{\sqrt{\beta }}E{V}^{-\mathsf{T}}{\phi }_{\mu }=\frac{1}{\sqrt{\beta }}E{W}_{\text{GH}}V{\phi }_{\mu }.$$

$${\tilde{x}}_n=\frac{1}{\sqrt{2}\alpha \beta }{x}_n^{\text{GH}}.$$

$${\tilde{Q}}_k(f):=\sum _{n=1}^N{\tilde{w}}_{k,n}f({\tilde{x}}_n),$$

Full text

††journal: BIT

1 Introduction

Let $\mu$ be the standard Gaussian measure on $\mathbb{R}$ and $f\colon\mathbb{R}\to\mathbb{R}$ a measurable function. We consider the problem of numerical computation of the integral with respect to $\mu$ of $f$ using a kernel quadrature rule (we reserve the term cubature for rules on higher dimensions) based on the Gaussian kernel

[TABLE]

with the length-scale $\ell>0$. Given any distinct nodes $x_{1},\ldots,x_{N}$, the kernel quadrature rule is an approximation of the form

[TABLE]

with its weights $w_{k}=(w_{k,1},\ldots,w_{k,N})\in\mathbb{R}^{N}$ solved from the linear system of equations

[TABLE]

where $[K]_{ij}\coloneqq k(x_{i},x_{j})$ and $[k_{\mu}]_{i}\coloneqq\int_{\mathbb{R}}k(x_{i},x)\operatorname{d\!}\mu(x)$. This is equivalent to uniquely selecting the weights such that the $N$ kernel translates $k(x_{1},\cdot),\ldots,k(x_{N},\cdot)$ are integrated exactly by the quadrature rule. Kernel quadrature rules can be interpreted as best quadrature rules in the reproducing kernel Hilbert space (RKHS) induced by a positive-definite kernel (Larkin, 1970), integrated kernel (radial basis function) interpolants (Bezhaev, 1991; Sommariva and Vianello, 2006), and posteriors to $\mu(f)$ under a Gaussian process prior on the integrand (Larkin, 1972; O’Hagan, 1991; Briol et al., 2019).

Recently, Fasshauer and McCourt (2012) have developed a method to circumvent the well-known problem that interpolation with the Gaussian kernel becomes often numerically unstable—in particular when $\ell$ is large—because the condition number of $K$ tends to grow with an exponential rate (Schaback, 1995). They do this by truncating the Mercer eigendecomposition of the Gaussian kernel after $M$ terms and replacing the interpolation basis $\{k(x_{n},\cdot)\}_{n=1}^{N}$ with the first $M$ eigenfunctions. In this article we show that application of this method with $M=N$ to kernel quadrature yields, when the nodes are selected by a suitable and fairly natural scaling of the nodes of the classical Gauss–Hermite quadrature rule, an accurate, explicit, and numerically stable approximation to the Gaussian kernel quadrature weights. Moreover, the proposed nodes appear to be a good and natural choice for the Gaussian kernel quadrature.

To be precise, Theorem 2.2 states that the quadrature rule $\widetilde{Q}_{k}$ that exactly integrates the first $N$ Mercer eigenfunctions of the Gaussian kernel and uses the nodes

[TABLE]

has the weights

[TABLE]

$\widetilde{w}_{k}=(\widetilde{w}_{k,1},\ldots,\widetilde{w}_{k,N})\in\mathbb{R}^{N}$, where $\alpha$ (for which the value $1/\sqrt{2}$ seems the most natural), $\beta$, and $\delta$ are constants defined in Equation 5, $\mathrm{H}_{n}$ are the probabilists’ Hermite polynomials (4), and $x_{n}^{\text{\tiny{GH}}}$ and $w_{n}^{\text{\tiny{GH}}}$ are the nodes and weights of the $N$-point Gauss–Hermite quadrature rule. We argue that these weights are a good approximation to $w_{k}$ and accordingly call them approximate Gaussian kernel quadrature weights. Although we derive no bounds for the error of this weight approximation, numerical experiments in Section 5 indicate that the approximation is accurate and that it appears that $\widetilde{w}_{k}\to w_{k}$ as $N\to\infty$. In Section 4 we extend the weight approximation for $d$-dimensional Gaussian tensor product kernel cubature rules of the form

[TABLE]

where $Q_{k,i}$ are one-dimensional Gaussian kernel quadrature rules. Since each weight of $Q_{k}^{d}$ is a product of weights of the univariate rules, an approximation for the tensor product weights is readily available.

It turns out that the approximate weight and the associated nodes $\tilde{x}_{n}$ have a number of desirable properties:

• •

We are not aware of any work on efficient selection of “good” nodes in the setting of this article. The Gauss–Hermite nodes (O’Hagan, 1991, Section 3) and random points (Rasmussen and Ghahramani, 2002) are often used, but one should clearly be able to do better, while computation of the optimal nodes (Oettershagen, 2017, Section 5.2) is computationally demanding. As such, given the desirable properties, listed below, of the resulting kernel quadrature rules, the nodes $\tilde{x}_{n}$ appear to be an excellent heuristic choice. These nodes also behave naturally when $\ell\to\infty$; see Section 2.5.

• •

Numerical experiments in Section 5.3 suggest that both $w_{k,n}$ (for the nodes $\tilde{x}_{n}$) and $\widetilde{w}_{k,n}$ are positive for any $N\in\mathbb{N}$ and every $n=1,\ldots,N$. Besides the optimal nodes, the weights for which are guaranteed to be positive when the Gaussian kernel is used (Richter-Dyn, 1971; Oettershagen, 2017), there are no node configurations that give rise to positive weights as far as we are aware of.

• •

Numerical experiments in Sections 5.1 and 5.3 demonstrate that computation of the approximate weights is numerically stable. Furthermore, construction of these weights only incurs a quadratic computational cost in the number of points, as opposed to the cubic cost of solving $w_{k}$ from Equation 2. See Section 2.6 for more details. Note that to obtain a numerically stable method it is not necessary to use the nodes $\tilde{x}_{n}$ as the method in (Fasshauer and McCourt, 2012) can be applied in a straightforward manner for any nodes. However, doing so one forgoes a closed form expression and has to use the QR decomposition.

• •

In Sections 3 and 4 we show that slow enough growth with $N$ of $\sum_{i=1}^{N}\mathinner{\lvert\widetilde{w}_{k,n}\rvert}$ (numerical evidence indicates this sum converges to one) guarantees that the approximate Gaussian kernel quadrature rule—as well as the corresponding tensor product version—converges with an exponential rate for functions in the RKHS of the Gaussian kernel. Convergence analysis is based on analysis of magnitude of the remainder of the Mercer expansion and rather explicit bounds on Hermite polynomials and their roots. Magnitude of the nodes $\tilde{x}_{n}$ is crucial for the analysis; if they were further spread out the proofs would not work as such.

• •

We find the connection to the Gauss–Hermite weights and nodes that the closed form expression for $\widetilde{w}_{k}$ provides intriguing and hope that it can be at some point used to furnish, for example, a rigorous proof of positivity of the approximate weights.

2 Approximate weights

This section contains the main results of this article. The main contribution is derivation, in Theorem 2.2, of the weights $\widetilde{w}_{k}$, that can be used to approximate the kernel quadrature weights. We also discuss positivity of these weights, the effect the kernel length-scale $\ell$ is expected to have on quality of the approximation, and computational complexity.

2.1 Eigendecomposition of the Gaussian kernel

Let $\nu$ be a probability measure on the real line. If the support of $\nu$ is compact, Mercer’s theorem guarantees that any positive-definite kernel $k$ admits an absolutely and uniformly convergent eigendecomposition

[TABLE]

for positive and non-increasing eigenvalues $\lambda_{n}$ and eigenfunctions $\varphi_{n}$ that are included in the RKHS $\mathcal{H}$ induced by $k$ and orthonormal in $L^{2}(\nu)$. Moreover, $\sqrt{\lambda_{n}}\varphi_{n}$ are $\mathcal{H}$-orthonormal. If the support of $\nu$ is not compact, the expansion (3) converges absolutely and uniformly on all compact subsets of $\mathbb{R}\times\mathbb{R}$ under some mild assumptions (Sun, 2005; Steinwart and Scovel, 2012). For the Gaussian kernel (1) and measure the eigenvalues and eigenfunctions are available analytically. For a collection of explicit eigendecompositions of some other kernels, see for instance (Fasshauer and McCourt, 2015, Appendix A)

Let $\mu_{\alpha}$ stand for the Gaussian probability measure,

[TABLE]

with variance $1/(2\alpha^{2})$ (i.e., $\mu=\mu_{1/\sqrt{2}}\,$) and

[TABLE]

for the (unnormalised) probabilists’ Hermite polynomial satisfying the orthogonality property $\langle\mathrm{H}_{n},\mathrm{H}_{m}\rangle_{L^{2}(\mu)}=n!\,\delta_{nm}$. Denote

[TABLE]

and note that $\beta>1$ and $\delta^{2}>0$. Then the eigenvalues and $L^{2}(\mu_{\alpha})$-orthonormal eigenfunctions of the Gaussian kernel are (Fasshauer and McCourt, 2012)

[TABLE]

and

[TABLE]

See (Fasshauer and McCourt, 2015, Section 12.2.1) for verification that these indeed are Mercer eigenfunctions and eigenvalues for the Gaussian kernel. The role of the parameter $\alpha$ is discussed in Section 2.4. The following result, also derivable from Equation 22.13.17 in (Abramowitz and Stegun, 1964), will be useful.

Lemma 1

The eigenfunctions (7) of the Gaussian kernel (1) satisfy

[TABLE]

for $m\geq 0$.

Proof

Since an Hermite polynomial of odd order is an odd function, $\mu(\varphi_{2m+1}^{\alpha})=0$. For even indices, use the explicit expression

[TABLE]

the Gaussian moment formula

[TABLE]

and the binomial theorem to conclude that

[TABLE]

∎

2.2 Approximation via QR decomposition

We begin by outlining a straightforward extension to kernel quadrature of the work of Fasshauer and McCourt (2012, 2015, Chapter 13) on numerically stable kernel interpolation. Recall that the kernel quadrature weights $w_{k}\in\mathbb{R}^{N}$ at distinct nodes $x_{1},\ldots,x_{N}$ are solved from the linear system $Kw_{k}=k_{\mu}$ with ${[K]_{ij}=k(x_{i},x_{j})}$ and ${[k_{\mu}]_{i}=\int_{\mathbb{R}}k(x_{i},x)\operatorname{d\!}\mu(x)}$. Truncation of the eigendecomposition (3) after $M\geq N$ terms111Low-rank approximations (i.e. $M<N$) are also possible (Fasshauer and McCourt, 2012, Section 6.1). yields the approximations $K\approx\Phi\Lambda\Phi^{\mathsf{T}}$ and $k_{\mu}\approx\Phi\Lambda\varphi_{\mu}$, where ${[\Phi]_{ij}\coloneqq\varphi_{j-1}^{\alpha}(x_{i})}$ is an $N\times M$ matrix, the diagonal $M\times M$ matrix $[\Lambda]_{ii}\coloneqq\lambda_{i-1}$ contains the eigenvalues in appropriate order, and ${[\varphi_{\mu}]_{i}\coloneqq\mu(\varphi_{i-1})}$ is an $M$-vector. The kernel quadrature weights $w_{k}$ can be therefore approximated by

[TABLE]

Equation 8 can be written in a more convenient form by exploiting the QR decomposition. The QR decomposition of $\Phi$ is

[TABLE]

for a unitary $Q\in\mathbb{R}^{N\times N}$, an upper triangular $R_{1}\in\mathbb{R}^{N\times N}$, and $R_{2}\in\mathbb{R}^{N\times(M-N)}$. Consequently,

[TABLE]

The decomposition

[TABLE]

of $\Lambda\in\mathbb{R}^{M\times M}$ into diagonal $\Lambda_{1}\in\mathbb{R}^{N\times N}$ and $\Lambda_{2}\in\mathbb{R}^{(M-N)\times(M-N)}$ allows for writing

[TABLE]

Therefore,

[TABLE]

where $I_{N}$ is the $N\times N$ identity matrix. If $\varepsilon^{2}/(\alpha^{2}+\delta^{2}+\varepsilon^{2})$ is small (i.e., $\ell$ is large), numerical ill-conditioning in Equation 9 for the Gaussian kernel is associated with the diagonal matrices $\Lambda_{1}^{-1}$ and $\Lambda_{2}$. Consequently, numerical stability can be significantly improved by performing the multiplications by these matrices in the terms $\Lambda_{1}^{-1}R_{1}^{-1}R_{2}\Lambda_{2}R_{2}^{\mathsf{T}}$ and $\Lambda_{1}^{-1}R_{1}^{-1}R_{2}\Lambda_{2}$ analytically; see (Fasshauer and McCourt, 2012, Sections 4.1 and 4.2) for more details.

Unfortunately, using the QR decomposition does not provide an attractive closed form solution for the approximate weights $\widetilde{w}_{k}^{M}$ for general $M$. Setting $M=N$ turns $\Phi$ into a square matrix, enabling its direct inversion and formation of an explicit connection to the classical Gauss–Hermite quadrature. The rest of the article is concerned with this special case.

2.3 Gauss–Hermite quadrature

Given a measure $\nu$ on $\mathbb{R}$, the $N$-point Gaussian quadrature rule is the unique $N$-point quadrature rule that is exact for all polynomials of degree at most $2N-1$. We are interested in Gauss–Hermite quadrature rules that are Gaussian rules for the Gaussian measure $\mu$:

[TABLE]

for every polynomial $p\colon\mathbb{R}\to\mathbb{R}$ with $\deg p\leq 2N-1$. The nodes $x_{1}^{\text{\tiny{GH}}},\ldots x_{N}^{\text{\tiny{GH}}}$ are the roots of the $N$th Hermite polynomial $\mathrm{H}_{N}$ and the weights $w_{1}^{\text{\tiny{GH}}},\ldots,w_{N}^{\text{\tiny{GH}}}$ are positive and sum to one. The nodes and the weights are related to the eigenvalues and eigenvectors of the tridiagonal Jacobi matrix formed out of three-term recurrence relation coefficients of normalised Hermite polynomials (Gautschi, 2004, Theorem 3.1).

We make use of the following theorem, a one-dimensional special case of a more general result due to Mysovskikh (1968). See also (Cools, 1997, Section 7). This result also follows from the Christoffel–Darboux formula (14).

Theorem 2.1

Let $\nu$ be a measure on $\mathbb{R}$. Suppose that $x_{1},\ldots,x_{N}$ and $w_{1},\ldots,w_{N}$ are the nodes and weights of the unique Gaussian quadrature rule. Let $p_{0},\ldots,p_{N-1}$ be the $L^{2}(\nu)$-orthonormal polynomials. Then the matrix $[P]_{ij}\coloneqq\sum_{n=0}^{N-1}p_{n}(x_{i})p_{n}(x_{j})$ is diagonal and has the diagonal elements $[P]_{ii}=1/w_{i}$.

2.4 Approximate weights at scaled Gauss–Hermite nodes

Let us now consider the approximate weights (8) with $M=N$. Assuming that $\Phi$ is invertible, we then have

[TABLE]

Note that the exponentially decaying Mercer eigenvalues, a major source of numerical instability, do not appear in the equation for $\widetilde{w}_{k}$. The weights $\widetilde{w}_{k}$ are those of the unique quadrature rule that is exact for the $N$ first eigenfunctions $\varphi_{0}^{\alpha},\ldots,\varphi_{N-1}^{\alpha}$. For the Gaussian kernel, we are in a position to do much more. Recalling the form of the eigenfunctions in Equation 7, we can write $\Phi=\sqrt{\beta}E^{-1}V$ for the diagonal matrix $[E]_{ii}\coloneqq\operatorname{e}^{\delta^{2}x_{i}^{2}}$ and the Vandermonde matrix

[TABLE]

of scaled and normalised Hermite polynomials. From this it is evident that $\Phi$ is invertible—which is just a manifestation of the fact that the eigenfunctions of a totally positive kernel constitute a Chebyshev system (Kellog, 1918; Pinkus, 1996). Consequently,

[TABLE]

Select the nodes

[TABLE]

Then the matrix $V$ defined in Equation 10 is precisely the Vandermonde matrix of the normalised Hermite polynomials and $VV^{\mathsf{T}}$ is the matrix $P$ of Theorem 2.1. Let $W_{\text{\tiny{GH}}}$ be the diagonal matrix containing the Gauss–Hermite weights. It follows that $V^{-\mathsf{T}}=W_{\text{\tiny{GH}}}V$ and

[TABLE]

Combining this equation with Lemma 1, we obtain the main result of this article.

Theorem 2.2

Let $x_{1}^{\text{\tiny{GH}}},\ldots,x_{N}^{\text{\tiny{GH}}}$ and $w_{1}^{\text{\tiny{GH}}},\ldots,w_{N}^{\text{\tiny{GH}}}$ stand for the nodes and weights of the $N$-point Gauss–Hermite quadrature rule. Define the nodes

[TABLE]

Then the weights $\widetilde{w}_{k}\in\mathbb{R}^{N}$ of the $N$-point quadrature rule

[TABLE]

defined by the exactness conditions $\widetilde{Q}_{k}(\varphi_{n}^{\alpha})=\mu_{\alpha}(\varphi_{n}^{\alpha})$ for $n=0,\ldots,N-1$, are

[TABLE]

where $\alpha$, $\beta$, and $\delta$ are defined in Equation 5 and $\mathrm{H}_{2m}$ are the probabilists’ Hermite polynomials (4).

Since the weights $\widetilde{w}_{k}$ are obtained by truncating of the Mercer expansion of $k$, it is to be expected that $\widetilde{w}_{k}\approx w_{k}$. This motivates our calling of these weights the approximate Gaussian kernel quadrature weights. We do not provide theoretical results on quality of this approximation, but the numerical experiments in Section 5.2 indicate that the approximation is accurate and that its accuracy increases with $N$. See (Fasshauer and McCourt, 2012) for related experiments.

An alternative non-analytical formula for the approximate weights can be derived using the Christoffel–Darboux formula (Gautschi, 2004, Section 1.3.3)

[TABLE]

From Equation 11 we then obtain (keep in mind that $x_{1}^{\text{\tiny{GH}}},\ldots,x_{N}^{\text{\tiny{GH}}}$ are the roots of $\mathrm{H}_{N}$)

[TABLE]

This formula is analogous to the formula

[TABLE]

for the Gauss–Hermite weights. Plugging this in, we get

[TABLE]

It appears that both $w_{k,n}$ and $\widetilde{w}_{k,n}$ of Theorem 2.2 are positive for many choices of $\alpha$; see Section 5.3 for experiments involving $\alpha=1/\sqrt{2}$. Unfortunately, we have not been able to prove this. In fact, numerical evidence indicates something slightly stronger. Namely that the even polynomial

[TABLE]

of degree $2\lfloor(N-1)/2\rfloor$ is positive for every $N\geq 1$ and (at least) every $0<\gamma\leq 1$. This would imply positivity of $\widetilde{w}_{k,n}$ since the Gauss–Hermite weights $w_{n}^{\text{\tiny{GH}}}$ are positive. For example, with $\alpha=1/\sqrt{2}$,

[TABLE]

As discussed in (Fasshauer and McCourt, 2012) in the context of kernel interpolation, the parameter $\alpha$ acts as a global scale parameter. While in interpolation it is not entirely clear how this parameter should be selected, in quadrature it seems natural to set $\alpha=1/\sqrt{2}$ so that the eigenfunctions are orthonormal in $L^{2}(\mu)$. This is the value that we use, though also other values are potentially of interest since $\alpha$ can be used to control the spread of the nodes independently of the length-scale $\ell$. In Section 3, we also see that this value leads to more natural convergence analysis.

2.5 Effect of the length-scale

Roughly speaking, magnitude of the eigenvalues

[TABLE]

determines how many eigenfunctions are necessary for an accurate weight approximation. We therefore expect that the approximation (13) is less accurate when the length-scale $\ell$ is small (i.e., $\varepsilon=1/(\sqrt{2}\ell)$ is large). This is confirmed by the numerical experiments in Section 5.

Consider then the case $\ell\to\infty$. This scenario is called the flat limit in scattered data approximation literature where it has been proved222It is interesting to note that the first published observation of analogous phenomenon is, as far as we are aware of, due to O’Hagan (1991, Section 3.3) in kernel quadrature literature, predating the work of Driscoll and Fornberg (2002). See also (Minka, 2000) for early quadrature-related work on the topic. that the kernel interpolant associated to an isotropic kernel with increasing length-scale converges to (i) the unique polynomial interpolant of degree $N-1$ to the data if the kernel is infinitely smooth (Larsson and Fornberg, 2005; Schaback, 2005; Lee et al., 2007) or (ii) to a polyharmonic spline interpolant if the kernel is of finite smoothness (Lee et al., 2014). In our case, $\ell\to\infty$ results in

[TABLE]

If the nodes are selected as in Equation 12, $\tilde{x}_{n}\to x_{n}^{\text{\tiny{GH}}}/(\sqrt{2}\alpha)$. That is, if $\alpha=1/\sqrt{2}$

[TABLE]

That the approximate weights convergence to the Gauss–Hermite ones can be seen, for example, from Equation 13 by noting that only the first term in the sum is retained at the limit. Based on the aforementioned results regarding convergence of kernel interpolants to polynomial ones at the flat limit, it is to be expected that also $w_{k,n}\to w_{n}^{\text{\tiny{GH}}}$ as $\ell\to\infty$ (we do not attempt to prove this). Because the Gauss–Hermite quadrature rule is the “best” for polynomials and kernel interpolants convergence to polynomials at the flat limit, the above observation provides another justification for the choice $\alpha=1/\sqrt{2}$ that we proposed the preceding section.

When it comes to node placement, the length-scale is having an intuitive effect if the nodes are selected according to Equation 12. For small $\ell$, the nodes are placed closer to the origin where most of the measure is concentrated as integrands are expected to converge quickly to zero as $\mathinner{\lvert x\rvert}\to\infty$, whereas for larger $\ell$ the nodes are more—but not unlimitedly—spread out in order to capture behaviour of functions that potentially contribute to the integral also further away from the origin.

2.6 On computational complexity

Because the Gauss–Hermite nodes and weights are related to the eigenvalues and eigenvectors of the tridiagonal Jacobi matrix (Gautschi, 2004, Theorem 3.1) they—and the points $\tilde{x}_{n}$—can be solved in quadratic time (in practice, these nodes and weights can be often tabulated beforehand). From Equation 13 it is seen that computation of each approximate weight is linear in $N$: there are approximately $(N-1)/2$ terms in the sum and the Hermite polynomials can be evaluated on the fly using the three-term recurrence formula $\mathrm{H}_{n+1}(x)=x\mathrm{H}_{n}(x)-n\mathrm{H}_{n-1}(x)$. That is, computational cost of obtaining $\tilde{x}_{n}$ and $\widetilde{w}_{k,n}$ for $n=1,\ldots,N$ is quadratic in $N$. Since the kernel matrix $K$ of the Gaussian kernel is dense, solving the exact kernel quadrature weights from the linear system (2) for the points $\tilde{x}_{n}$ incurs a more demanding cubic computational cost. Because computational cost of a tensor product rule does not depend on the nodes and weights after these have been computed, the above discussion also applies to the rules presented in Section 4.

3 Convergence analysis

In this section we analyse convergence in the reproducing kernel Hilbert space $\mathcal{H}\subset C^{\infty}(\mathbb{R})$ induced by the Gaussian kernel of quadrature rules that are exact for the Mercer eigenfunctions. First, we prove a generic result (Theorem 3.1) to this effect and then apply this to the quadrature rule with the nodes $\tilde{x}_{n}$ and weights $\widetilde{w}_{k,n}$. If $\sum_{n=1}^{N}\mathinner{\lvert\widetilde{w}_{k,n}\rvert}$ does not grow too fast with $N$, we obtain exponential convergence rates.

Recall some basic facts about reproducing kernel Hilbert spaces spaces (Berlinet and Thomas-Agnan, 2004): (i) $\langle f,k(x,\cdot)\rangle_{\mathcal{H}}=f(x)$ for any $f\in\mathcal{H}$ and $x\in\mathbb{R}$ and (ii) $f=\sum_{n=0}^{\infty}\lambda_{n}^{\alpha}\langle f,\varphi_{n}^{\alpha}\rangle\varphi_{n}^{\alpha}$ for any $f\in\mathcal{H}$. The worst-case error $e(Q)$ of a quadrature rule $Q(f)=\sum_{n=1}^{N}w_{n}f(x_{n})$ is

[TABLE]

Crucially, the worst-case error satisfies

[TABLE]

for any $f\in\mathcal{H}$. This justifies calling a sequence $\{Q_{N}\}_{N=1}^{\infty}$ of $N$-point quadrature rules convergent if $e(Q_{N})\to 0$ as $N\to\infty$. For given nodes $x_{1},\ldots,x_{N}$, the weights $w_{k}=(w_{k,1},\ldots,w_{k,N})$ of the kernel quadrature rule $Q_{k}$ are unique minimisers of the worst-case error:

[TABLE]

It follows that a rate of convergence to zero for $e(Q)$ also applies to $e(Q_{k})$.

A number of convergence results for kernel quadrature rules on compact spaces appear in (Bezhaev, 1991; Kanagawa et al., 2019; Briol et al., 2019). When it comes to the RKHS of the Gaussian kernel, characterised in (Steinwart et al., 2006; Minh, 2010), Kuo and Woźniakowski (2012) have analysed convergence of the Gauss–Hermite quadrature rule. Unfortunately, it turns out that the Gauss–Hermite rule converges in this space if and only if $\varepsilon^{2}<1/2$. Consequently, we believe that the analysis below is the first to establish convergence, under the assumption (supported by our numerical experiments) that the sum of $\mathinner{\lvert\widetilde{w}_{k,n}\rvert}$ does not grow too fast, of an explicitly constructed sequence of quadrature rules in the RKHS of the Gaussian kernel with any value of the length-scale parameter. We begin with two simple lemmas.

Lemma 2

The eigenfunctions $\varphi_{n}^{\alpha}$ admit the bound

[TABLE]

for a constant $K\leq 1.087$ and every $x\in\mathbb{R}$.

Proof

For each $n\geq 0$, the Hermite polynomials obey the bound

[TABLE]

for a constant $K\leq 1.087$ (Erdélyi, 1953, p. 208). See (Bonan and Clark, 1990) for other such bounds333In particular, the factor $n^{-1/6}$ could be added on the right-hand side. This would make little difference in convergence analysis of Theorem 3.1.. Thus

[TABLE]

∎

Lemma 3

Let $\alpha=1/\sqrt{2}$. Then

[TABLE]

for every $\ell>0$ if and only if $\rho\leq 2$.

Proof

The function

[TABLE]

satisfies $\gamma(0)=0$ and $\gamma(\varepsilon^{2})\to 1$ as $\varepsilon^{2}\to\infty$. The derivative

[TABLE]

is positive when $\rho\leq 2$. For $\rho>2$, the derivative has a single root at $\varepsilon_{0}^{2}=1/(4(\rho-2))$ so that $\gamma(\varepsilon_{0}^{2})>1$. That is, $\gamma(\varepsilon^{2})\in(0,1)$, and consequently $\gamma(\varepsilon^{2})^{1/2}\in(0,1)$, if and only if $\rho\leq 2$. ∎

Theorem 3.1

Let $\alpha=1/\sqrt{2}$. Suppose that the nodes $x_{1},\ldots,x_{N}$ and weights $w_{1},\ldots,w_{N}$ of an $N$-point quadrature rule $Q_{N}$ satisfy

$\sum_{n=1}^{N}\mathinner{\lvert w_{n}\rvert}\leq W_{N}$* for some $W_{N}\geq 0$;* 2. 2.

$Q_{N}(\varphi_{n}^{\alpha})=\mu(\varphi_{n}^{\alpha})$* for each $n=0,\ldots,M_{N}-1$ for some $M_{N}\geq 1$;* 3. 3.

$\sup_{1\leq n\leq N}\mathinner{\lvert x_{n}\rvert}\leq 2\sqrt{M_{N}}/\beta$.

Then there exist constants $C_{1},C_{2}>0$, independent of $N$ and $Q_{N}$, and $0<\eta<1$ such that

[TABLE]

Explicit forms of these constants appear in Equation 18.

Proof

For notational convenience, denote

[TABLE]

and $\varphi_{n}=\varphi_{n}^{\alpha}$. Because every $f\in\mathcal{H}$ admits the expansion $f=\sum_{n=0}^{\infty}\lambda_{n}\langle f,\varphi_{n}\rangle_{\mathcal{H}}\varphi_{n}$ and $Q_{N}(\varphi_{n})=\mu(\varphi_{n})$ for $n<M_{N}$, it follows from the Cauchy–Schwarz inequality and $\mathinner{\lVert\varphi_{n}\rVert}_{\mathcal{H}}=1/\sqrt{\lambda_{n}}$ that

[TABLE]

From Lemma 2 we have $\mathinner{\lvert\varphi_{n}(x)\rvert}\leq K\sqrt{\beta}\operatorname{e}^{x^{2}/4}$ for a constant $K\leq 1.087$. Consequently, the assumption $\sup_{1\leq m\leq N}\mathinner{\lvert x_{m}\rvert}\leq 2\sqrt{M_{N}}/\beta$ yields

[TABLE]

Combining this with Hölder’s inequality and $L^{2}(\mu)$-orthonormality of $\varphi_{n}$, that imply $\mu(\varphi_{n})\leq\mu(\varphi_{n}^{2})^{1/2}=1$, we obtain the bound

[TABLE]

Inserting this into Equation 16 produces

[TABLE]

Noticing that $\sqrt{\lambda}\operatorname{e}^{1/\beta^{2}}<1$ by Lemma 3 concludes the proof. ∎

Remark 1

From Lemma 3 we observe that the proof does not yield $\eta<1$ (for every $\ell$) if the assumption $\sup_{1\leq n\leq N}\mathinner{\lvert x_{n}\rvert}\leq 2\sqrt{M_{N}}/\beta$ on placement of the nodes is relaxed by replacing the constant $2$ on the right-hand side with $C>2$.

Consider now the $N$-point approximate Gaussian kernel quadrature rule ${\widetilde{Q}_{k,N}=\sum_{n=1}^{N}\widetilde{w}_{k,n}f(\tilde{x}_{n})}$ whose nodes and weights are defined in Theorem 2.2 and set $\alpha=1/\sqrt{2}$. The nodes $x_{n}^{\text{\tiny{GH}}}$ of the $N$-point Gauss–Hermite rule admit the bound (Area et al., 2004)

[TABLE]

for every $N\geq 1$. That is,

[TABLE]

Since the rule $\widetilde{Q}_{k,N}$ is exact for the first $N$ eigenfunctions, $M_{N}=N$. Hence the assumption on placement of the nodes in Theorem 3.1 holds. As our numerical experiments indicate that the weights $\widetilde{w}_{k,n}$ are positive and $\sum_{n=1}^{N}\mathinner{\lvert\widetilde{w}_{k,n}\rvert}\to 1$ as $N\to\infty$, it seems that the exponential convergence rate of Theorem 3.1 is valid for $\widetilde{Q}_{k,N}$ (as well as for the corresponding kernel quadrature rule $Q_{k,N}$) with $M_{N}=N$. Naturally, this result is valid whenever the growth of the absolute weight sum is, for example, polynomial in $N$.

Theorem 3.2

Let $\alpha=1/\sqrt{2}$ and suppose that $\sup_{N\geq 1}\sum_{n=1}^{N}\mathinner{\lvert\widetilde{w}_{k,n}\rvert}<\infty$. Then the quadrature rules $\widetilde{Q}_{k,N}(f)=\sum_{n=1}^{N}\widetilde{w}_{k,n}f(\tilde{x}_{n})$ and $Q_{k,N}(f)=\sum_{n=1}^{N}w_{k,n}f(\tilde{x}_{n})$ satisfy

[TABLE]

for $0<\eta<1$.

Another interesting case are the generalised Gaussian quadrature rules444Note that the cited results are for kernels and functions on compact intervals. However, generalisations for the whole real line are possible (Karlin and Studden, 1966, Chapter VI). for the eigenfunctions. As the eigenfunctions constitute a complete Chebyshev system (Kellog, 1918; Pinkus, 1996), there exists a quadrature rule $Q^{*}_{N}$ with positive weights $w_{1}^{*},\ldots,w_{N}^{*}$ such that $Q_{N}^{*}(\varphi_{n})=\mu(\varphi_{n})$ for every $n=0,\ldots,2N-1$ (Barrow, 1978). Appropriate control of the nodes of these quadrature rules would establish an exponential convergence result with the “double rate” $M_{N}=2N$.

4 Tensor product rules

Let $Q_{1},\ldots,Q_{d}$ be quadrature rules on $\mathbb{R}$ with nodes $X_{i}=\{x_{i,1},\ldots,x_{i,N_{i}}\}$ and weights $w_{1}^{i},\ldots,w_{N_{i}}^{i}$ for each $i=1,\ldots,d$. The tensor product rule on the Cartesian grid ${X\coloneqq X_{1}\times\cdots\times X_{d}\subset\mathbb{R}^{d}}$ is the cubature rule

[TABLE]

where $\mathcal{I}\in\mathbb{N}^{d}$ is a multi-index, $\mathcal{N}\coloneqq(N_{1},\ldots,N_{d})\in\mathbb{N}^{d}$, and the nodes and weights are

[TABLE]

We equip $\mathbb{R}^{d}$ with the $d$-variate standard Gaussian measure

[TABLE]

The following proposition is a special case of a standard result on exactness of tensor product rules (Oettershagen, 2017, Section 2.4).

Proposition 1

Consider the tensor product rule (19) and suppose that, for each $i=1,\ldots,d$, $Q_{i}(\varphi^{i}_{n})=\mu(\varphi^{i}_{n})$ for some functions $\varphi^{i}_{1},\ldots,\varphi_{N_{i}}^{i}\colon\mathbb{R}\to\mathbb{R}$. Then

[TABLE]

When a multivariate kernel is separable, this result can be used in constructing kernel cubature rules out of kernel quadrature rules. We consider $d$-dimensional separable Gaussian kernels

[TABLE]

where $\ell_{i}$ are dimension-wise length-scales. For each $i=1,\ldots,d$, the kernel quadrature rule $Q_{k,i}$ with nodes $X_{i}=\{x_{i,1},\ldots,x_{i,N_{i}}\}$ and weights $w_{k,1}^{i},\ldots,w_{k,N_{i}}^{i}$ is, by definition, exact for the $N_{i}$ kernel translates at the nodes:

[TABLE]

for each $n=1,\ldots,N_{i}$. Proposition 1 implies that the $d$-dimensional kernel cubature rule $Q_{k}^{d}$ at the nodes $X=X_{1}\times\cdots\times X_{d}$ is a tensor product of the univariate rules:

[TABLE]

with the weights being products of univariate Gaussian kernel quadrature weights, ${w_{k,\mathcal{I}}=\prod_{i=1}^{d}w_{k,\mathcal{I}(i)}}$. This is the case because each kernel translate $k^{d}(x,\cdot)$, $x\in X$, can be written as

[TABLE]

by separability of $k^{d}$.

We can extend Theorem 2.2 to higher dimensions if the node set is a Cartesian product of a number of scaled Gauss–Hermite node sets. For this purpose, for each $i=1,\ldots,d$ we use the $L(\mu_{\alpha_{i}})^{2}$-orthonormal eigendecomposition of the Gaussian kernel $k_{i}$. The eigenfunctions, eigenvalues, and other related constants from Section 2.1 for the eigendecomposition of the $i$th kernel are assigned an analogous subscript. Furthermore, use the notation

[TABLE]

Theorem 4.1

For $i=1,\ldots,d$, let $x_{i,1}^{\text{\tiny{GH}}},\ldots,x_{i,N_{i}}^{\text{\tiny{GH}}}$ and $w_{i,1}^{\text{\tiny{GH}}},\ldots w_{i,N_{i}}^{\text{\tiny{GH}}}$ stand for the nodes and weights of the $N_{i}$-point Gauss–Hermite quadrature rule and define the nodes

[TABLE]

Then the weights of the tensor product quadrature rule

[TABLE]

that is defined by the exactness conditions $\widetilde{Q}_{k}^{d}(\varphi_{\mathcal{I}})=\mu^{d}(\varphi_{\mathcal{I}})$ for every $\mathcal{I}\leq\mathcal{N}$, are ${\widetilde{w}_{k,\mathcal{I}}=\prod_{i=1}^{d}\widetilde{w}_{k,\mathcal{I}(i)}^{i}}$ for

[TABLE]

where $\alpha$, $\beta$, and $\delta$ are defined in Equation 5 and $\mathrm{H}_{2m}$ are the probabilists’ Hermite polynomials (4).

As in one dimension, the weights $\widetilde{w}_{k,\mathcal{I}}$ are supposed to approximate $w_{k,\mathcal{I}}$. Moreover, convergence rates can be obtained: a tensor product analogues of Theorems 3.1 and 3.2 follow from noting that every function $f\colon\mathbb{R}^{d}\to\mathbb{R}$ in the RKHS $\mathcal{H}^{d}$ of $k^{d}$ admits the multivariate Mercer expansion

[TABLE]

See (Kuo et al., 2017) for similar convergence analysis of tensor product Gauss–Hermite rules in $\mathcal{H}^{d}$.

Theorem 4.2

Let $\alpha_{1}=\cdots=\alpha_{d}=1/\sqrt{2}$. Suppose that the nodes $x_{i,1},\ldots,x_{i,N_{i}}$ and weights $w_{1}^{i},\ldots,w_{N_{i}}^{i}$ of the $N_{i}$-point quadrature rules $Q_{1,N_{1}},\ldots,Q_{d,N_{d}}$ satisfy

$\sup_{1\leq i\leq d}\sum_{n=1}^{N_{i}}\mathinner{\lvert w_{n}^{i}\rvert}\leq W_{\mathcal{N}}$* for some $W_{\mathcal{N}}\geq 1$;* 2. 2.

$Q_{i,N_{i}}(\varphi_{n}^{\alpha})=\mu(\varphi_{n}^{\alpha})$* for each $n=0,\ldots,M_{N_{i}}-1$ and $i=1,\ldots,d$ for some $M_{N_{i}}\geq 1$;* 3. 3.

$\sup_{1\leq n\leq{N_{i}}}\mathinner{\lvert x_{i,n}\rvert}\leq 2\sqrt{M_{N_{i}}}/\beta$* for each $i=1,\ldots,d$.*

Define the tensor product rule

[TABLE]

Then there exist constants $C>0$, independent of $\mathcal{N}$ and $Q_{\mathcal{N}}^{d}$, and $0<\eta<1$ such that

[TABLE]

where $M=\min(M_{N_{1}},\ldots,M_{N_{d}})$. Explicit forms of $C$ and $\eta$ appear in Equation 28.

Proof

The proof is largely analogous to that of Theorem 3.1. Since $f\in\mathcal{H}^{d}$ can be written as

[TABLE]

by defining the index set

[TABLE]

we obtain

[TABLE]

Consequently, the Cauchy–Schwarz inequality yields

[TABLE]

where we again use the notation

[TABLE]

Since $\mu(\varphi_{n})\leq 1$ for any $n\geq 0$, integration error for the eigenfunction $\varphi_{\mathcal{I}}$ satisfies

[TABLE]

Define the index sets $\mathcal{B}_{\mathcal{M}}^{j}(\mathcal{I})=\{j\leq i\leq d\,\mathrel{\mathop{\ordinarycolon}}\,\mathcal{I}(i)\geq M_{N_{i}}\}$ and their cardinalities $b_{\mathcal{M}}^{j}(\mathcal{I})=\#\mathcal{B}_{\mathcal{M}}^{j}(\mathcal{I})\leq d-j+1$ for $j\geq 1$. Because $\mathinner{\lvert\mu(\varphi_{\mathcal{I}(i)})-Q_{i,N_{i}}(\varphi_{\mathcal{I}(i)})\rvert}=0$ and $\mathinner{\lvert Q_{i,N_{i}}(\varphi_{\mathcal{I}(i)})\rvert}=\mathinner{\lvert\mu(\varphi_{\mathcal{I}(i)})\rvert}\leq 1$ if $\mathcal{I}(i)<M_{N_{i}}$, expansion of the recursive inequality (25) gives

[TABLE]

Equation 17 provides the bounds ${\mathinner{\lvert\mu(\varphi_{\mathcal{I}(i)})-Q_{i,N_{i}}(\varphi_{\mathcal{I}(i)})\rvert}\leq 1+K\sqrt{\beta}W_{\mathcal{N}}\operatorname{e}^{M_{N_{i}}/\beta^{2}}}$ and $\mathinner{\lvert Q_{i,N_{i}}(\varphi_{\mathcal{I}(i)})\rvert}\leq K\sqrt{\beta}W_{\mathcal{N}}\operatorname{e}^{M_{N_{i}}/\beta^{2}}$ for the constant $K=1.087$ that, when plugged in Equation 26, yield

[TABLE]

where the last inequality is based on the facts that $i\in\mathcal{B}_{\mathcal{M}}^{i}(\mathcal{I})$ if ${i\in\mathcal{B}_{\mathcal{M}}^{1}(\mathcal{I})}$ and ${1+K\sqrt{\beta}W_{\mathcal{N}}\operatorname{e}^{M_{N_{i}}/\beta^{2}}\leq 2K\sqrt{\beta}W_{\mathcal{N}}\operatorname{e}^{M_{N_{i}}/\beta^{2}}}$, a consequence of ${K,\beta,W_{\mathcal{N}}\geq 1}$. Equations 27 and 24, together with Lemma 3, now yield

[TABLE]

The claim therefore holds with

[TABLE]

∎

A multivariate version of Theorem 3.2 is obvious.

5 Numerical experiments

This section contains numerical experiments on properties and accuracy of the approximate Gaussian kernel quadrature weights defined in Theorems 2.2 and 4.1. The experiments have been implemented in MATLAB, and they are available at https://github.com/tskarvone/gauss-mercer. The value $\alpha=1/\sqrt{2}$ is used in all experiments. The experiments indicate that

Computation of the approximate weights in Equation 13 is numerically stable. 2. 2.

The weight approximation is quite accurate, its accuracy increasing with the number of nodes and the length-scale, as predicted in Section 2.5. 3. 3.

The weights $w_{k,n}$ and $\widetilde{w}_{k,n}$ are positive for every $N$ and $n=1,\ldots,N$ and their sums converge to one exponentially in $N$. 4. 4.

The quadrature rule $\widetilde{Q}_{k}$ converges exponentially, as implied by Theorem 3.2 and empirical observations on the behaviour of its weights. 5. 5.

In numerical integration of specific functions, the approximate kernel quadrature rule $\widetilde{Q}_{k}$ can achieve integration accuracy almost indistinguishable from that of the corresponding Gaussian kernel quadrature rule $Q_{k}$ and superior to some more traditional alternatives.

This suggest Equation 13 can be used as an accurate and numerically stable surrogate for computing the Gaussian kernel quadrature weights when the naive approach based on solving the linear system (2) is precluded by ill-conditioning of the kernel matrix. Furthermore, the choice (12) of the nodes by scaling the Gauss–Hermite nodes appears to yield an exponentially convergent kernel quadrature rule that has positive weights.

5.1 Numerical stability and distribution of weights

We have not encountered any numerical issues when computing the approximate weights (13). In this example we set $N=99$ and examine the distribution of approximate weights $\widetilde{w}_{k,n}$ for $\ell=0.05$, $\ell=0.4$ and $\ell=4$. Figure 1 depicts (i) approximate weights $\widetilde{w}_{k,n}$, (ii) absolute kernel quadrature weights $\mathinner{\lvert w_{k,n}\rvert}$ obtained by solving the linear system (2) for the points $\tilde{x}_{n}$ and, for $\ell=4$, (iii) Gauss–Hermite weights $w_{n}^{\text{\tiny{GH}}}$. The approximate weights $\widetilde{w}_{k,n}$ display no signs of numerical instabilities; their magnitudes vary smoothly and all of them are positive. That $\widetilde{w}_{k,1}>\widetilde{w}_{k,2}$ for $\ell=0.05$ appears to be caused by the sum in Equation 13 having not converged yet: the constant $2\alpha^{2}\beta^{2}/(1+2\delta^{2})-1$, that controls the rate of convergence of this sum, converges to $1$ as $\ell\to 0$ (in this case its value is $0.9512$) and $H_{2m}(x_{1}^{\text{\tiny{GH}}})>0$ for every $m=1,\ldots,49$ while $H_{2m}(x_{n}^{\text{\tiny{GH}}})<0$ for $m=46,47,48,49$. This and further experiments in Section 5.2 merely illustrates that quality of the weight approximation deteriorates when $\ell$ is small—as predicted in Section 2.5. Behaviour of $\widetilde{w}_{k,n}$ is in stark contrast to the naively computed weights $w_{k,n}$ that display clear signs of numerical instabilities for $\ell=0.4$ and $\ell=4$ (condition numbers of the kernel matrices were roughly $2.66\times 10^{16}$ and $3.59\times 10^{18}$). Finally, the case $\ell=4$ provides further evidence for numerical stability of Equation 13 since, based on Section 2.5, $\widetilde{w}_{k,n}\to w_{n}^{\text{\tiny{GH}}}$ as $\ell\to\infty$ and, furthermore, there is reason to believe that $w_{k,n}$ would share this property if they were computed in arbitrary-precision arithmetic. Section 5.3 and the experiments reported by Fasshauer and McCourt (2012) provide additional evidence for numerical stability of Equation 13.

5.2 Accuracy of the weight approximation

Next we assess quality of the weight approximation $\widetilde{w}_{k}\approx w_{k}$. Figure 2 depicts the results for a number of different length-scales in terms of norm of the relative weight error,

[TABLE]

As the kernel matrix quickly becomes ill-conditioned, computation of the kernel quadrature weights $w_{k}$ is challenging, particularly when the length-scale is large. To partially mitigate the problem we replaced the kernel quadrature weights with their QR decomposition approximations $\widetilde{w}_{k}^{M}$ derived in Section 2.2. The truncation length $M$ was selected based on machine precision; see (Fasshauer and McCourt, 2012, Section 4.2.2) for details. Yet even this does not work for large enough $N$. Because kernel quadrature rules on symmetric point sets have symmetric weights (Karvonen and Särkkä, 2018; Oettershagen, 2017, Section 5.2.4), breakdown in symmetricity of the computed kernel quadrature weights was used as a heuristic proxy for emergence of numerical instability: for each length-scale, relative errors are presented in Figure 2 until the first $N$ such that $\mathinner{\lvert 1-w_{k,N}/w_{k,1}\rvert}>10^{-6}$, ordering of the nodes being from smallest to the largest so that $w_{k,N}=w_{k,1}$ in absence of numerical errors.

5.3 Properties of the weights

Figure 3 shows the minimal weights $\min_{n=1,\ldots,N}\widetilde{w}_{k,n}$ and convergence to one of $\sum_{n=1}^{N}\mathinner{\lvert\widetilde{w}_{k,n}\rvert}$ for a number of different length-scales. These results provide strong numerical evidence for the conjecture that $\widetilde{w}_{k,n}$ remain positive and that the assumptions of Theorem 3.2 hold. Exact weights, as long as they can be reliably computed (see Section 5.2), exhibit behaviour practically indistinguishable from the approximate ones and are not therefore depicted separately in Figure 3.

5.4 Worst-case error

The worst-case error $e(Q)$ of a quadrature rule $Q(f)=\sum_{n=1}^{N}w_{n}f(x_{n})$ in a reproducing kernel Hilbert space induced by the kernel $k$ is explicitly computable:

[TABLE]

Figure 4 compares the worst-case errors in the RKHS of the Gaussian kernel for six different length-scales of (i) the classical Gauss–Hermite quadrature rule, (ii) the quadrature $\widetilde{Q}_{k}(f)=\sum_{n=1}^{N}\widetilde{w}_{k,n}f(\tilde{x}_{n})$ of Theorem 2.2, and (iii) the kernel quadrature rule with its nodes placed uniformly between the largest and smallest of $\tilde{x}_{n}$. We observe that $\widetilde{Q}_{k}$ is, for all length-scales, the fastest of these rules to converge (the kernel quadrature rule at $\tilde{x}_{n}$ yields WCEs practically indistinguishable from those of $\widetilde{Q}_{k}$ and is therefore not included). It also becomes apparent that the convergence rates derived in Theorems 3.1 and 3.2 for $\widetilde{Q}_{k}$ are rather conservative. For example, for $\ell=0.2$ and $\ell=1$ the empirical rates are $e(\widetilde{Q}_{k})=\mathcal{O}(\operatorname{e}^{-cN})$ with $c\approx 0.21$ and $c\approx 0.98$, respectively, whereas Equation 18 yields the theoretical values $c\approx 0.00033$ and $c\approx 0.054$, respectively.

5.5 Numerical integration

Set $\ell=1.2$ and consider the integrand

[TABLE]

When $0<c_{i}<4$ and $m_{i}\in\mathbb{N}$ for each $i=1,\ldots,d$, the function is in $\mathcal{H}$ (Minh, 2010, Theorems 1 and 3). Furthermore, the Gaussian integral of this function is available in closed form:

[TABLE]

when $m_{i}$ are even (when they are not even, the integral is obviously zero). Figure 5 shows integration error of the three methods (or, in higher dimensions, their tensor product versions) used in Section 5.4 and the kernel quadrature rule based on the nodes $\tilde{x}_{n}$ for (i) $d=1$, $m_{1}=6$, $c_{1}=3/2$ and (ii) $d=3$, $m_{1}=6$, $m_{2}=4$, $m_{3}=2$, $c_{1}=3/2$, $c_{2}=3$, $c_{3}=1/2$. As expected, there is little difference between $\widetilde{Q}_{k}$ and $Q_{k}$.