On the Numerical Approximation of the Karhunen-Lo\`{e}ve Expansion for Lognormal Random Fields

TL;DR

This paper provides a detailed error analysis for the numerical approximation of the Karhunen-Loève expansion applied to lognormal random fields, crucial for PDEs with random coefficients, balancing truncation, sampling, and discretization errors.

Contribution

It introduces new error estimates that optimize the approximation of the KL expansion for lognormal fields, integrating truncation, sampling, and numerical errors.

Findings

Derived error bounds for KL expansion approximation

Quantified the impact of truncation, sampling, and discretization errors

Validated theoretical results with numerical experiments in multiple stochastic dimensions

Abstract

The Karhunen-Lo\`{e}ve (KL) expansion is a popular method for approximating random fields by transforming an infinite-dimensional stochastic domain into a finite-dimensional parameter space. Its numerical approximation is of central importance to the study of PDEs with random coefficients. In this work, we analyze the approximation error of the Karhunen-Lo\`eve expansion for lognormal random fields. We derive error estimates that allow the optimal balancing of the truncation error of the expansion, the Quasi Monte-Carlo error for sampling in the stochastic domain and the numerical approximation error in the physical domain. The estimate is given in the number $M$ of terms maintained in the KL expansion, in the number of sampling points $N$, and in the discretization mesh size $h$ in the physical domain employed in the numerical solution of the eigenvalue problems during the expansion. The result is used to quantify the error in PDEs with random coefficients. We complete the theoretical analysis with numerical experiments in one and multiple stochastic dimensions.

Tables

Table 1. Table 1: The root mean square error between log ⁡ κ 𝜅 \log\kappa and log ⁡ κ M N , h superscript subscript 𝜅 𝑀 𝑁 ℎ \log\kappa_{M}^{N,h} for different numbers of truncation terms M 𝑀 M and different mesh sizes h ℎ h . Here, the number of sampling points is N = 1009 𝑁 1009 N=1009 and the dimension d ′ = 1 superscript 𝑑 ′ 1 d^{\prime}=1 with the optimal parameters being ( M , h ) := ( 5 , 1 / 101 ) assign 𝑀 ℎ 5 1 101 (M,h):=(5,1/101) and the corresponding root mean square error being 8.0493e-3.

$h\backslash M$	2	4	8
1/16	4.1217e-2	1.1933e-2	3.6194e-3
1/64	4.1216e-2	1.1931e-2	3.6060e-3
1/128	4.1216e-2	1.1931e-2	3.6059e-3
1/256	4.1216e-2	1.1931e-2	3.6059e-3

Table 2. Table 2: The root mean square error between log ⁡ κ 𝜅 \log\kappa and log ⁡ κ M N , h superscript subscript 𝜅 𝑀 𝑁 ℎ \log\kappa_{M}^{N,h} for different numbers of truncation terms M 𝑀 M and different mesh sizes h ℎ h . Here, the number of sampling points is N = 1009 𝑁 1009 N=1009 , dimension d ′ = 10 superscript 𝑑 ′ 10 d^{\prime}=10 and the shifted-averaged worse-case error is 3.0987e-3.

$h\backslash M$	2	4	8
1/16	6.0723e-3	1.8676e-4	1.7078e-4
1/64	7.2336e-3	1.0267e-4	1.0837e-5
1/128	6.2803e-3	7.9009e-5	2.7146e-6
1/256	6.1652e-3	6.6978e-5	2.7102e-6

Table 3. Table 3: The root mean square error between log ⁡ κ 𝜅 \log\kappa and log ⁡ κ M N , h superscript subscript 𝜅 𝑀 𝑁 ℎ \log\kappa_{M}^{N,h} for different numbers of truncation terms M 𝑀 M and different mesh sizes h ℎ h . Here, the number of sampling points N = 1009 𝑁 1009 N=1009 , dimension d ′ = 100 superscript 𝑑 ′ 100 d^{\prime}=100 and the shifted-averaged worse-case error is 3.1045e-3.

$h\backslash M$	2	4	8
1/16	2.8446e-3	2.8332e-3	2.8335e-3
1/64	3.7890e-4	2.9770e-4	2.9743e-4
1/128	2.8832e-4	7.8182e-5	7.8268e-5
1/256	2.2743e-4	1.9821e-5	1.9772e-5

Equations

$$\left\|{\kappa-\kappa_{M}^{N,h}}\right\|_{L^{p}(\Omega,L^{2}(D))}\lesssim M^{-\frac{s}{d}}+M^{\frac{s}{d}+\frac{3}{2}}h^{s}+\Big{(}\frac{M}{N}\Big{)}^{1/2}.$$

$$W^{m,p}(D)=\{u\in L^p(D):D^{\alpha }u\in L^p(D)\text{ for }0\le |\alpha |\le m\}.$$

$$\|u\|_{W^{m,p}(D)}=\left\{\begin{aligned} \Big{(}\sum\limits_{0\leq|\alpha|\leq m}\left\|D^{\alpha}u\right\|_{L^{p}(D)}^{p}\Big{)}^{\frac{1}{p}},&\text{ if }1\leq p<\infty,\\ \max\limits_{0\leq|\alpha|\leq m}\left\|D^{\alpha}u\right\|_{L^{\infty}(D)},&\text{ if }p=\infty.\end{aligned}\right.$$

$$\log \mathrm{d}\mathcal{P}(y)=\mathbf{\rho }\mathrm{d}y:=\prod _{j=1}^{d^{\prime }}\rho (y_j)\mathrm{d}(y_j)\text{ with }\rho (y):=\frac{1}{\sqrt{2\pi }}\exp (-\frac{y^2}{2}),$$

$$(\mathcal{S}v)(y)={\int }_D\log \kappa (y,x)v(x)\mathrm{d}x,$$

$$({\mathcal{S}}^{\ast }v)(x)={\int }_{\Omega }\log \kappa (y,x)v(y)\mathbf{\rho }\mathrm{d}y.$$

$$R(x,x^{\prime })={\int }_{\Omega }\log \kappa (y,x)\log \kappa (y,x^{\prime })\mathbf{\rho }\mathrm{d}y.$$

$$\mathcal{R}v(x)={\int }_{D}R(x,x^{\prime })v(x^{\prime })\mathrm{d}{x}^{\prime }={\int }_D{\int }_{\Omega }\log \kappa (y,x)\log \kappa (y,x^{\prime })v(x^{\prime })\mathbf{\rho }\mathrm{d}y\mathrm{d}{x}^{\prime }.$$

$${\psi }_n(y)=\frac{1}{\sqrt{{\lambda }_n}}{\int }_D\log \kappa (y,x){\varphi }_n(x)\mathrm{d}x.$$

$${\lambda }_n=\inf \{{\Vert \mathcal{R}-L\Vert }_{\mathcal{B}(L^2(D))}:L\in \mathfrak{F}(L^2(D)),\text{rank}(L)<n\}$$

$$\log \kappa (y,x)=\sum _{n=1}^{\infty }\sqrt{{\lambda }_n}{\varphi }_n(x){\psi }_n(y),$$

$${\lambda }_n$$

$$\displaystyle\Big{\|}{\sum\limits_{n>M}\sqrt{\lambda_{n}}\phi_{n}(x)\psi_{n}(y)}\Big{\|}_{L^{2}(\Omega\times D)}$$

$${\Vert {\varphi }_n\Vert }_{H^{\theta s}\left(D\right)}\le C(D,d,s)n^{\frac{\theta s}{d}}\text{ when n is sufficiently large.}$$

$${\partial }_x^{\mathbf{\alpha }}\log \kappa (y,x)=\sum _{n=1}^{\infty }\sqrt{{\lambda }_n}{\partial }_x^{\mathbf{\alpha }}{\varphi }_n(x){\psi }_n(y).$$

$$\Vert {\partial }_x^{\mathbf{\alpha }}\log \kappa {\Vert }_{L^2(\Omega \times D)}^2=\sum _{n=1}^{\infty }{\lambda }_n\Vert {\partial }_x^{\mathbf{\alpha }}{\varphi }_n{\Vert }_{L^2(D)}^2.$$

$$\Vert \log \kappa {\Vert }_{L^2(\Omega ,H^s(D))}^2=\sum _{n=1}^{\infty }{\lambda }_n{\Vert {\varphi }_n\Vert }_{H^s\left(D\right)}^2.$$

$$+\infty >\sum _{n=1}^{\infty }{\lambda }_n{\Vert {\varphi }_n\Vert }_{H^s\left(D\right)}^2$$

$$=C_{\text{ref}thm:truncationError}\sum _{n=1}^{\infty }{n}^{-1-ϵ}\cdot n^{-\frac{2s}{d}+ϵ}{\Vert {\varphi }_n\Vert }_{H^s\left(D\right)}^2$$

$${\Vert {\varphi }_n\Vert }_{H^s\left(D\right)}\le C(D,d,s)n^{\frac{s}{d}}.$$

$$\underset{\begin{array}{c}\{c_n(x){\}}_{n=1}^M\subset L^2(D)\\ \{c_n(x){\}}_{n=1}^M\text{ orthonormal}\end{array}}{\min }{\Vert \log \kappa (y,x)-\sum _{n=1}^M\left({\int }_D\log \kappa (y,x)c_n(x)\mathrm{d}x\right)c_n(x)\Vert }_{L^2\left(\Omega \times D\right)}=\sqrt{\sum _{n>M}{\lambda }_n}.$$

$$R_N(x,x^{\prime }):=\sum _{n=1}^N{\omega }_n\log \kappa (y_n,x)\log \kappa (y_n,x^{\prime }).$$

$${\psi }_n^N(y)=\frac{1}{\sqrt{{\lambda }_n^N}}{\int }_D\log \kappa (y,x){\varphi }_n^N(x)\mathrm{d}x.$$

$$V_h:=\{v\in H^1(D):v{|}_K\in P^k(K)\text{ for all }K\in {\mathcal{T}}_h\}.$$

$${\Vert v-I_{h}v\Vert }_{L^2\left(D\right)}$$

$$h^{d/2}{\Vert v-I_{h}v\Vert }_{C\left(D\right)}$$

$$\displaystyle I_{h}\Big{(}\lambda_{n}^{N,h}\phi_{n}^{N,h}-\mathcal{R}_{N}\phi_{n}^{N,h}\Big{)}=0.$$

$${\mathcal{R}}_{N,h}:=I_h{\mathcal{R}}_N{I}_h.$$

$$\log {\kappa }_M^{N,h}(y,x):=\sum _{n=1}^M\sqrt{{\lambda }_n^{N,h}}{\varphi }_n^{N,h}(x){\psi }_n(y).$$

$${\psi }_n^{N,h}(y):=\frac{1}{\sqrt{{\lambda }_n^{N,h}}}{h}^2\sum _{n=1}^N\sum _{K\in {\mathcal{T}}_h}\sum _{x_j\in I_K}\log \kappa (x_j,y_n){\varphi }_n^{N,h}(x_j)L_n(y).$$

$${\tilde{\psi }}_n^{N,h}(y):=\frac{1}{\sqrt{{\lambda }_n^{N,h}}}{\int }_D\log \kappa (y,x){\varphi }_n^{N,h}(x)\mathrm{d}x.$$

$${\Vert \log \kappa -\log {\kappa }_M^{N,h}\Vert }_{L^2\left(\Omega \times D\right)}$$

Full text

Keywords: Karhunen-Loève expansion, eigenvalue decay, approximation of bivariate functions, error estimates, lognormal random field.

1 Introduction

Partial differential equations (PDEs) with random coefficient have been widely employed to describe applications that are affected by a certain amount of uncertainty arising from imperfect/insufficient information about the problem, e.g., in the input data. The range of applications is broad and diverse and includes, e.g., oil field modelling, quantum mechanics and finance [9, 16, 23]. The dimension of the random coefficient can be huge or even infinite, which poses enormous computational challenge. To reduce its dimensionality, one can parameterize the random coefficient by means of the Karhunen-Loève (KL) expansion or the polynomial chaos (PC) expansion [15, 22], which greatly facilitates the subsequent numerical treatment, e.g., by the stochastic Galerkin method or the stochastic collocation method. Alternatively, one may expand the random field with respect to the hierarchical Faber basis or some wavelet type basis; see [5, 12] for details. In this paper, we will focus on the KL expansion, which is known to be optimal in the sense of the mean square error.

To formulate the problem, let $D\subset\mathbb{R}^{d}$ be an open bounded domain with a strong local Lipschitz boundary and let $(\widetilde{\Omega},\Sigma,\mathcal{P})$ be a complete separable probability space with $\sigma$-field $\Sigma\subset 2^{\widetilde{\Omega}}$ and probability measure $\mathcal{P}$. We will denote $\Omega:=(\widetilde{\Omega},\Sigma,\mathcal{P})$ for notational simplicity. Now, we consider a stochastic field $\kappa(y,x)\in L^{\infty}(\Omega,L^{2}(D))$ with its logarithm being a centered Gaussian field. The lognormal random field is frequently used in stochastic PDEs as a random diffusion coefficient.

In practical computation, its numerical approximation usually proceeds in three steps. In the first step, the centered random field $\log\kappa(y,x)$ is approximated by its $M$-term KL expansion for some $M\in\mathbb{N}_{+}$. The truncation error relies on the regularity of the bivariate function $\log\kappa(y,x)$ in the physical variable $x$, see [18] for details. In the second step, the covariance function $R(x,x^{\prime})$ of the centered Gaussian random field $\log\kappa(y,x)$ is approximated via a sampling method. By its very definition, the covariance function involves an integral over the stochastic domain $\Omega$, which is often of very high dimensional. For its approximation, various quadrature-type sampling methods, e.g., Monte-Carlo methods, (Quasi) Monte-Carlo (QMC) methods and sparse grids [10, 13, 14] can be applied, say with $N$ sampling points. These methods essentially require boundedness of the variation, the first or higher mixed derivatives of $\log\kappa(\cdot,x)$ for fixed $x\in D$ and then yield a corresponding order of convergence. In this paper, we focus on the QMC method, which has only a low regularity requirement on $\Omega$, namely that the first mixed derivative of $\log\kappa(\cdot,x)$ is bounded. The outcome of this second step is a function $R_{N}(x,x^{\prime})\in L^{2}(D\times D)$ that approximates the covariance function $R(x,x^{\prime})$. The associated self-adjoint operators are denoted as $\mathcal{R}_{N}$ and $\mathcal{R}$, respectively. Note that $\mathcal{R}_{N}$ is a finite rank operator with rank not larger than $N$. We shall prove in Proposition 3.1 that only the first $\lfloor N^{\frac{1}{2s/d+1}}\rfloor$ terms in the KL expansion of $\mathcal{R}_{N}$ are relevant to approximate the spectrum of $R(x,x^{\prime})$. Here, the nonnegative parameter $s$ denotes the regularity of the bivariate function $\log\kappa(y,x)$ in the physical variable $x$. This result implies that the number of KL truncation terms satisfies $M\leq\lfloor N^{\frac{1}{2s/d+1}}\rfloor$. The third step is to approximate the eigenvalue problem of the self-adjoint operator $\mathcal{R}_{N}$ by means of a conforming Galerkin finite element method (FEM) over a regular mesh with a mesh size $h$. Now, to estimate the error between $\kappa(y,x)$ and its numerical approximation $\kappa_{M}^{N,h}(y,x)$ with $M$ being the number of truncation terms, $N$ being the number of sampling points and $h$ being the mesh size, the eigenvalue approximation error is derived. Moreover, to balance the decay of the eigenvalues of the covariance kernel $R(x,x^{\prime})$ and the numerical approximation error, we need to take $h\ll N^{-1/s}$ in order to ensure convergence in the first place. Otherwise, no convergence rate is guaranteed when solving the eigenvalue problems numerically.

The main contribution of this work is threefold. First, we present the spectral analysis of the finite rank operator $\mathcal{R}_{N}$, which allows us to specify the number of truncation terms. Second, we recall the error rate of QMC quadrature and provide an error estimate of the numerical approximation to the eigenvalue problem associated with the operator $\mathcal{R}_{N}$ in terms of mesh size $h$. Third, we derive an error estimate of both, $\log\kappa-\log\kappa_{M}^{N,h}$ and $\kappa-\kappa_{M}^{N,h}$ in various norms. Our final estimates of $\|\kappa-\kappa_{M}^{N,h}\|_{L^{p}(\Omega,L^{2}(D))}$ and $\|\kappa-\kappa_{M}^{N,h}\|_{L^{p}(\Omega,C(D))}$ with $1\leq p<2$ are presented in Theorem 5.6. For example, we obtain the bound

[TABLE]

Moreover, we discuss the example of an elliptic PDE with lognormal random diffusion coefficient. There, using our previous results on the approximation of the lognormal random field, we can deduce bounds of the error between the solution $u$ of the PDE and its induced approximation $u_{M}^{N,h}$.

The remainder of the paper is organized as follows. We formulate in Section 2 the approximation of $\log\kappa$ by the KL expansion, explain the general sampling method and discuss the Galerkin approximation. In Section 3, we analyze the Quasi Monte-Carlo method to approximate the covariance kernel $R(x,x^{\prime})$, and derive a spectral estimate for $\mathcal{R}$ and $\mathcal{R}_{N}$ by means of the maximin principle and an eigenvalue decay estimate. In Section 4, we discuss the conforming Galerkin approximation of the eigenvalue problems of $\mathcal{R}_{N}$ and derive spectral estimates. The main error estimates between $\kappa$ and $\kappa_{M}^{N,h}$ in the $L^{p}(\Omega,L^{2}(D))$-norm and the $L^{p}(\Omega,C(D))$-norm with $p\in[1,2)$, respectively, are established in Section 5. Furthermore, we present an application of our results for an elliptic operator with lognormal random coefficients in Section 6. Two numerical tests are provided in Section 7 to verify our findings. Finally, we give some concluding remarks in Section 8.

2 Preliminaries

This section collects elementary facts on the KL expansion and its numerical approximation. To this end, the overall numerical approximation error is divided into three parts: the truncation error, the sampling error and the resulting approximation error of the eigenvalue problems.

We start with some notation. Let two Banach spaces $V_{1}$ and $V_{2}$ be given. Then, $\mathcal{B}(V_{1},V_{2})$ stands for the Banach space composed of all continuous linear operators from $V_{1}$ to $V_{2}$ and $\mathcal{B}(V_{1})$ stands for $\mathcal{B}(V_{1},V_{1})$. The set of nonnegative integers is denoted by $\mathbb{N}$. For any index $\alpha\in\mathbb{N}^{d}$, $|\alpha|$ is the sum of all components. The letters $M$, $N$ and $h$ are reserved for the truncation number of the KL modes, the number of sampling points and the mesh size. We write $A\lesssim B$ if $A\leq cB$ for some absolute constant $c$ which is independent of $M$, $N$ and $h$, and we likewise write $A\gtrsim B$. Moreover, for any $m\in\mathbb{N}$, $1\leq p\leq\infty$, we follow [1] and define the Sobolev space $W^{m,p}(D)$ by

[TABLE]

It is equipped with the norm

[TABLE]

The space $W_{0}^{m,p}(D)$ is the closure of $C^{\infty}_{0}(D)$ in $W^{m,p}(D)$. Its dual space is $W^{-m,q}(D)$, with ${1}/{p}+{1}/{q}=1$. Also we use $H^{m}(D)=W^{m,p}(D)$ for $p=2$. $(\cdot,\cdot)$ denotes the inner product in $L^{2}(D)$.

2.1 Karhunen-Loève expansion: continuous level

In this work, we consider a stochastic field $\kappa(y,x)\in L^{2}(\Omega\times D)$ with its logarithm being a centered Gaussian field, i.e.,

[TABLE]

where $\mathcal{P}$ is the probability measure on $\Omega$ introduced in Section 1. We denote the associated integral operator $\mathcal{S}:L^{2}(D)\rightarrow L^{2}(\Omega)$ by

[TABLE]

whereas its adjoint operator $\mathcal{S}^{*}:L^{2}(\Omega)\rightarrow L^{2}(D)$ is defined by

[TABLE]

Let $\mathcal{R}:L^{2}(D)\rightarrow L^{2}(D)$ be defined by $\mathcal{R}:=\mathcal{S}^{*}\mathcal{S}$. Then $\mathcal{R}$ is a nonnegative self-adjoint Hilbert-Schmidt operator with kernel $R\in L^{2}(D\times D):D\times D\to\mathbb{R}$ given by

[TABLE]

This is just the covariance function of the stochastic process $\log\kappa(x,{y})$. Moreover, for any $v\in L^{2}(D)$, we have

[TABLE]

The standard spectral theory for compact operators [25] implies that the operator $\mathcal{R}$ has at most countably many discrete eigenvalues, with zero being the only accumulation point, and each non-zero eigenvalue has only finite multiplicity. Let $\{\lambda_{n}\}_{n=1}^{\infty}$ be the sequence of eigenvalues (with multiplicity counted) associated to $\mathcal{R}$, which are ordered nonincreasingly, and let $\{\phi_{n}\}_{n=1}^{\infty}$ be the corresponding eigenfunctions that are orthonormal in $L^{2}(D)$. Furthermore, for any $\lambda_{n}\neq 0$, define

[TABLE]

One can verify that the sequence $\{\psi_{n}\}_{n=1}^{\infty}$ is uncorrelated and orthonormal in $L^{2}(\Omega)$, and therefore, $\{\psi_{n}\}_{n=1}^{\infty}$ are i.i.d normal random functions.

Note that the sequence $\{\lambda_{n}\}_{n=1}^{\infty}$ can be characterized by the so-called approximation numbers (cf. [21, Section 2.3.1]). They are defined by

[TABLE]

where $\mathfrak{F}(L^{2}(D))$ denotes the set of the finite rank operators on $L^{2}(D)$. This equivalency is frequently employed to estimate eigenvalues by constructing finite rank approximation operators to $\mathcal{R}$.

The KL expansion of the bivariate function $\log\kappa(y,x)$ then refers to the expression

[TABLE]

where the series converges in $L^{2}(\Omega\times D)$.

2.2 Karhunen-Loève expansion: $M$-term truncation

Now, we will truncate the KL expansion and discuss the resulting error. The studies on the $M$-term KL approximation to random fields are extensive. In [22], the authors derived the eigenvalue decay rates for random fields with their corresponding covariance kernels possessing certain regularity and considered the generalized fast multipole methods to solve the associated eigenvalue problems. Robust eigenvalue computation for smooth covariance kernels was studied in [24]. A comparison of $M$-term KL truncation and the sparse grids approximation was given in [17].

The result of this section is based on our recent paper [18], which proves a sharp eigenvalue decay rate under a mild assumption on the regularity of the bivariate function $\log\kappa(y,x)$ in the physical domain. To this end, we make the following assumption.

Assumption 2.1** (Regularity of $\log\kappa(y,x)$).**

There exists some $s\geq 0$ such that $\log\kappa(y,x)\in L^{\infty}(\Omega,H^{s}(D))$.

Under Assumption 2.1, by the definition of the kernel $R(x,x^{\prime})$, we have $R(x,x^{\prime})\in H^{s}(D)\times H^{s}(D)$.

The following eigenvalue decay estimate [18, Theorems 3.2, 3.3 and 3.4] will be used repeatedly.

Theorem 2.1**.**

Let Assumption 2.1 hold. Then, for any $M\in\mathbb{N}$ sufficiently large, there holds

[TABLE]

with the constant $C_{\ref{thm:truncationError}}:=\mathrm{diam}(D)^{2s}C_{\rm em}(d,s)C_{\rm{ext}}(D,s)\left\|\log\kappa\right\|_{L^{2}(\Omega,H^{s}\left(D\right))}^{2}$. Here, $C_{\rm{em}}(d,s)$ denotes the embedding constant between the Lorentz sequence spaces $\ell_{\frac{d}{d+2s},1}\hookrightarrow\ell_{\frac{d}{d+2s},\infty}$ and $C_{\rm{ext}}(D,s)$ is a constant depending only on $D$ and $s$.

The next lemma gives the regularity of the eigenfunctions $\{\phi_{n}\}_{n=1}^{\infty}$.

Lemma 2.1** (Regularity of the eigenfunctions $\{\phi_{n}\}_{n=1}^{\infty}$).**

Let Assumption 2.1 be valid. Then for all $0\leq\theta\leq 1$, there holds

[TABLE]

Here, $C(D,d,s)$ denotes a positive constant depending only on $D$, $d$ and $s$.

Proof.

We will only prove the result for $s\in\mathbb{N}_{+}$. The case for $s\in\mathbb{R}^{+}$ can be obtained by the interpolation method.

Let $\boldsymbol{\alpha}=[\alpha_{1},\cdots,\alpha_{d}]\in\mathbb{N}^{d}$ with $|\boldsymbol{\alpha}|:=\sum\limits_{i=1}^{d}\alpha_{i}\leq s$. The combination of Assumption 2.1 and decomposition (2.5) and an application of Lebesgue’s dominated convergence theorem lead to the expansion

[TABLE]

After taking the squared $L^{2}(\Omega\times D)$-norm on both sides, we arrive at

[TABLE]

Now we sum over all $\boldsymbol{\alpha}\in\mathbb{N}^{d}$ with $|\boldsymbol{\alpha}|\leq s$, and obtain by the definition of the Sobolev space $H^{s}(D)$ that

[TABLE]

At last, an application of Theorem 2.1 gives

[TABLE]

for any positive parameter $\epsilon$. With $0<\epsilon\to 0$ we obtain from the relation (2.7) that, when $n$ is sufficiently large, there holds

[TABLE]

This verifies (2.6) for $\theta=1$. By noting that $\left\|\phi_{n}\right\|_{L^{2}\left(D\right)}=1$, an application of [2, Theorem 3.3] yields then the desired estimate. ∎

It is worth to emphasize the optimality of the eigenfunctions $\{\phi_{n}\}_{n=1}^{\infty}$ in the sense that the mean-square error resulting from a finite-rank approximation of $\kappa(y,x)$ is minimized [15]. Thus, the eigenfunctions indeed minimize the truncation error in the $L^{2}$-sense, i.e.

[TABLE]

2.3 Sampling estimate of the continuous Karhunen-Loève approximation

Clearly, any numerical computation of the covariance function $R(x,x^{\prime})$ by a conventional quadrature method quickly becomes expensive and impractical when the dimensionality $d^{\prime}$ of the random domain $\Omega$ is large. This is due to the curse of dimensionality. To this end, depending on the regularity prerequisites with respect to the stochastic variable $y$, the Monte Carlo method, the Quasi-Monte Carlo (QMC) methods or the sparse grid method may be employed in approximating $R(x,x^{\prime})$. In this paper, we will focus on QMC.

Anyway, a numerical quadrature gives $R_{N}(x,x^{\prime})$, which is defined by

[TABLE]

Here, $N\in\mathbb{N}$ denotes the number of quadrature points and $\{{y}_{1},\cdots,{y}_{N}\}$ and $\{\omega_{1},\cdots,\omega_{N}\}$ are the corresponding quadrature points and weights. Clearly, $R_{N}\in L^{2}(D\times D)$ and $R_{N}:D\times D\to\mathbb{R}$.

Analogously, we denote by $\mathcal{R}_{N}$ the nonnegative self-adjoint Hilbert-Schmidt operator with kernel $R_{N}$. The operator $\mathcal{R}_{N}$ is of rank no greater than $N$ and hence compact. Analogously, we can define in nondecreasing order its eigenvalues and its normalized eigenfunctions in $L^{2}(D)$ as $\{\lambda_{n}^{N}\}_{n=1}^{N}$ and $\{\phi_{n}^{N}\}_{n=1}^{N}$, respectively.

Note at this point the following: If we are interested in a specific approximate realization of $\log\kappa(y,\cdot)$ for some $y\in\Omega$, then we have to consider the function $\psi_{n}^{N}(y)$ defined by

[TABLE]

To estimate the error between $\psi_{n}$ and $\psi_{n}^{N}$, we can apply finite elements $\mathcal{T}_{h}$ over $D$ as introduced in Subsection 2.4. This error depends on the regularity of $\log\kappa(y,\cdot)$ for given $y\in\Omega$. On the other hand, if we are only interested in certain statistical quantities of the Gaussian random field $\log\kappa$, then there is no need to calculate $\{\psi_{n}^{N}\}_{n=1}^{N}$, and we can take directly i.i.d normal random functions, e.g., $\{\psi_{n}\}_{n=1}^{N}$. This is indeed the situation many articles are concerned with, see e.g., [5, 11, 19].

2.4 Galerkin discretization of the sampled, truncated continuous Karhunen-Loève approximation

Now we describe the conforming Galerkin approximation of the eigenvalue problem on $\mathcal{R}_{N}$. To this end, let $\mathcal{T}_{h}$ be a regular quasi-uniform triangulation over the physical domain $D$ with a maximal mesh size $h$ and let $k:=\lceil s\rceil$. The associated finite element space $V_{h}$ is defined by

[TABLE]

Let $Q$ be the dimension of $V_{h}$. We then have Q=\mathcal{O}\Big{(}(\frac{h}{k})^{d}\Big{)}. The $L^{2}$-projection $I_{h}:L^{2}(D)\to V_{h}$ has the approximation property [8, Theorem 4.4.20]

[TABLE]

Here, the positive constant $C_{I_{h}}$ depends only on the regularity parameter of $\mathcal{T}_{h}$ and is independent of the mesh size $h$.

The conforming Galerkin approximation of the eigenvalue problem of $\mathcal{R}_{N}$ is to find $\{\lambda_{n}^{N,h},\phi_{n}^{N,h}\}_{n=1}^{Q}\subset\mathbb{R}\times V_{h}$ such that

[TABLE]

This is equivalent to the eigenvalue problem of the finite-rank operator on $L^{2}(D)$ defined by

[TABLE]

Let $\{\lambda_{n}^{N,h},\phi_{n}^{N,h}\}_{n=1}^{Q}$ be the corresponding eigenpairs with eigenvalues in nonincreasing order and eigenvectors orthonormal in $L^{2}(D)$. Then the $M$-term truncated KL expansion, denoted by $\kappa_{M}^{N,h}(y,x)$, is defined by

[TABLE]

Note at this point the following: Again, if we are mainly concerned with the approximation to a specific bivariate function $\log\kappa$ via the expression (2.14), then we have to replace $\psi_{n}$ with its numerical approximation

[TABLE]

Here, $I_{K}$ represents the quadrature points on each finite element $K\in\mathcal{T}_{h}$ and $\{L_{n}(y)\}_{n=1}^{N}$ denotes the Legendre polynomials of order $N$. Note that ${\psi}_{n}^{N,h}(y)$ is the numerical approximation by interpolation with sampling points $\{y_{n}\}_{n=1}^{N}$ to $\tilde{\psi}_{n}^{N,h}(y)$ defined as

[TABLE]

In view of the KL expansion (2.5) and the $M$-term truncation estimate (2.14), an application of the triangle inequality yields

[TABLE]

A main goal of this paper is to derive a sharp estimate of $\left\|\log\kappa-\log\kappa_{M}^{N,h}\right\|_{L^{2}\left(\Omega\times D\right)}$ in (2.17). To this end, it suffices to analyze the three terms on the right hand side of (2.17). Here, the first term represents the truncation error that can be estimated by Theorem 2.1, the second term is due to sampling of the KL approximation and the third term is induced by the Galerkin approximation error.

3 QMC method approximation error

In this section, we apply the QMC method based on the randomly shifted lattice rule and derive the sampling error corresponding to the second term in (2.17). To this end, we map the quadrature points $\Xi_{N}:=\{\xi_{1},\xi_{2},\cdots,\xi_{N}\}\subset[0,1]^{d^{\prime}}$ to $\mathbb{R}^{d^{\prime}}$ by the inverse of the cumulative distribution function of the standard normal distribution. The cumulative distribution function $\boldsymbol{\phi}(y)$ is defined by

[TABLE]

and its inverse is $\boldsymbol{\phi}^{-1}({y}):(0,1)^{d^{\prime}}\to\mathbb{R}^{d^{\prime}}.$ Upon changing variables, we obtain

[TABLE]

Then by taking ${y}_{i}:=\boldsymbol{\phi}^{-1}(\xi_{i})$ and $\omega_{i}:=\frac{1}{N}$ for $i=1,2,\cdots,N$ in (2.9), we get an approximation to the covariance function $R(x,x^{\prime})$, which is denoted by $R_{N}(x,x^{\prime})$.

To this end, we introduce the construction of the quadrature points $\Xi_{N}$, which is based on the fast CBC construction of randomly shifted lattice rules in the unanchored space [20]. The unanchored space $\mathcal{F}(\mathbb{R}^{d^{\prime}})$ is defined by

[TABLE]

Here, the positive function $\nu$ controls the boundary behavior of the functions in $\mathcal{F}(\mathbb{R}^{d^{\prime}{}})$. The collection of parameters $\gamma_{\boldsymbol{\alpha}}$ for all $\boldsymbol{\alpha}\subset\{1,\cdots,d^{\prime}{}\}$ controls the relative importance of various groups of variables, and

[TABLE]

Note that we will choose the weight function $\nu$ and the weight parameters $\gamma_{\boldsymbol{\alpha}}$, such that the bivariate function $\log\kappa(\cdot,x)$ belongs to $\mathcal{F}(\mathbb{R}^{d^{\prime}})$ for all $x\in D$.

We apply the CBC approach [20, Algorithm 6] to derive the generating vector $\boldsymbol{z}\in[0,1)^{d^{\prime}}$ with the number of sampling points being $N$. To this end, let the shift $\boldsymbol{\Delta}\in[0,1]^{d^{\prime}}$ be an i.i.d uniformly distributed vector. Then we obtain the randomly shifted (rank-1) lattice rule by

[TABLE]

Now, $R_{N}(x,x^{\prime})$ in (2.9) can be approximated by taking $y_{i}:=\boldsymbol{\phi}^{-1}(\xi_{i})$ for $i=1,2,\cdots,N$.

The error $e_{d^{\prime},N}(\boldsymbol{z})$ between $R(x,x^{\prime})$ and $R_{N}(x,x^{\prime})$ is measured by the shifted-averaged worse-case error defined by

[TABLE]

Thus, using the CBC Algorithm to calculate $R_{N}(x,x^{\prime})\in L^{2}(D\times D)$ defined in (2.9) yields a shifted-averaged worse-case error of $\mathcal{O}(N^{-1+\delta})$ for any $\delta>0$ with the construction cost of $\mathcal{O}(d^{\prime}N\log(N))$. Therefore, we start with the following setting.

Assumption 3.1** (Assumption on the sampling error).**

For some $\delta\in(0,1)$, there holds

[TABLE]

To approximate a bivariate function or a specific realization of the random field $\log\kappa(y,\cdot)$, we have introduced in the last section the quantities $\{\psi_{n}^{N}\}_{n=1}^{\infty}$, cf. (2.10), which are not orthonormal in $L^{2}(\Omega)$. Nevertheless, they are very close to an orthonormal basis when the approximation error between $\mathcal{R}$ and $\mathcal{R}_{N}$ is very small.

Lemma 3.1** (Near orthonormality of $\{\psi^{N}_{n}\}_{n=1}^{N}$).**

Let ${\psi}^{N}_{n}$ be defined as in (2.10). There holds

[TABLE]

Proof.

This is a direct consequence of the definition (2.10) and the eigenvalue problem for $\mathcal{R}_{N}$. ∎

Next, we give some estimates on the finite-rank approximation $\mathcal{R}_{N}$ and its spectrum.

Proposition 3.1** (Spectral estimate for $\mathcal{R}_{N}$).**

Let Assumption 2.1 hold, let $N\in\mathbb{N}_{+}$ be sufficiently large and let $M:=\lfloor N^{\frac{1}{2s/d+1}}\rfloor$. Then $\mathcal{R}_{N}\in\mathcal{B}(L^{2}(D),H^{s}(D))$ with

[TABLE]

For $1\leq n\leq N$, there holds

[TABLE]

Furthermore, let $\lambda_{k_{i}}$ be an eigenvalue of $\mathcal{R}$ with multiplicity $q_{i}$ for $i=1,2,\cdots$ and $k_{I-1}<N\leq k_{I}$ for some $I\in\mathbb{N}_{+}$. Assume that for sufficiently large $N$, there holds

[TABLE]

Then, for $1\leq n\leq N$, there holds

[TABLE]

In addition,

[TABLE]

Proof.

We can obtain from the definition (2.9) and the triangle inequality

[TABLE]

Then Assumption 2.1 leads to (3.5). The relation (3.6) is derived from the definition. To prove (3.8), fix $1\leq n\leq N$ and let $V_{n}=\text{span}\{\phi_{1},\cdots,\phi_{n}\}$ be a $n$-dimension subspace. Since $\mathcal{R}$ and $\mathcal{R}_{N}$ are nonnegative and self-adjoint, we obtain

[TABLE]

Next we estimate the lower bound of the minimum on the right hand side of (3.10). To this end, note that any $v\in V_{n}$ admits the expression $v=\sum\limits_{i=1}^{n}c_{i}\phi_{i}$ for some $\{c_{i}\}_{i=1}^{n}\subset\mathbb{R}^{n}$. Let $(v,v):=1$, then $\sum\limits_{i=1}^{n}c_{i}^{2}=1$. For any $\delta>0$, plugging in the expression for $v$ and applying (3.4) lead to

[TABLE]

The lower bound of the minimum can now be estimated using Lagrange multipliers. To this end, let $\mu\in\mathbb{R}$ and define

[TABLE]

Let $(\boldsymbol{c^{*}},\mu^{*})=(c_{1}^{*},\cdots,c_{n}^{*},\mu^{*})$ be the optimal point to the unconstrained minimization problem associated to $F(c_{1},\cdots,c_{n};\mu)$. Then $c_{1}^{*},\cdots,c_{n}^{*}$ have the same sign by the definition of $f$. Let $c_{i}^{*}\geq 0$ for all $1\leq i\leq n$. The optimality conditions read $\frac{\partial F}{\partial c_{i}}|_{(\boldsymbol{c^{*}},\mu^{*})}=0$ for $i=1,\cdots,n$, and $\frac{\partial F}{\partial\mu}|_{(\boldsymbol{c^{*}},\mu^{*})}=0$. This immediately implies

[TABLE]

The second relation in (3.12) implies for all $1\leq i\leq n$ that there holds

[TABLE]

Summing over $i=1,\cdots,n$ yields

[TABLE]

Recall that $c_{i}^{*}\geq 0$ for all $i=1,\cdots,n$. Together with (3.12), this implies $\mu^{*}<\lambda_{n}.$ Now combining (3.7) and (3.13) results in $\frac{1}{\lambda_{n}-\mu^{*}}\gtrsim N,$ and therefore, $\mu^{*}\geq\lambda_{n}-C_{1}N^{-1}$ for some positive constant $C_{1}$ independent of $N$. This, together with (3.12), (3.11) and (3.10), gives $\lambda_{n}-\lambda_{n}^{N}\lesssim N^{-1}.$ Analogously, by changing the roles of $\mathcal{R}$ and $\mathcal{R}_{N}$, we can show

[TABLE]

Consequently, (3.8) follows by Theorem 2.1.

It remains to prove (3.9). We only present the proof for $n=1$. For $n>1$, (3.9) can be shown similarly to [4, Theorem 9.1]. Since the whole space $L^{2}(D)$ is orthogonally decomposed as the direct sum of the range of $\mathcal{R}$ and its kernel, $\phi_{1}^{N}$ can be split into

[TABLE]

for some $v\in L^{2}(D)$ satisfying $\mathcal{R}v=0$. Recall that $\lambda_{1}=\cdots\lambda_{q_{1}}>\lambda_{q_{1}+1}$. This leads to

[TABLE]

which, combined with (3.4) and (3.14), gives

[TABLE]

By redefining $\phi_{1}$ to be \big{(}\sum\limits_{i=1}^{q_{1}}c_{i}^{2}\big{)}^{-\frac{1}{2}}\sum\limits_{i=1}^{q_{1}}c_{i}\phi_{i}, (3.9) is proved due to the spectral gap assumption (3.7). ∎

Remark 3.1**.**

If the number of sampling points $N$ is not sufficiently large then the spectral gap assumption (3.7) is not fulfilled. Then one can show that for $1\leq n\leq N$, there holds

[TABLE]

Thus, to make $\lambda_{n}^{N}$ an accurate approximation to $\lambda_{n}$ for $1\leq n\leq M$, we have to impose a much more stringent restriction on $M:=\lfloor N^{\frac{1}{2s/d+2}}\rfloor$. In this manner, we can show that $\frac{n}{N}\ll n^{-\frac{2s}{d}-1}$ for $1\leq n\leq M$.

Finally, we can give an estimate to the second term in (2.17):

Proposition 3.2** (Root mean square error to the second term in (2.17)).**

Let $N$ be sufficiently large and $M\leq\lfloor N^{\frac{1}{2s/d+1}}\rfloor$. Furthermore, let (3.7) be satisfied. Then there holds

[TABLE]

Proof.

Employing the triangle inequality yields

[TABLE]

Then the desired result follows from (3.14) and (3.9). ∎

Next we give an estimate on $\|\psi_{n}-\psi_{n}^{N}\|_{L^{2}(\Omega)}$.

Lemma 3.2** (Estimate on $\|\psi_{n}-\psi_{n}^{N}\|_{L^{2}(\Omega)}$).**

Let $N$ be sufficiently large and $M\leq\lfloor N^{\frac{1}{2s/d+1}}\rfloor$. Furthermore, let (3.7) be satisfied. Then there holds

[TABLE]

Proof.

By (2.3), we obtain

[TABLE]

Therefore, taking the $L^{2}(\Omega)$-norm on both sides yields

[TABLE]

where, in the last inequality, we have applied (3.9) and the inequality $|\sqrt{a}-\sqrt{b}|\leq\sqrt{|a-b|}$ for all $a,b\geq 0$. ∎

In order to numerically approximate each realization of the Gaussian random field $\log\kappa(y,\cdot)$ for given $y\in\Omega$, we need the following result.

Proposition 3.3**.**

Let $N$ be sufficiently large and $M\leq\lfloor N^{\frac{1}{2s/d+1}}\rfloor$. Furthermore, let (3.7) be satisfied. Then

[TABLE]

Proof.

The triangle inequality yields

[TABLE]

Then the desired result follows from (3.14), (3.9) and Lemma 3.2. ∎

4 Conforming Galerkin approximation estimate

In this section we derive an estimate for the third term in (2.17) by means of the approximation theory of conforming finite element methods. To this end, let

[TABLE]

Then $\mathcal{E}_{N,h}$ is a self-adjoint operator on $L^{2}(D)$ and we have the following error representation.

Lemma 4.1**.**

The error operator $\mathcal{E}_{N,h}$ has the property

[TABLE]

Proof.

For given $v\in L^{2}(D)$ and since $\mathcal{R}_{N,h}=I_{h}\mathcal{R}_{N}I_{h}$ and $(\mathcal{R}_{N}I_{h}v-I_{h}\mathcal{R}_{N}v,v)=0$, we obtain

[TABLE]

∎

A direct consequence of Lemma 4.1, together with the approximation property (2.12) and Proposition 3.1, is the upper bound estimate for the operator norm of $\mathcal{E}_{N,h}$

[TABLE]

Finally, we are ready to present the main result in this section.

Proposition 4.1** (Conforming Galerkin approximation estimate).**

Let Assumption 2.1 hold and let $N\in\mathbb{N}_{+}$ be sufficiently large and $M\leq\lfloor N^{\frac{1}{2s/d+1}}\rfloor$. Assume that the spectral gap condition (3.7) is valid. Then there are constants $C_{1}$, $C_{2}$ and $h_{0}\ll N^{-1/s}$ such that

[TABLE]

Furthermore, the eigenvectors $\{\phi^{N}_{n}\}_{n=1}^{N}$ can be selected such that

[TABLE]

Here, the constants $C_{1}$ and $C_{2}$ are independent of $h$ and $N$ and $h_{0}>0$ is sufficiently small.

Proof.

The proof follows from [4, Theorem 9.1], where the following identity plays a crucial role. We have

[TABLE]

for all $v\in L^{2}(D)$ satisfying $(v,v)=1$. This identity can be derived by definition directly. Together with Proposition 3.1 and the estimate (4.2), this completes the proof. ∎

Recall that if we want to approximate a certain realization of the random field $\log\kappa(y,\cdot)$ for some $y\in\Omega$, then we have to estimate the error between ${\psi}_{n}^{N}$ and $\psi_{n}^{N,h}$. To this end, we make the following assumption.

Assumption 4.1**.**

Let ${\psi}_{n}^{N,h}$ and $\tilde{\psi}_{n}^{N,h}$ be as defined in (2.15) and (2.16), respectively. Assume that, for some $N_{0}$ sufficiently large and a sufficiently small $h_{0}\ll N_{0}^{-1/s}$, there holds

[TABLE]

Note that Assumption 4.1 requires a certain regularity of the bivariate function $\log\kappa(y,x)\in L^{2}(\Omega\times D)$ over the stochastic domain $\Omega$ since the computation of $\tilde{\psi}_{n}^{N,h}$ involves the polynomial interpolation of $\psi_{n}^{N,h}$ over $\Omega$.

We then get the following result.

Lemma 4.2**.**

Let Assumption 4.1 be valid. Let ${\psi}_{n}^{N}$ and $\psi_{n}^{N,h}$ be as defined in (2.10) and (2.15), respectively. Then, for some $N_{0}$ sufficiently large and a sufficiently small $h_{0}\ll N_{0}^{-1/s}$, there holds

[TABLE]

Proof.

An application of the triangle inequality leads to

[TABLE]

The second term can be estimated by Assumption 4.1. The definitions (2.3) and (2.16) imply

[TABLE]

Taking the $L^{2}(\Omega)$-norm on both sides leads to

[TABLE]

Then the desired result follows from Lemma 3.1 and Proposition 4.1. ∎

5 Main estimates

In this section, we present the main estimate of the error between the lognormal random field $\kappa$ and its $M$-term numerical approximation $\kappa_{M}^{N,h}$ in the $L^{p}(\Omega,L^{2}(D))$-norm and in the $L^{p}(\Omega,C(D))$-norm for $p<2$. The overall procedure is as follows: We first derive in Sections 5.1 and 5.2 an estimate on $|\log\kappa-\log\kappa_{M}^{N,h}|$ with respect to the $L^{2}(\Omega\times D)$-norm and the $L^{2}(\Omega,C(D))$-norm, cf. Theorems 5.1 and 5.4. Then we establish the final results on $|\kappa-\kappa_{M}^{N,h}|$ in Theorem 5.6 by employing Fernique’s Theorem.

5.1 $L^{2}$ error estimate

First, we give an estimate for the third term in (2.17).

Proposition 5.1** (Galerkin approximation estimate in (2.17)).**

Let Assumption 2.1 hold, let $N\in\mathbb{N}_{+}$ be sufficiently large and let $M\leq\lfloor N^{\frac{1}{2s/d+1}}\rfloor$ and $h\ll N^{-1/s}$. Assume the spectral gap condition (3.7) to be valid. Then there holds

[TABLE]

Proof.

Due to the orthogonality of the basis functions $\{\psi_{n}\}_{n=1}^{\infty}$ in $L^{2}(\Omega)$, an application of the triangle inequality leads to

[TABLE]

Here, we have used the orthogonality of $\{\phi_{n}^{N}\}_{n=1}^{M}$ over $L^{2}(D)$ in the last step. Then, an application of Proposition 4.1 and Theorem 2.1 gives the desired result. ∎

Now, using Theorem 2.1, Propositions 3.2 and 5.1, we are finally ready to present an estimate for (2.17).

Theorem 5.1** (Root mean square error for $M$-term KL expansion).**

Let Assumption 2.1 hold. Let $N\in\mathbb{N}_{+}$ be large, let $M\leq\lfloor N^{\frac{1}{2s/d+1}}\rfloor$ and $h\ll N^{-1/s}$. Assume the spectral gap condition (3.7) to be valid. Let $\log\kappa_{M}^{N,h}$ be given as in (2.14). Then there holds

[TABLE]

Remark 5.1** (Complexity).**

According to Theorem 5.1, in order to approximate the Gaussian random field $\log\kappa$ by formula (2.14) with root mean square error of $\epsilon$ for a certain threshold accuracy $\epsilon>0$, we need to take the number of sampling points $N:=\mathcal{O}(\epsilon^{-2-d/s})$, the number of truncation terms $M:=\mathcal{O}(\epsilon^{-d/s})$ and the mesh size $h:=\mathcal{O}(\epsilon^{2/s+3d/2s^{2}})$.

To numerically approximate the realization of $\log\kappa(y,\cdot)$ for given $y\in\Omega$ by the $M$-term truncation formula (2.14), we can replace the i.i.d normal random functions $\{\psi_{n}(y)\}_{n=1}^{M}$ with $\{\psi_{n}^{N,h}(y)\}_{n=1}^{M}$ defined in (2.15). The error in this process can be estimated as follows.

Theorem 5.2** (Root mean square error for $M$-term KL expansion of a bivariate function).**

Let Assumptions 2.1 and 4.1 hold. Let $N\in\mathbb{N}_{+}$ be large, let $M\leq\lfloor N^{\frac{1}{2s/d+1}}\rfloor$ and $h\ll N^{-1/s}$. Assume the spectral gap condition (3.7) to be valid. Let

[TABLE]

Then there holds

[TABLE]

Proof.

An application of the triangle inequality together with the KL expansion (2.5) implies

[TABLE]

Now, the first term and the third term above can be bounded by Proposition 3.3 and Theorem 2.1, respectively. We only need to estimate the second term. The triangle inequality gives

[TABLE]

By Lemma 3.1 and (3.4), $\|\psi_{n}^{N}\|_{L^{2}(\Omega)}\lesssim 1$ for all $1\leq n\leq M.$ Then, Proposition 4.1, Lemma 4.2 and Theorem 2.1 together show the desired result. ∎

Remark 5.2** (Complexity).**

According to Theorem 5.2, in order to approximate a specific realization of the Gaussian random field $\log\kappa$ by formula (5.1) with root mean square error for a certain threshold accuracy $\epsilon>0$, we need to choose the number of sampling points $N:=\mathcal{O}(\epsilon^{-2-2d/s})$, the number of truncation terms $M:=\mathcal{O}(\epsilon^{-d/s})$ and the mesh size $h:=\mathcal{O}(\epsilon^{3/s+3d/2s^{2}})$.

5.2 Uniform error estimate

In order to derive a uniform error estimate of the Gaussian random field $\log\kappa(y,x)$, we require a further regularity assumption on $\log\kappa$ to guarantee that $\log\kappa\in L^{2}(\Omega,C(D))$. To this end, we make the following assumption.

Assumption 5.1** (Regularity of $\log\kappa(y,x)$).**

Let Assumption 2.1 hold. Furthermore, assume $s>d/2$.

Then, the following estimate is valid.

Theorem 5.3** (Uniform estimate on the eigenfunctions).**

Let Assumption 5.1 be satisfied. Then there holds

[TABLE]

Proof.

Due to Assumption 5.1, an application of (2.6) with $\theta s>d/2$ together with the Sobolev embedding implies the desired result. ∎

Remark 5.3** (Optimality of the uniform estimate (5.2)).**

In [6, Section 3] the uniform estimate of eigenfunctions was studied under certain assumptions on the stationary covariance kernel. There, the authors derived a similar uniform estimate as (5.2) by utilizing essentially the regularity of the Fourier transform of the covariance kernel. They showed the sharpness of their uniform estimate in the case when $D=[0,1]$ and for the stationary covariance kernel $R(x,x^{\prime})=R(x-x^{\prime})$ with its Fourier transform $\hat{R}=\chi_{[-F,F]}$. Then, the uniform estimate of the n-th eigenfunction is $\mathcal{O}(n^{1/2})$, see [7].

Proposition 5.2** (Uniform truncation estimate).**

Let Assumption 5.1 be satisfied. Then, for any $1\leq M\in\mathbb{N}$, there holds

[TABLE]

Proof.

By the Sobolev embedding theorem, Assumption 5.1 implies $R(x,x^{\prime})\in C(D\times D)$. Then one can obtain

[TABLE]

This and Theorem 2.1 yield the desired estimate. ∎

Proposition 5.3**.**

Let Assumption 5.1 be fulfilled. Then $\mathcal{R}$ and $\mathcal{R}_{N}$ are bounded from $L^{2}(D)$ to $C(D)$, i.e., $\mathcal{R}$ and $\mathcal{R}_{N}\in\mathcal{B}(L^{2}(D),C(D))$. In addition, it holds

[TABLE]

Furthermore, the eigenfunctions $\phi_{n}^{N}\in H^{s}(D)\hookrightarrow C(D)$. Let $N\in\mathbb{N}_{+}$ be sufficiently large and $M:=\lfloor N^{\frac{1}{2s/d+1}}\rfloor$. Assume the spectral gap condition (3.7) to be valid. When $0<h\leq h_{0}\ll N^{-\frac{1}{s}}$ for some sufficiently small $h_{0}$, there holds

[TABLE]

Proof.

The proof of (5.3) and (5.4) follows directly from basic operator theory. The bound (5.5) is a result of (2.13) and (5.4). By the definitions of $\phi_{n}^{N}$ and $\phi_{n}^{N,h}$, we obtain

[TABLE]

Together with Proposition 4.1 and (5.5), this yields

[TABLE]

Since $h\ll N^{-\frac{1}{s}}$, the second term can be bounded from above by the third term, and this completes the proof. ∎

Remark 5.4**.**

In a similar manner as in the proof to (5.6), if $N\in\mathbb{N}_{+}$ is sufficiently large and $M:=\lfloor N^{\frac{1}{2s/d+1}}\rfloor$, one can show that

[TABLE]

Combining (5.7) and Theorem 5.3, analogously to Proposition 3.1, we can derive the following estimate.

Proposition 5.4**.**

Let $N$ be sufficiently large and $M\leq\lfloor N^{\frac{1}{2s/d+1}}\rfloor$. Furthermore, let (3.7) be satisfied. Then there holds

[TABLE]

Proof.

This result follows from an application of the triangle inequality and (5.7). ∎

Proposition 5.5**.**

Let $N$ be sufficiently large and $M\leq\lfloor N^{\frac{1}{2s/d+1}}\rfloor$. There holds

[TABLE]

Proof.

An application of the triangle inequality leads to

[TABLE]

Then, an application of the inequalities (5.2) and (5.7), Propositions 4.1 and 5.3 and Theorem 2.1 reveals the desired result. ∎

Finally, the uniform estimate between $\log\kappa$ and $\log\kappa_{M}^{N,h}$ can be derived from Propositions 5.4, 5.5 and 5.2. We then obtain the following result.

Theorem 5.4** (Uniform estimate on $M$-term KL truncation of $\log\kappa$).**

Let Assumption 5.1 hold and let $N\in\mathbb{N}_{+}$ be large and $M\leq\lfloor N^{\frac{1}{2s/d+1}}\rfloor$. Assume the spectral gap condition (3.7) to be valid. Then there exists $h_{0}\ll N^{-\frac{1}{s}}$ sufficiently small, such that

[TABLE]

for all $\quad 0<h\leq h_{0}$.

5.3 Numerical estimate for the error between $\kappa$ and $\kappa_{M}^{N,h}$

In this section, by utilizing the preceding results on $|\log\kappa-\log\kappa_{M}^{N,h}|$ together with the mean value theorem, we will derive an error estimate between $\kappa$ and $\kappa_{M}^{N,h}$. Note at this point that the results in this part can be only applied to the case when $\log\kappa$ is a normal random field. One crucial tool which we will employ repeatedly below is Fernique’s theorem. For convenience, we recall it in the following.

Theorem 5.5** (Fernique’s theorem).**

Let $E$ be a real, separable Banach space and suppose that $X$ is an $E$-valued random variable which is a centered and Gaussian in the sense that, for each $x^{*}\in E^{*}$, $\langle X,x^{*}\rangle$ is a centered, $\mathbb{R}$-valued Gaussian random variable. If R=\inf\Big{\{}r\in[0,\infty):\mathbb{P}(\|X\|_{E}\leq r)\geq\frac{3}{4}\Big{\}}, then

[TABLE]

First, we give a priori bounds on $\kappa$ and $\kappa_{M}^{N,h}$.

Proposition 5.6**.**

Let Assumption 5.1 hold and let $N\in\mathbb{N}_{+}$ be large and $M\leq\lfloor N^{\frac{1}{2s/d+1}}\rfloor$. Assume the spectral gap condition (3.7) to be valid. Then there exists $h_{0}\ll N^{-\frac{1}{s}}$ sufficiently small such that, for all $\quad 0<h\leq h_{0}$, there holds

[TABLE]

Proof.

Note that $\log\kappa$ is a symmetric Gaussian random variable defined on $\Omega$ and valued in $C(D)$. By Fernique’s theorem, there exists $\alpha>0$ such that

[TABLE]

Hence, by Young’s inequality, we obtain

[TABLE]

and (5.8) leads to

[TABLE]

This shows the first assertion. The second one can be obtained in a similar manner. ∎

Now we can state the main result of this section.

Theorem 5.6**.**

Let Assumption 5.1 hold. Let $N\in\mathbb{N}_{+}$ be sufficiently large, let $M\leq\lfloor N^{\frac{1}{2s/d+1}}\rfloor$ and $h\ll N^{-\frac{1}{s}}$. Assume the spectral gap condition (3.7) to be valid. Then, for all $p<2$, there holds

[TABLE]

Proof.

The mean value theorem indicates

[TABLE]

This, combined with Hölder’s inequality, leads to

[TABLE]

where $1/p=1/2+1/q$. In view of Theorem 5.4 and Proposition 5.6, this proves (5.9). The second assertion (5.10) can be shown similarly using Theorem 5.1. This completes the proof. ∎

6 Application to elliptic PDEs with random diffusion coefficient

In this section, we use the results of Theorem 5.6 to analyze a model order reduction algorithm for a class of elliptic PDEs with lognormal random coefficient in the multi-query context. In the algorithm, we apply the Karhunen-Loève approximation to the stochastic diffusion coefficient $\kappa(y,x)$ to arrive at a truncated model with finite-dimensional noise. We shall provide an error analysis below. Throughout this section, we assume that the conditions of Theorem 5.6 are satisfied.

Let $D$ be an open bounded domain in $\mathbb{R}^{d}$ with a strong local Lipchitz boundary and let $(\Omega,\Sigma,\mathcal{P})$ be a given probability space. Consider the elliptic PDE with random coefficient

[TABLE]

for a.e. $y\in\Omega$, where the elliptic operator $\mathcal{L}$ is defined by

[TABLE]

and $\nabla$ denotes the derivative with respect to the spatial variable $x$. We assume the force term $f$ to be in $H^{-1}(D)$. In the model problem (6.1), the dependence of the diffusion coefficient $\kappa(y,x)$ on a stochastic variable $y\in\Omega$ reflects imprecise knowledge or lack of information.

The extra-coordinate $y$ poses significant computational challenges. One popular approach is the stochastic Galerkin method [3]. There, one often approximates the stochastic diffusion coefficient $\kappa(y,x)$ by a finite sum of products of deterministic and stochastic orthogonal bases (with respect to a certain probability measure). This gives a computationally more tractable finite-dimensional noise model. There, the choice of the employed orthogonal basis is crucial for the accurate and efficient approximation to $\kappa(y,x)$. In this work, we consider the KL approximation $\kappa_{M}^{N,h}(y,x)$ of the random field $\kappa(y,x)$ in (2.14).

First, we specify the functional analytic setting. Let $V=H^{1}_{0}(D)$ and let $H^{-1}(D)$ be its dual space. Then, for any given $y\in\Omega$, the weak formulation of problem (6.1) is to find $u(y,\cdot)\in V$ such that

[TABLE]

We first discuss the well-posedness of problem (6.2) for each $y\in\Omega$, which was proven in [11, Theorem 2.2]. By Assumption 5.1, $\kappa(y,x)\in C(D)$ a.e.. Let $\kappa_{\text{min}}(y):=\min\limits_{x\in\bar{D}}\kappa(y,x)$ and $\kappa_{\text{max}}(y):=\max\limits_{x\in\bar{D}}\kappa(y,x)$ for all $y\in\Omega$.

Proposition 6.1** (Ellipticity and boundedness of $\kappa$).**

The following bounds hold:

[TABLE]

Proposition 6.1, together with the Lax-Milgram theorem, guarantees that the weak formulation (6.2) is well-posed. Furthermore,

[TABLE]

After substituting the numerical KL approximation $\kappa_{M}^{N,h}(y,x)$ of the diffusion coefficient $\kappa(y,x)$ into problem (6.1), we arrive at a truncated problem with finite-dimensional noise: For a.e. $y\in\Omega$

[TABLE]

where $\mathcal{L}_{M}^{N,h}$ is the elliptic differential operator with the diffusion coefficient $\kappa_{M}^{N,h}$. The corresponding weak formulation is then to find $u_{M}^{N,h}(y,\cdot)\in V$ such that

[TABLE]

for any given $y\in\Omega$. Analogous to the continuous case, let $\kappa_{\text{min}}^{M,N,h}(y):=\min\limits_{x\in\bar{D}}\kappa_{M}^{N,h}(y,x)$ and $\kappa_{\text{max}}^{M,N,h}(y):=\max\limits_{x\in\bar{D}}\kappa_{M}^{N,h}(y,x)$ for all $y\in\Omega$. We can then state the well-posedness of problem (6.5).

Proposition 6.2** (Ellipticity and boundedness of $\kappa_{M}^{N,h}$).**

The following bounds hold:

[TABLE]

Proof.

This follows from the proof of [11, Theorem 2.2]. ∎

Due to Proposition 6.2 and the Lax-Milgram theorem, we obtain the well-posedness of problem (6.5). Furthermore, (6.5) and Proposition 6.2, together with Poincarè’s inequality, give the following a priori estimate

[TABLE]

The next result quantifies the effect of the perturbation of the coefficient $\kappa(y,x)$ on the solution $u(y,x)$.

Theorem 6.1**.**

Let Assumption 5.1 hold. Let $N\in\mathbb{N}_{+}$ be large, let $M\leq\lfloor N^{\frac{1}{2s/d+1}}\rfloor$ and $h\ll N^{-\frac{1}{s}}$. Assume the spectral gap condition (3.7) to be valid. Let $u$ and $u_{M}^{N,h}$ be solutions to (6.1) and (6.4), respectively. Then for all $p<2$, there holds

[TABLE]

Proof.

From the weak formulations for $u(y,x)$ and $u_{M}^{N,h}(y,x)$, cf. (6.2) and (6.5), we obtain for any $y\in\Omega$

[TABLE]

By setting $v=u-u_{M}^{N,h}\in V$ in the weak formulation (6.7) and using Proposition 6.1 and the generalized Hölder inequality, we have

[TABLE]

Consequently, we arrive at

[TABLE]

Finally, taking the $L^{p}(\Omega)$-norm on both sides and employing the generalized Hölder’s inequality, combined with Theorem 5.6, Proposition 6.1 and the *a priori *estimate (6.6), shows the desired result. ∎

Remark 6.1**.**

Note that this work is mainly concerned with the numerical approximation of the lognormal random coefficient. Therefore, we refrain from discussing the important issue of the numerical approximation of the associated elliptic problems, i.e., Problems (6.1) and (6.4). We refer to [19] for related results in this direction.

7 Numerical simulation

In this section, we provide numerical tests to verify the theoretical results presented in Section 5. Recall that the three parameters $M$, $N$ and $h$ denote the number of terms in the KL approximation, the number of sampling points and the mesh size. These parameters determine directly the computational cost involved.

We take $\Omega:=\mathbb{R}^{d^{\prime}}$ in the following simulation. In order to obtain the $M$-term KL truncation estimate to $\log\kappa(y,x)$ in the form of (2.14), we employ the fast CBC construction of randomly shifted lattice rules in the unanchored space [20] to estimate the kernel $\mathcal{R}_{N}(x,x^{\prime})\in L^{2}(D\times D)$ defined in (2.9). To this end, we employ the unanchored space $\mathcal{F}(\mathbb{R}^{d^{\prime}})$ as defined by (3), where the weight function and the weight parameters are

[TABLE]

We apply the CBC method [20, Algorithm 6] to derive the generating vector $\boldsymbol{z}\in[0,1)^{d^{\prime}}$ with the number of sampling points being $N=1009$. To this end, let the shift $\boldsymbol{\Delta}\in[0,1]^{d^{\prime}}$ be an i.i.d uniformly distributed vector. Then we obtain the randomly shifted (rank-1) lattice rule by formula (3.3). Now, $R_{N}(x,x^{\prime})$ in (2.9) can be approximated by taking $y_{i}:=\boldsymbol{\phi}^{-1}(\xi_{i})$ for $i=1,2,\cdots,N$.

The bivariate functions $\log\kappa(\cdot,x)$ employed in the following examples belong to $\mathcal{F}(\mathbb{R}^{d^{\prime}})$ for all $x\in D$. Thus, using the CBC Algorithm to calculate $\mathcal{R}_{N}(x,x^{\prime})\in L^{2}(D\times D)$ as defined in (2.9) yields a shifted-averaged worse-case error of $\mathcal{O}(N^{-1+\delta})$ for any $\delta>0$ with the construction cost of $\mathcal{O}(d^{\prime}N\log(N))$.

7.1 Example 1: $d^{\prime}=1$ and $N=1009$

Let

[TABLE]

see Fig. 1 for an illustration. One can verify that $\log\kappa(y,x)\in L^{2}(\Omega,H^{3/2-\delta}(D))$ for any $\delta>0$ with the physical domain $D:=[0,1]$ and the stochastic domain $\Omega:=\mathbb{R}$. Thus, we have $s:=3/2-\delta$ in this case. According to the definition of the finite element space $V_{h}$, we will use conforming quadratic finite element. Now as in Remark 5.1, let the tolerance be chosen as $\epsilon:=0.1$. Note that we always fix the sampling points $N=1009$. Then we can take the number of truncation terms $M\approx 5$ and the mesh size $h\approx 1/101$.

Since the dimension of the stochastic domain equals to one, we can choose the generating vector $\boldsymbol{z}:=1$. We then can derive the sampling points $\{y_{i}\}_{i=1}^{N}:=\{\boldsymbol{\phi}^{-1}(\xi_{i})\}_{i=1}^{N}$ from formula (3.3). The shifted-averaged worse-case error is 8.0171e-6. We present in Table 1 the root mean square error between $\log\kappa$ and $\log\kappa_{M}^{N,h}$ for different numbers of truncation terms $M=2,4,8$ and different mesh sizes $h:=1/16,1/64,1/128,1/256$.

Now, let us compare these computed results with the values that were predicted from our theory. To this end, we plug the fixed number of sampling points $N=1009$ into Remark 5.1 and derive that we can take the accuracy $\epsilon:=0.1$, the number of truncation terms $M:=5$ and the mesh size $h:=1/101$. Indeed, for $(M,h):=(5,1/101)$, we also obtain the optimal error in Table 1. This shows that our estimates are quite sharp and involve just small constants.

7.2 Example 2: $d^{\prime}=10,100$ and $N=1009$

Let the bivariate function $\kappa(y,x)$ be given by

[TABLE]

Then one can verify that $\log\kappa(y,x)\in L^{2}(D,H^{t}(\Omega))$ for any $t>0$, with the physical domain $D:=[0,1]$ and the stochastic domain $\Omega:=\mathbb{R}^{d^{\prime}}$. According to the definition of the finite element space $V_{h}$, we will use spectral element up to degree 10. The basis functions are Lagrange interpolation polynomials through the local Gauss-Lobatto integration points defined per element. Due to Theorem 2.1, the eigenvalues decay very fast since $s=\infty$ in this case.

We present the root mean square errors between $\log\kappa$ and $\log\kappa_{M}^{N,h}$ in Tables 2 and 3 for different mesh sizes $h:=1/16,1/64,1/128\text{ and }1/256$ and different numbers of truncation terms $M:=2,4\text{ and }8$ with $d^{\prime}=10$ and $d^{\prime}=100$, respectively.

Furthermore, for our fixed number of sampling points $N:=1009$ and for the accuracy $\epsilon:=0.1$, we can compare our computed results with the predicted ones due to Remark 5.1. We see that our estimates are again qualitatively quite sharp and involve just small constants.

8 Concluding remarks

In this work, we have analyzed the numerical approximation error in the Karhunen-Loève expansion to log normal random coefficients. We derived the numerical error in terms of the number $M$ of terms in the Karhunen-Loève expansion, the number $N$ of QMC sampling points to estimate the covariance function and the mesh size $h$ for the conforming Galerkin approximation to the eigenvalue problem. Our results show the basic relation

[TABLE]

among those three parameters, where $d$ is the dimension of the physical domain and $s$ denotes the regularity of the bivariate function in the physical domain. These results are also useful for the study of stochastic elliptic problems. We presented numerical results for one and multiple stochastic dimensions to support our theory.

The QMC method can be replaced by some properly adapted sparse grid method, if there is higher mixed regularity in $\log\kappa$ present with respect to the stochastic variables. Analogously, if the physical problem possesses higher regularity, then a more suitable FEM of higher order can be utilized. Then, of course, the sampling estimate, the Galerkin estimate and the resulting error estimates have to be modified accordingly. This would lead to a different balancing of the terms that in Remark 5.1.

Acknowledgements

The authors were supported by the Hausdorff Center for Mathematics in Bonn and the Sonderforschungsbereich 1060 The Mathematics of Emergent Effects funded by the Deutsche Forschungsgemeinschaft. Guanglian Li acknowledges the support from the Royal Society via the Newton International fellowship. Part of this work was done during her visit to IPAM in the Long program: Computational Issues in Oil Field Applications. We thank Prof. Dr. Markus Bachmayr for fruitful discussions.