# A filtered multilevel Monte Carlo method for estimating the expectation of
discretized random fields111 Distributed under a CC-BY 4.0 licence. This
project has received financial support from the CNRS (Centre National de la
Recherche Scientifique) through the 80—Prime program and the French national
program LEFE (Les Enveloppes Fluides et l’Environnement).

Jérémy Briant  IRIT, UMR 5505 CNRS-INP, Toulouse, France
<EMAIL_ADDRESS><EMAIL_ADDRESS>ehouarn.simon@toulouse-inp.fr).  CECI, UMR 5318 CNRS-CERFACS, Toulouse, France
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>weaver@cerfacs.fr).
Paul Mycek  CECI, UMR 5318 CNRS-CERFACS, Toulouse, France<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>weaver@cerfacs.fr).  Mayeul Destouches  Met
Office, Exeter, UK (mayeul.destouches@metoffice.gov.uk).  Olivier Goux  CECI,
UMR 5318 CNRS-CERFACS, Toulouse, France<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>weaver@cerfacs.fr).  Serge Gratton  IRIT, UMR 5505 CNRS-INP,
Toulouse, France<EMAIL_ADDRESS>serge.gratton@toulouse-
inp.fr, ehouarn.simon@toulouse-inp.fr).  ANITI, Toulouse, France.  Selime
Gürol  CECI, UMR 5318 CNRS-CERFACS, Toulouse, France<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>weaver@cerfacs.fr).  Ehouarn Simon  IRIT,
UMR 5505 CNRS-INP, Toulouse, France<EMAIL_ADDRESS><EMAIL_ADDRESS>ehouarn.simon@toulouse-inp.fr).  Anthony T.
Weaver  CECI, UMR 5318 CNRS-CERFACS, Toulouse, France<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>weaver@cerfacs.fr).

###### Abstract

In this paper, we investigate the use of multilevel Monte Carlo (MLMC) methods
for estimating the expectation of discretized random fields. Specifically, we
consider a setting in which the input and output vectors of the numerical
simulators (models) have inconsistent dimensions across the multilevel
hierarchy. This requires the introduction of grid transfer operators borrowed
from multigrid methods. Starting from a simple one-dimensional illustration,
we demonstrate numerically that the resulting MLMC estimator deteriorates the
estimation of high-frequency (small-scale) components of the discretized
expectation field compared to a (single-level) Monte Carlo (MC) estimator. By
adapting mathematical tools initially developed for multigrid methods, we
perform a theoretical spectral analysis of the MLMC estimator of the
expectation of discretized random fields, in the specific case of linear,
symmetric and circulant simulators. This analysis provides a spectral
decomposition of the variance into contributions associated with each scale
component of the discretized field. We then propose improved MLMC estimators
using a filtering mechanism similar to the smoothing process of multigrid
methods. The filtering operators improve the estimation of both the small- and
large-scale components of the variance, resulting in a reduction of the total
variance of the estimator. These improvements are quantified for the specific
class of simulators considered in our spectral analysis. The resulting
filtered MLMC (F-MLMC) estimator is applied to the problem of estimating the
discretized variance field of a diffusion-based covariance operator, which
amounts to estimating the expectation of a discretized random field. The
numerical experiments support the conclusions of the theoretical analysis even
with non-linear simulators, and demonstrate the improvements brought by the
proposed F-MLMC estimator compared to both a crude MC and an unfiltered MLMC
estimator.

#### Keywords:

Multilevel Monte Carlo, multigrid method, random field, spectral analysis,
filtering, diffusion operator.

##  1 Introduction

Monte Carlo (MC) estimation refers to a class of statistical methods that rely
on the sampling of random quantities to construct so-called MC estimators.
Such estimators are popular and widely-used owing to their simplicity and
flexibility. However, their slow convergence in terms of their root mean
square error (RMSE) with respect to the sample size makes them inefficient or
even unaffordable for the estimation of statistics of outputs of
computationally expensive numerical simulators. Alternative techniques often
rely on the use of surrogate models, which typically consist of an approximate
functional representation of the mapping between the inputs and the output of
interest of a numerical simulator. Commonly used surrogate models include,
among others, polynomial chaos expansions [42, 34, 62, 37], Gaussian processes
[46, 47] or neural networks [54]. Once constructed, such surrogate models are
inexpensive to evaluate, so that accurate MC estimators of their output can be
obtained at limited computational cost. However, the resulting estimators
suffer from model error (bias) coming from the approximation error of the
surrogate model with respect to the actual simulator, which may be significant
for complex, highly non-linear simulators. In addition, surrogate models
generally suffer from the so-called curse of dimensionality, which makes their
construction cost grow exponentially with the stochastic input dimension. This
curse can, in principle, be mitigated by resorting to model order reduction
techniques to lower the dimensionality of the input [51, 45], but this comes
at the expense of increasing model error even further.

Another avenue consists of modifying the original, inefficient MC methods to
reduce their RMSE. A first class of techniques relies on improving the space-
covering of the stochastic input space. Such techniques include Latin
hypercube sampling [40] and quasi Monte Carlo methods [35, Chapter 5]. A
second class is that of so-called variance reduction techniques [35, Chapter
4], which consist of rewriting the original MC estimator into a transformed
estimator with smaller variance (and hence RMSE, provided the bias remains
unchanged) for an equivalent computational cost. In recent years, many
multilevel/multifidelity estimation methods relying on variance reduction
techniques have been proposed in the literature [18, 43, 20, 50]. Among them,
the multilevel Monte Carlo (MLMC) method [29, 19, 18] relies on a collection
of simulators of different fidelities. In MLMC methods, such fidelities
typically correspond to different discrete resolutions, so that the collection
of fidelities can be structured as a hierachy of so-called levels, consisting
of simulators of increasing discrete resolution, and hence of increasing
accuracy and computational cost. MLMC estimators are constructed by combining
samples (or ensembles) from different levels in a particular manner that,
under certain assumptions, reduces the variance of the resulting multilevel
estimator, while leaving the bias unchanged. The closely-related multilevel
best linear unbiased estimators (MLBLUEs) were recently proposed to combine
more general multifidelity samples to provide an optimal reduction of the
variance [50, 49]. Originally designed for the estimation of the expected
value of scalar, real-valued random variables, the MLMC methodology has since
been extended to the estimation of higher-order statistical moments [2, 3] and
variance-based global sensitivity measures [41]. The analysis of MLMC
estimators has also been extended to the estimation of statistics of random
variables with values in separable Hilbert spaces [2, 3]. Likewise, strategies
for extending MLBLUEs to the estimation of statistical moments of multiple
simulator outputs were proposed in [13, 10, 14].

The estimation of a covariance matrix is a prominent example where MLMC
approaches have started to emerge, especially in the field of ensemble data
assimilation [30, 24, 13]. Another example, which will constitute the
motivating example of this paper, is embedded in the general approach of using
a discretized linear differential operator to represent the application of a
parametric form of a covariance matrix. In particular, we focus on an approach
that uses a discretized diffusion operator to represent a covariance operator
with a parametric kernel from the Matérn family [53]. Diffusion operators are
commonly used for modelling spatial covariances in ocean data assimilation
[58] and are closely related to other techniques for modelling spatial
covariances in atmospheric data assimilation [44], geostatistical modelling
[36], inverse problems [8] and uncertainty quantification [15]. Central to the
approach is the need to extract the diagonal elements (intrinsic variances) of
the diffusion-modelled covariance matrix so that they can be used to normalize
the matrix, to transform it into a (unit-diagonal) correlation matrix. Once
the covariance matrix is properly normalized, a desired variance field,
different from the intrinsic variance field, can be imposed. The standard
method for estimating the intrinsic variances of the diffusion-based
covariance matrix is the randomization method [58, 59]. This method relies on
the MC estimation of the expectation of a discretized random field, which, in
a more abstract formulation, may be viewed as the output of a numerical
simulator, whose input is also a discretized random field.

In this paper, we are interested in applying the MLMC methodology to improve
the efficiency of estimating the expectation of discretized random fields in
the abstract setting described above. To achieve this, certain adaptations of
the original MLMC formulation must be made. First, the input and output of the
abstract numerical simulator considered in our study are discretized fields,
typically represented as $n$-dimensional vectors for practical, numerical
reasons. While multilevel approaches have been proposed in the literature to
tackle the case of numerical simulators with vector inputs and outputs [10,
14], such vectors are simply used as a means to collect scalar quantities.
However, the components of discretized fields should _not_ be considered as
individual, separate scalar quantities, precisely because they encapsulate
physical and spatial information on the underlying continuous fields that are
discretized. We stress that the vectors at hand represent discretized fields
and should be thought of as discrete signals. Second, because the input and
output discretized fields have dimensions that depend on the discrete
resolution of the numerical simulator, designing a multilevel estimator
requires the introduction of grid transfer operators, which may be borrowed
from multigrid methods [6, 26, 55, 7]. These considerations raise crucial
questions about the effect of grid transfer operators on the different scales
(or, equivalently, frequencies) of the discretized fields (or signals) that
are transferred between grids of different resolutions, and their impact on
the resulting multilevel estimator. We address these important questions that,
to the best of our knowledge, have not been previously considered in the
literature. In particular, we present a numerical illustration of the effects
of grid transfer operators on the variance of the MLMC estimator, based on the
spectral decomposition of the estimator’s variance in a Fourier-like basis.
This leads us to propose a novel, filtered MLMC (F-MLMC) estimator that
mitigates the negative effects of the grid transfer operators on the
multilevel estimator. The introduction of filtering is supported by a
theoretical spectral analysis of the MLMC and F-MLMC estimators in a
simplified setting where the numerical simulators are considered to be linear,
symmetric and circulant operators. Finally, the proposed F-MLMC approach is
successfully tested on a 2D instance of our motivating example.

The remainder of this paper is organized as follows. Section 2 presents the
description of the motivating problem of estimating the normalization
coefficients of a diffusion-based covariance operator. This problem will be
used as an illustration and test case throughout the paper. Section 3
introduces the MLMC estimator of the expectation of discretized random fields,
and illustrates its capabilities and limitations on a one-dimensional (1D)
test case of the normalization problem. Following the conclusions of this
experiment, an F-MLMC estimator is introduced in section 4 and compared to the
unfiltered MLMC estimator on the same 1D test problem. A spectral analysis of
the MLMC and F-MLMC estimators is then carried out in section 5 to investigate
the effects of grid transfer and filtering operators on their variance, in the
specific case of linear, symmetric and circulant numerical simulators. In
section 6, we apply the MLMC and F-MLMC estimators to the 2D problem of
estimating the normalization coefficients of a two-dimensional (2D),
heterogeneous, diffusion-based covariance operator. General conclusions are
drawn in section 7, along with prospective avenues for future work.

##  2 Motivation and problem description

As discussed in the introduction, one area where the estimation of the
expectation of discretized random fields arises is covariance modelling,
specifically when estimating the intrinsic variances of a diffusion-based
covariance operator using a randomization method [58, 59]. This particular
problem, which is outlined in this section, will be our motivating example
throughout the paper.

Let $u\colon\mathcal{D}\to\mathbb{R}$ and $b\colon\mathcal{D}\to\mathbb{R}$ be
square-integrable functions on the domain $\mathcal{D}\subset\mathbb{R}^{d}$
where $d\in\\{1,2,3\\}$ is the spatial dimension. We consider numerical
solutions of the following elliptic equation, subject to application-dependent
boundary conditions (BCs):

$(I-\nabla\cdot\mathbf{K}\nabla)^{m}u=b,$ (1)

where $m$ is a positive integer, $I$ is the identity operator, and
$\mathbf{K}\colon\mathcal{D}\to\mathbb{R}^{d\times d}$ is a symmetric,
positive-definite (SPD) tensor field with entries $[K_{ij}]_{i,j=1,\ldots,d}$.
Equation 1 can be interpreted as a semi-discrete representation of a diffusion
equation integrated over $m$ time-steps where the temporal derivative is
discretized with a backward Euler (implicit) scheme, the time-step is equal to
unity, $\mathbf{K}$ is a diffusivity tensor, and the initial condition is $b$
[56]. If $\mathbf{K}$ is constant then the integral solution on
$\mathbb{R}^{d}$ defines a covariance operator whose kernel is a covariance
function from the Matérn class [60, 25].

We assume that the operator in eq. 1 is discretized in space on a (not
necessarily structured) grid of $n$ cells. We can then deduce the covariance
matrix associated with the numerical solution of eq. 1 as

$\mathbf{L}:=(\mathbf{I}-\bm{\Delta})^{-m}\mathbf{W}^{-1},$ (2)

where $\bm{\Delta}$ is the matrix representing a spatial discretization of the
differential operator ${\nabla\cdot\mathbf{K}\nabla}$, and
${\mathbf{W}\in\mathbb{R}^{n\times n}}$ is an SPD Gram matrix that encodes the
geometrical and structural information related to the discrete approximation
of the diffusive term on the grid. Specifically, $\mathbf{W}$ is such that
$\bm{\Delta}$ is self-adjoint with respect to the inner product whose
weighting matrix is $\mathbf{W}$, i.e.,
$\mathbf{W}\bm{\Delta}=\bm{\Delta}^{\mathrm{T}}\mathbf{W}$. Consequently, the
matrix $\mathbf{L}$ is self-adjoint (symmetric) with respect to the canonical
inner product. In the experiments, we consider only a diagonal diffusivity
tensor $K_{ij}=K_{ij}\delta_{ij}$ where $\delta_{ij}$ is the Kronecker delta.
Specifically, we define the diagonal elements according to the relation
${K_{ii}(\bm{x})=(2m-d-2)^{-1}(D_{ii}(\bm{x}))^{2}}$ where the elements
$[D_{ii}(\bm{x})]_{i=1,\ldots,d}$ correspond to the directional correlation
length-scales at the spatial location $\bm{x}$, and $m>d/2+1$ [56, section 3].

The matrix $\mathbf{L}$ is SPD but does not define a covariance matrix with
meaningful variances for applications like data assimilation. As such,
$\mathbf{L}$ must be normalized by its diagonal so that the desired variances
can be applied. Thus, we define the covariance matrix of interest as
$\mathbf{B}=\bm{\Sigma}\bm{\Gamma}\mathbf{L}\bm{\Gamma}\bm{\Sigma}$, where
$\bm{\Gamma}=\operatorname{Diag}(\operatorname{diag}(\mathbf{L}))^{-1/2}$ is a
normalizing diagonal matrix such that
$\operatorname{Diag}(\operatorname{diag}(\bm{\Gamma}\mathbf{L}\bm{\Gamma}))=\mathbf{I}_{n}$,
and $\bm{\Sigma}^{2}=\bm{\Sigma}\bm{\Sigma}$ is the diagonal matrix with
entries corresponding to the desired variances, i.e.,
$\operatorname{diag}(\mathbf{B})=\operatorname{diag}(\bm{\Sigma}^{2})=(\sigma_{1}^{2},\ldots,\sigma_{n}^{2})$.
In these expressions, the operator $\operatorname{diag}(\cdot)$ maps a matrix
to the vector consisting of the diagonal elements of that matrix, while the
operator $\operatorname{Diag}(\cdot)$ maps a vector to the diagonal matrix
whose diagonal consists of the entries of that vector. For large-scale
problems, the matrix $\mathbf{L}$ is not assembled, and only applications of
$\mathbf{L}$ to vectors are accessible. Thus, its diagonal entries
$\operatorname{diag}(\mathbf{L})=\bm{\theta}=(\theta_{1},\ldots,\theta_{n})$
are not explicitly stored and need to be determined differently [59]. A direct
way would be to recover these by applying $\mathbf{L}$ to the canonical basis
vectors of $\mathbb{R}^{n}$, i.e.,
$\theta_{k}=(\mathbf{L}\mathbf{e}_{k})_{k}$, for $k=1,\ldots,n$. For large
$n$, this approach is not computationally tractable.

An alternative strategy is to approximate $\bm{\theta}$ by randomization.
Taking $m=2q$ and introducing the factorization
${\mathbf{W}=\mathbf{V}\mathbf{V}^{\mathrm{T}}}$, then $\mathbf{L}$ can be
subsequently factored as $\mathbf{L}=\mathbf{A}\mathbf{A}^{\mathrm{T}}$, where

$\mathbf{A}:=(\mathbf{I}-\bm{\Delta})^{-q}(\mathbf{V}^{-1})^{\mathrm{T}},$ (3)

which, as for $\mathbf{L}$ itself, cannot be explicitly assembled in large-
scale applications. This decomposition of $\mathbf{L}$ implies that, for any
random vector $\mathbf{X}$ with $\mathbb{E}[\mathbf{X}]=\mathbf{0}_{n}$ and
$\mathbb{E}[\mathbf{X}\mathbf{X}^{\mathrm{T}}]=\mathbf{I}_{n}$,

$\bm{\theta}=\operatorname{diag}(\mathbf{L})=\operatorname{diag}(\mathbf{A}\mathbb{E}[\mathbf{X}\mathbf{X}^{\mathrm{T}}]\mathbf{A}^{\mathrm{T}})=\operatorname{diag}(\mathbb{C}[\mathbf{A}\mathbf{X}])=\mathbb{V}[\mathbf{A}\mathbf{X}]=\mathbb{E}[\mathbf{A}\mathbf{X}\odot\mathbf{A}\mathbf{X}],$
(4)

where $\odot$ denotes the Schur product (a.k.a. the Hadamard or element-wise
product).

A classical method for estimating the expectation $\bm{\theta}$ of the
$\mathbb{R}^{n}$-valued random vector
${\mathbf{Y}:=\mathbf{A}\mathbf{X}\odot\mathbf{A}\mathbf{X}}$ is through Monte
Carlo (MC) random sampling. Given a random $M$-sample
${\mathcal{X}=(\mathbf{X}^{(i)})_{i=0}^{M}}$ of $\mathbf{X}$, an unbiased
estimator $\hat{\bm{\theta}}$ of $\bm{\theta}$ is the sample mean

$\hat{\bm{\theta}}=\frac{1}{M}\sum_{i=1}^{M}(\mathbf{A}\mathbf{X}^{(i)})\odot(\mathbf{A}\mathbf{X}^{(i)}).$
(5)

We remark that the estimator $\hat{\bm{\theta}}$ only requires $M$
applications of $\mathbf{A}$ to a vector, typically with ${M\ll n}$ in large-
scale applications. Furthermore, by construction, the MC estimator defined by
eq. 5 yields non-negative estimates, which is a fundamental requirement for
the problem under consideration. Alternative MC estimators have been proposed
in the literature for estimating the diagonal of a general (not necessarily
symmetric) matrix [1, 27], which do not employ a factored form of the matrix
and hence do not guarantee non-negative estimates. In all cases, the MC sample
mean estimator is known to converge slowly with respect to the sample size
$M$; specifically, its root-mean-square error (RMSE) is
$\mathcal{O}(M^{-1/2})$ [48]. Here, we investigate the use of the MLMC
methodology [19, 18] to improve (in terms of RMSE) the estimation of
$\bm{\theta}$ and hence the efficiency of determining accurate normalization
coefficients for defining $\bm{\Gamma}$.

##  3 MLMC estimation of the expectation of discretized random fields

In this section, we introduce the MLMC ingredients used in the subsequent
spectral analysis, focusing on the estimation of the expectation of
discretized random fields. The notations introduced here are then used
throughout the remainder of this paper. First, general definitions and
notations are introduced, then, grid transfer operators are presented as well
as their use in the MLMC expectation estimator for discretized fields.
Finally, a 1D illustration showcases the MLMC estimator and its benefits
compared to a standard MC estimator.

The MLMC method aims to reduce the variance, or sampling error, of MC
estimators by combining samples of different fidelities. In favorable cases,
many cheap, low-fidelity samples are used to improve the sampling, while fewer
are required at finer and more expensive resolutions to correct the bias.
Under certain assumptions, [18, Theorem 1] and [9, Theorem 1] ensure that
there exists an allocation of samples on a finite number of levels such that
the computational cost of the MLMC estimator decreases at a faster rate as a
function of the mean square error (MSE), than that of the crude MC estimator.
In practice, MLMC is typically implemented as a sequential algorithm, whose
main idea is to start with a limited number of coarse fidelity levels, and add
as many finer levels as needed to reach a target MSE, with a prescribed
variance/bias balance. In the present work, however, we adopt a multilevel
approach that is closer to multifidelity methods. Specifically, a fine, high-
fidelity level is fixed (and thus so is the bias), while coarser, low-fidelity
levels are considered to reduce the variance, for a prescribed computational
budget.

