remarkRemark exampleExample assumptionAssumption Asymptotic analysis of
multilevel best linear unbiased estimatorsDaniel Schaden and Elisabeth Ullmann

# Asymptotic Analysis of Multilevel Best Linear Unbiased Estimators††thanks:
Submitted to the editors DATE. The authors gratefully acknowledge the support
by the Deutsche Forschungsgemeinschaft (DFG) through the International
Research Training Group IGDK 1754 “Optimization and Numerical Analysis for
Partial Differential Equations with Nonsmooth Structures” Projektnummer
188264188/GRK1754.

Daniel Schaden Department of Mathematics, Technical University of Munich,
Boltzmannstr. 3, 85748 Garching b. München, Germany<EMAIL_ADDRESS><EMAIL_ADDRESS>Elisabeth Ullmann22footnotemark: 2

###### Abstract

We study the computational complexity and variance of multilevel best linear
unbiased estimators introduced in [D. Schaden and E. Ullmann, SIAM/ASA J.
Uncert. Quantif., (2020)]. We specialize the results in this work to PDE-based
models that are parameterized by a discretization quantity, e.g., the finite
element mesh size. In particular, we investigate the asymptotic complexity of
the so-called sample allocation optimal best linear unbiased estimators
(SAOBs). These estimators have the smallest variance given a fixed
computational budget. However, SAOBs are defined implicitly by solving an
optimization problem and are difficult to analyze. Alternatively, we study a
class of auxiliary estimators based on the Richardson extrapolation of the
parametric model family. This allows us to provide an upper bound for the
complexity of the SAOBs, showing that their complexity is optimal within a
certain class of linear unbiased estimators. Moreover, the complexity of the
SAOBs is not larger than the complexity of Multilevel Monte Carlo. The
theoretical results are illustrated by numerical experiments with an elliptic
PDE.

###### keywords:

Uncertainty quantification, partial differential equation, Richardson
extrapolation, Monte Carlo, Multilevel Monte Carlo

35R60, 62J05, 62F12, 65N30, 65C05

##  1 Introduction

A common model in the field of uncertainty quantification (UQ) is a partial
differential equation (PDE) equipped with random coefficients, and other
inputs whose uncertainty is modeled by a probability distribution on a
suitable function space. An important building block in UQ is the estimation
of expected values of output quantities of interest linked with such random
PDEs. Monte Carlo (MC) estimators are often infeasible in this situation due
to the high cost per sample. In the last decade multilevel estimators have
been developed to address this problem and provide estimates with a much
smaller computational cost. This is achieved by variance reduction and by
working with a collection of PDE models with different resolutions or
fidelities. Typically, multilevel estimators couple an expensive, high
resolution PDE model with cheap, low resolution PDE models. Examples for
multilevel estimators are multilevel Monte Carlo (MLMC) [7, 8], multifidelity
Monte Carlo (MFMC) [14, 15] and approximate control variates (ACVs) [9]. In
this work we revisit the multilevel best linear unbiased estimator (BLUE)
introduced in [19]. The multilevel BLUE is a linear, unbiased combination of
MC estimators associated with different PDE resolutions. Importantly, the
multilevel BLUE selects the optimal linear combination of samples with respect
to the estimator variance.

The analysis in [19] is independent of the underlying models used for the
multilevel BLUE. Now we specialize the results in [19], and assume that we
work with PDE-based models with random coefficients. For discretized PDEs it
is often the case that the fidelity of a model can be linked with a
discretization parameter, e.g., the mesh size of a finite element space. Thus,
it is natural to assume that the models are parameterized by the mesh size,
where a small parameter value gives accurate results which might however be
expensive to compute. We make the form of the parameterization precise in the
following sections of the paper.

In [19] we introduced a sample allocation optimal estimator termed SAOB which
is a special case of the multilevel BLUE. The SAOB selects so called model
groups and the number of samples for each group such that the estimator
variance is minimal given a fixed computational budget. The SAOB is implicitly
defined in terms of the solution of an optimization problem. This complicates
its analysis, even if we make more assumptions on the models we work with. As
an auxiliary tool we now introduce and study a class of estimators termed RE
estimators. These are based on Richardson extrapolation (RE) [4, 17], and on
the telescoping sum approach in MLMC [7, 8]. It turns out that the classical
MLMC estimator in [7, 8] is a special case of an RE estimator.

RE is a well known technique in numerical analysis. It employs linear
combinations of a family of approximations to improve the accuracy of the
individual approximations within the family. Perhaps the most widely known
application of RE is the Romberg method [18] for numerical quadrature. In
addition, RE has been used for ordinary differential equations [5], stochastic
differential equations (SDEs) [21], and partial differential equations (see
e.g. [1, 2, 16]). By combining RE with the MLMC complexity theory [6, 8] we
provide an upper bound on the computational complexity of the SAOB and show
that its asymptotic complexity is not larger than the complexity of MLMC. We
further show that the coefficients of the SAOB converge to the coefficients in
the linear combination of the RE estimator for an academic toy problem. Hence
we can expect close links between the RE estimator and the SAOB.

The combination of RE and multilevel estimators has been discussed in a few
places in the literature. However, these are typically restricted to bias
errors of a specific form and SDE discretizations. Multilevel
Richardson–Romberg extrapolation has been explored already by Giles in the
pioneering MLMC paper [7] where one level of RE was used. In [12] the authors
combine RE and MLMC for specific discretizations of the Langevin equation.
Lemaire and Pagès [10] introduce a multilevel Richardson–Romberg estimator
termed ML2R for SDE discretizations where the bias error w.r.t. the
discretization parameter $h$ has the form $h^{\gamma^{k}}$ with linearly
growing exponents $\gamma^{k}=k\alpha$. An antithetic extension of the ML2R
estimator is studied in [11]. In our work we study RE estimators to analyze
the complexity of the SAOB. In fact, this idea grew out of numerical
experiments where we observed that the coefficients of the multilevel BLUE
approached RE coefficients. This highlights a novel application of RE in the
study of multilevel estimators and uncertainty quantification.

The main contributions of this paper are as follows: $(i)$ a general upper
bound on the complexity of the SAOB, $(ii)$ a specific upper bound in terms of
RE estimators for parametric model families and $(iii)$ a complete complexity
and variance analysis of the RE estimators.

The remainder of this work is structured as follows. In Section 2 we review
the necessary definitions and results introduced in [19], in particular, the
multilevel BLUE and the SAOB. We further present an asymptotic cost bound for
the SAOB, however, this bound is difficult to verify in practice. In Section 3
we provide upper bounds for the cost of the SAOB which are easier to verify.
These bounds are based on MC, MLMC and RE estimators. In Section 4 we verify
the theoretical results by numerical experiments for a PDE-based quantity of
interest in two space dimensions. In Section 5 we show that the SAOBs converge
to RE estimators for an academic toy problem. In Section 6 we compare the
SAOBs and RE estimators with the ACV estimators. Section 7 provides concluding
remarks.

##  2 Problem formulation

Let $Z$ be a scalar-valued random variable whose expectation $\mathbb{E}[Z]$
we want to estimate. We assume that exact sampling from $Z$ is not possible.
Hence we work with a family of approximations $Z_{1},\dots,Z_{L}$,
$L\in\mathbb{N}$. Typical scenarios we have in mind are finite element based
PDE discretizations where $\ell$ denotes the level of mesh refinement and
$Z_{\ell}(\omega)$ is an output quantity of interest which requires solving
the discretized PDE with random inputs depending on the event $\omega$. In
this case, for a large discretization parameter $\ell$ the approximation of
$Z(\omega)$ by $Z_{\ell}(\omega)$ is accurate yet computationally expensive.

Let us define the expectation and covariance for the random variables
$Z_{1},\dots,Z_{L}$ as follows,

$\displaystyle\mu_{\ell}:=\mathbb{E}[Z_{\ell}],\quad
c_{\ell,j}:=\operatorname{Cov}(Z_{\ell},Z_{j}),\quad\text{for all
}\ell,j\in\\{1,\dots,L\\}.$

We further define the mean vector $\mu:=(\mu_{1},\dots,\mu_{L})^{T}$ and the
model covariance matrix $C:=(c_{\ell,j})_{\ell,j=1}^{L}$. We assume throughout
this paper that the expectations and variances of $Z,Z_{1},\dots,Z_{L}$ exist
and are finite. For every subset $Q\subseteq\\{1,\dots,L\\}$ we define the
principal submatrix of $C$ as

$C_{Q,Q}:=(c_{\ell,j})_{\ell,j\in Q}\in\mathbb{R}^{|Q|\times|Q|}.$

We use a similar notation for vectors $\alpha\in\mathbb{R}^{L}$ with
$\alpha_{Q}:=(\alpha_{\ell})_{\ell\in Q}\in\mathbb{R}^{|Q|}$. We further
introduce the notation $\widehat{\mu}_{\alpha}$ for unbiased estimators of
$\alpha^{T}\mu$, i.e., estimators which satisfy
$\mathbb{E}[\widehat{\mu}_{\alpha}]=\alpha^{T}\mu$. In the special case
$\alpha=e_{\ell}$, that is, $\alpha$ is the $\ell$th unit vector, we write
$\widehat{\mu}_{\ell}$ for an unbiased estimator of
$e_{\ell}^{T}\mu=\mu_{\ell}$. Finally, we use the generic constant $c$ for
estimates, i.e. if it holds

$\phi(\ell)\leq c2^{-2\ell},$

for some function $\phi$, then $c$ is independent of $\ell$. Moreover, the
value of the generic constant $c$ may change from equation to equation.

###  2.1 Multilevel best linear unbiased estimator

Our goal is to construct a variance minimal, linear, unbiased estimator for
$\alpha^{T}\mu$, where the vector $\alpha\in\mathbb{R}^{L}$ is given. Recall
that the estimator $\widehat{\mu}_{\alpha}$ is linear and unbiased if there
exist coefficients $\beta_{1},\dots,\beta_{N}\in\mathbb{R}$ and samples
$Z_{\ell_{1}}(\omega_{j_{1}}),\dots,Z_{\ell_{N}}(\omega_{j_{N}})$ for
$\ell_{1},\dots,\ell_{N}\in\\{1,\dots,L\\}$ such that

$\widehat{\mu}_{\alpha}=\sum_{i=1}^{N}\beta_{i}Z_{\ell_{i}}(\omega_{j_{i}}),\quad\mathbb{E}[\widehat{\mu}_{\alpha}]=\alpha^{T}\mu.$

The events $\omega_{j_{i}}$ determine the correlation structure between the
samples and thus the variance of the estimator. In particular, not all events
are necessarily distinct and the same event might be used with different
output quantities. Throughout this paper we use the framework established in
[19]. Let $K:=2^{L}-1$, and let $S^{1},\dots,S^{K}\in
2^{\\{1,\dots,L\\}}\setminus\\{\emptyset\\}$ denote the non-empty subsets of
the index set $\\{1,\dots,L\\}$. In particular,

$S^{j}\not=S^{k}\quad\text{ for }j\not=k.$

Each collection of indices $S^{k}$ is a so-called model group. For an event
$\omega$ we define the corresponding model group evaluation, the covariance
matrix, the restriction and prolongation matrices, respectively, as follows,

$\displaystyle Z^{k}(\omega)$ $\displaystyle:=(Z_{\ell}(\omega))_{\ell\in
S^{k}}\in\mathbb{R}^{|S^{k}|},\quad$ $\displaystyle C^{k}$
$\displaystyle:=C_{S^{k},S^{k}}=\operatorname{Cov}(Z^{k},Z^{k})\in\mathbb{R}^{|S^{k}|\times|S^{k}|},$
$\displaystyle R^{k}v$ $\displaystyle:=v_{S^{k}}\in\mathbb{R}^{|S^{k}|}\text{
for all }v\in\mathbb{R}^{L},\quad$ $\displaystyle P^{k}$
$\displaystyle:=(R^{k})^{T}\in\mathbb{R}^{L\times|S^{k}|},\quad k=1,\dots,K.$

In the remainder of this paper we assume that each matrix $C^{k}$ defined
above is regular. Next we consider $m_{k}\in\mathbb{N}_{0}$ independent
samples associated with each model group $S^{k}$, denoted by $\omega_{i}^{k}$,
$i=1,\dots,m_{k}$. We further assume that $\omega_{i}^{k}$ and
$\omega_{i}^{j}$ are independent for $k\not=j$. We construct a linear unbiased
estimator using the samples

(2.1)
$(Z^{1}(\omega^{1}_{i}))_{i=1}^{m_{1}},\dots,(Z^{K}(\omega^{K}_{i}))_{i=1}^{m_{K}}.$

Note that the model group $S^{k}$ describes the statistical coupling of models
$Z_{\ell}$ with $\ell\in S^{k}$ and that $m_{k}$ is the number of independent
evaluations of $Z^{k}$. We organize the samples in (2.1) in a block vector and
define a block linear model as follows,

