# On the Numerical Approximation of the Karhunen-Loève Expansion for Random
Fields with Random Discrete Data

Michael Griebel and Guanglian Li and Christian Rieger Institut für Numerische
Simulation, Universität Bonn, Friedrich-Hirzebruch-Allee 7, 53115 Bonn,
Germany and Fraunhofer SCAI, Schloss Birlinghoven, 53754 Sankt Augustin,
Germany, E-mail<EMAIL_ADDRESS>Department of Mathematics, The
University of Hong Kong, Pokfulam Road, Hong Kong. E-mail:
<EMAIL_ADDRESS>Marburg, Fachbereich Mathematik und
Informatik, AG Numerik, Hans-Meerwein-Straße 6, 35032 Marburg.
<EMAIL_ADDRESS>
###### Abstract

Many physical and mathematical models involve random fields in their input
data. Examples are ordinary differential equations, partial differential
equations and integro–differential equations with uncertainties in the
coefficient functions described by random fields. They also play a dominant
role in problems in machine learning. In this article, we do not assume to
have knowledge of the moments or expansion terms of the random fields but we
instead have only given discretized samples for them. We thus model some
measurement process for this discrete information and then approximate the
covariance operator of the original random field. Of course, the true
covariance operator is of infinite rank and hence we can not assume to get an
accurate approximation from a finite number of spatially discretized
observations. On the other hand, smoothness of the true (unknown) covariance
function results in effective low rank approximations to the true covariance
operator. We derive explicit error estimates that involve the finite rank
approximation error of the covariance operator, the Monte-Carlo-type errors
for sampling in the stochastic domain and the numerical discretization error
in the physical domain. This permits to give sufficient conditions on the
three discretization parameters to guarantee that an error below a prescribed
accuracy $\varepsilon$ is achieved.

Keywords: covariance operators, eigenvalue decay, approximation of Gaussian-
type random fields, error estimates, Galerkin methods for eigenvalues, finite
elements, tapering estimators for sample covariance

Subject classification: 41A25, 41A35, 60F10, 65D40

##  1 Introduction

Mathematical models with random coefficients or random input data have been
widely employed to describe applications that are affected by a certain amount
of uncertainty arising from imperfect or insufficient information about the
problem. The range of applications is broad and diverse and includes
uncertainty quantification with ordinary differential equations, partial
differential equations and integro-differential equations or problems in
machine learning and data analysis in, e.g. oil field modeling, quantum
mechanics or finance.

For instance in oil field modeling, a common source of uncertainty stems from
the unknown soil parameters. Often, one assumes a statistical model for the
soil, for instance a Gaussian process with specified mean and covariance, and
fits the possibly remaining hyper-parameters of the process to the measured
data. Moreover in machine learning, a common source of randomness stems from
noisy observed data. This renders all subsequent quantities to be random
fields. Again one often assumes a statistical model for the noise, for
instance a Gaussian model, which perfectly fits to the linear Bayesian
framework as then the posterior is also Gaussian. Anyway, let now a
probability space $(\Omega,\mathcal{F},\mathbb{P})$ be given, where
$\mathcal{F}$ denotes a $\sigma$-algebra of measurable sets and $\mathbb{P}$
is a probability measure on this $\sigma$-algebra. On this probability space,
we consider a Gaussian random field

$\displaystyle\kappa:\Omega\times\mathcal{D}\to\mathbb{R}.$

Having such a Gaussian random field at hand, one usually parametrizes it by
means of the Karhunen-Loève (KL) expansion or a polynomial chaos (PC)
expansion [9], which fits to the $L^{2}$ setting, i.e.,

$\displaystyle\kappa(\boldsymbol{\omega},\boldsymbol{x})=\sum_{\ell=1}^{\infty}\psi_{\ell}(\boldsymbol{\omega})\phi_{\ell}(\boldsymbol{x}).$

This greatly facilitates the subsequent treatment of problems with $\kappa$
modeling their random input or their random coefficient functions, e.g., by
the stochastic Galerkin method or the stochastic collocation method.
Alternatively, in the Banach space setting, one may expand the random field
with respect to the hierarchical Faber basis or some wavelet type basis; see
[3] for details.

To further fix notation, a realization of the associated stochastic process is
the function

$\displaystyle\kappa(\boldsymbol{\omega},\cdot):\mathcal{D}\to\mathbb{R}\quad\boldsymbol{\omega}\in\Omega\text{
fixed}.$

We now can interpret the stochastic process $\kappa$ as a
$\mathbb{R}^{\mathcal{D}}$-valued random variable, i.e.,

$\displaystyle
K:\Omega\to\mathbb{R}^{\mathcal{D}},\quad\omega\mapsto\kappa(\boldsymbol{\omega},\cdot):\mathcal{D}\to\mathbb{R},$

where $\mathbb{R}^{\mathcal{D}}$ denotes the set of all maps
$\mathcal{D}\to\mathbb{R}$. Hence, using the notion of a push forward measure,
we define a distribution of this random variable $K$, and thus the random
process $\kappa$, as

$\displaystyle\mu=\mathbb{P}\circ K^{-1}.$

Any Gaussian random field $\kappa$ is determined completely by its first two
moments, i.e.,

$\displaystyle\mathbb{E}\left[\kappa\right]$
$\displaystyle:\mathcal{D}\to\mathbb{R},\quad\boldsymbol{x}\mapsto\mathbb{E}[\kappa]=\mathbb{E}\left[\kappa(\cdot,\boldsymbol{x})\right]\quad\text{and}$
$\displaystyle\text{Cov}_{\kappa}$
$\displaystyle:\mathcal{D}\times\mathcal{D}\to\mathbb{R},\quad(\boldsymbol{x},\boldsymbol{x}^{\prime})\mapsto\mathbb{E}\left[\left(\kappa-\mathbb{E}\left[\kappa(\cdot,\boldsymbol{x})\right]\right)\left(\kappa(\cdot,\boldsymbol{x}^{\prime})-\mathbb{E}\left[\kappa(\cdot,\boldsymbol{x}^{\prime})\right]\right)\right].$

For the sake of simplicity, we will assume that $\kappa$ is a centered
Gaussian field, i.e., $\mathbb{E}\left[K\right]\equiv 0$. This yields

$\displaystyle\text{Cov}_{\kappa}(\boldsymbol{x},\boldsymbol{x}^{\prime})=:R(\boldsymbol{x},\boldsymbol{x}^{\prime})=\int_{\Omega}\kappa(\boldsymbol{\omega},\boldsymbol{x})\kappa(\boldsymbol{\omega},\boldsymbol{x}^{\prime})\,\mathrm{d}\mathbb{P}(\boldsymbol{\omega})=\sum_{\ell=1}^{\infty}\lambda_{\ell}\phi_{\ell}(\boldsymbol{x})\phi_{\ell}(\boldsymbol{x}^{\prime}),$
(1.1)

where the series is a Mercer-type expansion. We will assume that the
eigenvalues are sorted in strictly decreasing order, i.e.,

$\displaystyle\lambda_{1}>\lambda_{2}>\dots\geq 0.$ (1.2)

In the case of multiple eigenvalues the projection onto the eigenspaces, i.e.,

$\displaystyle
R(\boldsymbol{x},\boldsymbol{x}^{\prime})=\sum_{\ell\in\mathbb{N}}\lambda_{\ell}\Pi_{\ell}(\boldsymbol{x},\boldsymbol{x}^{\prime})=\sum_{\ell\in\mathbb{N}}\lambda_{\ell}\sum_{k=1}^{N(\ell,k)}\phi_{\ell,k}(\boldsymbol{x})\phi_{\ell,k}(\boldsymbol{x^{\prime}}),$
(1.3)

is to be employed of course, but for reasons of simplicity we will stick here
to the notationally simpler Mercer series (1.1) instead. By repeating multiple
eigenvalues and modifying our approach accordingly, the more general case of
multiple eigenvalues can be covered in an analogous way.

Usually, such an expansion will be truncated yielding a so-called finite noise
approximation. To be precise, the approximation

$\displaystyle R(\boldsymbol{x},\boldsymbol{x}^{\prime})\approx
R^{L}(\boldsymbol{x},\boldsymbol{x}^{\prime}):=\sum_{\ell=1}^{L}\lambda_{\ell}\phi_{\ell}(\boldsymbol{x})\phi_{\ell}(\boldsymbol{x}^{\prime})$
(1.4)

is considered, where the quality of the approximation relies on sharp
eigenvalue decay estimates as given in Theorem 3.1 later on. This approach
however needs the set $\\{\phi_{\ell}\\}$ of continuous eigenfunctions of the
covariance operator $R$ in the first place, which is in general unknown. While
its knowledge is often a priorily assumed in the literature for further
numerical analysis, such complete information is mostly not available in
practice and such an approach is thus often not viable.

Therefore, in the following, we consider a different situation: We do not
assume to have access to the moments of the random field $\kappa$ in our
analysis. Thus we can not compute its Karhunen Loève expansion directly as an
eigenvalue problem for the integral operator with the covariance functions as
kernel. Instead, we will only assume

$\displaystyle
C_{\infty}:=\left(\int_{\Omega}\left\|\kappa\left(\boldsymbol{\omega},\cdot\right)\right\|_{C(\mathcal{D})}^{2}\,\mathrm{d}\mathbb{P}(\boldsymbol{\omega})\right)^{\frac{1}{2}}<\infty.$
(1.5)

Note that (1.5) implies uniform bounds on mean and covariance functions of the
random field $\kappa$. Indeed, the definition of the mean function and an
application of the Cauchy-Schwarz inequality leads to

$\displaystyle\left|\mathbb{E}\left[\kappa(\cdot,\boldsymbol{x})\right]\right|=\left|\int_{\Omega}\kappa(\cdot,\boldsymbol{x})\,\mathrm{d}\mathbb{P}\right|\leq\int_{\Omega}\left|\kappa(\cdot,\boldsymbol{x})\right|\,\mathrm{d}\mathbb{P}\leq\left(\int_{\Omega}\left|\kappa(\cdot,\boldsymbol{x})\right|^{2}\,\mathrm{d}\mathbb{P}\right)^{\frac{1}{2}}\left(\int_{\Omega}1\,\mathrm{d}\mathbb{P}\right)^{\frac{1}{2}}\leq
C_{\infty}.$ (1.6)

Moreover, we can obtain by the definition of the covariance function and an
application of the Cauchy-Schwarz inequality the bound

$\displaystyle\text{Cov}(\boldsymbol{x},\boldsymbol{x}^{\prime})$
$\displaystyle=\int_{\Omega}\left(\kappa(\cdot,\boldsymbol{x})-\mathbb{E}\left[\kappa(\cdot,\boldsymbol{x})\right]\right)\left(\kappa(\cdot,\boldsymbol{x}^{\prime})-\mathbb{E}\left[\kappa(\cdot,\boldsymbol{x}^{\prime})\right]\right)\,\mathrm{d}\mathbb{P}$
$\displaystyle\leq\left(\int_{\Omega}\left(\kappa(\cdot,\boldsymbol{x})-\mathbb{E}\left[\kappa(\cdot,\boldsymbol{x})\right]\right)^{2}\,\mathrm{d}\mathbb{P}\right)^{\frac{1}{2}}\left(\int_{\Omega}\left(\kappa(\cdot,\boldsymbol{x}^{\prime})-\mathbb{E}\left[\kappa(\cdot,\boldsymbol{x}^{\prime})\right]\right)^{2}\,\mathrm{d}\mathbb{P}\right)^{\frac{1}{2}}$
$\displaystyle\leq
2\left(\int_{\Omega}\kappa(\cdot,\boldsymbol{x})^{2}+\mathbb{E}\left[\kappa(\cdot,\boldsymbol{x})\right]^{2}\,\mathrm{d}\mathbb{P}\right)^{\frac{1}{2}}\left(\int_{\Omega}\kappa(\cdot,\boldsymbol{x}^{\prime})^{2}+\mathbb{E}\left[\kappa(\cdot,\boldsymbol{x}^{\prime})\right]^{2}\,\mathrm{d}\mathbb{P}\right)^{\frac{1}{2}}\leq
4C_{\infty}^{2}.$ (1.7)

Furthermore, we define
$C_{H^{s}}:=\left(\int_{\Omega}\left\|\kappa\left(\boldsymbol{\omega},\cdot\right)\right\|_{H^{s}\left(\mathcal{D}\right)}^{2}\,\mathrm{d}\mathbb{P}(\boldsymbol{\omega})\right)^{\frac{1}{2}}<\infty$
for later use as well. We refer to [16, Theorem 5.2] for conditions on
$\kappa$ which allow for such types of bounds.

Next we assume some spatial regularity of the random field $\kappa$.

######  Assumption 1.1.

Let $\kappa:\Omega\times\mathcal{D}\to\mathbb{R}$ and thus
$K:\Omega\to\mathbb{R}^{\mathcal{D}}$ be contained in the Bochner space

$\displaystyle L^{\infty}(\Omega,H^{s}(\mathcal{D}))\quad\text{for some }s\geq
0,$ (1.8)

where, for $s=0$, we have $H^{0}(\mathcal{D})=L^{2}(\mathcal{D})$.
Alternatively, let $\kappa:\Omega\times\mathcal{D}\to\mathbb{R}$ and thus
$K:\Omega\to\mathbb{R}^{\mathcal{D}}$ be contained in the Bochner space

$\displaystyle L^{\infty}(\Omega,H^{s}(\mathcal{D}))\quad\text{for some
}s>\frac{d}{2}.$ (1.9)

Now we introduce a discretization of the function space
$\mathbb{R}^{\mathcal{D}}:=\\{f:\mathcal{D}\to\mathbb{R}\\}$. To be precise,
we consider an underlying finite element space of the form

$\displaystyle\mathcal{V}_{h}:=\\{v\in H^{1}(\mathcal{D}):v|_{K}\in P^{\lceil
s\rceil}(K)\text{ for all
}K\in\mathcal{T}_{h}\\}\subset\mathbb{R}^{\mathcal{D}},$ (1.10)

where $\mathcal{T}_{h}$ is a regular quasi-uniform triangulation over the
physical domain $\mathcal{D}$ with a maximal mesh size $h$ and $\lceil
s\rceil\in\mathbb{N}$ will depend on the spatial regularity of $\kappa$. We
denote the dimension of $\mathcal{V}_{h}$ by $Q_{h}<\infty$. Moreover, we
expect $Q_{h}=\mathcal{O}(h^{-d}s^{d})$. For the rest of the paper, we will
fix some arbitrary basis of $\mathcal{V}_{h}$

$\displaystyle\mathcal{V}_{h}=\text{span}\left\\{\theta_{k}^{(h)}\ :\ 1\leq
k\leq Q_{h}\right\\}.$

Here we may employ a nodal basis as well as a discrete orthonormal basis later
on. Let $\boldsymbol{\theta}^{(h)}$ be

$\displaystyle\boldsymbol{\theta}^{(h)}:\mathcal{D}\to\mathbb{R}^{Q_{h}},\quad\boldsymbol{x}\mapsto\left(\theta^{(h)}_{1},\dots,\theta^{(h)}_{Q_{h}}\right)^{T}.$

Any function $v^{(h)}\in\mathcal{V}_{h}$ then admits the expression

$\displaystyle
v^{(h)}=\boldsymbol{V}^{(h)}\cdot\boldsymbol{\theta}^{(h)}\quad\text{with}\quad\boldsymbol{V}^{(h)}=\left(V^{(h)}_{1},\dots,V^{(h)}_{Q_{h}}\right)^{T}\in\mathbb{R}^{Q_{h}},$

and each function $v^{(h)}\in\mathcal{V}_{h}$ can be identified with its
coefficient vector $\boldsymbol{V}^{(h)}\in\mathbb{R}^{Q_{h}}$.

Next we consider the projection map
$\Pi_{\mathcal{V}_{h}}:L^{2}(\mathcal{D})\to\mathcal{V}_{h}$. Thus $s=0$ in
(1.8) from Assumption 1.1, which is here the weakest assumption possible.
Given such regularity, we recall the approximation power of the
$L^{2}$-projection
$\Pi_{\mathcal{V}_{h};L^{2}}:L^{2}(\mathcal{D})\to\mathcal{V}_{h}$, which is
given in [5, Theorem 4.4.20] as

$\displaystyle\left\|v-\Pi_{\mathcal{V}_{h}}v\right\|_{L^{2}\left(\mathcal{D}\right)}$
$\displaystyle\leq
C_{\Pi_{\mathcal{V}_{h};L^{2}}}h^{s}\left\|v\right\|_{H^{s}\left(\mathcal{D}\right)}\text{
for all }v\in H^{s}(\mathcal{D}),$ (1.11)
$\displaystyle\left\|v-\Pi_{\mathcal{V}_{h}}v\right\|_{C\left(\mathcal{D}\right)}$
$\displaystyle\leq
C_{\Pi_{\mathcal{V}_{h};C}}h^{s-\frac{d}{2}}\left\|v\right\|_{H^{s}\left(\mathcal{D}\right)}\text{
for all }v\in H^{s}(\mathcal{D})\text{ for }s>d/2.$ (1.12)

Here, the positive constants $C_{\Pi_{\mathcal{V}_{h};L^{2}}}$ and
$C_{\Pi_{\mathcal{V}_{h};C}}$ depend only on the shape regularity parameter of
$\mathcal{T}_{h}$ and are independent of the mesh size $h$.

In contrast to other approaches where the continuous random field
$\kappa:\Omega\times\mathcal{D}\to\mathbb{R}$ and thus
$K:\Omega\to\mathbb{R}^{\mathcal{D}}$ are assumed to be completely known, we
in the following only assume to have access to

$\displaystyle\kappa^{(h)}(\boldsymbol{\omega},\boldsymbol{x})=\boldsymbol{K}^{(h)}(\boldsymbol{\omega})\cdot\boldsymbol{\theta}^{(h)}(\boldsymbol{x}),$
(1.13)

where $\boldsymbol{K}^{(h)}(\boldsymbol{\omega})$ is a multivariate random
variable defined by

$\displaystyle\boldsymbol{K}^{(h)}:\Omega\to\mathbb{R}^{Q_{h}},\quad\boldsymbol{\omega}\mapsto\boldsymbol{K}^{(h)}(\boldsymbol{\omega})=\left(K_{1}^{(h)}(\boldsymbol{\omega}),\dots,K_{Q_{h}}^{(h)}(\boldsymbol{\omega})\right)^{T}.$
(1.14)

Furthermore we assume that the coefficients, i.e. the random variables
$\\{K^{(h)}_{j}(\boldsymbol{\omega})\\}_{j=1,\cdots,Q_{h}}$, are independent
identically distributed with a common distribution $\mu^{(h)}$. As we work
with centered Gaussian random fields, we have
$\mathbb{E}\left[\boldsymbol{K}^{(h)}\right]=\boldsymbol{0}\in\mathbb{R}^{Q_{h}}$.

This resembles more closely the practical situation in real applications where
in general the full random field $K$ is never completely accessible but only
given at a finite number of sample points which stem from independent
measurements in the first place. The random vector $\boldsymbol{K}^{(h)}$ is
indeed the information we practically have on the random input.

Note that an alternative approach would be to consider directly the _Carleman
operator_ $R$ (see (1.1)) for centered fields and to approximate it by
$\bar{R}^{(h)}$, i.e.,

$\displaystyle R(\boldsymbol{x},\boldsymbol{x}^{\prime})$
$\displaystyle\approx
R^{(h)}(\boldsymbol{x},\boldsymbol{x}^{\prime}):=\int_{\Omega}\kappa^{(h)}(\boldsymbol{\omega},\boldsymbol{x})\kappa^{(h)}(\boldsymbol{\omega},\boldsymbol{x}^{\prime})\,\mathrm{d}\mathbb{P}(\boldsymbol{\omega})$
$\displaystyle\approx
M^{-1}\sum_{m=1}^{M}\kappa^{(h)}(\boldsymbol{\omega}_{m},\boldsymbol{x})\kappa^{(h)}(\boldsymbol{\omega}_{m},\boldsymbol{x}^{\prime})=:\bar{R}^{(h)}(\boldsymbol{x},\boldsymbol{x}^{\prime}).$

This kernel $\bar{R}^{(h)}$ gives rise to an integral operator
$f\mapsto\int_{\mathcal{D}}f(\boldsymbol{x})\bar{R}^{(h)}(\boldsymbol{x},\cdot)d\boldsymbol{x}$
for which the stochastics is integrated out. Our aim is however to derive a
kernel via sampling estimation instead. But this sampling is not involved in
the above reference kernels ${R}$ and $\bar{R}^{(h)}$ at all. Therefore, we
consider the kernel $R^{(h)}$ as integral kernel on the discretized space
$\mathcal{V}_{h}$. In particular, we will have a Mercer-type expansion of the
form

$\displaystyle
R^{(h)}(\boldsymbol{x},\boldsymbol{x}^{\prime})=\sum_{\ell=1}^{Q_{h}}\lambda^{(h)}_{\ell}\phi_{\ell}^{(h)}(\boldsymbol{x})\phi_{\ell}^{(h)}(\boldsymbol{x}^{\prime})\quad\text{with}\quad\phi^{(h)}_{\ell}(\boldsymbol{x})=\boldsymbol{\Phi}^{(h)}_{\ell}\cdot\boldsymbol{\theta}^{(h)}(\boldsymbol{x}).$

This gives rise to a generalized matrix eigenvalue problem and we merely
estimate the associated stiffness matrix from samples.

Now, due to (1.9) from Assumption 1.1, we obtain the estimate

$\displaystyle\left\|\Pi_{\mathcal{V}_{h}}\kappa(\boldsymbol{\omega},\cdot)\right\|_{C(\mathcal{D})}$
$\displaystyle\leq\left\|\kappa\left(\boldsymbol{\omega},\cdot\right)\right\|_{C(\mathcal{D})}+\left\|\kappa\left(\boldsymbol{\omega},\cdot\right)-\Pi_{\mathcal{V}_{h}}\kappa(\boldsymbol{\omega},\cdot)\right\|_{C(\mathcal{D})}$
$\displaystyle\leq\left\|\
\kappa(\boldsymbol{\omega},\cdot)\right\|_{C(\mathcal{D})}+C_{\Pi_{\mathcal{V}_{h};C}}h^{s-\frac{d}{2}}\left\|\kappa\left(\boldsymbol{\omega},\cdot\right)\right\|_{H^{s}\left(\mathcal{D}\right)}\quad\text{for
}\mathbb{P}\text{ almost all }\boldsymbol{\omega}\in\Omega.$ (1.15)

Hence, after squaring and integrating over $\Omega$, we have

$\displaystyle\int_{\Omega}\left\|\Pi_{\mathcal{V}_{h}}\kappa(\boldsymbol{\omega},\cdot)\right\|_{C(\mathcal{D})}^{2}\,\mathrm{d}\mathbb{P}(\boldsymbol{\omega})\leq\int_{\Omega}\left(\left\|\
\kappa(\boldsymbol{\omega},\cdot)\right\|_{C(\mathcal{D})}+C_{\Pi_{\mathcal{V}_{h};C}}h^{s-\frac{d}{2}}\left\|\kappa\left(\boldsymbol{\omega},\cdot\right)\right\|_{H^{s}\left(\mathcal{D}\right)}\right)^{2}\,\mathrm{d}\mathbb{P}(\boldsymbol{\omega})$
$\displaystyle\leq 2\int_{\Omega}\left\|\
\kappa(\boldsymbol{\omega},\cdot)\right\|_{C(\mathcal{D})}^{2}+C_{\Pi_{\mathcal{V}_{h};C}}^{2}h^{2s-d}\left\|\kappa\left(\boldsymbol{\omega},\cdot\right)\right\|_{H^{s}\left(\mathcal{D}\right)}^{2}\,\mathrm{d}\mathbb{P}(\boldsymbol{\omega})\leq
2C_{\infty}^{2}+2C_{\Pi_{\mathcal{V}_{h};C}}^{2}h^{2s-d}C_{H^{s}}^{2}.$

