MDFEM: Multivariate decomposition finite element method for elliptic PDEs with lognormal diffusion coefficients using higher-order QMC and FEM

TL;DR

This paper presents the MDFEM, a novel method combining multivariate decomposition, higher-order QMC, and FEM to efficiently approximate expectations of solutions to elliptic PDEs with lognormal coefficients, achieving higher-order convergence.

Contribution

The paper introduces the MDFEM, integrating multivariate decomposition with higher-order QMC and FEM, for improved efficiency in high-dimensional elliptic PDEs with lognormal coefficients.

Findings

Achieves higher-order convergence rates in error versus cost.

Develops higher-order QMC rules for Gaussian integrals.

Provides complexity estimates for the MDFEM algorithm.

Abstract

We introduce the multivariate decomposition finite element method for elliptic PDEs with lognormal diffusion coefficient $a=\exp(Z)$ where $Z$ is a Gaussian random field defined by an infinite series expansion $Z(\boldsymbol{y}) = \sum_{j\ge1} y_j\,\phi_j$ with $y_j\sim\mathcal{N}(0,1)$ and a given sequence of functions $\{\phi_j\}_{j\ge1}$. We use the MDFEM to approximate the expected value of a linear functional of the solution of the PDE which is an infinite-dimensional integral over the parameter space. The proposed algorithm uses the multivariate decomposition method (MDM) to compute the infinite-dimensional integral by a decomposition into finite-dimensional integrals, which we resolve using quasi-Monte Carlo (QMC) methods, and for which we use the finite element method (FEM) to solve different instances of the PDE. We develop higher-order quasi-Monte Carlo rules for integration over the finite-dimensional Euclidean space with respect to the Gaussian distribution by use of a truncation strategy. By linear transformations of interlaced polynomial lattice rules from the unit cube to a multivariate box of the Euclidean space we achieve higher-order convergence rates for functions belonging to a class of anchored Gaussian Sobolev spaces, taking into account the truncation error. Under appropriate conditions, the MDFEM achieves higher-order convergence rates in term of error versus cost, i.e., to achieve an accuracy of $O(\epsilon)$ the computational cost is $O(\epsilon^{-1/\lambda-d'/\lambda}) = O(\epsilon^{-(p^*+d'/\tau)/(1-p^*)})$ where $\epsilon^{-1/\lambda}$ and $\epsilon^{-d'/\lambda}$ are respectively the cost of the quasi-Monte Carlo cubature and the finite element approximations, with $d' = d \, (1+\delta')$ for some $\delta' \ge 0$ and $d$ the physical dimension, and $0 < p^* \le (2+d'/\tau)^{-1}$ is a parameter representing the sparsity of $\{\phi_j\}_{j\ge1}$.

Equations

$$-\nabla \cdot (a(\mathbf{x},\mathbf{y})\nabla u(\mathbf{x},\mathbf{y}))$$

$$u(\mathbf{x},\mathbf{y})$$

$$a(\mathbf{x},\mathbf{y})$$

$$Z(\mathbf{x},\mathbf{y})$$

$$\mathbb{E}[G(u)]=\mathcal{I}(G(u))$$

$$:=\underset{s\to \infty }{\lim }{\int }_{{\mathbb{R}}^s}G(u(\cdot ,y_1,y_2,\dots ,y_s,0,0,\dots ))\mathrm{d}\mu ({\mathbf{y}}_{\{1:s\}}),$$

$$\mathrm{d}\mu ({\mathbf{y}}_{\mathfrak{u}})$$

$${\int }_{D}a(\mathbf{x},\mathbf{y})\nabla u(\mathbf{x},\mathbf{y})\cdot \nabla v(\mathbf{x})\mathrm{d}\mathbf{x}$$

$$\Vert v{\Vert }_{L_{p,\rho }({\mathbb{R}}^{\mathbb{N}};X)}$$

$$\kappa$$

$$\{b_j{\}}_{j\ge 1}\in {\ell }^{p^{\ast }}(\mathbb{N})\text{ for some }{p}^{\ast }\in (0,\infty ),$$

$$a_{\max }(\mathbf{y})$$

$$a_{\min }(\mathbf{y})$$

$$\Vert u(\cdot ,\mathbf{y}){\Vert }_V$$

$$\Vert u{\Vert }_{L_{p,\rho }({\mathbb{R}}^{\mathbb{N}};V)}\le {\Vert \frac{1}{a_{\min }}\Vert }_{L_{p,\rho }({\mathbb{R}}^{\mathbb{N}})}\Vert f{\Vert }_{V^{\ast }}<\infty ,$$

$$\Vert Z(\cdot ,\mathbf{y}){\Vert }_{L^{\infty }(D)}$$

$$\le \underset{\mathbf{x}\in D}{\sup }\sum _{j\ge 1}|y_j|\frac{b_j}{b_j}|{\varphi }_j(\mathbf{x})|\le \left(\underset{j\ge 1}{\sup }|y_j|b_j\right)\left(\underset{\mathbf{x}\in D}{\sup }\sum _{j\ge 1}\frac{|{\varphi }_j(\mathbf{x})|}{b_j}\right)<\infty ,$$

$$a_s(\mathbf{x},\mathbf{y})$$

$$\mathcal{I}(\mathcal{F})={\int }_{{\mathbb{R}}^{\mathbb{N}}}\mathcal{F}(\mathbf{y})\mathrm{d}\mu (\mathbf{y})$$

$$\mathcal{F}(y_1,\dots ,y_s,0,0,\dots )$$

$$\mathcal{F}(y_1,y_2,\dots )$$

$${\mathcal{F}}_{{\mathfrak{U}}_{ϵ}}(\mathbf{y})$$

$$\left|\mathcal{I}(\mathcal{F})-\sum _{\mathfrak{u}\in {\mathfrak{U}}_{ϵ}}{Q}_{\mathfrak{u},n_{\mathfrak{u}}}(F_{\mathfrak{u}}({\mathbf{y}}_{\mathfrak{u}}))\right|$$

$${\mathcal{A}}_{\mathfrak{u}}^{\mathbb{N}}$$

$$F_{\mathfrak{u}}({\mathbf{y}}_{\mathfrak{u}})=0\text{when }\exists j\in \mathfrak{u}:y_j=0.$$

$$({\partial }_{{\mathbf{y}}_{\mathfrak{u}}}^{{\mathbf{\omega }}_{\mathfrak{u}}}{F}_{\mathfrak{u}})({\mathbf{y}}_{\mathfrak{u}})$$

$${\int }_{D}a(\mathbf{x},{\mathbf{y}}_{\mathfrak{v}})\nabla u(\mathbf{x},{\mathbf{y}}_{\mathfrak{v}})\cdot \nabla v(\mathbf{x})\mathrm{d}\mathbf{x}$$

$$a(\cdot ,{\mathbf{y}}_{\mathfrak{v}})$$

$$u(\cdot ,\mathbf{y})$$

$$u_s(\cdot ,y_1,\dots ,y_s,0,\dots )$$

$$I_{\mathfrak{u}}(F_{\mathfrak{u}})$$

Full text

We develop higher-order quasi-Monte Carlo rules for integration over the finite-dimensional Euclidean space with respect to the Gaussian distribution by use of a truncation strategy. By linear transformations of interlaced polynomial lattice rules from the unit cube to a multivariate box of the Euclidean space we achieve higher-order convergence rates for functions belonging to a class of anchored Gaussian Sobolev spaces while taking into account the truncation error. These cubature rules are then used in the MDFEM algorithm.

Under appropriate conditions, the MDFEM achieves higher-order convergence rates in terms of error versus cost, i.e., to achieve an accuracy of $O(\epsilon)$ the computational cost is $O(\epsilon^{-1/\lambda-{d^{\prime}}/\lambda})=O(\epsilon^{-(p^{*}+{d^{\prime}}/\tau)/(1-p^{*})})$ where $\epsilon^{-1/\lambda}$ and $\epsilon^{-{d^{\prime}}/\lambda}$ are respectively the cost of the quasi-Monte Carlo cubature and the finite element approximations, with ${d^{\prime}}=d\,(1+{\delta^{\prime}})$ for some ${\delta^{\prime}}\geq 0$ and $d$ the physical dimension, and $0<p^{*}\leq(2+{d^{\prime}}/\tau)^{-1}$ is a parameter representing the sparsity of $\{\phi_{j}\}_{j\geq 1}$.

Keywords: elliptic PDE, stochastic diffusion coefficient, lognormal case, infinite-dimensional integration, multivariate decomposition method, finite element method, higher-order quasi-Monte Carlo, high dimensional quadrature/cubature, complexity bounds.

1 Introduction

In this paper we are concerned with the application of higher-order quasi-Monte Carlo (QMC) rules and multivariate decomposition methods (MDM) to elliptic PDEs with random diffusion coefficients. We focus on the lognormal diffusion coefficient, the logarithm of which is a Gaussian random field. The goal is to compute the expected value of some functional of the solution. This method was motivated by the need for new techniques for elliptic PDEs with smooth lognormal diffusion coefficients where the classical QMC approaches can only achieve first-order convergence, see, e.g., [14, 15, 18, 17].

The MDFEM was already analysed in the case of a uniform diffusion coefficient in [28], but in this case higher-order QMC rules are readily available for integration over the unit cube. For lognormal diffusions we cope with a more challenging problem since the expectation is taken with respect to the Gaussian distribution over an unbounded domain. Consequently, existing higher-order QMC algorithms are not directly applicable. To solve this problem we propose using a truncation method recently developed in [7], see also [27]. By exploiting the fast decay of the Gaussian distribution toward infinity, the Euclidean domain is truncated and the resulting integral is transformed to the unit cube using a linear transformation where suitable higher-order QMC rules can be applied. The proposed algorithm allows us to achieve higher-order convergence for sufficiently smooth integrands.

Let $D\subset\mathbb{R}^{d}$ be a bounded polygonal domain in $\mathbb{R}^{d}$, with typically $d=1,2$ or $3$, with boundary $\partial D$. We consider the following elliptic Dirichlet problem

[TABLE]

Here, the gradient operator $\nabla$ is taken with respect to ${\bm{x}}$ and $a:D\times\Omega^{\mathbb{N}}\to\mathbb{R}$ for some $\Omega\subseteq\mathbb{R}$.

We consider the case when ${\bm{y}}=\{y_{j}\}_{j\geq 1}$ is a sequence of parameters distributed on $\mathbb{R}^{\mathbb{N}}$ according to the product Gaussian measure $\mu=\bigotimes_{j\geq 1}\mathcal{N}(0,1)$, and the diffusion coefficient takes the form

[TABLE]

with

[TABLE]

where $\{\phi_{j}\}_{j\geq 1}$ is a suitable system of real-valued, bounded, and measurable functions.

Let us denote the natural numbers by $\mathbb{N}:=\{1,2,\ldots\}$, and $\mathbb{N}_{0}:=\{0,1,2,\ldots\}$. For any $s\in\mathbb{N}$ we use the shorthand notation $\{1:s\}$ to denote the set of indices $\{1,2,\ldots,s\}$. Let $G$ be a linear and bounded functional of the solution $u$. We are interested in computing the expected value of $G(u)$ with respect to the probability distribution $\mu$, i.e.,

[TABLE]

where for $\mathfrak{u}\subseteq\mathbb{N}$

[TABLE]

with $\rho_{\mathfrak{u}}({\bm{y}}_{\mathfrak{u}}):=\prod_{j\in\mathfrak{u}}\rho(y_{j})$. We note that, depending on what is most natural, we write $\mathrm{d}\mu({\bm{y}}_{\mathfrak{u}})$ or $\rho_{\mathfrak{u}}({\bm{y}}_{\mathfrak{u}})\,\mathrm{d}{\bm{y}}_{\mathfrak{u}}$ with the understanding that the product probability measure is always on all the variables of the integral.

The weak formulation of problem (1) is to find for a given ${\bm{y}}\in\Omega^{\mathbb{N}}$ the solution $u(\cdot,{\bm{y}})\in V:=H_{0}^{1}(D)$ such that

[TABLE]

The space $V=H_{0}^{1}(D)$ is equipped with the norm $\|v\|_{V}:=\|\nabla v\|_{L^{2}(D)}$.

Under some assumptions on the system $\{\phi_{j}\}_{j\geq 1}$ we have existence, uniqueness and an a priori estimate of the solution of the weak formulation by means of the Lax–Milgram lemma. For this we need to show a lower bound and an upper bound on $a({\bm{x}},{\bm{y}})$ for ${\bm{x}}\in D$, but since we have a lognormal field, it might be that $a({\bm{x}},{\bm{y}})$ is not bounded for certain values of ${\bm{y}}\in\mathbb{R}^{\mathbb{N}}$. However, under the following assumptions, which are standard, see, e.g., [1, 2, 18, 21], we can claim lower and upper bounds to hold for $\mu$ almost every ${\bm{y}}$. Let us first define, for some given space $X$ and $p\in[1,\infty)$, the space $L_{p,\rho}(\mathbb{R}^{\mathbb{N}};X)$ of all strongly measurable mappings $v:\mathbb{R}^{\mathbb{N}}\to X$ such that

[TABLE]

The following result is implied by [1, Section 2] and [18, Theorem 2, Proposition 3 and Corollary 6].

Proposition 1**.**

If there exists a positive sequence $\{b_{j}\}_{j\geq 1}$, with $0<b_{j}\leq 1$ for all $j$, such that

[TABLE]

and

[TABLE]

then it holds that the Gaussian field $Z$ and the lognormal field $a=\exp(Z)$ are elements of $L_{p,\rho}(\mathbb{R}^{\mathbb{N}};L^{\infty}(D))$ for every $p\in[1,\infty)$. Therefore, for $\mu$ almost every ${\bm{y}}\in\mathbb{R}^{\mathbb{N}}$ we have

[TABLE]

and, by the Lax–Milgram lemma, the solution $u(\cdot,{\bm{y}})$ then exists and is unique and satisfies

[TABLE]

and additionally, $a^{-1}_{\min}$ is an element of $L_{p,\rho}(\mathbb{R}^{\mathbb{N}})$ for every $p\in[1,\infty)$, such that

[TABLE]

for any $f\in V^{*}$ and $p\in[1,\infty)$.

Proposition 1 shows that for $\mu$ almost every ${\bm{y}}\in\mathbb{R}^{\mathbb{N}}$ there exists a unique solution to (5) where $\mu$ is the product Gaussian measure, and this stems from the measurability of $Z$ and $a=\exp(Z)$ in $L_{p,\rho}(\mathbb{R}^{\mathbb{N}},L^{\infty}(D))$, in particular for $p=1$. To make the “for $\mu$ almost every ${\bm{y}}$” more tangible we note that if we define the set $\Xi:=\{{\bm{y}}\in\mathbb{R}^{\mathbb{N}}:\sup_{j\geq 1}|y_{j}|\,b_{j}<\infty\}\subseteq\mathbb{R}^{\mathbb{N}}$, then for any ${\bm{y}}\in\Xi$

[TABLE]

under condition (6). From here we can show the wanted results in Proposition 1: $a_{\max}({\bm{y}})<\infty$ and $a_{\min}({\bm{y}})>0$ for all ${\bm{y}}\in\Xi$. It can be shown, see, e.g., [1], that the set $\Xi$ has full measure $\mu(\Xi)=1$ for $\{b_{j}\}_{j\geq 1}\in\ell^{p^{*}}$, hence “for $\mu$ almost every ${\bm{y}}$”.

To approximate the solution of (5) the standard approach is to truncate the infinite series (3) to $s$ terms. Let us denote the random field truncated to the first $s$ terms by $Z_{s}$ and the so-obtained truncated lognormal field by $a_{s}$, i.e.,

[TABLE]

and by $u_{s}$ the solution to (5) with $a=a_{s}$, then it is shown, under the assumptions of Proposition 1, in [1, 18] that $Z_{s}\to Z$ in $L^{\infty}(D)$, $a_{s}\to a$ in $L^{\infty}(D)$ and $u_{s}\to u$ in $V$ when $s\to\infty$ for almost every ${\bm{y}}\in\mathbb{R}^{\mathbb{N}}$. We note that all statements of Proposition 1 also hold for the truncation to $s$ terms, see, e.g., [18].

To compute $\mathbb{E}[G(u)]$, see (4), we cope with three computational challenges: the infinite number of variables of the integrand, the unboundedness of the integration domain, and the integrand involves solutions of a PDE. Let us discuss how to solve these problems in more detail.

First, to approximate the infinite-dimensional integral we will make use of the multivariate decomposition method which originated from the changing dimension method, see, e.g., [25, 31, 23]. The goal is to decompose the infinite-dimensional problem into multiple finite-dimensional ones. Let us illustrate the MDM for calculating the integral

[TABLE]

We use the anchored decomposition [26], which, for any finite $s$, states

[TABLE]

where ${\bm{y}}_{\mathfrak{u}}$ takes the values $y_{j}$ for $j\in\mathfrak{u}$ and [math] for those $j\notin\mathfrak{u}$. (Conditions and properties for $\mathcal{F}$ and the $F_{\mathfrak{u}}$ will be provided later, here we just want to illustrate the concept.) We note that $\mathcal{F}({\bm{y}}_{\mathfrak{v}})$ is often written as $\mathcal{F}([{\bm{y}}_{\mathfrak{v}};\bm{0}]_{\mathfrak{v}})$ or $\mathcal{F}(({\bm{y}}_{\mathfrak{v}};\bm{0}))$ or $\mathcal{F}(({\bm{y}}_{\mathfrak{v}},\bm{0}_{-\mathfrak{v}}))$ or $\mathcal{F}({\bm{y}}_{\mathfrak{v}},\bm{0}_{-\mathfrak{v}})$, but we use our convention to not clutter further expressions more than necessary. Thus, to obtain the projections $F_{\mathfrak{u}}$ we combine $\mathfrak{v}$-truncated evaluations of the function $\mathcal{F}$ in points where all components outside of $\mathfrak{v}\subseteq\mathfrak{u}$ are set to zero. Taking the limit, we write

[TABLE]

Instead of just truncating to the first $s$ dimensions, this decomposition allows us to be more selective and approximate the function $\mathcal{F}$, or its integral $\mathcal{I}(\mathcal{F})$ in our case, up to a certain accuracy by considering a set of important subsets $\mathfrak{U}_{\epsilon}\subset\mathbb{N}^{\mathbb{N}}$, called the active set, as follows

[TABLE]

This form is interesting when a bounded linear operation on $\mathcal{F}$, which would be costly in the number of “active variables”, could be approximated on $\mathcal{F}_{\mathfrak{U}_{\epsilon}}$ by wrapping the linear operation inside the sum where hopefully only sets $\mathfrak{u}$ of small cardinality appear. Note that a plain truncation to the first $s$ dimensions could be represented like this as well, but then the active set would contain all subsets of $\{1:s\}$, including the set $\{1:s\}$ itself, and it would be more efficient to just use $\mathcal{F}(y_{1},\ldots,y_{s},0,\ldots)$ in this case. The strength of the MDM is exactly that the active set can be chosen as to satisfy

[TABLE]

for a certain error request $\epsilon>0$ and appropriate cubature rules $Q_{\mathfrak{u},n_{\mathfrak{u}}}$. By assessing the relative contribution of specific sets of variables, for a given desired error, the MDM will decide which $\mathfrak{u}$ to include in the active set to approximate the infinite MDM sum. In Proposition 7 we will show that the sets in the active set have relatively low cardinalities, provided certain conditions on $\{\phi_{j}\}_{j\geq 1}$ are satisfied, and hence only relatively low-dimensional problems remain, which can be solved at small cost. As a comparison, for a certain error request $\epsilon$, the truncation strategy might need a truncation dimension of $s=1907$, since $y_{1907}$ still has large enough influence, thus needing a $1907$-dimensional cubature rule, while the analysis of the MDM approach might show it is sufficient to include the sets $\{1,1907\}$ and $\{2,1907\}$, along with a lot of other small sets, but never surpassing more than $7$-dimensional subproblems, see [10, Table 1] for examples. This is particularly useful when error bounds grow exponentially with the number of dimensions.

Next, to compute integrals over the Euclidean space, we exploit the fast decay of the Gaussian distribution to truncate the unbounded domain to bounded boxes. We then use a linear transformation to map the truncated integral into the unit cube. Finally, we apply existing higher-order quasi-Monte Carlo rules, in particular we apply interlaced polynomial lattice rules, see, e.g., [13], to approximate the truncated integral. Interlaced polynomial lattice rules have been used before in the PDE context to achieve higher-order convergence but with uniform diffusion in [8]. In contrast, most existing QMC methods for integration with respect to the normal density, and hence most existing QMC based methods for estimating expected values in the PDE context with lognormal diffusion, map the integral to the unit cube by using the inverse of the Gaussian cumulative distribution function and then use randomly shifted lattice rules to approximate this integral, see, e.g., [24, 29], and [14, 15, 18, 17, 22]. The aim of all these methods is to obtain dimension-independent convergence by making use of weighted function spaces. However, since using the inverse of the Gaussian cumulative distribution function might damage the smoothness of the integrand, they make use of first order QMC rules, namely, randomly shifted lattice rules, and limit themselves to first order convergence. The proposed QMC method in this paper avoids damaging the smoothness of the integrand by using a truncation and mapping strategy. As a result, it allows us to achieve higher-order convergence rates for sufficiently smooth integrands on the Euclidean space, since they retain their smoothness on the unit cube. We do not include weights in our function spaces and therefore the constants depend exponentially on the number of dimensions. This is by design, since by using the MDM we only need to tackle relatively low-dimensional integrals instead of approximating the truncated high-dimensional integral directly. The importance of the subproblems will determine which ones to include and how well to approximate them.

Lastly, for each variable ${\bm{y}}$ sampled by the QMC method the original stochastic PDE becomes a deterministic one. To solve each such problem we use the finite element method (FEM). We call the combination of the multivariate decomposition method (MDM) with the finite element method (FEM) the multivariate decomposition finite element method or MDFEM in short.

For the further discussion we need some notation and properties of the anchored decomposition. In general our multi-indices and variables are infinite-dimensional, e.g., ${\bm{\omega}}\in\mathbb{N}^{\mathbb{N}}$ and ${\bm{y}}\in\mathbb{R}^{\mathbb{N}}$. By a subscript $\mathfrak{u}\subseteq\mathbb{N}$ we mean to only retain those components in $\mathfrak{u}$ and set the other components to zero. We extend this notation to sets and for a set $\mathcal{A}$ we define

[TABLE]

In the case that we want the components of a multi-index or vector to range over all values of a set $\mathcal{A}$ for the components in $\mathfrak{u}$, and be zero otherwise, we write ${\bm{y}}\in\mathcal{A}^{\mathbb{N}}_{\mathfrak{u}}$, and in this case ${\bm{y}}_{\mathfrak{u}}={\bm{y}}$. When we only need the components $y_{j}$ for $j\in\mathfrak{u}$ we then write, with a slight abuse of notation, ${\bm{y}}_{\mathfrak{u}}\in\mathcal{A}^{|\mathfrak{u}|}$ or even ${\bm{y}}_{\mathfrak{u}}\in\mathcal{A}^{\mathbb{N}}_{\mathfrak{u}}$ instead of ${\bm{y}}\in\mathcal{A}^{\mathbb{N}}_{\mathfrak{u}}$. If there is possible confusion then we will indicate the meaning explicitly or use a more explicit notation. Note that, contrary to typical usage, we occasionally also write ${\bm{\omega}}_{\mathfrak{u}}\in\mathbb{N}_{0}^{|\mathfrak{u}|}$, with $\omega_{j}\in\mathbb{N}_{0}=\{0,1,2,\ldots\}$ for $j\in\mathfrak{u}$, which means that we are working with the indices in $\mathfrak{u}$, but also allow $\omega_{j}=0$ for $j\in\mathfrak{u}$ in this case.

We write $\partial^{{\bm{\omega}}_{\mathfrak{u}}}_{{\bm{y}}_{\mathfrak{u}}}:=\partial^{|{\bm{\omega}}_{\mathfrak{u}}|}/\prod_{j\in\mathfrak{u}}\partial^{\omega_{j}}_{y_{j}}$ with $|{\bm{\omega}}_{\mathfrak{u}}|:=\sum_{j\in\mathfrak{u}}\omega_{j}$ and by $(\partial^{{\bm{\omega}}_{\mathfrak{u}}}_{{\bm{y}}_{\mathfrak{u}}}F)({\bm{y}})$ we mean the value of such a partial derivative at ${\bm{y}}$. Whenever it is clear w.r.t. which variables the derivatives are taken then we just write $F^{({\bm{\omega}}_{\mathfrak{u}})}$ or $F^{({\bm{\omega}})}$.

Lemma 1**.**

Given the anchored decomposition of a function $\mathcal{F}$ (with anchor at zero) by (13) and (14), then $F_{\mathfrak{u}}$ depends only on the variables listed in $\mathfrak{u}$ and satisfies

[TABLE]

Furthermore, if $\mathcal{F}$ has continuous partial derivatives up to $\partial^{{\bm{\omega}}_{\mathfrak{u}}}_{{\bm{y}}_{\mathfrak{u}}}$ for some ${\bm{\omega}}_{\mathfrak{u}}\in\mathbb{N}_{0}^{|\mathfrak{u}|}$ then

[TABLE]

and

[TABLE]

Proof.

See Appendix A.1. ∎

In the following the function $\mathcal{F}$ will actually be $G(u)$ and the projections $F_{\mathfrak{u}}$ will be $G(u_{\mathfrak{u}})$. We remark that the properties of Lemma 1 will hold for the decompositions of both $u$ and $G(u)$. When we describe function spaces for the functions $F_{\mathfrak{u}}$ we will identify them with $|\mathfrak{u}|$-variate functions on $\mathbb{R}^{|\mathfrak{u}|}$.

The structure of this paper is as follows. Section 2 presents the main steps of the MDFEM. Section 3 introduces higher-order QMC rules for multivariate integration over the Euclidean space with respect to the Gaussian distribution based on a truncation strategy. A novel anchored Gaussian Sobolev function space is introduced and QMC rules are developed for that specific space. Section 4 discusses the parametric regularity of the solution of the PDE and shows that $u_{\mathfrak{u}}$ and $G(u_{\mathfrak{u}})$ live in a Bochner space based on the anchored Gaussian Sobolev space, and provides bounds on their norms. Section 5 states our assumptions on the higher-order convergence of the FE approximations. Section 6 presents the main contribution of this paper where the cost model, the construction and the complexity of the MDFEM algorithm are presented. Finally, a comparison with two existing methods, the QMCFEM [18] and MLQMCFEM [17], shows the benefit of the MDFEM.

2 Applying the MDM to PDEs

We now explain the ingredients of the MDFEM. Our main building block will be $u(\cdot,{\bm{y}}_{\mathfrak{v}})$. By $u(\cdot,{\bm{y}}_{\mathfrak{v}})$ we mean the weak solution of (1) with ${\bm{y}}={\bm{y}}_{\mathfrak{v}}$, and we call this the $\mathfrak{v}$-truncated solution. That is, for given ${\bm{y}}_{\mathfrak{v}}\in\mathbb{R}^{\mathbb{N}}_{\mathfrak{v}}=\{{\bm{y}}\in\mathbb{R}^{\mathbb{N}}:y_{j}=0\text{ for }j\notin\mathfrak{v}\}$, with $|\mathfrak{v}|<\infty$, the $\mathfrak{v}$-truncated solution $u(\cdot,{\bm{y}}_{\mathfrak{v}})\in V$ is the solution of the problem

[TABLE]

where

[TABLE]

This is not unlike the commonly used truncation strategy where one solves the problem (1) by truncating to the first $s$ dimensions, i.e., solving the problem for one (relatively large) set $\mathfrak{v}=\{1:s\}$, see (12). In contrast, the MDFEM will solve the problem for multiple (relatively small) sets $\mathfrak{v}$ and combine those results to obtain the $\mathfrak{u}$-projected solution $u_{\mathfrak{u}}$ (which combines all $\mathfrak{v}$-truncated solutions $u(\cdot,\cdot_{\mathfrak{v}})$ for all $\mathfrak{v}\subseteq\mathfrak{u}$), see (20) below, to approximate the expected value of a functional $G$ of the solution of (5) up to a given accuracy. A quantitative statement about the number of PDEs which need to be solved and the sizes of the sets $\mathfrak{u}$ to achieve a certain error is given later in Proposition 7.

Using the truncated solutions we obtain a multivariate decomposition of the full solution $u(\cdot,{\bm{y}})$ by means of the anchored decomposition, cf. (13) and (14),

[TABLE]

with $u(\cdot,{\bm{y}}_{\mathfrak{v}})$ being $\mathfrak{v}$-truncated solutions. Note that the full solution consists of $\mathfrak{u}$-projected solutions which in their turn consist of $\mathfrak{v}$-truncated solutions. The infinite sum in (20) will be truncated to the active set $\mathfrak{U}_{\epsilon}$. Compared to the more standard solution method of truncating the problem to the first $s$ dimensions, see (12), we remind the reader that, cf. (13),

[TABLE]

and hence the limit in (20) can be read as $u=\lim_{s\to\infty}u_{s}$. We will show in Section 4 that the functions $u_{\mathfrak{u}}$ and $G(u_{\mathfrak{u}})$ belong to the spaces $H_{\alpha,0,\rho,|\mathfrak{u}|}(\mathbb{R}^{|\mathfrak{u}|};V)$ and $H_{\alpha,0,\rho,|\mathfrak{u}|}(\mathbb{R}^{|\mathfrak{u}|})$, see Lemma 4, where the spaces $H_{\alpha,0,\rho,|\mathfrak{u}|}(\mathbb{R}^{|\mathfrak{u}|})$ are reproducing kernel Hilbert spaces which we introduce in Section 3 and we formally identify $\mathbb{R}^{\mathbb{N}}_{\mathfrak{u}}$ with $\mathbb{R}^{|\mathfrak{u}|}$ after a relabelling of the components. In Remark 1 we explain that they form an orthogonal decomposition of an infinite-variate reproducing kernel Hilbert space. We will denote integration for $F_{\mathfrak{u}}\in H_{\alpha,0,\rho,|\mathfrak{u}|}(\mathbb{R}^{|\mathfrak{u}|})$ by

[TABLE]

Convergence of the MDM decomposition (20) and equivalence of the infinite-dimensional integral with

[TABLE]

is then guaranteed under the conditions from [12] in the setting of applying the MDM to our PDE, see Remarks 1 and 4 for further details, and the conditions we ask in Section 6. Therefore, the convergence statements from [1, 18] for solving the $s$-truncated problem as in (12) will also hold in our case. Note that the active set $\mathfrak{U}_{\epsilon}$, which will be defined in (69), grows for $\epsilon\to 0$ and will include all subsets of $\{1:s_{\epsilon}\}$ for some $s_{\epsilon}$ for which $s_{\epsilon}\to\infty$ as $\epsilon\to 0$. (Indeed, for any finite set $\mathfrak{u}$ there is an $\epsilon$ such that $\mathfrak{u}\in\mathfrak{U}_{\epsilon}$.)

The above discussion was using the exact weak solution for each $\mathfrak{v}$-truncated problem which is typically not available and a numerical approximation will be computed using the FE method of which the approximation error will have to be taken into account in the error analysis. Let us define a family of finite-dimensional subspaces $V^{h}\subset V$, where $h>0$ is the mesh diameter, i.e., the largest element diameter over all elements of the FE mesh, and such that $V^{h}\subset V^{h^{\prime}}\subset V$ for $h^{\prime}<h$. The finite element approximation of the weak formulation of the $\mathfrak{v}$-truncated problem (18) for a given ${\bm{y}}_{\mathfrak{v}}$ is to find $u^{h}(\cdot,{\bm{y}}_{\mathfrak{v}})\in V^{h}$ such that the following equation holds

[TABLE]

We can now piece together the different parts of the MDFEM. To approximate (4) we make use of (22), replacing the integrals by cubature formulas and the solutions to the PDEs by FE approximations, and hence the MDFEM takes the form

[TABLE]

where $\mathfrak{U}_{\epsilon}$ is the active set, $Q_{\mathfrak{u},n_{\mathfrak{u}}}$ are cubature rules with nodes ${\bm{y}}_{\mathfrak{u}}^{(i)}$ and weights $w_{\mathfrak{u}}^{(i)}$, and

[TABLE]

We emphasize that in calculating $G(u_{\mathfrak{u}}^{h_{\mathfrak{u}}}(\cdot,{\bm{y}}_{\mathfrak{u}}))$ we are solving $2^{|\mathfrak{u}|}$ PDEs, all with the same FE mesh, indicated by the FE mesh diameter $h_{\mathfrak{u}}$ for all the $\mathfrak{v}$-truncated solutions in the formula above.

The computational cost of the MDFEM algorithm is given as

[TABLE]

There are three sources of error in the MDFEM approximation (24): the truncation error from truncating the infinite sum, the cubature errors in approximating the integrals, and the FEM errors in approximating the solutions to the PDEs. They are gathered into two terms as follows

[TABLE]

A sufficient condition to achieve an approximation error of at most $\epsilon$ is that both of these terms are less than $\epsilon/2$. This forces us to construct the active set $\mathfrak{U}_{\epsilon}$ such that the first term is bounded by $\epsilon/2$. For each $\mathfrak{u}\in\mathfrak{U}_{\epsilon}$ the FE space $V^{h_{\mathfrak{u}}}$ and the cubature rule $Q_{\mathfrak{u},n_{\mathfrak{u}}}$ are then chosen such that the second term is bounded by $\epsilon/2$ while minimizing the computational cost (27). This will be the strategy we follow in Section 6.

3 Higher-order quasi-Monte Carlo rules for finite-dimensional integration with respect to the Gaussian distribution using truncation for anchored integrand functions

In this section we consider quasi-Monte Carlo rules for approximating integrals over $\mathbb{R}^{s}$ with respect to the Gaussian distribution. Particularly, we are interested in computing $s$-dimensional integrals of the form

[TABLE]

where, with a slight abuse of notation, in this section $\rho$ is the product Gaussian distribution $\prod_{j=1}^{s}\rho(y_{j})$. In further usage we are calculating $I_{\mathfrak{u}}(F_{\mathfrak{u}})=I_{s}(F)$, where $I_{\mathfrak{u}}$ is the integral w.r.t. $\rho_{\mathfrak{u}}({\bm{y}}_{\mathfrak{u}})\,\mathrm{d}{\bm{y}}_{\mathfrak{u}}$, and the $s$ dimensions here will be a relabelling of the variables in ${\bm{y}}_{\mathfrak{u}}$ with the function $F=F_{\mathfrak{u}}=G(u^{h_{\mathfrak{u}}}_{\mathfrak{u}}(\cdot,{\bm{y}}_{\mathfrak{u}}))$ given by the anchored decomposition for some $\mathfrak{u}\subset\mathbb{N}$ such that $s=|\mathfrak{u}|<\infty$. We will show in Section 4, see Lemma 4, that the functions $G(u^{h_{\mathfrak{u}}}_{\mathfrak{u}}(\cdot,{\bm{y}}_{\mathfrak{u}}))$ belong to the function spaces which we introduce below.

To approximate the integral $I_{s}(F)$, we first truncate the Euclidean domain to a multidimensional bounded box, then use a linear mapping to transform the truncated integral into one over the unit cube, and finally approximate the integral over the unit cube using suitable cubature rules. More precisely, the truncated and transformed integrals have the following form

[TABLE]

for some $T>0$ and we define the mapping $\bm{T}:[0,1]^{s}\to[-T,T]^{s}$ by

[TABLE]

The resulting integral (30) is then approximated using an $n$-point QMC rule of the form

[TABLE]

where $\{{\bm{y}}^{(i)}\}_{i=0}^{n-1}$ are well chosen cubature points on the unit cube.

We note that it would be possible to truncate the box differently for each dimension, defining a box $[-T_{1},T_{1}]\times\cdots\times[-T_{s},T_{s}]$ as was done, e.g., in [27]. We do not pursue such a strategy here since in the application of the MDM we are not immediately making use of the different importances of the dimensions for each projected $F_{\mathfrak{u}}$.

3.1 Two reproducing kernel Hilbert spaces

We will introduce two reproducing kernel Hilbert spaces. The first function space is an anchored Sobolev space on the Euclidean space $\mathbb{R}^{s}$ with Gaussian measure, cf. (29), for which we will take into account that our integrand function $F$ is anchored. We will show in Section 4 that under appropriate conditions the functions $G(u_{\mathfrak{u}})$ belong to this first function space. The second function space is an unanchored Sobolev space on the unit cube to analyse our truncated and mapped integral, cf. the right hand side of (30). For this second function space we can then adjust techniques from [8] such that interlaced polynomial lattice rules will achieve higher-order convergence in approximating the integral. Contrary to most modern QMC results we do not introduce weights in the function spaces since we rely on the MDM to keep the number of dimensions limited. Therefore our error bounds will contain exponential factors in the number of dimensions. However, these weights will eventually show up when we handle the infinite-variate case, see Remark 1.

3.1.1 Anchored Gaussian Sobolev space for anchored functions $H_{\alpha,0,\rho,s}(\mathbb{R}^{s})$

We begin with introducing the univariate function space. For $\alpha\in\mathbb{N}$ the space $H_{\alpha,0,\rho}(\mathbb{R})$ consists of integrable functions over $\mathbb{R}$ with respect to the Gaussian distribution having absolutely continuous derivatives up to order $\alpha-1$ for any bounded interval and square integrable derivative of order $\alpha$ over $\mathbb{R}$ with respect to the Gaussian distribution and are anchored at [math]. Note that this allows us to use the Lebesgue version of the Fundamental Theorem of Calculus, and hence the Taylor theorem with integral remainder up to order $\alpha$. For $F,G\in H_{\alpha,0,\rho}(\mathbb{R})$ the anchored inner product is defined as

[TABLE]

Note that, because we use the anchored decomposition for the MDM, we only consider functions such that $F(0)=0$ and hence the term for $\tau=0$ which is normally there in the first sum is zero here. The associated norm is given by $\|\cdot\|_{H_{\alpha,0,\rho}(\mathbb{R})}:=\langle\cdot,\cdot\rangle^{1/2}_{H_{\alpha,0,\rho}(\mathbb{R})}$.

We note that another common choice of a Gaussian Sobolev space is the unanchored Gaussian Sobolev space or Hermite space, see, e.g., [20, 19, 7]. Instead of anchoring the values of the function and its derivatives up to order $\alpha-1$ at [math] as in (32), they are integrated out against the Gaussian distribution over $\mathbb{R}$. However, here we want to benefit from the anchored decomposition and therefore will use an anchored Sobolev space.

The anchored Gaussian Sobolev space for anchored functions $H_{\alpha,0,\rho}(\mathbb{R})$ is a reproducing kernel Hilbert space with kernel

[TABLE]

where $\mathds{1}\{X\}$ is the indicator function on $X$. Such a kernel for $\alpha=1$ was given in [29, Section 3.3], but we note that in this case the weight function in the formula for the inner-product and the reproducing kernel cannot be taken the same as in the integral (29) for $\alpha=1$, see [29, Table 1]; here we are in fact interested in $\alpha\geq 2$, but see also Remark 2. A slightly different kernel for an anchored Sobolev space over the Euclidean space with higher order smoothness, although without taking any weight function into account, was given in [30, Section 11.5.1]. For completeness we provide the full derivation of the reproducing kernel for inner products like (32) for general $\rho$ in Appendix A.2.

We define the multivariate space as the tensor product of the univariate spaces. The kernel of our space of anchored functions is then given by

[TABLE]

Note that this kernel itself also has the anchored property: if there is a $j\in\{1:s\}$ for which $x_{j}=0$ or $y_{j}=0$ then $K_{\alpha,0,\rho,s}({\bm{x}},{\bm{y}})=0$. The corresponding inner product is

[TABLE]

where, in the first line, $-\mathfrak{v}=\{1:s\}\setminus\mathfrak{v}$ and $({\bm{y}}_{\mathfrak{v}},\bm{0}_{-\mathfrak{v}})$ is a vector of $s$ variables such that $({\bm{y}})_{j}=y_{j}$ for $j\in\mathfrak{v}$ and [math] otherwise. The second line follows by our convention that ${\bm{y}}_{\mathfrak{v}}$ is a vector of the appropriate size which takes the value zero outside of $\mathfrak{v}$. In many references, see, e.g., [4], such inner products are usually written using a double sum, which here then takes the following form

[TABLE]

where $({\bm{\alpha}}_{\mathfrak{v}},\bm{\tau}_{-\mathfrak{v}})$ is a vector of $s$ variables such that the $j$th component equals $\tau_{j}$ for $j\in-\mathfrak{v}$ and $\alpha$ otherwise. The corresponding norm is given by $\|\cdot\|_{H_{\alpha,0,\rho,s}(\mathbb{R}^{s})}:=\langle\cdot,\cdot\rangle^{1/2}_{H_{\alpha,0,\rho,s}(\mathbb{R}^{s})}$. In Lemma 5 in Appendix A.2 we show how anchored functions in anchored spaces like $H_{\alpha,0,\rho,s}(\mathbb{R}^{s})$ can be represented by a Taylor series with integral remainder.

In analysing the error of the MDM algorithm we will need to be able to get estimates on integrals of functions $F_{\mathfrak{u}}\in H_{\alpha,0,\rho,|\mathfrak{u}|}(\mathbb{R}^{|\mathfrak{u}|})$. We will make use of the following properties of the reproducing kernel.

Proposition 2**.**

For $\alpha\in\mathbb{N}$ we have

[TABLE]

and, for $\mathfrak{u}\subset\mathbb{N}$, we have

[TABLE]

Hence for $F_{\mathfrak{u}}\in H_{\alpha,0,\rho,|\mathfrak{u}|}(\mathbb{R}^{|\mathfrak{u}|})$ we have

[TABLE]

such that $F_{\mathfrak{u}}\in H_{\alpha,0,\rho,|\mathfrak{u}|}(\mathbb{R}^{|\mathfrak{u}|})\subset L^{1}(\mathbb{R}^{|\mathfrak{u}|},\rho)$ are integrable w.r.t. the Gaussian measure $\rho_{\mathfrak{u}}$.

Moreover, if $2\leq\alpha<\infty$ then we also have

[TABLE]

Proof.

For $y>0$ we have

[TABLE]

The same bound holds for any $y<0$. Thus, using $\left(\sum_{j}|a_{j}|\right)^{1/2}\leq\sum_{j}|a_{j}|^{1/2}$ we have

[TABLE]

with the maximum for $\alpha\in\mathbb{N}$ achieved for $\alpha=3$. The bound on $I_{\mathfrak{u}}(F_{\mathfrak{u}})$ follows by using the reproducing property of the kernel

[TABLE]

To show the last claim we write

[TABLE]

with ${}_{2}F_{2}$ a generalized hypergeometric function. The sum over $r$ will stay finite when integrating against the normal density, similar like above. For the second part we have

[TABLE]

which is finite when $2\leq\alpha<\infty$. ∎

Remark 1*.*

Our discussion so far is on finite-variate spaces $H_{\alpha,0,\rho,|\mathfrak{u}|}(\mathbb{R}^{|\mathfrak{u}|})$ of anchored functions. We now show how they form an orthogonal decomposition for a space of $s$-variate functions (not just anchored functions) and we will then extend this to infinite-variate functions. The kernel for the $s$-variate weighted anchored Sobolev space can be expressed as, see, [26, Example 4.4],

[TABLE]

where $\{\gamma_{\mathfrak{u}}\}_{|\mathfrak{u}|<\infty}$ is a sequence of positive weights, where for $\mathfrak{u}=\emptyset$ we set $\gamma_{\emptyset}:=1$. For $\mathfrak{u}=\emptyset$ we have the space of constant functions and we define $K_{\alpha,0,\rho,|\emptyset|}:=1$ and $\|f_{\emptyset}\|_{K_{\alpha,0,\rho,|\emptyset|}}:=|f_{\emptyset}|$. The weights $\gamma_{\mathfrak{u}}$ in (39) model the importance of our subspaces and we will show in Section 6 that they can be taken of product form $\gamma_{\mathfrak{u}}=\prod_{j\in\mathfrak{u}}\gamma_{j}$ with $\gamma_{j}=\sqrt{2}\,b_{j}$. In Section 6 we will eventually demand that $\{b_{j}\}_{j\geq 1}\in\ell^{p^{*}}(\mathbb{N})$ for $p^{*}\in(0,1)$. Hence we know that $\sum_{j\geq 1}\gamma_{j}<\infty$. By taking the limit for $s\to\infty$ we obtain the reproducing kernel of the infinite-variate reproducing kernel Hilbert space $\mathcal{H}_{\alpha,\rho,\bm{\gamma}}(\mathbb{R}^{\mathbb{N}})$, see [12],

[TABLE]

with inner product

[TABLE]

and where $F_{\mathfrak{u}}$ and $G_{\mathfrak{u}}$ are obtained by the anchored decomposition of the functions $\mathcal{F}$ and $\mathcal{G}$. Since our kernels are the tensor products of a univariate kernel, the subspaces are orthogonal and their intersection only contains the zero function for $\mathfrak{u}\neq\mathfrak{v}$, see [12]. Technically, the kernel $\mathcal{K}$ is a reproducing kernel only when ${\bm{x}},{\bm{y}}\in\mathcal{Y}$ with $\mathcal{Y}:=\{{\bm{y}}\in\mathbb{R}^{\mathbb{N}}:\mathcal{K}({\bm{y}},{\bm{y}})<\infty\}$, but we can use [12, Condition (C3)],

[TABLE]

to show that $\mu(\mathcal{Y})=1$ for $\alpha\geq 2$ under the condition of product weights $\gamma_{\mathfrak{u}}=\prod_{j\in\mathfrak{u}}\gamma_{j}$ with $\sum_{j\geq 1}\gamma_{j}<\infty$. This can be shown using (38) and the same reasoning as in the proof of Proposition 6. This means $\mathcal{Y}$ can be replaced by $\mathbb{R}^{\mathbb{N}}$ in the “almost everywhere” sense for the infinite-variate reproducing kernel Hilbert space and we can write $\mathbb{R}^{|\mathfrak{u}|}$ as the domain for each subspace instead of the formal domain $\mathcal{Y}$. Finally, this condition allows us to claim convergence of $s$-truncated functions to the infinite-variate function in the Hilbert space and to claim equality for the MDM form of the infinite-dimensional integral

[TABLE]

see again [12]. We refer to Remark 4 for the verification that our integrand function $\mathcal{F}=G(u)$, with $u(\cdot,{\bm{y}})$ being the solution of the PDE, belongs to $\mathcal{H}_{\alpha,\rho,\bm{\gamma}}(\mathbb{R}^{\mathbb{N}})$ under the studied conditions.

3.1.2 Unanchored Sobolev space on the unit cube $H_{\alpha,s}([0,1]^{s})$

We need a function space over the unit cube to analyse the error for our cubature method of choice, which will be interlaced polynomial lattice rules. For this, let us define the unanchored Sobolev space over the unit cube $H_{\alpha,s}([0,1]^{s})$. This is the same space as was used in, e.g., [4]. For $\alpha\in\mathbb{N}$ the univariate space $H_{\alpha}([0,1])$ consists of integrable functions over $[0,1]$ having absolutely continuous derivatives up to order $\alpha-1$ and square integrable derivative of order $\alpha$. For the univariate case the unanchored inner product is

[TABLE]

with norm $\|\cdot\|_{H_{\alpha}([0,1])}:=\langle\cdot,\cdot\rangle^{1/2}_{H_{\alpha}([0,1])}$. Note that contrary to our Gaussian function space which is taking benefit of the functions being anchored here we must include the typical $\tau=0$ term in the first sum.

The multivariate space is the tensor product of the univariate spaces with inner product

[TABLE]

and norm $\|\cdot\|_{H_{\alpha,s}([0,1]^{s})}:=\langle\cdot,\cdot\rangle^{1/2}_{H_{\alpha,s}([0,1]^{s})}$.

Let us denote integration over the unit cube for functions $F\in H_{\alpha,s}([0,1]^{s})$ by

[TABLE]

and a QMC rule using a point set $P_{n}=\{{\bm{y}}^{(0)},\ldots,{\bm{y}}^{(n-1)}\}$ over the unit cube by

[TABLE]

For this space we can construct interlaced polynomial lattice rules which achieve the almost optimal order of convergence. The full derivation of the following result is given in Appendix A.4.

Theorem 1**.**

For $\alpha\in\mathbb{N}$, with $\alpha\geq 2$, let $F\in H_{\alpha,s}([0,1]^{s})$. For any $m\in\mathbb{N}$ an interlaced polynomial lattice rule of order $\alpha$ with point set $P_{n,\alpha}$ with $n=2^{m}$ points can be constructed with cost $O(\alpha sn\log(n))$ such that

[TABLE]

where

[TABLE]

Proof.

This result follows from combining Propositions 10 and 11 in the appendix. The construction cost follows from the algorithm for product weights, which we set all equal to $1$, in [8]. ∎

3.1.3 Norm in $H_{\alpha,s}([0,1]^{s})$ after mapping and truncating from $H_{\alpha,0,\rho,s}(\mathbb{R}^{s})$

To complete the analysis of mapping and truncating the integral (29) to (30) and using the result from Theorem 1 we need to show a bound on the norm in $H_{\alpha,s}([0,1]^{s})$ of the mapped function in terms of the norm in $H_{\alpha,0,\rho,s}(\mathbb{R}^{s})$ of the original function. Since the proof is quite long it is presented in the appendix.

Proposition 3**.**

For any $F\in H_{\alpha,0,\rho,s}(\mathbb{R}^{s})$ with $\alpha\in\mathbb{N}$ and $T\geq 1/(2\sqrt{2})$, the function $(F\rho)\circ\bm{T}:[0,1]^{s}\to\mathbb{R}^{s}$ belongs to $H_{\alpha,s}([0,1]^{s})$ and

[TABLE]

where

[TABLE]

and $I_{0}(\cdot)$ is the modified Bessel function of the first kind of order [math] with $I_{0}(1/2)\approx 1.06348$.

Proof.

See Appendix A.3. ∎

3.2 Higher-order quasi-Monte Carlo for integration over $\mathbb{R}^{s}$

We are now ready to state the error bound of the method described in the beginning of this section.

Theorem 2**.**

For $s\in\mathbb{N}$, $\alpha\in\mathbb{N}$, with $\alpha\geq 2$, let $F\in H_{\alpha,0,\rho,s}(\mathbb{R}^{s})$. Let $P_{n,\alpha}=\{{\bm{y}}^{(0)},\ldots,{\bm{y}}^{(n-1)}\}$ be the point set of an interlaced polynomial lattice rule of order $\alpha$ with $n=2^{m}$ points, with $m\in\mathbb{N}$, according to Theorem 1. Then, for any $\lambda\in[1,\alpha)$ and by taking $T=2+2\sqrt{\lambda\ln(n)}$, the QMC rule $Q_{s,n}$, defined in (31), using the point set $P_{n,\alpha}$ has an error bounded as

[TABLE]

where $C_{\alpha,\lambda,s}$ is a constant, independent of $F$ and $n$, defined below by (53).

Proof.

The error splits into two terms

[TABLE]

The first term is the domain truncation error. Similar as in Proposition 2, using the reproducing property of the kernel $K_{\alpha,0,\rho,s}$, we obtain

[TABLE]

To obtain a bound on the truncation error we will bound the integral on the right hand side. We have

[TABLE]

We will estimate each of the above integrals. Using the same argument as in the proof of Proposition 2 we have for any $T\geq 2$

[TABLE]

where we used $\alpha/\sqrt{2\alpha-1}\leq\sqrt{\alpha}$ for $\alpha\geq 1$, $\sqrt{\rho(y)}\leq 1/\sqrt{y}\leq 1/\sqrt{T}$ for $y\geq T$, and in the last step filled in $T=2$. Moreover, we have for $T>2$, and with the substitution $t=y/2-1$,

[TABLE]

Inserting this into (3.2) yields

[TABLE]

Applying this inequality together with (36) and (47) to (46) we obtain

[TABLE]

with

[TABLE]

We move to the second term of the total error which is the cubature error. Using the result of Theorem 1 and Proposition 3 we have for any $\lambda\in[1,\alpha)$

[TABLE]

where

[TABLE]

with $C_{1,\alpha}$ and $\widetilde{C}_{\alpha,\lambda,s}$, respectively, defined in (43) and in (42).

Combining (3.2), (49) and (51) leads to

[TABLE]

To balance the dominating terms in the square brackets we choose $T=2+2\sqrt{\lambda\ln(n)}$ such that $\mathrm{e}^{-(T/2-1)^{2}}=n^{-\lambda}$. Hence,

[TABLE]

where in the second inequality we used $2+2\sqrt{\lambda\ln(n)}\leq\left(\frac{2}{\sqrt{\ln 2}}+2\sqrt{\lambda}\right)\sqrt{\ln(n)}$ for $n\geq 2$ and

[TABLE]

with $C_{3,\alpha,s}$ and $C_{4,\alpha,\lambda,s}$, respectively, defined in (50) and (52). ∎

Remark 2*.*

A similar result as that of Theorem 2 can be shown for the anchored Gaussian Sobolev space with first order smoothness $H_{\alpha,0,\rho,s}(\mathbb{R}^{s})$ with $\alpha=1$ using randomly digitally shifted polynomial lattice rules which achieve the optimal convergence rate in $H_{\alpha,s}([0,1]^{s})$ with $\alpha=1$, see [9]. However, in this case we have to be slightly careful as (38) does not hold for $\alpha=1$. We note however that we can change from $y_{j}\sim\mathcal{N}(0,1)$ to $y_{j}\sim\mathcal{N}(0,c^{2})$ with $0<c<1$ in (3), in effect changing the variance of the normal distribution that we integrate against. To keep the law of the random field unchanged we will have to divide each $y_{j}$ by $c$ in (3). If we keep the weight in the kernel unchanged then the integral (38) w.r.t. $\mathcal{N}(0,c^{2})$ will be finite, see also [29, Table 1]. The effect of dividing $y_{j}$ by $c$ can now be moved into the basis functions $\phi_{j}$ and eventually ends up as multiplying each $b_{j}$ with $c$. Now take $c=1-\delta$ for arbitrarily small $\delta>0$. For condition (6) we then obtain $\kappa_{c}=\kappa/c\leq(1+\delta)\,\kappa<\infty$ if before we had $\kappa<\infty$. The same remark holds for condition (57) which asks $\kappa<\ln(2)/\alpha$ and which we will need in the next section. We here obtain $\kappa_{c}\leq(1+\delta)\,\kappa<\ln(2)/\alpha$ if before we had $\kappa<\ln(2)/\alpha$. Combining such randomized cubature rules with a suitable truncation of the Euclidean domain gives a similar convergence rate as in (44), however, of order $\lambda\in[\frac{1}{2},1)$.

Remark 3*.*

An alternative approach is to embed the function $(F\rho)\circ\bm{T}$ into the anchored Sobolev space over the unit cube and then use higher-order polynomial lattice rules, see, e.g., [6, 28]. However, a non-trivial result similar to Proposition 3 is then needed to obtain an explicit formula for the embedding constant. One of the reasons we choose our approach of mapping to the unanchored Sobolev space is that then the technique of [8] to make use of interlaced polynomial lattice rules could be used with a construction cost of $O(\alpha sn\log(n))$ as given in Theorem 1. The construction cost of higher-order polynomial lattice rules on the other hand grows exponentially in $n$ with respect to the smoothness, see [5]. We remark that this excessive construction cost could still be avoided by making use of yet another embedding of the anchored Sobolev space over the unit cube into the unanchored Sobolev space on the unit cube, see [11, Example 2.1]. In that case the interlaced polynomial lattice rules from Theorem 1 could then also be used.

4 Parametric regularity of the PDE solution

We next derive bounds for mixed derivatives of the solution $u(\cdot,{\bm{y}})$ with respect to ${\bm{y}}$. For $\alpha\in\mathbb{N}$ and $\mathfrak{u}\subset\mathbb{N}$ let us define the Bochner norm based on the spaces $H_{\alpha,0,\rho,|\mathfrak{u}|}(\mathbb{R}^{|\mathfrak{u}|})$, with inner product (3.1.1), and $V$ by

[TABLE]

Note that $\partial^{\bm{\tau}_{\mathfrak{u}}}_{{\bm{y}}_{\mathfrak{u}}}u_{\mathfrak{u}}$ is a function which only depends on the variables in $\mathfrak{u}$ and inside the integral we evaluate this function at ${\bm{y}}_{\mathfrak{v}}$ with $\mathfrak{v}\subseteq\mathfrak{u}$ setting all $y_{j}=0$ for $j\notin\mathfrak{v}$. The following result was also used in [28, Lemma 2] for the analysis of the MDFEM in the uniform case and allows us to use the regularity analysis on $u(\cdot,\cdot_{\mathfrak{u}})$ instead of on $u_{\mathfrak{u}}$, since their norms in $H_{\alpha,0,\rho,|\mathfrak{u}|}(\mathbb{R}^{|\mathfrak{u}|};V)$ coincide. Note however that $u(\cdot,\cdot_{\mathfrak{u}})\notin H_{\alpha,0,\rho,|\mathfrak{u}|}(\mathbb{R}^{|\mathfrak{u}|};V)$ since it does not satisfy the anchored properties, except maybe in exceptional cases.

Lemma 2**.**

For $\alpha\in\mathbb{N}$ and $\mathfrak{u}\subset\mathbb{N}$, let $u_{\mathfrak{u}}$ be obtained by the anchored decomposition (20) and let $u(\cdot,\cdot_{\mathfrak{u}})$ be the $\mathfrak{u}$-truncated solution, cf. (18). Let $G$ be a bounded linear functional on $V$ such that $|G(v)|\leq\|G\|_{V^{*}}\,\|v\|_{V}$ for all $v\in V$. Then it holds that

[TABLE]

Proof.

This follows directly from property (16) of the anchored decomposition and (54). ∎

For a given ${\bm{y}}_{\mathfrak{u}}\in\mathbb{R}^{\mathbb{N}}_{\mathfrak{u}}$ and with $v(\cdot,{\bm{y}}_{\mathfrak{u}})\in V$ let us introduce the notation

[TABLE]

Note that $\|\cdot\|_{V,a_{{\bm{y}}_{\mathfrak{u}}}}$ depends on ${\bm{y}}_{\mathfrak{u}}$. It is easy to see that for every $v(\cdot,{\bm{y}}_{\mathfrak{u}})\in V$

[TABLE]

and additionally, when $v(\cdot,{\bm{y}}_{\mathfrak{u}})=u(\cdot,{\bm{y}}_{\mathfrak{u}})\in V$ is the $\mathfrak{u}$-truncated solution, cf. (18), we obtain

[TABLE]

where $f\in V^{*}$ is the right hand side of the PDE and where we used (10) for ${\bm{y}}={\bm{y}}_{\mathfrak{u}}$.

The following result is modified from [21, Proposition 3.1], see also [1, Theorem 4.1], and accounts for the truncation to an arbitrary set $\mathfrak{u}\subset\mathbb{N}$.

Proposition 4**.**

Given $\alpha\in\mathbb{N}$, if there exists a sequence $\{b_{j}\}_{j\geq 1}$ with $0<b_{j}\leq 1$ for all $j$ and the constant $\kappa$, defined in (6), satisfies

[TABLE]

then, for $\mu$ almost every ${\bm{y}}\in\mathbb{R}^{\mathbb{N}}$, such that $u(\cdot,{\bm{y}}_{\mathfrak{u}})$ is the solution of the $\mathfrak{u}$-truncated problem, it holds that

[TABLE]

with

[TABLE]

and $C_{\kappa,\alpha}:=\sum_{k=0}^{\infty}\delta_{\kappa,\alpha}^{k}<\infty$, with $0<\delta_{\kappa,\alpha}<1$ being a constant depending on $\kappa$ and $\alpha$ such that $\kappa<\delta_{\kappa,\alpha}\ln(2)/\alpha$. This implies

[TABLE]

Proof.

The last inequality in the statement follows by using $0<b_{j}\leq 1$ which implies $\bm{b}_{\mathfrak{u}}^{-2}\leq\bm{b}_{\mathfrak{u}}^{-2\bm{\tau}_{\mathfrak{u}}}$ for every $\bm{\tau}_{\mathfrak{u}}\in\{1:\alpha\}^{|\mathfrak{u}|}$ together with (55) and (56). ∎

The previous result can now be used to show a bound on the norm of $u(\cdot,\cdot_{\mathfrak{u}})$ and as a consequence also on the norm of $u_{\mathfrak{u}}$.

Lemma 3**.**

Assume the sequence $\{b_{j}\}_{j\geq 1}$ satisfies the assumptions of Proposition 4 for a given $\alpha\in\mathbb{N}$, and, additionally that $\{b_{j}\}_{j\geq 1}\in\ell^{p^{*}}(\mathbb{N})$ for some $p^{*}\in(0,1]$. Then it holds that

[TABLE]

with $\bm{b}_{\mathfrak{u}}=\prod_{j\in\mathfrak{u}}b_{j}$ and

[TABLE]

Proof.

Using the definition of the Bochner norm (54), but with the $H_{\alpha,0,\rho,|\mathfrak{u}|}$ norm written as a double sum, cf. (3.1.1), we have

[TABLE]

where we applied Proposition 4 with ${\bm{y}}_{\mathfrak{u}}=({\bm{y}}_{\mathfrak{v}},\bm{0}_{\mathfrak{u}\setminus\mathfrak{v}})\in\mathbb{R}^{\mathbb{N}}_{\mathfrak{v}}\subseteq\mathbb{R}^{\mathbb{N}}_{\mathfrak{u}}$ for each $\mathfrak{v}\subseteq\mathfrak{u}$. Now we estimate the sum in the last expression. Our strategy is similar to that in [18, proof of Theorem 13]. Note that for any ${\bm{y}}_{\mathfrak{v}}\in\mathbb{R}^{\mathbb{N}}_{\mathfrak{v}}$ we have

[TABLE]

Therefore, we have

[TABLE]

where we used that the integral in the second line can be interpreted as the cumulative standard normal distribution evaluated at $2\kappa b_{j}$, and this can be bounded by $\exp(2(2\kappa b_{j})/\sqrt{2\pi})/2$, see, e.g.,  [14, p. 355]. Inserting (4) into (4) we obtain

[TABLE]

Since $\{b_{j}\}_{j\geq 1}\in\ell^{p^{*}}(\mathbb{N})\subseteq\ell^{1}(\mathbb{N})\subset\ell^{2}(\mathbb{N})$ for $p^{*}\in(0,1]$ the sum in the last expression is finite. Taking the square root on both sides finishes the proof. ∎

Combining Lemma 2 and Lemma 3 we obtain bounds for the norms of $u_{\mathfrak{u}}$ and $G(u_{\mathfrak{u}})$. Note that the arguments used to arrive at these bounds are based on the weak formulation with $u(\cdot,{\bm{y}}_{\mathfrak{u}})\in V$. They remain valid for the approximation $u_{\mathfrak{u}}^{h_{\mathfrak{u}}}$, in which we combine the FE approximations $u^{h_{\mathfrak{u}}}(\cdot,{\bm{y}}_{\mathfrak{v}})\in V^{h_{\mathfrak{u}}}$ for all $\mathfrak{v}\subseteq\mathfrak{u}$, since $V^{h_{\mathfrak{u}}}\subset V$, with constants independent of $h_{\mathfrak{u}}$. Note that this requires us to use the same FE mesh diameter $h_{\mathfrak{u}}$ for all the $\mathfrak{v}$-truncated solutions which we use to calculate $u_{\mathfrak{u}}^{h_{\mathfrak{u}}}$, see (25).

Lemma 4**.**

Assume the sequence $\{b_{j}\}_{j\geq 1}$ satisfies the assumptions of Proposition 4 for a given $\alpha\in\mathbb{N}$, and, additionally that $\{b_{j}\}_{j\geq 1}\in\ell^{p^{*}}(\mathbb{N})$ for some $p^{*}\in(0,1]$. Let $G\in V^{*}$ be a bounded linear functional. Then it holds that $u_{\mathfrak{u}}\in H_{\alpha,0,\rho,|\mathfrak{u}|}(\mathbb{R}^{|\mathfrak{u}|};V)$ with

[TABLE]

with $C_{\kappa,\alpha}^{\prime}$ given in (58). Moreover, for $V^{h_{\mathfrak{u}}}\subset V$, we have $u_{\mathfrak{u}}^{h_{\mathfrak{u}}}\in H_{\alpha,0,\rho,|\mathfrak{u}|}(\mathbb{R}^{|\mathfrak{u}|};V)$ and $G(u_{\mathfrak{u}}^{h_{\mathfrak{u}}})\in H_{\alpha,0,\rho,|\mathfrak{u}|}(\mathbb{R}^{|\mathfrak{u}|})$ with the same bounds on their norms as given for $u_{\mathfrak{u}}$ and $G(u_{\mathfrak{u}})$ up to a multiplicative constant independent of $h_{\mathfrak{u}}$.

5 Finite element approximation error

For bounding the error we also need results on the FE approximation error, which is the error between the true solution $u(\cdot,{\bm{y}}_{\mathfrak{v}})\in V$ and the discretized solution $u^{h_{\mathfrak{u}}}(\cdot,{\bm{y}}_{\mathfrak{v}})\in V^{h_{\mathfrak{u}}}\subset V$ of the weak form, cf. (18) and (23), where (25) dictates which FE approximations we need to take. The FE approximation error depends on the spatial regularity of the solution which depends on the smoothness of the domain $D$, the right hand side of the PDE $f$ and the spatial regularity of the diffusion coefficient $a(\cdot,{\bm{y}}_{\mathfrak{v}})=\exp(Z(\cdot,{\bm{y}}_{\mathfrak{v}}))$.

Our assumptions are quite standard, see, e.g., [18, Proposition 15] and also [8, Theorem 2.4 and Theorem 2.5]. We will assume that, given

[TABLE]

there exists a $t\in(0,\infty)$ such that

[TABLE]

and for which there exists a positive sequence $\{b_{j}\}_{j\geq 1}$, with $0<b_{j}\leq 1$ for all $j$, satisfying condition (6) and (7), which will be further strengthened, see Theorem 3, as well as

[TABLE]

then, for all $p\in[1,\infty)$ and $\tau<2t$ there is a constant $C>0$ such that, for any $h_{\mathfrak{v}}>0$, we have an FEM, using a continuous piecewise polynomial basis of total degree $r\geq\lceil\tau/2\rceil$, for which

[TABLE]

Here $C$ depends on the spatial regularity of $D$, $a$, $f$ and $G$ through (61)–(63), but is independent of $h_{\mathfrak{v}}$ and $\mathfrak{v}$. We note that to handle singularities, due to, e.g., reentrant corners for $d=2$, one either has to change the norms to incorporate a weight function as in [17, Proposition 2.3 and Remark 2.4] or use local mesh refinement as, e.g., in [15]. We refer to [17] for the definitions of the norms and omit such details here.

6 MDFEM

In the next subsections we explain the cost model for the MDFEM algorithm and select the active set, cubature rules and finite element approximations based on a priori error estimates. The main complexity result will be presented in Theorem 3, after which we compare the MDFEM with the QMCFEM and MLQMCFEM algorithms.

6.1 Computational cost

The total cost of the MDFEM algorithm (24) is comprised of the costs of computing $Q_{\mathfrak{u},n_{\mathfrak{u}}}(G(u^{h_{\mathfrak{u}}}_{\mathfrak{u}}))$ for all $\mathfrak{u}\in\mathfrak{U}(\epsilon)$, where $\mathfrak{U}(\epsilon)$ is the “active set” determined to reach a given error request $\epsilon>0$. For each $\mathfrak{u}\in\mathfrak{U}(\epsilon)$ the cost of computing $Q_{\mathfrak{u},n_{\mathfrak{u}}}(G(u^{h_{\mathfrak{u}}}_{\mathfrak{u}}))$ is given by $n_{\mathfrak{u}}$ times the cost of computing $G(u^{h_{\mathfrak{u}}}_{\mathfrak{u}})$, see (27). Based on the decomposition (26) the cost of evaluating $G(u^{h_{\mathfrak{u}}}_{\mathfrak{u}})$ is bounded by $2^{|\mathfrak{u}|}$ times the cost of evaluating the FE approximation of the $\mathfrak{u}$-truncated solution, where we assume the cost of approximating $u^{h_{\mathfrak{u}}}(\cdot,{\bm{y}}_{\mathfrak{u}})$ to be dominating those of $u^{h_{\mathfrak{u}}}(\cdot,{\bm{y}}_{\mathfrak{v}})$ for all $\mathfrak{v}\subseteq\mathfrak{u}$ since the stiffness matrix involves calculating (19) at a cost which we will estimate at $O(|\mathfrak{v}|)$. Technically this cost could be avoided by using a Gray code ordering of enumerating the sets $\mathfrak{v}\subseteq\mathfrak{u}$. Hence, we have

[TABLE]

For each ${\bm{y}}_{\mathfrak{u}}$, and hence for each ${\bm{y}}_{\mathfrak{v}}$ obtained from this ${\bm{y}}_{\mathfrak{u}}$ for $\mathfrak{v}\subseteq\mathfrak{u}$, the cost of evaluating the FE approximation $u^{h_{\mathfrak{u}}}(\cdot,{\bm{y}}_{\mathfrak{v}})$ is given by the cost of assembling the stiffness matrix plus the cost of solving the linear system. Due to the locality of the $O(h_{\mathfrak{u}}^{-d})$ basis functions of $V^{h_{\mathfrak{u}}}$, the stiffness matrix is sparse and has $O(h_{\mathfrak{u}}^{-d})$ nonzero entries. Each entry in turn needs at most $O(|\mathfrak{u}|)$ operations to evaluate the diffusion parameter, cf. (19), which could be avoided, see the previous remark, but we will leave this cost in, since there are other terms more dominating, cf. (80) forthcoming. We assume

[TABLE]

E.g., in [18, Section 10] it is assumed that the cost is nearly linear and that ${\delta^{\prime}}>0$ can be chosen arbitrarily small, making use of [16, Corollary 17]. Thus, using $(m+1)\leq 2m$ for $m\in\mathbb{N}$, we have

[TABLE]

To simplify the further exposition, we will write $\pounds_{\mathfrak{u}}:=2^{|\mathfrak{u}|}|\mathfrak{u}|$ and ${d^{\prime}}:=d\,(1+{\delta^{\prime}})$. Therefore, the total computational cost of the MDFEM is

[TABLE]

6.2 Error analysis

We will now give an a priori estimate of the total error of the MDFEM algorithm. We will use the higher-order convergence of our QMC based cubature rules from Section 3, together with the required parametric regularity results of the PDE solution from Section 4, and the convergence of the FE approximations from Section 5. For the convergence of the FE approximation we remind the reader that similar assumptions as given in (61)–(63), possibly extended with weighted norms or local mesh refinements, are needed to allow the convergence of (64) to hold.

The next result holds for any choice of $\mathfrak{U}(\epsilon)$ and with the MDFEM algorithm $Q_{\epsilon}$ making use of this active set. The actual choice of the active set, as well as the cubature rules and FE approximations, to reach a certain error request $\epsilon$ will be shown in the next sections based on the error bound in the next result.

We remark that for $\mathfrak{u}=\emptyset$ there is no integral to approximate since $I_{\emptyset}(F_{\emptyset})=F_{\emptyset}$ with $F_{\emptyset}=\mathcal{F}(\bm{0})$, which requires a single function evaluation to be computed. Hence, for $n_{\emptyset}=0$ the absolute value of the cubature error is $\|F_{\emptyset}\|_{H_{\alpha,0,\rho,0}}=|F_{\emptyset}|$, while for $n_{\emptyset}\geq 1$, we set $Q_{\emptyset,n_{\emptyset}}(F_{\emptyset})=F_{\emptyset}=\mathcal{F}(\bm{0})$ and the cubature error is [math]. To cover this case easily we define $0^{0}=1$ and $C_{\alpha,\lambda,0}:=1$.

Proposition 5**.**

Let the set $\mathfrak{U}(\epsilon)\subset\mathbb{N}^{\mathbb{N}}$ be given and assume that the conditions of Lemma 4 are satisfied for a given $\alpha\in\mathbb{N}$, and hence there is a sequence $\{b_{j}\}_{j\geq 1}\in\ell^{p^{*}}(\mathbb{N})$ for some $p^{*}\in(0,1]$, with $0<b_{j}\leq 1$ and $\kappa<\ln(2)/\alpha$. Let $G\in V^{*}$ be a bounded linear functional. Assume $D$, $a$, $f$ and $G$ have sufficient spatial regularity such that the application of $G$ to the FE approximations of the truncated solutions converge like $O(h^{\tau})$ as in (64). Let the cubature rules be defined as the transformed interlaced polynomial lattice rules with interlacing factor $\alpha$ as in Theorem 2, and hence $\alpha\geq 2$, and, with $m_{\mathfrak{u}}\in\mathbb{N}$, using $n_{\mathfrak{u}}=2^{m_{\mathfrak{u}}}\geq 2$ points. For $n_{\mathfrak{u}}=0$ and $n_{\mathfrak{u}}=1$ we take the zero approximation. Then the error of the MDFEM algorithm $Q_{\epsilon}$, see (24), based on the given set $\mathfrak{U}(\epsilon)$ can be bounded as

[TABLE]

for any $\lambda\in[1,\alpha)$, with the corresponding truncation points for the cubature rules chosen in accordance with this $\lambda$, i.e., $T_{\mathfrak{u}}=2+2\sqrt{\lambda\ln(n_{\mathfrak{u}})}$, and where $\alpha_{1}=\alpha/2+1/4$, $C_{\alpha,\lambda,|\mathfrak{u}|}$ is given by (53), $M_{\mathfrak{u}}=M^{|\mathfrak{u}|}$ with $M$ given by (36) and $\gamma_{\mathfrak{u}}:=\prod_{j\in\mathfrak{u}}\gamma_{j}$ with

[TABLE]

and $\gamma_{\emptyset}:=1$. The first term in (67) corresponds to the truncation error, while the second term corresponds to the combined error due to the FE approximations and QMC cubatures. The hidden constant in (67) is independent of the choice of the active set $\mathfrak{U}(\epsilon)$ and the choices of $n_{\mathfrak{u}}$ and $h_{\mathfrak{u}}$.

Proof.

The error splits into three terms:

[TABLE]

The first term is the truncation error from truncating to only a finite number of decomposed elements. The second term stems from the spatial discretization of the FE approximations. The last term is the cubature error arising from using cubature rules to approximate the integrals.

1.

Since we know that $G(u_{\mathfrak{u}})\in H_{\alpha,0,\rho,|\mathfrak{u}|}(\mathbb{R}^{|\mathfrak{u}|})$, see Lemma 4, we can make use of (37) in Proposition 2. Therefore, with $\{\gamma_{\mathfrak{u}}\}_{|\mathfrak{u}|<\infty}$ a sequence of positive weights, the truncation error can be bounded as

[TABLE]

2.

For the error due to the FE approximations we have

[TABLE]

Due to the linearity of $G$ we have

[TABLE]

Thus, using the assumption of an FE approximation error bound as in (64), we have

[TABLE]

with the hidden constant independent of $h_{\mathfrak{u}}$ and $\mathfrak{u}$ but dependent on the spatial regularity conditions of $D$, $a$, $f$ and $G$, see, e.g., (61)–(63). Hence the error incurred by the FE approximation can be bounded as

[TABLE]

with the hidden constant independent of $h_{\mathfrak{u}}$.

3.

From Lemma 4 we also know $G(u_{\mathfrak{u}}^{h_{\mathfrak{u}}})\in H_{\alpha,0,\rho,|\mathfrak{u}|}(\mathbb{R}^{|\mathfrak{u}|})$. For $\mathfrak{u}\neq\emptyset$ we obtain from Theorem 2, which holds for $n_{\mathfrak{u}}=2^{m_{\mathfrak{u}}}\geq 2$,

[TABLE]

for any $\lambda\in[1,\alpha)$, and where $\alpha_{1}=\alpha/2+1/4$ and $C_{\alpha,\lambda,|\mathfrak{u}|}$ is as defined in (53). The second bound also holds for $n_{\mathfrak{u}}=0$ and $n_{\mathfrak{u}}=1$, for which in both cases we take the zero approximation, since $M_{\mathfrak{u}}=M^{|\mathfrak{u}|}\leq C_{\alpha,\lambda,|\mathfrak{u}|}$ which follows from $C_{\alpha,\lambda,|\mathfrak{u}|}\geq(2/\sqrt{\ln(2)}+2\sqrt{\lambda})^{(\alpha+1/2)|\mathfrak{u}|}\geq(2/\sqrt{\ln(2)}+2)^{|\mathfrak{u}|}\approx 4.40224^{|\mathfrak{u}|}$ and $M<2.767$, and we used $C_{4,\alpha,\lambda,s}\geq 1$, see Proposition 2 and the constants referenced from (53). Therefore, with $\{\gamma_{\mathfrak{u}}\}_{|\mathfrak{u}|<\infty}$ a sequence of positive weights, the cubature error is bounded as

[TABLE]

where for $\mathfrak{u}=\emptyset$ we interpret $0^{0}$ as $1$. We remark that we deliberately pulled out $|\mathfrak{u}|^{-\alpha_{1}|\mathfrak{u}|}$ to control the logarithmic factor in $n_{\mathfrak{u}}$ later in Theorem 3. This technique was also used in [28] for the MDFEM in the uniform case.

To show that the above formula also holds for $\mathfrak{u}=\emptyset$ we recall from Remark 1 that we have $\|F_{\emptyset}\|_{H_{\alpha,0,\rho,0}}=|F_{\emptyset}|=|\mathcal{F}(\bm{0})|$. Thus, for $n_{\emptyset}=0$ the absolute value of the cubature error is $|\mathcal{F}(\bm{0})|$ while, if for $n_{\emptyset}\geq 1$ we set $Q_{\emptyset,n_{\emptyset}}(F_{\emptyset})=F_{\emptyset}$, the error is zero for $n_{\emptyset}\geq 1$. Hence for any $n\in\mathbb{N}_{0}$ we have

[TABLE]

We have now bounded all three contributions to the error. For the truncation error and the cubature error we still want to choose the weights $\gamma_{\mathfrak{u}}$. For the truncation error we obtain from Lemma 4 that

[TABLE]

with $\bm{b}_{\mathfrak{u}}=\prod_{j\in\mathfrak{u}}b_{j}$ and $C_{\kappa,\alpha}^{\prime}<\infty$, see (58), under the assumptions of $\{b_{j}\}_{j\geq 1}\in\ell^{p^{*}}(\mathbb{N})$ for $p^{*}\in(0,1]$. For the cubature error we obtain, also from Lemma 4, that such a bound holds with $u_{\mathfrak{u}}^{h_{\mathfrak{u}}}$ in place of $u_{\mathfrak{u}}$ and hence we also have

[TABLE]

By choosing $\gamma_{j}=\sqrt{2}\,b_{j}$, and with $\gamma_{\mathfrak{u}}=\prod_{j\in\mathfrak{u}}\gamma_{j}$, we have

[TABLE]

Combining all three errors we obtain the claimed bound for the total error. ∎

Remark 4*.*

In the previous proof we made use of a supremum-norm over all subspaces of the infinite-variate space, namely $\sup_{|\mathfrak{u}|<\infty}\gamma_{\mathfrak{u}}^{-1}\,\|F_{\mathfrak{u}}\|_{H_{\alpha,0,\rho,|\mathfrak{u}|}(\mathbb{R}^{|\mathfrak{u}|})}$. Since we know from Lemma 4 that $\|F_{\mathfrak{u}}\|_{H_{\alpha,0,\rho,|\mathfrak{u}|}(\mathbb{R}^{|\mathfrak{u}|})}\lesssim\bm{b}_{\mathfrak{u}}\,2^{|\mathfrak{u}|/2}$, where $\bm{b}_{\mathfrak{u}}=\prod_{j\in\mathfrak{u}}b_{j}$, we have chosen product weights $\gamma_{\mathfrak{u}}=\prod_{j\in\mathfrak{u}}\gamma_{j}$ with $\gamma_{j}=\sqrt{2}\,b_{j}$ such that this supremum-norm is finite. This same choice of weights $\gamma_{\mathfrak{u}}$ also shows that our infinite-variate function $\mathcal{F}=\sum_{|\mathfrak{u}|<\infty}F_{\mathfrak{u}}$ has finite norm in the infinite-variate reproducing kernel Hilbert space $\mathcal{H}_{\alpha,\rho,\bm{\gamma}}(\mathbb{R}^{\mathbb{N}})$ which we introduced in Remark 1, since

[TABLE]

which is finite when $\sum_{j\geq 1}b_{j}<\infty$ (using the technique as in the proof of Proposition 6). The summability of the $b_{j}$ is implied by our assumption in Proposition 5 which demands $\{b_{j}\}_{j\geq 1}\in\ell^{p^{*}}(\mathbb{N})$ for some $p^{*}\in(0,1]$.

6.3 Selection of the active set

Based on the expression of the truncation error in Proposition 5 we can choose the active set to reach a truncation error upper bounded by $\epsilon/2$ up to multiplicative constants.

Proposition 6**.**

Under the conditions of Proposition 5 with $p^{*}\in(0,1)$, for which $\{b_{j}\}_{j\geq 1}\in\ell^{p^{*}}(\mathbb{N})$, let the MDFEM active set be chosen by

[TABLE]

with $\gamma_{j}=\sqrt{2}\,b_{j}$ for all $j$ as in (68) and with $\gamma_{\mathfrak{u}}=\prod_{j\in\mathfrak{u}}\gamma_{j}$. Then

[TABLE]

and the MDFEM truncation error is bounded as

[TABLE]

Proof.

For the first claim we have the implications, with $\gamma_{\mathfrak{u}}=\prod_{j\in\mathfrak{u}}\gamma_{j}$ and $M_{\mathfrak{u}}=M^{|\mathfrak{u}|}$,

[TABLE]

where we used $\ln(1+x)\leq x$ for $x>-1$ and the last line is true since $\{b_{j}\}_{j\geq 1}\in\ell^{p^{*}}(\mathbb{N})$. For the second claim we have from Proposition 5

[TABLE]

The following result states that both the cardinality of the active set, as well as the cardinalities of each of the individual sets in the active set, increase very slowly with decreasing $\epsilon$. The (upper bound of the) cardinality of the active set also gets smaller for decreasing $p^{*}$.

Proposition 7**.**

Given $\gamma_{\mathfrak{u}}=\prod_{j\in\mathfrak{u}}\gamma_{j}$ with $\{\gamma_{j}\}\in\ell^{p^{*}}(\mathbb{N})$ for some $p^{*}\in(0,1)$ and $\mathfrak{U}(\epsilon)=\mathfrak{U}(\epsilon,p^{*})$ chosen as in (69), then for any $\epsilon>0$ it holds that

[TABLE]

and

[TABLE]

Proof.

See [31, 28]. ∎

6.4 Selection of the cubature rules and FEMs

By the combined error of the FE approximations and the cubature errors from Proposition 5 we can now choose the mesh diameters $h_{\mathfrak{u}}$ for the FE approximations and the number of cubature points $n_{\mathfrak{u}}$ for the cubature formulas to balance the combined error and obtain an upper bound of $\epsilon/2$ up to multiplicative constants using the method of Lagrange multipliers to minimize the associated cost. The maximum appearing in the bound of the next result will be dealt with later in Theorem 3.

Proposition 8**.**

Assume the conditions of Proposition 5 hold with $0<p^{*}\leq(2+{d^{\prime}}/\tau)^{-1}<\frac{1}{2}$ and the cost of solving the linear systems for the FEMs are $O(h_{\mathfrak{u}}^{-{d^{\prime}}})$ as in (65) and the application of $G$ to the FE approximations converge like $O(h_{\mathfrak{u}}^{\tau})$ as in (64). Take

[TABLE]

and take the FE mesh diameters $h_{\mathfrak{u}}$ and the number of cubature points $n_{\mathfrak{u}}$ as the solution to an optimization problem (specified as (74) in the proof of this statement, with solutions (76) and (78) and with $n_{\mathfrak{u}}=2^{\log_{2}(\lfloor k_{\mathfrak{u}}\rfloor)}\in\mathbb{N}_{0}$). Then the combined error from the FE approximations and the cubature approximations of the MDFEM algorithm, using transformed polynomial lattice rules with interlacing factor $\alpha$, is bounded as

[TABLE]

where $\alpha_{1}=\alpha/2+1/4$, and the computational cost is bounded as

[TABLE]

Proof.

From Proposition 5 we obtain

[TABLE]

Note that compared to Proposition 5 we have (possibly) increased the upper bound by multiplying the numerator by $2^{\lambda}$ while only multiplying the denominator by $2^{\lambda}$ for each $\mathfrak{u}$ with $n_{\mathfrak{u}}\in\mathbb{N}_{0}$ when $n_{\mathfrak{u}}\geq 1$. This is such that the optimization problem we formulate next will give us an upper bound for this error. We will deal with the max term later in Theorem 3 and now put an upper bound of $\epsilon/2$ on the sum over $\mathfrak{u}$. We are looking for positive real numbers $k_{\mathfrak{u}}$ and $h_{\mathfrak{u}}$ which are the solutions of the following minimization problem:

[TABLE]

Note that since we set $n_{\mathfrak{u}}=2^{\log_{2}(\lfloor k_{\mathfrak{u}}\rfloor)}\leq k_{\mathfrak{u}}$, the objective function is an upper bound on the cost (66) while the constraint is an upper bound on the sum over $\mathfrak{u}$ in the error bound (71) since $\max\{1,2\,n_{\mathfrak{u}}\}=\max\{1,2^{1+\log_{2}(\lfloor k_{\mathfrak{u}}\rfloor)}\}\geq k_{\mathfrak{u}}$. See also, e.g., [23, Section 4.3] for a similar technique.

Define $a_{\mathfrak{u}}:=\gamma_{\mathfrak{u}}\,2^{\lambda}\,C_{\alpha,\lambda,|\mathfrak{u}|}\,|\mathfrak{u}|^{\alpha_{1}|\mathfrak{u}|}$. The Lagrangian is then given by

[TABLE]

with $\xi$ the Lagrange multiplier. For each $\mathfrak{u}\in\mathfrak{U}(\epsilon)$ we obtain the equations

[TABLE]

From the second equation we obtain (after multiplying with $h_{\mathfrak{u}}$)

[TABLE]

which we insert into the first equation to obtain

[TABLE]

of which we substitute the first form into the constraint to obtain

[TABLE]

From (75) and (76) we also find

[TABLE]

where we have set

[TABLE]

Inserting this expression for $k_{\mathfrak{u}}$ into (77) we obtain the following expression for the Lagrange multiplier

[TABLE]

We require the sum $K_{\epsilon}$ in (79) to be uniformly bounded while $\epsilon\to 0$, that is, we require

[TABLE]

Since $\gamma_{\mathfrak{u}}=2^{|\mathfrak{u}|/2}\prod_{j\in\mathfrak{u}}b_{j}$ with $\{b_{j}\}_{j\geq 1}\in\ell^{p^{*}}(\mathbb{N})$ and, both $\pounds_{\mathfrak{u}}$ and $C_{\alpha,\lambda,|\mathfrak{u}|}$ are at most exponential in $|\mathfrak{u}|$, by applying [28, Lemma 1] the sum (80) is bounded if the following two conditions are satisfied:

[TABLE]

and

[TABLE]

Note that it is required in Theorem 2 that $\lambda\geq 1$, and together with (82) this restricts us to the case when $p^{*}$ is sufficiently small, that is,

[TABLE]

Making use of (75) and the first expression in (76), and then (77) and (79), the computational cost (66) is bounded like

[TABLE]

Using the stated conditions we know that $K_{\epsilon}\leq K<\infty$. Hence we can write

[TABLE]

To optimize the speed of convergence we have to pick $\lambda$ and $\tau$ as large as possible while satisfying (82). Hence we choose

[TABLE]

and we have $\lambda\geq 1$ due to (83). By choosing the interlacing order of the interlaced polynomial lattice rule to be $\alpha=\lfloor\lambda\rfloor+1$ we satisfy $\alpha\geq 2$ and $\lambda\in[1,\alpha)$, and we also satisfy (81), since

[TABLE]

for any $\lambda\geq 0$. This finishes the proof. ∎

Remark 5*.*

Since Proposition 8 is making use of interlaced polynomial lattice rules to achieve higher order convergence, i.e., $\lambda\geq 1$, we end up with the condition $0<p^{*}\leq(2+{d^{\prime}}/\tau)^{-1}$ there. Following the proof we see that if we make use of first order cubature rules with error bounds for $\tfrac{1}{2}\leq\lambda<1$, we obtain the condition $p^{*}\leq(3/2+{d^{\prime}}/(2\tau))^{-1}$ by using (82). In particular, we can use transformed randomly digitally shifted polynomial lattice rules from Remark 2 as cubature rules in the MDFEM algorithm. A full description and analysis of using such randomized cubature rules in the context of the MDFEM algorithm for the uniform case is given in [28]. We will write $\mathbb{E}^{\Delta(\epsilon)}[|\mathcal{I}(G(u))-Q_{\epsilon}(G(u))|^{2}]$ to denote the total mean square error of the randomized MDFEM algorithm using randomly digitally shifted polynomial lattice rules. The expected value is taken over a set $\Delta(\epsilon):=\{\bm{\Delta}^{(\mathfrak{u})}\}_{\mathfrak{u}\in\mathfrak{U}(\epsilon)}$ of independent random digital shifts, one for each subproblem $\mathfrak{u}$, and where each such shift $\bm{\Delta}^{(\mathfrak{u})}$ is uniformly distributed over $[0,1)^{|\mathfrak{u}|}$.

6.5 Main result

We are now able to analyze the complexity of the MDFEM algorithm. The analysis is under the conditions of Propositions 5–8 extended to the randomized setting using Remarks 2 and 5.

Theorem 3**.**

Assume that for a given $\alpha\in\mathbb{N}$ there is a sequence $\{b_{j}\}_{j\geq 1}\in\ell^{p^{*}}(\mathbb{N})$ for some $p^{*}\in(0,1)$, with $0<b_{j}\leq 1$ and

[TABLE]

Let $G\in V^{*}$ be a bounded linear functional. Assume that $D$, $a$, $f$ and $G$ have sufficient spatial regularity such that the application of $G$ to the FE approximations of the truncated solutions converge like $O(h_{\mathfrak{u}}^{\tau})$ as in (64) and assume that the cost of solving the linear systems for the FEMs are $O(h_{\mathfrak{u}}^{-{d^{\prime}}})$ as in (65).

For a given requested error $\epsilon>0$ take the active set $\mathfrak{U}(\epsilon)=\mathfrak{U}(\epsilon,p^{*})$ as in (69). Set

[TABLE]

Further, take the FE mesh diameters $h_{\mathfrak{u}}$ and the number of cubature points $n_{\mathfrak{u}}$ as the solution to the optimization problem (74) (see Proposition 8 and Remark 5, with solutions (76) and (78) and with $n_{\mathfrak{u}}=2^{\log_{2}(\lfloor k_{\mathfrak{u}}\rfloor)}$), such that the convergence of the cubature rules over the unit cube can be bounded by $O(n_{\mathfrak{u}}^{-\lambda})$, see Theorems 1 and 2 and Remark 2.

Then, for the MDFEM algorithm $Q_{\epsilon}$, given in (24), the following hold.

• •

If $0<p^{*}\leq(2+{d^{\prime}}/\tau)^{-1}$, i.e., $\lambda\geq 1$ and $\alpha\geq 2$, then, using transformed interlaced polynomial lattice rules with interlacing factor $\alpha$, we have

[TABLE]

• •

If $(2+{d^{\prime}}/\tau)^{-1}<p^{*}\leq(3/2+{d^{\prime}}/(2\tau))^{-1}$, i.e., $\frac{1}{2}\leq\lambda<1$ and $\alpha=1$, then, using transformed randomly digitally shifted polynomial lattice rules, we have

[TABLE]

In both cases the computational cost is bounded as

[TABLE]

Proof.

We first show the statement for the case $0<p^{*}\leq(2+{d^{\prime}}/\tau)^{-1}$. It follows from Propositions 5, 6 and 8 that the error of the MDFEM is then bounded by

[TABLE]

with computational cost

[TABLE]

Using (70) from Proposition 7 and $\ln(n_{\mathfrak{u}})\lesssim\ln(\epsilon^{-1/\lambda})\lesssim\ln(\epsilon^{-1})$, see (78) in combination with (79), we can use the same argument as in [28, Theorem 1], see also [31, Lemma 1], and we have

[TABLE]

where $\delta(\epsilon)=O\left(\ln(\ln(\ln(\epsilon^{-1})))/\ln(\ln(\epsilon^{-1}))\right)=o(1)$ as $\epsilon\to 0$. Hence, we can write

[TABLE]

For the case $(2+{d^{\prime}}/\tau)^{-1}<p^{*}\leq(3/2+{d^{\prime}}/(2\tau))^{-1}$ the statement follows using similar arguments, a full proof is provided in [28]. ∎

Finally we would like to compare the complexity of the MDFEM with the quasi-Monte Carlo finite element method (QMCFEM) [18] and with the multilevel quasi-Monte Carlo finite element method (MLQMCFEM) [17], both in the setting of using product weights and a wavelet expansion for the lognormal field.

The QMCFEM truncates the parameter vector ${\bm{y}}$ to some dimension $s$ and then approximates the $s$-truncated PDE for different samples ${\bm{y}}^{(i)}\in\mathbb{R}^{s}$, $i=0,\ldots,N-1$, obtained by a QMC method. Because $s$ might be arbitrarily large, the QMCFEM requires QMC rules with convergence independent of the dimension of the integrand. Such QMC rules over the Euclidean space $\mathbb{R}^{s}$ with the Gaussian distribution were developed in [24] by mapping randomly shifted lattice rules over the unit cube $[0,1]^{s}$ to $\mathbb{R}^{s}$ by the inverse of the normal cumulative distribution. However, since this mapping might damage the smoothness of the integrand, the convergence rate of these QMC rules was limited to first order (with respect to the number of QMC points $N$). Particularly, the QMCFEM was analysed in [18, Section 10] under the same conditions as those of Theorem 3 (marked with the subscripts “PROD” for product weights and “wav” for wavelet expansion, and making use of “Gaussian weight functions” for the QMC rules) to achieve a root-mean-square error bound of the form

[TABLE]

where $\delta>0$ is a parameter that might be chosen arbitrarily small (but then increases the hidden constant towards infinity), $N$ is the number of cubature points, $h$ is the FE mesh diameter and $s$ is the truncation dimension. (Note that [18, Section 10] considers the specific choice of a Gaussian random field with Matérn covariance for which we (optimistically) set $p^{*}=d/\nu$, with $\nu$ the smoothness parameter of the Matérn covariance function. Note additionally that in [3] it is proven that $p^{*}$ can be chosen arbitrarily close to $d/\nu$ in the case of a wavelet-type expansion.) The cost for the QMCFEM can be estimated by $\mathrm{cost}(Q^{\mathrm{QMCFEM}})\lesssim N(s+h^{-{d^{\prime}}})$. Since we are interested in the asymptotic best possible rates, we will ignore the hidden constants and focus on the rate. This means we formally set $\delta=0$ in (85) as a proxy of being arbitrarily close to the optimal rate and thereby obtaining a slightly optimistic bound on the work. By balancing all three contributions to be $\epsilon/3$ we find

[TABLE]

A smaller exponent means less work. For both exponents, $a_{\mathrm{QMCFEM}}$ as just defined, and $a_{\mathrm{MDFEM}}$ as given in (84), we can recognize two parts: the first part is the convergence order of the QMC cubature w.r.t. its cost and the second part is the convergence order of the linear functional applied to the FE approximation w.r.t. the cost of the FE approximation.

We now look at the case $0<p^{*}\leq(2+{d^{\prime}}/\tau)^{-1}<1/2<2/3$ when the MDFEM can employ a higher-order QMC cubature rule. We then have $\lambda\geq 1$ and hence $a_{\mathrm{MDFEM}}$ has the factor $1/\lambda$ for the QMC cubature part. For the QMCFEM the convergence of the QMC cubature is limited to $1$ and hence, for the case $0<p^{*}\leq 2/3$ we obtain

[TABLE]

It is easy to see that if ${d^{\prime}}/\tau\geq 1/4$ then we have $0<p^{*}\leq(2+{d^{\prime}}/\tau)^{-1}\leq(4{d^{\prime}}/\tau)(2+{d^{\prime}}/\tau)^{-1}$. Hence $a_{\mathrm{MDFEM}}=1/\lambda+{d^{\prime}}/\tau\leq a_{\mathrm{QMCFEM}}=1+{d^{\prime}}/\tau$ and the MDFEM and the QMCFEM perform the same when $\lambda=1$ and the MDFEM will perform better when $\lambda>1$. On the other hand, if ${d^{\prime}}/\tau<1/4$ then there are two cases to be considered. When $0<p^{*}\leq(4{d^{\prime}}/\tau)(2+{d^{\prime}}/\tau)^{-1}<(2+{d^{\prime}}/\tau)^{-1}$ we have exactly the same conclusion as for ${d^{\prime}}/\tau\geq 1/4$. Lastly, when $0<(4{d^{\prime}}/\tau)(2+{d^{\prime}}/\tau)^{-1}<p^{*}\leq(2+{d^{\prime}}/\tau)^{-1}$ the QMCFEM is in the unfortunate case that the term for the FEM convergence is larger than ${d^{\prime}}/\tau$ and so it will always lose.

Hence, for any $p^{*}$ such that $0<p^{*}\leq(2+{d^{\prime}}/\tau)^{-1}$ the MDFEM performs better or similar compared to the QMCFEM, with the possibility of higher-order convergence for the MDFEM.

As a second comparison we look at the multilevel variant of the QMCFEM which is the MLQMCFEM algorithm given in [17]. Since the QMC rules for the MLQMCFEM are the same kind as those for the QMCFEM, they can achieve at most order $1$ convergence and this requires $0<p^{*}<2/3$. The error bounds in [17, Theorem 6.3 and 6.5] take more complicated forms in which $\bar{\xi}$ there is the convergence of the QMC rule on each level, comparable with our $\lambda$, and $\eta$ there is our ${\delta^{\prime}}$. They are written in terms of the dominating cost, being either $\epsilon^{-1/\lambda}$, with $\lambda<1$ due to the QMC rules there, or $\epsilon^{-{d^{\prime}}/\tau}$ with extra log factors depending on the situation. So also in this case it is clear that in the case $0<p^{*}\leq(2+{d^{\prime}}/\tau)^{-1}<1/2<2/3$ if we take $\lambda\geq 1$ for the MDFEM and we can take ${d^{\prime}}/\tau\leq 1/\lambda$ then the MDFEM performs better or similar compared to the MLQMCFEM, with the possibility of higher-order convergence for the MDFEM.

Acknowledgements

The authors would like to thank the three anonymous referees for very careful and detailed commenting on the manuscript. The authors would also like to thank Frances Y. Kuo and Kien V. Nguyen for discussion on parts of the manuscript. We gratefully acknowledge the financial support from the Research Foundation Flanders (FWO) under grant G091920N.

Appendix A Appendix

We collect some results here which would otherwise disturb the flow of the paper.

A.1 Proof of Lemma 1

Proof of Lemma 1.

The first property, (15), is a well known property of the anchored decomposition and we omit the proof. We now prove the second property, (16), which is,

[TABLE]

i.e., for ${\bm{\omega}}_{\mathfrak{u}}\in\mathbb{N}^{|\mathfrak{u}|}$. First we move the derivative operator inside the explicit form for $F_{\mathfrak{u}}$, see (13),

[TABLE]

and realize that $\mathcal{F}(\cdot_{\mathfrak{v}})$ is a function which is constant in all variables not in $\mathfrak{v}$, and as such if we have an $\omega_{j}\geq 1$ for which $j\notin\mathfrak{v}$ then the derivative $(\partial^{{\bm{\omega}}_{\mathfrak{u}}}_{{\bm{y}}_{\mathfrak{u}}}\mathcal{F}(\cdot_{\mathfrak{v}}))$ is zero. Since ${\bm{\omega}}_{\mathfrak{u}}\in\mathbb{N}^{|\mathfrak{u}|}$ we have $\omega_{j}\geq 1$ for all $j\in\mathfrak{u}$ and thus the only remaining term for $\mathfrak{v}\subseteq\mathfrak{u}$ is the one for which $\mathfrak{v}=\mathfrak{u}$.

To show the third property, (17), which is

[TABLE]

we first note that here ${\bm{\omega}}_{\mathfrak{u}}\in\mathbb{N}_{0}^{|\mathfrak{u}|}$, i.e., $\omega_{j}$ is allowed to be zero for $j\in\mathfrak{u}$. Therefore, suppose there is a $j\in\mathfrak{u}$ such that $\omega_{j}=0$, then we are actually not taking the partial derivative w.r.t. the $j$th variable, and by the definition of the partial derivative this means we keep the $j$th variable fixed while we take the derivatives w.r.t. the variables in $\mathop{\mathrm{supp}}({\bm{\omega}}_{\mathfrak{u}})$. Then, in the case that $\omega_{j}=0$ we can do the evaluation at $y^{*}_{j}$ before taking the other partial derivatives, i.e., $(\partial^{{\bm{\omega}}_{\mathfrak{u}}}_{{\bm{y}}_{\mathfrak{u}}}F_{\mathfrak{u}})({\bm{y}}^{*}_{\mathfrak{u}})=(\partial^{{\bm{\omega}}_{\mathfrak{u}}}_{{\bm{y}}_{\mathfrak{u}}}(F_{\mathfrak{u}}|_{y_{j}=y^{*}_{j}}))({\bm{y}}^{*}_{\mathfrak{u}\setminus\{j\}})$. Since $j\in\mathfrak{u}$ we can use the first property (15) to deduce that $F_{\mathfrak{u}}|_{y_{j}=y^{*}_{j}}=0$ when $y^{*}_{j}=0$. ∎

A.2 Derivation of anchored Sobolev kernels and Taylor representation

We provide the derivation of the reproducing kernel of a univariate anchored Sobolev space w.r.t. a positive weight function $\rho$ for anchored functions. In particular this shows how to obtain the kernel given in (33) for the anchored Gaussian Sobolev space for anchored functions. For $\alpha\in\mathbb{N}$ the space consists of functions that have absolutely continuous derivatives up to order $\alpha-1$ for any bounded interval and have square integrable derivative of order $\alpha$ w.r.t. the weight function $\rho$ and are anchored at [math]. The domain is the support of $\rho$. These spaces are meant to be used for functions $F_{\mathfrak{u}}$ obtained by the anchored decomposition (13).

In accordance with [26, Example 4.4] and [12], the only constant function in our univariate anchored space is the zero function, and we can use the tensor product of this univariate kernel to create a multivariate kernel which represents multivariate anchored functions. A reproducing kernel $K$ for a reproducing kernel Hilbert space $\mathcal{H}(K)$ has the property that $K(x,y)$ as a function of $y$ is a function in $\mathcal{H}(K)$ for any $x$ in the domain. To satisfy this property, the kernel for the anchored Sobolev space on the unit cube for anchored functions in [26, Example 4.2] and [6, Section 5.2], should be amended to have the sum over the derivatives from $1$ to $\alpha-1$ where the “[math] otherwise” appears. Then it agrees with the kernel given here.

Our proof uses similar arguments as in [32, Section 1.2] where the kernel of the anchored Sobolev space over the unit cube is derived, but we modify the techniques for the case when the inner product contains a positive weight function $\rho$ in the $L^{2}$ inner product of the $\alpha$th derivatives, and when the functions are anchored.

Proposition 9**.**

The reproducing kernel of the anchored Sobolev space for anchored functions $H_{\alpha,0,\rho}(\mathbb{R})$, with $\rho(t)=\rho(-t)$, e.g., the anchored Gaussian Sobolev space, with inner product

[TABLE]

is given by

[TABLE]

where $\mathds{1}\{X\}$ is the indicator function on $X$.

Proof.

By Taylor’s theorem, and using $F(0)=0$, we have

[TABLE]

where we have written both forms as this helps to interpret how our kernel will operate in connection to the integral over $\mathbb{R}$ in the inner product. By the reproducing property we require

[TABLE]

where the derivatives of $K_{\alpha,0,\rho}$ are taken with respect to the first variable.

Comparing the two representations of $F$ leads us to choose the kernel $K_{\alpha,0,\rho}$ such that

[TABLE]

and

[TABLE]

Since $K_{\alpha,0,\rho}(\cdot,y)$ itself needs to be a function from the space $H_{\alpha,0,\rho}(\mathbb{R})$ we require that $K_{\alpha,0,\rho}(0,y)=0$ for any $y$ and its Taylor expansion with respect to the first variable is given by

[TABLE]

Therefore, inserting (86) and (87) into (88) we find for $y\geq 0$ that

[TABLE]

which can be written as

[TABLE]

Similarly, for $y\leq 0$ we have

[TABLE]

Hence, under the assumption that $\rho(t)=\rho(-t)$, the explicit formula for $K_{\alpha,0,\rho}$ can be written as

[TABLE]

We have stated Proposition 9 in a form which is relevant to the paper, but in fact the proof is more general. As stated, we do not require $\rho$ to be the Gaussian density. The condition $\rho(t)=\rho(-t)$ in the proposition is only stated to obtain an easy final expression and could be removed. The resulting kernel is valid for other positive weight functions $\rho$ and we assume the domain of the integrals are then truncated to the support of $\rho$. E.g., $\rho$ could be the uniform density on the unit cube and we then recover the anchored Sobolev space for anchored functions on $[0,1]$. If we want to drop the requirement that $F(0)=0$ then we add the constant one to the kernel and include the term for $\tau=0$ in the inner product. Such a kernel for $\rho\equiv 1$ over $\mathbb{R}$ is given in [30, Section 11.5.1]. The statements are also easily generalizable for an arbitrary anchor point, similar as the kernels in [26, Example 4.2] and [6, Section 5.2].

We take the multivariate space to be the tensor product of the univariate spaces and therefore the kernel is obtained as the product, see, e.g., [26, Example 4.4]. For anchored functions in such a multivariate anchored Sobolev space (not necessarily just the Gaussian case as we use in this paper) we have the following representation as a Taylor series with integral remainder term.

Lemma 5**.**

For $\mathfrak{u}\subset\mathbb{N}$ with $1\leq|\mathfrak{u}|<\infty$, assume $F_{\mathfrak{u}}\in H_{\alpha,0,\rho,|\mathfrak{u}|}(\Omega^{|\mathfrak{u}|})$, with $\Omega$ the support of the positive weight function $\rho$, is obtained from an anchored decomposition (13) of a function $\mathcal{F}$. Then we have the representation

[TABLE]

where $F_{\mathfrak{u}}^{(\bm{\nu}_{\mathfrak{u}})}(\bm{t}_{\mathfrak{v}})$ means to evaluate the function $F_{\mathfrak{u}}^{(\bm{\nu}_{\mathfrak{u}})}$ in the point $\bm{t}_{\mathfrak{v}}\in\mathbb{R}^{\mathbb{N}}_{\mathfrak{v}}$, i.e., setting all other arguments to zero, or, if viewing $F_{\mathfrak{u}}^{(\bm{\nu}_{\mathfrak{u}})}$ as a $|\mathfrak{u}|$-variate function then $F_{\mathfrak{u}}^{(\bm{\nu}_{\mathfrak{u}})}(\bm{t}_{\mathfrak{v}})=F_{\mathfrak{u}}^{(\bm{\nu}_{\mathfrak{u}})}(\bm{t}_{\mathfrak{v}},\bm{0}_{\mathfrak{u}\setminus\mathfrak{v}})$. For ${\bm{\omega}}_{\mathfrak{u}}\in\{0:\alpha\}^{|\mathfrak{u}|}$ we have

[TABLE]

The functions $F_{\mathfrak{u}}^{(\bm{\nu}_{\mathfrak{u}})}$ under the integrals could be replaced with $\mathcal{F}^{(\bm{\nu}_{\mathfrak{u}})}$ in both (89) and (90).

Proof.

Functions in $H_{\alpha,0,\rho,|\mathfrak{u}|}(\Omega^{|\mathfrak{u}|})$ satisfy all requirements to use the Taylor theorem with integral remainder up to order $\alpha$ in each direction successively, i.e., they have absolutely continuous derivatives up to order $\alpha-1$ for any bounded interval. Hence

[TABLE]

where we remark that we wrote $\bm{\nu}_{\mathfrak{u}}\in\{0:\alpha\}^{|\mathfrak{u}|}$, i.e., $\nu_{j}$ is allowed to be zero also for $j\in\mathfrak{u}$. We can now make use of Lemma 1. In particular, property (17) implies that any term for which there is at least one $j\in\mathfrak{u}$ for which $\nu_{j}=0$ will vanish. This proves (89).

To obtain the expression for $F_{\mathfrak{u}}^{({\bm{\omega}}_{\mathfrak{u}})}({\bm{y}}_{\mathfrak{u}})$ for ${\bm{\omega}}_{\mathfrak{u}}\in\{0:\alpha\}^{|\mathfrak{u}|}$ we work in a similar way by applying the Taylor theorem with integral remainder term to $F_{\mathfrak{u}}^{({\bm{\omega}}_{\mathfrak{u}})}$ for each $\omega_{j}<\alpha$. For $\omega_{j}=\alpha$ we cannot apply the Taylor theorem anymore, so we need to take care that in our expression we can just recover the $\alpha$th derivatives, i.e., we should not integrate those components for which $\omega_{j}=\alpha$. Therefore we introduce the set $\mathfrak{w}:=\{j\in\mathfrak{u}:\omega_{j}=\alpha\}$. As in the expression for $F_{\mathfrak{u}}$ we gather the indices for which we need the integral in the set $\mathfrak{v}$ which is now modified to exclude those indices in the set $\mathfrak{w}$. This proves (90).

By property (16) we also know that $F_{\mathfrak{u}}^{(\bm{\nu}_{\mathfrak{u}})}({\bm{y}}_{\mathfrak{u}})=\mathcal{F}^{(\bm{\nu}_{\mathfrak{u}})}({\bm{y}}_{\mathfrak{u}})$ for $\bm{\nu}_{\mathfrak{u}}\in\{1:\alpha\}^{|\mathfrak{u}|}$ and any ${\bm{y}}_{\mathfrak{u}}$ so also for ${\bm{y}}_{\mathfrak{u}}=(\bm{t}_{\mathfrak{v}},\bm{0}_{\mathfrak{u}\setminus\mathfrak{v}})\in\mathbb{R}^{\mathbb{N}}_{\mathfrak{u}}$ or ${\bm{y}}_{\mathfrak{u}}=(\bm{t}_{\mathfrak{v}},{\bm{y}}_{\mathfrak{w}},\bm{0}_{\mathfrak{z}})\in\mathbb{R}^{\mathbb{N}}_{\mathfrak{u}}$ as in (89) and (90) respectively. ∎

A.3 Proof of Proposition 3: norm embedding after mapping

We deferred the proof of Proposition 3 due to its length. We first show two lemmas that we need for the proof. The following result is taken from [7, Lemma 3].

Lemma 6**.**

For $F:\mathbb{R}^{s}\to\mathbb{R}$ having mixed partial derivatives for all $\bm{\tau}\in\{0:\alpha\}^{s}$ we have for any $\bm{\tau}\in\{0:\alpha\}^{s}$ that

[TABLE]

where the sum is over all ${\bm{\omega}}\in\mathbb{N}_{0}^{s}$ for which $0\leq\omega_{j}\leq\tau_{j}$ for all $j\in\{1,2,\ldots,s\}$, $c(\bm{\tau},{\bm{\omega}}):={\bm{\tau}\choose\bm{\tau}-{\bm{\omega}}}\sqrt{(\bm{\tau}-{\bm{\omega}})!}:=\prod_{j=1}^{s}{\tau_{j}\choose\tau_{j}-\omega_{j}}\sqrt{(\tau_{j}-\omega_{j})!}$ and $H_{\bm{\tau}}({\bm{y}}):=\prod_{j=1}^{s}H_{\tau_{j}}(y_{j})$, with $H_{\tau}$, for $\tau\in\mathbb{N}_{0}$, the $\tau$th normalized probabilistic Hermite polynomial given by

[TABLE]

Lemma 7**.**

For $\alpha\in\mathbb{N}$, $\tau\in\{0:\alpha\}$ and $\eta\in[0,2\alpha-1]$ we have

[TABLE]

with $I_{0}$ the modified Bessel function of the first kind of order [math].

Proof.

Using (91), the Cauchy–Schwarz inequality and $\int_{\mathbb{R}}|y|^{\eta}\rho(y)\,\mathrm{d}y=2^{\eta/2}\,\Gamma((\eta+1)/2)/\sqrt{\pi}$ for any $\eta\geq 0$, we have

[TABLE]

where in the second inequality we use $\tau\leq\alpha$, $\eta\leq 2\alpha-1$, $\tau-2k\geq 0$, $1/2\leq\tau-2k+(\eta+1)/2\leq 2\alpha$ for any $k\in\{0:\lfloor\tau/2\rfloor\}$ and $\Gamma(1/2)\leq\Gamma(2\alpha)$ for $\alpha\geq 1.43258$. For $\alpha=1$ the result follows by directly comparing the stated upper bound with the closed form solution of the integral for $\tau=0$ (given above) and $\tau=1$ (which is the expression above with $\eta$ replaced by $\eta+2$). ∎

We are now ready to give the proof of Proposition 3 which states that for $F\in H_{\alpha,0,\rho,s}(\mathbb{R}^{s})$ with $\alpha\in\mathbb{N}$ and $T\geq 1/(2\sqrt{2})$, the function $(F\rho)\circ\bm{T}:[0,1]^{s}\to\mathbb{R}^{s}$ belongs to $H_{\alpha,s}([0,1]^{s})$ and

[TABLE]

with $C_{1,\alpha}$ defined in (43). The proof starts along the lines of [7, proof of Lemma 4] but we need some extra work to arrive at the norm in the space $H_{\alpha,0,\rho,s}(\mathbb{R}^{s})$ for which we will make use of Lemma 7. We remind the reader that $F(y_{1},\ldots,y_{s})$ is actually some relabelling of $F_{\mathfrak{u}}({\bm{y}}_{\mathfrak{u}})$ with $|\mathfrak{u}|=s$ coming from an anchored decomposition and that $H_{\alpha,0,\rho,s}(\mathbb{R}^{s})$ is the anchored Gaussian Sobolev space specifically for such functions having the properties listed in Lemma 1.

Proof of Proposition 3.

Using the inner product (41), the Cauchy–Schwarz inequality and the chain rule $((F\rho)\circ T)^{(\bm{\tau})}({\bm{y}})=(2T)^{|\bm{\tau}|}(F\rho)^{(\bm{\tau})}({\bm{y}})$ we have

[TABLE]

Applying Lemma 6 to the last inequality we then write

[TABLE]

Applying the Cauchy–Schwarz inequality to the last integral leads to

[TABLE]

We will show below that for any $\bm{\tau}\in\{0:\alpha\}^{s}$ and ${\bm{\omega}}\in\mathbb{N}_{0}^{s}$ such that ${\bm{\omega}}\leq\bm{\tau}$ we have the uniform bound

[TABLE]

where $C_{*,\alpha}:=\alpha\,C_{\diamond,\alpha}$ with $C_{\diamond,\alpha}$ defined in (92). Inserting (95) into (A.3) then leads to

[TABLE]

Moreover, we have, for $\bm{\tau}\in\{0:\alpha\}^{s}$ and ${\bm{\omega}}\leq\bm{\tau}$,

[TABLE]

Inserting this into (96) implies

[TABLE]

where we used $(4T)^{2\alpha+2}/((4T)^{2}-1)\leq 2\,(4T)^{2\alpha}$ for $T\geq 1/(2\sqrt{2})>1/4$, and where

[TABLE]

We have now arrived at the claim of our statemement (93).

To complete the proof we still need to show (95). We are going to use the Taylor representation from Lemma 5 for the derivatives of $F=F_{\mathfrak{u}}$, where in our current exposition $\{1:s\}$ is a relabeling of $\mathfrak{u}$. Thus, using (90), and with the understanding that $F=F_{\mathfrak{u}}$, ${\bm{\omega}}={\bm{\omega}}_{\mathfrak{u}}\in\{0:\alpha\}^{|\mathfrak{u}|}$ and ${\bm{y}}={\bm{y}}_{\mathfrak{u}}\in\mathbb{R}^{|\mathfrak{u}|}$, we have, with some slight abuse of notation,

[TABLE]

Note that inside the integral we evaluate the function $F=F_{\mathfrak{u}}$ at $(\bm{t}_{\mathfrak{v}},{\bm{y}}_{\mathfrak{w}})=(\bm{t}_{\mathfrak{v}},{\bm{y}}_{\mathfrak{w}},\bm{0}_{\mathfrak{z}})\in\mathbb{R}^{|\mathfrak{u}|}$ with $\mathfrak{z}=\mathfrak{u}\setminus(\mathfrak{v}\cup\mathfrak{w})$ and $\mathfrak{v},\mathfrak{w},\mathfrak{z}$ are pairwise disjoint with $\mathfrak{u}=\mathfrak{v}\cup\mathfrak{w}\cup\mathfrak{z}$ since the above definitions are equivalent to

[TABLE]

where we used that $\omega_{j}\leq\nu_{j}\leq\alpha$ for the set $\mathfrak{w}$.

In what follows, for any $y\in\mathbb{R}$ we write $[0,y]^{*}=[0,y]$ if $y\geq 0$ and $[0,y]^{*}=[y,0]$ if $y<0$. For ${\bm{y}}\in\mathbb{R}^{s}$ we write $[\bm{0},{\bm{y}}]^{*}=[0,y_{1}]^{*}\times\cdots\times[0,y_{s}]^{*}$. So, squaring the above expression for $F^{({\bm{\omega}})}$ and applying the Cauchy–Schwarz inequality twice we obtain

[TABLE]

Note that in the second and third step we used that for $j\in\mathfrak{v}$ we know that $\omega_{j}\leq\alpha-1$, and in the last step we used that for $j\in\mathfrak{z}$ we have $\nu_{j}\leq\alpha-1$. We want to use this bound for $|F^{({\bm{\omega}})}({\bm{y}})|^{2}$ in $\mathcal{Y}^{2}(\bm{\tau},{\bm{\omega}})$, cf. (95), which is multiplying with $|H_{\bm{\tau}-{\bm{\omega}}}({\bm{y}})|^{2}\,\rho^{2}({\bm{y}})$ and taking the integral over $[-T,T]^{s}$. We can move the integral inside of the sum over $\bm{\nu}$ in the above expression to obtain

[TABLE]

with $B(\bm{\nu},{\bm{\omega}},{\bm{y}})$ defined in (97). Then

[TABLE]

We now split the integral over ${\bm{y}}$ into a product of three integrals over the pairwise disjoint sets $\mathfrak{v}\cup\mathfrak{w}\cup\mathfrak{z}=\{1:s\}$. We obtain

[TABLE]

where we used Lemma 7 and the constant $C_{\diamond,\alpha}$ defined there in (92), and $H_{\bm{\tau}}({\bm{y}})\sqrt{\rho({\bm{y}})}\leq 1$ for any $\bm{\tau}\in\mathbb{N}_{0}^{s}$ and any ${\bm{y}}\in\mathbb{R}^{s}$, see [7, Lemma 1]. Inserting (99) into (98) leads to

[TABLE]

where in the last line we obtain the norm based on the inner product (3.1.1) of the space $H_{\alpha,0,\rho,s}(\mathbb{R}^{s})$. This shows (95). The proof is now complete. ∎

A.4 Interlaced polynomial lattice rules for $H_{\alpha,s}([0,1]^{s})$

The aim of this section is to show that interlaced polynomial lattice rules can achieve the almost optimal order of convergence for integration in the space $H_{\alpha,s}([0,1]^{s})$ defined in Section 3.1.2. Interlaced polynomial lattice rules were also used in the setting of PDEs with random diffusion coefficient, but for the uniform case, i.e., with integrals directly expressible over the unit cube, in [8]. Here we map our integrals over the full space into the unit cube by the strategy described in Section 3. But the unanchored Sobolev space here is different from the unanchored Sobolev space considered in [8]. We adjust the analysis of [8] and [13] to show that the fast component-by-component construction algorithm as in [8] can also construct optimal interlaced polynomial lattice rules for our space $H_{\alpha,s}([0,1]^{s})$. We remind the reader that the inner product of our space was already given in (40) and (41). As explained in the introduction, we do not consider weighted function spaces, since the MDM already takes care to limit the number of dimensions for each subproblem.

Interlaced polynomial lattice rules are a modification of polynomial lattice rules to achieve higher-order convergence for integration over the unit cube in classes of Walsh spaces and weighted unanchored Sobolev spaces, see, e.g., [8, 13]. The aim is to approximate multivariate integrals over the $s$-dimensional unit cube

[TABLE]

by a quasi-Monte Carlo rule of the form

[TABLE]

where $P_{n}:=\{{\bm{y}}^{(i)}\}_{i=0}^{n-1}$ is the cubature point set. The worst-case error of the QMC rule $Q_{[0,1]^{s},P_{n}}$ in the normed space $H_{\alpha,s}([0,1]^{s})$ is defined by

[TABLE]

Hence, for any $F\in H_{\alpha,s}([0,1]^{s})$

[TABLE]

We need to introduce some necessary definitions. For simplicity we restrict ourselves to polynomial lattice rules over the finite field $\mathbb{Z}_{2}$. Let $\mathbb{Z}_{2}[\chi]$ denote the set of all polynomials over $\mathbb{Z}_{2}$ and $\mathbb{Z}_{2}[\chi^{-1}]$ denote the set of all formal Laurent series over $\mathbb{Z}_{2}$. For any $m\in\mathbb{N}$ let us define a mapping $\vartheta_{m}:\mathbb{Z}_{2}[\chi^{-1}]\to[0,1)$ by

[TABLE]

In the following we will identify any integer $k\in\{0,\ldots,2^{m}-1\}$, having binary expansion $k=\kappa_{0}+\kappa_{1}2+\cdots+\kappa_{m-1}2^{m-1}$, with the polynomial $k(\chi)=\kappa_{0}+\kappa_{1}\chi+\cdots+\kappa_{m-1}\chi^{m-1}\in\mathbb{Z}_{2}[\chi]$ and vice versa.

Definition 1** (polynomial lattice rule).**

For $m,s\in\mathbb{N}$ let $p\in\mathbb{Z}_{2}[\chi]$ be an irreducible polynomial such that $\deg(p)=m$ and let ${\bm{q}}=(q_{1},\ldots,q_{s})\in\mathscr{G}_{m}^{s}$ with

[TABLE]

A polynomial lattice point set $P_{p,m,s}({\bm{q}})$ is a set of $n=2^{m}$ points ${\bm{y}}^{(0)},\ldots,{\bm{y}}^{(2^{m}-1)}\in[0,1)^{s}$ where

[TABLE]

A QMC rule using this point set is called a polynomial lattice rule with generating vector ${\bm{q}}$ and modulus $p$.

The convergence for a polynomial lattice rule is typically that of a normal QMC rule, i.e., $O(n^{-1+\delta})$, $\delta>0$, under appropriate conditions and modulo $\log$-factors. By making use of interlacing we can obtain higher order convergence $O(n^{-\alpha+\delta})$, $\delta>0$, for $\alpha>1$, again, under appropriate conditions and modulo $\log$-factors. Interlacing is the process of combining $\alpha\in\mathbb{N}$, $\alpha\geq 2$, base-$b$ elements into one. I.e., we can interlace a tuple of $\alpha$ base-$2$ numbers in the interval $[0,1)$ and combine them into one base-$2$ number in the interval $[0,1)$. This is the aim of the digit interlacing function which will take a polynomial lattice point set in $\alpha s$ dimensions and interlace the points with a factor $\alpha$ to obtain an interlaced polynomial lattice point set in $s$ dimensions.

Definition 2** (interlaced polynomial lattice rule).**

Define the digit interlacing function $\mathscr{D}_{\alpha}:[0,1)^{\alpha}\to[0,1)$ with interlacing factor $\alpha\in\mathbb{N}$ by

[TABLE]

where $y_{j}=\xi_{1,j}2^{-1}+\xi_{2,j}2^{-2}+\cdots$ for $j=1,\ldots,\alpha$ and, in case the number of arguments is a multiple of $\alpha$, define $\mathscr{D}_{\alpha}:[0,1)^{\alpha s}\to[0,1)^{s}$ by

[TABLE]

For $m,s\in\mathbb{N}$ let $p\in\mathbb{Z}_{2}[\chi]$ be an irreducible polynomial such that $\deg(p)=m$ and let ${\bm{q}}=(q_{1},\ldots,q_{\alpha s})\in\mathscr{G}_{m}^{\alpha s}$. An interlaced polynomial lattice point set (of order $\alpha$) $\mathscr{D}_{\alpha}(P_{p,m,\alpha s}({\bm{q}}))$ is a set of $n=2^{m}$ points ${\bm{y}}^{(0)},\ldots,{\bm{y}}^{(2^{m}-1)}\in[0,1)^{s}$ such that

[TABLE]

where $\{{\bm{x}}^{(k)}\}_{k=0}^{2^{m}-1}\in[0,1)^{\alpha s}$ are the points of a polynomial lattice point set $P_{p,m,\alpha s}({\bm{q}})$. A QMC rule using this point set is called an interlaced polynomial lattice rule (of order $\alpha$) with generating vector ${\bm{q}}$ and modulus $p$.

To analyse the error we will make use of the dual of the point set. First we need to define vectors which have a specified support. Therefore, define for an integer vector ${\bm{k}}$ the function $\mathop{\mathrm{supp}}({\bm{k}}):=\{j:k_{j}\neq 0\}$ where the index $j$ ranges over the dimensions of ${\bm{k}}$. To range over all $s$-dimensional vectors with support on the set $\mathfrak{u}$ we write ${\bm{k}}_{\mathfrak{u}}\in\mathbb{N}^{s}_{\mathfrak{u}}$.

Definition 3** (dual of polynomial lattice point set).**

Given $k\in\mathbb{N}_{0}$ with binary expansion $k=\kappa_{0}+\kappa_{1}2+\cdots+\kappa_{a-1}2^{a-1}$ define the associated truncated polynomial

[TABLE]

where $\kappa_{a}=\cdots=\kappa_{m-1}=0$ if $a<m$. For ${\bm{k}}\in\mathbb{N}^{s}_{0}$ define $\operatorname{tr}_{m}({\bm{k}}):=(\operatorname{tr}_{m}(k_{1}),\ldots,\operatorname{tr}_{m}(k_{s}))$. The dual of the polynomial lattice point set $P$ with modulus $p$ and $\deg(p)=m$ and generating vector ${\bm{q}}\in\mathscr{G}_{m}^{s}$ is defined by

[TABLE]

and the dual with support $\mathfrak{u}$ by

[TABLE]

Since interlacing reduces $\alpha$ dimensions to a single dimension, the definition of the dual with support $\mathfrak{u}$ for an interlaced polynomial lattice point set needs to keep track of its source dimensions. Therefore, for $\mathfrak{v}\subseteq\{1:\alpha s\}$, we define, see also [8, Equation (3.26)],

[TABLE]

which tells us where the source dimensions end up in the interlaced point set.

Definition 4** (dual of interlaced polynomial lattice point set).**

Define the digit interlacing function for non-negative integers $\mathscr{E}_{\alpha}:\mathbb{N}_{0}^{\alpha}\to\mathbb{N}_{0}$ with interlacing factor $\alpha\in\mathbb{N}$ by

[TABLE]

where $k_{j}=\kappa_{0,j}+\kappa_{1,j}2+\kappa_{2,j}2^{2}+\cdots$ for $j=1,\ldots,\alpha$ and, in case the number of arguments is a multiple of $\alpha$, define $\mathscr{E}_{\alpha}:\mathbb{N}_{0}^{\alpha s}\to\mathbb{N}_{0}^{s}$ by

[TABLE]

The dual of the interlaced polynomial lattice point set $\mathscr{D}_{\alpha}(P_{p,m,\alpha s}({\bm{q}}))$ is defined by

[TABLE]

where $(P_{p,m,\alpha s}({\bm{q}}))^{\perp}$ is the dual of $P_{p,m,\alpha s}({\bm{q}})$ as given in Definition 3. The dual with support $\mathfrak{u}$ is defined by

[TABLE]

In order to state a bound on the worst-case error we still need to introduce a weight function which measures the importance of the $k$th Walsh basis functions and which will provide a link to the space $H_{\alpha,s}([0,1]^{s})$. For $k\in\mathbb{N}$ with binary expansion $k=\kappa_{1}2^{m_{1}-1}+\kappa_{2}2^{m_{2}-1}+\cdots+\kappa_{v}2^{m_{v}-1}$ such that $m_{1}>m_{2}>\cdots>m_{v}>0$ define

[TABLE]

and $\mu_{\alpha}(0)=0$. For ${\bm{k}}\in\mathbb{N}_{0}^{s}$ we define $\mu_{\alpha}({\bm{k}}):=\sum_{j=1}^{s}\mu_{\alpha}(k_{j})$. We can now state a first bound on the worst-case error for interlaced polynomial lattice rules in the space $H_{\alpha,s}([0,1]^{s})$.

Proposition 10**.**

For any $\alpha\in\mathbb{N}$, with $\alpha\geq 2$, we have

[TABLE]

with

[TABLE]

where $C_{1}:=2^{-1}$ and $C_{\tau}:=(5/3)^{\tau-2}\,2^{-\tau}$ for $\tau\geq 2$.

Proof.

The reproducing kernel for the space $H_{\alpha,s}([0,1]^{s})$ is well known and can be found, e.g., in [4]. For our unweighted tensor product space it is as follows

[TABLE]

where $B_{\tau}$ is the Bernoulli polynomial of degree $\tau\in\mathbb{N}$. We already gave the inner product for this space in (40) and (41). We will expand the kernel in a double Walsh series, see, e.g., [4, 13]. For ${\bm{k}},{\bm{\ell}}\in\mathbb{N}_{0}^{s}$ the $({\bm{k}},{\bm{\ell}})$th Walsh coefficient is defined by

[TABLE]

From, e.g., [4, Theorem 13], we have, with $\{{\bm{y}}^{(i)}\}_{i=0}^{2^{m}-1}\in[0,1)^{s}$ the points of the interlaced polynomial lattice rule,

[TABLE]

Using the “character property”, see, e.g., [13, Lemma 1], we have

[TABLE]

Hence,

[TABLE]

Using [4, Lemma 14 and Equation (13) together with Proposition 20] we have $\widehat{K}_{\alpha,s}({\bm{k}},{\bm{\ell}})=0$ if $\mathop{\mathrm{supp}}({\bm{k}})\neq\mathop{\mathrm{supp}}({\bm{\ell}})$. When both ${\bm{k}}$ and ${\bm{\ell}}$ are $\bm{0}$ we have $\widehat{K}_{\alpha,s}(\bm{0},\bm{0})=1$. Otherwise, for $\mathfrak{u}\neq\emptyset$ and ${\bm{k}}_{\mathfrak{u}},{\bm{\ell}}_{\mathfrak{u}}\in\mathbb{N}^{s}_{\mathfrak{u}}=\{{\bm{k}}\in\mathbb{N}_{0}^{s}:\mathop{\mathrm{supp}}({\bm{k}})=\mathfrak{u}\}$ we have

[TABLE]

with $\widehat{C}_{\alpha}$ defined in (100). Hence

[TABLE]

Thus, taking the square root on both sides and using $\left(\sum_{j}a_{j}\right)^{1/2}\leq\sum_{j}|a_{j}|^{1/2}$ on the right hand side, we obtain

[TABLE]

Using [8, Lemma 3.8 and the subsequent equation] we have

[TABLE]

from which the result follows. ∎

We are now in a similar situation as [8, Equation (3.30)] where a function $E_{d}({\bm{q}})$ is defined which is equal to the upper bound in Proposition 10 with “modified weights”, which in our case would be $\widetilde{\gamma}_{\mathfrak{v}}:=(2^{\alpha(\alpha-1)}\widehat{C}_{\alpha})^{|\mathfrak{u}_{\alpha}(\mathfrak{v})|/2}$, and after which a fast component-by-component construction algorithm is presented. In [8] the weights of the function space $\gamma_{\mathfrak{u}_{\alpha}(\mathfrak{v})}$ are also present in $\widetilde{\gamma}_{\mathfrak{v}}$, but in the unweighted setting here they are all $1$. Hence we can pull out the modified weights using $\widetilde{\gamma}_{\mathfrak{v}}\leq 2^{\alpha(\alpha-1)s/2}$ which holds for all $\mathfrak{v}\subseteq\{1:\alpha s\}$ since $\widehat{C}_{\alpha}<1$, see [4, Table 1 for $q=2$]. The following proposition now follows immediately from [8, Theorem 3.9] by using $\widetilde{\gamma}_{\mathfrak{v}}=1$ for all $\mathfrak{v}$. We note that using the actual weights $\widetilde{\gamma}_{\mathfrak{v}}$ would improve the result, but would not change the complexity for the MDFEM so we prefer this simpler result.

Proposition 11**.**

For any $\alpha\in\mathbb{N}$, with $\alpha\geq 2$, let $p$ be an irreducible polynomial with $\deg(p)=m$. For $d\in\mathbb{N}$ and ${\bm{q}}\in\mathscr{G}_{m}^{d}$ define

[TABLE]

A generating vector ${\bm{q}}^{*}=(q_{1}^{*},q_{2}^{*},\ldots,q_{d}^{*})\in\mathscr{G}_{m}^{d}$ can be constructed using a CBC approach for $d=1,2,\ldots$, minimizing $E_{d}({\bm{q}})$ in each step, such that

[TABLE]

for all $\lambda\in[1,\alpha)$.

Combining Propositions 10 and 11 for $d=\alpha s$ we obtain Theorem 1 in the main text.