Constructing QMC finite element methods for elliptic PDEs with random coefficients by a reduced CBC construction

TL;DR

None

Contribution

None

Abstract

In the analysis of using quasi-Monte Carlo (QMC) methods to approximate expectations of a linear functional of the solution of an elliptic PDE with random diffusion coefficient the sensitivity w.r.t. the parameters is often stated in terms of product-and-order-dependent (POD) weights. The (offline) fast component-by-component (CBC) construction of an $N$-point QMC method making use of these POD weights leads to a cost of $\mathcal{O}(s N \log(N) + s^2 N)$ with $s$ the parameter truncation dimension. When $s$ is large this cost is prohibitive. As an alternative Herrmann and Schwab introduced an analysis resulting in product weights to reduce the construction cost to $\mathcal{O}(s N \log(N))$. We here show how the reduced CBC method can be used for POD weights to reduce the cost to $\mathcal{O}(\sum_{j=1}^{\min\{s,s^*\}} (m-w_j+j) \, b^{m-w_j})$, where $N=b^m$ with prime $b$, $w_1 \le \cdots \le w_s$ are nonnegative integers and $s^*$ can be chosen much smaller than $s$ depending on the regularity of the random field expansion as such making it possible to use the POD weights directly. We show a total error estimate for using randomly shifted lattice rules constructed through the reduced CBC construction.

Equations

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

$$a(\mathbf{x},\mathbf{y})=a_0(\mathbf{x})+\sum _{j\ge 1}{y}_j{\psi }_j(\mathbf{x}),\mathbf{x}\in D,\mathbf{y}\in U.$$

$$\Vert v{\Vert }_V:={\left({\int }_D\sum _{j=1}^d|{\partial }_{x_j}v(\mathbf{x}){|}^2\mathrm{d}\mathbf{x}\right)}^{\frac{1}{2}}=\Vert \nabla v{\Vert }_{L^2(D)}.$$

$$A(\mathbf{y};u(\cdot ,\mathbf{y}),v)=⟨f,v{⟩}_{V^{\ast }\times V}={\int }_{D}f(\mathbf{x})v(\mathbf{x})\mathrm{d}\mathbf{x}\text{mbox}forallv\in V,$$

$$A(\mathbf{y};w,v):={\int }_{D}a(\mathbf{x},\mathbf{y})\nabla w(\mathbf{x})\cdot \nabla v(\mathbf{x})\mathrm{d}\mathbf{x}\text{mbox}forallw,v\in V,$$

$$\mathbb{E}[G(u)]:={\int }_{U}G(u(\cdot ,\mathbf{y}))\mu (\mathrm{d}\mathbf{y})={\int }_{U}G(u(\cdot ,\mathbf{y}))\mathrm{d}\mathbf{y}\approx Q_N(G(u_h^s)),$$

$$\mathbb{E}[G(u)]={\int }_{U}G(u(\cdot ,\mathbf{y}))\mathrm{d}\mathbf{y}:=\underset{s\to \infty }{\lim }{\int }_{{\left[-\frac{1}{2},\frac{1}{2}\right]}^s}G(u(\cdot ,(y_1,\dots ,y_s,0,0,\dots )))\mathrm{d}{y}_1\cdots \mathrm{d}{y}_s$$

$$Q_N(f)=\frac{1}{N}\sum _{k=1}^{N}F({\mathbf{t}}_k),$$

$$a_{0,\min }\le a_0(\mathbf{x})\le a_{0,\max },$$

$${\Vert \sum _{j\ge 1}\frac{|{\psi }_j|}{2a_0}\Vert }_{L^{\infty }(D)}\le \overline{\kappa }<1.$$

$$0<(1-\overline{\kappa })a_{0,\min }\le a_0(\mathbf{x})+\sum _{j\ge 1}{y}_j{\psi }_j(\mathbf{x})=a(\mathbf{x},\mathbf{y})\le (1+\overline{\kappa })a_{0,\max }.$$

$$|\mathbf{\nu }|:=\sum _{j\ge 1}{\nu }_j\text{and}\mathrm{supp}(\mathbf{\nu }):=\{j\ge 1:{\nu }_j>0\}$$

$$\mathcal{F}:=\{\mathbf{\nu }\in {\mathbb{N}}_0^{\mathbb{N}}:\mathrm{supp}(\mathbf{\nu })<\infty \}\text{and}{\mathcal{F}}_1:=\{\mathbf{\nu }\in \{0,1{\}}^{\mathbb{N}}:\mathrm{supp}(\mathbf{\nu })<\infty \},$$

$${\partial }^{\mathbf{\nu }}=\frac{{\partial }^{|\mathbf{\nu }|}}{\partial y_1^{{\nu }_1}\partial y_2^{{\nu }_2}\cdots },$$

$$\Vert v{\Vert }_{a_{\mathbf{y}}}^2:={\int }_{D}a(\mathbf{x},\mathbf{y})|\nabla v(\mathbf{x}){|}^2\mathrm{d}\mathbf{x}$$

$$(1-\overline{\kappa })a_{0,\min }\Vert v{\Vert }_V^2\le \Vert v{\Vert }_{a_{\mathbf{y}}}^2\text{for all}v\in V.$$

$$\Vert u(\cdot ,\mathbf{y}){\Vert }_{a_{\mathbf{y}}}^2$$

$$=⟨f,u(\cdot ,\mathbf{y}){⟩}_{V^{\ast }\times V}\le \Vert f{\Vert }_{V^{\ast }}\Vert u(\cdot ,\mathbf{y}){\Vert }_V\le \frac{\Vert f{\Vert }_{V^{\ast }}\Vert u(\cdot ,\mathbf{y}){\Vert }_{a_{\mathbf{y}}}}{\sqrt{(1-\overline{\kappa })a_{0,\min }}}$$

$$\Vert u(\cdot ,\mathbf{y}){\Vert }_{a_{\mathbf{y}}}^2\le \frac{\Vert f{\Vert }_{V^{\ast }}^2}{(1-\overline{\kappa })a_{0,\min }}.$$

$$\kappa$$

$$\sum _{\begin{array}{c}\mathbf{\omega }\le \mathbf{\nu }\\ |\mathbf{\omega }|=k\end{array}}{\mathbf{b}}^{-2\mathbf{\omega }}\Vert {\partial }^{\mathbf{\omega }}u(\cdot ,\mathbf{y}){\Vert }_V^2\le {\left({\left(\frac{2\tilde{\kappa }(\mathbf{\nu })}{1-\overline{\kappa }}\right)}^k\frac{\Vert f{\Vert }_{V^{\ast }}}{(1-\overline{\kappa })a_{0,\min }}\right)}^2,$$

$$\Vert {\partial }^{\mathbf{\nu }}u(\cdot ,\mathbf{y}){\Vert }_V\le {\mathbf{b}}^{\mathbf{\nu }}{\left(\frac{2\tilde{\kappa }(\mathbf{\nu })}{1-\overline{\kappa }}\right)}^{|\mathbf{\nu }|}\frac{\Vert f{\Vert }_{V^{\ast }}}{(1-\overline{\kappa })a_{0,\min }}.$$