(2.2) $\displaystyle Y$ $\displaystyle=H\mu+\eta,$ $\displaystyle Y$
$\displaystyle=(Y^{k})_{k=1}^{K},\quad$ $\displaystyle H$
$\displaystyle=(H^{k})_{k=1}^{K},\quad$ $\displaystyle\eta$
$\displaystyle=(\eta^{k})_{k=1}^{K},$ $\displaystyle Y^{k}$
$\displaystyle=(Z^{k}(\omega^{k}_{i}))_{i=1}^{m_{k}},\quad$ $\displaystyle
H^{k}$ $\displaystyle=(R^{k})_{i=1}^{m_{k}},\quad$ $\displaystyle\eta^{k}$
$\displaystyle=(Z^{k}(\omega^{k}_{i})-R^{k}\mu)_{i=1}^{m_{k}},$

where $Y$ is the vector of observations, $H$ is the design matrix describing a
linear relationship between the unknown $\mu$ and observations, and $\eta$ is
a noise vector with mean zero. This is a generalized linear model, where the
covariance matrix of the noise $\eta$ is block diagonal due to the assumed
independence structure of the samples,

$\displaystyle\operatorname{Cov}(\eta,\eta)$
$\displaystyle=\operatorname{diag}((\operatorname{Cov}(\eta^{k},\eta^{k}))_{k=1}^{K})=\operatorname{diag}(\operatorname{diag}((\operatorname{Cov}(Z^{k}(\omega^{k}_{i}),Z^{k}(\omega^{k}_{i})))_{i=1}^{m_{k}})_{k=1}^{K})$
$\displaystyle=\operatorname{diag}(((C^{k})_{i=1}^{m_{k}})_{k=1}^{K}).$

Finally, let $m=(m_{1},\dots,m_{K})^{T}$ denote the vector collecting the
number of samples in each model group, and define the matrix and vector

$\Psi(m):=\sum_{k=1}^{K}m_{k}P^{k}(C^{k})^{-1}R^{k}\in\mathbb{R}^{L\times
L},\quad
y(m):=\sum_{k=1}^{K}P^{k}(C^{k})^{-1}\sum_{i=1}^{m_{k}}Z^{k}(\omega_{i}^{k})\in\mathbb{R}^{L}.$

Then, the best linear unbiased estimator (BLUE)
$\widehat{\mu}^{\operatorname{B}}(m)$ for the mean vector $\mu$ in (2.2)
associated with the number of samples $m$ is the solution of the normal
equations

(2.3) $\Psi(m)\widehat{\mu}^{\operatorname{B}}(m)=y(m).$

Under certain assumptions the scalar value
$\widehat{\mu}^{\operatorname{B}}_{\alpha}(m):=\alpha^{T}\widehat{\mu}^{\operatorname{B}}(m)$
is also the BLUE for $\alpha^{T}\mu$ and the estimator variance is
$\operatorname{Var}(\widehat{\mu}^{\operatorname{B}}_{\alpha}(m))=\alpha^{T}\Psi(m)^{-1}\alpha$
(see [19, Theorem 2.7]).

###  2.2 Sample allocation optimal BLUE

Observe that the estimator $\widehat{\mu}^{\operatorname{B}}_{\alpha}(m)$
depends on the number of samples $m_{k}$ for each model group $S^{k}$. Now we
want to select $m$ optimally given a fixed computational budget. We assume
costs $w_{\ell}>0$ to compute a sample of $Z_{\ell}$, $\ell=1,\dots,L$. This
incurs the cost

$W_{k}:=\sum_{\ell\in S^{k}}\operatorname{Cost}(Z_{\ell})=\sum_{\ell\in
S^{k}}w_{\ell}$

to evaluate all models in the group $S^{k}$, $k=1,\dots,K$. For a fixed
computational budget $p>0$ we then solve the following sample allocation
problem:

(2.4)
$\begin{cases}\quad\min_{m\in\mathbb{N}_{0}^{K}}\quad\operatorname{Var}(\widehat{\mu}^{\operatorname{B}}_{\alpha}(m))&=\alpha^{T}\Psi(m)^{-1}\alpha\\\
\quad\text{s.t.}\quad\operatorname{Cost}(\widehat{\mu}^{\operatorname{B}}_{\alpha}(m))&=\sum_{k=1}^{K}m_{k}W_{k}\leq
p,\\\ \quad m_{k}&=0,\quad\text{if }|S^{k}|>q.\end{cases}$

In (2.4), the coupling number $q>0$ defines the maximal number of models that
are evaluated for the same input $\omega$. We then define the estimator termed
$\operatorname{SAOB},q$ as

(2.5)
$\widehat{\mu}^{\operatorname{SAOB,q}}_{\alpha}=\alpha^{T}\widehat{\mu}^{\operatorname{B}}(m^{*}),$

where $m^{*}$ is a minimizer of (2.4). If the coupling number $q=+\infty$, we
drop $q$ in the notation above. It can be proved that the
$\operatorname{SAOB}$ is variance minimal in the class of linear unbiased
estimators with costs bounded by $p$.

######  Theorem 2.1 ([19, Theorem 3.2]).

Let the model covariance matrix $C$ be positive definite. Then any linear
unbiased estimator $\widehat{\mu}_{\alpha}$ that uses the samples in (2.1) for
any $m$ such that $\operatorname{Cost}(\widehat{\mu}_{\alpha})\leq p$
satisfies

$\operatorname{Var}(\widehat{\mu}_{\alpha})\geq\operatorname{Var}(\widehat{\mu}^{\operatorname{SAOB}}_{\alpha}).$

Note that the $\operatorname{SAOB}$ is difficult to analyze since it depends
on a minimizer $m^{*}$ of the problem (2.4), and is constructed implicitly
with the solution of a linear regression problem. To make progress we now
present some general results on the complexity of linear unbiased estimators,
keeping in mind that the $\operatorname{SAOB}$ is a special case.

###  2.3 Asymptotic analysis of linear unbiased estimators

We are interested in linear unbiased estimators of $\alpha^{T}\mu$ that form a
linear combination of correlated MC estimators

(2.6) $\widehat{\mu}_{\alpha}=\sum_{k=1}^{K}\sum_{\ell\in
S^{k}}\beta_{\ell}^{k}\frac{1}{m_{k}}\sum_{i=1}^{m_{k}}Z_{\ell}(\omega_{i}^{k}).$

If the model group $S^{k}$ is not used and thus $m_{k}=0$, we define
$\beta^{k}=0$ and $\beta_{\ell}^{k}/m_{k}=0$. This is equivalent to excluding
the $k$th summand in (2.6). We further define $\beta^{k}_{\ell}=0$ for
$\ell\not\in S^{k}$. Combining the definitions in [19, Equ. (2.11)], [19, Equ.
(2.7)], and rearranging, we see that the coefficients $\beta^{k}$ in (2.6)
associated with $\widehat{\mu}_{\alpha}^{\text{SAOB}}$ can be written as

(2.7) $\beta^{k}=m_{k}P^{k}(C^{k})^{-1}R^{k}\Psi(m)^{-1}\alpha$

provided that $\Psi(m)$ is regular. Recall that $Z_{\ell}(\omega_{i}^{k})$ and
$Z_{\ell}(\omega_{j}^{n})$ are statistically independent if $i\neq j$ or
$k\neq n$. Thus the estimator variance is equal to

(2.8)
$\operatorname{Var}(\widehat{\mu}_{\alpha})=\sum_{k=1}^{K}\frac{1}{m_{k}}\operatorname{Var}\left(\sum_{\ell\in
S^{k}}\beta^{k}_{\ell}Z_{\ell}\right)=\sum_{k=1}^{K}\frac{1}{m_{k}}(\beta^{k})^{T}C\beta^{k}.$

Now, if $(\beta^{k})^{T}C\beta^{k}$ is asymptotically small, then only a small
number of samples $m_{k}$ is required to achieve a small variance
contribution. Therefore we want to obtain an expression for the optimal sample
allocation $m$ by minimizing (2.8) given a budget constraint

(2.9)
$\begin{cases}\quad\min_{m\in\mathbb{R}^{K}}J(m)&=\sum_{k=1}^{K}\frac{1}{m_{k}}(\beta^{k})^{T}C\beta^{k}\\\
\quad\text{s.t.}\sum_{k=1}^{K}m_{k}W_{k}&\leq p,\\\ \quad m_{k}&\geq
0,\quad\text{for }k=1,\dots,K.\end{cases}$

We now follow the classical MLMC approach [8, Section 1.3] to compute the
unique minimizer of (2.9). Note that we replace
$\operatorname{Var}(Z_{\ell}-Z_{\ell-1})$ by the more general term
$(\beta^{k})^{T}C\beta^{k}$.

######  Lemma 2.2 (Optimal sample allocation).

Let $C$ be positive definite, let $\alpha\not=0$ and the coefficients
$\beta^{k}_{\ell}$ such that the bias constraint is satisfied, i.e.
$\alpha=\sum_{k=1}^{K}\beta^{k}$. Then there exists a unique minimizer $m^{*}$
of (2.9) of the form

(2.10)
$m^{*}_{k}=\frac{p}{\sum_{k=1}^{K}(\left[(\beta^{k})^{T}C\beta^{k}\right]W_{k})^{1/2}}\left(\left[(\beta^{k})^{T}C\beta^{k}\right]/W_{k}\right)^{1/2}$

with associated variance

(2.11)
$J(m^{*})=\frac{1}{p}\left(\sum_{k=1}^{K}\left(\left[(\beta^{k})^{T}C\beta^{k}\right]W_{k}\right)^{1/2}\right)^{2}.$

In particular, the cost $p$ to achieve the variance $J(m^{*})=\varepsilon^{2}$
is equal to

$p=\varepsilon^{-2}\left(\sum_{k=1}^{K}\left(\left[(\beta^{k})^{T}C\beta^{k}\right]W_{k}\right)^{1/2}\right)^{2}.$

######  Proof 2.3.

The proof is provided in Appendix A.

Note that we assume real-valued numbers in (2.10) which is not practical.
However, rounding $m^{*}_{k}$ in (2.10) for all $K=2^{L}-1$ possible values of
$k$ may increase the total cost significantly, since $2^{L-1}$ model groups
contain the finest model $Z_{L}$. Fortunately, the result [19, Theorem 3.6]
tells us that we may choose $\beta^{k}\not=0$ for at most $L$ different
indices $k\in\\{1,\dots,K\\}$ without increasing the variance. Without loss of
generality this allows us to set $\beta^{L+1}=\beta^{L+2}=\dots=\beta^{K}=0$
by suitably renumbering $S^{1},\dots,S^{K}$. Since the estimator
$\widehat{\mu}_{\alpha}$ is an unbiased estimator for $\alpha^{T}\mu$, the
$\beta^{k}$ then satisfy

(2.12)
$\alpha=\sum_{k=1}^{L}\beta^{k},\quad\beta^{k}_{\ell}=0,\quad\ell\not\in
S^{k}.$

Let us now comment on the bias which is determined by the vector $\alpha$.
Observe that in practise it is not advisable to fix $\alpha$ independently of
the finest level $L$. For example, for $L=1$ we may choose $\alpha=1$,
however, for $L=2$ the choice $\alpha=(0,1)^{T}$ is suitable provided that
$Z_{2}$ is a better approximation of $Z$ compared to $Z_{1}$. Hence we
introduce a sequence of bias vectors $(\alpha^{L})_{L=1}^{\infty}$ that in
turn defines the sequence of estimators
$(\widehat{\mu}_{\alpha^{L}})_{L=1}^{\infty}$. The coefficients $\beta^{k}$
and model groups $S^{k}$ clearly depend on $L$, however for ease of notation
we drop this dependence. Furthermore, for $(\alpha^{L})^{T}\mu$ we assume that
$\alpha^{L}\in\mathbb{R}^{L}$ and
$\mu=(\mathbb{E}[Z_{1}],\dots,\mathbb{E}[Z_{L}])^{T}\in\mathbb{R}^{L}$ where
we again drop the dependence of $\mu$ on $L$. We now estimate the asymptotic
complexity of the sequence of estimators
$(\widehat{\mu}_{\alpha^{L}})_{L=1}^{\infty}$.

######  Theorem 2.4 (Asymptotic cost of linear unbiased estimators).

Assume that there exist positive constants $\gamma^{\operatorname{Bias}}$,
$\gamma^{\operatorname{Var}}$ and $\gamma^{\operatorname{Cost}}$, such that
the following statements hold for all $L\in\mathbb{N}$:

(M0) $\displaystyle\alpha^{L}$
$\displaystyle=\sum_{k=1}^{L}\beta^{k},\quad\quad\beta^{k}_{j}=0\text{ if
}j\not\in S^{k}$ $\displaystyle\quad\text{ for all }k\leq L,$ (M1)
$\displaystyle|(\alpha^{L})^{T}\mu-\mathbb{E}[Z]|$ $\displaystyle\leq
c2^{-L\gamma^{\operatorname{Bias}}},$ (M2)
$\displaystyle(\beta^{k})^{T}C\beta^{k}$ $\displaystyle\leq
c2^{-k\gamma^{\operatorname{Var}}}$ $\displaystyle\quad\text{ for all }k\leq
L,$ (M3) $\displaystyle W_{k}$ $\displaystyle\leq
c2^{k\gamma^{\operatorname{Cost}}}$ $\displaystyle\quad\text{ for all }k\leq
L,$