###  3.1 Preliminary notations and definitions

We first describe the MLMC setting for discretized random fields. Let
$\\{n_{\ell}\\}^{L}_{\ell=0}$ be an increasing sequence of positive integers
corresponding to the size of $L+1$ grids. Those grids of different resolutions
define the hierarchy of levels needed for the MLMC approach. Specifically,
level 0 corresponds to the coarsest level, while level $L$ is the finest. To
each level $\ell$, we associate the SPD Gram matrix
$\mathbf{W}_{\ell}=\mathbf{V}_{\ell}\mathbf{V}_{\ell}^{\mathrm{T}}\in\mathbb{R}^{n_{\ell}\times
n_{\ell}}$ that encodes the structural information related to the
approximation of continuous functions as a vector of $\mathbb{R}^{n_{\ell}}$
on this specific discrete grid. For each level $\ell$, we define the weighted
inner product $\langle\cdot,\cdot\rangle_{\mathbf{W}_{\ell}}$ between elements
of $\mathbb{R}^{n_{\ell}}$ as

$\forall\mathbf{u},\mathbf{v}\in\mathbb{R}^{n_{\ell}},\quad\langle\mathbf{u},\mathbf{v}\rangle_{\mathbf{W}_{\ell}}=\mathbf{u}^{\mathrm{T}}\mathbf{W}_{\ell}\mathbf{v}=\langle\mathbf{V}_{\ell}^{\mathrm{T}}\mathbf{u},\mathbf{V}_{\ell}^{\mathrm{T}}\mathbf{v}\rangle_{\mathbf{I}_{n_{\ell}}},$
(6)

where $\langle\cdot,\cdot\rangle_{\mathbf{I}_{n_{\ell}}}$ denotes the
canonical (Euclidean) dot product between vectors of $\mathbb{R}^{n_{\ell}}$.
The norm induced by $\langle\cdot,\cdot\rangle_{\mathbf{W}_{\ell}}$ is denoted
by $\|\cdot\|_{\mathbf{W}_{\ell}}$. The inner product space
$H_{\ell}:=(\mathbb{R}^{n_{\ell}},\langle\cdot,\cdot\rangle_{\mathbf{W}_{\ell}})$
is a separable Hilbert space.

Let $\mathbf{X}_{L}$ be the $\mathbb{R}^{n_{L}}$-valued random vector
corresponding to a discretized random field at discrete resolution $n_{L}$. We
define abstract numerical models
$f_{\ell}\colon\mathbb{R}^{n_{L}}\to\mathbb{R}^{n_{L}}$ of increasing
resolution (fidelity) $n_{\ell}$ and hence computational cost. We are
interested in estimating the expectation of the output of the finest (highest-
fidelity) model, $\mathbb{E}[f_{L}(\mathbf{X}_{L})]$. The MLMC estimator
$\hat{\bm{\mu}}^{\textnormal{MLMC}}_{L}$ of
$\mathbb{E}[f_{L}(\mathbf{X}_{L})]$ using $L+1$ independent $M_{\ell}$-samples
of $\mathbf{X}_{L}$,
$\\{\mathcal{X}^{(\ell)}_{L}=(\mathbf{X}_{L}^{(\ell,i)})^{M_{\ell}}_{i=1}\\}_{\ell=0}^{L}$,
is defined as

$\hat{\bm{\mu}}^{\textnormal{MLMC}}_{L}=\frac{1}{M_{0}}\sum_{i=1}^{M_{0}}f_{0}(\mathbf{X}_{L}^{(0,i)})+\sum_{\ell=1}^{L}\frac{1}{M_{\ell}}\sum_{i=1}^{M_{\ell}}\big{[}f_{\ell}(\mathbf{X}_{L}^{(\ell,i)})-f_{\ell-1}(\mathbf{X}_{L}^{(\ell,i)})\big{]}.$
(7)

However, in many applications, the different fidelity models are defined on
grids of different resolutions so that their input and output are vectors of
$\mathbb{R}^{n_{\ell}}$, instead of $\mathbb{R}^{n_{L}}$ as required by the
definition of $f_{\ell}$. We denote such abstract numerical models by
$\tilde{f}_{\ell}\colon\mathbb{R}^{n_{\ell}}\to\mathbb{R}^{n_{\ell}}$, for
$\ell=0,\dots,L$, and we are interested in the expectation of the finest
output, $\mathbb{E}[\tilde{f}_{L}(\mathbf{X}_{L})]$.

###  3.2 Grid transfer operators

Because the models $\tilde{f}_{\ell}$ have different domains and codomains
depending on the resolution, it is not possible to use the MLMC estimator eq.
7 directly. Grid transfer operators are needed to transfer the inputs and
outputs between the coarse grids of resolution $n_{\ell}$,
$\ell=0,\ldots,L-1$, and the finest grid of resolution $n_{L}$. In order to
ease the spectral analysis presented in this paper, we restrict ourselves to
linear grid transfer operators, although in principle non-linear transfer
operators may be used.

Operators that transfer a discretized field from a coarse grid onto a finer
grid are classically referred to as prolongation operators. We denote by
$P_{\ell}^{\ell^{\prime}}\colon\mathbb{R}^{n_{\ell}}\to\mathbb{R}^{n_{\ell^{\prime}}}$
the prolongation operator from level $\ell$ to a finer level
$\ell^{\prime}>\ell$ (i.e., with $n_{\ell^{\prime}}>n_{\ell}$). For linear
operators, we may identify $P_{\ell}^{\ell^{\prime}}$ with the matrix
$\mathbf{P}_{\ell}^{\ell^{\prime}}\in\mathbb{R}^{n_{\ell^{\prime}}\times
n_{\ell}}$ such that
$P_{\ell}^{\ell^{\prime}}(\mathbf{x}_{\ell})=\mathbf{P}_{\ell}^{\ell^{\prime}}\mathbf{x}_{\ell}$
for any $\mathbf{x}_{\ell}\in\mathbb{R}^{n_{\ell}}$. Besides, so-called
restriction operators
$R_{\ell^{\prime}}^{\ell}\colon\mathbb{R}^{n_{\ell^{\prime}}}\to\mathbb{R}^{n_{\ell}}$
are used to perform the “reverse” operation of transferring discretized fields
defined on level $\ell^{\prime}$ to a coarser level $\ell<\ell^{\prime}$.
Again, a linear restriction operator $R_{\ell^{\prime}}^{\ell}$ may be
identified with the appropriate matrix
$\mathbf{R}_{\ell^{\prime}}^{\ell}\in\mathbb{R}^{n_{\ell}\times
n_{\ell^{\prime}}}$. We may now define $f_{\ell}$ from $\tilde{f}_{\ell}$ as

$f_{\ell}=P_{\ell}^{L}\circ\tilde{f}_{\ell}\circ R_{L}^{\ell},$ (8)

for $\ell=0,\ldots,L$, so that the MLMC expression eq. 7 can be applied. These
numerical models involve transfer operators between a given level and the
finest level $L$. In practice, these operators may be defined from transfer
operators between successive levels, as

$\mathbf{P}_{\ell}^{L}:=\mathbf{P}^{L}_{L-1}\cdots\mathbf{P}^{\ell+2}_{\ell+1}\mathbf{P}^{\ell+1}_{\ell}\mbox{
and
}\mathbf{R}^{\ell}_{L}:=\mathbf{R}^{\ell}_{\ell+1}\cdots\mathbf{R}^{L-2}_{L-1}\mathbf{R}^{L-1}_{L},\quad\forall\ell=0,\ldots,L-1,$
(9)

and $\mathbf{R}_{L}^{L}=\mathbf{P}_{L}^{L}=\mathbf{I}_{n_{L}}$, so that
$\tilde{f}_{L}=f_{L}$.

###  3.3 The MLMC estimator

Considering a separable Hilbert space $H$ equipped with the inner product
$\langle\cdot,\cdot\rangle_{H}$ and induced norm $\|\cdot\|_{H}$, the space of
second-order $H$-valued random variables,

$L^{2}(\Omega,H):=\\{\xi\colon\Omega\to
H\mid\textstyle\int_{\Omega}\|\xi(\omega)\|^{2}_{H}\,\mathrm{d}\mathbb{P}(\omega)<+\infty\\},$
(10)

is a Hilbert space when equipped with the inner product
$\langle\cdot,\cdot\rangle_{L^{2}(\Omega,H)}$ defined as

$\forall\xi,\eta\in
L^{2}(\Omega,H),\quad\langle\xi,\eta\rangle_{L^{2}(\Omega,H)}:=\int_{\Omega}\langle\xi(\omega),\eta(\omega)\rangle_{H}\,\mathrm{d}\mathbb{P}(\omega).$
(11)

The induced norm is denoted by $\|\cdot\|_{L^{2}(\Omega,H)}$, and, for
convenience, we will from now on use the shorthand notation
$\|\cdot\|_{L^{2}(\Omega,H)}=\mathbb{E}[\|\cdot\|^{2}_{H}]^{1/2}$. In this
paper, we are particularly interested in the spaces $L^{2}(\Omega,H_{\ell})$,
for $\ell=0,\ldots,L$, with $H_{\ell}$ defined as in section 3.1. In what
follows, we assume that $\mathbf{X}_{L}\in L^{2}(\Omega,H_{L})$ and that, for
any $\mathbf{X}\in L^{2}(\Omega,H_{\ell})$,

  * •

$\tilde{f}_{\ell}(\mathbf{X})\in L^{2}(\Omega,H_{\ell})$, for
$\ell=0,\ldots,L$;

  * •

$P_{\ell-1}^{\ell}(\mathbf{X})\in L^{2}(\Omega,H_{\ell})$, for
$\ell=1,\ldots,L$;

  * •

$R_{\ell+1}^{\ell}(\mathbf{X})\in L^{2}(\Omega,H_{\ell})$, for
$\ell=0,\ldots,L-1$.

It follows that $\mathbf{Y}_{\ell}:=f_{\ell}(\mathbf{X}_{L})\in
L^{2}(\Omega,H_{L})$ for $\ell=0,\ldots,L$. A convergence analysis of the MLMC
method in such Hilbert spaces was proposed in [2], from which we reuse here
certain definitions and properties.

In particular, we are interested in this work in the normwise MSE of the MLMC
estimator with respect to the exact expectation $\bm{\mu}$, defined by

$\operatorname{MSE}(\hat{\bm{\mu}}^{\textnormal{MLMC}}_{L},\bm{\mu}):=\|\hat{\bm{\mu}}^{\textnormal{MLMC}}_{L}-\bm{\mu}\|^{2}_{L^{2}(\Omega,H_{L})}.$
(12)

As shown in [2, Theorem 3.1]), the MSE admits the decomposition

$\operatorname{MSE}(\hat{\bm{\mu}}^{\textnormal{MLMC}}_{L},\bm{\mu})=\mathcal{V}(\hat{\bm{\mu}}^{\textnormal{MLMC}}_{L})+\|\mathbb{E}[\mathbf{Y}_{L}]-\bm{\mu}\|^{2}_{\mathbf{W}_{L}},$
(13)

where

$\forall\mathbf{Y}\in
L^{2}(\Omega,H_{L}),\quad\mathcal{V}(\mathbf{Y}):=\|\mathbf{Y}-\mathbb{E}[\mathbf{Y}]\|^{2}_{L^{2}(\Omega,H_{L})}=\mathbb{E}[\|\mathbf{Y}\|^{2}_{\mathbf{W}_{L}}]-\|\mathbb{E}[\mathbf{Y}]\|^{2}_{\mathbf{W}_{L}}.$
(14)

The first term in eq. 13 is referred to as the variance of the multilevel
estimator, while the second term corresponds to the squared bias. We note that
the multilevel estimator is unbiased with respect to the expectation of the
output at the finest level $\bm{\mu}=\mathbb{E}[\mathbf{Y}_{L}]$. Furthermore,
the variance part of the MSE,
$\mathcal{V}(\hat{\bm{\mu}}^{\textnormal{MLMC}}_{L})$, can be further
decomposed level-wise into

$\mathcal{V}(\hat{\bm{\mu}}^{\textnormal{MLMC}}_{L})=\frac{1}{M_{0}}\mathcal{V}(\mathbf{Y}_{0})+\sum^{L}_{\ell=1}\frac{1}{M_{\ell}}\mathcal{V}(\mathbf{Y}_{\ell}-\mathbf{Y}_{\ell-1}).$
(15)

###  3.4 1D illustration

The variance estimation problem described in section 2 is considered here in a
1D setting. We consider cell-centered discretizations of the diffusion
operator on the domain $\mathcal{D}:=[0,1]$, with periodic boundary
conditions. Specifically, we consider uniform grids composed of $n_{\ell}$
cells of size $n_{\ell}^{-1}$, with the associated Gram matrices
$\mathbf{W}_{\ell}=n_{\ell}^{-1}\mathbf{I}_{n_{\ell}}$, corresponding to a
piecewise constant approximation of the solution. It follows immediately that
$\mathbf{V}_{\ell}=n_{\ell}^{-1/2}\mathbf{I}_{n_{\ell}}=\mathbf{W}_{\ell}^{1/2}$.
We define hierarchies of different depth, corresponding to
$L\in\\{0,\ldots,5\\}$, with a fixed finest discretization corresponding to
$n_{L}=512$. For $\ell=0,\dots,L-1$, a constant refinement factor
$n_{\ell+1}/n_{\ell}=2$ is defined, so that $n_{\ell}=2^{\ell-L}n_{L}$. The
inter-level prolongation and restriction operators are defined as

$\mathbf{P}_{\ell-1}^{\ell}:=\begin{bmatrix}1&&\\\ 1&&\\\ &1&\\\ &1&\\\
&&\ddots\end{bmatrix}\in\mathbb{R}^{n_{\ell}\times
n_{\ell-1}},\qquad\mathbf{R}^{\ell-1}_{\ell}:=\mathbf{V}_{\ell-1}^{-1}(\mathbf{P}_{\ell-1}^{\ell})^{\mathrm{T}}\mathbf{V}_{\ell}=\frac{1}{\sqrt{2}}(\mathbf{P}_{\ell-1}^{\ell})^{\mathrm{T}},$
(16)

respectively, and grid transfer operators between an arbitrary level and the
finest level are defined as in eq. 9. This particular choice of grid transfer
operators implies that

$\forall\mathbf{x}_{\ell}\in\mathbb{R}^{n_{\ell}},\;\|\mathbf{P}_{\ell}^{\ell^{\prime}}\mathbf{x}_{\ell}\|_{\mathbf{W}_{\ell^{\prime}}}=\|\mathbf{x}_{\ell}\|_{\mathbf{W}_{\ell}},\quad\mbox{and
}\forall\mathbf{X}_{\ell^{\prime}}\sim\mathcal{N}(\mathbf{0}_{n_{\ell^{\prime}}},\mathbf{I}_{n_{\ell^{\prime}}}),\;\mathbf{R}_{\ell^{\prime}}^{\ell}\mathbf{X}_{\ell^{\prime}}\sim\mathcal{N}(\mathbf{0}_{n_{\ell}},\mathbf{I}_{n_{\ell}}),$
(17)

for $\ell^{\prime}>\ell$. Further detail on the motivation for this choice is
given in appendix A.

The diffusion tensor field $\mathbf{K}$ (which actually reduces to a scalar
field in 1D) is taken to be constant, by setting ${D_{11}(\bm{x})=D}$ for all
${\bm{x}\in\mathcal{D}}$ (see section 2). The scalar value $D\in\mathbb{R}$
will be referred to as the length-scale and, unless stated otherwise, it will
be set to $D=0.06\approx 30n_{L}^{-1}$. Furthermore, we choose a fixed value
of ${m=10}$, implying that ${q=5}$. As described in section 2, the statistic
of interest is
$\mathbb{E}[f_{L}(\mathbf{X}_{L})]=\mathbb{E}[(\mathbf{A}_{L}\mathbf{X}_{L})\odot(\mathbf{A}_{L}\mathbf{X}_{L})]$,
and where $\mathbf{X}_{L}$ is a standard normal random vector, i.e.,
$\mathbf{X}_{L}\sim\mathcal{N}(\mathbf{0}_{n_{L}},\mathbf{I}_{n_{L}})$. The
abstract numerical models $f_{\ell}$ are defined by eq. 8, with
$\tilde{f}_{\ell}\colon\mathbf{X}_{\ell}\mapsto(\mathbf{A}_{\ell}\mathbf{X}_{\ell})\odot(\mathbf{A}_{\ell}\mathbf{X}_{\ell})$,
and $\mathbf{A}_{\ell}$ defined as in eq. 3 on the appropriate grid of level
$\ell$. For this experiment, we rely on the optimal sample allocation that
minimizes, for a given computational budget $\eta$, the variance of the
multilevel estimator [41],

$M_{\ell}=\left\lfloor\frac{\eta}{\mathcal{S}_{L}}\sqrt{\frac{\mathcal{V}_{\ell}}{\mathcal{C}_{\ell}+\mathcal{C}_{\ell-1}}}\right\rfloor^{+},\qquad\mathcal{V}_{\ell}:=\mathcal{V}(\mathbf{Y}_{\ell}-\mathbf{Y}_{\ell-1}),\qquad\mathcal{S}_{L}:=\sum_{\ell=0}^{L}\sqrt{\mathcal{V}_{\ell}(\mathcal{C}_{\ell}+\mathcal{C}_{\ell-1})},$
(18)

where $\mathcal{C}_{\ell}$ denotes the computational cost of one evaluation of
$f_{\ell}$ and $\lfloor\cdot\rfloor^{+}:=\max(1,\lfloor\cdot\rfloor)$, where
$\lfloor\cdot\rfloor$ denotes the floor function. It should be noted that, by
convention, the quantities indexed by $\ell=-1$ vanish, so that
$\mathcal{C}_{-1}=0$ and $\mathcal{V}_{0}=\mathcal{V}(\mathbf{Y}_{0})$. In
addition, we remark that, with eq. 18, the prescribed computational budget
$\eta$ may actually be exceeded by an additional cost lower than
$\mathcal{C}_{L}$. However, this extra cost should be negligible provided
$\eta\gg\mathcal{C}_{L}$. In what follows, we assume a linear computational
cost model, specifically $\mathcal{C}_{\ell}=\mathcal{O}(n_{\ell})$,
consistent with a fixed number of sparse (banded) matrix-vector applications
required for the evaluation of $f_{\ell}$. The variances $\mathcal{V}_{\ell}$
are estimated in a preprocessing stage with a pilot sample of size $1000$, and
the multilevel estimation is conducted with a computational budget
$\eta=100\mathcal{C}_{L}$.

(a) Sample allocation $\\{M_{\ell}\\}_{\ell=0}^{L}$ across levels.

(b) Total variance $\mathcal{V}(\hat{\bm{\mu}}^{\textnormal{MLMC}}_{L})$.

Figure 1: Optimal sample allocation across levels and resulting total variance
of the MC estimator ($L=0$) and different MLMC estimators ($L\in\\{1,3,5\\}$),
for the estimation problem described in section 2, with length-scale $D=0.06$.
The finest level $L$ always corresponds to a discretization with $n_{L}=512$
cells, and the total budget is set to $\eta=100\mathcal{C}_{L}$. The variance
is estimated from $1000$ estimators.

The sample allocation obtained with these parameters and the total variance of
different MLMC estimators is shown in fig. 1(a). Instead of requiring 100
fine-grid evaluations like a single-level Monte Carlo estimator, the 2-, 4-
and 6-level MLMC estimators only require a few evaluations on the finest
level, while more than 900 evaluations are used on the coarsest grid for the
6-level MLMC. The optimal allocation allows for MLMC estimators to reach a
lower total variance than the crude MC estimator, as shown in fig. 1(b). The
minimum variance is achieved with a 4-level MLMC ($L=3$), with a variance
reduction of about 70% compared to MC ($L=0$). The addition of two coarser
levels with $32$ and $16$ cells, respectively, actually deteriorates the total
variance.

The results above suggest that the MLMC estimators in eq. 7 have a lower
variance than the single-level MC estimator. However, the variance of the
estimator, according to the definition eq. 14, stands here in a normwise
(integral) sense, and we cannot infer any information for a specific scale
(frequency). In fact, the variance of the estimator associated with a given
frequency could provide valuable information for the analysis of multilevel
estimators involving discretized fields. In particular, the small-scale (high-
frequency) components of a signal on the fine grid cannot be represented on
coarser grids.