$$\Vert {\partial }^{\mathbf{\nu }}u(\cdot ,\mathbf{y}){\Vert }_{a_{\mathbf{y}}}^2=-\sum _{j\in \mathrm{supp}(\mathbf{\nu })}{\int }_D{\psi }_j(\mathbf{x})\nabla {\partial }^{\mathbf{\nu }-{\mathbf{e}}_j}u(\mathbf{x},\mathbf{y})\cdot \nabla {\partial }^{\mathbf{\nu }}u(\mathbf{x},\mathbf{y})\mathrm{d}\mathbf{x}.$$

$$⟨\mathbf{f},\mathbf{g}{⟩}_{\mathbf{\nu },k}:=\sum _{\begin{array}{c}\mathbf{\omega }\le \mathbf{\nu }\\ |\mathbf{\omega }|=k\end{array}}{\int }_D\sum _{j\in \mathrm{supp}(\mathbf{\omega })}{f}_{\mathbf{\omega },j}(\mathbf{x})g_{\mathbf{\omega },j}(\mathbf{x})\mathrm{d}\mathbf{x}.$$

$$\sum _{\begin{array}{c}\mathbf{\omega }\le \mathbf{\nu }\\ |\mathbf{\omega }|=k\end{array}}{\mathbf{b}}^{-2\mathbf{\omega }}\Vert {\partial }^{\mathbf{\omega }}u(\cdot ,\mathbf{y}){\Vert }_{a_{\mathbf{y}}}^2$$

$$=-\sum _{\begin{array}{c}\mathbf{\omega }\le \mathbf{\nu }\\ |\mathbf{\omega }|=k\end{array}}{\int }_D\sum _{j\in \mathrm{supp}(\mathbf{\omega })}{\mathbf{b}}^{-{\mathbf{e}}_j}{\mathbf{b}}^{-(\mathbf{\omega }-{\mathbf{e}}_j)}{\mathbf{b}}^{-\mathbf{\omega }}{\psi }_j(\mathbf{x})\nabla {\partial }^{\mathbf{\omega }-{\mathbf{e}}_j}u(\mathbf{x},\mathbf{y})\cdot \nabla {\partial }^{\mathbf{\omega }}u(\mathbf{x},\mathbf{y})\mathrm{d}\mathbf{x}$$

$$\le {\left({\int }_D\sum _{\begin{array}{c}\mathbf{\omega }\le \mathbf{\nu }\\ |\mathbf{\omega }|=k\end{array}}\sum _{j\in \mathrm{supp}(\mathbf{\omega })}{\mathbf{b}}^{-{\mathbf{e}}_j}|{\psi }_j(\mathbf{x})|{\left|{\mathbf{b}}^{-(\mathbf{\omega }-{\mathbf{e}}_j)}\nabla {\partial }^{\mathbf{\omega }-{\mathbf{e}}_j}u(\mathbf{x},\mathbf{y})\right|}^2\mathrm{d}\mathbf{x}\right)}^{\frac{1}{2}}$$

$$\times {\left({\int }_D\sum _{\begin{array}{c}\mathbf{\omega }\le \mathbf{\nu }\\ |\mathbf{\omega }|=k\end{array}}\sum _{j\in \mathrm{supp}(\mathbf{\omega })}{\mathbf{b}}^{-{\mathbf{e}}_j}|{\psi }_j(\mathbf{x})|{\left|{\mathbf{b}}^{-\mathbf{\omega }}\nabla {\partial }^{\mathbf{\omega }}u(\mathbf{x},\mathbf{y})\right|}^2\mathrm{d}\mathbf{x}\right)}^{\frac{1}{2}}.$$

$${\int }_D\sum _{\begin{array}{c}\mathbf{\omega }\le \mathbf{\nu }\\ |\mathbf{\omega }|=k\end{array}}\sum _{j\in \mathrm{supp}(\mathbf{\omega })}{\mathbf{b}}^{-{\mathbf{e}}_j}|{\psi }_j(\mathbf{x})|{\left|{\mathbf{b}}^{-(\mathbf{\omega }-{\mathbf{e}}_j)}\nabla {\partial }^{\mathbf{\omega }-{\mathbf{e}}_j}u(\mathbf{x},\mathbf{y})\right|}^2\mathrm{d}\mathbf{x}$$

$$\displaystyle\quad=\int_{D}\sum_{\begin{subarray}{c}{\bm{\omega}}\leq{\bm{\nu}}\\ |{\bm{\omega}}|=k-1\end{subarray}}\Bigg{(}\sum_{\begin{subarray}{c}j\in{\mathrm{supp}}({\bm{\nu}})\\ {\bm{\omega}}+{\bm{e}}_{j}\leq{\bm{\nu}}\end{subarray}}{\bm{b}}^{-{\bm{e}}_{j}}|\psi_{j}({\bm{x}})|\Bigg{)}\left|{\bm{b}}^{-{\bm{\omega}}}\nabla\partial^{{\bm{\omega}}}u({\bm{x}},{\bm{y}})\right|^{2}\,\mathrm{d}{\bm{x}}$$

$$\le {\Vert \sum _{j\in \mathrm{supp}(\mathbf{\nu })}\frac{|{\psi }_j|/b_j}{a(\cdot ,\mathbf{y})}\Vert }_{L^{\infty }(D)}\sum _{\begin{array}{c}\mathbf{\omega }\le \mathbf{\nu }\\ |\mathbf{\omega }|=k-1\end{array}}{\mathbf{b}}^{-2\mathbf{\omega }}{\int }_{D}a(\mathbf{x},\mathbf{y}){\left|\nabla {\partial }^{\mathbf{\omega }}u(\mathbf{x},\mathbf{y})\right|}^2\mathrm{d}\mathbf{x}$$

$$={\Vert \sum _{j\in \mathrm{supp}(\mathbf{\nu })}\frac{|{\psi }_j|/b_j}{a(\cdot ,\mathbf{y})}\Vert }_{L^{\infty }(D)}\sum _{\begin{array}{c}\mathbf{\omega }\le \mathbf{\nu }\\ |\mathbf{\omega }|=k-1\end{array}}{\mathbf{b}}^{-2\mathbf{\omega }}\Vert {\partial }^{\mathbf{\omega }}u(\cdot ,\mathbf{y}){\Vert }_{a_{\mathbf{y}}}^2,$$

$${\Vert \sum _{j\in \mathrm{supp}(\mathbf{\nu })}\frac{|{\psi }_j|/b_j}{a(\cdot ,\mathbf{y})}\Vert }_{L^{\infty }(D)}\le \frac{1}{1-\overline{\kappa }}{\Vert \sum _{j\in \mathrm{supp}(\mathbf{\nu })}\frac{|{\psi }_j|/b_j}{a_0}\Vert }_{L^{\infty }(D)}:=\frac{2\tilde{\kappa }(\mathbf{\nu })}{1-\overline{\kappa }},$$

$$\sum _{\begin{array}{c}\mathbf{\omega }\le \mathbf{\nu }\\ |\mathbf{\omega }|=k\end{array}}{\mathbf{b}}^{-2\mathbf{\omega }}\Vert {\partial }^{\mathbf{\omega }}u(\cdot ,\mathbf{y}){\Vert }_{a_{\mathbf{y}}}^2$$

Full text