Thus, with
$\kappa^{(h)}\left(\boldsymbol{\omega},\cdot\right):=\Pi_{\mathcal{V}_{h}}\kappa(\boldsymbol{\omega},\cdot)$,
we obtain

$\displaystyle
C_{\infty}^{(h)}:=\left(\int_{\Omega}\left\|\kappa^{(h)}\left(\boldsymbol{\omega},\cdot\right)\right\|_{C(\mathcal{D})}^{2}\,\mathrm{d}\mathbb{P}(\boldsymbol{\omega})\right)^{\frac{1}{2}}<\infty.$
(1.16)

Analogously by employing the $H^{s}$ norm instead of the $C(\mathcal{D})$
norm, we have

$\displaystyle
C_{H^{s}}^{(h)}:=\left(\int_{\Omega}\left\|\kappa^{(h)}\left(\boldsymbol{\omega},\cdot\right)\right\|_{H^{s}\left(\mathcal{D}\right)}^{2}\,\mathrm{d}\mathbb{P}(\boldsymbol{\omega})\right)^{\frac{1}{2}}\leq\left(\int_{\Omega}\left\|\kappa\left(\boldsymbol{\omega},\cdot\right)\right\|_{H^{s}\left(\mathcal{D}\right)}^{2}\,\mathrm{d}\mathbb{P}(\boldsymbol{\omega})\right)^{\frac{1}{2}}<\infty.$

As in the continuous case (1.6) and (1) , this implies

$\displaystyle\mathbb{E}\left[\kappa^{(h)}(\cdot,\boldsymbol{x})\right]\leq
C_{\infty}^{(h)}\quad\text{and}\quad\text{Cov}^{(h)}(\boldsymbol{x},\boldsymbol{x}^{\prime})\leq
4(C_{\infty}^{(h)})^{2}\quad\text{ for all
}\boldsymbol{x},\boldsymbol{x}^{\prime}\in\mathcal{D}.$ (1.17)

Furthermore, with the representation (1.13) and since we deal with (centered)
random vectors, the covariance is given as

$\displaystyle\boldsymbol{\Sigma}_{\boldsymbol{K}^{(h)}}:=\mathbb{E}\left[\boldsymbol{K}^{(h)}\otimes\boldsymbol{K}^{(h)}\right]-\mathbb{E}\left[\boldsymbol{K}^{(h)}\right]\otimes\mathbb{E}\left[\boldsymbol{K}^{(h)}\right]=\mathbb{E}\left[\boldsymbol{K}^{(h)}\otimes\boldsymbol{K}^{(h)}\right]\in\mathbb{R}^{Q_{h}\times
Q_{h}}.$ (1.18)

In contrast to the vast majority of the literature, where the covariance is
assumed to be known, we only assume to have given $M$ independent identically
distributed discrete samples
$\boldsymbol{K}^{(h)}_{1},\dots,\boldsymbol{K}^{(h)}_{M}$ of the random
variable $\kappa$. From these finitely many values we merely aim to derive an
estimate of the true covariance function. To this end, we have to deal with
the high-dimensional covariance matrix
$\boldsymbol{\Sigma}_{\boldsymbol{K}^{(h)}}\in\mathbb{R}^{Q_{h}\times Q_{h}}$.
Since $Q_{h}\sim h^{-d}$, this involves the well-known curse of dimension.
Consequently the estimation of the covariance matrix gets more challenging the
smaller the values of $h$ get. Therefore, we will assume that the covariance
matrix $\boldsymbol{\Sigma}_{\boldsymbol{K}^{(h)}}\in\mathbb{R}^{Q_{h}\times
Q_{h}}$ has a certain off-diagonal decay, that is

$\boldsymbol{\Sigma}_{\boldsymbol{K}^{(h)}}\in\mathcal{F}_{\alpha}=\mathcal{F}_{\alpha}(C_{\mathcal{F};1},C_{\mathcal{F};2}),$

where

$\displaystyle\mathcal{F}_{\alpha}:=\left\\{\boldsymbol{\Sigma}\in\mathbb{R}^{Q_{h}\times
Q_{h}}\ :\
\max_{k=1,\dots,Q_{h}}\sum_{\genfrac{}{}{0.0pt}{}{k^{\prime}=1}{\left|k^{\prime}-k\right|>c}}^{Q_{h}}\left|\boldsymbol{\Sigma}_{k,k^{\prime}}\right|\leq
C_{\mathcal{F};1}c^{-\alpha}\text{ for all }1\leq c\leq Q_{h}\text{ and
}\lambda_{\max}(\boldsymbol{\Sigma})\leq C_{\mathcal{F};2}\right\\}.$ (1.19)

Here, $\alpha$ modulates the speed of decay and
$C_{\mathcal{F};1},C_{\mathcal{F};2}$ are positive constants. This specific
class of matrices contains several relevant examples of discrete covariance
functions, see [4] and the overview [6]. Moreover, the choice of this class of
matrices will influence the approximation of continuous covariance functions
via so-called tapering weights, see (4.39). We will make use of the optimal
rates derived in [4]. We furthermore will work with generalized eigenvalue
problems. i.e., we do not compute the eigensystem of the covariance matrix
$\boldsymbol{\Sigma}_{K^{(h)}}$ but instead that of a certain transformed
matrix, see (4.14). The eigensystem of this system gives rise to a computable
Mercer-like expansion

$\displaystyle
R^{(h;M)}(\boldsymbol{x},\boldsymbol{x}^{\prime}):=\sum_{\ell=1}^{Q_{h}}\lambda^{(h;M)}_{\ell}\phi_{\ell}^{(h;M)}(\boldsymbol{x})\phi_{\ell}^{(h;M)}(\boldsymbol{x}^{\prime})\approx
R^{(L;h;M)}(\boldsymbol{x},\boldsymbol{x}^{\prime}):=\sum_{\ell=1}^{L}\lambda^{(h;M)}_{\ell}\phi_{\ell}^{(h;M)}(\boldsymbol{x})\phi_{\ell}^{(h;M)}(\boldsymbol{x}^{\prime}),$
(1.20)

where the $(h;M)$–notation indicates that we use an estimator which is based
on the $M$ finite samples
$\boldsymbol{K}^{(h)}_{1},\dots,\boldsymbol{K}^{(h)}_{M}$ and makes use of
certain appropriate tapering weights. Moreover $L$ indicates the finite-noise
approximation.

The main aim of this article to derive a bound on the approximation error
$\left\|{R}-R^{(L;h;M)}\right\|$ and of its expectation
$\mathbb{E}\left[\left\|{R}-R^{(L;h;M)}\right\|\right]$. To this end, we
couple the discretization parameters $L,h,M$ and consider their limits $h\to
0,L\to Q_{h}\to\infty,M\to\infty$ to recover the true continuous covariance
operator. For the approximation step with respect to $L$, we will make use of
the optimal decay estimates for the eigenvalues of the true covariance
operator as given [11]. For the two approximations steps involving $h$ and
$M$, we will derive corresponding estimates for the spatial discretization
error and the sampling error. Finally we properly combine the three
approximation steps. This coupling is our central result, see Theorem 5.1.
There we show that the bound

$\displaystyle\left\|{R}-R^{(L;h;M)}\right\|_{L^{2}(\mathcal{D}\times\mathcal{D})}$
$\displaystyle\lesssim
L^{-\frac{2s}{d}-\frac{1}{2}}+\left(L^{\frac{1}{2}}+G(L)\right)\left\|\widetilde{\boldsymbol{S}}^{(h)}-\widetilde{\boldsymbol{S}}^{(h;M)}\right\|_{2\to
2}\ ,$

holds with high probability, where $G(L)$ is a function of the truncation
parameter $L$ that depends on spectral properties of the true operator
associated to $R$ as given in (5.2). Moreover
$\widetilde{\boldsymbol{S}}^{(h)}$ denotes the finite element discretization
of the true operator and $\widetilde{\boldsymbol{S}}^{(h;M)}$ denotes its
sampling approximation, c.f. (4.7) and (4.13).

Since sampling involves randomness we are also interested in estimates in
expectation. Our second main result is Theorem 5.2. It states that

$\displaystyle\mathbb{E}\left[\left\|{R}-R^{(L;h;M)}\right\|_{L^{2}(\mathcal{D}\times\mathcal{D})}\right]$
$\displaystyle\lesssim
L^{-\frac{2s}{d}-\frac{1}{2}}+\left(L^{\frac{1}{2}}+G(L)\right)\tilde{\rho}^{\frac{1}{2}}_{h}(M)\lambda_{\max}\left(\boldsymbol{M}^{(h)}\right)$
$\displaystyle+L^{\frac{1}{2}}h^{-d}\exp\left(-M\rho_{1}H(L)\lambda_{\max}^{-2}\left(\boldsymbol{M}^{(h)}\right)\right),$

where $H(L)$ is again a function of the truncation parameter $L$ that depends
on spectral properties (5.29) of the true operator, $\tilde{\rho}_{h}(M)$
measures the approximation quality of our employed (tapering) estimator,
$\boldsymbol{M}^{(h)}$ is a classical finite element mass matrix and
$\rho_{1}$ denotes a constant given in (4.54) which is determined by the sub-
Gaussian property of the involved random variables. This permits to give
sufficient conditions on the three discretization parameters to guarantee that
an error below a prescribed accuracy $\varepsilon$ is achieved.

A novel contribution is the combination of our recent sharp bound (see [11])
on spectra of covariance operators with optimal statistical covariance
estimation methods, i.e., we derive an error analysis for the reconstruction
of the continuous covariance operator from finite measurements. We believe
that the presented framework is useful if only finite measured information on
uncertain problem parameters is available and no a priori modeling assumption
for the random field can be given. Such a situation is indeed often
encountered in many practical problems in uncertainty quantification and
machine learning.

The remainder of this paper is organized as follows: In Section 2 we introduce
notation and review some basic facts. In Section 3 we recall the KL-expansion
and sharp eigenvalue estimates. In Section 4, we present the reconstruction of
the covariance matrix as stiffness matrix with respect to the finite element
discretization. We focus on the discrete spatial approximation in Section 4.1
and on the statistical approximation in Section 4.2. In Section 4.3 we review
some tail estimates for Gaussian random variables. In Section 4.4 we give
bounds on the sampling covariance error. In Section 5 we present our final
error estimates for the covariance operator reconstruction and give sufficient
conditions on the three discretization parameters to guarantee that an error
below a prescribed accuracy $\sim\varepsilon$ is achieved. We give some
concluding remarks in Section 7.

##  2 Notation and basic facts

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 non-
negative integers is denoted by $\mathbb{N}$. For any index
$\alpha\in\mathbb{N}^{d}$, $|\alpha|$ is the sum of its components. The
letters $L$, $M$ 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 $L$, $M$
and $h$, and we likewise write $A\gtrsim B$. Moreover, for any
$s\in\mathbb{N}$, $1\leq p\leq\infty$, we follow [1] and define the Sobolev
space $W^{s,p}(\mathcal{D})$ by

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

It is equipped with the norm

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

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

We now recall the classical Monte–Carlo algorithm. To this end, we consider a
probability space $(\Omega,\mathcal{F},\mathbb{P})$ and a given probability
distribution $\mu$. Let $X_{1},\cdots,X_{M}$ for $M\in\mathbb{N}$ be
independent identically distributed real-valued random variables on
$(\Omega,\mathcal{F},\mathbb{P})$ with the same associated probability
distribution $\mu$, i.e,

$\displaystyle\mathbb{P}(X^{-1}_{m}(B))=\mu(B)\quad\text{for all
}B\in\mathcal{B}(\mathbb{R})\text{ and all }1\leq m\leq M.$

Here, $\mathcal{B}(\mathbb{R})$ denotes the Borel $\sigma$-algebra. Each
random variable $X_{m}$ induces a sub $\sigma$-algebra, defined by

$\displaystyle\sigma(X_{m})=\left\\{X^{-1}_{m}(B)\ :\
B\in\mathcal{B}(\mathbb{R})\right\\}.$

Note that the $\\{\sigma(X_{m})\\}_{m=1}^{M}$ are $\mathbb{P}$-independent
because the random variables $\\{X_{m}\\}_{m=1}^{M}$ are pairwise independent,
i.e., for every finite set $\\{m_{1},\dots m_{k}\\}\subset\mathbb{N}$ and
every $B_{m_{j}}\in\sigma(X_{m_{j}})$ for $1\leq j\leq k$, we have

$\displaystyle\mathbb{P}\left(\bigcap_{j=1}^{k}B_{m_{j}}\right)=\prod_{j=1}^{k}\mathbb{P}\left(B_{m_{j}}\right).$

Let $f:\mathbb{R}\to\mathbb{R}$ be a Borel function. Then $Y_{m}=f(X_{m})$ is
also a random variable, with its associated probability distribution being

$\displaystyle\mathbb{P}((f\circ
X_{m})^{-1}(B))=\mathbb{P}(X^{-1}_{m}(f^{-1}(B)))=\mu(f^{-1}(B))\quad\text{for
all }B\in\mathcal{B}(\mathbb{R})\text{ and all }m\in\mathbb{N}.$

Therefore, the composite random variables $\\{Y_{m}\\}$ are also identically
distributed, with their generated $\sigma$-algebras defined by

$\displaystyle\sigma(Y_{m})=\left\\{(f\circ X_{m})^{-1}(B)\ :\
B\in\mathcal{B}(\mathbb{R})\right\\}=\left\\{X^{-1}_{m}(f^{-1}(B))\ :\
B\in\mathcal{B}(\mathbb{R})\right\\}\subset\sigma(X_{m}).$

Consequently, the random variables $\\{Y_{m}\\}$ are pairwise independent as a
result of the independence of the $\\{X_{m}\\}$.

Next we introduce a real-valued random variable $X$ equipped with the
probability $\mu$. Let $Y$ be a random variable defined by $Y=f(X)$. Assume
that the first two moments of $Y$ exist. Then its mean and variance are given
as

$\displaystyle\mathbb{E}[Y]$
$\displaystyle=\int_{\Omega}Y(\omega)\,\mathrm{d}\mathbb{P}(\omega)=\int_{\Omega}(f\circ
X)(\omega)\,\mathrm{d}\mathbb{P}(\omega)=\int_{\mathbb{R}}f(x)\,\mathrm{d}\mu(x)<\infty,$
$\displaystyle\mathbb{V}[Y]$
$\displaystyle=\mathbb{E}[(Y-\mathbb{E}[Y])^{2}]=\int_{\mathbb{R}}\left(f(x)-\mathbb{E}[Y]\right)^{2}\,\mathrm{d}\mu(x)<\infty.$

By definition, $Y$ and $\\{Y_{m}\\}$ have the same mean

$\displaystyle\mathbb{E}[Y_{m}]=\int_{\mathbb{R}}f(x)\,\mathrm{d}\mu(x)=\mathbb{E}[Y]$

and hence, by linearity, we get

$\displaystyle\mathbb{E}\left[\frac{1}{M}\sum_{m=1}^{M}Y_{m}\right]=\mathbb{E}[Y].$

Furthermore, we obtain

$\displaystyle\mathbb{V}\left[\frac{1}{M}\sum_{m=1}^{M}Y_{m}\right]$
$\displaystyle=\mathbb{E}\left[\left(\frac{1}{M}\sum_{m=1}^{M}Y_{m}-\mathbb{E}[Y]\right)^{2}\right]=\text{Cov}\left(\frac{1}{M}\sum_{m=1}^{M}Y_{m},\frac{1}{M}\sum_{m=1}^{M}Y_{m}\right)$
$\displaystyle=\frac{1}{M^{2}}\sum_{m,\tilde{m}=1}^{M}\text{Cov}\left(Y_{m},Y_{\tilde{m}}\right)=\frac{1}{M^{2}}\sum_{m=1}^{N}\text{Cov}\left(Y_{m},Y_{m}\right)=\frac{1}{M}\mathbb{V}[Y],$
(2.1)

where we used the definition of covariance and the fact that mixed terms
vanish due to independence. The idea of Monte–Carlo is to approximate the mean
$\mathbb{E}[Y]$ by the estimator $\hat{\mathbb{E}}_{M}[Y]$ defined by

$\displaystyle\hat{\mathbb{E}}_{M}[Y]:=\frac{1}{M}\sum_{m=1}^{M}Y_{m}=\frac{1}{M}\sum_{m=1}^{M}f(X_{m}).$

Then an application of (2.1) yields the corresponding RMSE error

RMSE
$\displaystyle:=\sqrt{\mathbb{E}\left[\left(\hat{\mathbb{E}}_{M}[Y]-\mathbb{E}[Y]\right)^{2}\right]}=\frac{1}{\sqrt{M}}\mathbb{V}\left[Y\right]^{\frac{1}{2}}.$

For any $\epsilon>0$, a combination of Chebyshev’s inequality and the equality
(2.1) shows

$\displaystyle\mathbb{P}\left(\left\\{\left|\hat{\mathbb{E}}_{M}[Y]-\mathbb{E}[Y]\right|>\epsilon\right\\}\right)\leq\epsilon^{-2}\mathbb{V}\left[\hat{\mathbb{E}}_{M}[Y]\right]=\frac{1}{M\epsilon^{2}}\mathbb{V}(Y).$

Taking
$\epsilon:=\frac{1}{\sqrt{M\delta}}\left(\mathbb{V}(Y)\right)^{\frac{1}{2}}$
for some $\delta>0$, we obtain

$\displaystyle\mathbb{P}\left(\left\\{\left|\hat{\mathbb{E}}_{M}[Y]-\mathbb{E}[Y]\right|>\frac{1}{\sqrt{M\delta}}\left(\mathbb{V}(Y)\right)^{\frac{1}{2}}\right\\}\right)\leq\delta.$

##  3 Karhunen-Loève expansion: Continuous level

This section is concerned with the Karhunen-Loève expansion of the centered
Gaussian random field $\kappa$. Let $\mathcal{L}^{(d)}(\boldsymbol{x})$ be the
Lebesgue measure on the physical domain $\mathcal{D}$. For the sake of
simplicity, $L^{2}(\mathcal{D})$ and $L^{2}(\Omega)$ are short for
$L^{2}(\mathcal{D};\mathrm{d}\mathcal{L}^{(d)}(\boldsymbol{x}))$ and
$L^{2}(\Omega;\mathrm{d}\mathbb{P})$. We denote the associated integral
operator $\mathcal{S}:L^{2}(\mathcal{D})\rightarrow L^{2}(\Omega)$ by

$\displaystyle\mathcal{S}v=\int_{\mathcal{D}}\kappa(\boldsymbol{\omega},\cdot)v\,\mathrm{d}\mathcal{L}^{(d)},$
(3.1)

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

$\displaystyle\mathcal{S}^{*}v=\int_{\Omega}\kappa(\boldsymbol{\omega},\cdot)v\,\mathrm{d}\mathbb{P}(\boldsymbol{\omega}).$
(3.2)

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

$\displaystyle
R(\boldsymbol{x},\boldsymbol{x}^{\prime})=\int_{\Omega}\kappa(\boldsymbol{\omega},\boldsymbol{x})\kappa(\boldsymbol{\omega},\boldsymbol{x}^{\prime})\,\mathrm{d}\mathbb{P}(\boldsymbol{\omega})=\mathbb{E}\left[\kappa(\cdot,\boldsymbol{x})\kappa(\cdot,\boldsymbol{x}^{\prime})\right]$
(3.3)

Moreover, for any $v\in L^{2}(\mathcal{D})$, we have

$\displaystyle\mathcal{R}v(\boldsymbol{x})=\int_{\mathcal{D}}R(\boldsymbol{x},\boldsymbol{x}^{\prime})v(\boldsymbol{x}^{\prime})\,\mathrm{d}\mathcal{L}^{(d)}(\boldsymbol{x}^{\prime})=\int_{\mathcal{D}}\int_{\Omega}\kappa(\boldsymbol{\omega},\boldsymbol{x})\kappa(\boldsymbol{\omega},\boldsymbol{x}^{\prime})\,\mathrm{d}\mathbb{P}(\boldsymbol{\omega})\,v(\boldsymbol{x}^{\prime})\mathrm{d}\mathcal{L}^{(d)}(\boldsymbol{x}^{\prime}).$
(3.4)

The standard spectral theory for compact operators [19] 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_{\ell}\\}_{\ell=1}^{\infty}$ be the
sequence of eigenvalues (with multiplicity counted) associated to
$\mathcal{R}$, which are ordered non-increasingly, and let
$\\{\phi_{\ell}\\}_{\ell=1}^{\infty}$ be the corresponding eigenfunctions that
are orthonormal in $L^{2}(\mathcal{D})$. Furthermore, for any
$\lambda_{\ell}\neq 0$, define

$\psi_{\ell}(\boldsymbol{\omega})=\frac{1}{\sqrt{\lambda_{\ell}}}\int_{\mathcal{D}}\kappa(\boldsymbol{\omega},\boldsymbol{x})\phi_{\ell}\,\mathrm{d}\mathcal{L}^{(d)}(\boldsymbol{x}).$
(3.5)

One can verify that the sequence $\\{\psi_{\ell}\\}_{{\ell}=1}^{\infty}$ is
uncorrelated and orthonormal in $L^{2}(\Omega)$ and therefore,
$\\{\psi_{\ell}\\}_{{\ell}=1}^{\infty}$ are i.i.d normal random functions.
Note that the sequence $\\{\lambda_{\ell}\\}_{{\ell}=1}^{\infty}$ can be
characterized by the so-called approximation numbers (cf. [14, Section
2.3.1]). They are defined by

$\displaystyle\lambda_{\ell}=\inf\\{\left\|\mathcal{R}-\mathcal{T}\right\|_{\mathcal{B}(L^{2}(\mathcal{D}))}:\mathcal{T}\in\mathfrak{F}(L^{2}(\mathcal{D})),{\text{rank}}(\mathcal{T})<\ell\\},$
(3.6)

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

The KL expansion of the bivariate function $\kappa$ then refers to the
expression

$\kappa(\boldsymbol{\omega},\boldsymbol{x})=\sum\limits_{\ell=1}^{\infty}\sqrt{\lambda_{\ell}}\phi_{\ell}(\boldsymbol{x})\psi_{\ell}(\boldsymbol{\omega}),$
(3.7)

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

###  3.1 $L$-term truncation in case of continuous Karhunen-Loève expansion

Now we will truncate the KL expansion and discuss the resulting error. The
studies on the $L$-term KL approximation to random fields are extensive. For
example, in [15], decay rates for the eigenvalues of covariance kernels
possessing certain regularity were considered and generalized fast multipole
methods to solve the associated eigenvalue problems were studied. The robust
computation of eigenvalues for smooth covariance kernels was treated in [17].
A comparison of $L$-term KL truncation and the sparse grids approximation was
given in [10].