To study the error at different scales, we decompose the MSE in the Hartley
basis [28, 5, 4]. The Hartley basis is a Fourier-like basis, commonly used in
circulant embedding techniques for generating stationary Gaussian random
fields with prescribed covariance structure [23, 22, 31]. Its main advantage
is that it consists of purely real basis vectors, as opposed to the Fourier
basis (whose basis vectors are complex), thus easing interpretation and
visualization. On level $\ell$, the $n_{\ell}$ cell-centered Hartley basis
vectors $\\{\mathbf{h}^{\ell}_{k}\\}_{k=0}^{n_{\ell}-1}$ correspond to the
columns of the Hartley matrix $\mathbf{H}_{\ell}$ with entries

$(\mathbf{H}_{\ell})_{j,k}:=\frac{1}{\sqrt{n_{\ell}}}\left(\cos{\frac{2(j+\frac{1}{2})k\pi}{n_{\ell}}}+\sin{\frac{2(j+\frac{1}{2})k\pi}{n_{\ell}}}\right),\quad\forall
j,k=0,\dots,n_{\ell}-1.$ (19)

The matrices $\mathbf{H}_{\ell}$ are orthogonal, i.e., for any
$\ell=0,\dots,L$,
$\mathbf{H}_{\ell}^{\mathrm{T}}\mathbf{H}_{\ell}=\mathbf{H}_{\ell}\mathbf{H}_{\ell}^{\mathrm{T}}=\mathbf{I}_{n_{\ell}}$
(see appendix B). The MSE eq. 12 can thus be decomposed into contributions of
the individual Hartley modes,

$\operatorname{MSE}(\hat{\bm{\mu}}^{\textnormal{MLMC}}_{L},\bm{\mu})=\mathbb{E}[\|\hat{\bm{\mu}}^{\textnormal{MLMC}}_{L}-\bm{\mu}\|_{\mathbf{W}_{L}}^{2}]=\mathbb{E}[\|\mathbf{H}_{L}^{\mathrm{T}}\mathbf{W}_{L}^{1/2}(\hat{\bm{\mu}}^{\textnormal{MLMC}}_{L}-\bm{\mu})\|_{\mathbf{I}_{n_{L}}}^{2}],$
(20)

which, exploiting the linearity of the expectation operator and the fact that
$\mathbf{W}_{L}=n_{L}^{-1}\mathbf{I}_{n_{L}}$, may be compactly written as

$\operatorname{MSE}(\hat{\bm{\mu}}^{\textnormal{MLMC}}_{L},\bm{\mu})=\|\bm{\nu}\|_{1}:=\sum_{k=0}^{n_{L}-1}\nu_{k},\qquad\nu_{k}:=n_{L}^{-1}\mathbb{E}[((\mathbf{h}^{L}_{k})^{\mathrm{T}}(\hat{\bm{\mu}}^{\textnormal{MLMC}}_{L}-\bm{\mu}))^{2}],$
(21)

where $\bm{\nu}=(\nu_{k})_{k=0}^{n_{L}-1}$ denotes the spectral MSE. In what
follows, we take $\bm{\mu}=\mathbb{E}[\mathbf{Y}_{L}]$ as a reference, so that
the MSE coincides with the variance
$\mathcal{V}(\hat{\bm{\mu}}^{\textnormal{MLMC}}_{L})$, and
$\nu_{k}=n_{L}^{-1}\mathbb{V}[(\mathbf{h}^{L}_{k})^{\mathrm{T}}\hat{\bm{\mu}}^{\textnormal{MLMC}}_{L}]$.
We thus refer to $\bm{\nu}$ as the spectral variance of the estimator.
Furthermore, we define the cumulative variance,
$\bm{\nu}^{\textnormal{cml}}=(\nu^{\textnormal{cml}}_{k})_{k=0}^{n_{L}-1}$,
such that
$\nu^{\textnormal{cml}}_{k}=\sum_{k^{\prime}=0}^{k}\nu_{k^{\prime}}$, implying
that the total variance is given by
$\nu^{\textnormal{cml}}_{n_{L}-1}=\|\bm{\nu}\|_{1}$.

In fig. 2, we plot the spectral variance $\bm{\nu}$, as well as the
corresponding cumulative variance $\bm{\nu}^{\textnormal{cml}}$, associated
with the 2-, 4- and 6-level MLMC estimators ($L\in\\{1,3,5\\}$) and with the
single-level MC estimator ($L=0$), for the estimation problem described in
section 2, with length-scale $D=0.06$. Note that, in this and the following
figures, the Hartley basis vectors
$\\{\mathbf{h}^{\ell}_{k}\\}_{k=0}^{n_{\ell}-1}$ are reordered by increasing
representable frequency (see appendix C). Figure 2(a) shows that, for the
single-level MC estimator, most of the error arises from the first Hartley
modes, which are associated with the large scales, or low frequencies, of the
discretized field. This is confirmed in fig. 2(b), where the cumulative
variance of the MC estimator rapidly increases before reaching a plateau,
showing that the variance is concentrated on the first few modes. The spectral
variances of the MLMC estimators exhibit a similar behavior, but its decay in
the high frequencies is less pronounced. For the 2-level MLMC estimator, the
spectral variance starts with a similar decay as that of the MC estimator in
the low frequencies, before increasing again in the high frequencies.
Nevertheless, the variance is still concentrated on the first few modes, and
the low-frequency components of the spectral variance are lower than those of
the single-level MC estimator. This translates into a lower plateau reached by
the cumulative variance, as shown in Figure 2(b). Moreover, the variance
increase in the high-frequencies translates into a noticeable increase in the
cumulative variance in the last few Hartley modes, which, in turn, results in
a non-negligible increase in the total variance. As more levels are added, the
spectral variance significantly deteriorates in the high-frequencies. While
for the 4-level MLMC estimator, this deterioration is compensated by a lower
variance in the low frequencies, this is no longer the case for the 6-level
MLMC estimator, whose cumulative variance eventually gets larger than that of
the 4-level MLMC estimator, thus resulting in a larger total variance, as
already evidenced in fig. 1.

(a) Spectral variance $\bm{\nu}$.

(b) Cumulative variance $\bm{\nu}^{\textnormal{cml}}$.

Figure 2: Spectral and cumulative variance of the MC estimator ($L=0$) and
different MLMC estimators ($L\in\\{1,3,5\\}$), for the estimation problem
described in section 2, with length-scale $D=0.06$. The finest level $L$
always corresponds to a discretization with $n_{L}=512$ cells, and the total
budget is set to $\eta=100\mathcal{C}_{L}$. The variance is estimated from
$1000$ estimators.

It should be noted that the length-scale $D=0.06$ used for the experiments
presented above induces output fields mostly composed of large scales (low
frequencies). Decreasing its value increases the frequencies of the output
field, thus introducing smaller scales. Figure 3 shows the spectral and
cumulative variance of the MC and MLMC estimators for a reduced length-scale
$D=0.01\approx 5n_{L}^{-1}$. We observe that the decay of the spectral
variance of the MC estimator as the frequency increases is slower than for
$D=0.06$, so that the variance is concentrated on a wider low-frequency range.
For the multilevel estimators, we see that the significant deterioration of
the variance in the high-frequencies is no longer compensated by a better
estimation in the low frequencies, except for the 2-level estimator, which
remains slightly better than the single-level estimator in terms of total
variance.

(a) Spectral variance $\bm{\nu}$.

(b) Cumulative variance $\bm{\nu}^{\textnormal{cml}}$.

Figure 3: Same as fig. 2 but with length-scale $D=0.01$.

##  4 The filtered MLMC method

The spectral variance of the MLMC estimator for the 1D variance estimation
problem shows a significant deterioration of the variance in the high
frequencies. Inspired by multigrid methods [55, 7, 61], we propose an
improvement of the MLMC estimator for discretized random fields by adding pre-
and post-filtering (or, in multigrid terminology, smoothing). The objective is
to filter out the smaller scales, which cannot be represented on the coarse
grids, before using a restriction operator and after using a prolongation
operator.

###  4.1 Filtered grid transfer operators

Filtering the small-scale components out of a signal on level $\ell$ is
achieved using a low-pass filtering operator
$S_{\ell}\colon\mathbb{R}^{n_{\ell}}\to\mathbb{R}^{n_{\ell}}$. In what
follows, we resort to linear filtering operators $S_{\ell}$, which can thus be
identified with matrices $\mathbf{S}_{\ell}$. Filtered grid transfer operators
$\bar{P}_{\ell}^{L}$ and $\bar{R}_{L}^{\ell}$ are defined through their
corresponding matrices $\bar{\mathbf{P}}_{\ell}^{L}$ and
$\bar{\mathbf{R}}^{\ell}_{L}$ by

$\bar{\mathbf{P}}_{\ell}^{L}:=\bar{\mathbf{P}}^{L}_{L-1}\cdots\bar{\mathbf{P}}^{\ell+2}_{\ell+1}\bar{\mathbf{P}}^{\ell+1}_{\ell}\mbox{
and
}\bar{\mathbf{R}}^{\ell}_{L}:=\bar{\mathbf{R}}^{\ell}_{\ell+1}\cdots\bar{\mathbf{R}}^{L-2}_{L-1}\bar{\mathbf{R}}^{L-1}_{L},\quad\forall\ell=0,\ldots,L-1$
(22)

where

$\bar{\mathbf{R}}_{\ell}^{\ell-1}:=\mathbf{R}_{\ell}^{\ell-1}\mathbf{S}_{\ell},\quad\text{and}\quad\bar{\mathbf{P}}_{\ell-1}^{\ell}:=\mathbf{S}_{\ell}\mathbf{P}_{\ell-1}^{\ell}\quad\forall\ell=1,\ldots,L,$
(23)

and $\bar{\mathbf{R}}_{L}^{L}=\bar{\mathbf{P}}_{L}^{L}=\mathbf{I}_{n_{L}}$.
The filtered MLMC estimator, hereafter referred to as the F-MLMC estimator,
then reads

$\hat{\bm{\mu}}^{\textnormal{F-MLMC}}_{L}=\frac{1}{M_{0}}\sum_{i=1}^{M_{0}}\bar{f}_{0}(\mathbf{X}_{L}^{(0,i)})+\sum_{\ell=1}^{L}\frac{1}{M_{\ell}}\sum_{i=1}^{M_{\ell}}[\bar{f}_{\ell}(\mathbf{X}_{L}^{(\ell,i)})-\bar{f}_{\ell-1}(\mathbf{X}_{L}^{(\ell,i)})],$
(24)

where
$\bar{f}_{\ell}=\bar{P}_{\ell}^{L}\circ\tilde{f}_{\ell}\circ\bar{R}_{L}^{\ell}$,
for $\ell=0,\dots,L$. In what follows, we resort to the second-order Shapiro
filter [52, 16, 17] defined as

$\mathbf{S}_{\ell}:=\frac{1}{4}\begin{bmatrix}2&1&&&1\\\ 1&2&1&&\\\
&\ddots&\ddots&\ddots&\\\ &&1&2&1\\\ 1&&&1&2\end{bmatrix}.$ (25)

With the specific grid transfer operators defined in eq. 16, the operator
$\bar{\mathbf{P}}_{\ell-1}^{\ell}$ corresponds to the linear interpolation
operator between the levels $\ell-1$ and $\ell$.

The idea of smoothing highly oscillatory signals within the MLMC framework was
also recently exploited in [31], in a setting where the input of the
considered simulator is a discretized random field, but where the output is a
scalar quantity. Consequently, prolongation operators are not needed, and only
restriction operations are required. In fact, because the Cartesian, node-
centered, finite element discretization considered in [31] produces nested
meshes, the restriction operator reduces to a selection operator, which
selects the vector entries associated with the desired nodal unknowns. The
pre-smoothing (or pre-filtering) operation proposed in [31] is based on a
spectral truncation of the high-frequency components of the input signal, and
designed for each MLMC level such that the truncation error matches the
discretization error, so that the MLMC properties are best exploited. The
combination of this pre-filtering technique with spectral post-filtering,
along the lines of [13, section 4.4], may be investigated in future work.

###  4.2 1D illustration

(a) Sample allocation $\\{M_{\ell}\\}_{\ell=0}^{L}$ across levels.

(b) Total variance $\mathcal{V}(\hat{\bm{\mu}}^{\textnormal{MLMC}}_{L})$ and
$\mathcal{V}(\hat{\bm{\mu}}^{\textnormal{F-MLMC}}_{L})$.

Figure 4: Optimal sample allocation across levels for the 6-level ($L=5$) MLMC
and F-MLMC estimators, and total variance of the MC estimator ($L=0$) and
different MLMC and F-MLMC estimators ($L=1,\dots,5$), for the estimation
problem described in section 2, with length-scale $D=0.06$. The finest level
$L$ always corresponds to a discretization with $n_{L}=512$ cells, and the
total budget is set to $\eta=100\mathcal{C}_{L}$. The variance is estimated
from $1000$ estimators.

The 1D variance estimation problem described in section 2 and used in section
3 is again used to illustrate the effects of filtering. The optimal sample
allocation, obtained from eq. 18, is reported in fig. 4(a) for the 6-level
($L=5$) MLMC and F-MLMC estimators, for a length-scale $D=0.06$. We observe
that the filtered estimator allocates smaller sample sizes to the finer grids
than the unfiltered one, while a more than 50% larger sample size is allocated
to the coarsest level. This may be explained by the fact that, in the absence
of filtering, more effort needs to be made on the finer levels, due to the
large high-frequency components of the variance. This change in the optimal
allocation has a significant impact on the total variance, as shown by fig.
4(b). The F-MLMC estimators achieve a lower total variance than their
unfiltered counterpart regardless of the hierarchy, with the 6-level F-MLMC
($L=5$) achieving the lowest variance. The addition of filters allows the two
coarsest grids to reduce the total variance of the estimator further, while
the variance of the unfiltered MLMC estimator deteriorates when including
these levels.

(a) Spectral variance $\bm{\nu}$.

(b) Cumulative variance $\bm{\nu}^{\textnormal{cml}}$.

Figure 5: Spectral and cumulative variance of the MC estimator ($L=0$) and of
different F-MLMC estimators ($L\in\\{1,3,5\\}$), for the estimation problem
described in section 2, with length-scale $D=0.06$. The finest level $L$
always corresponds to a discretization with $n_{L}=512$ cells, and the total
budget is set to $\eta=100\mathcal{C}_{L}$. The variance is estimated from
1000 estimators.

Figure 5 shows the spectral variance $\bm{\nu}$ and the cumulative variance
$\bm{\nu}^{\textnormal{cml}}$ of the F-MLMC estimators with different grid
hierarchies corresponding to $L\in\\{1,\ldots,3\\}$. The effects of the
filters are especially visible on the spectral variance (fig. 5(a)), which is
significantly reduced, not only in the high frequencies, but also in the
lowest ones (i.e., corresponding to the first few Hartley modes), as compared
to the unfiltered MLMC (fig. 2(a)). The reduced error in the high frequencies
(small scales) allows the cumulative variance not to be impacted in the last
few Hartley modes (fig. 5(b)), as opposed to the unfiltered MLMC (fig. 2(b)).
Furthermore, the reduced error in the low frequencies (large scales)
translates into a lower plateau of the cumulative variance, hence a lower
total variance, than for the unfiltered MLMC. Specifically, the addition of
filters leads to a 90% reduction in total variance of the F-MLMC estimator
compared to the single-level MC estimator, corresponding to a 70% reduction
compared to the best, 4-level unfiltered MLMC estimator. The same quantities
are plotted in fig. 6 for a length-scale $D=0.01$. Again, the addition of
filters improves the multilevel estimation in both the low and the high
frequencies, which in turn benefits the cumulative and thus the total
variance. Filtering is here even more beneficial than for $D=0.06$, in the
sense that the F-MLMC estimator has a significantly lower variance than the
single-level MC estimator, which the unfiltered MLMC estimator failed to
achieve. Specifically, the variance of the 2-level F-MLMC estimator is reduced
by about 40% compared to the single-level MC estimator, and that of the
4-level F-MLMC is reduced by more than 50%, while its unfiltered counterpart
actually produced a higher variance than the single-level MC estimator.

(a) Spectral variance $\bm{\nu}$.

(b) Cumulative variance $\bm{\nu}^{\textnormal{cml}}$.

Figure 6: Same as Figure 5 but with length-scale $D=0.01$.

These experiments highlight that the effect of spurious high frequencies
caused by grid transfer operations is detrimental to the MLMC estimator. The
addition of pre- and post-filtering operations to mitigate these effects is
necessary for the multilevel estimator to reach its full potential.

##  5 Spectral analysis

In this section, we conduct a spectral analysis of the MLMC estimator eq. 7
and F-MLMC estimator eq. 24 to study more closely the effects of grid transfer
operators and filters on the variance of the multilevel estimator at different
scales. The analysis is conducted in a 1D setting similar to that of section
3.4, with Gram matrices
$\mathbf{W}_{\ell}=n_{\ell}^{-1}\mathbf{I}_{n_{\ell}}$, for $\ell=0,\dots,L$,
corresponding to a cell-centered, piecewise constant approximation of the
solution on uniform grids. The setting is further simplified to the case where
the numerical simulators are linear, and that their associated matrices are
symmetric and circulant.

###  5.1 MLMC with linear simulators

We consider linear numerical simulators of the form
$\tilde{f}_{\ell}\colon\mathbf{x}_{\ell}\mapsto\tilde{\mathbf{F}}_{\ell}\mathbf{x}_{\ell}$,
where $\tilde{\mathbf{F}}_{\ell}\in\mathbb{R}^{n_{\ell}\times n_{\ell}}$, so
that, for a given random vector $\mathbf{X}_{L}$ corresponding to a
discretized random field on the finest level $L$,

$\mathbf{Y}_{\ell}=\mathbf{F}_{\ell}\mathbf{X}_{L},\qquad\mathbf{F}_{\ell}:=\mathbf{P}^{L}_{\ell}\tilde{\mathbf{F}}_{\ell}\mathbf{R}^{\ell}_{L}\in\mathbb{R}^{n_{L}\times
n_{L}},\quad\forall\ell=0,\dots,L,$ (26)