1 Introduction and Problem Setting

We consider the parametric elliptic Dirichlet problem given by

[TABLE]

for $D\subset{\mathbb{R}}^{d}$ a bounded, convex Lipschitz polyhedron domain with boundary $\partial D$ and fixed spatial dimension $d\in\{1,2,3\}$. The function $f$ lies in $L^{2}(D)$, the parametric variable ${\bm{y}}=(y_{j})_{j\geq 1}$ belongs to a domain $U$, and the differential operators are understood to be with respect to the physical variable ${\bm{x}}\in D$. Here we study the “uniform case”, i.e., we assume that ${\bm{y}}$ is uniformly distributed on $U:=\left[-\frac{1}{2},\frac{1}{2}\right]^{{\mathbb{N}}}$ with uniform probability measure $\mu({\mathrm{d}}{\bm{y}})=\bigotimes_{j\geq 1}{\mathrm{d}}y_{j}={\mathrm{d}}{\bm{y}}$. The parametric diffusion coefficient $a({\bm{x}},{\bm{y}})$ is assumed to depend linearly on the parameters $y_{j}$ in the following way,

[TABLE]

For the variational formulation of (1.1), we consider the Sobolev space $V=H_{0}^{1}(D)$ of functions $v$ which vanish on the boundary $\partial D$ with norm

[TABLE]

The corresponding dual space of bounded linear functionals on $V$ with respect to the pivot space $L^{2}(D)$ is further denoted by $V^{*}=H^{-1}(D)$. Then, for given $f\in V^{*}$ and ${\bm{y}}\in U$, the weak (or variational) formulation of (1.1) is to find $u(\cdot,{\bm{y}})\in V$ such that

[TABLE]

with parametric bilinear form $A:U\times V\times V\to{\mathbb{R}}$ given by

[TABLE]

and duality pairing $\langle\cdot,\cdot\rangle_{V^{*}\times V}$ between $V^{*}$ and $V$. We will often identify elements $\varphi\in V$ with dual elements $L_{\varphi}\in V^{*}$. Indeed, for $\varphi\in V$ and $v\in V$, a bounded linear functional is given via $L_{\varphi}(v):=\int_{D}\varphi({\bm{x}})v({\bm{x}})\,\mathrm{d}{\bm{x}}=\langle\varphi,v\rangle_{L^{2}(D)}$ and by the Riesz representation theorem there exists a unique representer $\widetilde{\varphi}\in V$ such that $L_{\varphi}(v)=\langle\widetilde{\varphi},v\rangle_{L^{2}(D)}$ for all $v\in V$. Hence, the definition of the canonical duality pairing yields that $\langle L_{\varphi},v\rangle_{V^{*}\times V}=L_{\varphi}(v)=\langle\varphi,v\rangle_{L^{2}(D)}$.

Our quantity of interest is the expected value, with respect to ${\bm{y}}\in U$, of a given bounded linear functional $G\in V^{*}$ applied to the solution $u(\cdot,{\bm{y}})$ of the PDE. We therefore seek to approximate this expectation by numerically integrating $G$ applied to a finite element approximation $u_{h}^{s}(\cdot,{\bm{y}})$ of the solution $u^{s}(\cdot,{\bm{y}})\in H_{0}^{1}(D)=V$ of (1.3) with truncated diffusion coefficient $a({\bm{x}},({\bm{y}}_{\{1:s\}};0))$ where $\{1:s\}:=\{1,\ldots,s\}$ and we write $({\bm{y}}_{\{1:s\}};0)=(\tilde{y}_{j})_{j\geq 1}$ with $\tilde{y}_{j}=y_{j}$ for $j\in\{1:s\}$ and $\tilde{y}_{j}=0$ otherwise; that is,

[TABLE]

with $Q_{N}(\cdot)$ a linear quadrature rule using $N$ function evaluations. The infinite-dimensional integral ${\mathbb{E}}[G(u)]$ in (1.5) is defined as

[TABLE]

such that our integrands of interest are of the form $F({\bm{y}})=G(u(\cdot,{\bm{y}}))$ with ${\bm{y}}\in U$. In this article, we will employ (randomized) QMC methods of the form

[TABLE]

i.e., equal-weight quadrature rules with (randomly shifted) deterministic points ${\bm{t}}_{1},\ldots,{\bm{t}}_{N}\in\left[-\frac{1}{2},\frac{1}{2}\right]^{s}$. This elliptic PDE is a standard problem considered in the numerical analysis of computational methods in uncertainty quantification, see, e.g., [1, 2, 6, 8, 9, 10, 11, 12].

1.1 Existence of solutions of the variational problem

To assure that a unique solution to the weak problem (1.3) exists, we need certain conditions on the diffusion coefficient $a$. We assume $a_{0}\in L^{\infty}(D)$ and $\operatorname*{ess\,inf}_{{\bm{x}}\in D}a_{0}({\bm{x}})>0$, which is equivalent to the existence of two constants $0<a_{0,\min}\leq a_{0,\max}<\infty$ such that a.e. on $D$ we have

[TABLE]

and that there exists a $\overline{\kappa}\in(0,1)$ such that

[TABLE]

Via (1.7), we obtain that $|\sum_{j\geq 1}y_{j}\psi_{j}({\bm{x}})|\leq\overline{\kappa}\,a_{0}({\bm{x}})$ and hence, using (1.6), almost everywhere on $D$ and for any ${\bm{y}}\in U$

[TABLE]

These estimates yield the continuity and coercivity of $A({\bm{y}},\cdot,\cdot)$ defined in (1.4) on $V\times V$, uniformly for all ${\bm{y}}\in U$. The Lax–Milgram theorem then ensures the existence of a unique solution $u(\cdot,{\bm{y}})$ of the weak problem in (1.3).

1.2 Parametric regularity

Having established the existence of unique weak parametric solutions $u(\cdot,{\bm{y}})$, we investigate their regularity in terms of the behaviour of their mixed first-order derivatives. Our analysis combines multiple techniques which can be found in the literature, see, e.g., [2, 8, 10, 9, 1]. In particular we want to point out that our POD form bounds can take advantage of wavelet like expansions of the random field, a technique introduced in [1] and used to the advantage of QMC constructions by [9] to deliver product weights to save on the construction compared to POD weights. Although we end up again with POD weights, we will save on the construction cost by making use of a special construction method, called the reduced CBC construction, which we will introduce in Section 2.4. Let ${\bm{\nu}}=(\nu_{j})_{j\geq 1}$ with $\nu_{j}\in{\mathbb{N}}_{0}:=\{0,1,2,\ldots\}$ be a sequence of positive integers which we will refer to as a multi-index. We define the order $|{\bm{\nu}}|$ and the support ${\mathrm{supp}}({\bm{\nu}})$ as

[TABLE]

and introduce the sets ${\mathcal{F}}$ and ${\mathcal{F}}_{1}$ of finitely supported multi-indices as

[TABLE]

where ${\mathcal{F}}_{1}\subseteq{\mathcal{F}}$ is the restriction containing only ${\bm{\nu}}$ with $\nu_{j}\in\{0,1\}$. Then, for ${\bm{\nu}}\in{\mathcal{F}}$ denote the ${\bm{\nu}}$-th partial derivative with respect to the parametric variables ${\bm{y}}\in U$ by