The result of this section is based on our recent paper [11], which proves a
sharp eigenvalue decay rate under a mild assumption on the regularity of the
bivariate function $\kappa$ in the physical domain. To this end, we recall
Assumption 1.1 concerning the regularity of $\kappa$. Its weaker variant (1.8)
states that $K\in L^{\infty}(\Omega,H^{s}(\mathcal{D}))$ for some $s\geq 0$.
Its stronger variant (1.9) implies that the kernel belongs to
$H^{s}(\mathcal{D})\times H^{s}(\mathcal{D})$. The following eigenvalue decay
estimate [11, Theorems 3.2, 3.3 and 3.4] will be used repeatedly.

######  Theorem 3.1.

Let (1.9) from Assumption 1.1 hold. Then, there holds

$\displaystyle{{\lambda_{\ell}}}$ $\displaystyle\leq
C_{\ref{thm:truncationError}}\ell^{-\frac{2s}{d}-1}\quad\text{ for all
$\ell\geq 1$ and }$ (3.8)
$\displaystyle\Big{\|}{\sum\limits_{\ell>L}\sqrt{\lambda_{\ell}}\phi_{\ell}\psi_{\ell}}\Big{\|}_{L^{2}(\Omega\times\mathcal{D})}$
$\displaystyle\leq
C_{\ref{thm:truncationError}}^{1/2}\sqrt{\frac{d}{2s}}(L+1)^{-\frac{s}{d}}\quad\text{
when $L$ is sufficiently large.}$ (3.9)

Here, we use the constant
$C_{\ref{thm:truncationError}}:=\mathrm{diam}(\mathcal{D})^{2s}C_{\rm
em}(d,s)C_{\rm{ext}}(D,s)\left\|\kappa\right\|_{L^{2}(\Omega,H^{s}\left(\mathcal{D}\right))}^{2}$,
where $C_{\rm{em}}(d,s)$ denotes an embedding constant between certain Lorentz
sequence spaces, and $C_{\rm{ext}}(D,s)$ is a constant which depends only on
$\mathcal{D}$ and $s$.

The next lemma provides a regularity result for the eigenfunctions
$\\{\phi_{\ell}\\}_{\ell=1}^{\infty}$ which follows with operator
interpolation from [11, Lemma 3.1 & Remark after Assumption 3.1]

######  Lemma 3.1 (Regularity of the eigenfunctions
$\\{\phi_{\ell}\\}_{\ell=1}^{\infty}$).

Let (1.9) from Assumption 1.1 be valid. Then for all $0\leq\beta\leq 1$, there
holds

$\displaystyle\left\|\phi_{\ell}\right\|_{H^{\beta
s}\left(\mathcal{D}\right)}\leq C(\mathcal{D},d,s)\ell^{\frac{\beta
s}{d}}\quad\text{ when $\ell$ is sufficiently large}.$ (3.10)

###### Proof.

Assumption (1.9) states that $\kappa\in
L^{\infty}(\Omega,H^{s}(\mathcal{D}))\text{ for some }s\geq 0$. As
$\int_{\Omega}1\,\mathrm{d}\mathbb{P}(\boldsymbol{\omega})=1$, we deduce
$\kappa\in L^{2}(\Omega,H^{s}(\mathcal{D}))$ and

$\displaystyle\left\|\kappa\right\|^{2}_{L^{2}(\Omega,H^{s}(\mathcal{D}))}=\sum_{\ell=1}^{\infty}\lambda_{\ell}\left\|\phi_{\ell}\right\|^{2}_{H^{s}(\mathcal{D})}\leq
C_{\ref{thm:truncationError}}\sum_{\ell=1}^{\infty}\ell^{-1-\epsilon}\ell^{-\frac{2s}{d}+\epsilon}\left\|\phi_{\ell}\right\|^{2}_{H^{s}(\mathcal{D})}$

for all $\epsilon>0$. Considering the limit $\epsilon\to 0^{+}$, we infer that
$\left\|\phi_{\ell}\right\|_{H^{s}(\mathcal{D})}\lesssim\ell^{-\frac{s}{d}}$
for $\ell$ large. Since
$\left\|\phi_{\ell}\right\|_{H^{0}(\mathcal{D})}=\left\|\phi_{\ell}\right\|_{L^{2}(\mathcal{D})}=1$,
the statement follows via an interpolation argument. ∎

Here, $C(\mathcal{D},d,s)$ denotes a positive constant depending only on
$\mathcal{D}$, $d$ and $s$. We also get a bound in the uniform norm by using
Sobolev’s embedding theorem for $s=\frac{d}{2}+\epsilon\geq\frac{d}{2}$, i.e.,
we have

$\displaystyle\left\|\phi_{\ell}\right\|_{C\left(\mathcal{D}\right)}\leq\left\|\phi_{\ell}\right\|_{H^{s}\left(\mathcal{D}\right)}\leq
C(\mathcal{D},d,s)\ell^{\frac{s}{d}}=C(\mathcal{D},d)\ell^{\frac{1}{2}+\epsilon}\quad\text{
when $\ell$ is sufficiently large},$ (3.11)

where we use $\epsilon>0$ for an arbitrary small number which may change from
line to line. The eigenfunctions $\\{\phi_{\ell}\\}_{\ell=1}^{\infty}$ are
optimal in the sense that the mean-square error resulting from a finite-rank
approximation of $\kappa$ is minimized [9]. Thus, the eigenfunctions indeed
minimize the truncation error in the $L^{2}$-sense, i.e.

$\displaystyle\min\limits_{\begin{subarray}{c}\\{c_{\ell}(\boldsymbol{x})\\}_{\ell=1}^{L}\subset
L^{2}(\mathcal{D})\\\ \\{c_{\ell}(\boldsymbol{x})\\}_{\ell=1}^{L}\text{
orthonormal}\end{subarray}}\left\|\kappa(\boldsymbol{\omega},\boldsymbol{x})-\sum\limits_{\ell=1}^{L}\left(\int_{\mathcal{D}}\kappa(\boldsymbol{\omega},\boldsymbol{x}^{\prime})c_{\ell}(\boldsymbol{x}^{\prime})\,\mathrm{d}\mathcal{L}^{(d)}(\boldsymbol{x}^{\prime})\right)c_{\ell}(\boldsymbol{x})\right\|_{L^{2}\left(\Omega\times\mathcal{D}\right)}=\sqrt{\sum\limits_{\ell>L}\lambda_{\ell}}.$
(3.12)

##  4 Approximation of the covariance matrix

###  4.1 Spatial approximation

In this section, we investigate the influence of the spatial discretization.
Recall that the random field $\kappa^{(h)}$ defined in (1.13) is an
approximation to the random field $\kappa$. We assume for the moment to have
access to its covariance function

$\displaystyle
R^{(h)}(\boldsymbol{x},\boldsymbol{x}^{\prime}):=\int_{\Omega}\kappa^{(h)}(\boldsymbol{\omega},\boldsymbol{x})\kappa^{(h)}(\boldsymbol{\omega},\boldsymbol{x}^{\prime})\,\mathrm{d}\mathbb{P}(\boldsymbol{\omega})=\mathbb{E}\left[\kappa^{(h)}(\cdot,\boldsymbol{x})\kappa^{(h)}(\cdot,\boldsymbol{x}^{\prime})\right]$
(4.1)

and the associated integral operator

$\displaystyle\mathcal{R}^{(h)}f(\boldsymbol{x}):=\int_{\mathcal{D}}f(\boldsymbol{x}^{\prime})R^{(h)}(\boldsymbol{x},\boldsymbol{x}^{\prime})\,\mathrm{d}\boldsymbol{x}^{\prime}.$

Then we have the following result:

######  Lemma 4.1 (Semi-discrete spatial approximation error estimate).

Let $R\in L^{2}(\mathcal{D}\times\mathcal{D})$ be defined as in (3.3) and let
its numerical approximation $R_{h}$ be defined as in (4.1). Then it holds

$\displaystyle\left\|R^{(h)}-R\right\|_{C(\mathcal{D}\times\mathcal{D})}\leq
C_{\Pi_{\mathcal{V}_{h};L^{2}}}C_{\infty}C_{H^{s}}h^{s-\frac{d}{2}}.$ (4.2)

Furthermore, for the associated integral operators $\mathcal{R}$ and
$\mathcal{R}_{h}$, it holds

$\displaystyle\left\|\mathcal{R}-\mathcal{R}^{(h)}\right\|_{\mathcal{B}(L^{2}(\mathcal{D}))}\leq
C_{\mathcal{V}_{h};L^{2}}(C_{\infty}C_{H^{s}}+C_{\infty}^{(h)}C_{H^{s}}^{(h)})h^{s}.$
(4.3)

###### Proof.

Using $ab-a_{h}b_{h}=(a-a_{h})b+a_{h}(b-b_{h})$, we obtain

$\displaystyle\left\|R^{(h)}-R\right\|_{C(\mathcal{D}\times\mathcal{D})}$
$\displaystyle\leq\sup_{\boldsymbol{x},\boldsymbol{x}^{\prime}\in\mathcal{D}}\left|\int_{\Omega}\left(\kappa^{(h)}(\boldsymbol{\omega},\boldsymbol{x})-\kappa(\boldsymbol{\omega},\boldsymbol{x})\right)\kappa^{(h)}(\boldsymbol{\omega},\boldsymbol{x}^{\prime})\,\mathrm{d}\mathbb{P}(\boldsymbol{\omega})\right|$
$\displaystyle+\sup_{\boldsymbol{x},\boldsymbol{x}^{\prime}\in\mathcal{D}}\left|\int_{\Omega}\left(\kappa^{(h)}(\boldsymbol{\omega},\boldsymbol{x}^{\prime})-\kappa(\boldsymbol{\omega},\boldsymbol{x}^{\prime})\right)\kappa(\boldsymbol{\omega},\boldsymbol{x})\,\mathrm{d}\mathbb{P}(\boldsymbol{\omega})\right|$
$\displaystyle\leq\int_{\Omega}\left\|\kappa^{(h)}(\boldsymbol{\omega},\cdot)-\kappa(\boldsymbol{\omega},\cdot)\right\|_{C(\mathcal{D})}\left\|\kappa^{(h)}(\boldsymbol{\omega},\cdot)\right\|_{C(\mathcal{D})}\,\mathrm{d}\mathbb{P}(\boldsymbol{\omega})$
$\displaystyle+\int_{\Omega}\left\|\kappa^{(h)}(\boldsymbol{\omega},\cdot)-\kappa(\boldsymbol{\omega},\cdot)\right\|_{C(\mathcal{D})}\left\|\kappa(\boldsymbol{\omega},\cdot)\right\|_{C(\mathcal{D})}\,\mathrm{d}\mathbb{P}(\boldsymbol{\omega})$
$\displaystyle\leq
C_{\Pi_{\mathcal{V}_{h};L^{2}}}h^{s-\frac{d}{2}}\int_{\Omega}\left\|\kappa(\boldsymbol{\omega},\cdot)\right\|_{H^{s}\left(\mathcal{D}\right)}\left\|\kappa^{(h)}(\boldsymbol{\omega},\cdot)\right\|_{C(\mathcal{D})}\,\mathrm{d}\mathbb{P}(\boldsymbol{\omega})$
$\displaystyle+C_{\Pi_{\mathcal{V}_{h};L^{2}}}h^{s-\frac{d}{2}}\int_{\Omega}\left\|\kappa(\boldsymbol{\omega},\cdot)\right\|_{H^{s}\left(\mathcal{D}\right)}\left\|\kappa(\boldsymbol{\omega},\cdot)\right\|_{C(\mathcal{D})}\,\mathrm{d}\mathbb{P}(\boldsymbol{\omega}),$

where we used (1) in the last inequality. This, together with the assumption
(1.5), leads to

$\displaystyle\left\|R^{(h)}-R\right\|_{C(\mathcal{D}\times\mathcal{D})}$
$\displaystyle\leq
C_{\Pi_{\mathcal{V}_{h};L^{2}}}h^{s-\frac{d}{2}}\left(\int_{\Omega}\left\|\kappa(\boldsymbol{\omega},\cdot)\right\|_{H^{s}\left(\mathcal{D}\right)}^{2}\,\mathrm{d}\mathbb{P}(\boldsymbol{\omega})\right)^{\frac{1}{2}}\left(\int_{\Omega}\left\|\kappa(\boldsymbol{\omega},\cdot)\right\|_{C(\mathcal{D})}^{2}\,\mathrm{d}\mathbb{P}(\boldsymbol{\omega})\right)^{\frac{1}{2}}$
$\displaystyle\leq
C_{\Pi_{\mathcal{V}_{h};L^{2}}}C_{H^{s}}C_{\infty}h^{s-\frac{d}{2}}.$

This proves the first assertion (4.2). Now note that

$\displaystyle\left\|\mathcal{R}-\mathcal{R}^{(h)}\right\|_{\mathcal{B}(L^{2}(\mathcal{D}))}:=\sup_{\genfrac{}{}{0.0pt}{}{v\in
L^{2}(\mathcal{D})}{\left\|v\right\|_{L^{2}\left(\mathcal{D}\right)}=1}}\left\|\mathcal{R}v-\mathcal{R}^{(h)}v\right\|_{L^{2}\left(\mathcal{D}\right)}\leq\left\|R-R^{(h)}\right\|_{L^{2}\left(\mathcal{D}\times\mathcal{D}\right)}.$

Using (3.3) and (4.1), we observe

$\displaystyle\left\|R-R^{(h)}\right\|_{L^{2}\left(\mathcal{D}\times\mathcal{D}\right)}^{2}=\int_{\mathcal{D}}\int_{\mathcal{D}}\left(R(\boldsymbol{x},\boldsymbol{x}^{\prime})-R^{(h)}(\boldsymbol{x},\boldsymbol{x}^{\prime})\right)^{2}\,\mathrm{d}(\mathcal{L}^{(d)}\otimes\mathcal{L}^{(d)})(\boldsymbol{x},\boldsymbol{x}^{\prime})$
$\displaystyle=\int_{\mathcal{D}}\int_{\mathcal{D}}\left(\int_{\Omega}\kappa(\boldsymbol{\omega},\boldsymbol{x})\kappa(\boldsymbol{\omega},\boldsymbol{x}^{\prime})-\kappa^{(h)}(\boldsymbol{\omega},\boldsymbol{x})\kappa^{(h)}(\boldsymbol{\omega},\boldsymbol{x}^{\prime})\,\mathrm{d}\mathbb{P}(\boldsymbol{\omega})\right)^{2}\,\mathrm{d}(\mathcal{L}^{(d)}\otimes\mathcal{L}^{(d)})(\boldsymbol{x},\boldsymbol{x}^{\prime})$
$\displaystyle=\int_{\mathcal{D}}\int_{\mathcal{D}}\left(\mathbb{E}\left[\kappa(\cdot,\boldsymbol{x})\kappa(\cdot,\boldsymbol{x}^{\prime})-\kappa^{(h)}(\cdot,\boldsymbol{x})\kappa^{(h)}(\cdot,\boldsymbol{x}^{\prime})\right]\right)^{2}\,\mathrm{d}(\mathcal{L}^{(d)}\otimes\mathcal{L}^{(d)})(\boldsymbol{x},\boldsymbol{x}^{\prime}).$

Thus, we obtain

$\displaystyle\left\|R-R^{(h)}\right\|_{L^{2}\left(\mathcal{D}\times\mathcal{D}\right)}=\left\|\mathbb{E}\left[\kappa\kappa-\kappa^{(h)}\kappa^{(h)}\right]\right\|_{L^{2}\left(\mathcal{D}\times\mathcal{D}\right)}.$

We get by an application of the Jensen’s inequality, i.e.,
$f(\mathbb{E}{X})\leq\mathbb{E}{f(X)}$ with the convex function
$f(X)=\left\|X\right\|_{L^{2}\left(\mathcal{D}\times\mathcal{D}\right)}$, the
bound

$\displaystyle\left\|{R}-{R}^{(h)}\right\|_{L^{2}\left(\mathcal{D}\times\mathcal{D}\right)}$
$\displaystyle\leq\mathbb{E}\left[\left\|{\kappa\kappa-\kappa^{(h)}\kappa^{(h)}}\right\|_{L^{2}\left(\mathcal{D}\times\mathcal{D}\right)}\right]$
$\displaystyle\leq\mathbb{E}\left[\left\|(\kappa-\kappa^{(h)})\kappa\right\|_{L^{2}\left(\mathcal{D}\times\mathcal{D}\right)}+\left\|\kappa^{(h)}(\kappa-\kappa^{(h)})\right\|_{L^{2}\left(\mathcal{D}\times\mathcal{D}\right)}\right]$
$\displaystyle\leq
C_{\mathcal{V}_{h};L^{2}}(C_{\infty}C_{H^{s}}+C_{\infty}^{(h)}C_{H^{s}}^{(h)})h^{s}$

where we used (1.11) in the last step. This proves the second assertion (4.3).
∎

We are interested in solving the eigenvalue problem for the operator
$\mathcal{R}^{(h)}$, i.e., we consider finding
$(\phi_{\ell}^{(h)},\lambda^{(h)}_{\ell})\in\mathcal{V}_{h}\times\mathbb{R}$
such that

$\displaystyle\left(\mathcal{R}^{(h)}\phi_{\ell}^{(h)},v^{(h)}\right)_{\mathcal{D}}=\lambda^{(h)}_{\ell}\left(\phi_{\ell}^{(h)},v^{(h)}\right)_{\mathcal{D}}\quad\text{for
all }v^{(h)}\in\mathcal{V}_{h}=\text{span}\left\\{\theta_{k}^{(h)}\ :\ 1\leq
k\leq Q_{h}\right\\}.$ (4.4)

Using
$\phi_{\ell}^{(h)}:=\boldsymbol{\Phi}^{(h)}_{\ell}\cdot\boldsymbol{\theta}^{(h)}=\sum_{k=1}^{Q_{h}}\phi_{\ell;k}\theta_{k}^{(h)}\in\mathcal{V}_{h}$,
we derive for fixed but arbitrary
$v^{(h)}=\boldsymbol{V}^{(h)}\cdot\boldsymbol{\theta}^{(h)}\in\mathcal{V}_{h}$,
i.,e, for fixed but arbitrary $\boldsymbol{V}^{(h)}\in\mathbb{R}^{Q_{h}}$, the
identity

$\displaystyle\left(\mathcal{R}^{(h)}\phi_{\ell}^{(h)},v^{(h)}\right)_{\mathcal{D}}=\int_{\mathcal{D}}\int_{\mathcal{D}}\int_{\Omega}\kappa^{(h)}(\boldsymbol{\omega},\boldsymbol{x})\kappa^{(h)}(\boldsymbol{\omega},\boldsymbol{x}^{\prime})\,\mathrm{d}\mathbb{P}(\boldsymbol{\omega})\phi_{\ell}^{(h)}(\boldsymbol{x})\,\mathrm{d}\mathcal{L}^{(d)}(\boldsymbol{x})v^{(h)}(\boldsymbol{x}^{\prime})\,\mathrm{d}\mathcal{L}^{(d)}(\boldsymbol{x}^{\prime})$
$\displaystyle=\int_{\mathcal{D}}\int_{\mathcal{D}}\int_{\Omega}\boldsymbol{K}^{(h)}(\boldsymbol{\omega})\cdot\boldsymbol{\theta}(\boldsymbol{x})\boldsymbol{K}^{(h)}(\boldsymbol{\omega})\cdot\boldsymbol{\theta}(\boldsymbol{x}^{\prime})\,\mathrm{d}\mathbb{P}(\boldsymbol{\omega})\boldsymbol{\Phi}_{\ell}^{(h)}\cdot\boldsymbol{\theta}(\boldsymbol{x})\,\mathrm{d}\mathcal{L}^{(d)}(\boldsymbol{x})\boldsymbol{V}^{(h)}\cdot\boldsymbol{\theta}(\boldsymbol{x}^{\prime})\,\mathrm{d}\mathcal{L}^{(d)}(\boldsymbol{x}^{\prime})$
$\displaystyle=\int_{\Omega}\left(\int_{\mathcal{D}}\boldsymbol{K}^{(h)}(\boldsymbol{\omega})\cdot\boldsymbol{\theta}(\boldsymbol{x})\boldsymbol{\Phi}_{\ell}^{(h)}\cdot\boldsymbol{\theta}(\boldsymbol{x})\,\mathrm{d}\mathcal{L}^{(d)}(\boldsymbol{x})\right)\left(\int_{\mathcal{D}}\boldsymbol{K}^{(h)}(\boldsymbol{\omega})\cdot\boldsymbol{\theta}(\boldsymbol{x}^{\prime})\boldsymbol{V}^{(h)}\cdot\boldsymbol{\theta}(\boldsymbol{x}^{\prime})\,\mathrm{d}\mathcal{L}^{(d)}(\boldsymbol{x}^{\prime})\right)\,\mathrm{d}\mathbb{P}(\boldsymbol{\omega})$
$\displaystyle=\int_{\Omega}\boldsymbol{K}^{(h)}(\boldsymbol{\omega})\cdot\boldsymbol{M}^{(h)}\boldsymbol{\Phi}_{\ell}^{(h)}\boldsymbol{K}^{(h)}(\boldsymbol{\omega})\cdot\boldsymbol{M}^{(h)}\boldsymbol{V}^{(h)}\,\mathrm{d}\mathbb{P}(\boldsymbol{\omega})$
$\displaystyle=\int_{\Omega}\left(\boldsymbol{V}^{(h)}\right)^{T}\boldsymbol{M}^{(h)}\boldsymbol{K}^{(h)}(\boldsymbol{\omega})\otimes\boldsymbol{K}^{(h)}(\boldsymbol{\omega})\boldsymbol{M}^{(h)}\boldsymbol{\Phi}_{\ell}^{(h)}\,\mathrm{d}\mathbb{P}(\boldsymbol{\omega})=\left(\boldsymbol{V}^{(h)}\right)^{T}\boldsymbol{M}^{(h)}\boldsymbol{\Sigma}_{\boldsymbol{K}^{(h)}}\boldsymbol{M}^{(h)}\boldsymbol{\Phi}_{\ell}^{(h)}.$

Here $\boldsymbol{M}^{(h)}$ and $\boldsymbol{S}^{(h)}$ are the mass matrix and
the stiffness matrix defined by

$\displaystyle\boldsymbol{M}^{(h)}=(\boldsymbol{\theta}^{(h)},{\boldsymbol{\theta}^{(h)}}^{T})_{\mathcal{D}}\quad\text{and}\quad\boldsymbol{S}^{(h)}=\boldsymbol{M}^{(h)}\boldsymbol{\Sigma}_{\boldsymbol{K}^{(h)}}\boldsymbol{M}^{(h)}.$
(4.5)

Note that the symmetry of the mass matrix $\boldsymbol{M}^{(h)}$ implies the
symmetry of $\boldsymbol{S}^{(h)}$.

Thus, the corresponding matrix form of the generalized eigenvalue problem
(4.4) is to seek the eigenpair
$(\boldsymbol{\Phi}^{(h)}_{\ell},\lambda^{(h)}_{\ell})\in\mathbb{R}^{Q_{h}}\times\mathbb{R}$
for $\ell=1,\cdots,Q_{h}$, satisfying

$\displaystyle\boldsymbol{S}^{(h)}\boldsymbol{\Phi}^{(h)}_{\ell}=\lambda^{(h)}_{\ell}\boldsymbol{M}^{(h)}\boldsymbol{\Phi}^{(h)}_{\ell}.$
(4.6)

This is a generalized eigenvalue problem. We now transform it to a
conventional eigenvalue problem. To this end, observe that the mass matrix
$\boldsymbol{M}^{(h)}$ is symmetric positive-definite, and thus possesses a
Cholesky factorization