where the constant $c>0$ is independent of $L$. Then there exists a level
$L\in\mathbb{N}$, model groups $S^{1},\dots,S^{L}$ and numbers of samples
$m_{1},\dots,m_{L}$ to achieve
$\mathbb{E}[|\widehat{\mu}_{\alpha^{L}}-\mathbb{E}[Z]|^{2}]\leq\varepsilon^{2}$
with a cost bounded by

(2.13) $\operatorname{Cost}(\widehat{\mu}_{\alpha^{L}})\leq
c\varepsilon^{-\gamma^{\operatorname{Cost}}/\gamma^{\operatorname{Bias}}}+c\begin{cases}\varepsilon^{-2},\quad&\text{if
}\gamma^{\operatorname{Cost}}<\gamma^{\operatorname{Var}},\\\
\varepsilon^{-2}(\log\varepsilon)^{2},\quad&\text{if
}\gamma^{\operatorname{Cost}}=\gamma^{\operatorname{Var}},\\\
\varepsilon^{-2-\frac{\gamma^{\operatorname{Cost}}-\gamma^{\operatorname{Var}}}{\gamma^{\operatorname{Bias}}}},\quad&\text{if
}\gamma^{\operatorname{Cost}}>\gamma^{\operatorname{Var}}.\end{cases}$

######  Proof 2.5.

The proof is analogous to the proof of [6, Theorem 1]. Note that we do not use
the assumption
$\gamma^{\operatorname{Bias}}\geq\min\\{\gamma^{\operatorname{Var}},\gamma^{\operatorname{Cost}}\\}/2$.
Thus we do not explicitly bound the term
$\varepsilon^{-\gamma^{\operatorname{Cost}}/\gamma^{\operatorname{Bias}}}$
accounting for ceiling the number of samples.

Recall that the coefficients $\beta^{k}$ in (2.7) and the $L$ model groups
$S^{k}$ in (2.6) of the $\operatorname{SAOB}$ are chosen to minimize the
variance and are only given implicitly. Therefore, we typically cannot verify
the assumptions (M2) and (M3) for the $\operatorname{SAOB}$. Furthermore, it
is not clear how to choose the sequence of bias vectors
$(\alpha^{L})_{L=1}^{\infty}$ such that (M1) holds. Hence Theorem 2.4 is of
limited practical use to determine the asymptotic cost of the
$\operatorname{SAOB}$.

##  3 Upper asymptotic cost bounds

We now present an alternative, constructive approach for the complexity
analysis of the $\operatorname{SAOB}$. A straightforward consequence of
Theorem 2.1 is that the SAOB has the optimal complexity in the class of linear
unbiased estimators.

######  Theorem 3.1 (Complexity bound for $\operatorname{SAOB}$).

Let the model covariance matrix $C$ be positive definite and let
$\widehat{\mu}_{\alpha}$ be a linear unbiased estimator that uses the samples
in (2.1). If $\widehat{\mu}_{\alpha}$ estimates $\mathbb{E}[Z]$ with a mean
square error (MSE) bounded by $\varepsilon^{2}$ and with costs bounded by
$\phi(\varepsilon)$, that is,

$\displaystyle\mathbb{E}[|\widehat{\mu}_{\alpha}-\mathbb{E}[Z]|^{2}]\leq\varepsilon^{2}\quad\text{with}\quad\operatorname{Cost}(\widehat{\mu}_{\alpha})\leq\phi(\varepsilon),$

then the estimator $\widehat{\mu}^{\operatorname{SAOB}}_{\alpha}$ achieves the
same MSE with the same or smaller costs

$\displaystyle\mathbb{E}[|\widehat{\mu}^{\operatorname{SAOB}}_{\alpha}-\mathbb{E}[Z]|^{2}]\leq\varepsilon^{2}\quad\text{with}\quad\operatorname{Cost}(\widehat{\mu}^{\operatorname{SAOB}}_{\alpha})\leq\operatorname{Cost}(\widehat{\mu}_{\alpha})\leq\phi(\varepsilon).$

In particular, we may choose $p=\operatorname{Cost}(\widehat{\mu}_{\alpha})$
in (2.4).

######  Proof 3.2.

We use a bias variance decomposition and the bound on the MSE to show

$\mathbb{E}[(\widehat{\mu}_{\alpha}-\mathbb{E}[Z])^{2}]=(\alpha^{T}\mu-\mathbb{E}[Z])^{2}+\operatorname{Var}(\widehat{\mu}_{\alpha})\leq\varepsilon^{2}.$

Now, the estimator $\widehat{\mu}^{\operatorname{SAOB}}_{\alpha}$ has the same
bias as $\widehat{\mu}_{\alpha}$. Therefore we only compare the variance of
$\widehat{\mu}^{\operatorname{SAOB}}_{\alpha}$ and $\widehat{\mu}_{\alpha}$,
respectively, as unbiased estimators of $\alpha^{T}\mu$. If we choose
$p=\operatorname{Cost}(\widehat{\mu}_{\alpha})$ in Theorem 2.1, then – by
construction – the estimator $\widehat{\mu}^{\operatorname{SAOB}}_{\alpha}$
has equal or smaller variance than $\widehat{\mu}_{\alpha}$.

######  Remark 3.3.

Similarly to the proof of [19, Theorem 3.2] it can be shown that any linear
unbiased estimator $\widehat{\mu}_{\alpha}$ that uses the samples in (2.1) and
couples at most $q$ models using the same event $\omega$, that is $m_{k}=0$ if
$|S^{k}|>q$, satisfies

$\operatorname{Var}(\widehat{\mu}_{\alpha})\geq\operatorname{Var}(\widehat{\mu}^{\operatorname{SAOB,q}}_{\alpha}).$

Then we again choose $p=\operatorname{Cost}(\widehat{\mu}_{\alpha})$ in (2.4)
and are able to conclude that $\operatorname{SAOB},q$ is complexity optimal in
the class of linear unbiased estimators that couple at most $q$ models.

The complexity bound in Theorem 3.1 is formulated as a comparison principle,
and does not give an explicit expression for $\phi(\varepsilon)$. However,
explicit expressions may be obtained by the following procedure:

  * 1.

Let $(\widehat{\mu}_{\alpha^{L}})_{L=1}^{\infty}$ be a sequence of linear
unbiased estimators.

  * 2.

Bound the asymptotic complexity of
$(\widehat{\mu}_{\alpha^{L}})_{L=1}^{\infty}$ in terms of an explicit
expression $\phi(\varepsilon)$.

  * 3.

Verify the assumptions of Theorem 3.1 for all estimators
$(\widehat{\mu}_{\alpha^{L}})_{L=1}^{\infty}$. Then the complexity of the
$\operatorname{SAOB}$ is bounded by $\phi(\varepsilon)$.

Of course, this approach gives only an upper cost bound for the
$\operatorname{SAOB}$ that is not necessarily sharp, especially if the
estimators $(\widehat{\mu}_{\alpha^{L}})_{L=1}^{\infty}$ are chosen poorly. In
the following we describe a construction that we found useful. First, we
choose the vector $\alpha^{L}$ such that the bias is small. This can often be
achieved by Richardson extrapolation. Next, since the estimator
$\widehat{\mu}_{\alpha^{L}}$ has to be linear and unbiased, we use Lemma 2.2
as a guideline. Observe that for a fixed vector $\alpha^{L}$ the coefficients
of the $\operatorname{SAOB}$ satisfy

(3.1)
$\begin{cases}\quad\min_{\beta^{1},\dots,\beta^{K}\in\mathbb{R}^{L}}J(\beta)&=\sum_{k=1}^{K}\left(\left[(\beta^{k})^{T}C\beta^{k}\right]W_{k}\right)^{1/2}\\\
\quad\text{s.t.}\quad&\alpha^{L}=\sum_{k=1}^{K}\beta^{k},\\\
\quad&\beta^{k}_{\ell}=0\quad\text{if }\ell\not\in S^{k},\end{cases}$

where $\varepsilon^{-2}J^{2}$ is the cost of $\operatorname{SAOB}$ to achieve
a variance of $\varepsilon^{2}$. Note that this does not include the cost for
ceiling the number of samples. We may then choose suboptimal coefficients
$\beta^{k}$ which increases the variance but gives an asymptotic expression or
upper bound for $J^{2}$. The goal of this section is to demonstrate this proof
strategy for different choices of $\beta^{k}$.

  * •

The MC estimator has typically a high asymptotic cost and thus gives an upper
bound that is not sharp. To demonstrate the general idea we write down the MC
complexity in Section 3.1.

  * •

An improvement is possible by making $(\beta^{k})^{T}C\beta^{k}$ small such
that its contribution to $J$ is small while ensuring that the bias constraint
is still satisfied. The MLMC estimator achieves this by using a telescoping
sum and coefficients $\beta^{k}$ such that
$(\beta^{k})^{T}C\beta^{k}=\operatorname{Var}(Z_{\ell}-Z_{\ell-1})$. The
variance of the difference is often asymptotically small with a known rate. We
summarize this result in Section 3.2.

  * •

We further improve the asymptotic complexity of the MLMC estimator by using RE
to obtain even smaller asymptotic expressions for $(\beta^{k})^{T}C\beta^{k}$.
We formally define the resulting RE estimators in Section 3.3.

###  3.1 Monte Carlo estimator

The standard MC estimator uses $\alpha^{L}=e_{L}$, $S^{1}=\\{L\\}$,
$\beta^{1}=e_{L}$ and $\beta^{k}=0$ for all other $k$. The estimator reads

$\widehat{\mu}^{\operatorname{MC}}_{L}=\frac{1}{m_{1}}\sum_{i=1}^{m_{1}}Z_{L}(\omega^{1}_{i}).$

We have the following well known result for the complexity of the MC
estimator.

######  Corollary 3.4 (MC asymptotic cost, [6, Subsection 2.1]).

Assume that there exist positive constants $\gamma^{\operatorname{Bias}}$ and
$\gamma^{\operatorname{Cost}}$ such that the following statements hold for all
$\ell\in\mathbb{N}$,

(MC, M1) $\displaystyle|\mathbb{E}[Z_{\ell}]-\mathbb{E}[Z]|$
$\displaystyle\leq c2^{-\ell\gamma^{\operatorname{Bias}}},$ (MC, M2)
$\displaystyle\operatorname{Var}(Z_{\ell})$ $\displaystyle\leq c,$ (MC, M3)
$\displaystyle\operatorname{Cost}(Z_{\ell})$ $\displaystyle\leq
c2^{\ell\gamma^{\operatorname{Cost}}},$

where the constant $c>0$ is independent of $\ell$. Then there exists a final
level $L>0$ and number of samples $m_{1}$ to achieve a MSE of
$\mathbb{E}[|\widehat{\mu}^{\operatorname{MC}}_{L}-\mathbb{E}[Z]|^{2}]\leq\varepsilon^{2}$
with costs bounded by

(3.2) $\operatorname{Cost}(\widehat{\mu}^{\operatorname{MC}}_{L})\leq
c\varepsilon^{-2-\gamma^{\operatorname{Cost}}/\gamma^{\operatorname{Bias}}}.$

Using Theorem 3.1 we may bound the complexity of the SAOB by the complexity of
MC.

######  Corollary 3.5 (SAOB cost upper bounded by MC).

Let the assumptions of Corollary 3.4 be true. Then the
$\widehat{\mu}^{\operatorname{SAOB}}_{L}$ and
$\widehat{\mu}^{\operatorname{SAOB,q}}_{L}$ with $q\geq 1$ achieve a MSE of
$\varepsilon^{2}$ with a cost not larger than the cost of the MC estimator in
(3.2).

###  3.2 Multilevel Monte Carlo estimator

The MLMC estimator [7, 8] uses the vector $\alpha^{L}=e_{L}$ together with

$\displaystyle S^{1}$ $\displaystyle=\\{1\\},$
$\displaystyle\quad\beta^{1}=e_{1},$ $\displaystyle S^{k}$
$\displaystyle=\\{k-1,k\\},$ $\displaystyle\quad\beta^{k}=e_{k}-e_{k-1},$
$\displaystyle\quad k=2,\dots,L.$

The estimator reads

$\widehat{\mu}^{\operatorname{MLMC}}_{L}=\sum_{k=2}^{L}\frac{1}{m_{k}}\sum_{i=1}^{m_{k}}(Z_{k}(\omega^{k}_{i})-Z_{k-1}(\omega^{k}_{i}))+\frac{1}{m_{1}}\sum_{i=1}^{m_{1}}Z_{1}(\omega^{1}_{i}).$

We recall the well known MLMC complexity theorem.

######  Corollary 3.6 (MLMC asymptotic cost, [6, Theorem 1]).

Assume that there exist positive constants
$\gamma^{\operatorname{Bias}},\gamma^{\operatorname{Var}}$ and
$\gamma^{\operatorname{Cost}}$ such that the following statements hold for all
$\ell\in\mathbb{N}$,