with $\mathbf{P}^{L}_{L}=\mathbf{R}^{L}_{L}=\mathbf{I}_{n_{L}}$. For now, we
make no further assumption on $\tilde{\mathbf{F}}_{\ell}$. Exploiting the
linearity of the simulators, and letting
$\mathbf{G}:=\mathbb{E}[(\mathbf{X}_{L}-\mathbb{E}(\mathbf{X}_{L})(\mathbf{X}_{L}-\mathbb{E}(\mathbf{X}_{L})^{\mathrm{T}}]$
be the covariance matrix of $\mathbf{X}_{L}$, the variance eq. 15 of the
multilevel estimator reduces to

$\mathcal{V}(\hat{\bm{\mu}}^{\textnormal{MLMC}}_{L})=\frac{1}{M_{0}}\|\mathbf{F}_{0}\mathbf{G}^{1/2}\|^{2}_{F,\mathbf{W}_{L}}+\sum^{L}_{\ell=1}\frac{1}{M_{\ell}}\|(\mathbf{F}_{\ell}-\mathbf{F}_{\ell-1})\mathbf{G}^{1/2}\|^{2}_{F,\mathbf{W}_{L}},$
(27)

where
$\|\cdot\|_{F,\mathbf{W}}\colon\mathbf{M}\mapsto\operatorname{tr}(\mathbf{M}^{\mathrm{T}}\mathbf{W}\mathbf{M})^{1/2}$
denotes the $\mathbf{W}$-weighted Frobenius norm, for any SPD weighting matrix
$\mathbf{W}$ (see appendix D for the derivation of eq. 27). The variance in
eq. 27 emphasizes that the variance reduction is closely related to the
similarity of successive fidelity models. Again, the orthogonality of the
Hartley matrix $\mathbf{H}_{L}$ can be invoked to decompose the variance into
contributions of the individual Hartley modes,

$\mathcal{V}(\hat{\bm{\mu}}^{\textnormal{MLMC}}_{L})=\frac{1}{M_{0}}\|\mathbf{H}_{L}^{\mathrm{T}}\mathbf{F}_{0}\mathbf{G}^{1/2}\|^{2}_{F,\mathbf{W}_{L}}+\sum^{L}_{\ell=1}\frac{1}{M_{\ell}}\|\mathbf{H}_{L}^{\mathrm{T}}(\mathbf{F}_{\ell}-\mathbf{F}_{\ell-1})\mathbf{G}^{1/2}\|^{2}_{F,\mathbf{W}_{L}}.$
(28)

###  5.2 Two-level MLMC with linear, symmetric, circulant simulators

For the remainder of this analysis, we focus on 2-level MLMC estimators
($L=1$), with $n_{0}$ cells on the coarse grid and $n_{1}=2n_{0}$ cells on the
fine grid. Figure 7 depicts the basis vectors $\mathbf{h}^{1}_{k}$ and
$\mathbf{h}^{0}_{k}$ of the fine- and coarse-grid Hartley bases, discretized
on grids with $n_{1}=16$ and $n_{0}=8$ cells, respectively. These plots
highlight that, because of aliasing, the basis vectors exhibit a discrete
frequency that is different from their continuous counterpart. Specifically,
for $\ell\in\\{0,1\\}$, the vectors indexed by $k$ close to 0 or $n_{\ell}-1$
are discrete signals with low frequency, while their frequency increases as
$k$ tends to $n_{\ell}/2$.

Figure 7: The vectors $\mathbf{h}^{0}_{k}$ (orange) and $\mathbf{h}^{1}_{k}$
(blue) of a coarse-grid Hartley basis $\mathbf{H}_{0}$ with $n_{0}=8$ and of a
fine-grid Hartley basis with $n_{1}=16$. They are represented on cell-centered
grids of size 8 and 16 that discretize the domain $[0,1]$.

We now examine the effects of the grid transfer operators on the Hartley basis
vectors. Such effects have been studied extensively for multigrid methods [7,
61, 55] using different bases. We succinctly present here results for the
specific, cell-centered Hartley basis eq. 19. With the prolongation operator
$\mathbf{P}:=\mathbf{P}_{0}^{1}$ defined in eq. 16, we have (see appendix E)

$\mathbf{P}\mathbf{h}^{0}_{k}=\sqrt{2}(c_{k}\mathbf{h}^{1}_{k}-c_{n_{0}+k}\mathbf{h}^{1}_{n_{0}+k}),\quad\forall
k=0,\dots,n_{0}-1,$ (29)

where the coefficients $c_{k}:=\cos(k\pi/(2n_{0}))$ are strictly decreasing
with $k=0,\dots,2n_{0}-1$, i.e.,

$1=c_{0}>c_{1}>\dots>(c_{n_{0}}=0)>\dots>c_{2n_{0}-1}>-1.$ (30)

Prolongating a coarse-grid basis vector $\mathbf{h}^{0}_{k}$, with
$k=0,\dots,n_{0}-1$, produces a (fine-grid) vector consisting of a linear
combination of two fine-grid basis vectors $\mathbf{h}^{1}_{k}$ and
$\mathbf{h}^{1}_{n_{0}+k}$. Figure 7 shows that, for $k\leq n_{0}/2$, the
fine-grid signal $\mathbf{h}^{1}_{k}$ has the same frequency as the original,
coarse-grid signal $\mathbf{h}^{0}_{k}$, while $\mathbf{h}^{1}_{n_{0}+k}$ has
higher frequency. Conversely, for $k>n_{0}/2$, the fine-grid signal
$\mathbf{h}^{1}_{n_{0}+k}$ has the same frequency as the coarse-grid signal
$\mathbf{h}^{0}_{k}$, while $\mathbf{h}^{1}_{k}$ has higher frequency. In both
cases, the prolongation of a coarse-grid signal introduces spurious high-
frequency (i.e., small-scale) components to the prolongated signal.
Fortunately, both $\mathbf{h}^{1}_{k}$ and $\mathbf{h}^{1}_{n_{0}+k}$ are
damped by a factor which is closer to zero for the spurious, high-frequency
signals than for the consistent, low-frequency signals. Specifically, $c_{k}$
tends to $1$ and $c_{n_{0}+k}$ tends to $0$ as $k$ tends to $0$, thus damping
more severely the spurious, high-frequency signals $\mathbf{h}^{1}_{n_{0}+k}$
than the consistent, low-frequency signals $\mathbf{h}^{1}_{k}$. Conversely,
$c_{k}$ tends to $0$ and $c_{n_{0}+k}$ tends to $1$ as $k$ tends to $n_{0}$,
thus damping more severely the spurious, high-frequency signals
$\mathbf{h}^{1}_{k}$ than the consistent, low-frequency signals
$\mathbf{h}^{1}_{n_{0}+k}$.

For the restriction operator $\mathbf{R}:=\mathbf{R}_{1}^{0}$ defined in eq.
16, we have (see appendix F)

$\mathbf{R}\mathbf{h}^{1}_{k}=c_{k}\mathbf{h}^{0}_{k},\quad\mbox{and
}\mathbf{R}\mathbf{h}^{1}_{n_{0}+k}=-c_{n_{0}+k}\mathbf{h}^{0}_{k},\quad\forall
k=0,\dots,n_{0}-1.$ (31)

Restricting a fine-grid basis vector $\mathbf{h}^{1}_{k}$, with
$k=0,\dots,n_{0}-1$, produces a (coarse-grid) vector proportional to the
corresponding coarse-grid basis vector $\mathbf{h}^{0}_{k}$, specifically,
reduced by a factor $c_{k}\leq 1$. For $n_{0}/2<k<n_{0}$, the restricted
signal has lower frequency than the original, fine-grid signal, as illustrates
fig. 7. Similarly, restricting a fine-grid basis vector $\mathbf{h}^{1}_{k}$,
with $k=n_{0},\dots,2n_{0}-1$, produces a (coarse-grid) vector proportional to
the complementary coarse-grid basis vector $\mathbf{h}^{0}_{k-n_{0}}$,
specifically reduced by a factor $-c_{k}<1$. Again, high-frequency fine-grid
signals $\mathbf{h}^{1}_{k}$ corresponding to $n_{0}\leq k<3n_{0}/2$ are
restricted to a signal with lower frequency. In conclusion, high-frequency
fine-grid basis vectors that cannot be represented on the coarse grid are thus
restricted to lower-frequency signals. Fortunately, such signals are the most
damped, since they correspond to ranges of $k$ where $c_{k}$ is closer to $0$.

######  Remark 1.

Denoting
$\mathbf{C}:=\begin{bmatrix}\operatorname{Diag}(\\{c_{k}\\}_{k=0}^{n_{0}-1})&\operatorname{Diag}(\\{-c_{n_{0}+k}\\}_{k=0}^{n_{0}-1})\end{bmatrix}\in\mathbb{R}^{n_{0}\times
2n_{0}}$, the identities eqs. 29 and 31 may be compactly recast as

$\displaystyle\mathbf{P}\mathbf{H}_{0}$
$\displaystyle=\sqrt{2}\mathbf{H}_{1}\mathbf{C}^{\mathrm{T}},$
$\displaystyle\mathbf{R}\mathbf{H}_{1}$
$\displaystyle=\mathbf{H}_{0}\mathbf{C},$
$\displaystyle\mathbf{H}_{1}^{\mathrm{T}}\mathbf{P}$
$\displaystyle=\sqrt{2}\mathbf{C}^{\mathrm{T}}\mathbf{H}_{0}^{\mathrm{T}}$
$\displaystyle\mathbf{H}_{0}^{\mathrm{T}}\mathbf{R}$
$\displaystyle=\mathbf{C}\mathbf{H}_{1}^{\mathrm{T}},$ (32)

where the last two identities follow from the first two by exploiting the
orthogonality of $\mathbf{H}_{0}$ and $\mathbf{H}_{1}$.

We now assume that the operators $\tilde{\mathbf{F}}_{\ell}$ are symmetric,
circulant matrices. Such matrices can be diagonalized in the Hartley basis [4]
(see appendix G for the proof with the cell-centered basis), i.e.,

$\tilde{\mathbf{F}}_{\ell}=\mathbf{H}_{\ell}\bm{\Lambda}_{\ell}\mathbf{H}_{\ell}^{\mathrm{T}},\quad\bm{\Lambda}_{\ell}:=\operatorname{Diag}(\\{\lambda^{\ell}_{k}\\}_{k=0}^{n_{L}-1}),\quad\mbox{for
}\ell\in\\{0,1\\}.$ (33)

This property, along with the identities in eq. 32, allows
$\mathbf{F}_{0}:=\mathbf{P}\tilde{\mathbf{F}}_{0}\mathbf{R}=\mathbf{P}\mathbf{H}_{0}\bm{\Lambda}_{0}\mathbf{H}_{0}^{\mathrm{T}}\mathbf{R}$
to be decomposed as
${\mathbf{F}_{0}=\mathbf{H}_{1}\mathbf{M}\mathbf{H}_{1}^{\mathrm{T}}}$, where
$\mathbf{M}:=\mathbf{C}^{\mathrm{T}}\bm{\Lambda}_{0}\mathbf{C}=\begin{bmatrix}\mathbf{M}_{11}&\mathbf{M}_{12}\\\
\mathbf{M}_{21}&\mathbf{M}_{22}\end{bmatrix}$, with

$\displaystyle\mathbf{M}_{11}$
$\displaystyle:=\sqrt{2}\operatorname{Diag}(\\{c_{k}^{2}\lambda_{k}^{0}\\}_{k=0}^{n_{0}-1}),$
(34) $\displaystyle\mathbf{M}_{22}$
$\displaystyle:=\sqrt{2}\operatorname{Diag}(\\{c_{n_{0}+k}^{2}\lambda_{k}^{0}\\}_{k=0}^{n_{0}-1}),$
(35) $\displaystyle\mathbf{M}_{12}$
$\displaystyle:=\sqrt{2}\operatorname{Diag}(\\{-c_{k}c_{n_{0}+k}\lambda_{k}^{0}\\}_{k=0}^{n_{0}-1})=\mathbf{M}_{21}.$
(36)

As a consequence, eq. 28 becomes

$\mathcal{V}(\hat{\bm{\mu}}^{\textnormal{MLMC}}_{1})=\frac{1}{M_{0}}\|\mathbf{M}\mathbf{H}_{1}^{\mathrm{T}}\mathbf{G}^{1/2}\|^{2}_{F,\mathbf{W}_{1}}+\frac{1}{M_{1}}\|(\bm{\Lambda}_{1}-\mathbf{M})\mathbf{H}_{1}^{\mathrm{T}}\mathbf{G}^{1/2}\|^{2}_{F,\mathbf{W}_{1}}.$
(37)

To reduce the variance of the correction term in eq. 37, the difference
between $\mathbf{M}$ and $\bm{\Lambda}_{1}$ needs to be as small as possible.
First, we note that the two off-diagonal blocks $\mathbf{M}_{12}$ and
$\mathbf{M}_{21}$, which are themselves diagonal matrices, contribute to
increasing this difference. On the main diagonal, i.e., in the diagonal blocks
$\mathbf{M}_{11}$ and $\mathbf{M}_{22}$, scaled eigenvalues of
$\tilde{\mathbf{F}}_{0}$ appear twice. To compare these diagonal blocks to the
eigenvalues of $\tilde{\mathbf{F}}_{1}$ in $\bm{\Lambda}_{1}$, further
assumptions on $\tilde{\mathbf{F}}_{0}$ are required. We thus introduce the
Galerkin coarse-grid operator, which is an algebraic way of constructing the
coarse-grid operator $\tilde{\mathbf{F}}_{0}$ from the fine-grid operator
$\tilde{\mathbf{F}}_{1}$, and which is widely used in multigrid methods and
their analysis [55]. Specifically, the Galerkin operator is defined as

$\tilde{\mathbf{F}}_{0}:=\frac{1}{2}\mathbf{R}\tilde{\mathbf{F}}_{1}\mathbf{P}\in\mathbb{R}^{n_{0}\times
n_{0}}.$ (38)

This Galerkin operator is the optimal operator in terms of minimizing
$\|\tilde{\mathbf{F}}_{1}-\mathbf{P}\tilde{\mathbf{F}}_{0}\mathbf{R}\|_{F,\mathbf{W}_{1}}^{2}$
for $\mathbf{W}_{1}=n_{1}^{-1}\mathbf{I}_{n_{1}}$ and the grid transfer
operators defined by eq. 16 (see appendix H). It follows from eqs. 38 and 32
that the Galerkin operator $\tilde{\mathbf{F}}_{0}$ can be diagonalized in the
Hartley basis as
$\tilde{\mathbf{F}}_{0}=\mathbf{H}_{0}\bm{\Lambda}_{0}\mathbf{H}_{0}^{\mathrm{T}}$
with
$\bm{\Lambda}_{0}=\frac{1}{\sqrt{2}}\mathbf{C}\bm{\Lambda}_{1}\mathbf{C}^{\mathrm{T}}$.
In other words, the eigenvalues $\\{\lambda_{k}^{0}\\}_{k=0}^{n_{0}-1}$ of the
Galerkin operator $\tilde{\mathbf{F}}_{0}$ can be expressed from the
eigenvalues $\\{\lambda_{k}^{1}\\}_{k=0}^{n_{1}-1}$ of
$\tilde{\mathbf{F}}_{1}$,

$\lambda_{k}^{0}=\frac{1}{\sqrt{2}}\left(c_{k}^{2}\lambda_{k}^{1}+c_{n_{0}+k}^{2}\lambda_{n_{0}+k}^{1}\right),\quad\forall
k=0,\dots,n_{0}-1.$ (39)

The resulting blocks of $\mathbf{M}$ then read

$\displaystyle\mathbf{M}_{11}$
$\displaystyle=\operatorname{Diag}(\\{c_{k}^{4}\lambda_{k}^{1}+c_{k}^{2}c_{n_{0}+k}^{2}\lambda_{n_{0}+k}^{1}\\}_{k=0}^{n_{0}-1}),$
(40) $\displaystyle\mathbf{M}_{22}$
$\displaystyle=\operatorname{Diag}(\\{c_{n_{0}+k}^{4}\lambda_{n_{0}+k}^{1}+c_{k}^{2}c_{n_{0}+k}^{2}\lambda_{k}^{1}\\}_{k=0}^{n_{0}-1})=\operatorname{Diag}(\\{c_{k}^{4}\lambda_{k}^{1}+c_{k-n_{0}}^{2}c_{k}^{2}\lambda_{k-n_{0}}^{1}\\}_{k=n_{0}}^{2n_{0}-1}),$
(41) $\displaystyle\mathbf{M}_{12}$
$\displaystyle=\operatorname{Diag}(\\{-c_{k}^{3}c_{n_{0}+k}\lambda_{k}^{1}-c_{k}c_{n_{0}+k}^{3}\lambda_{n_{0}+k}^{1}\\}_{k=0}^{n_{0}-1})=\mathbf{M}_{21}.$
(42)

Using elementary trigonometric identities, we remark that
$c_{k-n_{0}}=c_{n_{0}+k-2n_{0}}=-c_{n_{0}+k}$. Therefore, the main diagonal of
$\mathbf{M}$ deviates from $\bm{\Lambda}_{1}$ by a multiplicative damping
factor $c_{k}^{4}\leq 1$ for $k=0,\ldots,2n_{0}-1$, on the one hand, and by
the addition of a spurious, complementary eigenvalue, though also damped by
$c_{k}^{2}c_{n_{0}+k}^{2}<1$ for $k=0,\ldots,2n_{0}-1$. The off-diagonal
blocks $\mathbf{M}_{12}$ and $\mathbf{M}_{21}$ introduce spurious terms that
contribute to the difference $\bm{\Lambda}_{1}-\mathbf{M}$ and thus increase
the variance of the 2-level estimator. Note that, because these terms are to
be compared with $0$, the comparison with the eigenvalues of the fine-grid
operator $\tilde{\mathbf{F}}_{1}$ are of little interest, which is why we
simply consider the damping factors $-\sqrt{2}c_{k}c_{n_{0}+k}$ with respect
to eigenvalues $\lambda^{0}_{k}$ of the coarse-grid operator
$\tilde{\mathbf{F}}_{0}$, for $k=0,\ldots,n_{0}-1$, given by eq. 36.

(a) Damping factors in $\mathbf{M}$.

(b) Damping factors in $\bar{\mathbf{M}}$.

Figure 8: Damping factors of the eigenvalues in $\mathbf{M}$ and
$\bar{\mathbf{M}}$ as functions of $k=0,\dots,n_{1}-1$ with $n_{1}=2n_{0}=32$,
when using the Galerkin operator eq. 38. The black curves correspond to the
factors of the correct eigenvalues on the main diagonal, while blue curves
represent the factors of the spurious eigenvalues. The orange curves
correspond to the factors of the off-diagonal blocks.

The evolution with $k=0,\ldots,n_{1}-1$ of the three damping factors is
presented in fig. 8(a). We observe that the eigenvalues $\lambda^{1}_{k}$
associated with fine, low-frequency Hartley basis vectors, i.e., for $k$ close
to $0$ and $n_{1}-1$ (see fig. 7), are well-represented on the main diagonal
of $\mathbf{M}$, since $c_{k}^{4}\approx 1$. At the same time, the damping
factors $c_{k}^{2}c_{n_{0}+k}^{2}$ corresponding to spurious components are
close to $0$, resulting in small values of the first and last few diagonal
entries of the difference $\bm{\Lambda}_{1}-\mathbf{M}$. In particular, we
note that the first entry of $\mathbf{M}_{11}$ exactly matches that of
$\bm{\Lambda}_{1}$, i.e., $\lambda^{1}_{0}$, because $c_{0}=1$ and
$c_{n_{0}}=0$. On the other hand, eigenvalues $\lambda^{1}_{k}$ associated
with fine, medium- to high-frequency Hartley basis vectors are severely damped
since $c_{k}^{4}$ quickly decreases to $0$ as $k$ approaches $n_{0}=n_{1}/2$.
Spurious diagonal components are also somewhat damped by a factor
$c_{k}^{2}c_{n_{0}+k}^{2}$, which is maximal for $k\in\\{n_{1}/4,3n_{1}/4\\}$.
Finally, the spurious, off-diagonal components are damped by a factor
$-\sqrt{2}c_{k}c_{n_{0}+k}$, which is maximal for
$k\in\\{n_{1}/4,3n_{1}/4\\}$.

###  5.3 Two-level F-MLMC with linear, symmetric, circulant simulators

We now consider the filtered grid transfer operators defined in section 4.1,
$\bar{\mathbf{P}}=\mathbf{S}_{1}\mathbf{P}$ and
$\bar{\mathbf{R}}=\mathbf{R}\mathbf{S}_{1}$, where $\mathbf{S}_{1}$ denotes
the second-order Shapiro filter defined in eq. 25. Similarly to eqs. 29 and
31, it is possible to study the effect the filtered transfer operators on the
Hartley basis vectors (see appendix I for the derivations). For the
prolongation $\bar{\mathbf{P}}$, we have

$\bar{\mathbf{P}}\mathbf{h}^{0}_{k}=\sqrt{2}(c_{k}^{3}\mathbf{h}^{1}_{k}-c_{n_{0}+k}^{3}\mathbf{h}^{1}_{n_{0}+k}),\quad\forall
k=0,\dots,n_{0}-1.$ (43)

The addition of the Shapiro filter raises the damping factors $c_{k}$ to the
power 3. The prolongated Hartley basis vectors are thus more severely damped
than they were without filtering. Again, the most damped fine-grid basis
vectors $\mathbf{h}^{1}_{k}$ are those corresponding to $k$ close to
$n_{0}=n_{1}/2$, i.e., high-frequency signals. For the restriction operator we
have

$\bar{\mathbf{R}}\mathbf{h}^{1}_{k}=c_{k}^{3}\mathbf{h}^{0}_{k},\quad\mbox{and
}\bar{\mathbf{R}}\mathbf{h}^{1}_{n_{0}+k}=-c_{n_{0}+k}^{3}\mathbf{h}^{0}_{k},\quad\forall
k=0,\dots,n_{0}-1.$ (44)

Similar conclusions can be drawn as for the unfiltered case, but with the
damping factors raised to the power of 3, thus increasing the damping, which
still affects more strongly the fine-grid, high-frequency signals that cannot
be represented on the coarse grid. Similarly to the unfiltered case (see
remark 1, identities eqs. 43 and 44 can be recast as

$\displaystyle\bar{\mathbf{P}}\mathbf{H}_{0}$
$\displaystyle=\sqrt{2}\mathbf{H}_{1}\mathbf{C}_{3}^{\mathrm{T}},$
$\displaystyle\bar{\mathbf{R}}\mathbf{H}_{1}$
$\displaystyle=\mathbf{H}_{0}\mathbf{C}_{3},$
$\displaystyle\mathbf{H}_{1}^{\mathrm{T}}\bar{\mathbf{P}}$
$\displaystyle=\sqrt{2}\mathbf{C}_{3}^{\mathrm{T}}\mathbf{H}_{0}^{\mathrm{T}}$
$\displaystyle\mathbf{H}_{0}^{\mathrm{T}}\bar{\mathbf{R}}$
$\displaystyle=\mathbf{C}_{3}\mathbf{H}_{1}^{\mathrm{T}},$ (45)

where
$\mathbf{C}_{3}:=\begin{bmatrix}\operatorname{Diag}(\\{c_{k}^{3}\\}_{k=0}^{n_{0}-1})&\operatorname{Diag}(\\{-c_{n_{0}+k}^{3}\\}_{k=0}^{n_{0}-1})\end{bmatrix}\in\mathbb{R}^{n_{0}\times
2n_{0}}$.

The impact of filtering on the total variance of the MLMC estimator is now
assessed, considering a 2-level MLMC estimator and assuming that
$\tilde{\mathbf{F}}_{1}$ and $\tilde{\mathbf{F}}_{0}$ are symmetric, circulant
matrices. From eqs. 33 and 45, we deduce the decomposition
$\bar{\mathbf{F}}_{0}:=\bar{\mathbf{P}}\tilde{\mathbf{F}}_{0}\bar{\mathbf{R}}=\mathbf{H}_{1}\bar{\mathbf{M}}\mathbf{H}_{1}^{\mathrm{T}}$,
where
$\bar{\mathbf{M}}:=\begin{bmatrix}\bar{\mathbf{M}}_{11}&\bar{\mathbf{M}}_{12}\\\
\bar{\mathbf{M}}_{21}&\bar{\mathbf{M}}_{22}\end{bmatrix}$, with

$\displaystyle\bar{\mathbf{M}}_{11}$
$\displaystyle=\sqrt{2}\operatorname{Diag}(\\{c_{k}^{6}\lambda_{k}^{0}\\}_{k=0}^{n_{0}-1}),$
(46) $\displaystyle\bar{\mathbf{M}}_{22}$
$\displaystyle=\sqrt{2}\operatorname{Diag}(\\{c_{n_{0}+k}^{6}\lambda_{k}^{0}\\}_{k=0}^{n_{0}-1}),$
(47) $\displaystyle\bar{\mathbf{M}}_{12}$
$\displaystyle=\sqrt{2}\operatorname{Diag}(\\{-c_{k}^{3}c_{n_{0}+k}^{3}\lambda_{k}^{0}\\}_{k=0}^{n_{0}-1})=\bar{\mathbf{M}}_{2,1}.$
(48)

Upon replacing $\mathbf{F}_{0}$ with its filtered counterpart
$\bar{\mathbf{F}}_{0}$, eq. 28 can be written as in eq. 37, but with
$\mathbf{M}$ replaced by $\bar{\mathbf{M}}$. The sparsity pattern of
$\bar{\mathbf{M}}$ is identical to that of $\mathbf{M}$, and its entries are
similar, but with damping factors raised to increased powers. We remark that
the off-diagonal blocks $\bar{\mathbf{M}}_{12}$ and $\bar{\mathbf{M}}_{21}$
have entries that are more strongly damped than those of their unfiltered
counterparts, $\mathbf{M}_{12}$ and $\mathbf{M}_{21}$, as can be visualized in
fig. 8 (orange plots), which contributes to reducing the off-diagonal entries
of $\bm{\Lambda}_{1}-\bar{\mathbf{M}}$ compared to those of
$\bm{\Lambda}_{1}-\mathbf{M}$.

To study more closely the diagonal entries of $\bar{\mathbf{M}}$ and compare
them to the eigenvalues of $\tilde{\mathbf{F}}_{1}$ in $\bm{\Lambda}_{1}$ we
resort to the Galerkin operator defined by
$\tilde{\mathbf{F}}_{0}:=\frac{1}{2}\bar{\mathbf{R}}\tilde{\mathbf{F}}_{1}\bar{\mathbf{P}}$.
This definition is inspired by the form of eq. 38, although there is no
guarantee of its optimality. Then, from eqs. 33 and 45, it follows that the
Galerkin operator $\tilde{\mathbf{F}}_{0}$ can be diagonalized in the Hartley
basis as
$\tilde{\mathbf{F}}_{0}=\mathbf{H}_{0}\bm{\Lambda}_{0}\mathbf{H}_{0}^{\mathrm{T}}$
with
$\bm{\Lambda}_{0}=\frac{1}{\sqrt{2}}\mathbf{C}_{3}\bm{\Lambda}_{1}\mathbf{C}_{3}^{\mathrm{T}}$,
or, equivalently,

$\lambda_{k}^{0}=\frac{1}{\sqrt{2}}\left(c_{k}^{6}\lambda_{k}^{1}+c_{n_{0}+k}^{6}\lambda_{n_{0}+k}^{1}\right),\quad\forall
k=0,\dots,n_{0}-1.$ (49)

The diagonal blocks of $\bar{\mathbf{M}}$ thus become

$\displaystyle\bar{\mathbf{M}}_{11}$
$\displaystyle=\operatorname{Diag}(\\{c_{k}^{12}\lambda_{k}^{1}+c_{n_{0}+k}^{6}c_{k}^{6}\lambda_{n_{0}+k}^{1}\\}_{k=0}^{n_{0}-1}),$
(50) $\displaystyle\bar{\mathbf{M}}_{22}$
$\displaystyle=\operatorname{Diag}(\\{c_{n_{0}+k}^{12}\lambda_{n_{0}+k}^{1}+c_{n_{0}+k}^{6}c_{k}^{6}\lambda_{k}^{1}\\}_{k=0}^{n_{0}-1})=\operatorname{Diag}(\\{c_{k}^{12}\lambda_{k}^{1}+c_{k-n_{0}}^{6}c_{k}^{6}\lambda_{k-n_{0}}^{1}\\}_{k=n_{0}}^{2n_{0}-1}).$
(51)

The conclusions are the same as for the diagonal blocks of $\mathbf{M}$, but
with damping factors raised to higher powers. Consequently, the spurious
eigenvalues on the main diagonal are damped more strongly than in the
unfiltered case, as can be visualized in fig. 8 (blue plots). Moreover, the
factors $c_{k}^{8}$ also damp more strongly the consistent eigenvalues
compared to the factors $c_{k}^{4}$ of $\mathbf{M}$, as can be visualized in
fig. 8 (black plots). Although the addition of filters induces a certain loss
of (mainly high-frequency) information, the spurious signals introduced by the
grid transfer operators are significantly reduced. These results suggest that,
for filtering to be beneficial, a tradeoff needs to be found between reducing
the detrimental effects induced by the grid transfer operators and degrading
of the information of the original signal. The second-order Shapiro filters
considered in this paper seem to offer a good compromise. Further endeavors to
study and improve the filtering process may be pursued in future work.

##  6 2D application

In this section, we again apply the filtered MLMC estimator to the variance
estimation problem described in section 2, but now considering a 2D diffusion
operator and using parameter settings more complex than those in the 1D
example presented earlier. Specifically, the diffusivity field $\mathbf{K}$ is
now specified to be non-uniform, making the variance field also non-uniform.
The considered domain is $\mathcal{D}=(0,1)\times(0,2)\subset\mathbb{R}^{2}$.
The boundary conditions are chosen periodic along both directions. The 2D
diffusivity tensor field $\mathbf{K}$ is chosen to be diagonal and
heterogeneous

$\mathbf{K}=\frac{1}{2m-4}\begin{bmatrix}D_{11}^{2}&0\\\
0&D_{22}^{2}\end{bmatrix},$ (52)

where, for $i=1,2$, $D_{ii}\colon\mathcal{D}\to\mathbb{R}$ represents a
length-scale field in the $i$-th direction. As previously, we let $m=2q=10$.
We model $D_{11}=\zeta(\omega_{1})$ and $D_{22}=\zeta(\omega_{2})$ as two
different realizations of a 2D, periodic Gaussian random field $\zeta$ over
$\mathcal{D}$ of uniform mean $\mu_{\zeta}$, and of quasi-Gaussian covariance
structure with uniform variance
$\sigma_{\zeta}^{2}=(\mu_{\zeta}/5)^{2}=0.04\mu_{\zeta}^{2}$ and uniform
length-scale $D_{\zeta}$. In the following experiments, two sets of parameters
are considered for $\zeta$, namely $(\mu_{\zeta}=0.12,D_{\zeta}=0.2)$ and
$(\mu_{\zeta}=0.02,D_{\zeta}=0.04)$. The corresponding realizations used in
the subsequent experiments are depicted in fig. 9.

(a) $\mu_{\zeta}=0.12$ and $D_{\zeta}=0.2$.

(b) $\mu_{\zeta}=0.02$ and $D_{\zeta}=0.04$.

Figure 9: Length-scale fields $D_{11}$ (top) and $D_{22}$ (bottom) used in the
reported experiments for the definition of $\mathbf{K}$ in eq. 52, with
$\mu_{\zeta}=0.12$ and $D_{\zeta}=0.2$ (left), and with $\mu_{\zeta}=0.02$ and
$D_{\zeta}=0.04$ (right).

The diffusion operator in eq. 1 is discretized on Cartesian grids of size
$n_{\ell}=n_{\ell}^{x}\times n_{\ell}^{y}$, with the length-scale fields
$D_{ij}$ discretized at edge centers, and the solution discretized at cell
centers. The finest grid considered here is composed of $n_{L}=256\times 128$
cells. Three coarser grids are used with a uniform coarsening factor of 4,
i.e., $n_{L-1}=128\times 64$, $n_{L-2}=64\times 32$ and $n_{L-3}=32\times 16$.
The resulting Gram matrix on level $\ell$ is given by
$\mathbf{W}_{\ell}=2n_{\ell}^{-1}\mathbf{I}_{n_{\ell}}$. The discrete
operators $\mathbf{A}_{\ell}\in\mathbb{R}^{n_{\ell}\times n_{\ell}}$ are
designed to apply to and return vectors of size $n_{\ell}$ whose entries are
associated with cell centers that are sorted by increasing $y$-coordinate
first, then by increasing $x$-coordinate. In other words, the entry indexed by
$k=in_{\ell}^{x}+j$ in such vectors is associated with a cell center located
at $(x_{j},y_{i})\in\mathbb{R}^{2}$, where $x_{j}:=(j+1/2)/n_{\ell}^{x}$ and
$y_{i}:=(i+1/2)/n_{\ell}^{y}$, for $j=1,\ldots,n_{\ell}^{x}$ and
$i=1,\ldots,n_{\ell}^{y}$. As in the 1D illustration, we are interested in the
multilevel estimation of the fine discretized field
$\bm{\theta}=\mathbb{E}[(\mathbf{A}_{L}\mathbf{X}_{L})\odot(\mathbf{A}_{L}\mathbf{X}_{L})]$.
With the Cartesian ordering of the unknowns described above, the 2D
prolongation and restriction operators may be constructed as the Kronecker
product of their 1D counterpart, defined in eq. 16, in the $x$ and $y$
directions. The 2D grid transfer operators defined in this manner are then
related by
$\mathbf{R}_{\ell}^{\ell-1}=\frac{1}{2}(\mathbf{P}_{\ell-1}^{\ell})^{\mathrm{T}}$,
for $\ell=1,\dots,L$. Likewise, the second-order 2D Shapiro filter is defined
as the Kronecker product of two 1D, second-order Shapiro filters defined in
eq. 25. The cost model is still considered linear in the number of cells,
implying here that one simulator evaluation on level $\ell$ is as
computationally expensive as 4 simulator evaluations on level $\ell-1$. The
baseline is a crude, single-level MC ($L=0$) computed with a sample of size
100; hence the total budget is $\eta=100\mathcal{C}_{L}$.

Figure 10: The top-left sub-figure shows the exact discretized field
$\bm{\theta}$ on level $L$, while the other sub-figures depict the expectation
of the single-level MC estimator (top-right), the expectation of the 2- and
4-level MLMC estimators (middle-left and middle-right), and the expectation of
the 2- and 4-level F-MLMC estimators (bottom-left and bottom-right). The
tensor field $\mathbf{K}$ corresponds to $\mu_{\zeta}=0.12$ and
$D_{\zeta}=0.2$ (fig. 9(a)). The expectation is approximated from 500
estimators, each constructed with a computational budget
$\eta=100\mathcal{C}_{L}$.

For a given diffusivity tensor field $\mathbf{K}$, the exact discretized field
$\bm{\theta}$ on the finest level $L$ is computed explicitly as
$\theta_{k}=(\mathbf{L}\mathbf{e}_{k})_{k}$, for $k=1,\ldots,n_{L}$, where
$\mathbf{L}$ denotes the discrete diffusion operator defined in eq. 2 on level
$L$ and $\mathbf{e}_{k}$ denotes the $k$-th canonical basis vector of
$\mathbb{R}^{n_{L}}$. The top-left sub-figure of fig. 10 presents the obtained
field $\bm{\theta}$ for a fixed tensor field $\mathbf{K}$ obtained as
described above with parameters $\mu_{\zeta}=0.12$ and $D_{\zeta}=0.2$. The
other sub-figures represent the expectation of the compared estimators, namely
the single-level MC estimator, and the 2- and 4-level (F-)MLMC estimators. The
expectation is approximated from 500 estimators, each constructed with a
computational budget $\eta=100\mathcal{C}_{L}$. These figures confirm a key
property of the MC and (F-)MLMC estimators, namely that they are unbiased.
Indeed, their expectations visually coincide (up to statistical error due to
the estimation) with the reference, i.e.,
$\mathbb{E}[\hat{\bm{\theta}}]=\bm{\theta}$. Consequently, the MSE solely
consists of the variance of the estimators.

Figure 11: Variance of the MC estimator (top), of the 2- and 4-level MLMC
estimators (middle-left and middle-right), and of the 2- and 4-level F-MLMC
estimators (bottom-left and bottom-right). The tensor field $\mathbf{K}$
corresponds to $\mu_{\zeta}=0.12$ and $D_{\zeta}=0.2$ (fig. 9(a)). The
variance is approximated from 500 estimators, each constructed with a
computational budget $\eta=100\mathcal{C}_{L}$.

Figure 11 shows the variance of the considered estimators. Although the
variance of the MLMC estimators (middle row) visually seems lower than that of
the MC estimator (top row), we clearly observe high-frequency fluctuations in
the variance field, both for the 2-level (left) and the 4-level (right)
unfiltered MLMC estimators. These high-frequency components of the variance
may have significant consequences on individual estimations of $\bm{\theta}$.
Indeed, while the MLMC estimations will, on average, match the desired field
(owing to the unbiasedness of the estimators), individual estimations will be
polluted by high-frequency error components because of the large high-
frequency components of the variance. The bottom row of fig. 11 demonstrates,
at least visually, that the addition of filtering effectively damps these
high-frequency components. These observations are confirmed by fig. 12, which
shows the spectral decomposition of the variance in the Hartley space,
allowing for the visualization of the contribution of each scale (or
frequency) to the variance. The Hartley matrix $\mathbf{H}_{L}$ used to
project the 2D variance (discretized) fields onto the Hartley spectral space
is defined as the Kronecker product of two 1D Hartley matrices defined in eq.
19. The 2D spectral variance is defined as
$\bm{\nu}=(\nu_{k})_{k=0}^{n_{L}-1}\in\mathbb{R}^{n_{L}}$, with

$\nu_{k}:=2n_{L}^{-1}\mathbb{E}[(\mathbf{h}^{L}_{k})^{\mathrm{T}}(\hat{\bm{\theta}}^{\textnormal{MLMC}}_{L}-\bm{\theta}))^{2}],\quad
k=0,\dots,n_{L}-1,$ (53)

where $\mathbf{h}^{L}_{k}$ denotes the $k$-th column of $\mathbf{H}_{L}$.
Again, the columns of $\mathbf{H}_{L}$ are re-ordered so that, in fig. 12, the
frequencies in the $x$ and $y$ directions increase along the associated axes,
starting from the lower-left corner, corresponding to low frequencies (large
scales). We observe that the variance of the unfiltered MLMC estimators
(middle row) clearly exhibits larger high-frequency components than the
single-level MC estimator (top row), with values of the same order of
magnitude as the low-frequency components. This not only translates into
noticeable high-frequency fluctuations of the variance field, as evidenced in
fig. 11, but it may also deteriorate the overall variance, as was the case in
some instances in the 1D illustration (see section 4.2).

Figure 12: Spectral variance $\bm{\nu}$ of the MC estimator (top), of the 2-
and 4-level MLMC estimators (middle-left and middle-right), and of the 2- and
4-level F-MLMC estimators (bottom-left and bottom-right). The tensor field
$\mathbf{K}$ corresponds to $\mu_{\zeta}=0.12$ and $D_{\zeta}=0.2$ (fig.
9(a)). The variance is approximated from 500 estimators, each constructed with
a computational budget $\eta=100\mathcal{C}_{L}$.

(a) $\mathbf{K}$ with $\mu_{\zeta}=0.12$ and $D_{\zeta}=0.2$ (fig. 9(a)).

(b) $\mathbf{K}$ with $\mu_{\zeta}=0.02$ and $D_{\zeta}=0.04$ (fig. 9(b)).

Figure 14: Cumulative variance of the MC estimator ($L=0$) and of different
MLMC estimators (left) and F-MLMC estimators (right) with $L\in\\{1,2,3\\}$.
The tensor field $\mathbf{K}$ corresponds to $\mu_{\zeta}=0.12$ and
$D_{\zeta}=0.2$ (top), and $\mu_{\zeta}=0.02$ and $D_{\zeta}=0.04$ (bottom).
The variance is approximated from 500 estimators, each constructed with a
computational budget $\eta=100\mathcal{C}_{L}$.

With the Cartesian indexing and the re-ordering of the Hartley basis vectors
described earlier, the 2D cumulative variance
$\bm{\nu}^{\textnormal{cml}}=(\nu^{\textnormal{cml}}_{k})_{k=0}^{n_{L}^{y}-1}\in\mathbb{R}^{n_{L}^{y}}$
is computed by adding the components of the 2D spectral variance eq. 53 shown
in fig. 12 in rectangle patterns starting with the bottom-left corner.
Specifically,

$\nu^{\textnormal{cml}}_{k}=\sum_{i=0}^{k}\sum_{j=0}^{2k+1}\nu_{in_{L}^{x}+j},\quad
k=0,\dots,n_{L}^{y}-1.$ (54)

As such, the first components of $\bm{\nu}^{\textnormal{cml}}$ represent the
cumulative variance associated with the low frequencies (large scales), while
the last component $\nu^{\textnormal{cml}}_{n_{L}^{y}-1}$ coincides with the
total variance of the estimator. Figure 14 presents the cumulative variance of
the MLMC and F-MLMC estimators for $\mathbf{K}$ corresponding to
$(\mu_{\zeta}=0.12,D_{\zeta}=0.2)$ and $(\mu_{\zeta}=0.02,D_{\zeta}=0.04)$. In
the first case (fig. 14(a)), for both the MC, the MLMC and the F-MLMC
estimators, most of the variance is concentrated in the lower frequencies.
Although the low-frequency contribution to the variance is significantly
reduced by the MLMC estimator, a non-negligible contribution of the high
frequencies to the variance is noticeable. Adding filters reduces the error in
the high frequencies, as was observed previously in the spectral decomposition
of the variance (fig. 12), but it also reduces the error on the lower
frequencies, thus leading to lower total variance. In the second case (fig.
14(b)), corresponding to a tensor field $\mathbf{K}$ with smaller scales (see
fig. 9(b)), the MLMC estimators deteriorate the variance compared to the crude
MC estimator. The addition of a coarser grid ($L=1$) and the corresponding
grid transfer operators induce significant variance in the high-frequency
components, leading to an increased total variance. The degradation is all the
more pronounced as coarser grids are added ($L=2$ and $L=3$). We observe that
the addition of filters mitigates these effects and allows the F-MLMC
estimators to produce an effective reduction of the variance. The addition of
coarser grids ($L=2$ and $L=3$) does not bring further improvement, but does
not deteriorate the total variance either. The low-frequency (large-scale)
components are already well-captured by the two finer grids. This can be
explained by the fact that the spectral content of the samples is shifted to
higher frequencies. In such a case, a different grid hierarchy, starting from
a finer grid, would be more appropriate, although finer discretizations may
not always be available or affordable in an operational context. These results
are summarized in terms of total variance in fig. 15.

(a) $\mathbf{K}$ with $\mu_{\zeta}=0.12$ and $D_{\zeta}=0.2$ (fig. 9(a)).

(b) $\mathbf{K}$ with $\mu_{\zeta}=0.02$ and $D_{\zeta}=0.04$ (fig. 9(b)).

Figure 15: Total variance of the MC estimator ($L=0$) and of different MLMC
estimators and F-MLMC estimators with $L\in\\{1,2,3\\}$. The tensor field
$\mathbf{K}$ corresponds to $\mu_{\zeta}=0.12$ and $D_{\zeta}=0.2$ (left), and
$\mu_{\zeta}=0.02$ and $D_{\zeta}=0.04$ (right). The variance is approximated
from 500 estimators, each constructed with a computational budget
$\eta=100\mathcal{C}_{L}$.

##  7 Conclusion

In this paper, we focused on the estimation of the expectation of a
discretized field using a multilevel, MLMC-like estimator. The different
fidelity levels considered are grids of different resolutions, which requires
the use of grid transfer operators in the estimator. The resulting MLMC
estimator can then be used to reduce the variance of the estimation compared
to a crude MC estimator, as confirmed with an idealized 1D problem of
estimating the discretized intrinsic variance field of a diffusion-based
covariance operator. However, projecting the variance of the MLMC estimator
onto a spectral space revealed some discrepancy in the estimation of the
different scales of the discretized field. In our experiments, the MLMC
estimator was still able to still achieve a lower total variance by improving
the estimation of the low-frequency (large-scale) components compared to the
MC estimator, but at the expense of degrading the estimation in the higher-
frequency (smaller-scale) components.

Inspired by multigrid methods, we proposed an improvement of the MLMC
estimator by adding filtering operators, resulting in the F-MLMC estimator.
Filtering out the high-frequency components of a discretized field before
restriction and after prolongation removes spurious features, thus yielding a
better estimation of both the small- and large-scale components. These
improvements significantly impact the total variance of the F-MLMC estimators,
which is also reduced compared to the MLMC estimators in our experiments. In
the specific case of linear, symmetric, circulant simulators, we quantified
the effects of grid transfer and filtering operators on the total variance of
(F-)MLMC estimators, which allowed us to improve our understanding of the
influence of each ingredient. The proposed F-MLMC estimators were applied to
the problem of estimating the discretized intrinsic variance field of a 2D
diffusion-based covariance operator with a non-uniform diffusivity field,
which relies on non-linear simulators
$\tilde{f}_{\ell}\colon\mathbf{X}_{\ell}\mapsto(\mathbf{A}_{\ell}\mathbf{X}_{\ell})\odot(\mathbf{A}_{\ell}\mathbf{X}_{\ell})$.
The conclusions of these experiments were consistent with the theoretical
results derived in the spectral analysis for linear, symmetric, circulant
simulators. Specifically, F-MLMC estimators do reduce the variance in both the
low and high frequencies compared to their unfiltered counterparts, thus
improving the total variance. Even in experiments where MLMC estimators were
not able to reduce (and, actually, did deteriorate) the variance compared to a
crude MC estimator, F-MLMC estimators still achieved lower total variance.

It should be noted that, for the particular problem of estimating discretized
variance fields considered in this paper, the proposed MLMC and F-MLMC
estimators are not guaranteed to be almost surely non-negative. Although
negative estimates were not encountered in our experiments, this is
nonetheless a serious limitation of the multilevel estimators, as already
pointed out for the MLMC estimation of the variance of random variables
(possibly with values in Hilbert spaces) [2]. A similar issue exists for the
estimation of covariance matrices, for which multilevel estimators are not
guaranteed to be almost surely SPD [30], although advanced (but
computationally expensive) approaches have been proposed to design multilevel
estimators that are SPD by construction [38, 39]. These crucial issues still
constitute an open research area.

Nonetheless, the investigations conducted in this paper and the proposed
F-MLMC estimator expand the range of use of MLMC-like methods to discretized
fields. In our study, the use of second-order Shapiro filters demonstrated the
benefits of applying pre- and post-smoothing at each level of the MLMC
estimators. A potential next step would be to investigate whether conditions
can be derived for selecting the grid transfer and filtering operators,
similar to the conditions stated in [12] for multigrid methods. Furthermore,
extensions of the MLBLUE techniques [50, 49] to discretized fields may allow
for the derivation of optimal spectral weights for each scale component of the
discretized field to act as a post-prolongation filter [13, section 4.4],
possibly in combination with the spectral pre-restriction smoothing technique
of [31]. Finally, for the specific variance estimation problem described in
section 2, sampling $\mathbf{A}\mathbf{X}$ requires solving (many) systems of
linear equations of the form $\mathbf{M}\mathbf{u}=\mathbf{b}$, involving an
SPD matrix $\mathbf{M}:=\mathbf{W}(\mathbf{I}-\bm{\Delta})$. In our
experiments, we solved these systems numerically using a Cholesky
decomposition of $\mathbf{M}$. However, for large-scale problems, computing
such a decomposition is not computationally tractable and sparse iterative
methods are typically used instead [57]. An interesting research avenue would
consist of exploiting further the hierarchical nature of MLMC techniques by
combining them with multigrid iterative methods for solving the systems of
linear equations, using the same grid hierarchy, along the lines of [32, 33].

##  Appendix A Choice of grid transfer operators

We detail here our choice of prolongation and restriction operators.
Specifically, the prolongation operator
$\mathbf{P}:=\mathbf{P}_{\ell-1}^{\ell}$ between a coarse level $\ell-1$ and
the next fine level $\ell$ is designed such that it satisfies the norm-
preserving property given by the first identity in eq. 17. By the definition
of the weighted inner product eq. 6, this property can be recast as
$\mathbf{x}_{\ell-1}^{\mathrm{T}}\mathbf{P}^{\mathrm{T}}\mathbf{W}_{\ell}\mathbf{P}\mathbf{x}_{\ell-1}=\mathbf{x}_{\ell-1}^{\mathrm{T}}\mathbf{W}_{\ell-1}\mathbf{x}_{\ell-1}$
for all $\mathbf{x}_{\ell}\in\mathbb{R}^{n_{\ell}}$, which, in turn, amounts
to $\mathbf{P}^{\mathrm{T}}\mathbf{W}_{\ell}\mathbf{P}=\mathbf{W}_{\ell-1}$.
For the particular choice
$\mathbf{W}_{\ell}=n_{\ell}^{-1}\mathbf{I}_{n_{\ell}}$, for $\ell=0,\ldots,L$,
this property becomes
$\mathbf{P}^{\mathrm{T}}\mathbf{P}=(n_{\ell}/n_{\ell-1})\mathbf{I}_{n_{\ell-1}}$.
We note that the prolongation operator defined in eq. 16 satisfies this
property with a refinement factor $n_{\ell}/n_{\ell-1}=2$. Furthermore, for
the variance estimation problem described in section 2 where normal random
vectors need to be restricted to coarser grids, we want the second identity in
eq. 17 to hold so that the restricted signals remain normally distributed with
the same mean and variance on the coarse grids. With a linear restriction
operator $\mathbf{R}:=\mathbf{R}_{\ell}^{\ell-1}$, a zero-mean normal random
vector trivially remains a zero-mean normal random vector when restricted to a
coarser grid. Importantly, we require the restricted normal random vector to
have a unit variance. To achieve this, we define
$\mathbf{R}=\mathbf{V}_{\ell-1}^{-1}\mathbf{P}^{\mathrm{T}}\mathbf{V}_{\ell}$
as in eq. 16, which implies that

$\mathbf{R}\mathbf{R}^{\mathrm{T}}=\mathbf{V}_{\ell-1}^{-1}\mathbf{P}^{\mathrm{T}}\mathbf{V}_{\ell}\mathbf{V}_{\ell}^{\mathrm{T}}\mathbf{P}(\mathbf{V}_{\ell-1}^{-1})^{\mathrm{T}}=\mathbf{V}_{\ell-1}^{-1}\mathbf{P}^{\mathrm{T}}\mathbf{W}_{\ell}\mathbf{P}(\mathbf{V}_{\ell-1}^{\mathrm{T}})^{-1}=\mathbf{V}_{\ell-1}^{-1}\mathbf{W}_{\ell-1}(\mathbf{V}_{\ell-1}^{\mathrm{T}})^{-1},$
(55)

where the last identity holds since
$\mathbf{P}^{\mathrm{T}}\mathbf{W}_{\ell}\mathbf{P}=\mathbf{W}_{\ell-1}$. The
decomposition
$\mathbf{W}_{\ell-1}=\mathbf{V}_{\ell-1}\mathbf{V}_{\ell-1}^{\mathrm{T}}$ then
implies that $\mathbf{R}\mathbf{R}^{\mathrm{T}}=\mathbf{I}_{n_{\ell-1}}$, and
hence that the unit variance is conserved on the coarser grids.

##  Appendix B Orthogonality of the Hartley matrix

We focus here on a Hartley matrix $\mathbf{H}$ of arbitrary size $n$,

$(\mathbf{H})_{j,k}:=\frac{1}{\sqrt{n}}\left(\cos\alpha_{jk}+\sin\alpha_{jk}\right),\quad\alpha_{jk}:=\frac{(2j+1)k\pi}{n},\quad\forall
j,k=0,\dots,n-1.$ (56)

We have, for $i,j=0,\ldots,n-1$,

$\displaystyle(\mathbf{H}\mathbf{H}^{\mathrm{T}})_{i,j}$
$\displaystyle=\frac{1}{n}\sum_{k=0}^{n-1}(\cos\alpha_{ik}+\sin\alpha_{ik})(\cos\alpha_{jk}+\sin\alpha_{jk})$
(57)
$\displaystyle=\frac{1}{n}\sum_{k=0}^{n-1}\left[\cos(\alpha_{ik}-\alpha_{jk})+\sin(\alpha_{ik}+\alpha_{jk})\right]$
(58)
$\displaystyle=\frac{1}{n}\sum_{k=0}^{n-1}\cos\frac{2(i-j)k\pi}{n}+\frac{1}{n}\sum_{k=0}^{n-1}\sin\frac{2(i+j+1)k\pi}{n}.$
(59)

From [21, 1.342 (1&2)], we deduce, for $i,j=0,\ldots,n-1$,

$\displaystyle\sum_{k=0}^{n-1}\cos\frac{2(i-j)k\pi}{n}$
$\displaystyle=n\delta_{ij},$
$\displaystyle\sum_{k=0}^{n-1}\sin\frac{2(i+j+1)k\pi}{n}$ $\displaystyle=0,$
(60)

where $\delta_{ij}$ denotes the Kronecker delta, thus proving that
$\mathbf{H}\mathbf{H}^{\mathrm{T}}=\mathbf{I}_{n}$. Similarly,

$(\mathbf{H}^{\mathrm{T}}\mathbf{H})_{i,j}=\frac{1}{n}\sum_{k=0}^{n-1}\cos\frac{(2k+1)(i-j)\pi}{n}+\frac{1}{n}\sum_{k=0}^{n-1}\sin\frac{(2k+1)(i+j)\pi}{n},$
(61)

and $\mathbf{H}^{\mathrm{T}}\mathbf{H}=\mathbf{I}_{n}$ follows from [21, 1.342
(3&4)].

##  Appendix C Reordering of the Hartley basis vectors

For the plots of the spectral variance, the columns of the 1D Hartley matrix
$\mathbf{H}_{L}$ used in eq. 20 are actually reordered so that they are sorted
by increasing representable frequency on the corresponding discrete grid (see
fig. 7, which depicts the Hartley basis vectors without reordering).
Specifically, the new matrix with reordered columns is obtained as
$\mathbf{H}_{L}\bm{\Pi}$, where $\bm{\Pi}=(\Pi_{j,k})_{j,k=0}^{n_{L}-1}$ is
the permutation matrix defined by $\Pi_{j,2k}=\delta_{j,k}$ and
$\Pi_{j,2k+1}=\delta_{j,n_{L}-k-1}$, for $j=0,\ldots,n_{L}-1$ and
$k=0,\ldots,n_{L}/2-1$ (assuming $n_{L}$ is even). The 2D Hartley matrix with
reordered columns is constructed as the Kronecker product of two 1D Hartley
matrices, each with reordered columns.

##  Appendix D Proof of eq. 27

By linearity of the expectation operator, and denoting
$\dot{\mathbf{X}}_{L}:=\mathbf{X}_{L}-\mathbb{E}[\mathbf{X}_{L}]$, we have
that
$\mathcal{V}(Y_{0}):=\mathbb{E}[\|Y_{0}-\mathbb{E}[Y_{0}]\|^{2}_{\mathbf{W}_{L}}]=\mathbb{E}[\|\mathbf{F}_{0}\dot{\mathbf{X}}_{L}\|^{2}_{\mathbf{W}_{L}}]$.
Likewise, for $\ell=1,\ldots,L$, we have
$\mathcal{V}(Y_{\ell}-Y_{\ell-1})=\mathbb{E}[\|(\mathbf{F}_{\ell}-\mathbf{F}_{\ell-1})\dot{\mathbf{X}}_{L}\|^{2}_{\mathbf{W}_{L}}]$.
Now, for any $\mathbf{M}\in\mathbb{R}^{n_{L}\times n_{L}}$,

$\displaystyle\mathbb{E}[\|\mathbf{M}\dot{\mathbf{X}}_{L}\|^{2}_{\mathbf{W}_{L}}]$
$\displaystyle=\mathbb{E}[(\mathbf{M}\dot{\mathbf{X}}_{L})^{\mathrm{T}}\mathbf{W}_{L}(\mathbf{M}\dot{\mathbf{X}}_{L})]=\operatorname{tr}\mathbb{E}[\dot{\mathbf{X}}_{L}^{\mathrm{T}}\mathbf{M}^{\mathrm{T}}\mathbf{W}_{L}\mathbf{M}\dot{\mathbf{X}}_{L}]$
(62)
$\displaystyle=\mathbb{E}[\operatorname{tr}(\dot{\mathbf{X}}_{L}^{\mathrm{T}}\mathbf{M}^{\mathrm{T}}\mathbf{W}_{L}\mathbf{M}\dot{\mathbf{X}}_{L})]=\mathbb{E}[\operatorname{tr}(\mathbf{M}^{\mathrm{T}}\mathbf{W}_{L}\mathbf{M}\dot{\mathbf{X}}_{L}\dot{\mathbf{X}}_{L}^{\mathrm{T}})]$
(63)
$\displaystyle=\operatorname{tr}(\mathbf{M}^{\mathrm{T}}\mathbf{W}_{L}\mathbf{M}\mathbb{E}[\dot{\mathbf{X}}_{L}\dot{\mathbf{X}}_{L}^{\mathrm{T}}])=\operatorname{tr}(\mathbf{M}^{\mathrm{T}}\mathbf{W}_{L}\mathbf{M}\mathbf{G})$
(64)
$\displaystyle=\operatorname{tr}((\mathbf{M}\mathbf{G}^{1/2})^{\mathrm{T}}\mathbf{W}_{L}(\mathbf{M}\mathbf{G}^{1/2}))=\|\mathbf{M}\mathbf{G}^{1/2}\|^{2}_{F,\mathbf{W}_{L}},$
(65)

where $\operatorname{tr}(\cdot)$ denotes the matrix trace operator and
$\mathbf{G}$ denotes the covariance matrix of $\mathbf{X}_{L}$. Furthermore,
assuming that $\mathbf{W}_{L}=w\mathbf{I}_{n_{L}}$ for some $w>0$, for any
$\mathbf{M}\in\mathbb{R}^{n_{L}\times n_{L}}$ and any orthogonal
$\mathbf{Q}\in\mathbb{R}^{n_{L}\times n_{L}}$, we have
$\|\mathbf{Q}^{\mathrm{T}}\mathbf{M}\|^{2}_{F,\mathbf{W}_{L}}=\operatorname{tr}(w\mathbf{M}^{\mathrm{T}}\mathbf{Q}\mathbf{Q}^{\mathrm{T}}\mathbf{M})=\operatorname{tr}(\mathbf{M}^{\mathrm{T}}\mathbf{W}_{L}\mathbf{M})=\|\mathbf{M}\|^{2}_{F,\mathbf{W}_{L}}$,
so that eq. 28 follows.

##  Appendix E Proof of eq. 29

For $j,k=0,\ldots,n_{0}-1$,

$(\mathbf{P}\mathbf{h}^{0}_{k})_{2j}=(\mathbf{P}\mathbf{h}^{0}_{k})_{2j+1}=(\mathbf{h}^{0}_{k})_{j}=\frac{\sqrt{2}}{\sqrt{n_{1}}}\left[\cos\frac{(4j+2)k\pi}{n_{1}}+\sin\frac{(4j+2)k\pi}{n_{1}}\right].$
(66)

Upon writing $(4j+2)k\pi=(4j+1)k\pi+k\pi$ and applying elementary
trigonometric identitites,

$(\mathbf{P}\mathbf{h}^{0}_{k})_{2j}=\sqrt{2}c_{k}(\mathbf{h}^{1}_{k})_{2j}+\frac{\sqrt{2}}{\sqrt{n_{1}}}\left[\sin\frac{k\pi}{n_{1}}\left(\cos\frac{(4j+1)k\pi}{n_{1}}-\sin\frac{(4j+1)k\pi}{n_{1}}\right)\right].$
(67)

Noticing that

$\displaystyle\sin\frac{k\pi}{n_{1}}$
$\displaystyle=\sin\mathopen{}\left(\frac{(n_{0}+k)\pi}{2n_{0}}-\frac{\pi}{2}\right)\mathclose{}=-c_{n_{0}+k},$
(68) $\displaystyle\cos\frac{(4j+1)k\pi}{n_{1}}$
$\displaystyle=\cos\mathopen{}\left(\frac{(4j+1)(n_{0}+k)\pi}{n_{1}}-\frac{\pi}{2}\right)\mathclose{}=\sin\frac{(4j+1)k\pi}{n_{1}},$
(69) $\displaystyle\sin\frac{(4j+1)k\pi}{n_{1}}$
$\displaystyle=\sin\mathopen{}\left(\frac{(4j+1)(n_{0}+k)\pi}{n_{1}}-\frac{\pi}{2}\right)\mathclose{}=-\cos\frac{(4j+1)k\pi}{n_{1}},$
(70)

we conclude that
$(\mathbf{P}\mathbf{h}^{0}_{k})_{2j}=\sqrt{2}(c_{k}(\mathbf{h}^{1}_{k})_{2j}-c_{n_{0}+k}(\mathbf{h}^{1}_{n_{0}+k})_{2j})$.
A similar derivation allows us to show that
$(\mathbf{P}\mathbf{h}^{0}_{k})_{2j+1}=\sqrt{2}(c_{k}(\mathbf{h}^{1}_{k})_{2j+1}-c_{n_{0}+k}(\mathbf{h}^{1}_{n_{0}+k})_{2j+1})$,
thus proving eq. 29.

##  Appendix F Proof of eq. 31

For $j=0,\ldots,n_{0}-1$ and $k=0,\ldots,2n_{0}-1$,

$\displaystyle(\mathbf{R}\mathbf{h}^{1}_{k})_{j}$
$\displaystyle=\frac{1}{\sqrt{2}}((\mathbf{h}^{1}_{k})_{2j}+(\mathbf{h}^{1}_{k})_{2j+1})$
(71)
$\displaystyle=\frac{1}{\sqrt{2n_{1}}}\left[\cos\frac{(4j+1)k\pi}{n_{1}}+\cos\frac{(4j+3)k\pi}{n_{1}}+\sin\frac{(4j+1)k\pi}{n_{1}}+\sin\frac{(4j+3)k\pi}{n_{1}}\right]$
(72)
$\displaystyle=\frac{1}{\sqrt{2n_{1}}}\left[2\cos\frac{(8j+4)k\pi}{2n_{1}}\cos\frac{2k\pi}{2n_{1}}+2\sin\frac{(8j+4)k\pi}{2n_{1}}\cos\frac{2k\pi}{2n_{1}}\right]$
(73)
$\displaystyle=\frac{c_{k}}{\sqrt{n_{0}}}\left[\cos\frac{(2j+1)k\pi}{n_{0}}+\sin\frac{(2j+1)k\pi}{n_{0}}\right].$
(74)

Hence, for $j,k=0,\ldots,n_{0}-1$, it follows immediately that
$\mathbf{R}\mathbf{h}^{1}_{k}=c_{k}\mathbf{h}^{0}_{k}$. Furthermore,

$(\mathbf{R}\mathbf{h}^{1}_{n_{0}+k})_{j}=\frac{c_{n_{0}+k}}{\sqrt{n_{0}}}\left[\cos\mathopen{}\left(\frac{(2j+1)k\pi}{n_{0}}+\pi\right)\mathclose{}+\sin\mathopen{}\left(\frac{(2j+1)k\pi}{n_{0}}+\pi\right)\mathclose{}\right]=-c_{n_{0}+k}\mathbf{h}^{0}_{k},$
(75)

thus proving eq. 31.

##  Appendix G Symmetric circulant matrices are diagonalizable in the Hartley
basis

In this appendix, we prove theorem 2 below, which states that symmetric,
circulant matrices can be diagonalized in the cell-centered Hartley basis
$\mathbf{H}$ defined by eq. 56. To do so, we start by recalling or proving
results on the node-centered Fourier basis
$\check{\mathbf{F}}\in\mathbb{C}^{n\times n}$ and Hartley bases
$\check{\mathbf{H}}^{\pm}\in\mathbb{R}^{n\times n}$ defined by

$\displaystyle\check{\mathbf{F}}$
$\displaystyle:=\check{\mathbf{H}}_{c}+i\check{\mathbf{H}}_{s},$
$\displaystyle\check{\mathbf{H}}^{\pm}$
$\displaystyle:=\check{\mathbf{H}}_{c}\pm\check{\mathbf{H}}_{s},$
$\displaystyle(\check{\mathbf{H}}_{c})_{j,k}$
$\displaystyle:=\frac{1}{\sqrt{n}}\cos\frac{2jk\pi}{n},$
$\displaystyle(\check{\mathbf{H}}_{s})_{j,k}$
$\displaystyle:=\frac{1}{\sqrt{n}}\sin\frac{2jk\pi}{n},$ (76)

where $i\in\mathbb{C}$ denotes the unit imaginary number such that $i^{2}=-1$.
It is clear that $\check{\mathbf{H}}_{c}$ and $\check{\mathbf{H}}_{s}$ are
real, symmetric matrices, and hence so are $\check{\mathbf{H}}^{\pm}$, while
$\check{\mathbf{F}}$ is a complex, symmetric (but not Hermitian) matrix.
Furthermore, $\check{\mathbf{F}}$ is unitary (see, e.g., [11]), i.e.,
$\check{\mathbf{F}}^{*}\check{\mathbf{F}}=\check{\mathbf{F}}\check{\mathbf{F}}^{*}=\mathbf{I}_{n}$,
where $\check{\mathbf{F}}^{*}=\check{\mathbf{H}}_{c}-i\check{\mathbf{H}}_{s}$
is the Hermitian transpose of $\check{\mathbf{F}}$.

######  Lemma 1.

$\check{\mathbf{H}}_{c}$ and $\check{\mathbf{H}}_{s}$ are such that

$\displaystyle\check{\mathbf{H}}_{c}^{2}+\check{\mathbf{H}}_{s}^{2}$
$\displaystyle=\mathbf{I}_{n},$
$\displaystyle\check{\mathbf{H}}_{c}\check{\mathbf{H}}_{s}$
$\displaystyle=\check{\mathbf{H}}_{s}\check{\mathbf{H}}_{c}=\bm{0}_{n},$
$\displaystyle\check{\mathbf{H}}_{c}\bm{1}_{n}$
$\displaystyle=\sqrt{n}\mathbf{e}_{1},$
$\displaystyle\check{\mathbf{H}}_{s}\bm{1}_{n}$ $\displaystyle=\bm{0}_{n},$
(77)

where $\bm{0}_{n}:=(0,\ldots,0)^{\mathrm{T}}\in\mathbb{R}^{n}$,
$\bm{1}_{n}:=(1,\ldots,1)^{\mathrm{T}}\in\mathbb{R}^{n}$, and $\mathbf{e}_{1}$
denotes the first column of $\mathbf{I}_{n}$.

###### Proof.

The first two identities follow from [4, Lemma 1], while the last two
identities follow from [21, 1.342 (1&2)]. ∎

######  Corollary 1.

$\check{\mathbf{H}}^{\pm}\bm{1}_{n}=\sqrt{n}\mathbf{e}_{1}$.

######  Corollary 2.

$\check{\mathbf{H}}^{\pm}$ are orthogonal, i.e.,
$(\check{\mathbf{H}}^{\pm})^{2}=\mathbf{I}_{n}$. Furthermore,
$\check{\mathbf{H}}^{+}\check{\mathbf{H}}^{-}=\check{\mathbf{H}}^{-}\check{\mathbf{H}}^{+}=\check{\mathbf{F}}^{2}$.

###### Proof.

By the definitions eq. 76 and lemma 1, we have

$\displaystyle\check{\mathbf{H}}^{\pm}\bm{1}_{n}$
$\displaystyle=\check{\mathbf{H}}_{c}\bm{1}_{n}\pm\check{\mathbf{H}}_{s}\bm{1}_{n}=\check{\mathbf{H}}_{c}\bm{1}_{n}=\sqrt{n}\mathbf{e}_{1},$
(78) $\displaystyle(\check{\mathbf{H}}^{\pm})^{2}$
$\displaystyle=(\check{\mathbf{H}}_{c}^{2}+\check{\mathbf{H}}_{s}^{2})\pm(\check{\mathbf{H}}_{c}\check{\mathbf{H}}_{s}+\check{\mathbf{H}}_{s}\check{\mathbf{H}}_{c})=\check{\mathbf{H}}_{c}^{2}+\check{\mathbf{H}}_{s}^{2}=\mathbf{I}_{n},$
(79) $\displaystyle\check{\mathbf{F}}^{2}$
$\displaystyle=\check{\mathbf{H}}_{c}^{2}-\check{\mathbf{H}}_{s}^{2}+i(\check{\mathbf{H}}_{c}\check{\mathbf{H}}_{s}+\check{\mathbf{H}}_{s}\check{\mathbf{H}}_{c})=\check{\mathbf{H}}_{c}^{2}-\check{\mathbf{H}}_{s}^{2},$
(80) $\displaystyle\check{\mathbf{H}}^{-}\check{\mathbf{H}}^{+}$
$\displaystyle=\check{\mathbf{H}}_{c}^{2}-\check{\mathbf{H}}_{s}^{2}+(\check{\mathbf{H}}_{c}\check{\mathbf{H}}_{s}-\check{\mathbf{H}}_{s}\check{\mathbf{H}}_{c})=\check{\mathbf{F}}^{2}=\check{\mathbf{H}}^{+}\check{\mathbf{H}}^{-}.$
(81)

∎

######  Lemma 2.

Let $\mathbf{A}=\operatorname{Circ}(\mathbf{a})\in\mathbb{R}^{n\times n}$ be a
symmetric, circulant matrix whose first column is
$\mathbf{a}=(a_{k})_{k=0}^{n-1}\in\mathbb{R}^{n}$, with $a_{n-i}=a_{i}$, for
$i=1,\ldots,n-1$. Then
$\check{\mathbf{H}}_{c}\mathbf{A}\check{\mathbf{H}}_{s}+\check{\mathbf{H}}_{s}\mathbf{A}\check{\mathbf{H}}_{c}=\check{\mathbf{H}}_{c}\mathbf{A}\check{\mathbf{H}}_{s}-\check{\mathbf{H}}_{s}\mathbf{A}\check{\mathbf{H}}_{c}=\bm{0}_{n}$.

###### Proof.

Proven in the intermediary steps of the proof of [4, Theorem 1]. ∎

######  Theorem 1.

Let $\mathbf{A}=\operatorname{Circ}(\mathbf{a})\in\mathbb{R}^{n\times n}$ be a
symmetric, circulant matrix whose first column is
$\mathbf{a}=(a_{k})_{k=0}^{n-1}\in\mathbb{R}^{n}$, with $a_{n-i}=a_{i}$, for
$i=1,\ldots,n-1$. Then
$\check{\mathbf{H}}^{+}\mathbf{A}\check{\mathbf{H}}^{+}=\check{\mathbf{H}}^{-}\mathbf{A}\check{\mathbf{H}}^{-}=\bm{\Lambda}$,
where
$\bm{\Lambda}=\sqrt{n}\operatorname{Diag}(\check{\mathbf{H}}^{+}\mathbf{a})=\sqrt{n}\operatorname{Diag}(\check{\mathbf{H}}^{-}\mathbf{a})=\sqrt{n}\operatorname{Diag}(\check{\mathbf{H}}_{c}\mathbf{a})$.

###### Proof.

The fact that
$\check{\mathbf{H}}^{+}\mathbf{A}\check{\mathbf{H}}^{+}=\sqrt{n}\operatorname{Diag}(\check{\mathbf{H}}^{+}\mathbf{a})=\bm{\Lambda}=\check{\mathbf{H}}_{c}\mathbf{A}\check{\mathbf{H}}_{c}+\check{\mathbf{H}}_{s}\mathbf{A}\check{\mathbf{H}}_{s}$
follows from [4, Theorem 1], since symmetric, circulant matrices belong to the
larger class of matrices considered in [4]. Then, from the definition of
$\check{\mathbf{H}}^{-}$ and lemma 2, we have

$\check{\mathbf{H}}^{-}\mathbf{A}\check{\mathbf{H}}^{-}=(\check{\mathbf{H}}_{c}\mathbf{A}\check{\mathbf{H}}_{c}+\check{\mathbf{H}}_{s}\mathbf{A}\check{\mathbf{H}}_{s})-(\check{\mathbf{H}}_{c}\mathbf{A}\check{\mathbf{H}}_{s}+\check{\mathbf{H}}_{s}\mathbf{A}\check{\mathbf{H}}_{c})=\check{\mathbf{H}}_{c}\mathbf{A}\check{\mathbf{H}}_{c}+\check{\mathbf{H}}_{s}\mathbf{A}\check{\mathbf{H}}_{s}=\bm{\Lambda}.$
(82)

Finally, from lemma 1 and corollary 1,

$\operatorname{diag}(\bm{\Lambda})=\bm{\Lambda}\bm{1}_{n}=\check{\mathbf{H}}^{-}\mathbf{A}\check{\mathbf{H}}^{-}\bm{1}_{n}=\check{\mathbf{H}}_{c}\mathbf{A}\check{\mathbf{H}}_{c}\bm{1}_{n}+\check{\mathbf{H}}_{s}\mathbf{A}\check{\mathbf{H}}_{s}\bm{1}_{n}=\sqrt{n}\check{\mathbf{H}}^{-}\mathbf{A}\mathbf{e}_{1}=\sqrt{n}\check{\mathbf{H}}_{c}\mathbf{A}\mathbf{e}_{1},$
(83)

which concludes the proof, since $\mathbf{A}\mathbf{e}_{1}=\mathbf{a}$. ∎

######  Corollary 3.

Let $\mathbf{A}=\operatorname{Circ}(\mathbf{a})\in\mathbb{R}^{n\times n}$ be a
symmetric, circulant matrix and let $\bm{\Lambda}$ be defined as in theorem 1.
Then
$\check{\mathbf{H}}^{+}\mathbf{A}\check{\mathbf{H}}^{-}=\check{\mathbf{H}}^{-}\mathbf{A}\check{\mathbf{H}}^{+}=\bm{\Lambda}\check{\mathbf{F}}^{2}=\check{\mathbf{F}}^{2}\bm{\Lambda}$.

###### Proof.

From theorem 1 and corollary 2, we have

$\displaystyle\check{\mathbf{H}}^{+}\mathbf{A}\check{\mathbf{H}}^{-}$
$\displaystyle=\check{\mathbf{H}}^{+}\mathbf{A}\check{\mathbf{H}}^{+}\check{\mathbf{H}}^{+}\check{\mathbf{H}}^{-}=\check{\mathbf{H}}^{+}\mathbf{A}\check{\mathbf{H}}^{+}\check{\mathbf{F}}^{2}=\bm{\Lambda}\check{\mathbf{F}}^{2},$
(84) $\displaystyle\check{\mathbf{H}}^{-}\mathbf{A}\check{\mathbf{H}}^{+}$
$\displaystyle=\check{\mathbf{H}}^{-}\mathbf{A}\check{\mathbf{H}}^{-}\check{\mathbf{H}}^{-}\check{\mathbf{H}}^{+}=\check{\mathbf{H}}^{-}\mathbf{A}\check{\mathbf{H}}^{-}\check{\mathbf{F}}^{2}=\bm{\Lambda}\check{\mathbf{F}}^{2},$
(85)
$\displaystyle(\check{\mathbf{H}}^{+}\mathbf{A}\check{\mathbf{H}}^{-})^{\mathrm{T}}$
$\displaystyle=(\bm{\Lambda}\check{\mathbf{F}}^{2})^{\mathrm{T}}=\check{\mathbf{F}}^{2}\bm{\Lambda}=\check{\mathbf{H}}^{-}\mathbf{A}\check{\mathbf{H}}^{+}=\bm{\Lambda}\check{\mathbf{F}}^{2},$
(86)

where the last row follows from the symmetry of $\mathbf{A}$,
$\check{\mathbf{H}}^{\pm}$ and $\check{\mathbf{F}}^{2}$. ∎

######  Lemma 3.

$\mathbf{H}=\check{\mathbf{H}}^{+}\mathbf{C}+\check{\mathbf{H}}^{-}\mathbf{S}$,
with $\mathbf{C}:=\operatorname{Diag}(\\{\cos\frac{k\pi}{n}\\}_{k=0}^{n-1})$
and $\mathbf{S}:=\operatorname{Diag}(\\{\sin\frac{k\pi}{n}\\}_{k=0}^{n-1})$.

###### Proof.

From elementary trigonometric identities, we have

$\displaystyle\cos\frac{(2j+1)k\pi}{n}$
$\displaystyle=\cos\frac{2jk\pi}{n}\cos\frac{k\pi}{n}-\sin\frac{2jk\pi}{n}\sin\frac{k\pi}{n},$
(87) $\displaystyle\sin\frac{(2j+1)k\pi}{n}$
$\displaystyle=\sin\frac{2jk\pi}{n}\cos\frac{k\pi}{n}-\cos\frac{2jk\pi}{n}\sin\frac{k\pi}{n},$
(88)

then lemma 3 follows. ∎

######  Lemma 4.

$\mathbf{C}\check{\mathbf{F}}^{2}\mathbf{S}+\mathbf{S}\check{\mathbf{F}}^{2}\mathbf{C}=\bm{0}_{n}$,
where $\mathbf{C}$ and $\mathbf{S}$ are defined as in lemma 3.

###### Proof.

From [4, Lemma 1], we have that, for $j,k=0,\dots,n-1$,
$(\check{\mathbf{F}}^{2})_{j,k}=1$ if $j=k=0$ or $j+k=n$, and
$(\check{\mathbf{F}}^{2})_{j,k}=0$ otherwise. Moreover,
$(\mathbf{C}\check{\mathbf{F}}^{2}\mathbf{S})_{j,k}=(\check{\mathbf{F}}^{2})_{j,k}\cos\frac{j\pi}{n}\sin\frac{k\pi}{n}$,
and, by symmetry,
$\mathbf{S}\check{\mathbf{F}}^{2}\mathbf{C}=(\mathbf{C}\check{\mathbf{F}}^{2}\mathbf{S})^{\mathrm{T}}$.
Now, if $j=k=0$,
$(\mathbf{C}\check{\mathbf{F}}^{2}\mathbf{S})_{j,k}=(\mathbf{S}\check{\mathbf{F}}^{2}\mathbf{C})_{j,k}=0$
trivially. Otherwise, if $j\neq 0$, $k\neq 0$, and $j+k\neq n$,
$(\mathbf{C}\check{\mathbf{F}}^{2}\mathbf{S})_{j,k}=(\mathbf{S}\check{\mathbf{F}}^{2}\mathbf{C})_{j,k}=0$.
Finally, if $j+k=n$, elementary trigonometric identities induce

$\displaystyle(\mathbf{C}\check{\mathbf{F}}^{2}\mathbf{S})_{j,k}$
$\displaystyle=\cos\frac{j\pi}{n}\sin\frac{(n-j)\pi}{n}=\cos\frac{j\pi}{n}\sin\mathopen{}\left(\pi-\frac{j\pi}{n}\right)\mathclose{}=\cos\frac{j\pi}{n}\sin\frac{j\pi}{n},$
(89) $\displaystyle(\mathbf{C}\check{\mathbf{F}}^{2}\mathbf{S})_{k,j}$
$\displaystyle=\sin\frac{j\pi}{n}\cos\frac{(n-j)\pi}{n}=\sin\frac{j\pi}{n}\cos\mathopen{}\left(\pi-\frac{j\pi}{n}\right)\mathclose{}=-\cos\frac{j\pi}{n}\sin\frac{j\pi}{n},$
(90)

so that
$(\mathbf{C}\check{\mathbf{F}}^{2}\mathbf{S}+\mathbf{S}\check{\mathbf{F}}^{2}\mathbf{C})_{j,k}=(\mathbf{C}\check{\mathbf{F}}^{2}\mathbf{S})_{j,k}+(\mathbf{C}\check{\mathbf{F}}^{2}\mathbf{S})_{k,j}=0$,
thus completing the proof. ∎

We are now ready to state and prove theorem 2.

######  Theorem 2.

Let $\mathbf{A}=\operatorname{Circ}(\mathbf{a})\in\mathbb{R}^{n\times n}$ be a
symmetric, circulant matrix and let $\bm{\Lambda}$ be defined as in theorem 1.
Then ${\mathbf{H}^{\mathrm{T}}\mathbf{A}\mathbf{H}=\bm{\Lambda}}$.

###### Proof.

Starting from the expression of lemma 3, and applying theorem 1, corollary 3
and lemma 4, we obtain

$\displaystyle\mathbf{H}^{\mathrm{T}}\mathbf{A}\mathbf{H}$
$\displaystyle=\mathbf{C}\check{\mathbf{H}}^{+}\mathbf{A}\check{\mathbf{H}}^{+}\mathbf{C}+\mathbf{S}\check{\mathbf{H}}^{-}\mathbf{A}\check{\mathbf{H}}^{-}\mathbf{S}+\mathbf{C}\check{\mathbf{H}}^{+}\mathbf{A}\check{\mathbf{H}}^{-}\mathbf{S}+\mathbf{S}\check{\mathbf{H}}^{-}\mathbf{A}\check{\mathbf{H}}^{+}\mathbf{C}$
(91)
$\displaystyle=\mathbf{C}\bm{\Lambda}\mathbf{C}+\mathbf{S}\bm{\Lambda}\mathbf{S}+\mathbf{C}\bm{\Lambda}\check{\mathbf{F}}^{2}\mathbf{S}+\mathbf{S}\bm{\Lambda}\check{\mathbf{F}}^{2}\mathbf{C}=\bm{\Lambda}(\mathbf{C}^{2}+\mathbf{S}^{2}+\mathbf{C}\check{\mathbf{F}}^{2}\mathbf{S}+\mathbf{S}\check{\mathbf{F}}^{2}\mathbf{C}),$
(92)

so that
$\mathbf{H}^{\mathrm{T}}\mathbf{A}\mathbf{H}=\bm{\Lambda}(\mathbf{C}^{2}+\mathbf{S}^{2})=\bm{\Lambda}$
since $\mathbf{C}^{2}+\mathbf{S}^{2}=\mathbf{I}_{n}$ by elementary
trigonometry. ∎

##  Appendix H Proof of the optimality of the Galerkin operator

Let $n_{1},n_{0}\in\mathbb{N}$ be such that $n_{1}>n_{0}>0$, and let
$\mathbf{A}_{1}\in\mathbb{R}^{n_{1}\times n_{1}}$ and
$\mathbf{A}_{0}\in\mathbb{R}^{n_{0}\times n_{0}}$ be two linear operators. Let
$\mathbf{R}\in\mathbb{R}^{n_{0}\times n_{1}}$ and
$\mathbf{P}\in\mathbb{R}^{n_{1}\times n_{0}}$ be a restriction operator and a
prolongation operator, respectively, such that
$\mathbf{R}\mathbf{R}^{\mathrm{T}}=\mathbf{I}_{n_{0}}$ and
$\mathbf{P}=\alpha\mathbf{R}^{\mathrm{T}}$, for some $\alpha\in\mathbb{R}$. We
seek the operator $\mathbf{A}_{0}^{*}$ that minimizes
$\|\mathbf{A}_{1}-\mathbf{P}\mathbf{A}_{0}\mathbf{R}\|^{2}_{\textrm{F},\mathbf{W}_{1}}$
among the linear operators $\mathbf{A}_{0}\in\mathbb{R}^{n_{0}\times n_{0}}$.
We restrict ourselves to the case where $\mathbf{W}_{1}=w\mathbf{I}_{n_{1}}$,
with $w>0$, so that

$\mathbf{A}_{0}^{*}:=\operatorname*{arg\,min}_{\mathbf{A}_{0}\in\mathbb{R}^{n_{0}\times
n_{0}}}\|\mathbf{A}_{1}-\mathbf{P}\mathbf{A}_{0}\mathbf{R}\|^{2}_{\textrm{F},\mathbf{W}_{1}}=\operatorname*{arg\,min}_{\mathbf{A}_{0}\in\mathbb{R}^{n_{0}\times
n_{0}}}\|\mathbf{A}_{1}-\mathbf{P}\mathbf{A}_{0}\mathbf{R}\|^{2}_{\textrm{F}},$
(93)

where $\|\cdot\|_{\textrm{F}}$ denotes the classical (unweighted) Frobenius
norm. We start by writing

$\displaystyle\|\mathbf{A}_{1}-\mathbf{P}\mathbf{A}_{0}\mathbf{R}\|^{2}_{\textrm{F}}$
$\displaystyle=\|\mathbf{A}_{1}-\mathbf{P}\mathbf{A}_{0}\mathbf{R}\|^{2}_{\textrm{F}}=\|\operatorname{vec}(\mathbf{A}_{1}-\mathbf{P}\mathbf{A}_{0}\mathbf{R})\|^{2}_{2}$
(94)
$\displaystyle=\|\operatorname{vec}(\mathbf{A}_{1})-\operatorname{vec}(\mathbf{P}\mathbf{A}_{0}\mathbf{R})\|^{2}_{2}=\|\operatorname{vec}(\mathbf{A}_{1})-(\mathbf{R}^{\mathrm{T}}\otimes\mathbf{P})\operatorname{vec}(\mathbf{A}_{1})\|^{2}_{2},$
(95)

so that eq. 93 is recast as an ordinary linear least squares problem whose
solution is

$\displaystyle\operatorname{vec}(\mathbf{A}_{0}^{*})$
$\displaystyle=\left((\mathbf{R}^{\mathrm{T}}\otimes\mathbf{P})^{\mathrm{T}}(\mathbf{R}^{\mathrm{T}}\otimes\mathbf{P})\right)^{-1}(\mathbf{R}^{\mathrm{T}}\otimes\mathbf{P})^{\mathrm{T}}\operatorname{vec}(\mathbf{A}_{1})$
(96)
$\displaystyle=\left((\mathbf{R}\otimes\mathbf{P}^{\mathrm{T}})(\mathbf{R}^{\mathrm{T}}\otimes\mathbf{P})\right)^{-1}(\mathbf{R}\otimes\mathbf{P}^{\mathrm{T}})\operatorname{vec}(\mathbf{A}_{1})$
(97)
$\displaystyle=(\mathbf{R}\mathbf{R}^{\mathrm{T}}\otimes\mathbf{P}^{\mathrm{T}}\mathbf{P})^{-1}\operatorname{vec}(\mathbf{P}^{\mathrm{T}}\mathbf{A}_{1}\mathbf{R}^{\mathrm{T}})=(\mathbf{I}_{n_{0}}\otimes\alpha^{2}\mathbf{I}_{n_{0}})^{-1}\operatorname{vec}((\alpha\mathbf{R})\mathbf{A}_{1}(\alpha^{-1}\mathbf{P}))$
(98)
$\displaystyle=\alpha^{-2}(\mathbf{I}_{n_{0}}\otimes\mathbf{I}_{n_{0}})\operatorname{vec}(\mathbf{R}\mathbf{A}_{1}\mathbf{P})=\alpha^{-2}\operatorname{vec}(\mathbf{R}\mathbf{A}_{1}\mathbf{P}),$
(99)

thus proving that
$\mathbf{A}_{0}^{*}=\alpha^{-2}\mathbf{R}\mathbf{A}_{1}\mathbf{P}$.

##  Appendix I Proof of eqs. 43 and 44

We start by looking at the effect of the Shapiro filter on the fine Hartley
basis vectors,
$(\mathbf{S}_{1}\mathbf{h}^{1}_{k})_{j}=\frac{1}{4}\left[(\mathbf{h}^{1}_{k})_{j-1}+2(\mathbf{h}^{1}_{k})_{j}+(\mathbf{h}^{1}_{k})_{j+1}\right]$.
Elementary trigonometric identities imply that, for $j,k=0,\ldots,n_{1}-1$,

$\displaystyle\cos\frac{(2j-1)k\pi}{n_{1}}+\cos\frac{(2j+1)k\pi}{n_{1}}$
$\displaystyle=2c_{k}\cos\frac{2jk\pi}{n_{1}},$ (100)
$\displaystyle\cos\frac{(2j+1)k\pi}{n_{1}}+\cos\frac{(2j+3)k\pi}{n_{1}}$
$\displaystyle=2c_{k}\cos\frac{(2j+2)k\pi}{n_{1}},$ (101)
$\displaystyle\sin\frac{(2j-1)k\pi}{n_{1}}+\sin\frac{(2j+1)k\pi}{n_{1}}$
$\displaystyle=2c_{k}\sin\frac{2jk\pi}{n_{1}},$ (102)
$\displaystyle\sin\frac{(2j+1)k\pi}{n_{1}}+\sin\frac{(2j+3)k\pi}{n_{1}}$
$\displaystyle=2c_{k}\sin\frac{(2j+2)k\pi}{n_{1}},$ (103)

leading to

$(\mathbf{S}_{1}\mathbf{h}^{1}_{k})_{j}=\frac{c_{k}}{2}\left[\cos\frac{2jk\pi}{n_{1}}+\cos\frac{(2j+2)k\pi}{n_{1}}+\sin\frac{2jk\pi}{n_{1}}+\sin\frac{(2j+2)k\pi}{n_{1}}\right]=c_{k}^{2}(\mathbf{h}^{1}_{k})_{j},$
(104)

where the last identity is obtained by applying the same trigonometric
identities as for eqs. 100 to 103, and thus showing that
$\mathbf{S}_{1}\mathbf{h}^{1}_{k}=c_{k}^{2}\mathbf{h}^{1}_{k}$, for
$k=0,\ldots,n_{1}-1$. Then, eq. 29 implies that, for $k=0,\ldots,n_{0}-1$,
$\bar{\mathbf{P}}\mathbf{h}^{0}_{k}=\mathbf{S}_{1}\mathbf{P}\mathbf{h}^{0}_{k}=\sqrt{2}\left[c_{k}\mathbf{S}_{1}\mathbf{h}^{1}_{k}-c_{n_{0}+k}\mathbf{S}_{1}\mathbf{h}^{1}_{n_{0}+k}\right]$,
from which eq. 43 follows. Furthermore, eq. 44 follows immediately from the
direct application of eq. 31 to
$\bar{\mathbf{R}}\mathbf{h}^{0}_{k}=\mathbf{R}\mathbf{S}_{1}\mathbf{h}^{0}_{k}=c_{k}^{2}\mathbf{R}\mathbf{h}^{0}_{k}$.

## Acknowledgments

This project has received financial support from the CNRS (Centre National de
la Recherche Scientifique) through the 80|Prime program and the French
national program LEFE (Les Enveloppes Fluides et l’Environnement).

For the purpose of Open Access, a CC-BY public copyright licence has been
applied by the authors to the present document and will be applied to all
subsequent versions up to the Author Accepted Manuscript arising from this
submission.

## References

  * [1] C. Bekas, E. Kokiopoulou and Y. Saad  “An Estimator for the Diagonal of a Matrix”  In _Applied Numerical Mathematics_ 57.11-12, 2007, pp. 1214–1229  DOI: 10.1016/j.apnum.2007.01.003
  * [2] Claudio Bierig and Alexey Chernov  “Convergence analysis of multilevel Monte Carlo variance estimators and application for random obstacle problems”  In _Numerische Mathematik_ 130.4, 2015, pp. 579–613  DOI: 10.1007/s00211-014-0676-3
  * [3] Claudio Bierig and Alexey Chernov  “Estimation of arbitrary order central statistical moments by the multilevel Monte Carlo method”  In _Stochastics and Partial Differential Equations Analysis and Computations_ 4.1, 2016, pp. 3–40  DOI: 10.1007/s40072-015-0063-9
  * [4] Dario Bini and Paola Favati  “On a Matrix Algebra Related to the Discrete Hartley Transform”  In _SIAM Journal on Matrix Analysis and Applications_ 14.2 Society for IndustrialApplied Mathematics, 1993, pp. 500–507  DOI: 10.1137/0614035
  * [5] R.. Bracewell  “Discrete Hartley Transform”  In _JOSA_ 73.12 Optica Publishing Group, 1983, pp. 1832–1835  DOI: 10.1364/JOSA.73.001832
  * [6] Achi Brandt  “Guide to Multigrid Development”  In _Multigrid Methods_ , Lecture Notes in Mathematics  Berlin, Heidelberg: Springer, 1982, pp. 220–312  DOI: 10.1007/BFb0069930
  * [7] William L. Briggs, Van Emden Henson and Steve F. McCormick  “A Multigrid Tutorial, Second Edition”  Society for IndustrialApplied Mathematics, 2000  DOI: 10.1137/1.9780898719505
  * [8] Tan Bui-Thanh, Omar Ghattas, James Martin and Georg Stadler  “A Computational Framework for Infinite-Dimensional Bayesian Inverse Problems Part I: The Linearized Case, with Application to Global Seismic Inversion”  In _SIAM Journal on Scientific Computing_ 35.6, 2013, pp. A2494–A2523  DOI: 10.1137/12089586X
  * [9] K.. Cliffe, M.. Giles, R. Scheichl and A.. Teckentrup  “Multilevel Monte Carlo Methods and Applications to Elliptic PDEs with Random Coefficients”  In _Computing and Visualization in Science_ 14.1, 2011, pp. 3–15  DOI: 10.1007/s00791-011-0160-x
  * [10] M. Croci, K.. Willcox and S.. Wright  “Multi-Output Multilevel Best Linear Unbiased Estimators via Semidefinite Programming”  In _Computer Methods in Applied Mechanics and Engineering_ 413, 2023, pp. 116130  DOI: 10.1016/j.cma.2023.116130
  * [11] Philip J. Davis  “Circulant Matrices: Second Edition”  American Mathematical Society, 2012  URL: https://bookstore.ams.org/chel-338
  * [12] Laurent Debreu, Emilie Neveu, Ehouarn Simon, François-Xavier Le Dimet and Arthur Vidard  “Multigrid solvers and multigrid preconditioners for the solution of variational data assimilation problems”  In _Quarterly Journal of the Royal Meteorological Society_ 142.694, 2016, pp. 515–528  DOI: 10.1002/qj.2676
  * [13] Mayeul Destouches, Paul Mycek and Selime Gürol  “Multivariate Extensions of the Multilevel Best Linear Unbiased Estimator for Ensemble-Variational Data Assimilation”, 2023  DOI: 10.48550/arXiv.2306.07017
  * [14] Thomas O. Dixon, James E. Warner, Geoffrey F. Bomarito and Alex A. Gorodetsky  “Covariance Expressions for Multi-Fidelity Sampling with Multi-Output, Multi-Statistic Estimators: Application to Approximate Control Variates”, 2023  DOI: 10.48550/arXiv.2310.00125
  * [15] D. Drzisga, B. Gmeiner, U. Rüde, R. Scheichl and B. Wohlmuth  “Scheduling Massively Parallel Multigrid for Multilevel Monte Carlo Methods”  In _SIAM Journal on Scientific Computing_ 39.5 Society for Industrial and Applied Mathematics, 2017, pp. S873–S897  DOI: 10.1137/16M1083591
  * [16] F. Falissard  “Genuinely multi-dimensional explicit and implicit generalized Shapiro filters for weather forecasting, computational fluid dynamics and aeroacoustics”  In _Journal of Computational Physics_ 253, 2013, pp. 344–367  DOI: 10.1016/j.jcp.2013.07.001
  * [17] F. Falissard  “Uneven-order decentered Shapiro filters for boundary filtering”  In _Journal of Computational Physics_ 292, 2015, pp. 168–175  DOI: 10.1016/j.jcp.2015.03.003
  * [18] Michael B. Giles  “Multilevel Monte Carlo methods”  In _Acta Numerica_ 24 Cambridge University Press, 2015, pp. 259–328  DOI: 10.1017/S096249291500001X
  * [19] Michael B. Giles  “Multilevel Monte Carlo Path Simulation”  In _Operations Research_ 56.3 INFORMS, 2008, pp. 607–617  DOI: 10.1287/opre.1070.0496
  * [20] Alex A. Gorodetsky, Gianluca Geraci, Michael S. Eldred and John D. Jakeman  “A Generalized Approximate Control Variate Framework for Multifidelity Uncertainty Quantification”  In _Journal of Computational Physics_ 408, 2020, pp. 109257  DOI: 10.1016/j.jcp.2020.109257
  * [21] I.. Gradshteyn and I.. Ryshik  “Table of Integrals, Series, and Products”  Elsevier, 2015  URL: https://linkinghub.elsevier.com/retrieve/pii/C20100648395
  * [22] I.. Graham, F.. Kuo, D. Nuyens, R. Scheichl and I.. Sloan  “Analysis of Circulant Embedding Methods for Sampling Stationary Random Fields”  In _SIAM Journal on Numerical Analysis_ 56.3, 2018, pp. 1871–1895  DOI: 10.1137/17M1149730
  * [23] I.G. Graham, F.Y. Kuo, D. Nuyens, R. Scheichl and I.H. Sloan  “Quasi-Monte Carlo Methods for Elliptic PDEs with Random Coefficients and Applications”  In _Journal of Computational Physics_ 230.10, 2011, pp. 3668–3694  DOI: 10.1016/j.jcp.2011.01.023
  * [24] A. Gregory, C.. Cotter and S. Reich  “Multilevel Ensemble Transform Particle Filtering”  In _SIAM Journal on Scientific Computing_ 38.3 Society for Industrial and Applied Mathematics, 2016, pp. A1317–A1338  DOI: 10.1137/15M1038232
  * [25] P. Guttorp and T. Gneiting  “Studies in the history of probability and statistics XLIX: On the Matérn correlation family”  In _Biometrika_ 93, 2006, pp. 989–995  DOI: 10.1093/biomet/93.4.989
  * [26] Wolfgang Hackbusch  “Multi-Grid Methods and Applications”  Berlin Heidelberg: Springer-Verlag, 1985  DOI: 10.1007/978-3-662-02427-0
  * [27] Eric Hallman, Ilse C.. Ipsen and Arvind K. Saibaba  “Monte Carlo Methods for Estimating the Diagonal of a Real Symmetric Matrix”  In _SIAM Journal on Matrix Analysis and Applications_ 44.1 Society for Industrial and Applied Mathematics, 2023, pp. 240–269  DOI: 10.1137/22M1476277
  * [28] R.V.L. Hartley  “A More Symmetrical Fourier Analysis Applied to Transmission Problems”  In _Proceedings of the IRE_ 30.3, 1942, pp. 144–150  DOI: 10.1109/JRPROC.1942.234333
  * [29] S. Heinrich  “Monte Carlo Complexity of Global Solution of Integral Equations”  In _Journal of Complexity_ 14.2, 1998, pp. 151–175  DOI: 10.1006/jcom.1998.0471
  * [30] Håkon Hoel, Kody J.. Law and Raul Tempone  “Multilevel ensemble Kalman filtering”  In _SIAM Journal on Numerical Analysis_ 54.3 Society for IndustrialApplied Mathematics, 2016, pp. 1813–1839  DOI: 10.1137/15M100955X
  * [31] Anastasia Istratuca and Aretha Teckentrup  “Smoothed Circulant Embedding with Applications to Multilevel Monte Carlo Methods for PDEs with Random Coefficients”  arXiv, 2023  DOI: 10.48550/arXiv.2306.13493
  * [32] Prashant Kumar, Cornelis W. Oosterlee and Richard P. Dwight  “A multigrid multilevel Monte Carlo method using high-order finite-volume scheme for lognormal diffusion problems”  In _International Journal for Uncertainty Quantification_ 7.1 Begel House Inc., 2017  DOI: 10.1615/Int.J.UncertaintyQuantification.2016018677
  * [33] Prashant Kumar, Carmen Rodrigo, Francisco J. Gaspar and Cornelis W. Oosterlee  “On Local Fourier Analysis of Multigrid Methods for PDEs with Jumping and Random Coefficients”  In _SIAM Journal on Scientific Computing_ 41.3 Society for Industrial and Applied Mathematics, 2019, pp. A1385–A1413  DOI: 10.1137/18M1173769
  * [34] O.. Le Maître and Omar M. Knio  “Spectral Methods for Uncertainty Quantification”, Scientific Computation  Dordrecht: Springer Netherlands, 2010  DOI: 10.1007/978-90-481-3520-2
  * [35] Christiane Lemieux  “Monte Carlo and Quasi-Monte Carlo Sampling”, Springer Series in Statistics  New York, NY: Springer, 2009  DOI: 10.1007/978-0-387-78165-5
  * [36] Finn Lindgren, Håvard Rue and Johan Lindström  “An explicit link between Gaussian fields and Gaussian Markov random fields: the stochastic partial differential equation approach”  In _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ 73.4, 2011, pp. 423–498  DOI: 10.1111/j.1467-9868.2011.00777.x
  * [37] Nora Lüthen, Stefano Marelli and Bruno Sudret  “Sparse Polynomial Chaos Expansions: Literature Survey and Benchmark”  In _SIAM/ASA Journal on Uncertainty Quantification_ 9.2 Society for Industrial and Applied Mathematics, 2021, pp. 593–649  DOI: 10.1137/20M1315774
  * [38] Aimee Maurais, Terrence Alsup, Benjamin Peherstorfer and Youssef Marzouk  “Multi-Fidelity Covariance Estimation in the Log-Euclidean Geometry”  arXiv, 2023  DOI: 10.48550/arXiv.2301.13749
  * [39] Aimee Maurais, Terrence Alsup, Benjamin Peherstorfer and Youssef Marzouk  “Multifidelity Covariance Estimation via Regression on the Manifold of Symmetric Positive Definite Matrices”  arXiv, 2023  DOI: 10.48550/arXiv.2307.12438
  * [40] M.. McKay, R.. Beckman and W.. Conover  “A Comparison of Three Methods for Selecting Values of Input Variables in the Analysis of Output from a Computer Code”  In _Technometrics_ 21.2 [Taylor & Francis, Ltd., American Statistical Association, American Society for Quality], 1979, pp. 239–245  DOI: 10.2307/1268522
  * [41] Paul Mycek and Matthias De Lozzo  “Multilevel Monte Carlo Covariance Estimation for the Computation of Sobol’ Indices”  In _SIAM/ASA Journal on Uncertainty Quantification_ 7.4 Society for IndustrialApplied Mathematics, 2019, pp. 1323–1348  DOI: 10.1137/18M1216389
  * [42] Habib N. Najm  “Uncertainty Quantification and Polynomial Chaos Techniques in Computational Fluid Dynamics”  In _Annual Review of Fluid Mechanics_ 41.1, 2009, pp. 35–52  DOI: 10.1146/annurev.fluid.010908.165248
  * [43] Benjamin Peherstorfer, Karen Willcox and Max Gunzburger  “Survey of Multifidelity Methods in Uncertainty Propagation, Inference, and Optimization”  In _SIAM Review_ 60.3 Society for IndustrialApplied Mathematics, 2018, pp. 550–591  DOI: 10.1137/16M1082469