$\displaystyle\boldsymbol{M}^{(h)}=:\boldsymbol{L}^{(h)}\left(\boldsymbol{L}^{(h)}\right)^{T}\in\mathbb{R}^{Q_{h}\times
Q_{h}}.$

Now let $\widetilde{\boldsymbol{S}}^{(h)}$ be a symmetric positive-definite
matrix defined by

$\displaystyle\widetilde{\boldsymbol{S}}^{(h)}:=\left(\boldsymbol{L}^{(h)}\right)^{T}\boldsymbol{\Sigma}_{\boldsymbol{K}^{(h)}}\boldsymbol{L}^{(h)}.$

Then we can solve for
$(\widetilde{\boldsymbol{\Phi}}^{(h)}_{\ell},\lambda^{(h)}_{\ell})\in\mathbb{R}^{Q_{h}}\times\mathbb{R}$
for $\ell=1,\cdots,Q_{h}$, satisfying

$\displaystyle\widetilde{\boldsymbol{S}}^{(h)}\widetilde{\boldsymbol{\Phi}}^{(h)}_{\ell}=\lambda^{(h)}_{\ell}\widetilde{\boldsymbol{\Phi}}^{(h)}_{\ell}.$
(4.7)

Suppose now that
$(\widetilde{\boldsymbol{\Phi}}^{(h)}_{\ell},\lambda^{(h)}_{\ell})\in\mathbb{R}^{Q_{h}}\times\mathbb{R}$
is an eigenpair of (4.7). Using

$\displaystyle\boldsymbol{S}^{(h)}=\boldsymbol{M}^{(h)}\boldsymbol{\Sigma}_{\boldsymbol{K}^{(h)}}\boldsymbol{M}^{(h)}=\boldsymbol{L}^{(h)}\widetilde{\boldsymbol{S}}^{(h)}\left(\boldsymbol{L}^{(h)}\right)^{T}$

we obtain

$\displaystyle\boldsymbol{S}^{(h)}\left(\boldsymbol{L}^{(h)}\right)^{-T}\widetilde{\boldsymbol{\Phi}}^{(h)}_{\ell}=\boldsymbol{L}^{(h)}\widetilde{\boldsymbol{S}}^{(h)}\widetilde{\boldsymbol{\Phi}}^{(h)}_{\ell}=\lambda^{(h)}\boldsymbol{L}^{(h)}\widetilde{\boldsymbol{\Phi}}^{(h)}_{\ell}=\lambda^{(h)}\boldsymbol{M}^{(h)}\left(\boldsymbol{L}^{(h)}\right)^{-T}\widetilde{\boldsymbol{\Phi}}^{(h)}_{\ell}.$

Thus, we have shown the equivalence

$\displaystyle(\widetilde{\boldsymbol{\Phi}}^{(h)}_{\ell},\lambda^{(h)}_{\ell})\in\mathbb{R}^{Q_{h}}\times\mathbb{R}\text{
is an eigenpair of \eqref{genalizedeigenvalue3}}\Longleftrightarrow$
$\displaystyle\left(\boldsymbol{\Phi}^{(h)}_{\ell}:=\left(\boldsymbol{L}^{(h)}\right)^{-T}\widetilde{\boldsymbol{\Phi}}^{(h)}_{\ell},\lambda^{(h)}_{\ell}\right)\in\mathbb{R}^{Q_{h}}\times\mathbb{R}\text{
is an eigenpair of \eqref{genalizedeigenvalue2}.}$

Note that, with
$\phi_{\ell}^{(h)}=\boldsymbol{\Phi}^{(h)}_{\ell}\cdot\boldsymbol{\theta}^{(h)}$,
we obtain a Mercer-like representation

$\displaystyle R^{(h)}(\boldsymbol{x},\boldsymbol{x}^{\prime})$
$\displaystyle=\sum_{k=1}^{Q_{h}}\lambda^{(h)}_{k}\phi_{k}^{(h)}(\boldsymbol{x})\phi_{k}^{(h)}(\boldsymbol{x}^{\prime})\approx
R^{(L;h)}(\boldsymbol{x},\boldsymbol{x}^{\prime}):=\sum_{k=1}^{L}\lambda^{(h)}_{k}\phi_{k}^{(h)}(\boldsymbol{x})\phi_{k}^{(h)}(\boldsymbol{x}^{\prime}),$
(4.8)

for $1\leq L\leq Q_{h}$.

We now give a classical error analysis by means of the approximation theory of
conforming finite element methods. To this end, consider the error operator

$\displaystyle\mathcal{E}^{(h)}=\mathcal{R}-\mathcal{R}^{(h)}:L^{2}(\mathcal{D})\to
L^{2}(\mathcal{D}),$ (4.9)

which is a self-adjoint operator on $L^{2}(\mathcal{D})$. We have the
following error representation.

######  Lemma 4.2.

The error operator $\mathcal{E}^{(h)}$ has the property

$\left(\mathcal{E}^{(h)}v,v\right)_{\mathcal{D}}=\left(v,(I-\Pi_{\mathcal{V}_{h};L^{2}})\mathcal{R}(I+\Pi_{\mathcal{V}_{h};L^{2}})v\right)_{\mathcal{D}}\quad\text{
for all }\quad v\in L^{2}(D),$

where we use the operator
$\Pi_{\mathcal{V}_{h};L^{2}}:L^{2}(\mathcal{D})\to\mathcal{V}_{h}$ from
(1.12).

###### Proof.

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

$\displaystyle\left(\mathcal{E}^{(h)}v,v\right)_{\mathcal{D}}$
$\displaystyle=\left(\left(\mathcal{R}-\Pi_{\mathcal{V}_{h};L^{2}}\mathcal{R}\Pi_{\mathcal{V}_{h};L^{2}}\right)v,v\right)_{\mathcal{D}}+\left(\mathcal{R}\Pi_{\mathcal{V}_{h};L^{2}}v-\Pi_{\mathcal{V}_{h};L^{2}}\mathcal{R}v,v\right)_{\mathcal{D}}$
$\displaystyle=\left(\left(I-\Pi_{\mathcal{V}_{h};L^{2}}\right)\mathcal{R}\left(I+\Pi_{\mathcal{V}_{h};L^{2}}\right)v,v\right)_{\mathcal{D}},$

which gives the assertion. ∎

A direct consequence of Lemma 4.2 is the upper bound for the operator norm of
$\mathcal{E}^{(h)}$

$\displaystyle\left\|\mathcal{E}^{(h)}\right\|_{\mathcal{B}(L^{2}(\mathcal{D}))}\leq
2C_{\Pi_{\mathcal{V}_{h};L^{2}}}h^{s}\left\|\mathcal{R}\right\|_{\mathcal{B}(L^{2}(\mathcal{D}),H^{s}(\mathcal{D}))},$
(4.10)

where we employed (1.12).

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

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

Let (1.8) from Assumption 1.1 hold. Then there are constants $C_{1}$, $C_{2}$
and $h_{0}$ such that

$\left|\lambda_{\ell}^{(h)}-\lambda_{\ell}\right|\leq
C_{1}\lambda_{\ell}^{-1}h^{2s}\quad\text{ for all }\quad 0<h\leq h_{0}.$

Furthermore, the eigenvectors $\\{\phi_{\ell}^{(h)}\\}_{\ell=1}^{Q_{h}}$ of
$\mathcal{R}^{(h)}$ can be selected such that

$\left\|\phi_{\ell}^{(h)}-\phi_{\ell}\right\|_{L^{2}\left(\mathcal{D}\right)}\leq
C_{2}\lambda_{\ell}^{-1}h^{s}\quad\text{ for all }\quad 0<h\leq h_{0}.$

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

###### Proof.

The proof follows from [2, Theorem 9.1]. ∎

A consequence of Proposition 4.1 is that, for any $\epsilon>0$ and for the
Lebesgue measure $\mathcal{L}^{(d)}$, there holds the estimate

$\displaystyle\mathcal{L}^{(d)}\left(\left\\{\
\left|\phi_{\ell}^{(h)}-\phi_{\ell}\right|\geq\epsilon\right\\}\right)\leq\epsilon^{-2}\int_{\mathcal{D}}\left|\phi_{\ell}^{(h)}-\phi_{\ell}\right|^{2}\,\mathrm{d}\mathcal{L}^{(d)}=\epsilon^{-2}\left\|\phi_{\ell}^{(h)}-\phi_{\ell}\right\|_{L^{2}\left(\mathcal{D}\right)}^{2}\leq
C_{2}\epsilon^{-2}\lambda_{\ell}^{-2}h^{2s}.$ (4.11)

###  4.2 Approximation of the covariance operator from samples

We will approximate the spectral decomposition (4.8), which is based on the
unknown stiffness matrix $\boldsymbol{S}^{(h)}$ by approximating the true
covariance matrix
$\boldsymbol{\Sigma}_{\boldsymbol{K}^{(h)}}\in\mathbb{R}^{Q_{h}\times Q_{h}}$.
Let us assume for the moment that we have such an estimate
$\boldsymbol{\Sigma}^{(h;M)}_{\boldsymbol{K}^{(h)}}\in\mathbb{R}^{Q_{h}\times
Q_{h}},$ where $\boldsymbol{\Sigma}^{(h;M)}_{\boldsymbol{K}^{(h)}}$ is also
symmetric. Then by repeating the procedure for the derivation (4.8) in an
analogous way, we obtain the approximations

$\displaystyle\boldsymbol{S}^{(h;M)}$
$\displaystyle:=\left(\boldsymbol{M}^{(h)}\right)^{T}\boldsymbol{\Sigma}^{(h;M)}_{\boldsymbol{K}^{(h)}}\boldsymbol{M}^{(h)}\approx\boldsymbol{S}^{(h)}=\left(\boldsymbol{M}^{(h)}\right)^{T}\boldsymbol{\Sigma}^{(h)}_{\boldsymbol{K}^{(h)}}\boldsymbol{M}^{(h)}\quad\text{and}\quad$
(4.12) $\displaystyle\widetilde{\boldsymbol{S}}^{(h;M)}$
$\displaystyle:=\left(\boldsymbol{L}^{(h)}\right)^{T}\boldsymbol{\Sigma}^{(h;M)}_{\boldsymbol{K}^{(h)}}\boldsymbol{L}^{(h)}\approx\widetilde{\boldsymbol{S}}^{(h)}=\left(\boldsymbol{L}^{(h)}\right)^{T}\boldsymbol{\Sigma}^{(h)}_{\boldsymbol{K}^{(h)}}\boldsymbol{L}^{(h)},$
(4.13)

and we encounter the following eigenvalue problem: Seek
$(\widetilde{\boldsymbol{\Phi}}^{(h;M)}_{\ell},\lambda^{(h;M)}_{\ell})\in\mathbb{R}^{Q_{h}}\times\mathbb{R}$
for all $\ell=1,\cdots,Q_{h}$ such that

$\displaystyle\widetilde{\boldsymbol{S}}^{(h;M)}\widetilde{\boldsymbol{\Phi}}^{(h;M)}_{\ell}=\lambda^{(h;M)}_{\ell}\widetilde{\boldsymbol{\Phi}}^{(h;M)}_{\ell}.$
(4.14)

As $\widetilde{\boldsymbol{S}}^{(h;M)}$ is a symmetric matrix, we have that
$\left(\widetilde{\boldsymbol{\Phi}}^{(h;M)}_{\ell},\widetilde{\boldsymbol{\Phi}}^{(h;M)}_{\ell^{\prime}}\right)_{\ell_{2}}=\delta_{\ell,\ell^{\prime}}$
after normalization. Moreover, we encounter the generalized eigenvalue problem

$\displaystyle\boldsymbol{S}^{(h;M)}\boldsymbol{\Phi}^{(h;M)}_{\ell}=\lambda^{(h;M)}_{\ell}\boldsymbol{M}^{(h)}\boldsymbol{\Phi}^{(h;M)}_{\ell}.$
(4.15)

Analogously, we have the equivalence

$\displaystyle(\widetilde{\boldsymbol{\Phi}}^{(h;M)}_{\ell},\lambda^{(h;M)}_{\ell})\in\mathbb{R}^{Q_{h}}\times\mathbb{R}\text{
is an eigenpair of \eqref{genalizedeigenvalue3-estimator}}\Longleftrightarrow$
$\displaystyle\left(\boldsymbol{\Phi}^{(h;M)}_{\ell}:=\left(\boldsymbol{L}^{(h)}\right)^{-T}\widetilde{\boldsymbol{\Phi}}^{(h;M)}_{\ell},\lambda^{(h;M)}_{\ell}\right)\in\mathbb{R}^{Q_{h}}\times\mathbb{R}\text{
is an eigenpair of \eqref{genalizedeigenvalue4-estimator}}.$

Thus we can derive approximations to $\widetilde{\boldsymbol{S}}^{(h)}$ and
$\boldsymbol{S}^{(h)}$ by

$\displaystyle\widetilde{\boldsymbol{S}}^{(h;M)}=\sum_{k=1}^{Q_{h}}\lambda^{(h;M)}_{k}\widetilde{\boldsymbol{\Phi}}^{(h;M)}_{k}\otimes\widetilde{\boldsymbol{\Phi}}^{(h;M)}_{k}\quad\text{and}\quad\boldsymbol{S}^{(h;M)}=\sum_{k=1}^{Q_{h}}\lambda^{(h;M)}_{k}\boldsymbol{\Phi}^{(h;M)}_{k}\otimes\boldsymbol{\Phi}^{(h;M)}_{k}\in\mathbb{R}^{Q_{h}\times
Q_{h}}.$ (4.16)

Then the eigensystem (4.16) gives rise to a computable Mercer-like expansion

$\displaystyle
R^{(h;M)}(\boldsymbol{x},\boldsymbol{x}^{\prime}):=\sum_{k=1}^{Q_{h}}\lambda^{(h;M)}_{k}\phi^{(h;M)}_{k}(\boldsymbol{x})\,\phi^{(h;M)}_{k}(\boldsymbol{x}^{\prime}).$
(4.17)

Here, we use analogously

$\displaystyle\phi_{k}^{(h;M)}:=\boldsymbol{\Phi}^{(h;M)}_{k}\cdot\boldsymbol{\theta}^{(h)}\text{
for all }k=1,\cdots,Q_{h}.$ (4.18)

We observe, using
$\boldsymbol{M}^{(h)}=\boldsymbol{L}^{(h)}\left(\boldsymbol{L}^{(h)}\right)^{T}\in\mathbb{R}^{Q_{h}\times
Q_{h}},$ that

$\displaystyle\left(\phi_{k}^{(h;M)},\phi_{k^{\prime}}^{(h;M)}\right)_{L_{2}(\mathcal{D})}$
$\displaystyle=\left(\boldsymbol{\Phi}^{(h;M)}_{k}\right)^{T}\boldsymbol{M}^{(h)}\boldsymbol{\Phi}^{(h;M)}_{k^{\prime}}=\left(\left(\boldsymbol{L}^{(h)}\right)^{T}\boldsymbol{\Phi}^{(h;M)}_{k},\left(\boldsymbol{L}^{(h)}\right)^{T}\boldsymbol{\Phi}^{(h;M)}_{k^{\prime}}\right)_{\ell_{2}}$
$\displaystyle=\left(\widetilde{\boldsymbol{\Phi}}^{(h;M)}_{k},\widetilde{\boldsymbol{\Phi}}^{(h;M)}_{k^{\prime}}\right)_{\ell_{2}}=\delta_{k,k^{\prime}}.$
(4.19)

To bound the sampling error for the eigenvalues, we can apply Weyl’s
inequality (see for instance [12, (Eq. 3.1)]) to any of the two eigenvalue
problems (4.14) or (4.15), respectively. We will focus on (4.14) as (4.15) is
a generalized eigenvalue problem, which makes it harder to derive bounds on
the eigenvectors. We obtain

$\displaystyle\left|\lambda^{(h)}_{\ell}-\lambda^{(h;M)}_{\ell}\right|\leq\left\|\widetilde{\boldsymbol{S}}^{(h)}-\widetilde{\boldsymbol{S}}^{(h;M)}\right\|_{2\to
2},\quad 1\leq\ell\leq\min\\{Q_{h},L\\}.$ (4.20)

Now we can take Proposition 4.1 into account to get that, under Assumption
1.1, there are constants $C_{1}$ and $h_{0}$ such that

$\displaystyle\left|\lambda_{\ell}-\lambda^{(h;M)}_{\ell}\right|$
$\displaystyle\leq\left|\lambda_{\ell}-\lambda^{(h)}_{\ell}\right|+\left|\lambda^{(h)}_{\ell}-\lambda^{(h;M)}_{\ell}\right|$
$\displaystyle\leq
C_{1}\lambda_{\ell}^{-1}h^{2s}+\left\|\widetilde{\boldsymbol{S}}^{(h)}-\widetilde{\boldsymbol{S}}^{(h;M)}\right\|_{2\to
2}\quad\text{ for all }\quad 0<h\leq h_{0}\text{ and
}1\leq\ell\leq\min\\{Q_{h},L\\}.$ (4.21)

For perturbation bounds on the eigenvectors, see [18], and for a recent
versions of the Davis-Kahan theorem including randomness, see [13]. We use the
variant presented in [20, Corollary 1]. To this end, recall the assumption
that the eigenvalues are ordered non-increasingly, i.e.,
$\lambda^{(h)}_{1}\geq\dots\geq\lambda^{(h)}_{Q_{h}}$ and
$\lambda^{(h;M)}_{1}\geq\dots\geq\lambda^{(h;M)}_{Q_{h}}$. We define the
discrete spectral gap

$\displaystyle\delta^{(h;M)}_{\ell}:=\min\left\\{\left|\lambda^{(h;M)}_{\ell-1}-\lambda^{(h)}_{\ell}\right|,\left|\lambda^{(h)}_{\ell}-\lambda^{(h;M)}_{\ell+1}\right|\right\\},\quad\text{with
}\lambda^{(h;M)}_{0}=\infty\text{ and }\lambda^{(h;M)}_{Q_{h}+1}=-\infty.$
(4.22)

The quantity $\delta_{\ell}^{(h;M)}$ is problematic as it contains
$\lambda^{(h)}_{\ell},\lambda^{(h;M)}_{\ell\pm 1}$. We therefore aim to
replace it by a quantity which only depends on the continuous problem. To this
end, we adopt the strategy presented in [20]: Let $\delta_{\ell}$ be the
continuous spectral gap defined by

$\displaystyle\delta_{\ell}:=\min\left\\{\lambda_{\ell-1}-\lambda_{\ell},\lambda_{\ell}-\lambda_{\ell+1}\right\\}\text{
with }\lambda_{0}=\infty.$ (4.23)

######  Assumption 4.1 (Spectral gap).

Assume that there is a sufficiently small positive parameter $h_{1}\leq h_{0}$
and assume that $\widetilde{\boldsymbol{S}}^{(h;M)}$ is a good approximation
to $\widetilde{\boldsymbol{S}}^{(h)}$, i.e., for $h\leq h_{1}$, there holds

$\displaystyle\delta_{\ell}\geq
4C_{1}h^{2s}\lambda_{\ell+1}^{-1}+4\left\|\widetilde{\boldsymbol{S}}^{(h)}-\widetilde{\boldsymbol{S}}^{(h;M)}\right\|_{2\to
2}.$ (4.24)

######  Theorem 4.1.

Let (4.24) be valid. It then holds

$\displaystyle\delta^{(h;M)}_{\ell}\geq\frac{1}{4}\delta_{\ell}.$ (4.25)

Furthermore it holds

$\displaystyle\left\|\widetilde{\boldsymbol{\Phi}}^{(h)}_{\ell}-\widetilde{\boldsymbol{\Phi}}^{(h;M)}_{\ell}\right\|_{2}\leq
C\frac{\left\|\widetilde{\boldsymbol{S}}^{(h)}-\widetilde{\boldsymbol{S}}^{(h;M)}\right\|_{2\to
2}}{\delta^{(h;M)}_{\ell}}\leq
4C\frac{\left\|\widetilde{\boldsymbol{S}}^{(h)}-\widetilde{\boldsymbol{S}}^{(h;M)}\right\|_{2\to
2}}{\delta_{\ell}},$ (4.26)

where we fixed a sign for $\widetilde{\boldsymbol{\Phi}}^{(h;M)}_{\ell}$ such
that
$\widetilde{\boldsymbol{\Phi}}^{(h)}_{\ell}\cdot\widetilde{\boldsymbol{\Phi}}^{(h;M)}_{\ell}\geq
0$. Moreover for
$\phi_{\ell}^{(h;M)}:=\left(\boldsymbol{L}^{(h)}\right)^{T}\widetilde{\boldsymbol{\Phi}}^{(h;M)}_{k}\cdot\boldsymbol{\theta}^{(h)}$
the bound

$\displaystyle\left\|\phi_{\ell}^{(h)}-\phi_{\ell}^{(h;M)}\right\|_{L_{2}(\mathcal{D})}\leq
C\frac{\left\|\widetilde{\boldsymbol{S}}^{(h)}-\widetilde{\boldsymbol{S}}^{(h;M)}\right\|_{2\to
2}}{\delta^{(h;M)}_{\ell}}\leq
4C\frac{\left\|\widetilde{\boldsymbol{S}}^{(h)}-\widetilde{\boldsymbol{S}}^{(h;M)}\right\|_{2\to
2}}{\delta_{\ell}}$ (4.27)

holds.

###### Proof.

The first assertion is a generalization of an argument presented in [20] to
the case of several discretization parameters. Recall the discrete spectral
gap $\delta^{(h;M)}_{\ell}$ in (4.22). We show (4.25) by the following
strategy: First, we observe

$\displaystyle\delta_{\ell}$
$\displaystyle\leq\left|\lambda_{\ell-1}-\lambda_{\ell}\right|=\left|\lambda_{\ell-1}-\lambda^{(h)}_{\ell-1}+\lambda^{(h)}_{\ell-1}-\lambda^{(h)}_{\ell}+\lambda^{(h)}_{\ell}-\lambda_{\ell}\right|$
$\displaystyle\leq\left|\lambda_{\ell-1}-\lambda^{(h)}_{\ell-1}\right|+\left|\lambda^{(h)}_{\ell-1}-\lambda^{(h)}_{\ell}\right|+\left|\lambda^{(h)}_{\ell}-\lambda_{\ell}\right|.$

Similarly, we obtain

$\displaystyle\delta_{\ell}$
$\displaystyle\leq\left|\lambda_{\ell}-\lambda_{\ell+1}\right|=\left|\lambda_{\ell}-\lambda^{(h)}_{\ell}+\lambda^{(h)}_{\ell}-\lambda^{(h)}_{\ell+1}+\lambda^{(h)}_{\ell+1}-\lambda_{\ell+1}\right|$
$\displaystyle\leq\left|\lambda_{\ell}-\lambda^{(h)}_{\ell}\right|+\left|\lambda^{(h)}_{\ell}-\lambda^{(h)}_{\ell+1}\right|+\left|\lambda^{(h)}_{\ell+1}-\lambda_{\ell+1}\right|.$

Using Proposition 4.1 and the fact that the eigenvalues are sorted, i.e.,
$\lambda_{\ell-1}\geq\lambda_{\ell}\geq\lambda_{\ell+1}$, we obtain for all
$0<h\leq h_{0}$ the bound