(MLMC, M1) $\displaystyle|\mathbb{E}[Z_{\ell}]-\mathbb{E}[Z]|\leq
c2^{-\ell\gamma^{\operatorname{Bias}}},$ (MLMC, M2)
$\displaystyle\operatorname{Var}(Z_{\ell}-Z_{\ell-1})\leq
c2^{-\ell\gamma^{\operatorname{Var}}},$ (MLMC, M3)
$\displaystyle\operatorname{Cost}(Z_{\ell})\leq
c2^{\ell\gamma^{\operatorname{Cost}}}.$

where the constant $c>0$ is independent of $\ell$. Then there exists a final
level $L>0$ and numbers of samples $m_{1},\dots,m_{L}$ to achieve a MSE of
$\mathbb{E}[|\widehat{\mu}^{\operatorname{MLMC}}_{L}-\mathbb{E}[Z]|^{2}]\leq\varepsilon^{2}$
with costs bounded by

(3.3) $\operatorname{Cost}(\widehat{\mu}^{\operatorname{MLMC}}_{L})\leq
c\varepsilon^{-\gamma^{\operatorname{Cost}}/\gamma^{\operatorname{Bias}}}+c\begin{cases}\varepsilon^{-2},\quad&\text{if
}\gamma^{\operatorname{Cost}}<\gamma^{\operatorname{Var}},\\\
\varepsilon^{-2}(\log\varepsilon)^{2},\quad&\text{if
}\gamma^{\operatorname{Cost}}=\gamma^{\operatorname{Var}},\\\
\varepsilon^{-2-\frac{\gamma^{\operatorname{Cost}}-\gamma^{\operatorname{Var}}}{\gamma^{\operatorname{Bias}}}},\quad&\text{if
}\gamma^{\operatorname{Cost}}>\gamma^{\operatorname{Var}}.\end{cases}$

Note that the complexity of the MLMC estimator is in general smaller compared
to the complexity of the MC estimator. However, (MLMC, M2) is a stronger
assumption compared to (MC, M2). Fortunately, this stronger assumption is
often satisfied if the model $Z_{\ell}$ is derived from the discretization of
an ODE or PDE and thus satisfies the error estimate

$\operatorname{Var}(Z_{\ell}-Z_{\ell-1})\leq\mathbb{E}[(Z_{\ell}-Z_{\ell-1})^{2}]\leq
c(\mathbb{E}[(Z_{\ell}-Z)^{2}]+\mathbb{E}[(Z_{\ell-1}-Z)^{2}])\leq
c2^{-\ell\gamma^{\operatorname{Var}}},$

where $\gamma^{\operatorname{Var}}>0$. We generalize this idea in the next
subsection, where we also obtain rates for $\gamma^{\operatorname{Bias}}$ and
$\gamma^{\operatorname{Var}}$. Using again Theorem 3.1 gives a complexity
bound for the SAOB.

######  Corollary 3.7 (SAOB cost upper bounded by MLMC).

Let the assumptions of Corollary 3.6 be true. Then the
$\widehat{\mu}^{\operatorname{SAOB}}_{L}$ and
$\widehat{\mu}^{\operatorname{SAOB,q}}_{L}$ with $q\geq 2$ achieve a MSE of
$\varepsilon^{2}$ with a cost not larger than the cost of the MLMC estimator
in (3.3).

###  3.3 Richardson extrapolation estimator

Richardson extrapolation (RE) [1, 10, 17, 18] is a well known technique to
improve the accuracy of a collection of parametric models . The idea is to
linearly combine models with a low accuracy, and to increase the convergence
rate to the truth $Z$. We remark that RE was already used in [7, 10] to
increase the bias rate $\gamma^{\operatorname{Bias}}$. The analysis in this
section is similar to the analysis of MLMC. We replace the assumptions (MLMC,
M1), (MLMC, M2) and (MLMC, M3) with the following assumptions.

[Pathwise Asymptotic Expansion]

  * $(i)$

There exists $q\in\mathbb{N}$, rates $0=\gamma^{1}<\dots<\gamma^{q}$ and
random variables $c_{2},\dots,c_{q}$ with bounded second moment such that it
holds

(3.4)
$Z_{\ell}(\omega)=Z(\omega)+\sum_{j=2}^{q-1}c_{j}(\omega)2^{-\ell\gamma^{j}}+\mathcal{O}(2^{-\ell\gamma^{q}}),\quad\text{as
}\ell\rightarrow+\infty.$

  * $(ii)$

The costs for an evaluation of $Z_{\ell}$ are bounded by

(3.5) $\operatorname{Cost}(Z_{\ell})\leq
c2^{\ell\gamma^{\operatorname{Cost}}}.$

The $\mathcal{O}$-notation in (3.4) has to be understood in the $L^{2}$-sense
for $\ell\rightarrow+\infty$, that is,

$\mathbb{E}\left[\left(Z_{\ell}(\omega)-Z(\omega)-\sum_{j=2}^{q-1}c_{j}(\omega)2^{-\ell\gamma^{j}}\right)^{2}\right]\leq
c2^{-2\ell\gamma^{q}}\quad\text{for all }\ell\text{ large enough}.$

Note that (3.5) is in fact (MLMC, M3). For a fixed event $\omega$ the
assumption in (3.4) is classical in the theory of Richardson extrapolation. It
has been adapted to our setting where we deal with realizations of random
variables in contrast to deterministic quantities.

The idea of RE is to linearly combine the models $Z_{\ell}$ to remove the term
$\sum_{j=2}^{q}c_{j}(\omega)2^{-\ell\gamma^{j}}$ in (3.4). To this end we
define the RE coefficients

(3.6) $v^{k,q}:=\begin{cases}0,&\quad\text{if }k=0,\\\ e_{1},&\quad\text{if
}k=1,\\\
(2^{\gamma^{k}}Dv^{k-1,q}-v^{k-1,q})/(2^{\gamma^{k}}-1),&\quad\text{if
}1<k<q,\\\ Dv^{k-1,q},&\quad\text{if }k\geq q,\end{cases}$

where $D$ is a down shift matrix

(3.7) $D:=\begin{pmatrix}0&0\\\
I_{L-1,L-1}&0\end{pmatrix}\in\mathbb{R}^{L\times L}.$

The idea of the RE estimator is analogous to MLMC in the sense that we also
use a telescoping sum. Let $\alpha^{L}=v^{L,q}$ and

(3.8) $\displaystyle S^{1}$ $\displaystyle=\\{1\\},$
$\displaystyle\;\beta^{1}=v^{1,q},$ $\displaystyle S^{k}$
$\displaystyle=\\{\max\\{k-q+1,1\\},\max\\{k-q+2,1\\},\dots,k\\},$
$\displaystyle\;\beta^{k}=v^{k,q}-v^{k-1,q},$ $\displaystyle\quad
k=2,\dots,L.$

The $\operatorname{RE,q}$ estimator is then defined as

$\widehat{\mu}^{\operatorname{RE,q}}_{v^{L,q}}=\sum_{k=1}^{L}\sum_{\ell\in
S^{k}}(v^{k,q}_{\ell}-v^{k-1,q}_{\ell})\frac{1}{m_{k}}\sum_{i=1}^{m_{k}}Z_{\ell}(\omega^{k}_{i}).$

A close inspection shows that for $q=2$ this estimator is actually the MLMC
estimator. The $\operatorname{RE},q$ couples at most $q$ models such that
$|S^{k}|\leq q$ for $k=1,\dots,L$. We have the following result.

######  Lemma 3.8.

Let Section 3.3 be true. Then the following estimates hold

(RE, M1) $\displaystyle|(v^{k,q})^{T}\mu-\mathbb{E}[Z]|$ $\displaystyle\leq
c2^{-k\gamma^{q}},$ (RE, M2)
$\displaystyle(v^{k,q}-v^{k-1,q})^{T}C(v^{k,q}-v^{k-1,q})$ $\displaystyle\leq
c2^{-2k\gamma^{q}},$ (RE, M3) $\displaystyle W_{k}$ $\displaystyle\leq
c2^{k\gamma^{\operatorname{Cost}}}.$

######  Proof 3.9.

See Appendix B.

We now state complexity rates for the RE estimator which can be proved using
Lemma 3.8 and standard arguments.

######  Corollary 3.10 (RE asymptotic cost).

Let Section 3.3 or (RE, M1), (RE, M2) and (RE, M3) be true. Then there exists
a final level $L>0$ and numbers of samples $m_{1},\dots,m_{L}$ to achieve a
MSE of
$\mathbb{E}[|\widehat{\mu}^{\operatorname{RE,q}}_{v^{L,q}}-\mathbb{E}[Z]|^{2}]\leq\varepsilon^{2}$
with costs bounded by

(3.9) $\operatorname{Cost}(\widehat{\mu}^{\operatorname{RE,q}}_{v^{L,q}})\leq
c\varepsilon^{-\gamma^{\operatorname{Cost}}/\gamma^{q}}+c\begin{cases}\varepsilon^{-2},\quad&\text{if
}\gamma^{\operatorname{Cost}}<2\gamma^{q},\\\
\varepsilon^{-2}(\log\varepsilon)^{2},\quad&\text{if
}\gamma^{\operatorname{Cost}}=2\gamma^{q},\\\
\varepsilon^{-2-\frac{\gamma^{\operatorname{Cost}}-2\gamma^{q}}{\gamma^{q}}},\quad&\text{if
}\gamma^{\operatorname{Cost}}>2\gamma^{q}.\end{cases}$

Observe that the complexity of the RE estimator is improved compared to the
MLMC estimator since we can leverage the additional smoothness of the models
$Z_{\ell}$ in assumption (3.4). We obtain an increased bias rate of
$\gamma^{q}$ instead of $\gamma^{2}$ and an increased variance rate of
$2\gamma^{q}$ instead of $2\gamma^{2}$. Using Theorem 3.1 we obtain another
complexity bound for the SAOB.

######  Corollary 3.11 (SAOB cost upper bounded by RE).

Let the assumptions of Corollary 3.10 be true. Then the
$\widehat{\mu}^{\operatorname{SAOB}}_{v^{L,q}}$ and
$\widehat{\mu}^{\operatorname{SAOB,s}}_{v^{L,q}}$ with $s\geq q$ achieve a MSE
of $\varepsilon^{2}$ with a cost not larger than the cost of the RE estimator
in (3.9).

###  3.4 Weighted Richardson extrapolation estimator

Finally, we generalize the RE estimators to obtain different bias and variance
rates. Note that for $s\in\mathbb{N}$ the definition (3.6) shows that the
vectors

$\\{v^{1,s}-v^{0,s},v^{2,s}-v^{1,s},\dots,v^{L,s}-v^{L-1,s}\\}$

form a basis of $\mathbb{R}^{L}$. Thus there exist weights $a_{1},\dots,a_{L}$
such that for $t\in\mathbb{N}$ it holds

(3.10) $v^{L,t}=\sum_{k=1}^{L}a_{k}(v^{k,s}-v^{k-1,s}).$

Now, we define the weighted RE estimator

$\widehat{\mu}^{\operatorname{RE,s}}_{v^{L,t}}=\sum_{k=1}^{L}a_{k}\sum_{\ell\in
S^{k}}(v^{k,s}_{\ell}-v^{k-1,s}_{\ell})\frac{1}{m_{\ell}}\sum_{i=1}^{m_{\ell}}Z_{\ell}(\omega^{k}_{i}),$

which is a linear unbiased estimator for $(v^{L,t})^{T}\mu$. A result similar
to Lemma 3.8 holds.

######  Lemma 3.12.

Let Section 3.3 be true and $t,s\leq q$. Then the following estimates hold:

(WRE, M1) $\displaystyle|(v^{\ell,t})^{T}\mu-\mathbb{E}[Z]|$
$\displaystyle\leq c2^{-\ell\gamma^{t}},$ (WRE, M2)
$\displaystyle[a_{k}(v^{k,s}-v^{k-1,s})]^{T}C[a_{k}(v^{k,s}-v^{k-1,s})]$
$\displaystyle\leq ca^{2}_{k}2^{-2k\gamma^{s}},$ (WRE, M3) $\displaystyle
W_{k}$ $\displaystyle\leq c2^{k\gamma^{\operatorname{Cost}}}.$

Furthermore, if $s\leq t$ then $|a_{k}|\leq c$ with a constant $c$ independent
of $k$.

######  Proof 3.13.

See Appendix C.

Note that in numerical experiments we observed that the property $|a_{k}|\leq
c$ with a constant $c$ independent of $k$ also holds for $s>t$, however, we do
not see a way to prove it. We directly obtain the analogous result of
Corollary 3.10.

######  Corollary 3.14 (Weighted RE asymptotic cost).

Let Section 3.3 or (WRE, M1), (WRE, M2) and (WRE, M3) be true with $s\leq
t\leq q$. Then there exists a final level $L>0$ and numbers of samples
$m_{1},\dots,m_{L}$ to achieve a MSE of
$\mathbb{E}[|\widehat{\mu}^{\operatorname{RE,s}}_{v^{L,t}}-\mathbb{E}[Z]|^{2}]\leq\varepsilon^{2}$
with costs bounded by

