Lattice Rules Meet Kernel Cubature
January 16, 2025
Abstract
Rank-1 lattice rules are a class of equally weighted quasi-Monte Carlo methods that achieve essentially linear convergence rates for functions in a reproducing kernel Hilbert space (RKHS) characterized by square-integrable first-order mixed partial derivatives. In this work, we explore the impact of replacing the equal weights in lattice rules with optimized cubature weights derived using the reproducing kernel. We establish a theoretical result demonstrating a doubled convergence rate in the one-dimensional case and provide numerical investigations of convergence rates in higher dimensions. We also present numerical results for an uncertainty quantification problem involving an elliptic partial differential equation with a random coefficient.
1 Introduction↩︎
Computing the expected value \(\mathbb{E}_{\mathbb{P}}[f]\) of a function \(f\colon D \to \mathbb{R}\) over a domain \(D \subseteq \mathbb{R}^s\) with respect to a probability distribution \(\mathbb{P}\) is a fundamental problem in fields such as uncertainty quantification, machine learning, statistics, financial mathematics, and statistical mechanics. Since these integrals are often intractable analytically, they are approximated numerically using an empirical mean: \[\label{KaKlScSu:equ:empirical95mean95approximation} \mathbb{E}_{\mathbb{P}}[f] = \int_D f(x)\,{\rm d}\mathbb{P}(x) \approx \sum_{k=0}^{n-1} w_k f({\boldsymbol{t}}_k) = {\mathbb{E}}_{\mathbb{P}_{{\boldsymbol{T}}}^{\boldsymbol{w}}}[f], \qquad {\mathbb{P}}_{{\boldsymbol{T}}}^{{\boldsymbol{w}}} := \sum_{k=0}^{n-1} w_{k} \delta_{{\boldsymbol{t}}_{k}}.\tag{1}\] The central challenge in constructing higher-order cubature methods lies in the careful selection of evaluation points \({\boldsymbol{T}}= ({\boldsymbol{t}}_k)_{k=0}^{n-1} \in D^n\) and weights \({\boldsymbol{w}}= (w_k)_{k=0}^{n-1} \in {\mathbb{R}}^n\) to ensure favorable approximation properties of the error \(|\mathbb{E}_{\mathbb{P}}[f] - \mathbb{E}_{\mathbb{P}_{{\boldsymbol{T}}}^{\boldsymbol{w}}}[f]|\).
To address this challenge for potentially high-dimensional integration problems, one can either use sampling-based methods or numerical cubature rules. Sampling-based approaches include methods like Markov chain Monte Carlo (MCMC), which construct a Markov chain with the target distribution \(\mathbb{P}\) as its stationary distribution, and importance sampling, which modifies the probability measure to reduce variance and often enables direct sampling. In contrast, numerical cubature methods such as sparse grids and quasi-Monte Carlo (QMC) methods construct the nodes \({\boldsymbol{T}}\) and weights \({\boldsymbol{w}}\) deterministically and can achieve faster convergence rates under sufficient smoothness assumptions on the integrand. While QMC is fundamentally deterministic, in our numerical experiments (Section 4) we use a randomized QMC method based on several random shifts of a fixed lattice which is common practice in QMC. We focus on a particularly simple QMC rule in which \({\boldsymbol{T}}\) is chosen to be a lattice, as introduced in Section 2.4.
In this work, we consider cubature rules with a known convergence rate in a reproducing kernel Hilbert space (RKHS) \(\mathcal{H}\). We interpret the cubature rule as an element of the subspace \(V_{{\boldsymbol{T}}} = \mathrm{span}(K({\boldsymbol{t}}_0,\setbox 0x to\wd 0{\hss\cdot\hss}), \ldots, K({\boldsymbol{t}}_{n-1},\setbox 0x to\wd 0{\hss\cdot\hss})) \subseteq \mathcal{H}\), spanned by the reproducing kernel with its first argument fixed at the cubature nodes \({\boldsymbol{t}}_k\).
Specifically, we focus on the kernel mean embeddings of \(\mathbb{P}\) and \(\mathbb{P}_{{\boldsymbol{T}}}^{\boldsymbol{w}^{\ast}}\): \[h = \int_D K({\boldsymbol{x}}, \setbox 0x to\wd 0{\hss\cdot\hss}) \, \mathrm{d}\mathbb{P}({\boldsymbol{x}}) \in \mathcal{H}, \qquad h_{{\boldsymbol{T}}}^{{\boldsymbol{w}}^{\ast}} = \sum_{k=0}^{n-1} w_k^* K({\boldsymbol{t}}_k, \setbox 0x to\wd 0{\hss\cdot\hss}) \in V_{{\boldsymbol{T}}} \subseteq \mathcal{H},\] and choose the weights \({\boldsymbol{w}}^*\) so that \(h_{{\boldsymbol{T}}}^{{\boldsymbol{w}}^{\ast}}\) becomes the \(\mathcal{H}\)-orthogonal projection (or equivalently, the kernel interpolant, cf.Lemma 1) of \(h\) onto \(V_{{\boldsymbol{T}}}\). We refer to the resulting cubature rule \(Q_{{\boldsymbol{T}}}^{{\boldsymbol{w}}^{\ast}} f := \sum_{k=0}^{n-1} w_{k}^{\ast} f({\boldsymbol{t}}_{k})\) as kernel cubature. By the reproducing-type properties (cf.@eq:KaKlScSu:equ:KME95reproducing95property below) of \(h\) and \(h_{{\boldsymbol{T}}}^{{\boldsymbol{w}}^{\ast}}\), \[\mathbb{E}_{\mathbb{P}}[f] = \langle f, h \rangle_{\mathcal{H}}, \qquad \mathbb{E}_{\mathbb{P}_{{\boldsymbol{T}}}^{{\boldsymbol{w}}^{\ast}}}[f] = \langle f, h_{{\boldsymbol{T}}}^{{\boldsymbol{w}}^{\ast}} \rangle_{\mathcal{H}},\] the choice of \({\boldsymbol{w}}^*\) ensures that these two expected values are close, yielding a cubature rule as in 1 with favorable approximation properties. In fact, it guarantees that the worst-case error \[\label{KaKlScSu:equ:WCE95general95form95in95RKHS} e(Q_{{\boldsymbol{T}}}^{{\boldsymbol{w}}} , {\mathcal{H}}) := \sup_{\| f \|_{{\mathcal{H}}}=1} | \mathbb{E}_{\mathbb{P}}[f] - Q_{{\boldsymbol{T}}}^{{\boldsymbol{w}}} f| = \sup_{\| f \|_{{\mathcal{H}}}=1} \langle h - h_{{\boldsymbol{T}}}^{{\boldsymbol{w}}} , f \rangle_{{\mathcal{H}}} = \| h - h_{{\boldsymbol{T}}}^{{\boldsymbol{w}}} \|_{{\mathcal{H}}}\tag{2}\] is minimized by \({\boldsymbol{w}}^{\ast}\), where \(Q_{{\boldsymbol{T}}}^{{\boldsymbol{w}}} f := {\mathbb{E}}_{{\mathbb{P}}_{{\boldsymbol{T}}}^{{\boldsymbol{w}}}}[f]\). While conventional QMC analysis primarily focuses on the worst-case error, we aim to achieve an additional improvement in the convergence of the approximation error, motivated by the observation (cf.Proposition [KaKlScSu:prop:approximation95error95additional95gain]) that, for a fixed \(f \in {\mathcal{H}}\), \[\label{KaKlScSu:equ:approximation95error95additional95gain} |\mathbb{E}_{\mathbb{P}}[f] - Q_{{\boldsymbol{T}}}^{{\boldsymbol{w}}^{\ast}}f| \leq e(Q_{{\boldsymbol{T}}}^{{\boldsymbol{w}}^{\ast}}, {\mathcal{H}})\, \mathrm{dist}_{{\mathcal{H}}}(f,V_{{\boldsymbol{T}}}),\tag{3}\] a result that is unique to the optimally weighted cubature rule \(Q_{{\boldsymbol{T}}}^{{\boldsymbol{w}}^{\ast}}\). This expectation arises from the intuition that the distance \(\mathrm{dist}_{{\mathcal{H}}}(f,V_{{\boldsymbol{T}}}) = \inf_{v \in V_{{\boldsymbol{T}}}} \|{f-v}\|_{{\mathcal{H}}}\) between \(f\) and \(V_{{\boldsymbol{T}}}\) should decrease as \(n\) increases and \(V_{{\boldsymbol{T}}}\) increasingly approximates \({\mathcal{H}}\). However, we cannot yet establish a precise rate for this decay.
The optimal weights \({\boldsymbol{w}}^{\ast}\) coincide with Bayesian cubature weights for the canonical choice of the prior [1]–[5]. The recent work of Hickernell and Jagadeeswaran [6]–[8] has investigated the construction of optimized cubature weights for fixed sequences of lattice points and Sobol\('\) nets, but only within the context of shift-invariant kernels and Walsh kernels. We note that our construction is related to recent studies on kernel interpolation over lattice point sets [9]–[11] in the sense that the kernel cubature of a function \(f \in {\mathcal{H}}\) is equivalent to computing the integral of its kernel interpolant over the cubature point set. However, our work addresses the non-periodic setting, while the works [9]–[11] only discuss kernel interpolation over lattice point sets in the periodic setting.
Contributions. We make the following contributions to weighted QMC methods:
We propose a weighted version of QMC cubature, where the weights minimize the distance between the kernel mean embeddings of the empirical and true distribution, leading to low approximation error, as observed in various numerical experiments across low and high dimensions.
While, by construction, the worst-case error \(e(Q_{{\boldsymbol{T}}}^{{\boldsymbol{w}}^{\ast}}, {\mathcal{H}})\) is minimal among all possible weights \({\boldsymbol{w}}\), it shows only a slight improvement over the equally weighted case. However, we observe a significant reduction in the approximation error \(|\mathbb{E}_{\mathbb{P}}[f] - \mathbb{E}_{\mathbb{P}_{{\boldsymbol{T}}}^{\boldsymbol{w}}}[f]|\), which we attribute to the factor \(\mathrm{dist}_{{\mathcal{H}}}(f,V_{{\boldsymbol{T}}})\) in the bound 3 . Unlike the constant factor \(\| f \|_{{\mathcal{H}}}\) in the classical bound for equally weighted QMC, \(\mathrm{dist}_{{\mathcal{H}}}(f,V_{{\boldsymbol{T}}})\) can be expected to decrease as \(n\) increases.
Although the explicit error convergence rates for optimally weighted lattice point sets remain an open problem in higher dimensions, we prove in the one-dimensional case that this approach leads to a doubled rate of convergence compared to the equally weighted case.
We numerically investigate the behavior of the worst-case error in Sobolev spaces of higher smoothness \(\alpha = 4\) using weights optimized for the less smooth setting \(\alpha = 2\). A significant improvement in convergence speed is observed. This experiment is conducted in low dimensions only (\(s=2\)) with tent-transformed lattices, which are known to enhance the convergence rate from first to second order.
Outline. This document is structured as follows. After introducing our setup and notation in Section 2, we provide our theoretical contributions in Section 3. In Section 4 we numerically demonstrate the improvement of kernel cubature over equally weighted lattice rules and provide a conclusion in Section 5.
2 Preliminaries and Notation↩︎
Throughout this manuscript, we will use the following general notation: \(D\subseteq {\mathbb{R}}^{s}\), \(s\in{\mathbb{N}}\), will be the domain of interest equipped with its Borel \(\sigma\)-algebra and a probability measure \({\mathbb{P}}\), typically \(D = [0,1]^{s}\) with uniform measure \({\mathbb{P}}= \mathsf{Unif}_{D}\). We denote by \(\{ v \} := (v_{j} - \lfloor v_{j} \rfloor)_{j=1,\dots,s}\) the componentwise fractional part of a vector \(v \in \mathbb{R}^{s}\), by \(\mathbf{1} = (1)_{k=0}^{n-1}\) the \(n\)-dimensional unit vector, and by \(\mathbb{1} \colon D \to {\mathbb{R}}\) the constant unit function. Further, for a function \(f \in {\mathcal{H}}\) we denote \[\begin{align} If &= \int_{D} f \, \mathrm{d} {\mathbb{P}}= {\mathbb{E}}_{{\mathbb{P}}}[f], && \\ Q_{{\boldsymbol{T}}} f &= n^{-1} \sum_{k=0}^{n-1} f({\boldsymbol{t}}_{k}) = {\mathbb{E}}_{{\mathbb{P}}_{{\boldsymbol{T}}}}[f], & {\mathbb{P}}_{{\boldsymbol{T}}} &= n^{-1} \sum_{k=0}^{n-1} \delta_{{\boldsymbol{t}}_{k}}, \\ Q_{{\boldsymbol{T}}}^{{\boldsymbol{w}}} f &= \sum_{k=0}^{n-1} w_{k} f({\boldsymbol{t}}_{k}) = {\mathbb{E}}_{{\mathbb{P}}_{{\boldsymbol{T}}}^{{\boldsymbol{w}}}}[f], & {\mathbb{P}}_{{\boldsymbol{T}}}^{{\boldsymbol{w}}} &= \sum_{k=0}^{n-1} w_{k} \delta_{{\boldsymbol{t}}_{k}}, \end{align}\] where the cubature rules \(Q_{{\boldsymbol{T}}}\) and \(Q_{{\boldsymbol{T}}}^{{\boldsymbol{w}}}\) are based on the evaluation points \({\boldsymbol{T}}= ({\boldsymbol{t}}_{k})_{k=0}^{n-1} \in D^{n}\), which in this work will be a (potentially shifted and tent-transformed, cf.Section 2.4) lattice, and cubature weights, \({\boldsymbol{w}}= (w_{k})_{k=0}^{n-1} \in {\mathbb{R}}^{n}\). Here, \({\mathbb{P}}_{{\boldsymbol{T}}}\) and \({\mathbb{P}}_{{\boldsymbol{T}}}^{{\boldsymbol{w}}}\) denote the corresponding (possibly signed) discrete measures on \(D\). Note that we make no assumptions on the cubature weights to be non-negative or to sum to one. This is in line with common practice in Bayesian cubature, where the primary objective is optimization rather than strictly enforcing a probabilistic interpretation of the weights. Consequently, the resulting measures, while being finite, may fail to be probability measures, can attain negative values, and may introduce a bias.
2.1 Reproducing Kernel Hilbert Spaces and Kernel Mean Embeddings↩︎
The (possibly signed) measures \({\mathbb{P}},{\mathbb{P}}_{{\boldsymbol{T}}},{\mathbb{P}}_{{\boldsymbol{T}}}^{{\boldsymbol{w}}}\) will be embedded into a reproducing kernel Hilbert space (RKHS; [12]) \({\mathcal{H}}\) corresponding to a symmetric and positive definite kernel \(K \colon D\times D \to {\mathbb{R}}\). Note that we work with strictly positive definite kernels rather than semi-positive definite ones in the sense that the Gram matrix \(G=(K({\boldsymbol{x}}_i,{\boldsymbol{x}}_j))_{i,j=1}^N\) is (strictly) positive definite, and thereby invertible, for all \(N\in\mathbb{N}\) and pairwise distinct \({\boldsymbol{x}}_i\in D\), \(i=1,\dots,N\). For a signed measure \(\mu\) on \(D\) its kernel mean embedding (\(\mathsf{KME}\)) is defined by \[\mathsf{KME}(\mu) := \int_{D} K({\boldsymbol{x}},\setbox 0x to\wd 0{\hss\cdot\hss}) \, \mathrm{d} \mu({\boldsymbol{x}}) \in {\mathcal{H}}.\] Strictly speaking, the KME is defined only for certain combinations of kernels and signed measures [12], in particular, the corresponding integral must be well defined. We omit these technical details here, as the assumptions are always fulfilled for the measures and kernels considered in this paper. Importantly, the KME satisfies a reproducing-type property [13] \[\label{KaKlScSu:equ:KME95reproducing95property} \langle \mathsf{KME}(\mu) , f \rangle_{{\mathcal{H}}} = {\mathbb{E}}_{\mu}[f], \qquad f \in {\mathcal{H}}.\tag{4}\] After defining \(V_{{\boldsymbol{T}}} = \mathrm{span}(K({\boldsymbol{t}}_0,\setbox 0x to\wd 0{\hss\cdot\hss}), \ldots, K({\boldsymbol{t}}_{n-1},\setbox 0x to\wd 0{\hss\cdot\hss})) \subseteq \mathcal{H}\) and denoting by \(P_{V_{{\boldsymbol{T}}}}\colon {\mathcal{H}}\to V_{{\boldsymbol{T}}}\) the corresponding \({\mathcal{H}}\)-orthogonal projection, the following embeddings will be crucial: \[\begin{align} h &= \mathsf{KME}({\mathbb{P}}) = \int_{D} K({\boldsymbol{x}},\setbox 0x to\wd 0{\hss\cdot\hss}) {\mathbb{P}}(\mathrm d {\boldsymbol{x}}) \in {\mathcal{H}}, \\ h_{{\boldsymbol{T}}} &= \mathsf{KME}({\mathbb{P}}_{{\boldsymbol{T}}}) = \int_{D} K({\boldsymbol{x}},\setbox 0x to\wd 0{\hss\cdot\hss}) {\mathbb{P}}_{{\boldsymbol{T}}}(\mathrm d{\boldsymbol{x}}) = n^{-1} \sum_{k=0}^{n-1} K({\boldsymbol{t}}_{k},\setbox 0x to\wd 0{\hss\cdot\hss}) \in V_{{\boldsymbol{T}}} \subseteq {\mathcal{H}}, \\ h_{{\boldsymbol{T}}}^{{\boldsymbol{w}}} &= \mathsf{KME}({\mathbb{P}}_{{\boldsymbol{T}}}^{{\boldsymbol{w}}}) = \int_{D} K({\boldsymbol{x}},\setbox 0x to\wd 0{\hss\cdot\hss}) {\mathbb{P}}_{{\boldsymbol{T}}}^{{\boldsymbol{w}}}(\mathrm d{\boldsymbol{x}}) = \sum_{k=0}^{n-1} w_{k} K({\boldsymbol{t}}_{k},\setbox 0x to\wd 0{\hss\cdot\hss}) \in V_{{\boldsymbol{T}}} \subseteq {\mathcal{H}}. \end{align}\]
2.2 Kernel Cubature↩︎
The basic idea of this paper is to view \(h\) and \(h_{{\boldsymbol{T}}}^{{\boldsymbol{w}}}\) as representing the integration operator and the cubature rule from 1 as elements in \({\mathcal{H}}\)—after all, by 4 , \(\mathbb{E}_{\mathbb{P}}[f] = \langle h , f \rangle_{{\mathcal{H}}}\) and \({\mathbb{E}}_{\mathbb{P}_{{\boldsymbol{T}}}^{\boldsymbol{w}}}[f] = \langle h_{{\boldsymbol{T}}}^{\boldsymbol{w}} , f \rangle_{{\mathcal{H}}}\). Hence, in order to reduce the approximation error in 1 , it seems natural to choose the weights \({\boldsymbol{w}}\) such that \(h_{{\boldsymbol{T}}}^{{\boldsymbol{w}}}\) is the best approximation of \(h\) in \({\mathcal{H}}\), that is, its orthogonal projection onto \(V_{{\boldsymbol{T}}}\). It is well-known [14] that such orthogonal projections within RKHSs correspond to interpolation:
Lemma 1. Let \({\boldsymbol{T}}= ({\boldsymbol{t}}_k)_{k=0}^{n-1} \in D^n\) be any point set in \(D\). The orthogonal projection \(\hat{g} = P_{V_{{\boldsymbol{T}}}} g\) of each \(g \in {\mathcal{H}}\) onto \(V_{{\boldsymbol{T}}} = \mathrm{span}(K({\boldsymbol{t}}_0,\setbox 0x to\wd 0{\hss\cdot\hss}), \ldots, K({\boldsymbol{t}}_{n-1},\setbox 0x to\wd 0{\hss\cdot\hss})) \subseteq \mathcal{H}\) coincides with the unique solution of the following interpolation problem: find \(\hat{g} \in V_{{\boldsymbol{T}}}\) such that \(\hat{g}({\boldsymbol{t}}_{k}) = g({\boldsymbol{t}}_{k})\) for \(k = 0,\dots,n-1\). In particular, \(h_{{\boldsymbol{T}}}^{{\boldsymbol{w}}^{\ast}} = P_{V_{{\boldsymbol{T}}}} h\) is given by the unique solution \({\boldsymbol{w}}^{\ast}\) of \[\label{KaKlScSu:eq:gramsystem} {\mathcal{K}}_{{\boldsymbol{T}}} {\boldsymbol{w}}^{\ast} = (h({\boldsymbol{t}}_{k}))_{k=0}^{n-1},\qquad{(1)}\] where \({\mathcal{K}}_{{\boldsymbol{T}}} = (K({\boldsymbol{t}}_{k},{\boldsymbol{t}}_{\ell}))_{k,\ell=0}^{n-1}\) is the Gram matrix.
Definition 1. We refer to the cubature weights \({\boldsymbol{w}}^{\ast}\) given by ?? as optimal weights* and to the corresponding weighted cubature rule \(Q_{{\boldsymbol{T}}}^{{\boldsymbol{w}}^{\ast}}\) as kernel cubature.*
Another advantage of the embeddings \(h,h_{{\boldsymbol{T}}}^{{\boldsymbol{w}}^{\ast}}\) is that the norm of their difference naturally describes the worst-case cubature error 2 . Since the optimal weights \({\boldsymbol{w}}^{\ast}\) stem from an orthogonal projection, the Pythagorean theorem implies \[\label{KaKlScSu:equ:optimal95wce95formula} e(Q_{{\boldsymbol{T}}}^{{\boldsymbol{w}}^{\ast}} , {\mathcal{H}})^{2} = \| h - h_{{\boldsymbol{T}}}^{{\boldsymbol{w}}^{\ast}} \|_{{\mathcal{H}}}^{2} = \| h \|_{{\mathcal{H}}}^{2} - \| h_{{\boldsymbol{T}}}^{{\boldsymbol{w}}^{\ast}} \|_{{\mathcal{H}}}^{2} = \| h \|_{{\mathcal{H}}}^{2} - (w^{\ast})^{\top} {\mathcal{K}}_{{\boldsymbol{T}}} w^{\ast}.\tag{5}\] We will see in the next subsection that, for the kernels \(K\) considered in this paper, \(h = \mathbb{1}\) turns out to be the constant unit function and 5 reduces to \(e(Q_{{\boldsymbol{T}}}^{{\boldsymbol{w}}} , {\mathcal{H}})^{2} = 1 - \sum_{k=0}^{n-1} w_{k}^{\ast}\), making it easily computable once the weights are established.
2.3 Specific Kernels and the Corresponding Sobolev Spaces↩︎
In this paper, we consider the Sobolev space of dominating-mixed smoothness, see [15] and [16].
Definition 2 (Sobolev space of dominating mixed smoothness). Let \(s\in\mathbb{Z}_+\) and let \(\boldsymbol{\gamma}=(\gamma_{\mathrm{\mathfrak{u}}})_{\mathrm{\mathfrak{u}}\subseteq\{1,\ldots,s\}}\) be a sequence of positive weights, termed coordinate weights. The weighted Sobolev space \({\mathcal{H}}^{\alpha}_{s,\boldsymbol{\gamma}}\) of order \(\alpha \in\mathbb{Z}_+\) is a reproducing kernel Hilbert space with inner product \[\begin{align} &\langle f , g \rangle _{{\mathcal{H}}^{\alpha}_{s,\boldsymbol{\gamma}}} := \left(\int_{[0,1]^s} f({\boldsymbol{x}}){\mathrm{d}}{\boldsymbol{x}}\right)\; \left(\int_{[0,1]^s} g({\boldsymbol{x}}){\mathrm{d}}{\boldsymbol{x}}\right) + \\ &\sum_{\emptyset\ne\mathrm{\mathfrak{u}}\subseteq\{1,\ldots,s\}} \kern-5mm \gamma^{-1}_{\mathrm{\mathfrak{u}}} \sum_{ \substack{ {\boldsymbol{\tau}}\in \{0,\ldots,\alpha \}^{|\mathfrak{u}|} \\ {\mathfrak{v}}:= \{j: \tau_j=\alpha \}} } \int_{[{\boldsymbol{0}}_{\mathfrak{v}},{\boldsymbol{1}}_{\mathfrak{v}}]} \left(\int_{[{\boldsymbol{0}}_{-{\mathfrak{v}}},{\boldsymbol{1}}_{-{\mathfrak{v}}}]} \kern-8mm f^{({\boldsymbol{\tau}},{\boldsymbol{0}}_{-\mathrm{\mathfrak{u}}})} ({\boldsymbol{x}}) {\mathrm{d}}{\boldsymbol{x}}_{-{\mathfrak{v}}} \right) \left(\int_{[{\boldsymbol{0}}_{-{\mathfrak{v}}},{\boldsymbol{1}}_{-{\mathfrak{v}}}]} \kern-8mm g^{({\boldsymbol{\tau}},{\boldsymbol{0}}_{-\mathrm{\mathfrak{u}}})} ({\boldsymbol{x}}) {\mathrm{d}}{\boldsymbol{x}}_{-{\mathfrak{v}}} \right) \kern-0.7mm {\mathrm{d}}{\boldsymbol{x}}_{\mathfrak{v}} , \end{align}\] where for \({\mathfrak{v}}\subseteq \{1,\ldots,s\}\) we write \([{\boldsymbol{a}}_{\mathfrak{v}},{\boldsymbol{b}}_{\mathfrak{v}}] := \prod_{j\in {\mathfrak{v}}} [a_j,b_j]\) and likewise for \(-{\mathfrak{v}}:= \{1,\ldots,s\} \setminus {\mathfrak{v}}\), \([{\boldsymbol{a}}_{-{\mathfrak{v}}},{\boldsymbol{b}}_{-{\mathfrak{v}}}] = \prod_{j \in -{\mathfrak{v}}} [a_j,b_j] := \prod_{j \in \{1,\ldots,s\} \setminus {\mathfrak{v}}} [a_j,b_j]\). The reproducing kernel is given by \[\begin{align} \label{KaKlScSu:eq:kern} K^{\alpha}_{s,\boldsymbol{\gamma}} ({\boldsymbol{x}},{\boldsymbol{y}}) &:= 1+ \sum_{\emptyset\ne\mathrm{\mathfrak{u}}\subseteq\{1,\ldots,s\}} \gamma_{\mathrm{\mathfrak{u}}}\prod_{j\in \mathrm{\mathfrak{u}}} \left( -1 + K^{\alpha}_{1,1}(x_j,y_j)\right) , \end{align}\qquad{(2)}\] with \[\begin{align} K^{\alpha}_{1,1}(x_j,y_j)= 1+\sum_{\tau=1}^\alpha \frac{B_{\tau}(x_j)}{\tau!} \, \frac{B_{\tau}(y_j)}{\tau!} + (-1)^{\alpha+1} \frac{\widetilde{B}_{2\alpha}(x_j-y_j)}{(2\alpha)!}, \end{align}\] where \(\widetilde{B}_{2\alpha}\) is the \(1\)-periodic Bernoulli polynomial of order \(2\alpha\), i.e., denoting by \(\{x-y\}\) the fractional part of \(x-y\) and using the standard Bernoulli polynomial \(B_{2\alpha}\) we define \[\widetilde{B}_{2\alpha}(x-y):=B_{2\alpha}(\{x-y\}).\]
The term “coordinate weights” is used to help distinguish the weights \(\boldsymbol{\gamma}\) from the cubature weights \({\boldsymbol{w}}\).
We remark that, since Bernoulli polynomials \(B_{\alpha}(x),\; \alpha\ge 1,\) integrate to \(0\) over the unit interval \([0,1]\), we have \[\int_0 ^1 K^{\alpha}_{1,1}(x_j,\setbox 0x to\wd 0{\hss\cdot\hss}) \, {\rm d} x_j = \mathbb{1}, \quad \text{ hence,} \quad \int_{[0,1]^s} K^{\alpha}_{s,\boldsymbol{\gamma}} ({\boldsymbol{x}},\setbox 0x to\wd 0{\hss\cdot\hss}) \, {\rm d} {\boldsymbol{x}} = \mathbb{1}.\] Due to this fact, we can represent the worst-case error 2 by \[\begin{align} \label{KaKlScSu:eq:wce95general} \begin{aligned} e(Q_{{\boldsymbol{T}}}^{{\boldsymbol{w}}} , {\mathcal{H}}^{\alpha}_{s,\boldsymbol{\gamma}})^{2} &= \langle \mathbb{1} , \mathbb{1} \rangle_{{\mathcal{H}}^{\alpha}_{s,\boldsymbol{\gamma}}} - 2 \sum_{k=0}^{n-1}w_k \langle \mathbb{1} , K^{\alpha}_{s,\boldsymbol{\gamma}}({\boldsymbol{t}}_k,\setbox 0x to\wd 0{\hss\cdot\hss}) \rangle_{{\mathcal{H}}^{\alpha}_{s,\boldsymbol{\gamma}}} \\ & + \sum_{k=0}^{n-1}\sum_{k'=0}^{n-1} w_{k}w_{k'} \langle K^{\alpha}_{s,\boldsymbol{\gamma}} ({\boldsymbol{t}}_k,\setbox 0x to\wd 0{\hss\cdot\hss}) , K^{\alpha}_{s,\boldsymbol{\gamma}} ({\boldsymbol{t}}_{k'},\setbox 0x to\wd 0{\hss\cdot\hss}) \rangle_{{\mathcal{H}}^{\alpha}_{s,\boldsymbol{\gamma}}} \\ &= 1-2\sum_{k=0}^{n-1}w_k+\sum_{k=0}^{n-1}\sum_{k'=0}^{n-1} w_{k}w_{k'} K^{\alpha}_{s,\boldsymbol{\gamma}} ({\boldsymbol{t}}_k,{\boldsymbol{t}}_{k'}), \end{aligned} \end{align}\tag{6}\] which for equal and optimal weights reduces to \[\label{KaKlScSu:eq:wce-weights} e(Q_{{\boldsymbol{T}}}^{{\boldsymbol{w}}} , {\mathcal{H}}^{\alpha}_{s,\boldsymbol{\gamma}})^{2} = \begin{cases} -1 + n^{-2} \sum_{k,k'=0}^{n-1} K^{\alpha}_{s,\boldsymbol{\gamma}} ({\boldsymbol{t}}_k,{\boldsymbol{t}}_{k'}) & \text{if } {\boldsymbol{w}}= (n^{-1})_{k=0}^{n-1}, \\ 1 - \sum_{k=0}^{n-1} w_{k}^{\ast} & \text{if } {\boldsymbol{w}}= {\boldsymbol{w}}^{\ast}. \end{cases}\tag{7}\]
2.4 Lattice Points and the Tent Transform↩︎
Let \(\tilde{{\boldsymbol{T}}} = (\tilde{{\boldsymbol{t}}}_k)_{k=0}^{n-1} \subset [0,1]^s\) be a lattice point set defined by \[\tilde{{\boldsymbol{t}}}_k = \left\{ \frac{k \boldsymbol{z}}{n} \right\}, \quad k = 0, \ldots, n - 1,\] where \(\{\setbox 0x to\wd 0{\hss\cdot\hss}\}\) denotes the componentwise fractional part and \(\boldsymbol{z}\) is the so-called generating vector, which consists of \(s\) elements of integers, \(\boldsymbol{z}\in\{1, \ldots, n-1\}^s\). In Section 4, our cubature rules \(Q_{{\boldsymbol{T}}}\) and \(Q_{{\boldsymbol{T}}}^{{\boldsymbol{w}}}\) will be based on a point set \({\boldsymbol{T}}= ({\boldsymbol{t}}_k)_{k=0}^{n-1} \subset [0,1]^s\) defined in one of the following ways:
as the unshifted lattice point set itself, \({\boldsymbol{T}}= \tilde{{\boldsymbol{T}}}\);
as a randomly shifted lattice point set, defined by \({\boldsymbol{t}}_k = \{ \tilde{{\boldsymbol{t}}}_k + \Delta \}\), where \(\Delta\) is a fixed random shift sampled uniformly from \([0,1]^s\);
or as a randomly shifted and tent-transformed lattice point set, defined by \({\boldsymbol{t}}_k = \phi(\{ \tilde{{\boldsymbol{t}}}_k + \Delta \})\), where the so-called baker’s transform is applied componentwise: \[\phi(\boldsymbol{t}) := (\phi(t_1), \ldots, \phi(t_s)), \quad \phi(t) := 1 - |2t - 1|, \quad t \in \mathbb{R}.\]
We refer to [17]–[20] for the general theory of numerical integration using lattice rules, and to [15], [21], [22] for the use of tent-transformed lattice rules in the context of non-periodic functions. Our motivation for using tent-transformed lattice rules stems from their ability to achieve second-order convergence in Sobolev spaces of smoothness \(\alpha = 2\); see [15].
3 Theoretical Considerations↩︎
Classical QMC theory typically bounds the approximation error by the inequality \[|I f - Q_{{\boldsymbol{T}}}^{{\boldsymbol{w}}} f| \leq e(Q_{{\boldsymbol{T}}}^{{\boldsymbol{w}}}, {\mathcal{H}})\, \| f \|_{{\mathcal{H}}}, \qquad f \in {\mathcal{H}}.\] Since the second factor \(\| f \|_{{\mathcal{H}}}\) is constant with respect to \(n\) for fixed \(f\) (often normalized to one for simplicity), the error is ultimately controlled by the worst-case error \(e(Q_{{\boldsymbol{T}}}^{{\boldsymbol{w}}}, {\mathcal{H}})\). While \({\boldsymbol{w}}^{\ast}\) minimizes the worst-case error among all possible weights, \[e(Q_{{\boldsymbol{T}}}^{{\boldsymbol{w}}^{\ast}}, {\mathcal{H}}) \leq e(Q_{{\boldsymbol{T}}}^{{\boldsymbol{w}}}, {\mathcal{H}}) \qquad \text{for all } {\boldsymbol{w}}\in {\mathbb{R}}^{n},\] its construction via an orthogonal projection further improves this bound by replacing the constant factor \(\| f \|_{{\mathcal{H}}}\) with the distance \(\mathrm{dist}_{{\mathcal{H}}}(f,V_{{\boldsymbol{T}}})\), which can be expected to decrease as \(n\) grows and \(V_{{\boldsymbol{T}}}\) increasingly approximates \({\mathcal{H}}\):
Let \({\boldsymbol{T}}= ({\boldsymbol{t}}_k)_{k=0}^{n-1} \in D^n\) be any point set in \(D\) and \({\mathcal{H}}\) be an RKHS with kernel \(K\colon D\times D \to \mathbb{R}\). Let \(V_{{\boldsymbol{T}}} = \mathrm{span}(K({\boldsymbol{t}}_0,\setbox 0x to\wd 0{\hss\cdot\hss}), \ldots, K({\boldsymbol{t}}_{n-1},\setbox 0x to\wd 0{\hss\cdot\hss})) \subseteq \mathcal{H}\) and let \({\boldsymbol{w}}^{\ast}\) satisfy ?? . Then \[|If - Q_{{\boldsymbol{T}}}^{{\boldsymbol{w}}^{\ast}} f| \leq e(Q_{{\boldsymbol{T}}}^{{\boldsymbol{w}}^{\ast}}, {\mathcal{H}})\, \mathrm{dist}_{{\mathcal{H}}}(f,V_{{\boldsymbol{T}}}),\] where \(\mathrm{dist}_{{\mathcal{H}}}(f,V_{{\boldsymbol{T}}})\) denotes the distance between \(f\) and \(V_{{\boldsymbol{T}}}\) in \({\mathcal{H}}\).
Proof. Since \(h_{{\boldsymbol{T}}}^{{\boldsymbol{w}}^{\ast}} = P_{V_{{\boldsymbol{T}}}} h\), we obtain \(h-h_{{\boldsymbol{T}}}^{{\boldsymbol{w}}^{\ast}} \perp V_{{\boldsymbol{T}}}\) and, by the Cauchy–Schwarz inequality \[\begin{align} |If - Q_{{\boldsymbol{T}}}^{{\boldsymbol{w}}^{\ast}} f| &= |\langle h-h_{{\boldsymbol{T}}}^{{\boldsymbol{w}}^{\ast}} , f \rangle_{{\mathcal{H}}}| \\ &= |\langle h-h_{{\boldsymbol{T}}}^{{\boldsymbol{w}}^{\ast}} , f - P_{V_{{\boldsymbol{T}}}} f \rangle_{{\mathcal{H}}}| \\ &\leq \| h-h_{{\boldsymbol{T}}}^{{\boldsymbol{w}}^{\ast}} \|_{{\mathcal{H}}} \, \| f - P_{V_{{\boldsymbol{T}}}} f \|_{{\mathcal{H}}} \\ &= e(Q_{{\boldsymbol{T}}}^{{\boldsymbol{w}}^{\ast}}, {\mathcal{H}})\, \mathrm{dist}_{{\mathcal{H}}}(f,V_{{\boldsymbol{T}}}), \end{align}\] proving the claim.0◻ ◻
We noted above that the distance \(\mathrm{dist}_{{\mathcal{H}}}(f, V_{{\boldsymbol{T}}})\) can be expected to decrease as \(n\) increases. This is intuitive, as the space \(V_{{\boldsymbol{T}}}\) becomes richer and more capable of approximating elements in \({\mathcal{H}}\). However, to the best of our knowledge, no theoretical result is currently available that quantifies the rate of this convergence across arbitrary RKHSs. Quantitative rates are known in specific cases, for example for periodic Sobolev spaces and certain shift-invariant kernels, but do not appear to extend directly to the non-periodic Sobolev spaces of dominating mixed smoothness considered in this paper. There appear to be connections to Gaussian process regression on lattices, which is an active area of research [23], [24].
3.1 Rate Doubling in the One-Dimensional Setting↩︎
In this subsection, we investigate the effect of using optimal weights for an (unshifted) lattice rule in dimension \(s=1\), which simply corresponds to a left-Riemann rule, that is, \[t_k=\frac{k}{n},\quad k=0,\ldots,n-1.\] If \(f\in \mathcal{H}_{1,\mathbf{1}}^1\), then it is a consequence of 7 and the identity \(\sum_{k,k'=0}^{n-1}K_{1,\mathbf{1}}^1(\frac{k}{n},\frac{k'}{n})=\frac{3n^2+1}{3}\) that the equally weighted quadrature rule \(Q_{\boldsymbol{T}} f =\frac{1}{n}\sum_{k=0}^{n-1}f(t_k)\) admits the error rate \[|I f - Q_{\boldsymbol{T}} f|=\mathcal{O}(n^{-1}).\]
Defining the sequence of weights \((w_k^*)_{k=0}^{n-1}\) as the solution to the system ?? with \(K_{1,\mathbf{1}}^1(x,y)=1+\frac{1}{2}B_2(|x-y|)+(x-1/2)(y-1/2)\), \(x,y\in[0,1]\), denoting the one-dimensional kernel corresponding to \(\mathcal{H}_{1,\mathbf{1}}^1\), we obtain \[w_0^*=\frac{1}{2n}\frac{12n^3}{12n^3+n+3}, \qquad w_k^*=2w_0^*,\quad k\in\{1,\ldots,n-2\}, \qquad w_{n-1}^*=3w_0^*.\] In this special case, the optimally weighted quadrature rule \(Q_{\boldsymbol{T}}^{{\boldsymbol{w}}^{\ast}} f = \sum_{k=0}^{n-1} w_{k}^{\ast} f(t_k)\) exhibits a quadratic error rate:
Lemma 2. Suppose that \(f\in \mathcal{H}_{1,\mathbf{1}}^2\). Then \[|I f - Q_{\boldsymbol{T}}^{\boldsymbol{w}^*} f |=\mathcal{O}(n^{-2}).\]
Proof. The quadrature error can be recast as \[\begin{align} \int_0^1 f(y)\,{\rm d}y-\sum_{k=0}^{n-1}w_k^*f(t_k)&=\int_0^1\langle f,K_{1,\mathbf{1}}^1(\setbox 0x to\wd 0{\hss\cdot\hss},y)\rangle_{\mathcal{H}_{1,\mathbf{1}}^1}\,{\rm d}y-\sum_{k=0}^{n-1}w_k^* \langle f,K_{1,\mathbf{1}}^1(\setbox 0x to\wd 0{\hss\cdot\hss},t_k)\rangle_{\mathcal{H}_{1,\mathbf{1}}^1}\\ &=\bigg\langle f,\int_0^1 K_{1,\mathbf{1}}^1(\setbox 0x to\wd 0{\hss\cdot\hss},y)\,{\rm d}y-\sum_{k=0}^{n-1}w_k^*K_{1,\mathbf{1}}^1(\setbox 0x to\wd 0{\hss\cdot\hss},t_k)\bigg\rangle_{\mathcal{H}_{1,\mathbf{1}}^1}\\ &=\langle f,h-h_{{\boldsymbol{T}}}^{{\boldsymbol{w}}^{\ast}}\rangle_{\mathcal{H}_{1,\mathbf{1}}^1}, \end{align}\] where \(h:=\int_0^1 K_{1,\mathbf{1}}^1(\setbox 0x to\wd 0{\hss\cdot\hss},y)\,{\rm d}y=\mathbb{1}\) and \(h_{{\boldsymbol{T}}}^{{\boldsymbol{w}}^{\ast}}:=\sum_{k=0}^{n-1}w_k^*K_{1,\mathbf{1}}^1(\setbox 0x to\wd 0{\hss\cdot\hss},t_k)\) is the kernel interpolant of \(h\). Using the definition of the inner product in the weighted Sobolev space of smoothness \(\alpha=1\) we have \[\begin{align} \langle f,h-h_{{\boldsymbol{T}}}^{{\boldsymbol{w}}^{\ast}}\rangle_{\mathcal{H}_{1,\mathbf{1}}^1}&=\bigg(\int_0^1f(y)\,{\rm d}y\bigg)\bigg(\int_0^1 (h(y)-h_{{\boldsymbol{T}}}^{{\boldsymbol{w}}^{\ast}}(y))\,{\rm d}y\bigg)\\ &\quad +\int_0^1f'(y)(h'(y)-(h_{{\boldsymbol{T}}}^{{\boldsymbol{w}}^{\ast}})'(y))\,{\rm d}y. \end{align}\] Using integration by parts, the absolute value of the latter integral can be estimated as \[\begin{align} & \bigg| \int_0^1f'(y)(h'(y)-(h_{{\boldsymbol{T}}}^{{\boldsymbol{w}}^{\ast}})'(y))\,{\rm d}y \bigg|\\ &= \bigg| f'(1)(h(1)-h_{{\boldsymbol{T}}}^{{\boldsymbol{w}}^{\ast}}(1))-f'(0)(h(0)-h_{{\boldsymbol{T}}}^{{\boldsymbol{w}}^{\ast}}(0))-\int_0^1 f''(y)(h(y)-h_{{\boldsymbol{T}}}^{{\boldsymbol{w}}^{\ast}}(y))\,{\rm d}y \bigg|\\ &\leq \frac{6n}{12n^3+n+3} |f'(1)| +\|f''\|_{L^2(0,1)}\|h-h_{{\boldsymbol{T}}}^{{\boldsymbol{w}}^{\ast}}\|_{L^2(0,1)}, \end{align}\] where we used the fact that \(h(0)-h_{\boldsymbol{T}}^{\boldsymbol{w}^*}(0)=0\), a consequence of the fact that \(h_{\boldsymbol{T}}^{\boldsymbol{w}^*}\) interpolates \(h\) at the lattice point \(t_0=0\), as well as the identity \[\begin{align} h(1)-h_{\boldsymbol{T}}^{\boldsymbol{w}^*}(1)&=1-\frac{12n^2}{12n^3+n+3}\bigg(\frac{1}{2} K_{1,\mathbf{1}}^1(1,0)+\frac{3}{2} K_{1,\mathbf{1}}^1\bigg(1,\frac{n-1}{n}\bigg)+\sum_{k=1}^{n-2}K_{1,\mathbf{1}}^1\bigg(1,\frac{k}{n}\bigg)\bigg)\\ &=1-\frac{12n^2}{12n^3+n+3}\bigg(\frac{29}{12} + \frac{3}{4n^2}-\frac{3}{2n}+\sum_{k=1}^{n-2}\frac{5n^2+3k^2}{6n^2}\bigg)\\ &=1-\frac{12n^2}{12n^3+n+3}\bigg(n-\frac{5}{12n}+\frac{1}{4n^2}\bigg)\\ &=\frac{6n}{12n^3+n+3}, \end{align}\] where we used \(K_{1,\mathbf{1}}^1\big(1,\frac{k}{n}\big)=\frac{5}{6}+\frac{k^2}{2n^2}\), \(k\in\{0,\ldots,n-1\}\), and \(\sum_{k=1}^{n-2}k^2=\frac{(2n-3)(n-1)(n-2)}{6}\). Therefore \[|\langle f,h-h_{{\boldsymbol{T}}}^{{\boldsymbol{w}}^{\ast}}\rangle_{\mathcal{H}_{1,\mathbf{1}}^1}|\leq \frac{6n}{12n^3+n+3}|f'(1)| + \|h-h_{{\boldsymbol{T}}}^{{\boldsymbol{w}}^{\ast}}\|_{L^2(0,1)}\big(\|f\|_{L^2(0,1)}+\|f''\|_{L^2(0,1)}\big).\] Since \(\frac{6n}{12n^3+n+3}|f'(1)|=\mathcal{O}(n^{-2})\), it remains to assess the convergence rate of \(\|h-h_{{\boldsymbol{T}}}^{{\boldsymbol{w}}^{\ast}}\|_{L^2(0,1)}\). Making use of the identities \[\begin{align} &\sum_{k=0}^{n-1}w_k^*=\frac{24n^2}{12n^3+n+3}+\frac{12n^2(n-2)}{12n^3+n+3},\\ &\sum_{k=0}^{n-1}\sum_{\ell=0}^{n-1}w_k^*w_\ell^* \bigg(\frac{46}{45}-\frac{1}{6} t_k^2+\frac{1}{12}t_k^3-\frac{1}{24}t_k^4+\frac{1}{4} t_k^2t_\ell -\frac{1}{6} t_\ell^2\\ &\quad\quad\quad+\frac{1}{4} t_kt_\ell^2-\frac{1}{4} t_k^2t_\ell^2+\frac{1}{12}t_\ell^3-\frac{1}{24}t_\ell^4+\frac{1}{12}|t_k-t_\ell|^3\bigg)=\frac{720n^6+n^2+60n-45}{5(12n^3+n+3)^2}, \end{align}\] where the latter identity is valid for \(n\geq 2\), we obtain, for \(n\geq 2\), \[\begin{align} \|h-&h_{{\boldsymbol{T}}}^{{\boldsymbol{w}}^{\ast}}\|_{L^2(0,1)}^2 = \int_0^1 \bigg(1-\sum_{k=0}^{n-1}w_k^*K_{1,\mathbf{1}}^1(x,t_k)\bigg)^2\,{\rm d}x\\ &=\int_0^1\bigg(1-2\sum_{k=0}^{n-1}w_k^*K_{1,\mathbf{1}}^1(x,t_k)+\sum_{k=0}^{n-1}\sum_{\ell=0}^{n-1}w_k^*w_\ell^* K_{1,\mathbf{1}}^1(x,t_k)K_{1,\mathbf{1}}^1(x,t_\ell)\bigg)\,{\rm d}x\\ &=1-2\sum_{k=0}^{n-1}w_k^*+\sum_{k=0}^{n-1}\sum_{\ell=0}^{n-1}w_k^*w_\ell^* \bigg(\frac{46}{45}-\frac{1}{6} t_k^2+\frac{1}{12}t_k^3-\frac{1}{24}t_k^4+\frac{1}{4} t_k^2t_\ell -\frac{1}{6} t_\ell^2\\ &\quad\quad\!\!+\frac{1}{4} t_kt_\ell^2-\frac{1}{4} t_k^2t_\ell^2+\frac{1}{12}t_\ell^3-\frac{1}{24}t_\ell^4+\frac{1}{12}|t_k-t_\ell|^3\bigg) \\ &=\frac{6n(n+15)}{5(12n^3+n+3)^2}. \end{align}\] Hence, \(|If-Q_{{\boldsymbol{T}}}^{\boldsymbol{w}^*}f|=\mathcal{O}(n^{-2})\) as claimed.0◻ ◻
Lemma 2 illustrates that using kernel cubature can improve the error convergence rate of a lattice rule in one dimension.
4 Numerical Experiments↩︎
We briefly describe the construction of the Gramian matrix \(\mathcal{K}_{\boldsymbol{T}}\) for common types of coordinate weights and analyze the computational complexity in Section 4.1. We then investigate the behavior of kernel cubature applied to QMC point sets in two experiments. In Section 4.2, we compare the performance of equally weighted lattice rules against kernel cubature constructed for the same point sets, applied to an elliptic partial differential equation (PDE) with a parametric input coefficient, while Section 4.3 compares the worst-case errors of equally weighted and optimally weighted QMC point sets in unweighted Sobolev spaces \(\mathcal{H}_{s,\mathbf{1}}^{\alpha}\) with varying smoothness parameters \(\alpha\).
4.1 Construction of the Gramian Matrix and Its Computational Complexity↩︎
In the context of information-based complexity and uncertainty quantification, it is often the case to consider the dominant cost to be the \(n\) evaluations of the integrand \(f\), see e.g., [25], [26]. By contrast, the construction of the lattice rule and the computation of the associated weights are typically negligible in comparison and can often be performed “offline”.
That said, for completeness, we provide here a brief analysis of the computational cost associated with constructing the cubature weights, which involves two main steps: (1) assembling the Gramian matrix \({\mathcal{K}}_{{\boldsymbol{T}}}\), and (2) solving the resulting linear system. The second step, solving the resulting linear system, has a computational complexity of \(\mathcal{O}(n^3)\) using standard direct methods such as Cholesky decomposition, which is appropriate here given that the kernel matrix is symmetric and positive definite. We now turn our attention to the first step: computing entries of the kernel matrix, i.e., evaluating \(K_{s,{\boldsymbol{\gamma}}}^{\alpha}({\boldsymbol{x}}, {\boldsymbol{y}})\) for given \({\boldsymbol{x}}, {\boldsymbol{y}}\in [0,1]^s\). Recall that \[K_{s,{\boldsymbol{\gamma}}}^{\alpha}({\boldsymbol{x}},{\boldsymbol{y}})=\sum_{{\mathfrak{u}}\subseteq\{1,\ldots,s\}}\gamma_{{\mathfrak{u}}}\prod_{j\in{\mathfrak{u}}}\eta_{\alpha}(x_j,y_j),\label{KaKlScSu:eq:generickernel}\tag{8}\] where \(\eta_{\alpha}(x,y)=\sum_{\tau=1}^{\alpha}\frac{1}{(\tau!)^2}B_{\tau}(x)B_{\tau}(y)+\frac{(-1)^{\alpha+1}}{(2\alpha)!}\widetilde{B}_{2\alpha}(x-y)\). The cost of evaluating this kernel depends on the structure of the coordinate weights \({\boldsymbol{\gamma}}=(\gamma_{{\mathfrak{u}}})_{{\mathfrak{u}}\subseteq\{1,\ldots,s\}}\). We refer to [9] for details and provide a brief overview below. In what follows, it is assumed that evaluating \(\eta_\alpha\) has a constant cost. We also use the convention that a product over an empty set is defined to be equal to 1.
Product weights \(\gamma_{{\mathfrak{u}}}=\prod_{j\in{\mathfrak{u}}}\widetilde{\gamma}_j\) are specified by a sequence of nonnegative numbers \((\widetilde{\gamma}_j)_{j=1}^s\), and the kernel can be equivalently written as \[K_{s,{\boldsymbol{\gamma}}}^{\alpha}({\boldsymbol{x}},{\boldsymbol{y}})=\prod_{j=1}^s (1+\widetilde{\gamma}_j\eta_\alpha(x_j,y_j)).\] This expression can be evaluated in \(\mathcal{O}(s)\) time for one pair \(({\boldsymbol{x}},{\boldsymbol{y}})\) and, for any point set \(\boldsymbol{T}=(\boldsymbol{t}_k)_{k=0}^{n-1}\) in \([0,1]^s\), the matrix \(\mathcal{K}_{\boldsymbol{T}}=(K_{s,\boldsymbol{\gamma}}^{\alpha}({\boldsymbol{t}}_k,{\boldsymbol{t}}_\ell))_{k,\ell=0}^{n-1}\) can be assembled in \(\mathcal{O}(sn^2)\) time.
Product-and-order dependent (POD) weights \(\gamma_{{\mathfrak{u}}}=\Gamma_{|{\mathfrak{u}}|}\prod_{j\in{\mathfrak{u}}}\widetilde{\gamma}_j\) are specified by two sequences of nonnegative numbers \((\Gamma_k)_{k=0}^s\) and \((\widetilde{\gamma}_j)_{j=1}^s\), and the kernel can be equivalently written as \[\begin{align} K_{s,{\boldsymbol{\gamma}}}^{\alpha}({\boldsymbol{x}},{\boldsymbol{y}})=\sum_{\ell=0}^s \Gamma_{\ell}P_{s,\ell},\label{KaKlScSu:eq:podkern} \end{align}\tag{9}\] where the sequence \((P_{k,\ell})_{k,\ell=0}^s\) can be computed recursively by \[\begin{align} &P_{k,0}=1\quad\text{for all}~k\in\{0,\ldots,s\},\\ &P_{k,\ell}=0\quad\text{for all}~k\in\{0,\ldots,s\}~\text{and}~\ell\in\{k+1,\ldots,s\},\\ &P_{k,\ell}=P_{k-1,\ell}+\widetilde{\gamma}_k \eta_\alpha(x_k,y_k)P_{k-1,\ell-1}~\text{for all}~k\in\{1,\ldots,s\}~\text{and}~\ell\in\{1,\ldots,k\}. \end{align}\] The cost to obtain \(K_{s,{\boldsymbol{\gamma}}}^\alpha({\boldsymbol{x}},{\boldsymbol{y}})\) using the expression 9 is \(\mathcal{O}(s^2)\) for one pair \(({\boldsymbol{x}},{\boldsymbol{y}})\) and, for any point set \(\boldsymbol{T}=(\boldsymbol{t}_k)_{k=0}^{n-1}\) in \([0,1]^s\), the matrix \(\mathcal{K}_{\boldsymbol{T}}=(K_{s,\boldsymbol{\gamma}}^{\alpha}({\boldsymbol{t}}_k,{\boldsymbol{t}}_\ell))_{k,\ell=0}^{n-1}\) can be assembled in \(\mathcal{O}(s^2n^2)\) time.
Since the system matrix \(\mathcal{K}_{\boldsymbol{T}}\) is symmetric, it has at most \(\frac{n(n+1)}{2}\) unique entries that need to be constructed. The cost of solving the matrix equation ?? is \(\mathcal{O}(n^3)\) independently of the weights \({\boldsymbol{\gamma}}\) and dimension \(s\).
It is possible to obtain recurrence formulas for other classes of coordinate weights: for example, using smoothness-driven product-and-order dependent (SPOD) weights with smoothness degree \(\sigma\in\mathbb{N}\), the cost to obtain \(K_{s,\boldsymbol{\gamma}}^{\alpha}(\boldsymbol{x},\boldsymbol{y})\) is \(\mathcal{O}(s^2\sigma^2)\) for one pair \((\boldsymbol{x},\boldsymbol{y})\) and, for any point set \(\boldsymbol{T}=(\boldsymbol{t}_k)_{k=0}^{n-1}\) in \([0,1]^s\), the matrix \(\mathcal{K}_{\boldsymbol{T}}=(K_{s,\boldsymbol{\gamma}}^{\alpha}({\boldsymbol{t}}_k,{\boldsymbol{t}}_\ell))_{k,\ell=0}^{n-1}\) can be assembled in \(\mathcal{O}(s^2\sigma^2n^2)\) time. For details, we refer to [9].
4.2 Application to PDE Uncertainty Quantification↩︎
Let \(\Omega=(0,1)^2\). We consider the elliptic PDE \[\begin{align} \begin{cases} -\nabla \setbox 0x to\wd 0{\hss\cdot\hss}(a(\boldsymbol{x},\boldsymbol{y})\nabla u(\boldsymbol{x},\boldsymbol{y}))=f(\boldsymbol{x}),&\boldsymbol{x}\in \Omega,~\boldsymbol{y}\in[-\tfrac12,\tfrac12]^s,\\ u(\boldsymbol{x},\boldsymbol{y})=0,&\boldsymbol{x}\in\partial \Omega,~\boldsymbol{y}\in[-\tfrac12,\tfrac12]^s, \end{cases}\label{KaKlScSu:eq:pde} \end{align}\tag{10}\] equipped with the parametric diffusion coefficient \[\begin{align} \label{KaKlScSu:eq:pde2} a(\boldsymbol{x},\boldsymbol{y})=\frac{1}{2}+\frac{1}{2}\sum_{j=1}^s j^{-2}y_j\sin(j\pi x_{1})\sin(j\pi x_{2}) \end{align}\tag{11}\] for \(\boldsymbol{x}=(x_1,x_2)\in \Omega\) and \(\boldsymbol{y}=(y_1,\ldots,y_s)\in[-\tfrac12,\tfrac12]^s\), where \(s\) is referred to as the truncation dimension. In this case, it can be shown (cf.[27]) that, using the POD coordinate weights \[\begin{align} \label{KaKlScSu:eq:podweights} \gamma_{{\mathfrak{u}}}=\bigg(|{\mathfrak{u}}|!\prod_{j\in{\mathfrak{u}}}\frac{b_j}{\sqrt{{2\zeta(2\lambda)}/{(2\pi^2)^\lambda}}}\bigg)^{\frac{2}{1+\lambda}}\quad \text{for all}~{\mathfrak{u}}\subseteq\{1,\ldots,s\}, \end{align}\tag{12}\] with \(b_j=(1-\frac{1}{2}\zeta(2))^{-1} j^{-2}\) and \(\lambda=\frac{1}{2-2\delta}\), \(\delta=0.05\), randomly shifted rank-1 lattice rules achieve dimension-independent convergence rates for the root-mean-square integration error of the expected value \[\begin{align} \mathbb{E}[G(u)]=\int_{[-\frac{1}{2},\frac{1}{2}]^s}G(u(\setbox 0x to\wd 0{\hss\cdot\hss},\boldsymbol{y}))\,{\mathrm d}\boldsymbol{y},\label{KaKlScSu:eq:expectedvalue} \end{align}\tag{13}\] where \(G\!:H_0^1(\Omega)\to\mathbb{R}\) is an arbitrary bounded linear functional (the quantity of interest). The PDE 10 was discretized using a first-order finite element method with mesh width \(h=2^{-5}\).
Figure 1: The cubature errors for the PDE example 10 –11 were computed using both equally weighted rank-1 lattice rules (“QMC PDE”) and kernel cubature rules (“Kernel PDE”), with both methods evaluated on the same lattice point sets. We used \(s\in\{1,5,20,100\}\) as the truncation dimensions. The weights for the kernel cubature were obtained subject to \(\mathcal{H}_{s,{\boldsymbol{\gamma}}}^1\) and we also illustrate the computed worst-case errors for both equally weighted lattice rules (“QMC WCE”) and the corresponding kernel cubature rules (“Kernel WCE”) in the space \(\mathcal{H}_{s,{\boldsymbol{\gamma}}}^1\). The cubature errors have been averaged over \(R=8\) random shifts..
We set \(s\in\{1,5,20,100\}\) as the truncation dimension in 11 , choose \(f(\boldsymbol{x})=x_1\) as the source term, and set \(G(v)=\int_{\Omega}v(\boldsymbol{x})\,{\mathrm d}\boldsymbol{x}\) as the quantity of interest. We used the fast component-by-component algorithm [28] to find an extensible generating vector for \(n=2^k\), \(k=1,\ldots,12\), corresponding to the POD coordinate weights 12 . We approximated the expected value 13 using two methods based on the same fixed lattice rule:
Equally weighted lattice QMC rules for \(n=2^k\), \(k=1,\ldots,10\);
Kernel cubature rules over the same lattice points for \(n=2^k\), \(k=1,\ldots,10\).
For error estimation, we applied \(R=8\) random shifts to each point set. As the reference solution, we used the numerical approximations corresponding to an \(n=2^{12}\) point QMC rule for experiment 1 and an \(n=2^{12}\) point kernel cubature point set for experiment 2. To obtain the weights \(\boldsymbol{w}^*\) of the randomly shifted kernel cubature rules, we solved the linear system ?? for each randomly shifted lattice point set corresponding to the kernel ?? with \(\alpha=1\) and POD coordinate weights 12 . In this case, it is necessary to assemble the elements of the Gram matrix \(\mathcal{K}_{{\boldsymbol{T}}}\) recursively; see Section 4.1 for details. In addition, we computed the worst-case errors for both methods using the square root of the expression 7 . The results are displayed in Figure 1.
While the kernel cubature does not improve the essentially linear worst-case error rate of the underlying rank-1 lattice point set, the numerical results seem to indicate that the kernel cubature rate is significantly better. For dimension \(s=1\), the kernel cubature rate for the PDE problem is double that of the equally weighted cubature rule (as indicated by Lemma 2) while for increased dimensions \(s\) the observed cubature convergence rate lies approximately between \(-1.6\) and \(-1.7\). This improvement may be attributed to the error decomposition presented in Proposition [KaKlScSu:prop:approximation95error95additional95gain]: if the distance between \(f\) and \(V_{\boldsymbol{T}}\) is decreasing as \(n\) grows, then the observed kernel cubature rate can exceed the worst-case error rate by a significant amount.
4.3 Worst-Case Errors for Sobolev Spaces of Varying Smoothness↩︎
In this subsection, we demonstrate improved convergence rates of lattice rules with optimized weights by directly calculating worst-case errors. We first optimize the weight for \({\mathcal{H}}^2_{s,{\boldsymbol{1}}}\) and calculate the worst-case error. We consider shifted and tent-transformed lattice rules \({\boldsymbol{T}}= \phi(\{\tilde{{\boldsymbol{T}}} + \Delta\})\) as described in Section 2.4, where a single random shift \(\Delta\) is added to the lattice before the transformation, preventing points from coinciding after the transform is applied (which would render the Gram matrix \({\mathcal{K}}_{{\boldsymbol{T}}}\) singular). It is shown in [15] that second order convergence can be achieved by a tent-transformed lattice rule in \({\mathcal{H}}^{2}_{s,\boldsymbol{\gamma}}\). In this numerical example we show that the optimally weighted tent-transformed lattice rule can achieve even faster convergence in \({\mathcal{H}}^{4}_{2,{\boldsymbol{1}}}\).
Figure 2 shows the two-dimensional case. We use the generating vector \({\boldsymbol{z}}=(1,182667)\) for \(n=2^2,\ldots,2^{10}\), from [29] (available online at Frances Kuo’s webpage [30]). We remind the reader that the worst-case errors can be easily calculated by 6 and 7 , i.e., \[\begin{align} e(Q_{{\boldsymbol{T}}}^{{\boldsymbol{w}}^*} , {\mathcal{H}}^{\alpha}_{s,\boldsymbol{\gamma}})^{2} &= 1-\sum_{k=0}^{n-1} w_k^{\ast}, \\ e(Q_{{\boldsymbol{T}}} , {\mathcal{H}}^{\alpha}_{s,\boldsymbol{\gamma}})^{2} &= -1+n^{-2}\kern-2mm\sum_{k,k'=0} ^{n-1} K^{\alpha}_{s,\boldsymbol{\gamma}} ({\boldsymbol{t}}_k,{\boldsymbol{t}}_{k'}), \end{align}\] and \[e(Q_{{\boldsymbol{T}}}^{{\boldsymbol{w}}^*} , {\mathcal{H}}^{2\alpha}_{s,\boldsymbol{\gamma}})^{2} = 1-2\sum_{k=0}^{n-1} w_k^{\ast}+\sum_{k=0}^{n-1}\sum_{k'=0}^{n-1} w^{\ast}_{k}w^{\ast}_{k'} K^{2\alpha}_{s,\boldsymbol{\gamma}} ({\boldsymbol{t}}_k,{\boldsymbol{t}}_{k'}),\] where in the last line, the weights \({\boldsymbol{w}}^{\ast}\) are optimized for \({\mathcal{H}}^{\alpha}_{s,\boldsymbol{\gamma}}\), not \({\mathcal{H}}^{2\alpha}_{s,\boldsymbol{\gamma}}\).
We observe that a convergence rate faster than second order is attained. Note that there are other ways to achieve convergence faster than \(\mathcal{O}(n^{-2})\) using lattice rules for non-periodic functions. One example is the periodization strategy, which uses a change of variables in order to obtain a periodic integrand from a non-periodic one. However, it is not known how to avoid the curse of dimensionality when using this strategy, see [31]. Another example is the symmetrized lattice rule [15], but this is also cursed by dimensionality because the required number of points grows exponentially with the dimension.
5 Concluding Remarks↩︎
In this paper, we introduced a weighted version of QMC cubature, referred to as kernel cubature, where the weights are chosen to minimize the distance between the kernel mean embeddings of the true and empirical distributions. We provided a theoretical result (Proposition [KaKlScSu:prop:approximation95error95additional95gain]) suggesting an improved convergence rate for kernel cubature compared to the equally weighted case, proved a corresponding statement for dimension \(s=1\), and presented numerical results in dimensions \(s = 1, 5, 20,\) and \(100\), focusing on lattice rules and an elliptic PDE problem with a random coefficient. Additionally, we explored the behavior of the worst-case error when weights optimized for a Sobolev space of dominating-mixed smoothness \(\alpha\) were applied to a space of higher smoothness. In dimension \(s=2\), this led to a significant acceleration in the convergence rate.
While this work established several theoretical insights, its primary focus was experimental. Future research will aim to establish improved convergence guarantees in arbitrary dimensions, as well as investigate constructions of lattices and other QMC point sets specifically designed for kernel cubature.
Acknowledgments↩︎
IK and CS were funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany’s Excellence Strategy (EXC-2046/1, project 390685689) through project EF1-19 of the Berlin Mathematics Research Center MATH+. This work of YS was supported by the Research Council of Finland (decisions 348503 and 359181). The work of VK was supported by the Research Council of Finland (Flagship of Advanced Mathematics for Sensing, Imaging and Modelling grant 359183). We thank Frances Kuo, Fred Hickernell and Robert Gruhlke for helpful collegial discussions.