$\displaystyle\max\left\\{\left|\lambda_{\ell-1}-\lambda^{(h)}_{\ell-1}\right|,\left|\lambda^{(h)}_{\ell}-\lambda_{\ell}\right|,\left|\lambda_{\ell+1}-\lambda^{(h)}_{\ell+1}\right|\right\\}\leq
C_{1}\lambda_{\ell+1}^{-1}h^{2s}\leq\frac{1}{4}\delta_{\ell},$

where we have used the spectral gap assumption (4.24) in the last step. Hence,
we infer

$\displaystyle\delta_{\ell}\leq\min\left\\{\left|\lambda^{(h)}_{\ell-1}-\lambda^{(h)}_{\ell}\right|,\left|\lambda^{(h)}_{\ell}-\lambda^{(h)}_{\ell+1}\right|\right\\}+\frac{1}{2}\delta_{\ell}=\delta^{(h)}_{\ell}+\frac{1}{2}\delta_{\ell}.$

This inequality implies

$\displaystyle\frac{1}{2}\delta_{\ell}\leq\delta^{(h)}_{\ell}=\min\left\\{\left|\lambda^{(h)}_{\ell-1}-\lambda^{(h)}_{\ell}\right|,\left|\lambda^{(h)}_{\ell}-\lambda^{(h)}_{\ell+1}\right|\right\\}.$
(4.28)

Moreover, we have

$\displaystyle\delta^{(h)}_{\ell}\leq\left|\lambda^{(h)}_{\ell-1}-\lambda^{(h)}_{\ell}\right|=\left|\lambda^{(h)}_{\ell-1}-\lambda^{(h;M)}_{\ell-1}+\lambda^{(h;M)}_{\ell-1}-\lambda^{(h)}_{\ell}\right|\leq\left|\lambda^{(h)}_{\ell-1}-\lambda^{(h;M)}_{\ell-1}\right|+\left|\lambda^{(h;M)}_{\ell-1}-\lambda^{(h)}_{\ell}\right|$

and also

$\displaystyle\delta^{(h)}_{\ell}\leq\left|\lambda^{(h)}_{\ell}-\lambda^{(h)}_{\ell+1}\right|=\left|\lambda^{(h)}_{\ell}-\lambda^{(h;M)}_{\ell}+\lambda^{(h;M)}_{\ell}-\lambda^{(h)}_{\ell+1}\right|\leq\left|\lambda^{(h)}_{\ell}-\lambda^{(h;M)}_{\ell}\right|+\left|\lambda^{(h;M)}_{\ell}-\lambda^{(h)}_{\ell+1}\right|.$

Using Weyl’s Theorem, i.e., (4.20), and the spectral gap assumption (4.24), we
obtain

$\displaystyle\max\left\\{\left|\lambda^{(h)}_{\ell-1}-\lambda^{(h;M)}_{\ell-1}\right|,\left|\lambda^{(h)}_{\ell}-\lambda^{(h;M)}_{\ell}\right|\right\\}\leq\left\|\widetilde{\boldsymbol{S}}^{(h)}-\widetilde{\boldsymbol{S}}^{(h;M)}\right\|_{2\to
2}\leq\frac{1}{4}\delta_{\ell}\leq\frac{1}{2}\delta^{(h)}_{\ell},$

where we have used (4.28) in the last step. Hence, we obtain

$\displaystyle\delta^{(h;M)}_{\ell}=\min\left\\{\left|\lambda^{(h)}_{\ell}-\lambda^{(h;M)}_{\ell-1}\right|,\left|\lambda^{(h)}_{\ell+1}-\lambda^{(h;M)}_{\ell}\right|\right\\}\geq\frac{1}{2}\delta^{(h)}_{\ell}\geq\frac{1}{4}\delta_{\ell}$

where have used again (4.28) in the last step.

Moreover, we observe

$\displaystyle\left\|\phi_{\ell}^{(h)}-\phi_{\ell}^{(h;M)}\right\|^{2}_{L_{2}(\mathcal{D})}=\left(\Phi^{(h)}_{\ell}-\Phi^{(h;M)}_{\ell}\right)^{T}\boldsymbol{M}^{(h)}\left(\Phi^{(h)}_{\ell}-\Phi^{(h;M)}_{\ell}\right)=\left\|\widetilde{\Phi}^{(h)}_{\ell}-\widetilde{\Phi}^{(h;M)}_{\ell}\right\|^{2}_{\ell_{2}}.$

The two other assertions follow directly with the Davis-Kahan theorem as in
[20, Corollary1]. ∎

Finally, we are left with bounding
$\left\|\widetilde{\boldsymbol{S}}^{(h)}-\widetilde{\boldsymbol{S}}^{(h;M)}\right\|_{2\to
2}$. To this end, we first observe

$\displaystyle\widetilde{\boldsymbol{S}}^{(h)}-\widetilde{\boldsymbol{S}}^{(h;M)}$
$\displaystyle:=\left(\boldsymbol{L}^{(h)}\right)^{T}\left(\boldsymbol{\Sigma}_{\boldsymbol{K}^{(h)}}-\boldsymbol{\Sigma}^{(h;M)}_{\boldsymbol{K}^{(h)}}\right)\boldsymbol{L}^{(h)}\in\mathbb{R}^{Q_{h}\times
Q_{h}}.$

Moreover, we have

$\displaystyle\left\|\widetilde{\boldsymbol{S}}^{(h)}-\widetilde{\boldsymbol{S}}^{(h;M)}\right\|_{2\to
2}$
$\displaystyle=\sup_{\boldsymbol{x}\in\mathbb{R}^{Q_{h}}\setminus\\{\boldsymbol{0}\\}}\frac{\left(\left(\widetilde{\boldsymbol{S}}^{(h)}-\widetilde{\boldsymbol{S}}^{(h;M)}\right)\boldsymbol{x},\boldsymbol{x}\right)_{\ell_{2}}}{\|\boldsymbol{x}\|^{2}}$
$\displaystyle=\sup_{\boldsymbol{y}\in\mathbb{R}^{Q_{h}}\setminus\\{\boldsymbol{0}\\}}\frac{\left(\left(\boldsymbol{\Sigma}_{\boldsymbol{K}^{(h)}}-\boldsymbol{\Sigma}^{(h;M)}_{\boldsymbol{K}^{(h)}}\right)\boldsymbol{y},\boldsymbol{y}\right)_{\ell_{2}}}{\|\left(\boldsymbol{L}^{(h)}\right)^{-1}\boldsymbol{y}\|_{\ell_{2}}^{2}},$
(4.29)

where the choice $\boldsymbol{y}=\boldsymbol{L}^{(h)}\boldsymbol{x}$ is
admissible since $\boldsymbol{L}^{(h)}$ is regular. Note here also that we
have

$\displaystyle\lambda_{\min}\left(\left(\boldsymbol{M}^{(h)}\right)^{-1}\right)\|\boldsymbol{y}\|^{2}\leq\left(\left(\boldsymbol{M}^{(h)}\right)^{-1}\boldsymbol{y},\boldsymbol{y}\right)\leq\lambda_{\max}\left(\left(\boldsymbol{M}^{(h)}\right)^{-1}\right)\|\boldsymbol{y}\|^{2}$

and we observe
$\lambda_{\min}\left(\left(\boldsymbol{M}^{(h)}\right)^{-1}\right)=\lambda^{-1}_{\max}\left(\boldsymbol{M}^{(h)}\right)$
and
$\lambda_{\max}\left(\left(\boldsymbol{M}^{(h)}\right)^{-1}\right)=\lambda^{-1}_{\min}\left(\boldsymbol{M}^{(h)}\right)$.
Consequently, we obtain the bound

$\lambda_{\min}\left(\boldsymbol{M}^{(h)}\right)\left\|\boldsymbol{\Sigma}_{\boldsymbol{K}^{(h)}}-\boldsymbol{\Sigma}^{(h;M)}_{\boldsymbol{K}^{(h)}}\right\|_{2\to
2}\leq\left\|\widetilde{\boldsymbol{S}}^{(h)}-\widetilde{\boldsymbol{S}}^{(h;M)}\right\|_{2\to
2}\leq\lambda_{\max}\left(\boldsymbol{M}^{(h)}\right)\left\|\boldsymbol{\Sigma}_{\boldsymbol{K}^{(h)}}-\boldsymbol{\Sigma}^{(h;M)}_{\boldsymbol{K}^{(h)}}\right\|_{2\to
2}.$ (4.30)

The matrix norm
$\left\|\boldsymbol{\Sigma}_{\boldsymbol{K}^{(h)}}-\boldsymbol{\Sigma}^{(h;M)}_{\boldsymbol{K}^{(h)}}\right\|_{2\to
2}$ will be estimated in probabilistic terms in Section 4.4.

###  4.3 Sub-Gaussian tails for observed random vectors

In order to invoke recent results for bounding the sampling error, we first
need to consider the sub-Gaussian property of the random variables (1.14),
i.e., of the random vector $\boldsymbol{K}^{(h)}$. This random vector consists
of the coefficients of the discrete random field
$\kappa^{(h)}:\Omega\times\mathcal{D}\to\mathbb{R}$, i.e.,

$\displaystyle\kappa^{(h)}(\boldsymbol{\omega},\boldsymbol{x})=\sum_{k=1}^{Q_{h}}K^{(h)}_{k}(\boldsymbol{\omega})\theta_{k}^{(h)}(\boldsymbol{x})=\boldsymbol{K}^{(h)}(\boldsymbol{\omega})\cdot\boldsymbol{\theta}(\boldsymbol{x}).$

We mainly focus on two different basis functions. First, we choose
$\\{\theta^{(h)}_{k}\\}_{1\leq k\leq Q_{h}}$ to be a nodal basis and, next, we
choose $\\{\theta^{(h)}_{k}\\}_{1\leq k\leq Q_{h}}$ to be an
$L^{2}$-orthonormal basis.

####  4.3.1 Nodal basis

To start, let us consider the setting of so-called standard information. To
this end, recall the point set
$X_{Q_{h}}=\left\\{\boldsymbol{x}_{1},\dots,\boldsymbol{x}_{Q_{h}}\right\\}\subset\mathcal{D}$,
which determines the nodal basis functions and thus the degrees of freedom for
the finite element space $\mathcal{V}_{h}$. In this case, it holds

$\displaystyle\kappa^{(h)}(\boldsymbol{\omega},\boldsymbol{x})=\boldsymbol{K}^{(h)}(\boldsymbol{\omega})\cdot\boldsymbol{\theta}(\boldsymbol{x})=\sum_{k=1}^{Q_{h}}\kappa^{(h)}_{k}(\boldsymbol{\omega},\boldsymbol{x}_{k})\theta_{k}^{(h)}(\boldsymbol{x}).$

Thus, we have

$\displaystyle\boldsymbol{K}^{(h)}:\Omega\to\mathbb{R}^{Q_{h}},\quad\boldsymbol{\omega}\mapsto\boldsymbol{K}^{(h)}(\boldsymbol{\omega})=\left(\kappa^{(h)}(\boldsymbol{\omega},\boldsymbol{x}_{1}),\dots,\kappa^{(h)}(\boldsymbol{\omega},\boldsymbol{x}_{Q_{h}})\right)^{T}.$
(4.31)

Since $\kappa^{(h)}$ is assumed to be a centered Gaussian random field, it
follows by definition that the random vector $\boldsymbol{K}^{(h)}$ is also
distributed according to a multi-variate Gaussian law, i.e.,

$\displaystyle\boldsymbol{K}^{(h)}\sim\mathcal{N}\left(\boldsymbol{0},\boldsymbol{\Sigma}_{\boldsymbol{K}^{(h)}}\right)$

with

$\displaystyle\boldsymbol{\Sigma}_{\boldsymbol{K}^{(h)}}:=\mathbb{E}\left[\left(\boldsymbol{K}^{(h)}-\mathbb{E}\left[\boldsymbol{K}^{(h)}\right]\right)\otimes\left(\boldsymbol{K}^{(h)}-\mathbb{E}\left[\boldsymbol{K}^{(h)}\right]\right)\right]=\mathbb{E}\left[\boldsymbol{K}^{(h)}\otimes\boldsymbol{K}^{(h)}\right].$

We can invoke (1.17) and (1.18), i.e.,
$\left(\boldsymbol{\Sigma}_{\boldsymbol{K}^{(h)}}\right)_{k,k^{\prime}}\leq
4(C_{\infty}^{(h)})^{2}$ . Finally, we obtain by Chernoff’s inequality [8] the
bound

$\displaystyle\mathbb{P}\left\\{\boldsymbol{v}\cdot\left(\boldsymbol{K}^{(h)}-\mathbb{E}\left[\boldsymbol{K}^{(h)}\right]\right)>t\right\\}\leq\exp\left(-\frac{t^{2}}{8(C_{\infty}^{(h)})^{2}}\right)$

for any $\boldsymbol{v}$ with $\|\boldsymbol{v}\|_{2}=1$. Note that this
implies $4(C_{\infty}^{(h)})^{2}=\rho^{-1}$ in the sense of [6, Definition 1].

####  4.3.2 $L^{2}$-orthonormal basis

Now let us consider the setting of so-called linear information. Here, we can
not directly rely on the definition of a centered Gaussian random field to
deduce that the coefficients are distributed with a multivariate normal law.
In this case, the random vector (1.14) gets encoded as

$\displaystyle\boldsymbol{K}^{(h)}:\Omega\to\mathbb{R}^{Q_{h}},\quad\boldsymbol{\omega}$
$\displaystyle\mapsto\left(\int_{\mathcal{D}}\kappa^{(h)}(\boldsymbol{\omega},\boldsymbol{x})\theta^{(h)}_{1}(\boldsymbol{x})\,\mathrm{d}\mathcal{L}^{(d)}(\boldsymbol{x}),\dots,\int_{\mathcal{D}}\kappa^{(h)}(\boldsymbol{\omega},\boldsymbol{x})\theta^{(h)}_{Q_{h}}(\boldsymbol{x})\,\mathrm{d}\mathcal{L}^{(d)}(\boldsymbol{x})\right)^{T}$
$\displaystyle=\int_{\mathcal{D}}\kappa^{(h)}(\boldsymbol{\omega},\boldsymbol{x})\boldsymbol{\theta}^{(h)}(\boldsymbol{x})\,\mathrm{d}\mathcal{L}^{(d)}(\boldsymbol{x}),$

where $\\{\theta_{j}^{(h)},j=1,\ldots,Q_{h}\\}$ is an $L_{2}$-orthonormal
basis. The individual random variables are given as

$\displaystyle
K^{(h)}_{k}:\Omega\to\mathbb{R}^{Q_{h}},\quad\boldsymbol{\omega}\mapsto\int_{\mathcal{D}}\kappa^{(h)}(\boldsymbol{\omega},\boldsymbol{x})\theta^{(h)}_{k}(\boldsymbol{x})\,\mathrm{d}\mathcal{L}^{(d)}(\boldsymbol{x}).$
(4.32)

Hence, for the expected value, we observe

$\displaystyle\mathbb{E}\left[K^{(h)}_{k}\right]$
$\displaystyle=\int_{\Omega}\int_{\mathcal{D}}\kappa^{(h)}(\boldsymbol{\omega},\boldsymbol{x})\theta^{(h)}_{k}(\boldsymbol{x})\,\mathrm{d}\mathcal{L}^{(d)}(\boldsymbol{x})\,\mathrm{d}\mathbb{P}(\boldsymbol{\omega})=\int_{\mathcal{D}}\int_{\Omega}\kappa^{(h)}(\boldsymbol{\omega},\boldsymbol{x})\,\mathrm{d}\mathbb{P}(\boldsymbol{\omega})\theta^{(h)}_{k}(\boldsymbol{x})\,\mathrm{d}\mathcal{L}^{(d)}(\boldsymbol{x})$
$\displaystyle=\int_{\mathcal{D}}\mathbb{E}\left[\kappa^{(h)}(\cdot,\boldsymbol{x})\right]\theta^{(h)}_{k}(\boldsymbol{x})\,\mathrm{d}\mathcal{L}^{(d)}(\boldsymbol{x})=0,$

which implies

$\displaystyle\mathbb{E}\left[\boldsymbol{K}^{(h)}\right]=\int_{\mathcal{D}}\mathbb{E}\left[\kappa^{(h)}(\cdot,\boldsymbol{x})\right]\boldsymbol{\theta}^{(h)}(\boldsymbol{x})\,\mathrm{d}\mathcal{L}^{(d)}(\boldsymbol{x})=\boldsymbol{0}.$

For the variance, we observe

$\displaystyle\mathbb{E}\left[\left(\boldsymbol{K}^{(h)}-\mathbb{E}\left[\boldsymbol{K}^{(h)}\right]\right)\otimes\left(\boldsymbol{K}^{(h)}-\mathbb{E}\left[\boldsymbol{K}^{(h)}\right]\right)\right]=\mathbb{E}\left[\boldsymbol{K}^{(h)}\otimes\boldsymbol{K}^{(h)}\right]=:\boldsymbol{\Sigma}_{\boldsymbol{K}^{(h)}},$

where

$\left(\boldsymbol{\Sigma}_{\boldsymbol{K}^{(h)}}\right)_{k,k^{\prime}}=\int_{\mathcal{D}}\int_{\mathcal{D}}\text{Cov}_{\kappa^{(h)}}(\boldsymbol{x},\boldsymbol{x}^{\prime})\theta_{k}^{h}(\boldsymbol{x})\theta_{k^{\prime}}^{h}(\boldsymbol{x}^{\prime})\mathrm{d}\boldsymbol{x}\mathrm{d}\boldsymbol{x}^{\prime}$
(4.33)

and

$\left(\boldsymbol{\Sigma}_{\boldsymbol{K}^{(h)}}\right)_{k,k^{\prime}}\leq\|\text{Cov}_{\kappa^{(h)}}\|_{L^{2}(\mathcal{D}\times\mathcal{D})}\leq
4(C_{\infty}^{(h)})^{2}\cdot|\mathcal{D}|.$ (4.34)

For all $\boldsymbol{v}\in\mathbb{R}^{Q_{h}}$ with $\|\boldsymbol{v}\|_{2}=1$,
we have

$\displaystyle\boldsymbol{v}\cdot\left(\boldsymbol{K}^{(h)}-\mathbb{E}\left[\boldsymbol{K}^{(h)}\right]\right)$
$\displaystyle=\sum_{k=1}^{Q_{h}}\left(\kappa^{(h)}(\omega_{j},\cdot),\theta_{j}^{(h)}\right)_{L^{2}(\mathcal{D})}v_{j}\sim\mathcal{N}(0,\boldsymbol{v}\cdot\boldsymbol{\Sigma}_{\boldsymbol{K}^{(h)}}\boldsymbol{v}).$

Furthermore, given $\boldsymbol{x}\in\mathcal{D}$, it holds

$\displaystyle\kappa^{(h)}(\cdot,\boldsymbol{x})\sim\mathcal{N}\left(0,\text{Cov}_{\kappa^{(h)}}(\boldsymbol{x},\boldsymbol{x})\right).$
(4.35)

Therefore we can invoke Chernoff’s inequality to obtain for all
$\boldsymbol{x}\in\mathcal{D}$ and $t>0$

$\displaystyle\mathbb{P}\left\\{\
\left|\kappa^{(h)}(\omega,\boldsymbol{x})\right|>t\right\\}$
$\displaystyle\leq\exp\left(-\frac{t^{2}}{2\text{Cov}_{\kappa^{(h)}}(\boldsymbol{x},\boldsymbol{x})}\right)\leq\exp\left(-\frac{t^{2}}{8(C_{\infty}^{(h)})^{2}|\mathcal{D}|}\right),$

which is uniform for $\boldsymbol{x}\in\mathcal{D}$, see (4.34). Consequently,
we have, for all $\boldsymbol{v}\in\mathbb{R}^{Q_{h}}$ with
$\|\boldsymbol{v}\|_{2}=1$, the estimate

$\displaystyle\mathbb{P}\left\\{\left|\boldsymbol{v}\cdot\left(\boldsymbol{K}^{(h)}-\mathbb{E}\left[\boldsymbol{K}^{(h)}\right]\right)\right|>t\right\\}\leq\exp\left(-\frac{t^{2}}{{8(C_{\infty}^{(h)}})^{2}|\mathcal{D}|}\right).$

Note that this implies $4(C_{\infty}^{(h)})^{2}|\mathcal{D}|=\rho^{-1}$ in the
sense of [6, Definition 1]. Hence, the random vector $\boldsymbol{K}^{(h)}$
obeys a sub-Gaussian bound in the sense of [7, Eq. (7)]. We will make use of
this estimate together with the results of [7] to bound the variance of the
sampled covariance matrix later on.

###  4.4 Covariance estimation from samples using tapering and decay
assumptions

In this subsection, we focus on the sampling error and specifically on the
induced variance of our estimator for the covariance matrix. We assume in this
subsection the parameter $h$ to be fixed. Thus the dimension of the covariance
matrix is given as $Q_{h}$. We will provide a constructive approach to obtain
an estimation to the covariance matrix
$\boldsymbol{\Sigma}_{\boldsymbol{K}^{(h)}}$. To this end, with given samples
$\boldsymbol{K}^{(m;h)}$ for $m=1,\cdots,M$, let the sample mean and the
maximum likelihood estimator for the covariance matrix
$\boldsymbol{\Sigma}_{\boldsymbol{K}^{(h)}}$ be [7, Eq. (2)]

$\displaystyle\bar{\boldsymbol{K}}^{(h;M)}$
$\displaystyle:=\frac{1}{M}\sum_{m=1}^{M}\boldsymbol{K}^{(m;h)}\in\mathbb{R}^{Q_{h}}\quad\text{and}$
(4.36) $\displaystyle\bar{\boldsymbol{\Sigma}}^{(h;M)}_{\boldsymbol{K}^{(h)}}$
$\displaystyle:=\frac{1}{M}\sum_{m=1}^{M}\left(\boldsymbol{K}^{(m;h)}-\bar{\boldsymbol{K}}^{(h)}\right)\otimes\left(\boldsymbol{K}^{(m;h)}-\bar{\boldsymbol{K}}^{(h)}\right)\in\mathbb{R}^{Q_{h}\times
Q_{h}}.$ (4.37)

We know from [6, Lemma 1] that the direct covariance estimator (4.37) suffers
from the curse of dimension with respect to high spatial resolution $Q_{h}$.
Consequently, we need a better estimator when $Q_{h}$ is large. To this end,
we follow [6] and assume an off diagonal decay of the covariance matrix
$\boldsymbol{\Sigma}_{\boldsymbol{K}^{(h)}}\in\mathbb{R}^{Q_{h}\times Q_{h}}$,
i.e., we assume
$\boldsymbol{\Sigma}_{\boldsymbol{K}^{(h)}}\in\mathcal{F}_{\alpha}=\mathcal{F}_{\alpha}(C_{\mathcal{F};1},C_{\mathcal{F};2})$,
where

$\displaystyle\mathcal{F}_{\alpha}:=\left\\{\boldsymbol{\Sigma}\in\mathbb{R}^{Q_{h}\times
Q_{h}}\ :\
\max_{k=1,\dots,Q_{h}}\sum_{\genfrac{}{}{0.0pt}{}{k^{\prime}=1}{\left|k^{\prime}-k\right|>c}}^{Q_{h}}\left|\boldsymbol{\Sigma}_{k,k^{\prime}}\right|\leq
C_{\mathcal{F};1}c^{-\alpha}\text{ for all }1\leq c\leq Q_{h}\text{ and
}\lambda_{\max}(\boldsymbol{\Sigma})\leq C_{\mathcal{F};2}\right\\}.$ (4.38)