[TABLE]

and for a sequence ${\bm{b}}=(b_{j})_{j\geq 1}\subset{\mathbb{R}}^{{\mathbb{N}}}$, set ${\bm{b}}^{\bm{\nu}}:=\prod_{j\geq 1}b_{j}^{\nu_{j}}$. We further write ${\bm{\omega}}\leq{\bm{\nu}}$ if $\omega_{j}\leq\nu_{j}$ for all $j\geq 1$ and denote by ${\bm{e}}_{i}\in{\mathcal{F}}_{1}$ the multi-index with components $e_{j}=\delta_{i,j}$. For a fixed ${\bm{y}}\in U$, we introduce the energy norm $\|\cdot\|_{a_{\bm{y}}}^{2}$ in the space $V$ via

[TABLE]

for which it holds true by (1.8) that

[TABLE]

Consequently, we have that $(1-\overline{\kappa})\,a_{0,\min}\|u(\cdot,{\bm{y}})\|_{V}^{2}\leq\|u(\cdot,{\bm{y}})\|_{a_{\bm{y}}}^{2}$ and hence the definition of the dual norm $\|\cdot\|_{V^{*}}$ yields the following initial estimate from (1.3) and (1.4),

[TABLE]

which gives in turn

[TABLE]

In order to exploit the decay of the norm sequence $(\|\psi_{j}\|_{L^{\infty}(D)})_{j\geq 1}$ of the basis functions, we extend condition (1.7) as follows. To characterize the smoothness of the random field, we assume that there exist a sequence of reals ${\bm{b}}=(b_{j})_{j\geq 1}$ with $0<b_{j}\leq 1$ for all $j$, a constant $\kappa\in(0,1)$ and therefore also constants $\widetilde{\kappa}({\bm{\nu}})\leq\kappa$ for all ${\bm{\nu}}\in{\mathcal{F}}_{1}$ such that

[TABLE]

We remark that condition (1.7) is included in this assumption by letting $b_{j}=1$ for all $j\geq 1$ and that $0<\overline{\kappa}\leq\kappa<1$. Using the above estimations we can derive the following theorem for the mixed first-order partial derivatives.

Theorem 1**.**

Let ${\bm{\nu}}\in{\mathcal{F}}_{1}$ be a multi-index of finite support and let $k\in\{0,1,\ldots,|{\bm{\nu}}|\}$. Then, for every $f\in V^{*}$ and every ${\bm{y}}\in U$,

[TABLE]

with $\widetilde{\kappa}({\bm{\nu}})$ as in (1.11). Moreover, for $k=|{\bm{\nu}}|$ we obtain

[TABLE]

Proof.

For the special case ${\bm{\nu}}={\bm{0}}$, the claim follows by combining (1.9) and (1.10). For ${\bm{\nu}}\in{\mathcal{F}}_{1}$ with $\left|{\bm{\nu}}\right|>0$, as is known from, e.g., [2] and [11, Appendix], the linearity of $a({\bm{x}},{\bm{y}})$ gives rise to the following identity for any ${\bm{y}}\in U$:

[TABLE]

For sequences of $L^{2}(D)$-integrable functions ${\bm{f}}=(f_{{\bm{\omega}},j})_{{\bm{\omega}}\in{\mathcal{F}},j\geq 1}$ with $f_{{\bm{\omega}},j}:D\to{\mathbb{R}}$, we define the inner product $\langle{\bm{f}},{\bm{g}}\rangle_{{\bm{\nu}},k}$ as follows,

[TABLE]

We can then apply the Cauchy–Schwarz inequality to ${\bm{f}}=(f_{{\bm{\omega}},j})$ and ${\bm{g}}=(g_{{\bm{\omega}},j})$ with $f_{{\bm{\omega}},j}={\bm{b}}^{-{\bm{e}}_{j}/2}|\psi_{j}|^{\frac{1}{2}}\,{\bm{b}}^{-({\bm{\omega}}-{\bm{e}}_{j})}\nabla\partial^{{\bm{\omega}}-{\bm{e}}_{j}}u(\cdot,{\bm{y}})$ and $g_{{\bm{\omega}},j}={\bm{b}}^{-{\bm{e}}_{j}/2}|\psi_{j}|^{\frac{1}{2}}\,{\bm{b}}^{-{\bm{\omega}}}\nabla\partial^{{\bm{\omega}}}u(\cdot,{\bm{y}})$ to obtain, with the help of (1.12),

[TABLE]

The first of the two factors above is then bounded as follows,

[TABLE]

while the other factor can be bounded trivially. Furthermore, using (1.8), we have for any ${\bm{y}}\in U$

[TABLE]

so that, combining these three estimates, we obtain

[TABLE]

Therefore, we finally obtain that

[TABLE]

which inductively gives

[TABLE]

where the last inequality follows from the initial estimate (1.10). The estimate (1.9) then gives

[TABLE]

which yields the first claim. The second claim follows since the sum over the ${\bm{\omega}}\leq{\bm{\nu}}$ with $|{\bm{\omega}}|=|{\bm{\nu}}|$ and ${\bm{\nu}}\in{\mathcal{F}}_{1}$ consists only of the term corresponding to ${\bm{\omega}}={\bm{\nu}}$. ∎

Corollary 1**.**

Under the assumptions of Theorem 1, there exists a number $\kappa(k)$ for each $k\in{\mathbb{N}}$, given by

[TABLE]

such that $\widetilde{\kappa}({\bm{\nu}})\leq\kappa(k)\leq\kappa<1$ for all ${\bm{\nu}}\in{\mathcal{F}}_{1}$ with $|{\bm{\nu}}|=k$. Then for ${\bm{\nu}}\in{\mathcal{F}}_{1}$, every $f\in V^{\ast}$, and every ${\bm{y}}\in U$, the solution $u(\cdot,{\bm{y}})$ satisfies

[TABLE]

Note that since $0<\overline{\kappa}\leq\kappa<1$, the results of Theorem 1 and Corollary 1 remain also valid for $\overline{\kappa}$ replaced by $\kappa$.

The obtained bounds on the mixed first-order derivatives turn out to be of product and order-dependent (so-called POD) form; that is, they are of the general form

[TABLE]

with a map $\Gamma:{\mathbb{N}}_{0}\to{\mathbb{R}}$, a sequence of reals ${\bm{b}}=(b_{j})_{j\geq 1}\in{\mathbb{R}}^{{\mathbb{N}}}$ and some constant $C\in{\mathbb{R}}_{+}$. This finding motivates us to consider this special type of bounds in the following error analysis.

2 Quasi-Monte Carlo finite element error

We analyze the error ${\mathbb{E}}[G(u)]-Q_{N}(G(u_{h}^{s}))$ obtained by applying QMC rules to the finite element approximation $u_{h}^{s}$ to approximate the expected value

[TABLE]

To this end, we introduce the finite element approximation $u_{h}^{s}({\bm{x}},{\bm{y}}):=u_{h}({\bm{x}},({\bm{y}}_{\{1:s\}};0))$ of a solution of (1.3) with truncated diffusion coefficient $a({\bm{x}},({\bm{y}}_{\{1:s\}};0))$, where $u_{h}$ is a finite element approximation as defined in (2.2) and $({\bm{y}}_{\{1:s\}};0)=(y_{1},\ldots,y_{s},0,0,\ldots)$. The overall absolute QMC finite element error is then bounded as follows