(3.11) $\operatorname{Cost}(\widehat{\mu}^{\operatorname{RE,s}}_{v^{L,t}})\leq
c\varepsilon^{-\gamma^{\operatorname{Cost}}/\gamma^{t}}+c\begin{cases}\varepsilon^{-2},\quad&\text{if
}\gamma^{\operatorname{Cost}}<2\gamma^{s},\\\
\varepsilon^{-2}(\log\varepsilon)^{2},\quad&\text{if
}\gamma^{\operatorname{Cost}}=2\gamma^{s},\\\
\varepsilon^{-2-\frac{\gamma^{\operatorname{Cost}}-2\gamma^{s}}{\gamma^{t}}},\quad&\text{if
}\gamma^{\operatorname{Cost}}>2\gamma^{s}.\end{cases}$

The smallest asymptotic cost is achieved if $s=t=q$ with maximal $q$ such that
Section 3.3 holds. In this case Corollary 3.14 is identical to Corollary 3.10.
The complexity bound for the SAOB follows again from Theorem 3.1.

######  Corollary 3.15 (SAOB cost upper bounded by weighted RE).

Let the assumptions of Corollary 3.14 be true. Then the
$\widehat{\mu}^{\operatorname{SAOB}}_{v^{L,t}}$ and
$\widehat{\mu}^{\operatorname{SAOB,s^{\prime}}}_{v^{L,t}}$ with
$s^{\prime}\geq s$ achieve a MSE of $\varepsilon^{2}$ with a cost not larger
than the cost of the weighted RE estimator in (3.11).

##  4 Numerical experiments with a PDE in 2D space

In this section we verify the main results of this paper numerically. We
revisit the example in [19, Sec. 6.3]. Recall that the random variable $Z$ is
the domain average with realizations

(4.1) $Z(\omega):=\frac{1}{|D_{obs}|}\int_{D_{obs}}y(x,\omega)dx,$

where $D_{obs}:=(3/4,7/8)\times(7/8,1)$ is a subset of the unit square domain
$D:=(0,1)^{2}$. The function $y$ is a weak solution of the elliptic boundary
value problem whose pathwise, strong form reads

(4.2) $\displaystyle-\operatorname{div}(a(x,\omega)\nabla y(x,\omega))$
$\displaystyle=1,\quad x\in D,$ $\displaystyle y(x,\omega)$
$\displaystyle=0,\quad x\in\partial D$

for almost all $\omega$. The diffusion coefficient is
$a(x,\omega):=\exp(b(x,\omega))$, where $b$ is a mean zero Gaussian random
field with Whittle–Matérn covariance function [20] with smoothness parameter
$\nu=3/2$, variance $\sigma^{2}=1$, and correlation length $\rho=0.5$. To
discretize $Z$ we use a standard piecewise linear finite element (FE)
discretization of (4.2) with $L=6$ different levels of accuracy. We obtain
$Z_{\ell}$ by uniform refinement of the previous mesh, starting with a regular
grid for $Z_{1}$. We define the bias as follows:

(4.3) $\operatorname{Bias}(Z_{\ell}):=|\mathbb{E}[Z_{\ell}]-\mathbb{E}[Z]|.$

The discretization data is summarized in Table 1. We use $10^{5}$ pilot
samples to estimate the $\operatorname{Bias}$ in (4.3), the model covariance
matrix $C$, and the average cost $w_{\ell}$ to compute a sample of $Z_{\ell}$.
The cost for the pilot samples is not included in the subsequent analysis. The
fourth column in Table 1 suggests the cost increase rate
$\gamma^{\operatorname{Cost}}=2$ which corresponds to a four-fold cost
increase as expected in 2D space.