Here, $\lambda_{\max}(\boldsymbol{\Sigma})$ is the maximum eigenvalue of
$\boldsymbol{\Sigma}$, and $\alpha$ modulates the speed of decay while
$C_{\mathcal{F};1}$ and $C_{\mathcal{F};2}$ are positive parameters. Now, let
weights for any even integer $1\leq\tau\leq Q_{h}$ (see [7, (Eq. 5)]) be given
as

$\displaystyle
w_{k,k^{\prime}}:=w_{k,k^{\prime}}(\tau):=\begin{cases}1,&\left|k-k^{\prime}\right|\leq\frac{\tau}{2},\\\
2\left(1-\frac{\left|k-k^{\prime}\right|}{\tau}\right),&\frac{\tau}{2}<\left|k-k^{\prime}\right|<\tau,\\\
0,&\left|k-k^{\prime}\right|\geq\tau.\end{cases}$ (4.39)

Then the associated tapering estimator (c.f. [7, (Eq. 4)]) is defined element-
wise by

$\displaystyle\left(\boldsymbol{\Sigma}^{(h;M;\tau)}_{\boldsymbol{K}^{(h)}}\right)_{{}_{k,k^{\prime}}}:=w_{k,k^{\prime}}(\tau)\left(\bar{\boldsymbol{\Sigma}}^{h;M}_{\boldsymbol{K}^{(h)}}\right)_{k,k^{\prime}},\quad
1\leq k,k^{\prime}\leq Q_{h}.$ (4.40)

Note that the estimator
$\boldsymbol{\Sigma}^{(h;M;\tau)}_{\boldsymbol{K}^{(h)}}$ is self-adjoint.
Moreover the bounds

$\displaystyle\mathbb{E}\left[\left\|\boldsymbol{\Sigma}^{(h;M;\tau)}_{\boldsymbol{K}^{(h)}}-\mathbb{E}\left[\boldsymbol{\Sigma}^{(h;M;\tau)}_{\boldsymbol{K}^{(h)}}\right]\right\|_{2\to
2}^{2}\right]$ $\displaystyle\lesssim\frac{\tau+\log(Q_{h})}{M},$ (4.41)
$\displaystyle\left\|\mathbb{E}\left[\boldsymbol{\Sigma}^{(h;M;\tau)}_{\boldsymbol{K}^{(h)}}\right]-\boldsymbol{\Sigma}_{\boldsymbol{K}^{(h)}}\right\|_{2\to
2}^{2}$ $\displaystyle\lesssim\tau^{-2\alpha}$ (4.42)

on the variance and the bias hold, see [7, (Eqs. 13 & 14)]. By the triangle
inequality we deduce the convergence of
$\boldsymbol{\Sigma}^{(h;M;\tau)}_{\boldsymbol{K}^{(h)}}\to\boldsymbol{\Sigma}_{\boldsymbol{K}^{(h)}}$
in expectation as

$\displaystyle\mathbb{E}\left[\left\|\boldsymbol{\Sigma}^{(h;M;\tau)}_{\boldsymbol{K}^{(h)}}-\boldsymbol{\Sigma}_{\boldsymbol{K}^{(h)}}\right\|_{2\to
2}^{2}\right]$ $\displaystyle\leq
2\left(\left\|\mathbb{E}\left[\boldsymbol{\Sigma}^{(h;M;\tau)}_{\boldsymbol{K}^{(h)}}\right]-\boldsymbol{\Sigma}_{\boldsymbol{K}^{(h)}}\right\|_{2\to
2}^{2}+\mathbb{E}\left[\left\|\boldsymbol{\Sigma}^{(h;M;\tau)}_{\boldsymbol{K}^{(h)}}-\mathbb{E}\left[\boldsymbol{\Sigma}^{(h;M;\tau)}_{\boldsymbol{K}^{(h)}}\right]\right\|_{2\to
2}^{2}\right]\right)$
$\displaystyle\lesssim\frac{\tau+\log(Q_{h})}{M}+\tau^{-2\alpha},$ (4.43)

which motivates the choice of the tapering parameter
$\tau=M^{\frac{1}{2\alpha+1}}$, c.f. [7, (Eq. 15)]. This gives the following
result, c.f. [7, Theorem 2 & (Eq. 31)].

######  Proposition 4.2.

For the tapering estimator (4.40) with $\tau=M^{\frac{1}{2\alpha+1}}$, i.e.,

$\displaystyle\left(\boldsymbol{\Sigma}^{(h;M)}_{\boldsymbol{K}^{(h)}}\right)_{{}_{k,k^{\prime}}}:=w_{k,k^{\prime}}(M^{\frac{1}{2\alpha+1}})\left(\bar{\boldsymbol{\Sigma}}^{(h;M)}_{\boldsymbol{K}^{(h)}}\right)_{k,k^{\prime}},\quad
1\leq k,k^{\prime}\leq Q_{h},$ (4.44)

there holds the error estimate

$\displaystyle\mathbb{E}\left[\left\|\boldsymbol{\Sigma}^{(h;M)}_{\boldsymbol{K}^{(h)}}-\boldsymbol{\Sigma}_{\boldsymbol{K}^{(h)}}\right\|_{2\to
2}^{2}\right]\lesssim M^{-\frac{2\alpha}{2\alpha+1}}+\frac{\log(Q_{h})}{M}$
(4.45)

if the condition $Q_{h}\geq M^{\frac{1}{2\alpha+1}}$ is satisfied. In the case
$Q_{h}<M^{\frac{1}{2\alpha+1}}$, we can directly use the estimator from (4.37)

$\displaystyle\bar{\boldsymbol{\Sigma}}^{(h;M)}_{\boldsymbol{K}^{(h)}}$
$\displaystyle:=\frac{1}{M}\sum_{m=1}^{M}\left(\boldsymbol{K}^{(m;h)}-\bar{\boldsymbol{K}}^{(h)}\right)\otimes\left(\boldsymbol{K}^{(m;h)}-\bar{\boldsymbol{K}}^{(h)}\right)\in\mathbb{R}^{Q_{h}\times
Q_{h}}$

and obtain the bound

$\displaystyle\mathbb{E}\left[\left\|\bar{\boldsymbol{\Sigma}}^{(h;M)}_{\boldsymbol{K}^{(h)}}-\boldsymbol{\Sigma}_{\boldsymbol{K}^{(h)}}\right\|_{2\to
2}^{2}\right]\lesssim\frac{Q_{h}}{M}.$ (4.46)

Note that the result is in fact _rate optimal_.

We introduce the notation

$\displaystyle\rho_{h}(M)$
$\displaystyle:=\begin{cases}M^{-\frac{2\alpha}{2\alpha+1}}+\frac{\log(Q_{h})}{M},&Q_{h}\geq
M^{\frac{1}{2\alpha+1}}\\\
\frac{Q_{h}}{M},&Q_{h}<M^{\frac{1}{2\alpha+1}}\end{cases}$ (4.47)
$\displaystyle\lesssim\tilde{\rho}_{h}(M):=\begin{cases}M^{-\frac{2\alpha}{2\alpha+1}}+\frac{d\log(h^{-1})}{M},&Q_{h}\geq
M^{\frac{1}{2\alpha+1}}\\\
\frac{h^{-d}}{M},&Q_{h}<M^{\frac{1}{2\alpha+1}}\end{cases},$ (4.48)

using $Q_{h}\sim s^{d}h^{-d}$ for the definition of $\tilde{\rho}_{h}(M)$.
Then, we obtain the following corollary.

######  Corollary 4.1.

For the tapering estimator (4.40) with $\tau=M^{\frac{1}{2\alpha+1}}$ the
bounds

$\displaystyle\mathbb{E}\left[\left\|\widetilde{\boldsymbol{S}}^{(h)}-\widetilde{\boldsymbol{S}}^{(h;M)}\right\|_{2\to
2}\right]$ $\displaystyle\lesssim\rho^{\frac{1}{2}}_{h}(M)\leavevmode\nobreak\
\lambda_{\max}\left(\boldsymbol{M}^{(h)}\right),$ (4.49)
$\displaystyle\mathbb{E}\left[\left\|\widetilde{\boldsymbol{S}}^{(h)}-\widetilde{\boldsymbol{S}}^{(h;M)}\right\|^{2}_{2\to
2}\right]$
$\displaystyle=\mathbb{E}\left[\left\|\left(\boldsymbol{\Sigma}^{(h;M)}_{\boldsymbol{K}^{(h)}}-\boldsymbol{\Sigma}_{\boldsymbol{K}^{(h)}}\right)\right\|^{2}_{2\to
2}\right]\lambda_{\max}^{2}\left(\boldsymbol{M}^{(h)}\right)\lesssim\rho_{h}(M)\leavevmode\nobreak\
\lambda_{\max}^{2}\left(\boldsymbol{M}^{(h)}\right)$ (4.50)

hold.

###### Proof.

First, we obtain by (4.30)

$\displaystyle\left\|\widetilde{\boldsymbol{S}}^{(h)}-\widetilde{\boldsymbol{S}}^{(h;M)}\right\|_{2\to
2}\leq\lambda_{\max}\left(\boldsymbol{M}^{(h)}\right)\left\|\boldsymbol{\Sigma}^{(h;M)}_{\boldsymbol{K}^{(h)}}-\boldsymbol{\Sigma}_{\boldsymbol{K}^{(h)}}\right\|_{2\to
2}.$ (4.51)

We can apply Jensen’s inequality for the convex function $\psi(x)=x^{2}$,

$\displaystyle\psi\left(\mathbb{E}\left[X\right]\right)\leq\mathbb{E}\left[\psi(X)\right].$

Plugging (4.51) into Jensen’s inequality yields

$\displaystyle\left(\mathbb{E}\left[\left\|\widetilde{\boldsymbol{S}}^{(h)}-\widetilde{\boldsymbol{S}}^{(h;M)}\right\|_{2\to
2}\right]\right)^{2}\leq\mathbb{E}\left[\left\|\left(\boldsymbol{\Sigma}^{(h;M)}_{\boldsymbol{K}^{(h)}}-\boldsymbol{\Sigma}_{\boldsymbol{K}^{(h)}}\right)\right\|^{2}_{2\to
2}\right]\lambda_{\max}^{2}\left(\boldsymbol{M}^{(h)}\right).$

Thus, we get from Proposition 4.2 and with the definition (4.47) the bound

$\displaystyle\mathbb{E}\left[\left\|\widetilde{\boldsymbol{S}}^{(h)}-\widetilde{\boldsymbol{S}}^{(h;M)}\right\|_{2\to
2}\right]\lesssim\rho^{\frac{1}{2}}_{h}(M)\leavevmode\nobreak\
\lambda_{\max}\left(\boldsymbol{M}^{(h)}\right).$

Moreover, we have

$\displaystyle\mathbb{E}\left[\left\|\widetilde{\boldsymbol{S}}^{(h)}-\widetilde{\boldsymbol{S}}^{(h;M)}\right\|^{2}_{2\to
2}\right]\leq\lambda_{\max}^{2}\left(\boldsymbol{M}^{(h)}\right)\mathbb{E}\left[\left\|\left(\boldsymbol{\Sigma}^{(h;M)}_{\boldsymbol{K}^{(h)}}-\boldsymbol{\Sigma}_{\boldsymbol{K}^{(h)}}\right)\right\|^{2}_{2\to
2}\right]\lesssim\rho_{h}(M)\leavevmode\nobreak\
\lambda_{\max}^{2}\left(\boldsymbol{M}^{(h)}\right).$

∎

Next, we need to verify that Assumption (4.24), i.e.,

$\displaystyle\delta_{\ell}\geq
4C_{1}h^{2s}\lambda_{\ell+1}^{-1}+4\left\|\widetilde{\boldsymbol{S}}^{(h)}-\widetilde{\boldsymbol{S}}^{(h;M)}\right\|_{2\to
2}\quad\text{for all }1\leq\ell\leq L,$

holds at least with high probability. To this end, we have the following
result:

######  Proposition 4.3.

Let the mesh size $h$ be sufficiently small and $M$ be sufficiently large such
that, for all $1\leq\ell\leq L$,

$\displaystyle\frac{\delta_{\ell}}{2}\geq
4C_{1}h^{2s}\lambda_{L+1}^{-1}\quad\text{and}\quad
CM^{-\frac{2\alpha}{2\alpha+1}}\leq\frac{\min_{1\leq\ell\leq
L}\delta_{\ell}}{16\leavevmode\nobreak\
\lambda_{\max}\left(\boldsymbol{M}^{(h)}\right)}$ (4.52)

holds. Then we have with probability

$\displaystyle
p_{0}:=1-2Q_{h}5^{\tau}\exp\left(-M\rho_{1}\left(\frac{\min_{1\leq\ell\leq
L}\delta_{\ell}}{48\leavevmode\nobreak\
\lambda_{\max}\left(\boldsymbol{M}^{(h)}\right)}\right)^{2}\right)\quad\text{where}\quad\tau=M^{\frac{1}{2\alpha+1}}$
(4.53)

that the condition (4.24) in our spectral gap assumption 4.1 is satisfied.

###### Proof.

Recall (4.30), i.e.,

$\displaystyle\left\|\widetilde{\boldsymbol{S}}^{(h)}-\widetilde{\boldsymbol{S}}^{(h;M)}\right\|_{2\to
2}\leq\lambda_{\max}\left(\boldsymbol{M}^{(h)}\right)\left\|\boldsymbol{\Sigma}^{(h;M)}_{\boldsymbol{K}^{(h)}}-\boldsymbol{\Sigma}_{\boldsymbol{K}^{(h)}}\right\|_{2\to
2}.$

Thus, it is sufficient to ensure

$\displaystyle\delta_{\ell}\geq
4C_{1}h^{2s}\lambda_{\ell+1}^{-1}+4\left\|\boldsymbol{\Sigma}^{(h;M)}_{\boldsymbol{K}^{(h)}}-\boldsymbol{\Sigma}_{\boldsymbol{K}^{(h)}}\right\|_{2\to
2}\lambda_{\max}\left(\boldsymbol{M}^{(h)}\right)\quad\text{for all
}1\leq\ell\leq L.$

Condition (4.52) reduces this to

$\displaystyle\frac{\min_{1\leq\ell\leq L}\delta_{\ell}}{8\leavevmode\nobreak\
\lambda_{\max}\left(\boldsymbol{M}^{(h)}\right)}\geq\left\|\boldsymbol{\Sigma}^{(h;M)}_{\boldsymbol{K}^{(h)}}-\boldsymbol{\Sigma}_{\boldsymbol{K}^{(h)}}\right\|_{2\to
2}.$

We aim to show that this inequality is only violated with small probability.
To this end, we derive by repeating the arguments of [7, Lemmata 2 & 3] and
[7, Remark 1] the bound

$\displaystyle\mathbb{P}\left\\{\left\|\boldsymbol{\Sigma}^{(h;M)}_{\boldsymbol{K}^{(h)}}-\mathbb{E}\left[\boldsymbol{\Sigma}^{(h;M)}_{\boldsymbol{K}^{(h)}}\right]\right\|_{2\to
2}>t\right\\}\leq
2Q_{h}5^{\tau}\exp\left(-\frac{1}{9}Mt^{2}\rho_{1}\right)\quad\text{for all
}0<t<\rho_{1},$ (4.54)

where $\rho_{1}$ is a constant, see [7, page 2142], and
$\tau=M^{\frac{1}{2\alpha+1}}$. Thus, we obtain with (4.42) and by setting
$\tau=M^{\frac{1}{2\alpha+1}}$ the bound

$\displaystyle\mathbb{P}\left\\{\left\|\boldsymbol{\Sigma}^{(h;M)}_{\boldsymbol{K}^{(h)}}-\boldsymbol{\Sigma}_{\boldsymbol{K}^{(h)}}\right\|_{2\to
2}>t\right\\}$
$\displaystyle\leq\mathbb{P}\left\\{\left\|\boldsymbol{\Sigma}^{(h;M)}_{\boldsymbol{K}^{(h)}}-\mathbb{E}\left[\boldsymbol{\Sigma}^{(h;M)}_{\boldsymbol{K}^{(h)}}\right]\right\|_{2\to
2}+\left\|\mathbb{E}\left[\boldsymbol{\Sigma}^{(h;M)}_{\boldsymbol{K}^{(h)}}\right]-\boldsymbol{\Sigma}_{\boldsymbol{K}^{(h)}}\right\|_{2\to
2}>t\right\\}$
$\displaystyle\leq\mathbb{P}\left\\{\left\|\boldsymbol{\Sigma}^{(h;M)}_{\boldsymbol{K}^{(h)}}-\mathbb{E}\left[\boldsymbol{\Sigma}^{(h;M)}_{\boldsymbol{K}^{(h)}}\right]\right\|_{2\to
2}+CM^{-\frac{2\alpha}{2\alpha+1}}>t\right\\}.$ (4.55)

Now, let $\bar{t}\in(0,\rho_{1})$ be fixed and let $M$ be large enough such
that $CM^{-\frac{2\alpha}{2\alpha+1}}\leq\frac{\bar{t}}{2}$. Then, it holds
that

$\displaystyle\mathbb{P}\left\\{\left\|\boldsymbol{\Sigma}^{(h;M)}_{\boldsymbol{K}^{(h)}}-\boldsymbol{\Sigma}_{\boldsymbol{K}^{(h)}}\right\|_{2\to
2}>\bar{t}\right\\}$
$\displaystyle\leq\mathbb{P}\left\\{\left\|\boldsymbol{\Sigma}^{(h;M)}_{\boldsymbol{K}^{(h)}}-\mathbb{E}\left[\boldsymbol{\Sigma}^{(h;M)}_{\boldsymbol{K}^{(h)}}\right]\right\|_{2\to
2}>\frac{\bar{t}}{2}\right\\}\leq
2Q_{h}5^{\tau}\exp\left(-\frac{1}{36}M\bar{t}^{2}\rho_{1}\right).$ (4.56)

Now, we set $\bar{t}:=\frac{1}{8}\min_{1\leq\ell\leq L}\delta_{\ell}$ and
observe that the second inequality in condition (4.52) ensures
$CM^{-\frac{2\alpha}{2\alpha+1}}\leq\frac{1}{2}\bar{t}=\frac{1}{16\leavevmode\nobreak\
\lambda_{\max}\left(\boldsymbol{M}^{(h)}\right)}\min_{1\leq\ell\leq
L}\delta_{\ell}$. Thus, we get

$\displaystyle\mathbb{P}\left\\{\left\|\boldsymbol{\Sigma}^{(h;M)}_{\boldsymbol{K}^{(h)}}-\boldsymbol{\Sigma}_{\boldsymbol{K}^{(h)}}\right\|_{2\to
2}>\frac{\min_{1\leq\ell\leq L}\delta_{\ell}}{8}\right\\}$
$\displaystyle\leq\mathbb{P}\left\\{\left\|\boldsymbol{\Sigma}^{(h;M)}_{\boldsymbol{K}^{(h)}}-\mathbb{E}\left[\boldsymbol{\Sigma}^{(h;M)}_{\boldsymbol{K}^{(h)}}\right]\right\|_{2\to
2}>\frac{\bar{t}}{2}\right\\}$ $\displaystyle\leq
2Q_{h}5^{\tau}\exp\left(-M\rho_{1}{\left(\frac{\min_{1\leq\ell\leq
L}\delta_{\ell}}{48\leavevmode\nobreak\
\lambda_{\max}\left(\boldsymbol{M}^{(h)}\right)}\right)^{2}}\right),$ (4.57)

using $36*8^{2}=9*2^{8}$ and thus $\sqrt{36*8^{2}}=3*2^{4}=48$. This proves
the claim. ∎

This establishes a further central ingredient for our later analysis, namely
the optimal bound for the sampling error.

##  5 Error analysis for covariance operator reconstruction

We are now in the position to derive an overall error bound for the
reconstruction of the continuous covariance operator from finite samples. To
this end recall the covariance kernel $R$ in (1.1), i.e.,

$\displaystyle
R(\boldsymbol{x},\boldsymbol{x}^{\prime})=\int_{\Omega}\kappa(\boldsymbol{\omega},\boldsymbol{x})\kappa(\boldsymbol{\omega}^{\prime},\boldsymbol{x}^{\prime})\,\mathrm{d}\mathbb{P}=\sum_{\ell=1}^{\infty}\lambda_{\ell}\phi_{\ell}(\boldsymbol{x})\phi_{\ell}(\boldsymbol{x}^{\prime}).$

Furthermore let $\phi_{\ell}^{(h;M)}$ be given as in (4.18) and let
$L\in\mathbb{N}$ be a truncation parameter with $1\leq L\leq Q_{h}$. We define

$\displaystyle
R^{(L;h;M)}(\boldsymbol{x},\boldsymbol{x}^{\prime}):=\sum_{\ell=1}^{L}\lambda^{(h;M)}_{\ell}\phi_{\ell}^{(h;M)}(\boldsymbol{x})\phi_{\ell}^{(h;M)}(\boldsymbol{x}^{\prime}).$
(5.1)

We are now interested in the approximation error $R-R^{(L;h;M)}$. To this end,
let us define

$\displaystyle
G(L):=\left(\sum_{\ell=1}^{L}\left(\frac{\lambda_{\ell}}{\delta_{\ell}}\right)^{2}\right)^{\frac{1}{2}},$
(5.2)

which is a function of the truncation parameter $L$ that is determined by
spectral properties of the unknown covariance kernel $R$ only.
However, in special cases, it is possible to compute $G(L)$ explicitly. As an
example, let us consider Brownian motion. In the univariate case we have the
kernel

$\displaystyle R^{(1)}_{B}:(0,1)^{2}\to\mathbb{R},\quad
R^{(1)}_{B}(x,x^{\prime}):=\min\\{x,x^{\prime}\\}.$ (5.3)

Then it is known that

$\displaystyle\lambda_{R^{(1)}_{B}}(\ell)=\pi^{-2}\left(\ell-\frac{1}{2}\right)^{-2}=\left(\frac{2}{\pi}\right)^{2}\left(2\ell-1\right)^{-2}\in\left(\frac{2}{\pi}\right)^{2}(2\mathbb{N}-1)^{-2}.$

Thus, we have

$\displaystyle\delta_{R^{(1)}_{B}}(\ell):=\frac{2\ell}{\pi^{2}(\ell^{2}-\frac{1}{4})^{2}}=\frac{1}{2}\left(\frac{2}{\pi}\right)^{2}\frac{\ell}{\left(\ell^{2}-\frac{1}{4}\right)^{2}}$

and consequently we obtain

$\displaystyle\frac{\lambda_{R^{(1)}_{B}}(\ell)}{\delta_{R^{(1)}_{B}}(\ell)}=\frac{2}{\ell}\left(\frac{\ell^{2}-\frac{1}{4}}{2\ell-1}\right)^{2}=2\frac{\left(2\ell+1\right)^{2}}{16\ell}=\frac{\left(2\ell+1\right)^{2}}{8\ell}.$

Finally, we get

$\displaystyle
G^{2}_{R^{(1)}_{B}}(L)=\sum_{\ell=1}^{L}\left(\frac{\lambda_{R^{(1)}_{B}}(\ell)}{\delta_{R^{(1)}_{B}}(\ell)}\right)^{2}=\frac{1}{64}\sum_{\ell=1}^{L}\frac{(2\ell+1)^{4}}{\ell^{2}}\approx
L^{3}\quad\text{for }L\text{ large.}$