[TABLE]

The first term on the right hand side of (2.1) will be referred to as (dimension) truncation error, the second term is the finite element discretization error and the last term is the QMC quadrature error for the integrand $u_{h}^{s}$. In the following sections we will analyze these different error terms separately.

2.1 Finite Element Approximation

Here, we consider the approximation of the solution $u(\cdot,{\bm{y}})$ of (1.3) by a finite element approximation $u_{h}(\cdot,{\bm{y}})$ and assess the finite element discretization error. More specifically, denote by $\{V_{h}\}_{h>0}$ a family of subspaces $V_{h}\subset V$ of finite dimension $M_{h}$ such that $V_{h}\to V$ as $h\to 0$. We define the parametric finite element (FE) approximation as follows: for $f\in V^{*}$ and given ${\bm{y}}\in U$, find $u_{h}(\cdot,{\bm{y}})\in V_{h}$ such that

[TABLE]

To establish convergence of the finite element approximations, we need some further conditions on $a({\bm{x}},{\bm{y}})$. To this end, we define the space $W^{1,\infty}(D)\subseteq L^{\infty}(D)$ endowed with the norm $\|v\|_{W^{1,\infty}(D)}=\max\{\|v\|_{L^{\infty}(D)},\|\nabla v\|_{L^{\infty}(D)}\}$ and require that

[TABLE]

Under these conditions and using that $f\in L^{2}(D)$, it was proven in [12, Theorems 7.1 and 7.2] that for any ${\bm{y}}\in U$ the approximations $u_{h}(\cdot,{\bm{y}})$ satisfy

[TABLE]

In addition, if (the representer of) the bounded linear functional $G\in V^{*}$ lies in $L^{2}(D)$ we have for any ${\bm{y}}\in U$, as $h\to 0$,

[TABLE]

where the constants $C_{1},C_{2},C_{3}>0$ are independent of $h$ and ${\bm{y}}$. Since the above statements hold true for any ${\bm{y}}\in U$, they remain also valid for $u^{s}({\bm{x}},{\bm{y}}):=u({\bm{x}},({\bm{y}}_{\{1:s\}};0))$ and $u_{h}^{s}({\bm{x}},{\bm{y}}):=u_{h}({\bm{x}},({\bm{y}}_{\{1:s\}};0))$.

2.2 Dimension Truncation

For every $s\in{\mathbb{N}}$ and ${\bm{y}}\in U$, we formally define the solution of the parametric weak problem (1.3) corresponding to the diffusion coefficient $a({\bm{x}},({\bm{y}}_{\{1:s\}};0))$ with sum truncated to $s$ terms as

[TABLE]

In [8, Proposition 5.1] it was shown that for the solution $u^{s}$ the following error estimates are satisfied.

Theorem 2**.**

Let $\overline{\kappa}\in(0,1)$ be such that (1.7) is satisfied and assume furthermore that there exists a sequence of reals ${\bm{b}}=(b_{j})_{j\geq 1}$ with $0<b_{j}\leq 1$ for all $j$ and a constant $\kappa\in[\overline{\kappa},1)$ as defined in (1.11). Then, for every ${\bm{y}}\in U$ and each $s\in{\mathbb{N}}$

[TABLE]

Moreover, if it holds for $\kappa$ that $\frac{\kappa\,a_{0,\max}}{(1-\overline{\kappa})\,a_{0,\min}}\sup_{j\geq s+1}b_{j}<1$, then for every $G\in V^{*}$ we have

[TABLE]

In the following subsection, we will discuss how to approximate the finite-dimensional integral of solutions of the form (2.5) by means of QMC methods.

2.3 Quasi-Monte Carlo Integration

For a real-valued function $F:[-\tfrac{1}{2},\tfrac{1}{2}]^{s}\to{\mathbb{R}}$ defined over the $s$-dimensional unit cube centered at the origin, we consider the approximation of the integral $I_{s}(F)$ by $N$-point QMC rules $Q_{N}(F)$, i.e.,

[TABLE]

with quadrature points ${\bm{t}}_{1},\ldots,{\bm{t}}_{N}\in[-\tfrac{1}{2},\tfrac{1}{2}]^{s}$. As a quality criterion of such a rule, we define the worst-case error for QMC integration in some Banach space ${\mathcal{H}}$ as

[TABLE]

In this article, we consider randomly shifted rank-1 lattice rules as randomized QMC rules, with underlying points of the form

[TABLE]

with generating vector ${\bm{z}}\in{\mathbb{Z}}^{s}$, uniform random shift ${\bm{\Delta}}\in[0,1]^{s}$ and component-wise applied fractional part, denoted by $\{{\bm{x}}\}$. For simplicity, we denote the worst-case error using a shifted lattice rule with generating vector ${\bm{z}}$ and shift ${\bm{\Delta}}$ by $e_{N,s}({\bm{z}},{\bm{\Delta}})$.

For randomly shifted QMC rules, the probabilistic error bound

[TABLE]

holds for all $F\in{\mathcal{H}}$, with shift-averaged worst-case error

[TABLE]

As function space ${\mathcal{H}}$ for our integrands $F$, we consider the weighted, unanchored Sobolev space ${\mathcal{W}}_{s,{\bm{\gamma}}}$, which is a Hilbert space of functions defined over $[-\frac{1}{2},\frac{1}{2}]^{s}$ with square integrable mixed first derivatives and general non-negative weights ${\bm{\gamma}}=(\gamma_{{\mathfrak{u}}})_{{\mathfrak{u}}\subseteq\{1:s\}}$. More precisely, the norm for $F\in{\mathcal{W}}_{s,{\bm{\gamma}}}$ is given by

[TABLE]

where $\{1:s\}:=\{1,\ldots,s\}$, $\frac{\partial^{|{\mathfrak{u}}|}F}{\partial{\bm{y}}_{\mathfrak{u}}}$ denotes the mixed first derivative with respect to the variables ${\bm{y}}_{\mathfrak{u}}=(y_{j})_{j\in{\mathfrak{u}}}$ and we set ${\bm{y}}_{-{\mathfrak{u}}}=(y_{j})_{j\in\{1:s\}\setminus{\mathfrak{u}}}$.

For the efficient construction of good lattice rule generating vectors, we consider the so-called reduced component-by-component (CBC) construction introduced in [5]. For $b\in{\mathbb{N}}$ and $m\in{\mathbb{N}}_{0}$, we define the group of units of integers modulo $b^{m}$ via

[TABLE]

and note that ${\mathbb{Z}}_{b^{0}}^{\times}={\mathbb{Z}}_{1}^{\times}=\{0\}$ since $\gcd(0,1)=1$. Henceforth, let $b$ be prime and recall that then, for $m\geq 1$, $|{\mathbb{Z}}_{b^{m}}^{\times}|=\varphi(b^{m})=b^{m-1}\varphi(b)$ and $|{\mathbb{Z}}_{b}^{\times}|=\varphi(b)=(b-1)$, where $\varphi$ is Euler’s totient function. Let ${\bm{w}}:=(w_{j})_{j\geq 1}$ be a non-decreasing sequence of integers in ${\mathbb{N}}_{0}$, the elements of which we will refer to as reduction indices. In the reduced CBC algorithm the components $\widetilde{z}_{j}$ of the generating vector $\widetilde{{\bm{z}}}$ of the lattice rule will be taken as multiples of $b^{w_{j}}$.