Model | #Nodes | $w_{\ell}$ | $w_{\ell}/w_{\ell-1}$ | $\operatorname{Var}(Z_{\ell})$ | $\operatorname{Bias}(Z_{\ell})$
---|---|---|---|---|---
$Z_{1}$ | 81 | 0.0019s | – | $1.52\cdot 10^{-2}$ | $3.80\cdot 10^{-3}$
$Z_{2}$ | 289 | 0.0025s | 1.33 | $2.35\cdot 10^{-2}$ | $1.02\cdot 10^{-3}$
$Z_{3}$ | 1089 | 0.0050s | 2.00 | $2.62\cdot 10^{-2}$ | $2.58\cdot 10^{-4}$
$Z_{4}$ | 4225 | 0.0185s | 3.70 | $2.69\cdot 10^{-2}$ | $6.17\cdot 10^{-5}$
$Z_{5}$ | 16641 | 0.0736s | 3.98 | $2.71\cdot 10^{-2}$ | $1.24\cdot 10^{-5}$
$Z_{6}$ | 66049 | 0.3149s | 4.28 | $2.71\cdot 10^{-2}$ | $2.98\cdot 10^{-6}$
Table 1: The number of finite element basis functions (#Nodes) characterizes
the discretization. $w_{\ell}$ is the average time to compute a realization of
$Z_{\ell}$, and the cost increase factor is denoted by $w_{\ell}/w_{\ell-1}$.
The last two columns contain the estimates for the variance and bias.

Since the random field $b$ is smooth and we integrate over $y$ in (4.1), we
expect a RE expansion of the form

(4.4)
$Z_{\ell}(\omega)=Z(\omega)+c_{2}(\omega)2^{-2\ell}+c_{3}(\omega)2^{-4\ell}+o(2^{-4\ell})\quad\text{as
}\ell\rightarrow+\infty.$

If the expansion in (4.4) holds, we have $q=3$ with rates $\gamma^{1}=0$,
$\gamma^{2}=2$ and $\gamma^{3}=4$. Therefore we expect to observe the
following rates according to Lemma 3.8:

(4.5) $\displaystyle\operatorname{Bias}(Z_{\ell})\simeq 2^{-2\ell},$
$\displaystyle\operatorname{Var}(Z_{\ell}-Z_{\ell-1})\simeq 2^{-4\ell},$
$\displaystyle\operatorname{Bias}((v^{\ell,3})^{T}(Z_{j})_{j=1}^{\ell})\simeq
2^{-4\ell},$
$\displaystyle\operatorname{Var}((v^{\ell,3}-v^{\ell-1,3})^{T}(Z_{j})_{j=1}^{\ell})\simeq
2^{-8\ell}$

and $\operatorname{Var}(Z_{\ell})\simeq 1$. The notation
$\phi_{1}(\ell)\simeq\phi_{2}(\ell)$ means that there exist constants
$c_{1},c_{2}>0$ such that

$\phi_{1}(\ell)\leq c_{1}\phi_{2}(\ell),\quad\phi_{2}(\ell)\leq
c_{2}\phi_{1}(\ell),\quad\text{for all }\ell\text{ sufficiently large}.$

Corollary 3.10 can be applied if (4.5) is true, which we assume in the
remainder of this section. We provide a numerical verification of the rates in
Figure 1. A formal proof of the rates in (4.5) or the expansion in (4.4) is
beyond the scope of this paper.

Figure 1: Biases, variances and work w.r.t. the refinement level $\ell$.

###  4.1 Complexity of estimators

We now study the computational complexity of various unbiased estimators to
approximate $\mu_{\ell}=\mathbb{E}[Z_{\ell}]$. First, we consider Monte Carlo
$\widehat{\mu}^{\operatorname{MC}}_{\ell}$, Multilevel Monte Carlo
$\widehat{\mu}^{\operatorname{MLMC}}_{\ell}$ and the SAOB
$\widehat{\mu}^{\operatorname{SAOB,2}}_{\ell}$,
$\widehat{\mu}^{\operatorname{SAOB,3}}_{\ell}$, and
$\widehat{\mu}^{\operatorname{SAOB}}_{\ell}$ with coupling numbers $q=2$,
$q=3$ and $q=+\infty$, respectively. We also consider the Multifidelity Monte
Carlo estimator $\widehat{\mu}^{\operatorname{MFMC}}_{\ell}$ in [14].

We measure the accuracy of $\widehat{\mu}_{\ell}$, an unbiased estimator of
$\mu_{\ell}=\mathbb{E}[Z_{\ell}]$, by the MSE

(4.6)
$\mathbb{E}[(\widehat{\mu}_{\ell}-\mathbb{E}[Z])^{2}]=\operatorname{Bias}(Z_{\ell})^{2}+\operatorname{Var}(\widehat{\mu}_{\ell}).$

For each level $\ell=1,\dots,L$ we ensure that
$\operatorname{Bias}(Z_{\ell})^{2}=\operatorname{Var}(\widehat{\mu}_{\ell})$
resulting in a MSE equal to $2\cdot\operatorname{Bias}(Z_{\ell})^{2}$. Let
$\varepsilon>0$ denote a given tolerance. To achieve a MSE
$\leq\varepsilon^{2}$ in (4.6) all estimators in this section have
computational costs bounded by an expression of the form

(4.7)
$\text{Cost}(\widehat{\mu}_{\ell})\leq\phi(\varepsilon)=c\varepsilon^{-\gamma^{\operatorname{Cost}}/\gamma^{\operatorname{Bias}}}+c\begin{cases}\varepsilon^{-2},\quad&\text{if
}\gamma^{\operatorname{Cost}}<\gamma^{\operatorname{Var}},\\\
\varepsilon^{-2}(\log\varepsilon)^{2},\quad&\text{if
}\gamma^{\operatorname{Cost}}=\gamma^{\operatorname{Var}},\\\
\varepsilon^{-2-({\gamma^{\operatorname{Cost}}-\gamma^{\operatorname{Var}}})/\gamma^{\operatorname{Bias}}},\quad&\text{if
}\gamma^{\operatorname{Cost}}>\gamma^{\operatorname{Var}}.\end{cases}$

Note that the cost bound $\phi(\varepsilon)$ in (4.7) consists of two terms.
The first summand arises from ceiling the optimal, but possibly fractional
number of samples, and the second summand arises from the constraint MSE
$\leq\varepsilon^{2}$.

In addition to the unbiased estimators for $\mu_{\ell}$ we also consider
$\widehat{\mu}^{\operatorname{RE,3}}_{v^{\ell,3}}$,
$\widehat{\mu}^{\operatorname{SAOB,2}}_{v^{\ell,3}}$,
$\widehat{\mu}^{\operatorname{SAOB,3}}_{v^{\ell,3}}$ and
$\widehat{\mu}^{\operatorname{SAOB}}_{v^{\ell,3}}$, which are unbiased
estimators for $(v^{\ell,3})^{T}\mu$. We record the complexity bounds for the
estimators in Table 2. For $\gamma^{\operatorname{Cost}}=2$ we observe that
all estimators have the optimal complexity $\varepsilon^{-2}$ except the MC
estimator with complexity $\varepsilon^{-3}$. The complexity bounds are
verified numerically in Figure 2 where we plot the computed costs. We
generated this plot by first computing the optimal fractional number of
samples for each estimator and then ceiling this number.

A more interesting setting is achieved by using the artificial, increased cost
per sample

(4.8) $\operatorname{Cost}(Z_{\ell})=10^{-6}\cdot 2^{6\ell},$

which corresponds to a rate of $\gamma^{\operatorname{Cost}}=6$. We list the
asymptotic upper bound for the cost $\phi$ of all estimators in Table 2. Note
that the bounds for MC, MLMC and the RE estimators follows from Section 3. The
bound for MFMC is proved in [13]. For the $\operatorname{SAOB}$ we use the
procedure outlined at the beginning of Section 3, that is, we consider a
sequence of estimators whose complexity is an upper bound for the SAOBs. In
particular, for $\widehat{\mu}^{\operatorname{SAOB}}_{L}$ we obtain the rate
$\varepsilon^{-3}$ by comparison with
$\widehat{\mu}^{\operatorname{MLMC}}_{L}$, and for
$\widehat{\mu}^{\operatorname{SAOB}}_{v^{L,3}}$ we obtain the rate
$\varepsilon^{-2}$ by comparison with
$\widehat{\mu}^{\operatorname{RE,3}}_{v^{L,3}}$. Remark 3.3 shows that
$\widehat{\mu}^{\operatorname{SAOB,q}}_{\alpha}$ has equal or smaller
complexity compared to $\widehat{\mu}^{\operatorname{RE,q}}_{\alpha}$ and thus
we use the bound for the latter. The special case
$\widehat{\mu}^{\operatorname{RE,2}}_{L}=\widehat{\mu}^{\operatorname{MLMC}}_{L}$
shows that $\widehat{\mu}^{\operatorname{SAOB,q}}_{L}$ has costs not exceeding
the costs of MLMC for $q\geq 2$.

The resulting costs of the estimators are plotted in Figure 3, thereby
confirming the claims made in Table 2. Note that the cost for
$\widehat{\mu}^{\operatorname{MLMC}}_{L}$,
$\widehat{\mu}^{\operatorname{SAOB,3}}_{L}$ and
$\widehat{\mu}^{\operatorname{SAOB}}_{L}$ is of order $\varepsilon^{-3}$ and
hence suboptimal. However, for $\widehat{\mu}^{\operatorname{SAOB,3}}_{L}$ and
$\widehat{\mu}^{\operatorname{SAOB}}_{L}$ the suboptimal cost is a result of
ceiling the number of samples: Since at least one high fidelity model has to
be evaluated, the total cost is lower bounded by $\varepsilon^{-3}$. Without
ceiling, the total cost would in fact be optimal with an order of
$\varepsilon^{-2}$. This is not the case for
$\widehat{\mu}^{\operatorname{MLMC}}_{L}$ where the suboptimal cost is of
order $\varepsilon^{-3}$ with or without ceiling the number of samples.

In Figure 3 we also plot $\operatorname{SAOB}$(*), which is the
$\operatorname{SAOB}$ without ceiling the number of samples. This estimator
requires the evaluation of $\approx 5\cdot 10^{-4}$ samples of the high-
fidelity model which is not a natural number and thus impossible in practice.
This shows that the complexity is $\varepsilon^{-2}$ if we drop the term
$\varepsilon^{-\gamma^{\operatorname{Cost}}/\gamma^{\operatorname{Bias}}}=\varepsilon^{-3}$
associated with ceiling. Importantly, we can achieve the optimal complexity by
using the bias vector $v^{\ell,3}$ instead of $v^{\ell,2}=e_{\ell}$, which is
shown in the right image of Figure 3. This change improves the overall
complexity to $\varepsilon^{-2}$, since the rounding cost does not dominate.
Finally, the $\widehat{\mu}^{\operatorname{SAOB,2}}_{v^{\ell,3}}$ does not
couple three models and thus only achieves a rate of $\varepsilon^{-2.5}$,
which is a consequence of the small rate $\gamma^{\operatorname{Var}}=4$.

Estimator |  |  |  | True $\gamma^{\text{Cost}}=2$ | Artificial $\gamma^{\text{Cost}}=6$ | Justification
---|---|---|---|---|---|---
| $q$ | $\gamma^{\text{Bias}}$ | $\gamma^{\text{Var}}$ | $\phi(\varepsilon)$ | $\phi(\varepsilon)$ |
$\widehat{\mu}^{\operatorname{MC}}_{L}$ | 1 | 2 | 0 | $\varepsilon^{-1}+\varepsilon^{-3}$ | $\varepsilon^{-3}+\varepsilon^{-5}$ | Corollary 3.4
$\widehat{\mu}^{\operatorname{MLMC}}_{L}$ | 2 | 2 | 4 | $\varepsilon^{-1}+\varepsilon^{-2}$ | $\varepsilon^{-3}+\varepsilon^{-3}$ | Corollary 3.6
$\widehat{\mu}^{\operatorname{MFMC}}_{L}$ | $L$ | 2 | 4 | $\varepsilon^{-1}+\varepsilon^{-2}$ | $\varepsilon^{-3}+\varepsilon^{-3}$ | [13]
$\widehat{\mu}^{\operatorname{SAOB,2}}_{L}$ | 2 | 2 | 4 | $\varepsilon^{-1}+\varepsilon^{-2}$ | $\varepsilon^{-3}+\varepsilon^{-3}$ | Corollary 3.7
$\widehat{\mu}^{\operatorname{SAOB,3}}_{L}$ | 3 | 2 | 8 | $\varepsilon^{-1}+\varepsilon^{-2}$ | $\varepsilon^{-3}+\varepsilon^{-2}$ |  Corollary 3.15 (*)
$\widehat{\mu}^{\operatorname{SAOB}}_{L}$ | $L$ | 2 | 8 | $\varepsilon^{-1}+\varepsilon^{-2}$ | $\varepsilon^{-3}+\varepsilon^{-2}$ |  Corollary 3.15 (*)
$\widehat{\mu}^{\operatorname{MC}}_{v^{L,3}}$ | 1 | 4 | 0 | $\varepsilon^{-0.5}+\varepsilon^{-2.5}$ | $\varepsilon^{-1.5}+\varepsilon^{-3.5}$ | Corollary 3.4
$\widehat{\mu}^{\operatorname{RE,2}}_{v^{L,3}}$ | 2 | 4 | 4 | $\varepsilon^{-0.5}+\varepsilon^{-2}$ | $\varepsilon^{-1.5}+\varepsilon^{-2.5}$ | Corollary 3.14
$\widehat{\mu}^{\operatorname{RE,3}}_{v^{L,3}}$ | 3 | 4 | 8 | $\varepsilon^{-0.5}+\varepsilon^{-2}$ | $\varepsilon^{-1.5}+\varepsilon^{-2}$ | Corollary 3.10
$\widehat{\mu}^{\operatorname{SAOB,2}}_{v^{L,3}}$ | 2 | 4 | 4 | $\varepsilon^{-0.5}+\varepsilon^{-2}$ | $\varepsilon^{-1.5}+\varepsilon^{-2.5}$ | Corollary 3.15
$\widehat{\mu}^{\operatorname{SAOB,3}}_{v^{L,3}}$ | 3 | 4 | 8 | $\varepsilon^{-0.5}+\varepsilon^{-2}$ | $\varepsilon^{-1.5}+\varepsilon^{-2}$ | Corollary 3.11
$\widehat{\mu}^{\operatorname{SAOB}}_{v^{L,3}}$ | $L$ | 4 | 8 | $\varepsilon^{-0.5}+\varepsilon^{-2}$ | $\varepsilon^{-1.5}+\varepsilon^{-2}$ | Corollary 3.11
Table 2:  Bias, variance, and cost rates together with the resulting
complexity bounds in (4.7) to achieve a MSE of order $\varepsilon^{2}$. The
cost bound $\phi$ is given as order of $\varepsilon$, where the left or right
term is associated with the corresponding term of (4.7). The coupling number
$q$ denotes the maximal number of models that are evaluated using the same
event $\omega$. The justification marked with (*) requires Corollary 3.14 and
Corollary 3.15 to hold with $s>t$, which we did not prove but we conjecture
that it is true. Alternatively, we obtain the exact same upper bound for the
total complexity, namely $\varepsilon^{-2}$ for $\gamma^{\text{Cost}}=2$ and
$\varepsilon^{-3}$ for $\gamma^{\text{Cost}}=6$ by using the bounds
$\varepsilon^{-1}+\varepsilon^{-2}$ for $\gamma^{\text{Cost}}=2$ and
$\varepsilon^{-3}+\varepsilon^{-3}$ for $\gamma^{\text{Cost}}=6$ which follow
from Corollary 3.7.

Figure 2: True cost $\gamma^{\text{Cost}}=2$: unbiased estimators of
$\mu_{\ell}$ on the left and $(v^{\ell,3})^{T}\mu$ on the right. Note that the
MC estimator on the left is not equal to the MC estimator on the right due to
the different bias. The RE, 2 estimator on the right is the weighted RE
estimator $\widehat{\mu}^{\operatorname{RE,2}}_{v^{\ell,3}}$.

Figure 3: Artificial cost $\gamma^{\text{Cost}}=6$: unbiased estimators of
$\mu_{\ell}$ on the left and $(v^{\ell,3})^{T}\mu$ on the right. Note that the
MC estimator on the left is not equal to the MC estimator on the right due to
the different bias. The RE, 2 estimator on the right is the weighted RE
estimator $\widehat{\mu}^{\operatorname{RE,2}}_{v^{\ell,3}}$.
$\operatorname{SAOB},3$ and $\operatorname{SAOB}$ have almost identical cost.
The costs for $\operatorname{SAOB}$(*) on the left is computed without ceiling
the number of samples.

######  Remark 4.1.

Let us describe an informal derivation of the rates
$\gamma^{\operatorname{Bias}}$ and $\gamma^{\operatorname{Var}}$ for (4.7) in
Table 2. For MC, MLMC and the RE estimators the rates follow from Section 3
assuming the rates in (4.5) are true. For MFMC the authors of [13] show that
the rates are the same as for MLMC. For the SAOB the rate
$\gamma^{\operatorname{Bias}}$ is equal to the bias rate for the RE
estimators. The rate $\gamma^{\operatorname{Var}}$ follows from the coupling
number $q$, i.e. using a single model results in
$\gamma^{\operatorname{Var}}=0$, coupling two models gives
$\gamma^{\operatorname{Var}}=4$ and coupling three or more models gives
$\gamma^{\operatorname{Var}}=8$. This reasoning is valid for all examined
estimators except for MFMC.

Notice that we can only prove upper bounds on the cost of
$\widehat{\mu}^{\operatorname{SAOB,q}}_{\alpha}$. However, our numerical
experiment suggests that these cost bounds are sharp in some cases. We examine
this finding from a different angle in the next section.

##  5 Convergence of SAOB to the RE estimator

In this section we conduct a numerical experiment to show that in some cases
the coefficients $\beta^{k}$ for the SAOB converge to the coefficients of RE
estimators. Recall that if the aforementioned estimators use the same model
groups $S^{k}$ they can be expressed as

$\widehat{\mu}_{L}=\sum_{k=1}^{L}\sum_{\ell\in
S^{k}}\beta^{k}_{\ell}\frac{1}{m_{k}}\sum_{i=1}^{m_{k}}Z_{\ell}(\omega^{k}_{i}),$

where the coefficients $\beta^{k}$ depend on the estimator. We compare
unbiased estimators for $\mathbb{E}[Z_{L}]$ and thus
$e_{L}=\sum_{k=1}^{L}\beta^{k}$. For $\operatorname{SAOB},q$ and
$\operatorname{RE},q$ with $q=2,3,4$ we define the difference in the
coefficients $\beta^{k}$ as follows

(5.1)
$r^{q}(\ell_{0})=\left(\sum_{k=1}^{L}\|\beta^{k,\operatorname{SAOB,q}}(\ell_{0})-\beta^{k,\operatorname{RE},q}\|^{2}\right)^{1/2},$

where $\ell_{0}$ is a discretization constant determining the initial
accuracy. We further examine the relative loss of the variance using the RE
estimator instead of the optimal $\operatorname{SAOB}$. To this end we define

(5.2)
$e^{q}(\ell_{0})=\frac{\operatorname{Var}(\widehat{\mu}^{\operatorname{RE,q}}_{L}(\ell_{0}))-\operatorname{Var}(\widehat{\mu}^{\operatorname{SAOB,q}}_{L}(\ell_{0}))}{\operatorname{Var}(\widehat{\mu}^{\operatorname{SAOB,q}}_{L}(\ell_{0}))}.$

We want to show the convergence of (5.1) and (5.2) to zero for
$\ell_{0}\rightarrow+\infty$. We achieve this by using an academic toy model
with $L=4$ and $\ell=1,\dots,4$, where

(5.3) $\displaystyle Z_{\ell}(\omega)$
$\displaystyle=Z(\omega)+\sum_{k=2}^{4}c_{k}(\omega)2^{-(k-1)(\ell+\ell_{0})}+0.1\xi_{\ell}(\omega)2^{-3(\ell+\ell_{0})},$
$\displaystyle(Z,c_{2},c_{3},c_{4})^{T}$ $\displaystyle\sim N(0,Q),\quad
Q_{ij}=\exp(-|i-j|),\quad i,j=1,\dots,4,$ $\displaystyle\xi_{\ell}$
$\displaystyle\sim N(0,1).$

We assume that $\xi_{i}$ and $\xi_{j}$ are independent for $i\not=j$ and that
the $\xi_{\ell}$ are independent of $Z,c_{2},c_{3},c_{4}$. The parameter
$\ell_{0}$ controls the accuracy of the coarsest model. This model satisfies
Section 3.3 with $q=4$, $\gamma^{2}=1$, $\gamma^{3}=2$ and $\gamma^{4}=3$.
Furthermore, we fix artificial costs of
$\operatorname{Cost}(Z_{\ell})=4^{\ell-1}$ and remark that the model
covariance matrix $C$ can be computed analytically.

Recall that $\operatorname{RE},2$ is the MLMC estimator if the bias
$\alpha=e_{\ell}$. We plot the computed values of (5.1) and (5.2) for
different values of $\ell_{0}=0,\dots,6$ in Figure 4. We conclude that both
quantities converge to zero showing that for large $\ell_{0}$ the estimators
$\operatorname{SAOB},q$ and $\operatorname{RE},q$ are almost identical for the
problem in (5.3). Notice that we allow fractional samples
$m_{k}\in\mathbb{R}_{\geq 0}$ and do not ceil. This does not change the
results of this section in a fundamental way, since we could also arbitrarily
increase the budget or the variance by scaling.

Figure 4: The left image shows the convergence of the coefficients of the
$\operatorname{SAOB},q$ to the coefficients of the $\operatorname{RE},q$
estimators w.r.t. $\ell_{0}$. The right image shows the relative variance to
the corresponding $\operatorname{SAOB},q$ estimator. A reference rate is drawn
dashed. Here $r^{2}$ and $e^{2}$ is the coefficient and relative variance
increment of $\operatorname{RE},2$, which is the MLMC estimator.

######  Remark 5.1.

We informally state an explanation for the convergence of the coefficients of
$\operatorname{SAOB},q$ to $\operatorname{RE},q$. The variance of
$\operatorname{RE},q$ satisfies

$\operatorname{Var}(\widehat{\mu}^{\operatorname{RE,q}}_{L})=\sum_{k=1}^{L}\frac{\operatorname{Var}\left(\sum_{\ell\in
S^{k}}\beta^{k}_{\ell}Z_{\ell}\right)}{m_{k}}=\sum_{k=1}^{L}\frac{\mathcal{O}(2^{-2\gamma^{q}\ell_{0}})}{m_{k}}.$

To achieve the rate $2^{-2\gamma^{q}\ell_{0}}$ for every variance term the
coefficients $\beta^{k}$ are often asymptotically uniquely determined. Scaling
these coefficients down is not allowed since we have the bias constraint
$e_{L}=\sum_{k=1}^{L}\beta^{k}$. The proof of Lemma 3.8 shows that the linear
combinations $(v^{k,q}-v^{k-1,q})$ needed to achieve an increased order are
uniquely determined with the exception of some degenerate cases, i.e., if the
$c_{k}$ in Section 3.3 are zero or linearly dependent.

######  Remark 5.2.

The numerical results in this section suggest that the estimators
$\operatorname{RE},q$ converge to $\operatorname{SAOB},q$ if the coarsest grid
is fine enough. We thus suspect that for hierarchical models
$\operatorname{RE},q$ and $\operatorname{SAOB},q$ often have the same
asymptotic cost. We also observed this numerically in Section 4.1.
$\operatorname{SAOB}$ may couple more models than $\operatorname{RE},q$,
however if the remainder term in Section 3.3 cannot be removed by linearly
combining more models we expect

$(\beta^{k})^{T}C\beta^{k}\geq c2^{-2k\gamma^{q}},$

and thus $\widehat{\mu}^{\operatorname{SAOB}}_{\alpha}$ and
$\widehat{\mu}^{\operatorname{RE,q}}_{\alpha}$ should have the same variance
rate of $\gamma^{\operatorname{Var}}=2\gamma^{q}$. They also have the same
bias rate $\gamma^{\operatorname{Bias}}$. For $\operatorname{SAOB}$ it is
reasonable to assume that $S^{k}=\\{1,\dots,k\\}$ for $k=1,\dots,L$, since
$\operatorname{SAOB}$ should use all models $Z_{1},\dots,Z_{L}$ if $Z_{1}$ is
fine enough and the asymptotic cost of $S^{k}$ is given by its finest model,
i.e., we may add all coarse models without an asymptotic cost increase.
Therefore both estimators have equal bias, variance and cost rates and thus
$\widehat{\mu}^{\operatorname{SAOB}}_{\alpha}$ has costs asymptotically equal
to $\widehat{\mu}^{\operatorname{RE,q}}_{\alpha}$ where $q$ is maximal such
that Section 3.3 is satisfied. The reasoning in this remark is however not a
formal proof.

##  6 Comparison with ACV estimators

We return to the example in Section 4 and compare the SAOBs and RE estimators
with the ACV estimators [9], in particular, the ACV-IS, ACV-MF and ACV-KL. For
brevity we provide only the definition of the ACV-MF estimator,

(6.1) $\widehat{\mu}^{\operatorname{ACV-
MF}}_{L}=\frac{1}{n_{L}}\sum_{i=1}^{n_{L}}Z_{L}(\omega^{i})+\sum_{\ell=1}^{L-1}\beta_{\ell}\left(\frac{1}{n_{L}}\sum_{i=1}^{n_{L}}Z_{\ell}(\omega^{i})-\frac{1}{n_{\ell}}\sum_{i=1}^{n_{\ell}}Z_{\ell}(\omega^{i})\right),$

and refer to [9, Def. 2] and [9, Def. 4] for the definition of the ACV-IS and
ACV-KL, respectively. The comparison is carried out in a separate section
since the ACV estimators have been developed very recently, and no asymptotic
complexity results are known to date. Recall that an important motivation for
the introduction of the ACV estimators is their increased variance reduction
capacity compared to MLMC and MFMC. In fact, the ACV estimators in [9] reach
the exact same lower variance bound in the infinite low fidelity data limit as
the multilevel BLUEs (cf. [19, Sec. 4]). However, the improved variance
reduction does not necessarily result in a smaller asymptotic cost for
parametric, PDE-based models which we study in this paper. We observed this
already in the numerical example in Section 4 for the SAOBs where the
asymptotic complexity is improved in some settings, but depends on the bias,
variance and cost rate in the considered example.

We plot the cost of the ACV estimators in Figure 5 for the artificial cost
rate $\gamma^{\operatorname{Cost}}=6$ without and with ceiling the number of
samples.

Figure 5: The left image shows the MSE for the computed costs not ceiling the
number of samples and the right image with ceiling for the artificial cost
$\gamma^{\operatorname{Cost}}=6$. For small $\varepsilon$ the cost for the
ACV-MF and ACV-KL estimator are identical. Both plots display the estimators
with a bias of $\mu_{\ell}$.

We can clearly see that in this example the ACV estimators have the asymptotic
cost of order $\varepsilon^{-3}$. This is the same cost as the other
estimators, and in particular, the MLMC estimator. The results in the right
plot of Figure 5 should be compared with the left plot in Figure 3.

We now offer a possible explanation for this observation by looking at the
coefficients $\beta^{k}$ of the various estimators. Since the ACV estimators
are also linear unbiased estimators for $\mu_{L}$ we compare them to the SAOB,
4 and the weighted RE, 4 estimator. Each of these three estimators can be
written as

$\widehat{\mu}_{L}=\sum_{k=1}^{6}\sum_{\ell\in
S^{k}}\beta^{k}_{\ell}\frac{1}{m_{k}}\sum_{i=1}^{m_{k}}Z_{\ell}(\omega^{i}_{k}),$

with different model groups $S^{1},\dots,S^{6}$ for every estimator. We plot
the resulting coefficients $\beta^{k}$ in Figure 6.

Figure 6: Coefficients in front of every model group for some estimators in
Figure 5 where the smallest value of $\varepsilon$ has been used. A column
represents the models in the model group. For example, $S^{5}=\\{2,\dots,5\\}$
for RE, 4 and SAOB, 4 and $S^{5}=\\{1,\dots,6\\}$ for ACV-MF. The coefficient
$\beta^{4}_{3}=-1.83$ for both SAOB, 4 and RE, 4. Blank entries are not
contained in the respective model group. Entries with values $0.00$ or $-0.00$
are in the respective model group but have a small coefficient. The bias
constraint ensures that the column vectors add to $e_{L}$ (up to small
rounding errors due to the display of only two decimal places).

Using the definition of the ACV-MF estimator in (6.1) and using $n_{\ell}\geq
n_{L}$ it is straightforward to verify that the ACV-MF coefficients satisfy
the sign pattern

(6.2)
$\operatorname{sign}(\beta^{L}_{\ell})=-\operatorname{sign}(\beta^{k}_{\ell})\quad\text{
for all }k\in\\{1,\dots,L-1\\}\text{ with }\ell\in S^{k}.$

However, the RE estimator exhibits a chequerboard pattern w.r.t. the sign of
the coefficients, which cannot be satisfied under the conditions (6.2). This
reduces the variance reduction capability of the ACV-MF estimator, since some
linear combinations of models are prohibited by construction. Thus, the
estimator variance

$\operatorname{Var}\left(\sum_{\ell\in
S^{k}}\beta^{k}_{\ell}Z_{\ell}\right)=(\beta^{k})^{T}C\beta^{k}$

may be larger than necessary. If the linear combinations that asymptotically
decrease the estimator variance require multiple sign changes, which seems to
be the case according to Section 5 and Figure 6, then the smaller complexity
of the RE estimators cannot be achieved by the ACV-MF estimator. Since for the
ACV-MF we can have a single sign change similar to the sign pattern in MLMC,
we believe that the asymptotic complexity of the ACV-MF estimator is in
general equal to the complexity of MLMC for the parametric PDE-based models in
this paper. We remark that the sign pattern (6.2) holds also true for the ACV-
IS and ACV-KL estimator. Hence, a performance similar to the ACV-MF can be
expected.

##  7 Conclusions

In this paper we study the asymptotic properties of the multilevel BLUE for
the expectation of scalar-valued, PDE-based random outputs. The main tool of
our analysis is a pathwise expansion of the random output in terms of a
discretization parameter, linking a collection of models associated with the
output. We use the idea of Richardson extrapolation (RE) combined with MLMC
theory to analyse the complexity of a certain RE estimator. This allows us to
bound the asymptotic complexity of the SAOB which is not worse than the
complexity of the RE estimators and in particular MLMC. Numerical experiments
with a smooth PDE-based output in 2D space suggest that a RE type of expansion
for the quantity of interest holds both in terms of the bias and variance. If
we use the true sample cost, then the SAOB, MLMC and MFMC estimators give the
optimal complexity with respect to the mean-square error. For an artificial,
increased cost, the SAOB and RE estimators have a smaller complexity compared
to MLMC and MFMC. However, it remains an open research question whether the
SAOB has an asymptotically smaller complexity than the RE estimator.

## Acknowledgements

The authors thank Michael Ulbrich for the finite element code that was used to
solve the PDE and to sample from the mean zero Gaussian random field with
Matern 3/2 covariance in Section 4.

##  Appendix A Proof of Lemma 2.2

First, we remark that (2.10) is valid if $\beta^{k}=0$, since then by our
convention $m^{*}_{k}=0$. On the other hand, if $\beta^{k}\not=0$ and since
$C$ is positive definite, we have

$\lim_{m_{k}\rightarrow 0}J(m)=+\infty.$

Furthermore, since (2.9) clearly has a feasible point, the number of samples
is lower bounded,

$m^{*}_{k}>0\quad\text{for all }k\in\\{1,...,K\\}\text{ with
}\beta^{k}\not=0.$

Similarly, since $W_{k}>0$, the cost constraint ensures the upper bound

$m^{*}_{k}<c\quad\text{for all }k\in\\{1,...,K\\}.$

Hence (2.9) is an optimization problem over a compact set and since $J$ is
convex and thus continuous, a minimizer $m^{*}$ exists. The function $J$ is
monotonically decreasing in $m$ and thus the cost constraint is satisfied with
equality at a minimizer

$\sum_{k=1}^{K}m^{*}_{k}W_{k}=p,$

which allows us to replace the inequality constraint with an equality
constraint. The Karush–Kuhn–Tucker conditions [3, Section 5.5] are necessary
and sufficient for the convex optimization problem (2.9). Thus, with Lagrange-
Multipliers $\lambda\in\mathbb{R}$ and $\xi_{1},\dots,\xi_{K}\in\mathbb{R}$,

(A.1) $\displaystyle-\frac{(\beta^{k})^{T}C\beta^{k}}{(m^{*}_{k})^{2}}+\lambda
W_{k}-\xi_{k}$ $\displaystyle=0\qquad\text{for all }k=\\{1,\dots,K\\}\text{
with }\beta^{k}\not=0,$ (A.2) $\displaystyle\sum_{k=1}^{K}m^{*}_{k}W_{k}$
$\displaystyle=p,$ (A.3) $\displaystyle m^{*}_{k}\geq 0,\quad\xi_{k}$
$\displaystyle\geq 0,\quad\xi_{k}m^{*}_{k}=0\qquad\text{for all
}k=\\{1,\dots,K\\}\text{ with }\beta^{k}\not=0$

and $m^{*}_{k}=0$ if $\beta^{k}=0$. For $\beta^{k}\not=0$ we have already
shown $m^{*}_{k}>0$, thus $\xi_{k}=0$ and (A.1) reads

(A.4) $m^{*}_{k}=\left(\frac{(\beta^{k})^{T}C\beta^{k}}{\lambda
W_{k}}\right)^{1/2}.$

This expression is well defined for $\lambda>0$, which is valid since
otherwise (A.1) cannot be satisfied at the minimizer $m^{*}$. We insert (A.4)
into the cost constraint (A.2) and arrive at

(A.5)
$\frac{1}{\lambda^{1/2}}=\frac{p}{\sum_{k=1}^{K}((\beta^{k})^{T}C\beta^{k}W_{k})^{1/2}}.$

Here the denominator is positive since $W_{k}>0$, $C$ is positive definite and
$\beta^{k}\not=0$ for some index $k$. The last statement follows from the bias
constraint $\alpha=\sum_{k=1}^{K}\beta^{k}$ and $\alpha\not=0$ from the
assumptions of this lemma. We now insert (A.5) into (A.4) to obtain the result
(2.10), which also shows the uniqueness of a minimizer. Inserting this $m^{*}$
into $J$ then shows (2.11).

##  Appendix B Proof of Lemma 3.8

We first show that the following statement is true for all $k=1,\dots,L$ and
$\ell_{0}\in\\{0,\dots,L-\ell\\}$:

(B.1)
$\sum_{\ell=1}^{L-\ell_{0}}v^{k,q}_{\ell}Z_{\ell+\ell_{0}}=Z+\sum_{j=k+1}^{q-1}c^{k}_{j}2^{-(k+\ell_{0})\gamma^{j}}+\mathcal{O}(2^{-(k+\ell_{0})\gamma^{q}}),$

where $c^{k}_{j}$ are random variables with bounded second moment. The
remainder term $\mathcal{O}(\cdot)$ also has a bounded second moment. The
statement for $k=1$ is Section 3.3 (i). Now let $1<k<q$ and assume the
induction hypothesis is true for $k-1$. Observe that $v^{k-1}_{L-\ell_{0}}=0$
if $\ell_{0}\in\\{0,\dots,L-k\\}$ due to the recursion (3.6). Defining
$v^{k-1,q}_{0}=0$ we have

$\displaystyle\sum_{\ell=1}^{L-\ell_{0}}v^{k,q}_{\ell}Z_{\ell+\ell_{0}}$
$\displaystyle=\frac{1}{2^{\gamma^{k}}-1}\sum_{\ell=1}^{L-\ell_{0}}(2^{\gamma^{k}}v^{k-1,q}_{\ell-1}-v^{k-1,q}_{\ell})Z_{\ell+\ell_{0}}$
$\displaystyle=\frac{1}{2^{\gamma^{k}}-1}\left(2^{\gamma^{k}}\sum_{\ell=1}^{L-\ell_{0}}v^{k-1,q}_{\ell}Z_{\ell+\ell_{0}+1}-\sum_{\ell=1}^{L-\ell_{0}}v^{k-1,q}_{\ell}Z_{\ell+\ell_{0}}\right).$

Now we apply the induction hypothesis, which leads to

$\displaystyle\sum_{\ell=1}^{L-\ell_{0}}v^{k,q}_{\ell}Z_{\ell+\ell_{0}}$
$\displaystyle=Z+\frac{2^{\gamma^{k}}}{2^{\gamma^{k}}-1}\left(\sum_{j=k}^{q-1}c^{k-1}_{j}2^{-(k+\ell_{0})\gamma^{j}}+\mathcal{O}(2^{-(k+\ell_{0})\gamma^{q}})\right)$
$\displaystyle-\frac{1}{2^{\gamma^{k}}-1}\left(\sum_{j=k}^{q-1}c^{k-1}_{j}2^{-(k+\ell_{0}-1)\gamma^{j}}+\mathcal{O}(2^{-(k+\ell_{0}-1)\gamma^{q}})\right).$

The random coefficients thus satisfy

$c^{k}_{j}=\frac{1}{2^{\gamma^{k}}-1}\left(2^{\gamma^{k}}c^{k-1}_{j}-2^{\gamma^{j}}c^{k-1}_{j}\right),$

where a similar expression for the remainder term is valid. Notice that
$c^{k}_{k}=0$ and that the $c^{k}_{j}$ have bounded second moment as well as
the remainder term. Therefore (B.1) holds for $1\leq k<q$. For $k\geq q$
observe $v^{k-1,q}_{L-\ell_{0}}=0$ and thus

$\sum_{\ell=1}^{L-\ell_{0}}v^{k,q}_{\ell}Z_{\ell+\ell_{0}}=\sum_{\ell=1}^{L-\ell_{0}}v^{k-1,q}_{\ell-1}Z_{\ell+\ell_{0}}=\sum_{\ell=1}^{L-\ell_{0}-1}v^{k-1,q}_{\ell}Z_{\ell+\ell_{0}+1}=\sum_{\ell=1}^{L-\ell_{0}}v^{k-1,q}_{\ell}Z_{\ell+\ell_{0}+1},$

which by repeatedly applying this process allows us to reduce the case of
$k\geq q$ to $k=q-1$ by increasing the value of $\ell_{0}$. This shows (B.1).

Let us now prove the bias estimate (RE, M1) using (B.1). For $\ell_{0}=0$ we
conclude

$|(v^{k,q})^{T}\mu-\mathbb{E}[Z]|=\left|\sum_{j=k+1}^{q-1}\mathbb{E}[c^{k}_{j}]2^{-k\gamma^{j}}+\mathcal{O}(2^{-k\gamma^{q}})\right|\leq
c\begin{cases}2^{-k\gamma^{k+1}},&\text{if }k<q,\\\
2^{-k\gamma^{q}},&\text{otherwise}.\end{cases}$

The bound (RE, M1) for the case $k<q$ is obtained using the crude estimate

$2^{-k\gamma^{k+1}}=2^{-k\gamma^{k+1}}2^{k\gamma^{q}}2^{-k\gamma^{q}}\leq
1\cdot 2^{q\gamma^{q}}2^{-k\gamma^{q}}\leq c2^{-k\gamma^{q}},$

where now $c$ is independent of $k$. The variance estimate (RE, M2) can be
derived similarly from (B.1). Here the key idea is that the difference
$v^{k,q}-v^{k-1,q}$ is used to remove $Z$ from (B.1), which we use with
$\ell_{0}=0$, $k$ and $k-1$,

$\displaystyle\sum_{\ell=1}^{L}(v^{k,q}_{\ell}-v^{k-1,q}_{\ell})Z_{\ell}$
$\displaystyle=\sum_{j=k+1}^{q-1}c^{k}_{j}2^{-k\gamma^{j}}+\mathcal{O}(2^{-k\gamma^{q}})-\sum_{j=k}^{q-1}c^{k-1}_{j}2^{-(k-1)\gamma^{j}}+\mathcal{O}(2^{-(k-1)\gamma^{q}})$
$\displaystyle=\sum_{j=k}^{q-1}\widetilde{c}^{k-1}_{j}2^{-(k-1)\gamma^{j}}+\mathcal{O}(2^{-(k-1)\gamma^{q}})$

for suitably defined random variables $\widetilde{c}^{k-1}_{j}$ with bounded
second moment. Thus taking the variance yields the desired result. Finally,
the estimate (RE, M3) follows from the definition of the model groups in (3.8)
and the geometric cost of $Z_{\ell}$ (3.5) in Section 3.3 (ii),

$\displaystyle W_{k}$
$\displaystyle=\operatorname{Cost}\left(S^{k}\right)=\operatorname{Cost}\left(Z_{\max\\{k-q+1,1\\}},Z_{\max\\{k-q+2,1\\}},\dots,Z_{k}\right)=\sum_{\ell=\max\\{k-q+1,1\\}}^{k}\operatorname{Cost}(Z_{\ell})$
$\displaystyle\leq c\sum_{\ell=1}^{k}2^{\ell\gamma^{\operatorname{Cost}}}\leq
c2^{k\gamma^{\operatorname{Cost}}}.$

##  Appendix C Proof of Lemma 3.12

We only have to prove the boundedness of $a_{k}$. For $t=s$ we have $a_{k}=1$
and the standard RE estimator. For $s<t$ observe from (3.6) that there exist
coefficients $d$ such that

$\sum_{k=1}^{t}d_{k}v^{k,s}=v^{t,t}.$

We use the property of the shift matrix $D$ to obtain

$\sum_{k=L-t+1}^{L}d_{k-L+t}v^{k,s}=\sum_{k=1}^{t}d_{k}v^{L-t+k,s}=D^{L-t}v^{t,t}=v^{L,t}.$

We combine this with the basis property of the differences $v^{k,s}-v^{k-1,s}$
in (3.10), we define $a_{L+1}=0$ and use $v^{0,s}=0$ to rewrite the sum in
terms of differences of the $a_{k}$,

$v^{L,t}=\sum_{k=1}^{L}a_{k}(v^{k,s}-v^{k-1,s})=\sum_{k=1}^{L}a_{k}(v^{k,s}-v^{k-1,s})=\sum_{k=1}^{L}(a_{k}-a_{k+1})v^{k,s}.$

We summarize the chain of equations

$\sum_{k=L-t+1}^{L}d_{k-L+t}v^{k,s}=\sum_{k=1}^{L}(a_{k}-a_{k+1})v^{k,s}$

and since the $v^{k,s}$ are linearly independent the last $t$ coefficients
$a_{L-t+1},\dots,a_{L}$ only depend on $d$ and are independent of $L$. The
remaining coefficients satisfy $a_{k}=a_{k+1}$ for $\ell=1,\dots,L-t$ and thus
$|a_{k}|\leq c$ with a constant $c$ independent of $k$ and the finest level
$L$.

## References

  * [1] M. Asadzadeh, A. H. Schatz, and W. Wendland, A new approach to Richardson extrapolation in the finite element method for second order elliptic problems, Mathematics of Computation, 78 (2009), pp. 1951–1973, https://doi.org/10.1090/S0025-5718-09-02241-8.
  * [2] H. Blum, Q. Lin, and R. Rannacher, Asymptotic error expansion and Richardson extrapolation for linear finite elements, Numer. Math., 49 (1986), pp. 11–37, https://doi.org/10.1007/BF01389427.
  * [3] S. Boyd and L. Vandenberghe, Convex Optimization, Cambridge University Press, Cambridge, UK, 2004, https://doi.org/10.1017/CBO9780511804441.
  * [4] C. Brezinski, Convergence acceleration during the 20th century, vol. 122, 2000, pp. 1–21, https://doi.org/10.1016/S0377-0427(00)00360-5.  Numerical analysis 2000, Vol. II: Interpolation and extrapolation.
  * [5] R. Bulirsch and J. Stoer, Numerical treatment of ordinary differential equations by extrapolation methods, Numerische Mathematik, 8 (1966), pp. 1–13, https://doi.org/10.1007/BF02165234.
  * [6] K. A. Cliffe, M. B. Giles, R. Scheichl, and A. L. Teckentrup, Multilevel Monte Carlo methods and applications to elliptic PDEs with random coefficients, Comput. Vis. Sci., 14 (2011), pp. 3–15, https://doi.org/10.1007/s00791-011-0160-x.
  * [7] M. B. Giles, Multi-level Monte Carlo path simulation, Operations Research, 56 (2008), pp. 607–617, https://doi.org/10.1287/opre.1070.0496.
  * [8] M. B. Giles, Multilevel Monte Carlo methods, Acta Numerica, 24 (2015), pp. 259–328, https://doi.org/10.1017/S096249291500001X.
  * [9] A. A. Gorodetsky, G. Geraci, M. S. Eldred, and J. D. Jakeman, A generalized approximate control variate framework for multifidelity uncertainty quantification, Journal of Computational Physics, 408 (2020), p. 109257, https://doi.org/10.1016/j.jcp.2020.109257.
  * [10] V. Lemaire and G. Pagès, Multilevel Richardson-Romberg extrapolation, Bernoulli, 20 (2017), pp. 1029–1067, https://doi.org/10.3150/16-BEJ822.
  * [11] C. Mbaye, G. Pagès, and F. Vrins, An antithetic approach of multilevel Richardson-Romberg extrapolation estimator for multidimensional SDEs, in Numerical analysis and its applications, vol. 10187 of Lecture Notes in Comput. Sci., Springer, Cham, 2017, pp. 482–491.
  * [12] E. H. Müller, R. Scheichl, and T. Shardlow, Improving multilevel Monte Carlo for stochastic differential equations with application to the Langevin equation, Proc. A., 471 (2015), pp. 20140679, 20, https://doi.org/10.1098/rspa.2014.0679.
  * [13] B. Peherstorfer, M. Gunzburger, and K. Willcox, Convergence analysis of multifidelity Monte Carlo estimation, Numerische Mathematik, 139 (2018), pp. 683–707, https://doi.org/10.1007/s00211-018-0945-7.
  * [14] B. Peherstorfer, K. Willcox, and M. Gunzburger, Optimal Model Management for Multifidelity Monte Carlo Estimation, SIAM Journal on Scientific Computing, 38 (2016), pp. A3163–A3194, https://doi.org/10.1137/15M1046472.
  * [15] B. Peherstorfer, K. Willcox, and M. Gunzburger, Survey of multifidelity methods in uncertainty propagation, inference, and optimization, SIAM Rev., 60 (2018), pp. 550–591, https://doi.org/10.1137/16M1082469.
  * [16] R. Rannacher, Extrapolation techniques in the finite element method (a survey), no. MATC7 in Proc. Summer School on Numerical Analysis, Helsinki, Univ. of Tech., 1988, pp. 80–113.
  * [17] L. F. Richardson, The approximate arithmetical solution by finite differences of physical problems involving differential equations, with an application to the stresses in a masonry dam, Philosophical Transactions of the Royal Society of London. Series A, 210 (1911), pp. 307–357, https://doi.org/10.1098/rsta.1911.0009.
  * [18] W. Romberg, Vereinfachte numerische Integration, Det Kongelige Norske Videnskabers Selskab Forhandlinger, 28 (1955), pp. 30–36.
  * [19] D. Schaden and E. Ullmann, On multilevel best linear unbiased estimators, SIAM/ASA Journal on Uncertainty Quantification, 8 (2020), pp. 601–635, https://doi.org/10.1137/19M1263534.
  * [20] M. L. Stein, Interpolation of spatial data, Springer Series in Statistics, Springer-Verlag, New York, 1999, https://doi.org/10.1007/978-1-4612-1494-6.  Some theory for Kriging.
  * [21] D. Talay and L. Tubaro, Expansion of the global error for numerical schemes solving stochastic differential equations, Stochastic Anal. Appl., 8 (1990), pp. 483–509 (1991), https://doi.org/10.1080/07362999008809220.