Moreover, for a two-dimensional Brownian motion, we have the eigenvalues

$\displaystyle\lambda_{R^{(2)}_{B}}(\ell_{1},\ell_{2}):=\lambda_{R^{(1)}_{B}}(\ell_{1})\lambda_{R^{(1)}_{B}}(\ell_{2})=\left(\frac{2}{\pi}\right)^{4}(2\ell_{1}-1)^{-2}(2\ell_{2}-1)^{-2},\quad\ell_{1},\ell_{2}\in\mathbb{N}.$

Hence, we observe

$\displaystyle\text{image}(\lambda_{R^{(2)}_{B}})=\lambda_{R^{(2)}_{B}}(\mathbb{N}\times\mathbb{N})=\left(\frac{2}{\pi}\right)^{4}(2\mathbb{N}-1)^{-2},$

i.e., the image of $\lambda_{R^{(2)}_{B}}$ consists up to the pre-factor
$\left(\frac{2}{\pi}\right)^{4}$ again of the squared inverses of the odd
numbers. Hence, we can sort the eigenvalues (ignoring multiplicities) by the
ordering in the odd natural numbers. Moreover, it holds that

$\displaystyle\lambda_{R^{(2)}_{B}}(\ell_{1},\ell_{2})=\lambda_{R^{(1)}_{B}}(1)\lambda_{R^{(1)}_{B}}\left(k(\ell_{1},\ell_{2})\right),\quad\text{where
}k(\ell_{1},\ell_{2}):=\frac{(2\ell_{1}-1)(2\ell_{2}-1)+1}{2}.$ (5.4)

In the general $d$-dimensional situation, we have

$\displaystyle\lambda_{R^{(d)}_{B}}(\ell_{1},\dots,\ell_{d})=\prod_{j=1}^{d}\lambda_{R^{(d)}_{B}}(\ell_{j})=\left(\lambda_{R^{(1)}_{B}}(1)\right)^{d-1}\lambda^{(1)}(k(\ell_{1},\dots,\ell_{d})),$

where

$\displaystyle
k(\ell_{1},\dots,\ell_{d}):=\frac{1}{2}\left(\prod_{j=1}^{d}(2\ell_{j}-1)+1\right).$

Furthermore, we have with $\ell:=k(\ell_{1},\dots,\ell_{d})$

$\displaystyle\delta_{R^{(d)}_{B}}(\ell)=\left(\lambda_{R^{(1)}_{B}}(1)\right)^{d-1}\delta_{R^{(1)}_{B}}(\ell).$

Thus, for the $d$-dimensional Brownian motion, we get

$\displaystyle G^{2}_{R^{(d)}_{B}}(L)$
$\displaystyle=\sum_{\ell=1}^{L}\left(\frac{\lambda_{R^{(d)}_{B}}(\ell)}{\delta_{R^{(d)}_{B}}(\ell)}\right)^{2}=\sum_{\ell=1}^{L}\left(\frac{\
\lambda_{R^{(1)}_{B}}(\ell)}{\delta_{R^{(1)}_{B}}(\ell)}\right)^{2}=G^{2}_{R^{(1)}_{B}}(L)=\frac{1}{64}\sum_{\ell=1}^{L}\frac{(2\ell+1)^{4}}{\ell^{2}}\approx
L^{3}\quad\text{for }L\text{ large.}$ (5.5)

Note here that $G_{R^{d}_{B}}(L)$ is now completely independent of the
dimension.

In the general situation, we have the following result:

######  Theorem 5.1.

Let $Q_{h}\geq L$. Under the condition of Proposition 4.3, i.e., (4.52) and
the assumption $C_{2}h^{s}\lambda_{L}^{-1}\leq 1$, there holds with
probability $p_{0}$ given in (4.53) that

$\displaystyle\left\|{R}-R^{(L;h;M)}\right\|_{L^{2}(\mathcal{D}\times\mathcal{D})}$
$\displaystyle\lesssim
L^{-\frac{2s}{d}-\frac{1}{2}}+L^{\frac{1}{2}}h^{s}+\left(L^{\frac{1}{2}}+G(L)\right)\left\|\widetilde{\boldsymbol{S}}^{(h)}-\widetilde{\boldsymbol{S}}^{(h;M)}\right\|_{2\to
2}$ (5.6)

with $G(L)$ from (5.2). If we use the assumption
$C_{2}h^{s}\lambda_{L}^{-1}\leq 1$ again, we obtain the simplified bound

$\displaystyle\left\|{R}-R^{(L;h;M)}\right\|_{L^{2}(\mathcal{D}\times\mathcal{D})}$
$\displaystyle\lesssim
L^{-\frac{2s}{d}-\frac{1}{2}}+\left(L^{\frac{1}{2}}+G(L)\right)\left\|\widetilde{\boldsymbol{S}}^{(h)}-\widetilde{\boldsymbol{S}}^{(h;M)}\right\|_{2\to
2}.$ (5.7)

###### Proof.

Since (4.52) holds, it is ensured that the condition (4.24) in our spectral
gap assumption 4.1 is satisfied with the probability $p_{0}$ from (4.53). Now
recall the definition of $R$ in (1.1), of $R^{(L)}$ in (1.4), of $R^{(L;h)}$
in (4.8) and of $R^{(L;h;M)}$ in (5.1). We will split the proof into three
parts. Observe that we can decompose the approximation error as

$\displaystyle R(\boldsymbol{x},\boldsymbol{x}^{\prime})$
$\displaystyle-R^{(L;h;M)}(\boldsymbol{x},\boldsymbol{x}^{\prime})=$ (5.8)
$\displaystyle\underbrace{R(\boldsymbol{x},\boldsymbol{x}^{\prime})-R^{(L)}(\boldsymbol{x},\boldsymbol{x}^{\prime})}_{E_{1}}+\underbrace{R^{(L)}(\boldsymbol{x},\boldsymbol{x}^{\prime})-R^{(L;h)}(\boldsymbol{x},\boldsymbol{x}^{\prime})}_{E_{2}}+\underbrace{R^{(L;h)}(\boldsymbol{x},\boldsymbol{x}^{\prime})-R^{(L;h;M)}(\boldsymbol{x},\boldsymbol{x}^{\prime})}_{E_{3}}.$

The first error term

$\displaystyle
E_{1}(\boldsymbol{x},\boldsymbol{x}^{\prime}):=R(\boldsymbol{x},\boldsymbol{x}^{\prime})-R^{(L)}(\boldsymbol{x},\boldsymbol{x}^{\prime})=\sum_{\ell=L+1}^{\infty}\lambda_{\ell}\phi_{\ell}(\boldsymbol{x})\phi_{\ell}(\boldsymbol{x}^{\prime})$
(5.9)

depends solely on the truncation parameter $L$. The second error term

$\displaystyle
E_{2}(\boldsymbol{x},\boldsymbol{x}^{\prime}):=R^{(L)}(\boldsymbol{x},\boldsymbol{x}^{\prime})-R^{(L;h)}(\boldsymbol{x},\boldsymbol{x}^{\prime})=\sum_{\ell=1}^{L}\lambda_{\ell}\phi_{\ell}(\boldsymbol{x})\phi_{\ell}(\boldsymbol{x}^{\prime})-\sum_{\ell=1}^{L}\lambda^{(h)}_{\ell}\phi^{(h)}_{\ell}(\boldsymbol{x})\phi^{(h)}_{\ell}(\boldsymbol{x}^{\prime})$
(5.10)

depends on the spatial discretization $h$ and on the truncation parameter $L$.
The third error term

$\displaystyle E_{3}(\boldsymbol{x},\boldsymbol{x}^{\prime})$
$\displaystyle:=R^{(L;h)}(\boldsymbol{x},\boldsymbol{x}^{\prime})-R^{(L;h;M)}(\boldsymbol{x},\boldsymbol{x}^{\prime})$
$\displaystyle=\sum_{\ell=1}^{L}\lambda^{(h)}_{\ell}\phi^{(h)}_{\ell}(\boldsymbol{x})\phi^{(h)}_{\ell}(\boldsymbol{x}^{\prime})-\sum_{\ell=1}^{L}\lambda^{(h;M)}_{\ell}\phi^{(h;M)}_{\ell}(\boldsymbol{x})\phi^{(h;M)}_{\ell}(\boldsymbol{x}^{\prime})$
(5.11)

depends on the sample size $M$, the spatial discretization $h$ and on the
truncation parameter $L$.

To derive a bound for (5.9) we use the orthonormality of the eigenfunctions
$\\{\phi_{\ell}\\}_{\ell=1}^{\infty}$ and obtain directly for the integral
operator

$\displaystyle\mathcal{E}_{1}:L^{2}(\mathcal{D})\to L^{2}(\mathcal{D}),\quad
v\mapsto\int_{\mathcal{D}}E_{1}(\boldsymbol{x},\cdot)v(\boldsymbol{x})\,\mathrm{d}\mathcal{L}^{(d)}(\boldsymbol{x})$
(5.12)

the identity

$\displaystyle\left\|\mathcal{E}_{1}\right\|_{L^{2}(\mathcal{D})\to
L^{2}(\mathcal{D})}=\left\|E_{1}\right\|_{L^{2}(\mathcal{D}\times\mathcal{D})}=\left(\sum_{\ell=L+1}^{\infty}\lambda_{\ell}^{2}\right)^{\frac{1}{2}}.$
(5.13)

Now, we can use inequality (3.8) of Theorem 3.1 and derive the bound

$\displaystyle\left(\sum_{\ell=L+1}^{\infty}\lambda_{\ell}^{2}\right)^{\frac{1}{2}}\leq
C_{\ref{thm:truncationError}}\left(\sum_{\ell=L+1}^{\infty}\ell^{-\frac{4s}{d}-2}\right)^{\frac{1}{2}}\leq
C_{\ref{thm:truncationError}}\left(\int_{\ell=L}^{\infty}x^{-\frac{4s}{d}-2}\,\mathrm{d}x\right)^{\frac{1}{2}}\leq
C_{\ref{thm:truncationError}}\left(\frac{4s}{d}+1\right)^{-\frac{1}{2}}L^{-\frac{2s}{d}-\frac{1}{2}}.$

For the second term (5.10), we have

$\displaystyle E_{2}(\boldsymbol{x},\boldsymbol{x}^{\prime})$
$\displaystyle:=\sum_{\ell=1}^{L}\lambda_{\ell}\phi_{\ell}(\boldsymbol{x})\phi_{\ell}(\boldsymbol{x}^{\prime})-\sum_{\ell=1}^{L}\lambda^{(h)}_{\ell}\phi^{(h)}_{\ell}(\boldsymbol{x})\phi^{(h)}_{\ell}(\boldsymbol{x}^{\prime})$
$\displaystyle=\sum_{\ell=1}^{L}\left(\lambda_{\ell}-\lambda^{(h)}_{\ell}\right)\phi_{\ell}(\boldsymbol{x})\phi_{\ell}(\boldsymbol{x}^{\prime})+\sum_{\ell=1}^{L}\lambda^{(h)}_{\ell}\left(\phi_{\ell}(\boldsymbol{x})-\phi^{(h)}_{\ell}(\boldsymbol{x})\right)\phi_{\ell}(\boldsymbol{x}^{\prime})$
$\displaystyle+\sum_{\ell=1}^{L}\lambda^{(h)}_{\ell}\phi^{(h)}_{\ell}(\boldsymbol{x})\left(\phi_{\ell}(\boldsymbol{x}^{\prime})-\phi^{(h)}_{\ell}(\boldsymbol{x}^{\prime})\right)$
$\displaystyle=:E_{2;1}(\boldsymbol{x},\boldsymbol{x}^{\prime})+E_{2;2}(\boldsymbol{x},\boldsymbol{x}^{\prime})+E_{2;3}(\boldsymbol{x},\boldsymbol{x}^{\prime}).$

The orthonormality of $\\{\phi_{\ell}\\}_{\ell=1}^{\infty}$ and an application
of Proposition 4.1 leads to

$\displaystyle\left\|{E}_{2;1}\right\|_{L^{2}\left(\mathcal{D}\times\mathcal{D}\right)}$
$\displaystyle\leq\left(\sum_{\ell=1}^{L}|\lambda_{\ell}-\lambda^{(h)}_{\ell}|^{2}\right)^{1/2}\leq\left(\sum_{\ell=1}^{L}C^{2}_{1}h^{4s}\lambda_{\ell}^{-2}\right)^{1/2}$
$\displaystyle\leq
C_{1}h^{s}\left(\sum_{\ell=1}^{L}h^{2s}\lambda_{\ell}^{-2}\right)^{1/2}\leq\frac{C_{1}}{C_{2}}L^{\frac{1}{2}}h^{s},$
(5.14)

where the last inequality holds since $C_{2}h^{s}\lambda_{L}^{-1}\leq 1$. This
bound implies

$\displaystyle C_{2}h^{s}\leq\lambda_{L}\leq L^{-\frac{2s}{d}-1}.$ (5.15)

Now we turn to $E_{2;2}$. The orthonormality of
$\\{\phi_{\ell}\\}_{\ell=1}^{\infty}$ and the bound
$\lambda^{(h)}_{\ell}\leq\lambda_{\ell}$ for $1\leq\ell\leq L$ due to the
Courant min-max principle for decreasing eigenvalues implies

$\displaystyle\left\|{E}_{2;2}\right\|_{L^{2}\left(\mathcal{D}\times\mathcal{D}\right)}^{2}\leq\sum_{\ell=1}^{L}\left(\lambda^{(h)}_{\ell}\right)^{2}\left\|\phi_{\ell}-\phi_{\ell}^{(h)}\right\|_{L^{2}\left(\mathcal{D}\right)}^{2}\leq\sum_{\ell=1}^{L}\lambda^{2}_{\ell}\left\|\phi_{\ell}-\phi_{\ell}^{(h)}\right\|_{L^{2}\left(\mathcal{D}\right)}^{2}\leq
C_{2}^{2}h^{2s}L,$

where we used
$\left\|\phi_{\ell}^{(h)}-\phi_{\ell}\right\|_{L^{2}\left(\mathcal{D}\right)}\leq
C_{2}\lambda_{\ell}^{-1}h^{s}$ in the last step, see Proposition 4.1. Thus, we
obtain

$\displaystyle\left\|{E}_{2;2}\right\|_{L^{2}\left(\mathcal{D}\times\mathcal{D}\right)}$
$\displaystyle\leq C_{2}h^{s}L^{\frac{1}{2}}.$ (5.16)

For the term $E_{2,3}$ the orthonormality of the series
$\\{\phi_{\ell}^{(h)}\\}_{\ell=1}^{Q_{h}}$, which stems from the eigenvalue
problem (4.4), leads to

$\displaystyle\left\|{E}_{2;3}\right\|_{L^{2}\left(\mathcal{D}\times\mathcal{D}\right)}^{2}$
$\displaystyle=\sum_{\ell=1}^{L}\left(\lambda^{(h)}_{\ell}\right)^{2}\left\|\phi_{\ell}-\phi_{\ell}^{(h)}\right\|_{L^{2}\left(\mathcal{D}\right)}^{2}.$

Thus, we obtain with verbatim the same computations as for (5.16), the bound

$\displaystyle\left\|{E}_{2;3}\right\|_{L^{2}\left(\mathcal{D}\times\mathcal{D}\right)}$
$\displaystyle\leq C_{2}h^{s}L^{\frac{1}{2}}.$ (5.17)

Altogether with the estimates (5), (5.16) and (5.17), we derive

$\displaystyle\left\|{E}_{2}\right\|_{L^{2}\left(\mathcal{D}\times\mathcal{D}\right)}$
$\displaystyle=\left\|{E}_{2;1}+{E}_{2;2}+{E}_{2;3}\right\|_{L^{2}\left(\mathcal{D}\times\mathcal{D}\right)}\leq\left\|{E}_{2;1}\right\|_{L^{2}\left(\mathcal{D}\times\mathcal{D}\right)}+\left\|{E}_{2;2}\right\|_{L^{2}\left(\mathcal{D}\times\mathcal{D}\right)}+\left\|{E}_{2;3}\right\|_{L^{2}\left(\mathcal{D}\times\mathcal{D}\right)}$
$\displaystyle\leq\left(\frac{C_{1}}{C_{2}}+2C_{2}\right)h^{s}L^{\frac{1}{2}}\leq\left(\frac{C_{1}}{C_{2}}+2C_{2}\right)C^{-1}_{2}L^{-\frac{2s}{d}-1}L^{\frac{1}{2}}=\left(\frac{C_{1}}{C_{2}}+2C_{2}\right)C^{-1}_{2}L^{-\frac{2s}{d}-\frac{1}{2}}.$
(5.18)

For the third term $E_{3}$ in (5), we proceed in an analogous fashion as in
the splitting of (5.10). We then have the decomposition

$\displaystyle
E_{3}(\boldsymbol{x},\boldsymbol{x}^{\prime})=\sum_{\ell=1}^{L}\lambda^{(h)}_{\ell}\phi_{\ell}^{(h)}(\boldsymbol{x})\phi_{\ell}^{(h)}(\boldsymbol{x}^{\prime})-\sum_{\ell=1}^{L}\lambda^{(h;M)}_{\ell}\phi^{(h;M)}_{\ell}(\boldsymbol{x})\phi^{(h;M)}_{\ell}(\boldsymbol{x}^{\prime})$
$\displaystyle=\sum_{\ell=1}^{L}\left(\lambda^{(h)}_{\ell}-\lambda^{(h;M)}_{\ell}\right)\phi_{\ell}^{(h)}(\boldsymbol{x})\phi_{\ell}^{(h)}(\boldsymbol{x}^{\prime})+\sum_{\ell=1}^{L}\lambda^{(h;M)}_{\ell}\left(\phi_{\ell}^{(h)}(\boldsymbol{x})-\phi^{(h;M)}_{\ell}(\boldsymbol{x})\right)\phi^{(h)}_{\ell}(\boldsymbol{x}^{\prime})$
$\displaystyle+\sum_{\ell=1}^{L}\lambda^{(h;M)}_{\ell}\left(\phi_{\ell}^{(h)}(\boldsymbol{x}^{\prime})-\phi_{\ell}^{(h;M)}(\boldsymbol{x}^{\prime})\right)\phi^{(h;M)}_{\ell}(\boldsymbol{x})$
$\displaystyle=:E_{3;1}+E_{3;2}+E_{3;3}.$

Now, for the term $E_{3,1}$, the orthonormality of the basis
$\\{\phi_{\ell}^{(h)}\\}_{\ell=1}^{Q_{h}}$ and an application of (4.20)
implies

$\displaystyle\left\|{E}_{3;1}\right\|_{L^{2}\left(\mathcal{D}\right)}$
$\displaystyle=\left(\sum_{\ell=1}^{L}\left(\lambda^{(h)}_{\ell}-\lambda^{(h;M)}_{\ell}\right)^{2}\right)^{\frac{1}{2}}\leq
L^{1/2}\left\|\widetilde{\boldsymbol{S}}^{(h)}-\widetilde{\boldsymbol{S}}^{(h;M)}\right\|_{2\to
2}.$ (5.19)

To derive an upper estimate for the term $E_{3,2}$, we will use our bound from
Theorem 4.1 for the eigenfunction approximation in (4.27). Thus we have with
probability $p_{0}$ that

$\displaystyle\left\|\phi^{(h)}_{\ell}-\phi^{(h;M)}_{\ell}\right\|_{L^{2}\left(\mathcal{D}\right)}^{2}\leq
4C^{2}\frac{\left\|\widetilde{\boldsymbol{S}}^{(h)}-\widetilde{\boldsymbol{S}}^{(h;M)}\right\|^{2}_{2\to
2}}{\delta_{\ell}^{2}},$ (5.20)

where we note at this point that this bound makes use of (4.24) from the
spectral gap assumption 4.1. Due to the orthonormality of the basis
$\\{\phi^{(h)}_{\ell}\\}_{\ell=1}^{Q_{h}}$, this yields

$\displaystyle\left\|{E}_{3;2}\right\|_{L^{2}\left(\mathcal{D}\times\mathcal{D}\right)}^{2}$
$\displaystyle=\sum_{\ell=1}^{L}(\lambda^{(h;M)}_{\ell})^{2}\left\|\phi^{(h)}_{\ell}-\phi^{(h;M)}_{\ell}\right\|_{L^{2}\left(\mathcal{D}\right)}^{2}\leq
4C^{2}\sum_{\ell=1}^{L}\left(\frac{\left\|\widetilde{\boldsymbol{S}}^{(h)}-\widetilde{\boldsymbol{S}}^{(h;M)}\right\|_{2\to
2}}{\delta_{\ell}}\right)^{2}\left(\lambda^{(h;M)}_{\ell}\right)^{2}.$

Hence, we obtain

$\displaystyle\left\|{E}_{3;2}\right\|_{L^{2}\left(\mathcal{D}\times\mathcal{D}\right)}\leq
2C\left\|\widetilde{\boldsymbol{S}}^{(h)}-\widetilde{\boldsymbol{S}}^{(h;M)}\right\|_{2\to
2}\left(\sum_{\ell=1}^{L}\left(\frac{\lambda^{(h;M)}_{\ell}}{\delta_{\ell}}\right)^{2}\right)^{\frac{1}{2}}.$

Moreover, we can use (4.20), i.e.,

$\displaystyle\left|\lambda^{(h)}_{\ell}-\lambda^{(h;M)}_{\ell}\right|\leq\left\|\widetilde{\boldsymbol{S}}^{(h)}-\widetilde{\boldsymbol{S}}^{(h;M)}\right\|_{2\to
2},\quad 1\leq\ell\leq\min\\{Q_{h},L\\}.$

to get

$\displaystyle\left|\lambda^{(h;M)}_{\ell}\right|\leq\left|\lambda^{(h)}_{\ell}-\lambda^{(h;M)}_{\ell}\right|+\left|\lambda^{(h)}_{\ell}\right|\leq\left\|\widetilde{\boldsymbol{S}}^{(h)}-\widetilde{\boldsymbol{S}}^{(h;M)}\right\|_{2\to
2}+\left|\lambda_{\ell}\right|,$ (5.21)

where we used again $\lambda^{(h)}_{\ell}\leq\lambda_{\ell}$ for
$1\leq\ell\leq L$ in the last step. This yields with probability $p_{0}$

$\displaystyle\left\|{E}_{3;2}\right\|_{L^{2}\left(\mathcal{D}\times\mathcal{D}\right)}\leq
2C\left\|\widetilde{\boldsymbol{S}}^{(h)}-\widetilde{\boldsymbol{S}}^{(h;M)}\right\|_{2\to
2}\left(\sum_{\ell=1}^{L}\left(\frac{\left\|\widetilde{\boldsymbol{S}}^{(h)}-\widetilde{\boldsymbol{S}}^{(h;M)}\right\|_{2\to
2}+\left|\lambda_{\ell}\right|}{\delta_{\ell}}\right)^{2}\right)^{\frac{1}{2}}.$
(5.22)

For the term $E_{3,3}$, we use orthonormality of the
$\\{\phi_{\ell}^{(h;M)}\\}$ to get

$\displaystyle\left\|E_{3;3}\right\|_{L^{2}\left(\mathcal{D}\times\mathcal{D}\right)}^{2}$
$\displaystyle=\sum_{\ell=1}^{L}(\lambda^{(h;M)}_{\ell})^{2}\left\|\phi^{(h)}_{\ell}-\phi^{(h;M)}_{\ell}\right\|_{L^{2}\left(\mathcal{D}\right)}^{2}=\left\|E_{3;2}\right\|_{L^{2}\left(\mathcal{D}\times\mathcal{D}\right)}^{2}.$
(5.23)