In [5], the reduced CBC construction was introduced to construct rank-1 lattice rules for $1$-periodic functions in a weighted Korobov space ${\mathcal{H}}(K_{s,\alpha,{\bm{\gamma}}})$ of smoothness $\alpha$ (see, e.g., [14]). We denote the worst-case error in ${\mathcal{H}}(K_{s,\alpha,{\bm{\gamma}}})$ using a rank-1 lattice rule with generating vector ${\bm{z}}$ by $e_{N,s}({\bm{z}})$. Following [5], the reduced CBC construction is then given in Algorithm 1.

The following theorem, proven in [5], states that the algorithm yields generating vectors with a small integration error for general weights $\gamma_{{\mathfrak{u}}}$ in the Korobov space.

Theorem 3**.**

For a prime power $N=b^{m}$ let $\widetilde{{\bm{z}}}=(b^{w_{1}}z_{1},\ldots,b^{w_{s}}z_{s})$ be constructed according to Algorithm 1 with integer reduction indices $0\leq w_{1}\leq\cdots\leq w_{s}$. Then for every $d\in\{1:s\}$ and every $\lambda\in(1/\alpha,1]$ it holds for the worst-case error in the Korobov space ${\mathcal{H}}(K_{s,\alpha,{\bm{\gamma}}})$ with $\alpha>1$ that

[TABLE]

This theorem can be extended to the weighted unanchored Sobolev space ${\mathcal{W}}_{s,{\bm{\gamma}}}$ using randomly shifted lattice rules as follows.

Theorem 4**.**

For a prime power $N=b^{m}$, $m\in{\mathbb{N}}_{0}$, and for $F\in{\mathcal{W}}_{s,{\bm{\gamma}}}$ belonging to the weighted unanchored Sobolev space defined over $[-\frac{1}{2},\frac{1}{2}]^{s}$ with weights ${\bm{\gamma}}=(\gamma_{{\mathfrak{u}}})_{{\mathfrak{u}}\subseteq\{1:s\}}$, a randomly shifted lattice rule can be constructed by the reduced CBC algorithm, see Algorithm 1, such that for all $\lambda\in(1/2,1]$,

[TABLE]

with integer reduction indices $0\leq w_{1}\leq\cdots\leq w_{s}$ and $\varrho(\lambda)=2\zeta(2\lambda)(2\pi^{2})^{-\lambda}$.

Proof.

Using Theorem 3 and the connection that the shift-averaged kernel of the Sobolev space equals the kernel of the Korobov space ${\mathcal{H}}(K_{s,\alpha,\widetilde{{\bm{\gamma}}}})$ with $\alpha=2$ and weights $\widetilde{\gamma}_{\mathfrak{u}}=\gamma_{\mathfrak{u}}/(2\pi^{2})^{|{\mathfrak{u}}|}$, see, e.g., [7, 13], the result follows from

[TABLE]

It follows that we can construct the lattice rule in the weighted Korobov space using the connection mentioned in the proof of the previous theorem.

2.4 Implementation of the reduced CBC algorithm

Similar to other variants of the CBC construction, we present a fast version of the reduced CBC method for POD weights in Algorithm 2 for which Theorems 3 and 4 still hold. The full derivation of Algorithm 2 is given in Section 5, here we only introduce the necessary notation. The squared worst-case error for POD weights ${\bm{\gamma}}=(\gamma_{{\mathfrak{u}}})_{{\mathfrak{u}}\subseteq\{1:s\}}$ with $\gamma_{{\mathfrak{u}}}=\Gamma(|{\mathfrak{u}}|)\prod_{j\in{\mathfrak{u}}}\gamma_{j}$ and $\gamma_{\emptyset}=1$ in the weighted Korobov space ${\mathcal{H}}(K_{s,\alpha,{\bm{\gamma}}})$ with $\alpha>1$ can be written as

[TABLE]

where $\omega(x)=\sum_{0\neq h\in{\mathbb{Z}}}{\mathrm{e}}^{2\pi{\mathrm{i}}\,hx}/\left|h\right|^{\alpha}$, see, e.g., [7, 13], and for $n\in{\mathbb{N}}$ we define $\Omega_{n}$ as

[TABLE]

We assume that the values of the function $\omega$ can be computed at unit cost. For integers $0\leq w^{\prime}\leq w^{\prime\prime}\leq m$ and given base $b$ we define the “fold and sum” operator, which divides a length $b^{m-w^{\prime}}$ vector into blocks of equal length $b^{m-w^{\prime\prime}}$ and sums them up, i.e.,

[TABLE]

where $\cramped{I_{b^{m-w^{\prime\prime}}}}$ is the identity matrix of size $\cramped{b^{m-w^{\prime\prime}}\times b^{m-w^{\prime\prime}}}$. The computational cost of applying $P_{w^{\prime\prime},w^{\prime}}^{m}$ is the length of the input vector $\cramped{{\mathcal{O}}(b^{m-w^{\prime}})}$. It should be clear that $\cramped{P_{w^{\prime\prime\prime},w^{\prime\prime}}^{m}\,P_{w^{\prime\prime},w^{\prime}}^{m}\,{\bm{v}}=P_{w^{\prime\prime\prime},w^{\prime}}^{m}\,{\bm{v}}}$ for $0\leq w^{\prime}\leq w^{\prime\prime}\leq w^{\prime\prime\prime}\leq m$. In step 4 of Algorithm 2 the notation $.*$ denotes the element-wise product of two vectors and $\Omega_{b^{m-w_{j}}}(z_{j},:)$ means to take the row corresponding to $z=z_{j}$ from the matrix. Furthermore, Algorithm 2 includes an optional step in which the reduction indices are adjusted in case $w_{1}>0$, the auxiliary variable $w_{0}=0$ is introduced to satisfy the recurrence relation.

The standard fast CBC algorithm for POD weights has a complexity of ${\mathcal{O}}(s\,N\log N+s^{2}N)$, see, e.g., [7, 13]. The cost of our new algorithm can be substantially lower as is stated in the following theorem. We stress that the presented algorithm is the first realization of the reduced CBC construction for POD weights. Our new algorithm improves upon the one stated in [5] which only considers product weights, but the same technique can be used there since POD weights are more general and include product weights.

Theorem 5**.**

Given a sequence of integer reduction indices $0\leq w_{1}\leq w_{2}\leq\cdots$, the reduced CBC algorithm for a prime power $N=b^{m}$ points in $s$ dimensions as specified in Algorithm 2 can construct a lattice rule with near optimal worst-case error as in Theorem 4 with an arithmetic cost of

[TABLE]

where $s^{*}$ is defined to be the largest integer such that $w_{s^{*}}<m$. The memory cost is ${\mathcal{O}}(\sum_{j=1}^{\min\{s,s^{\ast}\}}b^{m-w_{j}})$. In case of product weights ${\mathcal{O}}(\sum_{j=1}^{\min\{s,s^{\ast}\}}(m-w_{j})\,b^{m-w_{j}})$ operations are required for the construction with memory ${\mathcal{O}}(b^{m-w_{1}})$.

Proof.

We refer to Algorithm 2. Step 1 can be calculated in ${\mathcal{O}}(j\,b^{m-w_{j-1}})$ operations (and we may assume $w_{0}=w_{1}$ since the case $w_{1}>0$ can be reduced to the case $w_{1}=0$). The matrix-vector multiplication in step 2 can be done by exploiting the block-circulant structure to obtain a fast matrix-vector product by FFTs at a cost of ${\mathcal{O}}((m-w_{j})\,b^{m-w_{j}})$, see, e.g., [3, 4]. We ignore the possible saving by pre-computation of FFTs on the first columns of the blocks in the matrices $\Omega_{b^{m-w_{j}}}$ as this has cost ${\mathcal{O}}((m-w_{1})\,b^{m-w_{1}})$ and therefore is already included in the cost of step 2. Finally, the vectors ${\bm{q}}_{j,\ell}$ for $\ell=1,\ldots,j$ in step 4 can be calculated in ${\mathcal{O}}(j\,b^{m-w_{j-1}})$. To obtain the total complexity we remark that the applications of the “fold and sum” operator, marked by the square brackets could be performed in iteration $j-1$ such that the cost of steps 1 and 4 in iteration $j$ are only ${\mathcal{O}}(j\,b^{m-w_{j}})$ instead of ${\mathcal{O}}(j\,b^{m-w_{j-1}})$. The cost of the additional fold and sum to prepare for iteration $j$ in iteration $j-1$, which can be performed after step 4, is then equal to the cost of step 4 in that iteration. Since we can assume $w_{0}=w_{1}$ we obtain the claimed construction cost. Note that the algorithm is written in such a way that the vectors ${\bm{q}}_{j,\ell-1}$ can be reused for storing the vectors ${\bm{q}}_{j,\ell}$ (which might be smaller). Similarly for the vectors $\overline{{\bm{q}}}_{j}$. Therefore the memory cost is ${\mathcal{O}}(\sum_{j=1}^{\min\{s,s^{\ast}\}}b^{m-w_{j}})$. The result for product weights can be obtained similarly, see, e.g., [13]. ∎

3 QMC finite element error analysis

We now combine the results of the previous subsections to analyze the overall QMC finite element error. We consider the root mean square error (RMSE) given by

[TABLE]

The error ${\mathbb{E}}[G(u)]-Q_{N}(G(u_{h}^{s}))$ can be written as

[TABLE]

such that due to the fact that ${\mathbb{E}}_{{\bm{\Delta}}}(Q_{N}(f))=I_{s}(f)$ for any integrand $f$ we obtain

[TABLE]

Then, noting that ${\mathbb{E}}[G(u)]-I_{s}(G(u_{h}^{s}))={\mathbb{E}}[G(u)]-I_{s}(G(u^{s}))+I_{s}(G(u^{s}))-I_{s}(G(u_{h}^{s}))$,

[TABLE]

and since for general $x,y\in{\mathbb{R}}$ it holds that $2xy\leq x^{2}+y^{2}$, we obtain furthermore

[TABLE]

From the previous subsections we can then use (2.6) for the truncation part, (2.1), which holds for general ${\bm{y}}\in U$ and thus also for $y_{\{1:s\}}$, for the finite element error, and Theorem 4 for the QMC integration error to obtain the following error bound for the mean square error ${\mathbb{E}}_{{\bm{\Delta}}}[|{\mathbb{E}}[G(u)]-Q_{N}(G(u_{h}^{s}))|^{2}]=:e^{\text{MSE}}_{N,s,h}(G(u))$,

[TABLE]

for some constants $K_{1},K_{2}\in{\mathbb{R}}_{+}$ and provided that $\frac{a_{0,\max}}{(1-\overline{\kappa})\,a_{0,\min}}\kappa\,\sup_{j\geq s+1}b_{j}<1$.

3.1 Derivative bounds of POD form

In the following we assume that we have general bounds on the mixed partial derivatives $\partial^{{\bm{\nu}}}u(\cdot,{\bm{y}})$ which are of POD form; that is,

[TABLE]

with a map $\Gamma:{\mathbb{N}}_{0}\to{\mathbb{R}}$, a sequence of reals $\widetilde{{\bm{b}}}=(\widetilde{b}_{j})_{j\geq 1}\in{\mathbb{R}}^{{\mathbb{N}}}$ and some constant $C\in{\mathbb{R}}_{+}$. Such bounds can be found in the literature and we provided a new derivation in Theorem 1 also leading to POD weights.

For bounding the norm $\|G(u_{h}^{s})\|_{{\mathcal{W}}_{s,{\bm{\gamma}}}}$, we can then use (3.2) and the definition in (2.7) to proceed as outlined in [11], to obtain the estimate

[TABLE]

Denoting ${\bm{w}}:=(w_{j})_{j\geq 1}$ and using (3.3), the contribution of the quadrature error to the mean square error $e_{N,h,s}^{\text{MSE}}(G(u))$ can be upper bounded by

[TABLE]

where we define

[TABLE]

The term $C_{{\bm{\gamma}},{\bm{w}},\lambda}$ can be bounded as

[TABLE]

Due to [12, Lemma 6.2] the latter term is minimized by choosing the weights $\gamma_{\mathfrak{u}}$ as

[TABLE]

Then we set

[TABLE]

and easily see that also

[TABLE]

which implies that $C_{{\bm{\gamma}},{\bm{w}},\lambda}\leq A_{\lambda}^{1+1/\lambda}$. We demonstrate how the term $A_{\lambda}$ can be estimated for the derivative bounds derived in Section 1.2.

In view of Theorem 1, assume in the following that

[TABLE]

Note that we could also choose $\Gamma(\left|{\mathfrak{u}}\right|)=\kappa(\left|{\mathfrak{u}}\right|)^{\left|{\mathfrak{u}}\right|}$ above, in which case the subsequent estimate of $A_{\lambda}$ can be done analogously, but to make the argument less technical, we consider the slightly coarser variant $\Gamma(\left|{\mathfrak{u}}\right|)=\kappa^{\left|{\mathfrak{u}}\right|}$ here. In this case,

[TABLE]

Note that, as $\lambda\leq 1$, it holds that $b^{\frac{-w_{\ell}}{2\lambda}}\leq b^{\frac{-w_{\ell}}{2}}$ and hence

[TABLE]

We now proceed similarly to the proof of Theorem 6.4 in [12]. Let $(\alpha_{j})_{j\geq 1}$ be a sequence of positive reals, to be specified below, which satisfies $\Sigma:=\sum_{j=1}^{\infty}\alpha_{j}<\infty$. Dividing and multiplying by $\prod_{j\in{\mathfrak{u}}}\alpha_{j}^{(2\lambda)/(1+\lambda)}$, and applying Hölder’s inequality with conjugate components $p=(1+\lambda)/(2\lambda)$ and $p^{*}=(1+\lambda)/(1-\lambda)$,

[TABLE]