Thus, we again obtain probability $p_{0}$

$\displaystyle\left\|{E}_{3;3}\right\|_{L^{2}\left(\mathcal{D}\times\mathcal{D}\right)}\leq
2C\left\|\widetilde{\boldsymbol{S}}^{(h)}-\widetilde{\boldsymbol{S}}^{(h;M)}\right\|_{2\to
2}\left(\sum_{\ell=1}^{L}\left(\frac{\left\|\widetilde{\boldsymbol{S}}^{(h)}-\widetilde{\boldsymbol{S}}^{(h;M)}\right\|_{2\to
2}+\left|\lambda_{\ell}\right|}{\delta_{\ell}}\right)^{2}\right)^{\frac{1}{2}}.$
(5.24)

Altogether, combining the inequalities (5.19), (5.22) and (5.24) we obtain
probability $p_{0}$

$\displaystyle\left\|{E}_{3}\right\|_{L^{2}\left(\mathcal{D}\times\mathcal{D}\right)}=\left\|{E}_{3;1}+{E}_{3;2}+{E}_{3;3}\right\|_{L^{2}\left(\mathcal{D}\times\mathcal{D}\right)}\leq\left\|{E}_{3;1}\right\|_{L^{2}\left(\mathcal{D}\times\mathcal{D}\right)}+\left\|{E}_{3;2}\right\|_{L^{2}\left(\mathcal{D}\times\mathcal{D}\right)}+\left\|{E}_{3;3}\right\|_{L^{2}\left(\mathcal{D}\times\mathcal{D}\right)}$
$\displaystyle\leq\left\|\widetilde{\boldsymbol{S}}^{(h)}-\widetilde{\boldsymbol{S}}^{(h;M)}\right\|_{2\to
2}\left(L^{\frac{1}{2}}+4C\left(\sum_{\ell=1}^{L}\left(\frac{\left\|\widetilde{\boldsymbol{S}}^{(h)}-\widetilde{\boldsymbol{S}}^{(h;M)}\right\|_{2\to
2}+\left|\lambda_{\ell}\right|}{\delta_{\ell}}\right)^{2}\right)^{\frac{1}{2}}\right).$
(5.25)

Thus, we get with (4.24) from the spectral gap assumption 4.1

$\displaystyle\delta_{\ell}\geq
4C_{1}h^{2s}\lambda_{\ell+1}^{-1}+4\left\|\widetilde{\boldsymbol{S}}^{(h)}-\widetilde{\boldsymbol{S}}^{(h;M)}\right\|_{2\to
2}\geq
4\left\|\widetilde{\boldsymbol{S}}^{(h)}-\widetilde{\boldsymbol{S}}^{(h;M)}\right\|_{2\to
2},$ (5.26)

which is satisfied with probability $p_{0}$, that

$\displaystyle\sum_{\ell=1}^{L}\left(\frac{\left\|\widetilde{\boldsymbol{S}}^{(h)}-\widetilde{\boldsymbol{S}}^{(h;M)}\right\|_{2\to
2}+\left|\lambda_{\ell}\right|}{\delta_{\ell}}\right)^{2}\leq
2\sum_{\ell=1}^{L}\left(\frac{\left\|\widetilde{\boldsymbol{S}}^{(h)}-\widetilde{\boldsymbol{S}}^{(h;M)}\right\|_{2\to
2}}{\delta_{\ell}}\right)^{2}+2\sum_{\ell=1}^{L}\left(\frac{\left|\lambda_{\ell}\right|}{\delta_{\ell}}\right)^{2}$
$\displaystyle\leq\frac{1}{8}L+2\sum_{\ell=1}^{L}\left(\frac{\lambda_{\ell}}{\delta_{\ell}}\right)^{2}=\frac{1}{8}L+2G^{2}(L).$
(5.27)

Thus, we obtain, using $\sqrt{a+b}\leq\sqrt{a}+\sqrt{b}$ for $a,b\geq 0$, that
the bound (5) reduces to

$\displaystyle\left\|{E}_{3}\right\|_{L^{2}\left(\mathcal{D}\times\mathcal{D}\right)}\leq\left\|\widetilde{\boldsymbol{S}}^{(h)}-\widetilde{\boldsymbol{S}}^{(h;M)}\right\|_{2\to
2}\left(L^{\frac{1}{2}}+\sqrt{2}CL^{\frac{1}{2}}+4\sqrt{2}CG(L)\right),$
(5.28)

which holds with probability $p_{0}$. Finally, combining the inequalities
(5.13), (5) and (5.28) we obtain

$\displaystyle\left\|R-R^{(L;h;M)}\right\|_{\mathcal{D}\times\mathcal{D}}\leq\left\|E_{1}+E_{2}+E_{3}\right\|_{\mathcal{D}\times\mathcal{D}}\leq\left\|E_{1}\right\|_{\mathcal{D}\times\mathcal{D}}+\left\|E_{2}\right\|_{\mathcal{D}\times\mathcal{D}}+\left\|E_{3}\right\|_{\mathcal{D}\times\mathcal{D}},$

which shows the assertion. ∎

The bound in Theorem 5.1 still depends on
$\left\|\widetilde{\boldsymbol{S}}^{(h)}-\widetilde{\boldsymbol{S}}^{(h;M)}\right\|_{2\to
2}$. This term however involves the $M$-dependent sampling approximation
(4.13) and is thus a random variable. In the following, we therefore are
interested in the expected value of the error. To this end, note that we can
not use Corollary 4.1 directly as we need the spectral gap assumption 4.1 to
be valid. But this is only the case with probability $p_{0}$, which needs to
be properly be taken care of. Moreover, we need to introduce the quantity

$\displaystyle H(L):={\left(\frac{1}{48}\min_{1\leq\ell\leq
L}\delta_{\ell}\right)^{2}},$ (5.29)

which is a function of the truncation parameter $L$ that depends on spectral
properties as given in (4.23) of the true kernel $R$ . For general operators
this quantity is not easy to obtain.

However, in the simple case of $d$-dimensional Brownian motion we can obtain a
precise value as

$\displaystyle H_{R^{(d)}_{B}}(L)$
$\displaystyle:=\frac{1}{{48^{2}}}\min_{\ell=1}^{L}\delta^{2}_{R^{(d)}_{B}}(\ell)=\frac{1}{{2304}}\min_{\ell=1}^{L}\left(\lambda_{R^{(1)}_{B}}(1)\right)^{2(d-1)}\left(\delta_{R^{(1)}_{B}}(\ell)\right)^{2}$
$\displaystyle=\frac{1}{{2304}}\left(\lambda_{R^{(1)}_{B}}(1)\right)^{2(d-1)}\frac{1}{4}\left(\frac{2}{\pi}\right)^{4}\min_{\ell=1}^{L}\left(\frac{\ell}{\left(\ell^{2}-\frac{1}{4}\right)^{2}}\right)^{2}=\frac{1}{{9216}}\left(\frac{2}{\pi}\right)^{4d}\min_{\ell=1}^{L}\left(\frac{\ell}{\left(\ell^{2}-\frac{1}{4}\right)^{2}}\right)^{2}$
$\displaystyle\approx L^{-6}\quad\text{for }L\text{ large.}$ (5.30)

Note here that, in contrast to the value $G_{R^{(d)}_{B}}(L)$ in (5.5), the
dimension enters the overall value of $H_{R^{(d)}_{B}}(L)$ in (5) now
exponentially.

Altogether, in the general situation, we have the following result.

######  Theorem 5.2.

Under the conditions of Proposition 4.3, i.e. (4.52), with condition (5.15)
and with definition (4.48), there holds

$\displaystyle\mathbb{E}\left[\left\|{R}-R^{(L;h;M)}\right\|_{L^{2}(\mathcal{D}\times\mathcal{D})}\right]$
$\displaystyle\lesssim
L^{-\frac{2s}{d}-\frac{1}{2}}+\left(L^{\frac{1}{2}}+G(L)\right)\tilde{\rho}^{\frac{1}{2}}_{h}(M)\lambda_{\max}\left(\boldsymbol{M}^{(h)}\right)$
$\displaystyle+L^{\frac{1}{2}}h^{-d}\exp\left(-M\rho_{1}H(L)\lambda_{\max}^{-2}\left(\boldsymbol{M}^{(h)}\right)\right),$
(5.31)

where $H(L)$ is defined in (5.29) and $\rho_{1}$ denotes a constant given in
(4.54).

###### Proof.

From Theorem 5.1, we have with condition (5.15) that

$\displaystyle\mathbb{E}\left[\left\|{R}-R^{(L;h;M)}\right\|_{L^{2}(\mathcal{D}\times\mathcal{D})}\right]\lesssim
L^{-\frac{2s}{d}-\frac{1}{2}}+h^{s}L^{\frac{1}{2}}+\mathbb{E}\left[\left\|E_{3}\right\|_{L^{2}\left(\mathcal{D}\times\mathcal{D}\right)}\right]$
$\displaystyle\lesssim
L^{-\frac{2s}{d}-\frac{1}{2}}+\mathbb{E}\left[\left\|E_{3}\right\|_{L^{2}\left(\mathcal{D}\times\mathcal{D}\right)}\right]$
$\displaystyle\lesssim
L^{-\frac{2s}{d}-\frac{1}{2}}+\mathbb{E}\left[\left\|E_{3;1}\right\|_{L^{2}\left(\mathcal{D}\times\mathcal{D}\right)}+\left\|E_{3;2}\right\|_{L^{2}\left(\mathcal{D}\times\mathcal{D}\right)}+\left\|E_{3;3}\right\|_{L^{2}\left(\mathcal{D}\times\mathcal{D}\right)}\right].$

The term
$\left\|{E}_{3;1}\right\|_{L^{2}\left(\mathcal{D}\times\mathcal{D}\right)}$
from (5.19) can be handled directly and we obtain with Corollary 4.1 the bound

$\displaystyle\mathbb{E}\left[\left\|{E}_{3;1}\right\|_{L^{2}\left(\mathcal{D}\times\mathcal{D}\right)}\right]$
$\displaystyle\leq
L^{1/2}\mathbb{E}\left[\left\|\widetilde{\boldsymbol{S}}^{(h)}-\widetilde{\boldsymbol{S}}^{(h;M)}\right\|_{2\to
2}\right]\lesssim L^{1/2}\rho^{\frac{1}{2}}_{h}(M)$ (5.32)

Next, we define the set, where the assumption 4.1 is violated, as bad and the
good set simply as the complement of the bad set, i.e.,

$\displaystyle M_{\text{bad}}$
$\displaystyle:=\left\\{\boldsymbol{\omega}\in\Omega:\frac{\min_{1\leq\ell\leq
L}\delta_{\ell}}{4}<\left\|\widetilde{\boldsymbol{S}}^{(h)}-\widetilde{\boldsymbol{S}}^{(h;M)}\right\|_{2\to
2}\right\\},$ $\displaystyle M_{\text{good}}$
$\displaystyle:=\left\\{\boldsymbol{\omega}\in\Omega:\frac{\min_{1\leq\ell\leq
L}\delta_{\ell}}{4}\geq\left\|\widetilde{\boldsymbol{S}}^{(h)}-\widetilde{\boldsymbol{S}}^{(h;M)}\right\|_{2\to
2}\right\\}.$

We note that Proposition 4.3 yields

$\displaystyle\left|M_{\text{bad}}\right|\leq
2Q_{h}5^{\tau}\exp\left(-M\rho_{1}H(L)\lambda_{\max}^{-2}\left(\boldsymbol{M}^{(h)}\right)\right)\quad\text{where}\quad\tau=M^{\frac{1}{2\alpha+1}}.$

Now, we split the term
$\mathbb{E}\left[\left\|{E}_{3;j}\right\|_{L^{2}\left(\mathcal{D}\times\mathcal{D}\right)}\right]$
as follows

$\displaystyle\mathbb{E}\left[\left\|{E}_{3;j}\right\|_{L^{2}\left(\mathcal{D}\times\mathcal{D}\right)}\right]=\int_{\boldsymbol{\omega}\in
M_{\text{bad}}}\left\|{E}_{3;j}\right\|_{L^{2}\left(\mathcal{D}\times\mathcal{D}\right)}\,\mathrm{d}\mathbb{P}(\boldsymbol{\omega})+\int_{\boldsymbol{\omega}\in
M_{\text{good}}}\left\|{E}_{3;j}\right\|_{L^{2}\left(\mathcal{D}\times\mathcal{D}\right)}\,\mathrm{d}\mathbb{P}(\boldsymbol{\omega}),$

and we also recall that
$\left\|{E}_{3;2}\right\|_{L^{2}\left(\mathcal{D}\times\mathcal{D}\right)}=\left\|{E}_{3;3}\right\|_{L^{2}\left(\mathcal{D}\times\mathcal{D}\right)}$
as we have seen in (5.23). In order to treat the first integral relating to
the bad set, we will show that the integrand is point-wise bounded by a
deterministic quantity on $M_{\text{bad}}$, i.e., that
$\left\|{E}_{3;j}\right\|_{L^{2}\left(\mathcal{D}\times\mathcal{D}\right)}\leq
C$ holds. First, we use (4.2) to get at least

$\displaystyle\left\|\phi^{(h)}_{\ell}-\phi^{(h;M)}_{\ell}\right\|_{L^{2}\left(\mathcal{D}\right)}^{2}$
$\displaystyle\leq
2\left\|\phi^{(h)}_{\ell}\right\|_{L^{2}\left(\mathcal{D}\right)}^{2}+2\left\|\phi^{(h;M)}_{\ell}\right\|_{L^{2}\left(\mathcal{D}\right)}^{2}\leq
4.$

Moreover, we use again (5.21), i.e.,

$\displaystyle\left|\lambda^{(h;M)}_{\ell}\right|\leq\left|\lambda^{(h)}_{\ell}-\lambda^{(h;M)}_{\ell}\right|+\left|\lambda^{(h)}_{\ell}\right|\leq\left\|\widetilde{\boldsymbol{S}}^{(h)}-\widetilde{\boldsymbol{S}}^{(h;M)}\right\|_{2\to
2}+\left|\lambda_{\ell}\right|.$

Thus, we obtain for $j=2,3$ the bound

$\displaystyle\left\|{E}_{3;j}\right\|_{L^{2}\left(\mathcal{D}\times\mathcal{D}\right)}^{2}$
$\displaystyle=\sum_{\ell=1}^{L}(\lambda^{(h;M)}_{\ell})^{2}\left\|\phi^{(h)}_{\ell}-\phi^{(h;M)}_{\ell}\right\|_{L^{2}\left(\mathcal{D}\right)}^{2}\leq
4\sum_{\ell=1}^{L}(\lambda^{(h;M)}_{\ell})^{2}$ $\displaystyle\leq
4\sum_{\ell=1}^{L}\left(\left\|\widetilde{\boldsymbol{S}}^{(h)}-\widetilde{\boldsymbol{S}}^{(h;M)}\right\|_{2\to
2}+\left|\lambda_{\ell}\right|\right)^{2}\leq
8\sum_{\ell=1}^{L}\left\|\widetilde{\boldsymbol{S}}^{(h)}-\widetilde{\boldsymbol{S}}^{(h;M)}\right\|^{2}_{2\to
2}+8\sum_{\ell=1}^{L}\left|\lambda_{\ell}\right|^{2}$ $\displaystyle\leq
8L\left\|\widetilde{\boldsymbol{S}}^{(h)}-\widetilde{\boldsymbol{S}}^{(h;M)}\right\|^{2}_{2\to
2}+8L\left|\lambda_{1}\right|^{2}$

and get

$\displaystyle\left\|{E}_{3;j}\right\|_{L^{2}\left(\mathcal{D}\times\mathcal{D}\right)}\leq
2\sqrt{2}\
L^{\frac{1}{2}}\left(\left\|\widetilde{\boldsymbol{S}}^{(h)}-\widetilde{\boldsymbol{S}}^{(h;M)}\right\|_{2\to
2}+\lambda_{1}\right)$

using $\sqrt{a+b}\leq\sqrt{a}+\sqrt{b}$ for $a,b\geq 0$. Hence, we obtain on
the bad set

$\displaystyle\int_{\boldsymbol{\omega}\in
M_{\text{bad}}}\left\|{E}_{3;j}\right\|_{L^{2}\left(\mathcal{D}\times\mathcal{D}\right)}\,\mathrm{d}\mathbb{P}(\boldsymbol{\omega})\leq
2\sqrt{2}\ L^{\frac{1}{2}}\int_{\boldsymbol{\omega}\in
M_{\text{bad}}}\left\|\widetilde{\boldsymbol{S}}^{(h)}-\widetilde{\boldsymbol{S}}^{(h;M)}\right\|_{2\to
2}\,\mathrm{d}\mathbb{P}(\boldsymbol{\omega})+2\sqrt{2}\
L^{\frac{1}{2}}\left|M_{\text{bad}}\right|\lambda_{1}.$

For the first term, we use (4.49) to get

$\displaystyle\int_{\boldsymbol{\omega}\in
M_{\text{bad}}}\left\|\widetilde{\boldsymbol{S}}^{(h)}-\widetilde{\boldsymbol{S}}^{(h;M)}\right\|_{2\to
2}\,\mathrm{d}\mathbb{P}(\boldsymbol{\omega})$
$\displaystyle\leq\int_{\boldsymbol{\omega}\in\Omega}\left\|\widetilde{\boldsymbol{S}}^{(h)}-\widetilde{\boldsymbol{S}}^{(h;M)}\right\|_{2\to
2}\,\mathrm{d}\mathbb{P}(\boldsymbol{\omega})\lesssim\rho^{\frac{1}{2}}_{h}(M)\lambda_{\max}\left(\boldsymbol{M}^{(h)}\right)$

and thus obtain

$\displaystyle\int_{\boldsymbol{\omega}\in
M_{\text{bad}}}\left\|{E}_{3;j}\right\|_{L^{2}\left(\mathcal{D}\times\mathcal{D}\right)}\,\mathrm{d}\mathbb{P}(\boldsymbol{\omega})\lesssim
2\sqrt{2}L^{\frac{1}{2}}\lambda_{\max}\left(\boldsymbol{M}^{(h)}\right)\rho^{\frac{1}{2}}_{h}(M)+2\sqrt{2}\left|M_{\text{bad}}\right|L^{\frac{1}{2}}\lambda_{1}.$

We are left with estimating the integral of
$\left\|E_{3}\right\|_{L^{2}\left(\mathcal{D}\times\mathcal{D}\right)}$ over
the set $M_{\text{good}}$. To this end, we employ the bound (5.28), which is
only valid on the good set as we have the condition (5.26). Recall now

$\displaystyle\left\|{E}_{3}\right\|_{L^{2}\left(\mathcal{D}\times\mathcal{D}\right)}\leq\left\|\widetilde{\boldsymbol{S}}^{(h)}-\widetilde{\boldsymbol{S}}^{(h;M)}\right\|_{2\to
2}\left(L^{\frac{1}{2}}+\sqrt{2}CL^{\frac{1}{2}}+4\sqrt{2}CG(L)\right)\quad\text{on
}M_{\text{good}}.$

Thus we get

$\displaystyle\int_{M_{\text{good}}}\left\|E_{3}\right\|_{L^{2}\left(\mathcal{D}\times\mathcal{D}\right)}\,\mathrm{d}\mathbb{P}(\boldsymbol{\omega})$
$\displaystyle\leq\left(L^{\frac{1}{2}}+\sqrt{2}CL^{\frac{1}{2}}+4\sqrt{2}CG(L)\right)\int_{M_{\text{good}}}\left\|\widetilde{\boldsymbol{S}}^{(h)}-\widetilde{\boldsymbol{S}}^{(h;M)}\right\|_{2\to
2}\,\mathrm{d}\mathbb{P}(\boldsymbol{\omega})$
$\displaystyle\leq\left(L^{\frac{1}{2}}+\sqrt{2}CL^{\frac{1}{2}}+4\sqrt{2}CG(L)\right)\mathbb{E}\left[\left\|\widetilde{\boldsymbol{S}}^{(h)}-\widetilde{\boldsymbol{S}}^{(h;M)}\right\|_{2\to
2}\right]$
$\displaystyle\lesssim\left(L^{\frac{1}{2}}+\sqrt{2}CL^{\frac{1}{2}}+4\sqrt{2}CG(L)\right)\rho^{\frac{1}{2}}_{h}(M)\lambda_{\max}\left(\boldsymbol{M}^{(h)}\right).$

Together with Proposition 4.3 and the notation from (4.48) this completes the
proof. ∎

##  6 Discussion

It now remains to put the result of Theorem 5.2 into context. To this end, let
us aim at an error bound

$\mathbb{E}\left[\|R-R^{(L,h,M)}\|_{L^{2}(\mathcal{D}\times\mathcal{D})}\right]\leq
c\varepsilon$

with a prescribed, fixed accuracy $\varepsilon>0$ and a small constant $c$.
The question is then: What are the conditions on the discretization parameters
$L,h,M$ to achieve this aim? For the sake of simplicity, we restrict ourselves
here to the situation
$\lambda_{\max}\left(\boldsymbol{M}^{(h)}\right)\approx\lambda_{\min}\left(\boldsymbol{M}^{(h)}\right)\approx
1$, which is satisfied for instance for any orthonormal basis.111For the nodal
basis, we would obtain, after proper scaling,
$\lambda_{\max}\left(\boldsymbol{M}^{(h)}\right)\approx\lambda_{\min}\left(\boldsymbol{M}^{(h)}\right)\approx
h^{d}$. The resulting discussion for this case is left to the reader.
Moreover, we assume the smoothness assumption (1.9), i.e., $s>d/2$. We proceed
in a term by term fashion as follows:

For the first term in our error bound (5.2) we obtain from
$L^{-2s/d-1/2}\stackrel{{\scriptstyle!}}{{\leq}}\varepsilon$ directly the
condition $L\geq\varepsilon^{-2d/(4s+d)}$. We also point out that choosing $L$
much larger will destroy the balance of the error contributions. Hence, we
will assume

$\displaystyle L_{\varepsilon}:=\lceil\varepsilon^{-\frac{2d}{4s+d}}\rceil.$
(6.1)

The coupling (6.1) influences the choices for $h$ in terms of $L$ via the
Proposition 4.3, i.e. (4.52) and the relation (5.15). We get as sufficient
condition on $h$

$\displaystyle
h^{2s}\lesssim\min\left\\{H^{\frac{1}{2}}(L)\lambda_{L+1},\lambda^{2}_{L}\right\\}\lesssim\min\left\\{\delta_{\ell}\lambda_{L+1},\lambda^{2}_{L}\
:\ 1\leq\ell\leq L\right\\}$ (6.2)

with the notation (5.29).

For the second error term we now assume $G(L)=L^{\tilde{\gamma}}$ for some
$\tilde{\gamma}\in\mathbb{R}$. Thus

$\displaystyle L^{\frac{1}{2}}+G(L)\lesssim
L^{\gamma},\quad\text{with}\quad\gamma:=\max\\{\frac{1}{2},\tilde{\gamma}\\}.$
(6.3)

Consequently, we have
$(L^{\frac{1}{2}}+G(L))\tilde{\rho}^{\frac{1}{2}}_{h}(M)\lesssim
L^{\gamma}\tilde{\rho}^{\frac{1}{2}}_{h}(M)$. At this point, we note that the