where we define

[TABLE]

For the first factor we estimate

[TABLE]

By the ratio test, the latter expression is finite if we choose $(\alpha_{j})_{j\geq 1}$ such that $L:=\sup_{k\in{\mathbb{N}}}\kappa\,b^{\frac{-w_{k}}{2}}(k+1)^{-1}=\kappa b^{\frac{-w_{1}}{2}}/2<1/\Sigma$. Hence we assume that $(\alpha_{j})_{j\geq 1}$ is chosen such that indeed $L<1/\Sigma$. Note that $L$ is small if $\kappa$ is small, which means that $\Sigma$ can be allowed to be large in this case. Consider now the term

[TABLE]

We require

[TABLE]

To this end, we choose $\alpha_{j}:=\frac{\left(b_{j}b^{w_{j}}\right)^{p}}{\theta}$, where $\frac{\theta}{L}>\sum_{j=1}^{\infty}\left(b_{j}b^{w_{j}}\right)^{p}$. Then,

[TABLE]

as long as we choose $\lambda$ such that

[TABLE]

We denote the upper bound in (3.8) by $\overline{A}(\lambda)$. Similarly to what is done in [12, Proof of Theorem 6.4], we see that Condition (3.9) is satisfied if $\lambda\geq\frac{p}{2-p}$. Again, similarly to [12, Proof of Theorem 6.4] we see that the latter can be achieved by choosing

[TABLE]

Hence by choosing $\lambda$ equal to $\lambda_{p}$, we get an efficient bound on $C_{{\bm{\gamma}},{\bm{w}},\lambda_{p}}=A_{\lambda_{p}}^{1+1/\lambda_{p}}$, as long as the $w_{j}$ are chosen to guarantee convergence of $\sum_{j=1}^{\infty}\left(b_{j}b^{w_{j}}\right)^{p}$.

4 Combined error bound

The derivation in the previous section leads to the following result.

Theorem 6**.**

Given the PDE in (1.1) for which we characterized the regularity of the random field by a sequence of $b_{j}$ with sparsity $p\in(0,1)$ and determined a sequence of $w_{j}$ such that $\sum_{j=1}^{\infty}(b_{j}\,b^{w_{j}})^{p}<\infty$, we can construct the generating vector for an $N$-point randomized lattice rule using the reduced CBC algorithm (Algorithm 2), at the cost of $\mathcal{O}(\sum_{j=1}^{\min\{s,s^{\ast}\}}(m-w_{j}+j)\,b^{m-w_{j}})$ operations, such that, assuming that (1.11), (2.3) and $\frac{\kappa\,a_{0,\max}}{(1-\overline{\kappa})\,a_{0,\min}}\sup_{j\geq s+1}b_{j}<1$ hold, we obtain an upper bound

[TABLE]

where the implied constant is independent of $s$, $h$ and $N$.

Observe that if the $w_{j}$ increase sufficiently fast, the construction cost of Algorithm 2 does not depend anymore on the increasing dimensionality. Further note that the first term on the right-hand side of (4.1) is small if $\sup_{j\geq s+1}b_{j}$ is small, and, since we assumed that $b_{j}$ must tend to zero by assumption (3.6), we can shrink the first summand by choosing $s$ sufficiently large. By choosing $h$ sufficiently small, and $N$ sufficiently large, we can also make the other two summands in the overall error bound small.

Note that $\sup_{j\geq s+1}b_{j}\leq\sum_{j\geq s+1}b_{j}$, and that (3.6) yields $\sum_{j=1}^{\infty}b_{j}^{p}<\infty$, which implies that one can use the machinery developed in [12] to obtain a cost analysis similar to [12, Theorem 8.1]. Note, in particular, that it is sufficient to choose $N$ of order $\mathcal{O}(\varepsilon^{-\lambda_{p}/2})$, independently of $s$, to meet an error threshold of $\varepsilon$.

5 Derivation of the fast reduced CBC algorithm

Finally in this last section the derivation of the fast reduced CBC algorithm for POD weights in Algorithm 2 is given. For prime $b$ and $m\in{\mathbb{N}}$ let $N=b^{m}$. Consider a generating vector $\widetilde{{\bm{z}}}=(b^{w_{1}}z_{1},\ldots,b^{w_{d}}z_{d})$ with $z_{j}\in{\mathbb{Z}}_{b^{m-w_{j}}}^{\times}$ and integer $0\leq w_{j}\leq m$ for each $j=1,\ldots,d$. Furthermore, for an integer $0\leq w^{\prime}\leq m$, the squared worst-case error can be written as

[TABLE]

We note that this holds for any integer $0\leq w^{\prime}\leq m$ and, in particular, for $w=0$ this is the vector being used in the normal fast CBC algorithm. We now write the error in terms of the previous error, as is standard for CBC algorithms, by splitting the expression into subsets ${\mathfrak{u}}\subseteq\{1:d\}$ for which $d\not\in{\mathfrak{u}}$ and $d\in{\mathfrak{u}}$, to obtain

[TABLE]

Since the choice of $z_{d}\in{\mathbb{Z}}_{b^{m-w_{d}}}^{\times}$ is modulo $b^{m-w_{d}}$, we can make a judicious choice for splitting up $k=k^{\prime}+t\,b^{m-w_{d}}$ for which the effect of dimension $d$ (for a choice of $z_{d}$) is then constant for all $t\in{\mathbb{Z}}_{b^{w_{d}}}$. We obtain

[TABLE]

Then we observe that for all $0\leq w_{d-1}\leq w_{d}\leq m$, with $k^{\prime}\in{\mathbb{Z}}_{b^{m-w_{d}}}$, writing $t=t^{\prime}+t^{\prime\prime}\,b^{w_{d}-w_{d-1}}\in{\mathbb{Z}}_{b^{w_{d}}}$ with $t^{\prime}\in{\mathbb{Z}}_{b^{w_{d}-w_{d-1}}}$ and $t^{\prime\prime}\in{\mathbb{Z}}_{b^{w_{d-1}}}$, leads to

[TABLE]

where $k^{\prime\prime}=k^{\prime}+t^{\prime}\,b^{m-w_{d}}\in{\mathbb{Z}}_{b^{m-w_{d-1}}}$ as required for $q_{d-1,\ell,w_{d-1}}(k^{\prime\prime})$. Note that this is the property of the “fold and sum” operator as introduced in (2.8) and mentioned there. Using matrix-vector notation, we rewrite the expression in (5.1) for all $z_{d}\in{\mathbb{Z}}_{b^{m-w_{d}}}^{\times}$ as

[TABLE]

where ${\bm{e}}_{b^{m},d}^{2}\in{\mathbb{R}}^{\varphi(b^{m-w_{d}})}$ is the vector with components $e_{b^{m},d}^{2}(b^{w_{1}}z_{1},\ldots,b^{w_{d}}z_{d})$ for all $z_{d}\in{\mathbb{Z}}_{b^{m-w_{d}}}^{\times}$. After $z_{d}$ has been selected we can calculate (for $\ell=1,\ldots,d$)

[TABLE]

In Algorithm 2 the vectors ${\bm{q}}_{j,\ell,w_{j}}$ are denoted by just ${\bm{q}}_{j,\ell}$.