# A filtered multilevel Monte Carlo method for estimating the expectation of cell-centered discretized random fields

Jérémy Briant, Paul Mycek, Mayeul Destouches, Olivier Goux, Serge Gratton, Selime Gürol, Ehouarn Simon, Anthony T. Weaver

► To cite this version:

Jérémy Briant, Paul Mycek, Mayeul Destouches, Olivier Goux, Serge Gratton, et al. A filtered multilevel Monte Carlo method for estimating the expectation of cell-centered discretized random fields. International Journal for Uncertainty Quantification, In press, 10.1615/Int.J.UncertaintyQuantification.2025057754. hal-04337114v3

HAL Id: hal-04337114

https://hal.science/hal-04337114v3

Submitted on 11 Jul 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
# A filtered multilevel Monte Carlo method for estimating the expectation of cell-centered discretized random fields*

Jérémy Briant¹,³, Paul Mycek²,³,⁴, Mayeul Destouches⁵, Olivier Goux²,³, Serge Gratton¹,⁶, Selime Gürol²,³, Ehouarn Simon¹, and Anthony T. Weaver²,³

¹IRIT, UMR 5505, Université de Toulouse, CNRS, Toulouse INP, Toulouse, France (jeremy.briant@toulouse-inp.fr, serge.gratton@toulouse-inp.fr, ehouarn.simon@toulouse-inp.fr).

²Cerfacs, Toulouse, France (mycek@cerfacs.fr, goux@cerfacs.fr, gurol@cerfacs.fr, weaver@cerfacs.fr).

³CECI, UMR 5318, Université de Toulouse, Cerfacs, CNRS, IRD, Toulouse, France.

⁴Concace, Airbus CR&T, Cerfacs, Inria, France.

⁵CNRM, UMR 3589, Université de Toulouse, Météo-France, CNRS, Toulouse, France (mayeul destouches@meteo.fr).

⁶ANITI, Toulouse, France.

## 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 numerical simulators have inconsistent dimensions across the multilevel hierarchy. This motivates the introduction of grid transfer operators borrowed from multigrid methods. By adapting mathematical tools from 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. We then propose filtered MLMC (F-MLMC) estimators based on a filtering mechanism similar to the smoothing process of multigrid methods, and we show that 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. Next, the conclusions of the spectral analysis are experimentally verified with a one-dimensional illustration. Finally, the proposed 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 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.

*CC-BY 4.0 licence.* 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).
# 1 Introduction

In recent years, multilevel and multifidelity Monte Carlo (MC) methods (see, e.g., [1, 2] for surveys on these methods) have grown in popularity as a means to accelerate traditional MC methods by combining simulators of different fidelities to leverage the lower computational cost of lower-fidelity simulators while preserving the accuracy of the high-fidelity simulator. Such lower fidelity numerical simulators can be obtained, e.g., from a coarser spatial and/or temporal discretization as is the case in multigrid methods [3–6]. In such instances, the different fidelities are referred to as levels, and the multifidelity method is then referred to as a multilevel method. Among multilevel methods, the multi-level Monte Carlo (MLMC) method [1, 7, 8], which is the focus of the present paper, combines samples (or ensembles) from different levels in such a way that, under certain assumptions, the variance of the resulting multilevel estimator is reduced, while leaving the bias unchanged. 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 [9, 10] and variance-based global sensitivity measures [11]. The analysis of MLMC estimators has also been extended to the estimation of statistics of random variables with values in separable Hilbert spaces [9, 10].

The estimation of covariance matrices is a prominent example where multilevel approaches have started to emerge, especially in the field of data assimilation [12–16]. 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 [17]. Diffusion operators are commonly used for modelling spatial covariances in ocean data assimilation [18] and are closely related to other techniques for modelling spatial covariances in atmospheric data assimilation [19], geostatistical modelling [20], inverse problems [21] and uncertainty quantification [22]. 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 [18, 23]. 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 cell-centered discretized random fields in the abstract setting described above. The main specificities of this setting are that, first, the considered simulators are based on a hierarchy of cell-centered discretizations on grids of different resolution; second, the input and output of the numerical simulators are discretized fields whose dimensions depend on the level; and third, we focus on the discrete representation (as opposed to a functional representation) of these discretized fields, arising from the discretization of the strong form of the differential operators at hand (e.g., in our example, diffusion operators). These specificities motivate the introduction of grid transfer operators, namely restriction and prolongation operators, which may be borrowed from multigrid methods [3–6].

Using MLMC with simulators whose inputs and/or outputs are discretized fields has been studied before [24–33], typically in a functional framework based on the finite element (FE) or finite volume (FV) method, which prompts the use of natural grid transfer operators, especially prolongation operators, which are then usually only implicitly defined. Therefore, their impact on the quality of the resulting multilevel estimator is generally not studied. Nonetheless, the introduction of grid transfer operators raises crucial questions about their effect on the different spatial scales (or, equivalently, frequencies) of the discretized fields (or signals) that are transferred between grids of different resolution. We note that the MLMC estimator proposed in [33] circumvents the need to use restriction operators
by explicitly exploiting the joint distribution of the pair of input fields defined on successive levels within each MLMC correction term, which can be exhibited for the specific problem addressed in [33].

The question of the representation of small scales on coarse grids within the MLMC framework has been explored previously by [30, 32] in an FE setting where the input of the considered simulator is a discretized random field. The authors proposed the idea of smoothing highly oscillatory input signals based on a spectral truncation of their high-frequency components, which is controlled for each MLMC level so that the truncation error matches the discretization error. In this paper, we investigate these important questions in the abstract setting described above. Our objective is to gain a better understanding of the impact of grid transfer operators on the quality of the MLMC estimator in order to propose remedies for their negative effects. This leads us to propose a novel, filtered MLMC (F-MLMC) estimator that mitigates the negative effects of the grid transfer operators on the quality of the multilevel estimator. In multigrid methods [3–6, 34], studying the effects of grid transfer operators on the different frequencies of a discrete signal, which is typically achieved through a spectral analysis, is core to determining the effectiveness of the method. Along the same lines, we present a spectral analysis of the MLMC and F-MLMC estimators in a simplified setting where the numerical simulators are assumed to be linear, symmetric and circulant operators.

The remainder of this paper is organized as follows. Section 2 introduces the MLMC and F-MLMC estimators of the expectation of random fields as adaptations of the framework in [9] to our specific setting, relying on (possibly filtered) grid transfer operators. A spectral analysis of the MLMC and F-MLMC estimators is then carried out in section 3 to investigate the effects of grid transfer on their variance, in the specific case of linear, symmetric and circulant numerical simulators. In section 4, we apply the MLMC and F-MLMC estimators, first, to a one-dimensional illustration that comply with the assumptions of the spectral analysis, and, second, to our target problem of estimating the intrinsic variances of a two-dimensional, heterogeneous, diffusion-based covariance operator. General conclusions are drawn in section 5, along with prospective avenues for future work.

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

The MLMC method aims to improve the accuracy of MC estimators by combining samples of different fidelities. In favorable cases, large low-fidelity samples (i.e., with many members that are cheap to generate on coarse levels) are used to reduce the sampling error (variance), while smaller samples are required at finer and more expensive levels to correct the bias. Under certain assumptions, [24, Theorem 1] ensures that there exists a sample allocation 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.

In this paper, we focus on the multilevel estimation of the expectation of cell-centered discretized random fields. We assume that such discretized random fields can be represented as random variables defined on some probability space $(\Omega, \Sigma, \mathbb{P})$, with values in separable Hilbert spaces and having finite second-order moment. We thus rely on the mathematical framework proposed in [9] for the design and analysis of MC and MLMC estimators in such a setting. Specifically, given 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 on $(\Omega, \Sigma, \mathbb{P})$,

$$L^2(\Omega, H) := \{\xi: \Omega \to H \mid \int_{\Omega} \|\xi(\omega)\|_H^2 d\mathbb{P}(\omega) < +\infty\}, \quad (1)$$
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 d\mathbb{P}(\omega), \qquad (2) $$

whose induced norm is denoted by $\|\cdot\|_{L^2(\Omega, H)}$. More specifically, we are interested in discretized random fields that correspond to the output of a numerical simulator whose input is also a discretized random field. Formally, we consider an abstract, deterministic numerical simulator as $f: H_{\text{in}} \to H_{\text{out}}$, where $H_{\text{in}}$ and $H_{\text{out}}$ are two separable Hilbert spaces, and we assume that $f$ is such that $f(X) \in L^2(\Omega, H_{\text{out}})$ for any $X \in L^2(\Omega, H_{\text{in}})$. Then, given an input $X \in L^2(\Omega, H_{\text{in}})$, the Bochner integral $\mathbb{E}[f(X)] := \int_{\Omega} f(X(\omega))d\mathbb{P}(\omega)$, which we shall refer to as the expectation of the output $f(X)$, is well-defined and we have $\mathbb{E}[f(X)] \in H_{\text{out}}$ [35, 36]. From here on, for the sake of exposition and without loss of generality, we restrict ourselves to the case where $H_{\text{in}} = H_{\text{out}} = H$.

We now turn to the multilevel estimation of $\mathbb{E}[f(X)]$, which relies on a hierarchy (sequence) of $L+1$ numerical simulators ($f_\ell: H \to H$)$_{\ell=0}^L$ of increasing fidelity. We are interested in estimating the expectation of the output of the finest (highest-fidelity) model, $\mathbb{E}[f_L(X)]$. The MLMC estimator $\hat{\mu}_L^{\text{MLMC}}$ of $\mathbb{E}[f_L(X)]$, using a sequence ($X^{(\ell,i)}$)$_{i=1}^{M_\ell}$)$_{\ell=0}^L$ of $L+1$ independent samples of $X$, is defined as

$$ \hat{\mu}_L^{\text{MLMC}} = \frac{1}{M_0} \sum_{i=1}^{M_0} f_0(X^{(0,i)}) + \sum_{\ell=1}^{L} \frac{1}{M_\ell} \sum_{i=1}^{M_\ell} [f_\ell(X^{(\ell,i)}) - f_{\ell-1}(X^{(\ell,i)})] \in L^2(\Omega, H). \qquad (3) $$

The accuracy of this MLMC estimator may be quantified by its normwise MSE with respect to some $\mu \in H$, defined by (see, e.g., [9])

$$ \text{MSE}(\hat{\mu}_L^{\text{MLMC}}, \mu) := \|\hat{\mu}_L^{\text{MLMC}} - \mu\|_{L^2(\Omega, H)}^2. \qquad (4) $$

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

$$ \text{MSE}(\hat{\mu}_L^{\text{MLMC}}, \mu) = \mathcal{V}(\hat{\mu}_L^{\text{MLMC}}) + \|\mathbb{E}[f_L(X)] - \mu\|_H^2, \qquad (5) $$

where

$$ \forall Y \in L^2(\Omega, H), \quad \mathcal{V}(Y) := \|Y - \mathbb{E}[Y]\|_{L^2(\Omega, H)}^2 = \mathbb{E}[\|Y\|_H^2] - \|\mathbb{E}[Y]\|_H^2, \qquad (6) $$

with the shorthand notation $\mathbb{E}[\| \cdot \|_H^2]^{1/2} = \| \cdot \|_{L^2(\Omega, H)}$. The first term in eq. (5) is referred to as the variance of the multilevel estimator, while the second term corresponds to the squared bias. Furthermore, [9, Theorem 3.1] shows that the variance $\mathcal{V}(\hat{\mu}_L^{\text{MLMC}})$ can be further decomposed level-wise
into

$$ \mathcal{V}(\hat{\mu}_L^{\text{MLMC}}) = \frac{1}{M_0} \mathcal{V}(f_0(X)) + \sum_{\ell=1}^{L} \frac{1}{M_\ell} \mathcal{V}(f_\ell(X) - f_{\ell-1}(X)). \quad (7) $$

We note that the multilevel estimator is unbiased with respect to the expectation of the output at the finest level $\mu_L := \mathbb{E}[f_L(X)]$, so that $\text{MSE}(\hat{\mu}_L^{\text{MLMC}}, \mu_L) = \mathcal{V}(\hat{\mu}_L^{\text{MLMC}})$.

## 2.2 MLMC with inconsistent input and output spaces

In typical applications, especially when the simulators are related to discretized partial differential equations (PDEs), and that their fidelity is dictated by the quality of the discretization, the input and output spaces $H_\ell$ are typically finite-dimensional and generally not the same across the simulators of the hierarchy. The multilevel hierarchy of numerical simulators then becomes $(\tilde{f}_\ell: H_\ell \to H_\ell)_{\ell=0}^L$, and we are interested in the expectation of the highest-fidelity output, $\mathbb{E}[\tilde{f}_L(X_L)] \in H_L$, for some input $X_L \in L^2(\Omega, H_L)$. As a consequence, the MLMC estimator defined in eq. (3) cannot be used directly, but transfer operators may be introduced to transfer inputs and outputs between the highest-fidelity space $H_L$ and lower-fidelity spaces $H_\ell$ with $\ell = 0, \dots, L-1$. Operators that transfer a discretized field from a level $\ell$ to a higher-fidelity level $\ell' > \ell$ will be referred to as prolongation operators. We denote by $P_\ell^{\ell'}: H_\ell \to H_{\ell'}$ the prolongation operator from level $\ell$ to a higher-fidelity level $\ell' > \ell$. In addition, restriction operators $R_\ell^{\ell'}: H_{\ell'} \to H_\ell$ are used to perform the “reverse” operation of transferring discretized fields defined on level $\ell'$ to a lower-fidelity level $\ell < \ell'$. We may now define $f_\ell$ from $\tilde{f}_\ell$ as

$$ f_\ell = P_\ell^L \circ \tilde{f}_\ell \circ R_L^\ell: H_L \to H_L, \quad (8) $$

for $\ell = 0, \dots, L$.

The numerical simulators in eq. (8) involve transfer operators between the highest-fidelity level $L$ and lower-fidelity levels. In practice, grid transfer operators between arbitrary levels can be defined as the composition of transfer operators between successive levels, as

$$ \begin{aligned} \forall \ell = 0, \dots, L-1, \quad \forall \ell' = \ell + 1, \dots, L, \quad P_\ell^{\ell'} &:= P_{\ell'-1}^{\ell'} \circ P_{\ell'-2}^{\ell'-1} \circ \dots \circ P_{\ell+1}^{\ell+2} \circ P_\ell^{\ell+1}, \\ R_{\ell'}^{\ell} &:= R_{\ell+1}^{\ell} \circ R_{\ell+2}^{\ell+1} \circ \dots \circ R_{\ell'-1}^{\ell'-2} \circ R_{\ell'}^{\ell'-1}, \end{aligned} \quad (9) $$

and, for $\ell = 0, \dots, L$, $R_\ell^\ell = P_\ell^\ell = I_\ell$, where $I_\ell: H_\ell \to H_\ell$ is the identity operator on $H_\ell$. We further assume that, for any $X \in L^2(\Omega, H_\ell)$,

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

* $P_\ell^{\ell+1}(X) \in L^2(\Omega, H_{\ell+1})$, for $\ell = 0, \dots, L-1$;

* $R_\ell^{\ell-1}(X) \in L^2(\Omega, H_{\ell-1})$, for $\ell = 1, \dots, L$.

It follows that, for any $X_L \in L^2(\Omega, H_L)$ and for $\ell = 0, \dots, L$, $f_\ell(X_L) \in L^2(\Omega, H_L)$. Consequently, the MLMC expression eq. (3) can now be used with $H = H_L$ and $X = X_L \in L^2(\Omega, H_L)$. The resulting MLMC formulation, i.e., based on transferred versions $f_\ell$ of the simulators $\tilde{f}_\ell$ as defined in eq. (8) and composed transfer operators as in eq. (9), constitutes the core of the proposed methodology, and will be the focus of the remainder of this paper, in a discrete, finite-dimensional setting described in the next section.
## 2.3 Discrete finite-dimensional spaces

In this paper, we focus on numerical simulators that arise from the discrete approximation of the strong form of elliptic PDEs, i.e., through the direct discretization of the differential operators involved in such PDEs. This is the case, for instance, in the finite difference (FD) or finite volume (FV) discretization of a PDE. This choice is motivated by the considered application, presented in section 4.1, and more importantly by the target ocean numerical model [37], which relies on an FD-like discretization. Consequently, input and output discretized fields of such simulators may only be evaluated at a finite number of discrete locations, and are thus represented as finite-dimensional real vectors, whose size depends on the number of discretization points in the spatial domain of interest $\mathcal{D}$. Formally, we consider a polytopal tessellation $\mathcal{T}$ of a polytopal domain $\mathcal{D} \subset \mathbb{R}^d$, with $d \in \{1, 2, 3\}$. We refer to elements of $\mathcal{T}$ as cells, and we denote by $|\mathcal{T}|$ the cardinality, or size, of $\mathcal{T}$, i.e., its number of cells.

We now describe the MLMC hierarchy of simulators (and corresponding Hilbert spaces) for this discretization setting. Let $(\mathcal{T}_\ell)_{\ell=0}^L$ be a sequence of $L+1$ tessellations of $\mathcal{D}$ of increasing size $n_\ell := |\mathcal{T}_\ell|$, for $\ell = 0, \dots, L$. These tessellations of different sizes define the geometric hierarchy of levels of the MLMC approach. Specifically, $\mathcal{T}_0$ is the coarsest tessellation corresponding to the lowest-fidelity level, while $\mathcal{T}_L$ is the finest tessellation corresponding to the highest-fidelity level. To each level $\ell = 0, \dots, L$, we associate the symmetric positive definite (SPD) Gram matrix $\mathbf{W}_\ell \in \mathbb{R}^{n_\ell \times n_\ell}$, which encodes the structural information related to the discrete approximation of scalar fields as vectors of $\mathbb{R}^{n_\ell}$ on the corresponding cell-centered discretization on $\mathcal{T}_\ell$. We then 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}}, \quad (10) $$

where $\mathbf{V}_\ell \in \mathbb{R}^{n_\ell \times n_\ell}$ is an invertible matrix such that $\mathbf{W}_\ell = \mathbf{V}_\ell \mathbf{V}_\ell^\mathrm{T}$, which exists since $\mathbf{W}_\ell$ is SPD, and where $\langle \cdot, \cdot \rangle_{\mathbf{I}_{n_\ell}}$ denotes the canonical (Euclidean) dot product between vectors of $\mathbb{R}^{n_\ell}$. The inner product space $H_\ell := (\mathbb{R}^{n_\ell}, \langle \cdot, \cdot \rangle_{\mathbf{W}_\ell})$ is a separable Hilbert space, so that the multilevel hierarchy of simulators ($\tilde{f}_\ell: H_\ell \to H_\ell$) $\mathcal{L}_{\ell=0}^L$ fits in the framework described in section 2.2. For clarity, we shall from here on simply write $\mathbb{R}^{n_\ell}$ for $(\mathbb{R}^{n_\ell}, \langle \cdot, \cdot \rangle_{\mathbf{W}_\ell})$, and, unless explicitly stated otherwise, use $\langle \cdot, \cdot \rangle_{\mathbf{W}_\ell}$ as the default inner product between vectors of $\mathbb{R}^{n_\ell}$. Furthermore, we shall denote by $\| \cdot \|_{\mathbf{W}_\ell}$ the norm induced by $\langle \cdot, \cdot \rangle_{\mathbf{W}_\ell}$. In this discrete setting, the MLMC estimator defined in eq. (3) becomes

$$ \hat{\mu}_L^{\text{MLMC}} = \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} [f_\ell(\mathbf{X}_L^{(\ell,i)}) - f_{\ell-1}(\mathbf{X}_L^{(\ell,i)})], \quad (11) $$

with the grid transferred version of the simulators defined as in eq. (8). It should be remarked that the $L+1$ independent samples $(\{\mathbf{X}_L^{(\ell,i)}\}_{i=1}^{M_\ell})_{\ell=0}^L$ of $\mathbf{X}_L$ are all generated on the finest grid.

In the functional framework that is typically employed for simulators whose inputs and/or outputs are discretized fields [24–32], grid transfer operators are implicitly employed, but usually not explicitly defined or even mentioned. This is especially the case for the prolongation operator, since discretized fields are defined everywhere in the spatial domain through their functional representation. Nevertheless, the practical computation of the multilevel correction terms in eq. (11) implicitly requires the introduction of prolongation operators. By focusing on the discrete framework, we explicitly exhibit the grid transfer operators and study their effect on the quality of the multilevel estimator. We shall see that the operators that seem the most natural in the functional framework, such as the canonical injection for the prolongation operator in settings with nested approximation spaces, actually yield multilevel estimators with poor properties. This is particularly true for cell-centered discretizations, which are the focus of this paper (see section 2.5), and which are also known in the multigrid community for causing similar difficulties [38–41].
## 2.4 The filtered MLMC method

In this paper, we are particularly concerned with the impact of the grid transfer operators on the representation of the different scales (or frequencies) of the multilevel correction fields. While it is clear that, owing to the telescoping correction mechanism, the bias of the multilevel estimator (11) is unaffected, the impact on its variance, and hence on the MSE, remains to be investigated. As will be illustrated in sections 3 and 4, the introduction of grid transfer operators induces a deterioration of the multilevel estimator's variance in the high frequencies. Inspired by multigrid methods [5, 6, 34], 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. The idea of filtering highly oscillatory input fields in an MLMC framework has been proposed in [30, 32], for simulators with scalar outputs. This filtering is achieved through spectral truncation in the Karhunen-Loève [42–44] or circulant embedding [45, 46] framework.

In this paper, we propose a general filtering formalism for the abstract, grid-transferred MLMC framework (11). Because the outputs of the simulators considered here are discretized fields, we also apply the filtering process to the outputs of coarser levels after they have been prolonged to the fine grid. In the proposed formalism, filtering the small-scale components out of a signal defined on level $ℓ$ is achieved using a low-pass pre- or post-filtering operator $S_{ℓ}^{\text{pre/post}}: \mathbb{R}^{n_{ℓ}} \to \mathbb{R}^{n_{ℓ}}$. These filters can then be composed with the grid transfer operators between successive levels to obtain pre-filtered restriction operators and post-filtered prolongation operators,

$$ \forall \ell = 1, \dots, L, \quad \bar{R}_{\ell}^{\ell-1} := R_{\ell}^{\ell-1} \circ S_{\ell}^{\text{pre}}, \quad \text{and} \quad \bar{P}_{\ell-1}^{\ell} := S_{\ell}^{\text{post}} \circ P_{\ell-1}^{\ell}. \qquad (12) $$

Similar to eq. (9), filtered grid transfer operators between arbitrary levels can be defined as the composition of transfer operators between successive levels, as

$$ \begin{aligned} \forall \ell = 0, \dots, L-1, \quad \forall \ell' = \ell + 1, \dots, L, \quad \bar{P}_{\ell}^{\ell'} &:= \bar{P}_{\ell'-1}^{\ell'} \circ \bar{P}_{\ell'-2}^{\ell'-1} \circ \dots \circ \bar{P}_{\ell+1}^{\ell+2} \circ \bar{P}_{\ell}^{\ell+1}, \\ \bar{R}_{\ell}^{\ell'} &:= \bar{R}_{\ell+1}^{\ell} \circ \bar{R}_{\ell+2}^{\ell+1} \circ \dots \circ \bar{R}_{\ell'-1}^{\ell'-2} \circ \bar{R}_{\ell'}^{\ell'-1}, \end{aligned} \qquad (13) $$

and, for $ℓ = 0, \dots, L$, $\bar{R}_{\ell}^{\ell} = \bar{P}_{\ell}^{\ell} = I_{\ell}$. A new (filtered) simulator $\bar{f}_{\ell}$ is then defined on each level $ℓ$ as

$$ \bar{f}_{\ell} = \bar{P}_{\ell}^{L} \circ \tilde{f}_{\ell} \circ \bar{R}_{L}^{\ell}, \qquad (14) $$

so that the filtered MLMC (F-MLMC) estimator is obtained by replacing $f_{\ell}$ with $\bar{f}_{\ell}$ in eq. (11):

$$ \hat{\mu}_{L}^{\text{F-MLMC}} = \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)})]. \qquad (15) $$

## 2.5 Practical considerations

While the proposed methodology solely relies on eqs. (8) and (11) for the unfiltered MLMC estimator, and on eqs. (12) to (15) for the F-MLMC estimator, we describe in this section further practical considerations that will be assumed for the remainder of this paper.

### 2.5.1 Discrete setting and MLMC hierarchy

In what follows, we consider cell-centered discretizations, meaning that the discretization points are located at the center of the cells of a tessellation $\mathcal{T}$. It should be noted that functional representations can be deduced a posteriori from such discretized fields. In particular, if $\mathbf{u} \in \mathbb{R}^n$ is a cell-centered
discretized field on a tessellation $\mathcal{T} := \{T_i\}_{i=1}^n$, a natural functional representation is a piecewise constant field $u \in \mathbb{P}^0$, where $\mathbb{P}^0 := \text{span}(\{1_{T_i}\}_{i=1}^n)$ and $1_T: \mathcal{D} \to \{0, 1\}$ denotes the indicator function associated with the subset $T \subseteq \mathcal{D}$. The piecewise constant representation can then be defined through $u \in \mathbb{R}^n$ as $u = u^T \phi$, where $\phi := [1_{T_1}, \dots, 1_{T_n}]^T$.

In the remainder of this paper, we prescribe a fixed finest level $L$. This is consistent with the fact that, in complex operational settings, the high-fidelity level is typically dictated by the application. We further assume that the tessellations of the grid hierarchy $(\mathcal{T}_\ell)_{\ell=0}^L$ are nested, i.e., for $\ell = 1, \dots, L$, for any $T \in \mathcal{T}_\ell$, there is a unique $T' \in \mathcal{T}_{\ell-1}$ such that $T \subset T'$. We note, however, that the sets of discretization points, i.e., the sets of cell centers, are not nested. It should be stressed that the nestedness of the tessellations is not strictly required, but facilitates the design of transfer operators. Non-nested grid hierarchies would call for the design of more complex, and possibly more expensive, transfer operators.

## 2.5.2 Optimal sample allocation

The sample sizes $M_\ell$ are optimally allocated across levels such that, for a given computational budget $C$, the variance of the multilevel estimator is minimal [11],

$$ M_{\ell} = \left\lfloor \frac{C}{S_{L}} \sqrt{\frac{\mathcal{V}_{\ell}}{C_{\ell} + C_{\ell-1}}} \right\rfloor^{+}, \quad \mathcal{V}_{\ell} := \mathcal{V}(\tilde{f}_{\ell}(\mathbf{X}_{L}) - \tilde{f}_{\ell-1}(\mathbf{X}_{L})), \quad S_{L} := \sum_{\ell=0}^{L} \sqrt{\mathcal{V}_{\ell}(C_{\ell} + C_{\ell-1})}, \quad (16) $$

where $C_\ell$ denotes the mean computational cost of evaluating $\tilde{f}_\ell(\mathbf{X}_L)$ and $\|\cdot\|^+ := \max(1, \|\cdot\|)$, where $\|\cdot\|$ denotes the floor function. By convention, the quantities indexed by $\ell = -1$ vanish, so that $C_{-1} = 0$ and $\mathcal{V}_0 = \mathcal{V}(\tilde{f}_0(\mathbf{X}_L))$. Furthermore, for conciseness, we have used the same notation $\tilde{f}_\ell$ to refer either to the filtered version of $f_\ell$ (i.e., $\tilde{f}_\ell$ as defined in section 2.4) when considering F-MLMC, or to its unfiltered counterpart (i.e., simply $f_\ell$) when dealing with unfiltered MLMC. As a matter of fact, the latter is a special case of the former, with identity operators as filters. For a fixed finest level $L$, it is easy to see from the optimal sample allocation (16) (ignoring the rounding of $M_\ell$) that the variance of the resulting multilevel estimator is $S_L^2/C$, while the variance of the single-level, high-fidelity MC estimator at equivalent cost is $C_L \mathcal{V}(f_L(\mathbf{X}_L))/C$. This shows that, while the variance of the multilevel estimator is reduced (regardless of the budget) provided that $S_L^2 < C_L \mathcal{V}(f_L(\mathbf{X}_L))$, the convergence rate remains unchanged compared to a standard single-level MC estimator, which is typical of multifidelity methods where $L$ is fixed (see, e.g., [47]). It should be noted that, because of the introduction of grid transfer operations, the level-wise cost $C_\ell$ can be decomposed into two contributions, specifically, the cost $\tilde{C}_\ell$ of evaluating the simulator $\tilde{f}_\ell$ on the one hand, and the cost $C_\ell^{\text{transf}}$ of applying the (possibly filtered) grid transfer operators on the other hand.

## 2.5.3 Unfiltered grid transfer operators

In practical applications, the transfer operators introduced in section 2.2 are typically chosen to be linear. A linear prolongation operator $P_\ell^{\ell'}$ may be identified with the matrix $\mathbf{P}_\ell^{\ell'} \in \mathbb{R}^{n_{\ell'} \times n_\ell}$ such that $P_\ell^{\ell'}(\mathbf{x}_\ell) = \mathbf{P}_\ell^{\ell'} \mathbf{x}_\ell$ for any $\mathbf{x}_\ell \in \mathbb{R}^{n_\ell}$, and, likewise, a linear restriction operator $R_\ell^{\ell'}$ may be identified with the appropriate matrix $\mathbf{R}_\ell^{\ell'} \in \mathbb{R}^{n_\ell \times n_{\ell'}}$. In the remainder of this paper, we shall only consider linear transfer operators, and hence employ the associated matrix notation. Then, eq. (9) implies

$$ \forall \ell = 0, \dots, L-1, \quad \mathbf{P}_\ell^L := \mathbf{P}_{L-1}^L \cdots \mathbf{P}_{\ell+1}^{\ell+2} \mathbf{P}_\ell^{\ell+1} \quad \text{and} \quad \mathbf{R}_\ell^L := \mathbf{R}_{\ell+1}^\ell \cdots \mathbf{R}_{L-1}^{L-2} \mathbf{R}_L^{L-1}, \quad (17) $$
and $\mathbf{R}_L^L = \mathbf{P}_L^L = \mathbf{I}_{n_L}$. Consequently, the MLMC estimator (11) can be written as

$$ \hat{\mu}_L^{\text{MLMC}} = \frac{1}{M_0} \sum_{i=1}^{M_0} \mathbf{P}_0^L \tilde{f}_0(\mathbf{R}_L^0 \mathbf{X}_L^{(0,i)}) + \sum_{\ell=1}^{L} \frac{1}{M_\ell} \sum_{i=1}^{M_\ell} \left[ \mathbf{P}_\ell^L \tilde{f}_\ell(\mathbf{R}_L^\ell \mathbf{X}_L^{(\ell,i)}) - \mathbf{P}_{\ell-1}^L \tilde{f}_{\ell-1}(\mathbf{R}_L^{\ell-1} \mathbf{X}_L^{(\ell,i)}) \right], \quad (18) $$

$$ = \frac{1}{M_0} \mathbf{P}_0^L \sum_{i=1}^{M_0} \tilde{f}_0(\mathbf{R}_L^0 \mathbf{X}_L^{(0,i)}) + \sum_{\ell=1}^{L} \frac{1}{M_\ell} \mathbf{P}_\ell^L \sum_{i=1}^{M_\ell} \left[ \tilde{f}_\ell(\mathbf{R}_L^\ell \mathbf{X}_L^{(\ell,i)}) - \mathbf{P}_{\ell-1}^\ell \tilde{f}_{\ell-1}(\mathbf{R}_L^{\ell-1} \mathbf{X}_L^{(\ell,i)}) \right], \quad (19) $$

where the second identity is obtained by exploiting the linearity of the prolongation operators and the composition definition in eq. (17). The latter is to be preferred in practical implementations, as it involves fewer applications of the prolongation operators. Finally, in the remainder of this paper, we shall consider prolongation operators satisfying

$$ \forall \ell = 0, \dots, L-1, \quad (\mathbf{P}_\ell^{\ell+1})^\mathrm{T} \mathbf{W}_{\ell+1} \mathbf{P}_\ell^{\ell+1} = \mathbf{W}_\ell. \quad (20) $$

A major consequence of eq. (20) is that

$$ \forall \ell = 0, \dots, L-1, \quad \forall \mathbf{x}_\ell \in \mathbb{R}^{n_\ell}, \quad \|\mathbf{P}_\ell^{\ell+1} \mathbf{x}_\ell\|_{\mathbf{W}_{\ell+1}} = \mathbf{x}_\ell^\mathrm{T} (\mathbf{P}_\ell^{\ell+1})^\mathrm{T} \mathbf{W}_{\ell+1} \mathbf{P}_\ell^{\ell+1} \mathbf{x}_\ell = \mathbf{x}_\ell^\mathrm{T} \mathbf{W}_\ell \mathbf{x}_\ell = \|\mathbf{x}_\ell\|_{\mathbf{W}_\ell}, \quad (21) $$

indicating that the prolongation operator is norm-preserving. We note that, by eq. (17), similar properties then hold for transfer operators between arbitrary levels.

### 2.5.4 Filtering operators and filtered grid transfer operators

Without loss of generality, we shall use the same operator $S_\ell = S_\ell^{\text{pre}} = S_\ell^{\text{post}}$ for both pre- and post-filtering on each level. Moreover, we shall use linear filtering operators $S_\ell$, which can thus be identified with matrices $\mathbf{S}_\ell$, so that the filtered grid transfer operators remain linear. The resulting filtered grid transfer operators $\bar{\mathbf{P}}_\ell^L$ and $\bar{\mathbf{R}}_\ell^L$ are then defined through their corresponding matrices $\bar{\mathbf{P}}_\ell^L$ and $\bar{\mathbf{R}}_\ell^L$ by

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

where

$$ \forall \ell = 1, \dots, L, \quad \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 (23) $$

and $\bar{\mathbf{R}}_L^L = \bar{\mathbf{P}}_L^L = \mathbf{I}_{n_L}$. The F-MLMC estimator then reads

$$ \hat{\mu}_L^{\text{F-MLMC}} = \frac{1}{M_0} \sum_{i=1}^{M_0} \bar{\mathbf{P}}_0^L \tilde{f}_0(\bar{\mathbf{R}}_L^0 \mathbf{X}_L^{(0,i)}) + \sum_{\ell=1}^{L} \frac{1}{M_\ell} \sum_{i=1}^{M_\ell} \left[ \bar{\mathbf{P}}_\ell^L \tilde{f}_\ell(\bar{\mathbf{R}}_L^\ell \mathbf{X}_L^{(\ell,i)}) - \bar{\mathbf{P}}_{\ell-1}^L \tilde{f}_{\ell-1}(\bar{\mathbf{R}}_L^{\ell-1} \mathbf{X}_L^{(\ell,i)}) \right], \quad (24) $$

$$ = \frac{1}{M_0} \bar{\mathbf{P}}_0^L \sum_{i=1}^{M_0} \tilde{f}_0(\bar{\mathbf{R}}_L^0 \mathbf{X}_L^{(0,i)}) + \sum_{\ell=1}^{L} \frac{1}{M_\ell} \bar{\mathbf{P}}_\ell^L \sum_{i=1}^{M_\ell} \left[ \tilde{f}_\ell(\bar{\mathbf{R}}_L^\ell \mathbf{X}_L^{(\ell,i)}) - \bar{\mathbf{P}}_{\ell-1}^\ell \tilde{f}_{\ell-1}(\bar{\mathbf{R}}_L^{\ell-1} \mathbf{X}_L^{(\ell,i)}) \right]. \quad (25) $$

Again, we note that the unfiltered MLMC estimator corresponds to a special case of the filtered MLMC estimator with identity filtering operators. However, we note that, in general, the filtered transfer operators defined in section 2.4 no longer satisfy the norm-preserving property of eq. (20). It should be stressed that the proposed methodology remains valid with non-linear (filtered or unfiltered) transfer operators; in particular, the telescoping sum property still guarantees the unbiasedness of the multilevel estimators in eqs. (11) and (15).
# 3 Spectral analysis

In this section, we conduct a spectral analysis of the MLMC estimator (18) and of the F-MLMC estimator (24) to study more closely the effects of grid transfer operators on their variance at different spatial scales. The analysis is performed in a one-dimensional spatial domain, with cell-centered discretized input and output fields on uniform, nested grids. The setting is further simplified to the case where the numerical simulators are linear, and that their associated matrices are symmetric and circulant.

## 3.1 Linear simulators

We start by considering linear simulators of the form $\tilde{f}_{\ell}: \mathbf{x}_{\ell} \mapsto \tilde{\mathbf{F}}_{\ell}\mathbf{x}_{\ell}$, for $\ell = 0, \dots, L$, where $\tilde{\mathbf{F}}_{\ell} \in \mathbb{R}^{n_{\ell} \times n_{\ell}}$, and, consistently with eqs. (8) and (18), we define $\mathbf{F}_{\ell} := \mathbf{P}_{\ell}^{L} \tilde{\mathbf{F}}_{\ell} \mathbf{R}_{L}^{\ell} \in \mathbb{R}^{n_{L} \times n_{L}}$, 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}$. Then, by linearity of the expectation operator, the variance $\mathcal{V}(\hat{\mu}_{L}^{\text{MLMC}})$ of the multilevel estimator (18), given by eq. (7), becomes

$$ \mathcal{V}(\hat{\mu}_{L}^{\text{MLMC}}) = \frac{1}{M_0} \mathbb{E}[\|\mathbf{F}_0 \dot{\mathbf{X}}_L\|_{\mathbf{W}_L}^2] + \sum_{\ell=1}^{L} \frac{1}{M_{\ell}} \mathbb{E}[\|(\mathbf{F}_{\ell} - \mathbf{F}_{\ell-1}) \dot{\mathbf{X}}_L\|_{\mathbf{W}_L}^2], \quad (26) $$

where $\dot{\mathbf{X}}_L := \mathbf{X}_L - \mathbb{E}[\mathbf{X}_L]$. Now, for any matrix $\mathbf{F} \in \mathbb{R}^{n_L \times n_L}$,

$$ \begin{aligned} \mathbb{E}[\|\mathbf{F}\dot{\mathbf{X}}_L\|_{\mathbf{W}_L}^2] &= \mathbb{E}[(\mathbf{F}\dot{\mathbf{X}}_L)^{\mathrm{T}}\mathbf{W}_L(\mathbf{F}\dot{\mathbf{X}}_L)] = \mathrm{tr}\mathbb{E}[\dot{\mathbf{X}}_L^{\mathrm{T}}\mathbf{F}^{\mathrm{T}}\mathbf{W}_L\mathbf{F}\dot{\mathbf{X}}_L] = \mathbb{E}[\mathrm{tr}(\dot{\mathbf{X}}_L^{\mathrm{T}}\mathbf{F}^{\mathrm{T}}\mathbf{W}_L\mathbf{F}\dot{\mathbf{X}}_L)] \\ &= \mathbb{E}[\mathrm{tr}(\mathbf{F}^{\mathrm{T}}\mathbf{W}_L\mathbf{F}\dot{\mathbf{X}}_L\dot{\mathbf{X}}_L^{\mathrm{T}})] = \mathrm{tr}(\mathbf{F}^{\mathrm{T}}\mathbf{W}_L\mathbf{F}\mathbb{E}[\dot{\mathbf{X}}_L\dot{\mathbf{X}}_L^{\mathrm{T}}]) = \mathrm{tr}(\mathbf{F}^{\mathrm{T}}\mathbf{W}_L\mathbf{F}\mathbf{G}) \end{aligned} \quad (27) \\ = \mathrm{tr}((\mathbf{F}\mathbf{G}^{1/2})^{\mathrm{T}}\mathbf{W}_L(\mathbf{F}\mathbf{G}^{1/2})) = \|\mathbf{F}\mathbf{G}^{1/2}\|_{F,\mathbf{W}_L}^2 \end{aligned} $$

where $\mathrm{tr}(\cdot)$ denotes the matrix trace operator, $\| \cdot \|_{F,\mathbf{W}}: \mathbf{F} \mapsto \mathrm{tr}(\mathbf{F}^{\mathrm{T}}\mathbf{W}\mathbf{F})^{1/2}$ denotes the $\mathbf{W}$-weighted Frobenius norm for any SPD weighting matrix $\mathbf{W}$, and $\mathbf{G} := \mathbb{E}[\dot{\mathbf{X}}_L\dot{\mathbf{X}}_L^{\mathrm{T}}] = \mathbf{G}^{1/2}(\mathbf{G}^{1/2})^{\mathrm{T}}$ denotes the covariance matrix of $\mathbf{X}_L$. Thus, eq. (26) reduces to

$$ \mathcal{V}(\hat{\mu}_{L}^{\text{MLMC}}) = \frac{1}{M_0} \|\mathbf{F}_0 \mathbf{G}^{1/2}\|_{F, \mathbf{W}_L}^2 + \sum_{\ell=1}^{L} \frac{1}{M_{\ell}} \|(\mathbf{F}_{\ell} - \mathbf{F}_{\ell-1}) \mathbf{G}^{1/2}\|_{F, \mathbf{W}_L}^2, \quad (28) $$

which emphasizes that the variance reduction is closely related to the similarity of successive fidelity models. Furthermore, for any orthogonal matrix $\mathbf{Q} \in \mathbb{R}^{n_L \times n_L}$ with respect to the weighted inner product $\langle \cdot, \cdot \rangle_{\mathbf{W}_L}$, i.e., such that $\mathbf{Q}\mathbf{W}_L\mathbf{Q}^{\mathrm{T}} = \mathbf{Q}^{\mathrm{T}}\mathbf{W}_L\mathbf{Q} = \mathbf{I}_{n_L}$, we have, for any $\mathbf{x} \in \mathbb{R}^{n_L}$,

$$ \|\mathbf{Q}^{\mathrm{T}}\mathbf{V}^{\mathrm{T}}\mathbf{x}\|_{\mathbf{W}_L} = \|\mathbf{Q}\mathbf{V}^{\mathrm{T}}\mathbf{x}\|_{\mathbf{W}_L} = \|\mathbf{x}\|_{\mathbf{W}_L}, \quad \text{and} \quad \|\mathbf{Q}^{\mathrm{T}}\mathbf{x}\|_{\mathbf{W}_L} = \|\mathbf{Q}\mathbf{x}\|_{\mathbf{W}_L} = \|\mathbf{x}\|_{\mathbf{I}_{n_L}}, \quad (29) $$

and, similarly, for any $\mathbf{F} \in \mathbb{R}^{n_L \times n_L}$,

$$ \|\mathbf{Q}^{\mathrm{T}}\mathbf{V}^{\mathrm{T}}\mathbf{F}\|_{F, \mathbf{W}_L} = \|\mathbf{Q}\mathbf{V}^{\mathrm{T}}\mathbf{F}\|_{F, \mathbf{W}_L} = \|\mathbf{F}\|_{F, \mathbf{W}_L}, \quad \text{and} \quad \|\mathbf{Q}^{\mathrm{T}}\mathbf{F}\|_{F, \mathbf{W}_L} = \|\mathbf{Q}\mathbf{F}\|_{F, \mathbf{W}_L} = \|\mathbf{F}\|_{F, \mathbf{I}_{n_L}}. \quad (30) $$

## 3.2 One-dimensional spatial setting

The analysis is conducted in a one-dimensional (1D) setting with periodic, cell-centered discretized input and output random fields defined on a unit spatial domain $\mathcal{D} := [0, 1) \subset \mathbb{R}$. Specifically, we consider a hierarchy of uniform grids $(\mathcal{T}_{\ell})_{\ell=0}^{L}$, where $\mathcal{T}_{\ell} := \{T_{\ell,i}\}_{i=1}^{n_{\ell}}$ and $T_{\ell,i} := [(i-1)/n_{\ell}, i/n_{\ell})$ on level $\ell$, along with the associated Gram matrices $\mathbf{W}_{\ell} = n_{\ell}^{-1}\mathbf{I}_{n_{\ell}}$. It follows immediately that
<figure><img src="image_2.png" /><figcaption>Figure 1: The vectors $\mathbf{h}_k^0$ (orange) of a coarse-grid Hartley basis on level $\ell = 0$, with $n_0 = 8$, and the vectors $\mathbf{h}_k^1$ (blue) of a fine-grid Hartley basis on level $\ell = 1$, with $n_1 = 16$. For $\ell \in \{0, 1\}$ and for each vector $k = 0, \dots, n_\ell - 1$, the entries of the vector $[\mathbf{h}_k^\ell]_j$ are plotted against $x_j := (2j+1)/(2n_\ell) \in (0, 1)$, for $j = 0, \dots, n_\ell - 1$.</figcaption></figure>

Prolongating a coarse-grid basis vector $\mathbf{h}_k^0$, 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}_k^1$ and $\mathbf{h}_{n_0+k}^1$. Figure 1 shows that, for $k \le n_0/2$, the fine-grid signal $\mathbf{h}_k^1$ has the same frequency as the original, coarse-grid signal $\mathbf{h}_k^0$, while $\mathbf{h}_{n_0+k}^1$ has higher frequency. Conversely, for $k > n_0/2$, the fine-grid signal $\mathbf{h}_{n_0+k}^1$ has the same frequency as the coarse-grid signal $\mathbf{h}_k^0$, while $\mathbf{h}_k^1$ has higher frequency. In both cases, the prolongation of a coarse-grid signal introduces spurious high-frequency (i.e., small-scale) components to the prolonged signal. Fortunately, both $\mathbf{h}_k^1$ and $\mathbf{h}_{n_0+k}^1$ 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}_{n_0+k}^1$ than the consistent, low-frequency signals $\mathbf{h}_k^1$. 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}_k^1$ than the consistent, low-frequency signals $\mathbf{h}_{n_0+k}^1$.

For the restriction operator $\mathbf{R} := \mathbf{R}_1^0$ defined in eq. (31), we have (see section C)

$$ \mathbf{R}\mathbf{h}_k^1 = 2c_k\mathbf{h}_k^0, \quad \text{and} \quad \mathbf{R}\mathbf{h}_{n_0+k}^1 = -2c_{n_0+k}\mathbf{h}_k^0, \quad \forall k = 0, \dots, n_0 - 1. \qquad (35) $$

Restricting a fine-grid basis vector $\mathbf{h}_k^1$, with $k = 0, \dots, n_0 - 1$, produces a (coarse-grid) vector proportional to the corresponding coarse-grid basis vector $\mathbf{h}_k^0$, specifically by a factor $2c_k \le 2$. For $n_0/2 < k < n_0$, the restricted signal has lower frequency than the original, fine-grid signal, as illustrates fig. 1. Similarly, restricting a fine-grid basis vector $\mathbf{h}_k^1$, with $k = n_0, \dots, 2n_0 - 1$, produces a (coarse-grid) vector proportional to the complementary coarse-grid basis vector $\mathbf{h}_{k-n_0}^0$, specifically by a factor $-2c_k < 2$. Again, high-frequency fine-grid signals $\mathbf{h}_k^1$ corresponding to $n_0 \le 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} := [\text{Diag}(\{c_k\}_{k=0}^{n_0-1}) \quad \text{Diag}(\{-c_{n_0+k}\}_{k=0}^{n_0-1})] \in \mathbb{R}^{n_0 \times 2n_0}$, the identities in eqs. (33) and (35) may be compactly recast as

$$ \mathbf{PH}_0 = \mathbf{H}_1\mathbf{C}^T, \quad \mathbf{RH}_1 = 2\mathbf{H}_0\mathbf{C}, \quad \mathbf{H}_1^T\mathbf{P} = 2\mathbf{C}^T\mathbf{H}_0^T \quad \mathbf{H}_0^T\mathbf{R} = \mathbf{C}\mathbf{H}_1^T, \qquad (36) $$

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

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

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

$$ \tilde{\mathbf{F}}_\ell = \mathbf{H}_\ell \mathbf{\Lambda}_\ell \mathbf{H}_\ell^T, \quad \mathbf{\Lambda}_\ell := \text{Diag}(\{\lambda_k^\ell\}_{k=0}^{n_\ell-1}), \quad \text{for } \ell \in \{0, 1\}. \qquad (37) $$

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

$$ \mathbf{M}_{11} := \text{Diag}(\{c_k^2\lambda_k^0\}_{k=0}^{n_0-1}), \qquad (38) $$

$$ \mathbf{M}_{22} := \text{Diag}(\{c_{n_0+k}^2\lambda_k^0\}_{k=0}^{n_0-1}), \qquad (39) $$

$$ \mathbf{M}_{12} := \text{Diag}(\{-c_k c_{n_0+k} \lambda_k^0\}_{k=0}^{n_0-1}) = \mathbf{M}_{21}. \qquad (40) $$

As a consequence, eq. (28) becomes

$$ \mathcal{V}(\hat{\mu}_1^{\text{MLMC}}) = \frac{1}{M_0} \|\mathbf{H}_1 \mathbf{M} \mathbf{H}_1^T \mathbf{G}^{1/2}\|_{F, \mathbf{W}_L}^2 + \frac{1}{M_1} \|\mathbf{H}_1 (\mathbf{\Lambda}_1 - \mathbf{M}) \mathbf{H}_1^T \mathbf{G}^{1/2}\|_{F, \mathbf{W}_L}^2, \qquad (41) $$

which, owing to the orthogonality of the Hartley matrix, and more particularly eq. (30), can be recast as

$$ \mathcal{V}(\hat{\mu}_1^{\text{MLMC}}) = \frac{1}{M_0} \|\mathbf{M}\mathbf{H}_1^T\mathbf{G}^{1/2}\|_{F, \mathbf{I}_{n_L}}^2 + \frac{1}{M_1} \|(\mathbf{\Lambda}_1 - \mathbf{M})\mathbf{H}_1^T\mathbf{G}^{1/2}\|_{F, \mathbf{I}_{n_L}}^2. \qquad (42) $$

To reduce the variance of the correction term in eq. (42), the difference between $\mathbf{M}$ and $\mathbf{\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 $\mathbf{\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{F}_0$ from the fine-grid operator $\tilde{F}_1$, and which is widely used in multigrid methods and their analysis [5]. Specifically, the Galerkin operator is defined as

$$ \tilde{F}_0 := \frac{1}{4} \mathbf{R} \tilde{F}_1 \mathbf{P} \in \mathbb{R}^{n_0 \times n_0}. \quad (43) $$

As shown in section E, this operator is the optimal operator in terms of minimizing $\|\tilde{F}_1 - \mathbf{P}\tilde{F}_0\|_{F,\mathbf{W}_1}^2$ for $\mathbf{W}_1 = n_1^{-1}\mathbf{I}_{n_1}$ and the grid transfer operators defined by eq. (31). It follows from eqs. (36) and (43) that the Galerkin operator $\tilde{F}_0$ can be diagonalized in the Hartley basis as $\tilde{F}_0 = \mathbf{H}_0\mathbf{\Lambda}_0\mathbf{H}_0^T$ with $\mathbf{\Lambda}_0 = \mathbf{C}\mathbf{\Lambda}_1\mathbf{C}^T$. In other words, the eigenvalues $\{\lambda_k^0\}_{k=0}^{n_0-1}$ of the Galerkin operator $\tilde{F}_0$ can be expressed from the eigenvalues $\{\lambda_k^1\}_{k=0}^{n_0-1}$ of $\tilde{F}_1$,

$$ \lambda_k^0 = c_k^2 \lambda_k^1 + c_{n_0+k}^2 \lambda_{n_0+k}^1, \quad \forall k = 0, \dots, n_0-1. \quad (44) $$

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

$$ \mathbf{M}_{11} = \mathrm{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}), \quad (45) $$

$$ \mathbf{M}_{22} = \mathrm{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}) = \mathrm{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}), \quad (46) $$

$$ \mathbf{M}_{12} = \mathrm{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}. \quad (47) $$

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 $\mathbf{\Lambda}_1$ by a multiplicative damping factor $c_k^4 \le 1$ for $k = 0, \dots, 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, \dots, 2n_0 - 1$. The off-diagonal blocks $\mathbf{M}_{12}$ and $\mathbf{M}_{21}$ introduce spurious terms that contribute to the difference $\mathbf{\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{F}_1$ is of little interest, which is why we simply consider the damping factors $-c_k c_{n_0+k}$ with respect to eigenvalues $\lambda_k^0$ of the coarse-grid operator $\tilde{F}_0$, for $k = 0, \dots, n_0 - 1$, given by eq. (40).

The evolution with $k = 0, \dots, n_1 - 1$ of the four damping factors is presented in fig. 2(a). We observe that the eigenvalues $\lambda_k^1$ associated with low-frequency vectors of the fine Hartley basis, i.e., for $k$ close to 0 and $n_1 - 1$ (see fig. 1), 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 $\mathbf{\Lambda}_1 - \mathbf{M}$. In particular, we note that the first entry of $\mathbf{M}_{11}$ exactly matches that of $\mathbf{\Lambda}_1$, i.e., $\lambda_0^1$, because $c_0 = 1$ and $c_{n_0} = 0$. On the other hand, eigenvalues $\lambda_k^1$ associated with medium- to high-frequency vectors of the fine Hartley basis 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 factors $-c_k^3 c_{n_0+k}$ and $-c_k c_{n_0+k}^3$, whose magnitudes are maximal for $k \in \{n_1/8 + 1, 7n_1/8 - 1\}$ and $k \in \{3n_1/8 - 1, 5n_1/8 + 1\}$.

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

We now turn to the F-MLMC estimator based on the second-order Shapiro filter [51–53] 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} \in \mathbb{R}^{n_{\ell} \times n_{\ell}}. \quad (48) $$
<figure><img src="image_3.png" /><figcaption>Figure 2: Damping factors of the eigenvalues in M and $\bar{M}$ as functions of $k = 0, \dots, n_1 - 1$ with $n_1 = 2n_0 = 32$, when using the Galerkin operator eq. (43). 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 and green curves correspond to the factors of the off-diagonal blocks.</figcaption></figure>

With the specific grid transfer operators defined in eq. (31), the operator $\bar{P}_{\ell-1}^\ell$ corresponds to the linear interpolation operator between the levels $\ell-1$ and $\ell$. We again restrict ourselves to a two-level analysis and define the filtered grid transfer operators as $\bar{P} = S_1P$ and $\bar{R} = RS_1$, where $P := P_0^1$ and $R := R_1^0$ as in section 3.3, and where $S_1$ denotes the second-order Shapiro filter defined in eq. (48). Similarly to eqs. (33) and (35), it is possible to study the effect the filtered transfer operators on the Hartley basis vectors (see section F). For the prolongation $\bar{P}$, we have

$$ \bar{P}h_k^0 = c_k^3 h_k^1 - c_{n_0+k}^3 h_{n_0+k}^1, \quad \forall k = 0, \dots, n_0 - 1. \qquad (49) $$

The addition of the Shapiro filter raises the damping factors $c_k$ to the power 3. The prolonged Hartley basis vectors are thus more severely damped than they were without filtering. Again, the most damped fine-grid basis vectors $h_k^1$ are those corresponding to $k$ close to $n_0 = n_1/2$, i.e., high-frequency signals. For the restriction operator we have

$$ \bar{R}h_k^1 = 2c_k^3h_k^0, \quad \text{and} \quad \bar{R}h_{n_0+k}^1 = -2c_{n_0+k}^3h_k^0, \quad \forall k = 0, \dots, n_0 - 1. \qquad (50) $$

Similar conclusions can be drawn as for the unfiltered case, but with the damping factors raised to the power of 3, thus increasing their effect, 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. (49) and (50) can be recast as

$$ \bar{P}H_0 = H_1C_3^T, \quad \bar{R}H_1 = 2H_0C_3, \quad H_1^T\bar{P} = 2C_3^T H_0^T, \quad H_0^T\bar{R} = C_3H_1^T, \qquad (51) $$

where $C_3 := [\text{Diag}(\{c_k^3\}_{k=0}^{n_0-1}) \text{Diag}(\{-c_{n_0+k}^3\}_{k=0}^{n_0-1})] \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{F}_1$ and $\tilde{F}_0$ are symmetric, circulant matrices. From eqs. (37) and (51), we deduce the decomposition $\bar{F}_0 := \bar{P}\bar{F}_0\bar{R} = H_1\bar{M}\bar{H}_1^T$, where $\bar{M} := \begin{bmatrix} \bar{M}_{11} & \bar{M}_{12} \\ \bar{M}_{21} & \bar{M}_{22} \end{bmatrix}$, with

$$ \bar{M}_{11} = \text{Diag}(\{c_k^6\lambda_k^0\}_{k=0}^{n_0-1}), \qquad (52) $$

$$ \bar{M}_{22} = \text{Diag}(\{c_{n_0+k}^6\lambda_k^0\}_{k=0}^{n_0-1}), \qquad (53) $$

$$ \bar{M}_{12} = \text{Diag}(\{-c_k^3 c_{n_0+k}^3 \lambda_k^0\}_{k=0}^{n_0-1}) = \bar{M}_{21}. \qquad (54) $$
Upon replacing $\mathbf{F}_0$ with its filtered counterpart $\tilde{\mathbf{F}}_0$, eq. (28) can be written as in eq. (42), but with $\mathbf{M}$ replaced by $\tilde{\mathbf{M}}$. The sparsity pattern of $\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 $\tilde{\mathbf{M}}_{12}$ and $\tilde{\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. 2 (orange plots), which contributes to reducing the off-diagonal entries of $\mathbf{\Lambda}_1 - \tilde{\mathbf{M}}$ compared to those of $\mathbf{\Lambda}_1 - \mathbf{M}$.

To study more closely the diagonal entries of $\tilde{\mathbf{M}}$ and compare them to the eigenvalues of $\tilde{\mathbf{F}}_1$ in $\mathbf{\Lambda}_1$ we resort to the Galerkin operator defined by $\tilde{\mathbf{F}}_0 := \frac{1}{4}\mathbf{R}\tilde{\mathbf{F}}_1\mathbf{P}$. This definition is inspired by the form of eq. (43), although there is no guarantee of its optimality. Then, from eqs. (37) and (51), 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\mathbf{\Lambda}_0\mathbf{H}_0^T$ with $\mathbf{\Lambda}_0 = \mathbf{C}_3\mathbf{\Lambda}_1\mathbf{C}_3^T$, or, equivalently,

$$ \lambda_k^0 = c_k^6 \lambda_k^1 + c_{n_0+k}^6 \lambda_{n_0+k}^1, \quad \forall k = 0, \dots, n_0-1. \qquad (55) $$

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

$$ \tilde{\mathbf{M}}_{11} = \mathrm{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}), \qquad (56) $$

$$ \tilde{\mathbf{M}}_{22} = \mathrm{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}) = \mathrm{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}), \qquad (57) $$

$$ \tilde{\mathbf{M}}_{12} = \mathrm{Diag}(\{-c_k^9 c_{n_0+k}^3 \lambda_k^1 - c_k^3 c_{n_0+k}^9 \lambda_{n_0+k}^1\}_{k=0}^{n_0-1}) = \tilde{\mathbf{M}}_{21}. \qquad (58) $$

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. 2 (blue plots). Moreover, the factors $c_k^{12}$ also damp more strongly the consistent eigenvalues compared to the factors $c_k^4$ of $\mathbf{M}$, as can be visualized in fig. 2 (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.

# 4 Numerical experiments

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 [18, 23]. This particular problem, which is briefly outlined below, will be the focus of our numerical experiments.

## 4.1 Problem description

Let $u: \mathcal{D} \to \mathbb{R}$ and $b: \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, \qquad (59) $$

where $m$ is a positive integer, $I$ is the identity operator, and $\mathbf{K}: \mathcal{D} \to \mathbb{R}^{d \times d}$ is a symmetric, positive-definite (SPD) tensor field with entries $[K_{ij}]_{i,j=1,...,d}$. If $\mathbf{K}$ is constant, the integral solution on $\mathbb{R}^d$
corresponds to the application of a covariance operator, whose kernel is a covariance function from the Matérn class, to b [54, 55].

We assume that the operator in eq. (59) 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. (59) as

$$L := (I - \Delta)^{-m} W^{-1}, \quad (60)$$

where $\Delta$ is the matrix representing a spatial discretization of the differential operator $\nabla \cdot K\nabla$, and $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, $W$ is such that $\Delta$ is self-adjoint with respect to $\langle \cdot, \cdot \rangle_W$, i.e., $W\Delta = \Delta^T W$. Consequently, the matrix $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}(x) = (2m - d - 2)^{-1}(D_{ii}(x))^2$ where the elements $[D_{ii}(x)]_{i=1,...,d}$ correspond to the directional correlation length-scales at the spatial location $x$, and $m > d/2 + 1$ [56, section 3].

The matrix $L$ is SPD but does not define a covariance matrix with meaningful variances for applications like data assimilation. As such, $L$ must be normalized by its diagonal so that the desired variances can be applied. Thus, we define the covariance matrix of interest as $B = \Sigma\Gamma L\Gamma\Sigma$, where $\Gamma = \text{Diag}(\text{diag}(L))^{-1/2}$ is a normalizing diagonal matrix such that $\text{Diag}(\text{diag}(\Gamma L\Gamma)) = I_n$, and $\Sigma^2 = \Sigma\Sigma$ is the diagonal matrix with entries corresponding to the desired variances, i.e., $\text{diag}(B) = \text{diag}(\Sigma^2) = (\sigma_1^2, ..., \sigma_n^2)$. In these expressions, the operator $\text{diag}(\cdot)$ maps a matrix to the vector consisting of the diagonal elements of that matrix, while the operator $\text{Diag}(\cdot)$ maps a vector to the diagonal matrix whose diagonal consists of the entries of that vector.

The main computational challenge lies in computing the diagonal $\theta := \text{diag}(L) = (\theta_1, ..., \theta_n)$. Indeed, for large-scale problems, the matrix $L$ is not assembled, and only applications of $L$ to vectors are accessible. Thus, its diagonal entries are not explicitly stored and need to be determined differently [23]. A direct way would be to recover these by applying $L$ to the canonical basis vectors of $\mathbb{R}^n$, i.e., $\theta_k = (Le_k)_k$, for $k = 1, ..., n$. For large $n$, this approach is not computationally tractable. Here, we focus on an alternative strategy, used in operational systems, consisting of approximating $\theta$ by randomization [23]. Taking $m = 2q$, $L$ can be subsequently factored (see section G) as $L = AWA^T$, where

$$A := (I - \Delta)^{-q} W^{-1}, \quad (61)$$

which, as for $L$ itself, cannot be explicitly assembled in large-scale applications. This decomposition of $L$ implies that, for any random vector $X$ with $E[X] = 0_n$ and $E[XX^T] = W$,

$$\theta = \text{diag}(L) = \text{diag}(A E[XX^T] A^T) = \text{diag}(C[AX]) = V[AX] = E[AX \odot AX], \quad (62)$$

where $E[\cdot]$ and $V[\cdot]$ denote the element-wise expectation and variance of a random vector, $C[\cdot]$ denotes the covariance matrix of a random vector, and where $\odot$ denotes the Schur product (a.k.a. the Hadamard or element-wise product) between two vectors or matrices of the same size. One way to construct such a vector $X$ is to define $X = VZ$, where $V$ arises from the factorization $W = VV^T$, and where $Z$ is a random vector with $E[Z] = 0_n$ and $E[ZZ^T] = I_n$. In the rest of this paper, we consider normal random vectors $X \sim N(0_n, W)$, which can be constructed by $X = VZ$, with $Z \sim N(0_n, I_n)$. The expectation $\theta$ of the $\mathbb{R}^n$-valued random vector $Y := AX \odot AX$ can then be estimated using MC sampling. Given a random $M$-sample $\{X^{(i)}\}_{i=0}^M$ of $X$, an unbiased estimator $\hat{\theta}$ of $\theta$ is the sample mean

$$\hat{\theta} = \frac{1}{M} \sum_{i=1}^{M} (AX^{(i)} \odot AX^{(i)}). \quad (63)$$
We remark that the estimator $\hat{\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. (63) 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 [57, 58], which do not employ a factored form of the matrix and hence do not guarantee non-negative estimates. Here, we investigate the use of the MLMC and F-MLMC estimators described in sections 2 and 2.4 to improve (in terms of RMSE) the estimation of $\theta$ and hence the efficiency of determining accurate normalization coefficients $\Gamma$. It should be noted that the proposed (F)-MLMC construction does not guarantee non-negative estimates either, although negative estimates were never encountered in our experiments.

In practice, the matrix-vector product $\mathbf{A}\mathbf{x}$ can be computed by solving the sequence of SPD systems of linear equations $(\mathbf{W}(\mathbf{I} - \mathbf{\Delta}))\mathbf{y}^{(j)} = \mathbf{z}^{(j)})_{j=1}^q$, with $\mathbf{z}^{(1)} := \mathbf{x}$ and $\mathbf{z}^{(j)} := \mathbf{W}\mathbf{y}^{(j-1)}$ for $j = 2, \dots, q$, so that $\mathbf{y}^{(q)} = \mathbf{A}\mathbf{x}$. In our experiments, the numerical solving of these systems is achieved by precomputing a Cholesky decomposition of $\mathbf{W}(\mathbf{I} - \mathbf{\Delta})$. However, for large-scale problems, it is not reasonably possible to compute or store such a decomposition, and sparse iterative methods are typically used instead [59, 60].

## 4.2 1D illustration

Before turning to the (F-)MLMC estimation of $\theta$, we briefly consider a simpler experiment to verify numerically the conclusions of the spectral analysis of section 3. Specifically, we consider a hierarchy of simulators $(f_\ell)_{\ell=0}^L$ defined through eq. (8) by $\tilde{f}_\ell: \mathbf{x}_\ell \mapsto \mathbf{A}_\ell\mathbf{x}_\ell$, where $\mathbf{A}_\ell$ is the factor in eq. (61) arising from a cell-centered discretization on the 1D domain $\mathcal{D} := [0, 1]$, with periodic boundary conditions. Specifically, we consider the same setting and grid transfer operators as in section 3.2, and we define hierarchies of different depths (i.e., numbers of levels), corresponding to $L \in \{0, \dots, 5\}$, with a fixed finest discretization corresponding to $n_L = 512$. The diffusion tensor field $\mathbf{K}$, which reduces here to a scalar field, is taken to be constant by setting $D_{11}(\mathbf{x}) = D \in \mathbb{R}$ for all $\mathbf{x} \in \mathcal{D}$. The scalar value $D$ will be referred to as the length-scale. We note that in this setting, the matrices $\mathbf{A}_\ell$ are symmetric and circulant [61].

The statistical parameter of interest is $\mathbb{E}[f_L(\mathbf{X}_L)] = \mathbb{E}[\mathbf{A}_L\mathbf{X}_L]$, with $\mathbf{X}_L \sim \mathcal{N}(\mathbf{0}_{n_L}, \mathbf{W}_L)$. It trivially follows that $\mathbb{E}[f_L(\mathbf{X}_L)] = \mathbf{0}_{n_L}$. The MLMC estimator $\hat{\mu}_L^{\text{MLMC}}$ is defined as in eqs. (18) and (19), where the transfer operators are defined as in eq. (31). With this choice, it follows from eqs. (17) and (20) that the MLMC estimator $\hat{\mu}_L^{\text{MLMC}}$ involves, at each correction level, input random vectors that are such that $\mathbf{R}_L^0\mathbf{X}_L \sim \mathcal{N}(\mathbf{0}_{n_0}, \mathbf{W}_0)$ and $[(\mathbf{R}_L^\ell\mathbf{X}_L)^\text{T}}(\mathbf{R}_L^{\ell-1}\mathbf{X}_L)^\text{T}]^\text{T} \sim \mathcal{N}(\mathbf{0}_{n_\ell+n_{\ell-1}}, \tilde{\mathbf{W}}_\ell)$ for $\ell = 1, \dots, L$, with

$$ \tilde{\mathbf{W}}_\ell := \begin{bmatrix} \mathbf{W}_\ell & \mathbf{W}_\ell \mathbf{P}_{\ell-1}^\ell \\ \mathbf{R}_\ell^{\ell-1} \mathbf{W}_\ell & \mathbf{W}_{\ell-1} \end{bmatrix} \in \mathbb{R}^{(n_\ell+n_{\ell-1}) \times (n_\ell+n_{\ell-1})}, \quad (64) $$

which has the same structure as [33, eq. (3.12)]. In fact, eq. (64) and [33, eq. (3.12)] are equivalent under the functional, piecewise constant representation of the input fields. The filtered estimator $\hat{\mu}_L^{\text{F-MLMC}}$ is defined as in eq. (15), with the same transfer operators and the 2nd-order Shapiro filters defined in eq. (48).

We now describe the cost model that will be prescribed in the remainder of this paper. First, the cost $\tilde{C}_\ell$ of evaluating $\tilde{f}_\ell$ is assumed to be proportional to the number of cells $n_\ell$ on level $\ell$, i.e., $\tilde{C}_\ell = \alpha n_\ell$. This is consistent with a fixed number $q \times T$ of matrix-free applications of $\mathbf{A}_\ell$, where $q$ represents the number of diffusion steps in eq. (61), i.e., the number of linear systems to be solved, and $T$ the fixed number of iterations used to solve each of these systems. Second, the cost $C_\ell^{\text{transf}}$ of applying the (possibly filtered) grid transfer operators is modeled as $C_\ell^{\text{transf}} = \beta \sum_{k=\ell+1}^L n_k$, where $\beta$ denotes the constant (level-independent) number of floating-point operations per cell on level $k$ for applying the (possibly filtered) restriction operator from level $k$ to level $k-1$. This model is consistent with the fact
that the restriction from the finest grid is defined as the composition of restriction operators between successive levels (see eqs. (17) and (22)). Exploiting the nestedness of the tessellations, and assuming a constant refinement factor $s := n_\ell/n_{\ell-1} > 1$ between two successive grids, the restriction cost reduces to $C_\ell^{\text{transf}} = \frac{s}{s-1}\beta(s^{L-\ell}-1)n_\ell$. Eventually, the cost $C_\ell = \tilde{C}_\ell + C_\ell^{\text{transf}}$ of evaluating $f_\ell$ is given by

$$C_{\ell} = s^{\ell-L} \left[ 1 + \gamma \left( s^{L-\ell} - 1 \right) \right] C_{L}, \quad \text{with } \gamma := \frac{s}{s-1} \frac{\beta}{\alpha} > 0. \qquad (65)$$

It is then easy to see from eq. (65) that $C_\ell/C_{\ell-1} \to s$ as $\gamma \to 0$.

In the target application, Chebyshev iterations (see, e.g., [62–64]) are typically used [60]. We note that each iteration involves one application of the operator $(\mathbf{I}_{n_\ell} - \mathbf{\Delta}_\ell)$, which, for a second-order discretization of the 1D Laplace operator represented by a 3-point stencil, requires 3 floating-point operations per cell, hence $\alpha = 3qT$. In what follows, we choose $m = 2q = 10$ and $T = 20$, which are typical values [60], so that $\alpha = 300$. For the unfiltered MLMC, the 1D restriction, as defined in eq. (31), from a given fine level $k$ to the immediately coarser level $k-1$ requires $s = 2$ operations per coarse cell, i.e., a total of $sn_{k-1} = n_k$ operations, hence $\beta = 1$. The 1D pre-filtering defined in eq. (48), represented by a 3-point stencil, requires 3 additional operations on level $k$, so that $\beta = 4$ for the filtered MLMC. Table 1 summarizes the normalized costs $C_\ell/C_L$ and cost reduction factors $C_\ell/C_{\ell-1}$ of the simulators $f_\ell$ (MLMC) and $\tilde{f}_\ell$ (F-MLMC) for this 1D example.

<table><caption>Table 1: Normalized costs $C_\ell/C_L$ and cost reduction factors $C_\ell/C_{\ell-1}$ of the simulators $f_\ell$ (MLMC) and $\tilde{f}_\ell$ (F-MLMC) corresponding to the 1D example of section 4.2.</caption><thead><tr><th rowspan="2">ℓ</th><th colspan="2">MLMC<br>(γ = 1/150)</th><th colspan="2">F-MLMC<br>(γ = 4/150)</th></tr><tr><th>C<sub>ℓ</sub>/C<sub>L</sub></th><th>C<sub>ℓ</sub>/C<sub>ℓ-1</sub></th><th>C<sub>ℓ</sub>/C<sub>L</sub></th><th>C<sub>ℓ</sub>/C<sub>ℓ-1</sub></th></tr></thead><tbody><tr><td>L</td><td>1.00</td><td>1.99</td><td>1.00</td><td>1.95</td></tr><tr><td>L - 1</td><td>0.50</td><td>1.97</td><td>0.51</td><td>1.90</td></tr><tr><td>L - 2</td><td>0.26</td><td>1.95</td><td>0.27</td><td>1.82</td></tr><tr><td>L - 3</td><td>0.13</td><td>1.90</td><td>0.15</td><td>1.70</td></tr><tr><td>L - 4</td><td>0.07</td><td>1.82</td><td>0.09</td><td>1.53</td></tr><tr><td>L - 5</td><td>0.04</td><td>N/A</td><td>0.06</td><td>N/A</td></tr></tbody></table>

For all the experiments conducted in the remainder of this paper, the optimal multilevel sample allocation is given by eq. (16) for a fixed computational budget $C = 100C_L$. The variances $V_\ell$ required for the computation of eq. (16) are estimated in a preprocessing stage from a pilot sample of size 1000.

In what follows, $\hat{\mu}_L$ denotes an unbiased estimator of $\mu_L := \mathbb{E}[f_L(\mathbf{X}_L)]$, typically, the plain MC estimator $\hat{\mu}_L^{\text{MLMC}}$ on the finest level $L$, the MLMC estimator $\hat{\mu}_L^{\text{F-MLMC}}$ or the F-MLMC estimator $\hat{\mu}_L^{\text{F-MLMC}}$. We recall that the unbiasedness of $\hat{\mu}_L$ implies that MSE($\hat{\mu}_L$, $\mu_L$) = $V(\hat{\mu}_L)$. In order to examine the different scales of the proposed estimators, we decompose the variance $V(\hat{\mu}_L)$ into contributions of the individual Hartley modes. Specifically, exploiting the orthogonality of $\mathbf{H}_L$, it follows from eq. (29) that

$$V(\hat{\mu}_L) = \mathbb{E}[\|\hat{\mu}_L - \mu_L\|_{\mathbf{W}_L}^2] = \mathbb{E}[\|\mathbf{H}_L^T \mathbf{V}_L^T (\hat{\mu}_L - \mu_L)\|_{\mathbf{W}_L}^2], \qquad (66)$$

which, owing to the linearity of the expectation operator and the fact that $\mathbf{W}_L = n_L^{-1}\mathbf{I}_{n_L}$, may be recast as

$$V(\hat{\mu}_L) = \|\nu\|_1 := \sum_{k=0}^{n_L-1} \nu_k, \quad \nu_k := n_L^{-2} \mathbb{E}[(\mathbf{h}_k^L)^\mathrm{T} (\hat{\mu}_L - \mathbb{E}[f_L(\mathbf{X}_L)]))^2] = n_L^{-2} \mathbb{V}[(\mathbf{h}_k^L)^\mathrm{T} \hat{\mu}_L], \qquad (67)$$
where $\nu = (\nu_k)_{k=0}^{n_L-1}$ will be referred to as the spectral variance. Furthermore, we define the cumulative spectral variance $\nu^{\text{cml}} = (\nu_k^{\text{cml}})_{k=0}^{n_L-1}$ such that $\nu_k^{\text{cml}} = \sum_{k'=0}^k \nu_{k'}$, implying that the total variance is given by $\nu_{n_L-1}^{\text{cml}} = \|\nu\|_1$. For better visualization and interpretation, the columns of the 1D Hartley matrix $\mathbf{H}_L$ used in eq. (66) are actually reordered so that they are sorted by increasing representable frequency on the corresponding discrete grid (see fig. 1, which depicts the Hartley basis vectors without reordering). Specifically, the new matrix with reordered columns is obtained as $\mathbf{H}_L \mathbf{\Pi}$, where $\mathbf{\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, \dots, n_L-1$ and $k = 0, \dots, n_L/2-1$ (assuming $n_L$ is even), and where $\delta_{i,j}$ denotes the Kronecker delta.

<figure><img src="image_4.png" /><figcaption>Figure 3: Spectral variance (top) and cumulative spectral variance (bottom) of the MC estimator ($L=0$) and different MLMC (left) and F-MLMC (right) estimators ($L \in \{1, 3, 5\}$), for the illustrative 1D estimation problem described in section 4.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 $C = 100C_L$. The variance is estimated from 1000 estimators.</figcaption></figure>

Figure 3 reports the spectral variance $\nu$, as well as the corresponding cumulative variance $\nu^{\text{cml}}$, associated with the 2-, 4- and 6-level MLMC and F-MLMC estimators ($L \in \{1, 3, 5\}$) and with the single-level MC estimator ($L=0$), for the estimation of $\mathbb{E}[f_L(\mathbf{X}_L)] = \mathbf{0}_{n_L}$, where the simulators correspond to length-scale $D = 0.06$. Figures 3(a) and 3(b) show 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 figs. 3(c) and 3(d), 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. For the 2-level MLMC estimator ($L = 1$), the spectral variance plotted in fig. 3(a) 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 spectral variance, as shown in fig. 3(c). 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 ($L = 3$), this deterioration is compensated by a lower variance in the low frequencies, this is no longer the case for the 6-level MLMC estimator ($L = 5$), whose cumulative spectral variance eventually gets larger than that of the 4-level MLMC estimator, thus resulting in a larger total variance.

Figures 3(b) and 3(d) show the spectral variance $\nu$ and the cumulative variance $\nu^{cml}$ of the F-MLMC estimators with different grid hierarchies, corresponding to $L \in \{1, 3, 5\}$. The effects of the filters are especially visible on the spectral variance plotted in fig. 3(b), 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), compared to that of the unfiltered MLMC plotted in fig. 3(a). The reduced error in the high frequencies (i.e., small scales) prevents the cumulative spectral variance, plotted in fig. 3(d), from being significantly impacted in the last few Hartley modes, as opposed to that of the unfiltered MLMC, plotted in fig. 3(c). Furthermore, the reduced error in the low frequencies (i.e., large scales) translates into a lower plateau of the cumulative variance, hence a lower total variance, than for the unfiltered MLMC. This is well visible in fig. 4(a), which summarizes the total variance of the MC estimator and the different MLMC and F-MLMC estimators for $L \in \{1, \dots, 5\}$. 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, and to a 55% reduction compared to the best, 5-level unfiltered MLMC estimator.

<figure><img src="image_5.png" /><figcaption>Figure 4: Total variance of the MC estimator ($L = 0$) and different MLMC and F-MLMC estimators ($L \in \{1, \dots, 5\}$), for the illustrative 1D estimation problem described in section 4.2, with length-scale $D = 0.06$ (left) and $D = 0.01$ (right). The finest level $L$ always corresponds to a discretization with $n_L = 512$ cells, and the total budget is set to $C = 100C_L$. The variance is estimated from 1000 estimators.</figcaption></figure>

Figure 5 compares F-MLMC estimators with pre-filtering only, post-filtering only, and both pre- and post-filtering, along with the (unfiltered) MLMC and MC estimators, on the experiment with $D = 0.06$. The computational cost model for these estimators is the same as in eq. (65), with $\beta = 4$ for the F-MLMC estimator with pre-filtering only, and $\beta = 1$ for the F-MLMC estimator with post-filtering only. The total variances reported on fig. 5(a) show that both the pre-filtering and the post-filtering
operations contribute to reducing the variance of the estimator. The cumulative variance of the 6-level estimators ($L = 5$) presented in fig. 5(b) explicitly shows the effects of the pre- and post-filtering on the different frequencies of the variance. Using an estimator with post-filtering alone leads to a reduction of the variance in the high frequencies, since the spurious high-frequencies components introduced by the prolongation operator are damped. In addition, we observe that the post-filtering also improves the estimation of the low frequencies. On the other hand, the pre-filtering reduces the frequencies that cannot be represented on the coarse grids and that lead to spurious low-frequency components, thus improving the variance in those low frequencies, corresponding to the very first entries of the cumulative variance vector $v^{\text{cm1}}$. The joint use of pre- and post-filtering combines both benefits, further reducing the variance over the entire range of frequencies.

<figure><img src="image_6.png" /><figcaption>Figure 5: Total variance (left) of the MC estimator ($L = 0$) and different MLMC and F-MLMC estimators, ($L \in \{1, \dots, 5\}$), and cumulative spectral variance (right) of the different estimators with $L = 5$. Both figures are for the illustrative 1D estimation problem described in section 4.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 $C = 100C_L$. The variance is estimated from 1000 estimators.</figcaption></figure>

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 6 shows the spectral and cumulative variance of the MC, MLMC and F-MLMC estimators for a smaller 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 MLMC 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. 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, as shown in fig. 4(b), which the unfiltered MLMC estimator fails to achieve. Specifically, the variance of the 4-level F-MLMC estimator is reduced by about 60% compared to the single-level MC estimator, and by more than 50% compared to the best, 2-level unfiltered MLMC estimator.

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 is necessary to mitigate these effects, so that the multilevel estimator can reach its full potential.
<figure><img src="image_7.png" /><figcaption>(a) Spectral variance $\nu$ of MC and MLMC estimators.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>(b) Spectral variance $\nu$ of MC and F-MLMC estimators.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>(c) Cumulative spectral variance $\nu^{\text{cm1}}$ of MC and MLMC estimators.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>(d) Cumulative spectral variance $\nu^{\text{cm1}}$ of MC and F-MLMC estimators.</figcaption></figure>

<caption>Figure 6: Same as fig. 3 but with length-scale $D = 0.01$.</caption>

## 4.3 2D application

We now apply the MLMC and F-MLMC methodology to the variance estimation problem described in section 4.1, i.e., the estimation of $\theta$ in eq. (62), in a 2D setting with more complex diffusion operators than those used in the 1D illustration of the previous section. 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, 2) \times (0, 1) \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}, \qquad (68) $$

where, for $i = 1, 2$, $D_{ii}: \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. 7.

The diffusion operator in eq. (59) is discretized on Cartesian grids of size $n_\ell = n_\ell^x \times n_\ell^y$, with the length-scale fields $D_{11}$ and $D_{22}$ discretized at edge centers, and the solution discretized at cell centers, as depicted in fig. 8. The coarse operators $\mathbf{A}_\ell$ are based on restrictions of $D_{11}$ and $D_{22}$, obtained by
<figure><img src="image_11.png" /><figcaption>Figure 7: Length-scale fields $D_{11}$ (top) and $D_{22}$ (bottom) used in the reported experiments for the definition of $\mathbf{K}$ in eq. (68), with $\mu_\zeta = 0.12$ and $D_\zeta = 0.2$ (left), and with $\mu_\zeta = 0.02$ and $D_\zeta = 0.04$ (right).</figcaption></figure>

recursively averaging, at each point of a given coarse level, the values of the two nearest points of the immediately finer length-scale field, from the finest level $L$ down to the desired level $\ell < L$. 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 $x$-coordinate first, then by increasing $y$-coordinate. In other words, the entry indexed by $k = jn_\ell^x + i$ in such vectors is associated with a cell center located at $(x_i, y_j) \in \mathbb{R}^2$, where $x_i := 2(i + 1/2)/n_\ell^x$ and $y_j := (j + 1/2)/n_\ell^y$, for $i = 0, \dots, n_\ell^x - 1$ and $j = 0, \dots, n_\ell^y - 1$.

The hierarchy of simulators ($f_\ell$)$_{\ell=0}^L$ is defined through eq. (8) by $\tilde{f}_\ell: \mathbf{x}_\ell \mapsto (\mathbf{A}_\ell\mathbf{x}_\ell) \odot (\mathbf{A}_\ell\mathbf{x}_\ell)$, and we are interested in the multilevel estimation of the fine discretized field $\boldsymbol{\theta}_L := \mathbb{E}[f_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. (31), in the $x$ and $y$ directions. The 2D grid transfer operators defined in this manner satisfy eq. (20). Likewise, the second-order 2D Shapiro filter is defined as the Kronecker product of two 1D, second-order Shapiro filters defined in eq. (48). The cost model is still defined as in eq. (65), with the 2D refinement factor $s = 4$. In two dimensions, the application of $(\mathbf{I}_{n_\ell} - \boldsymbol{\Delta}_\ell)$ now requires 5 floating-point operations per cell on level $\ell$, consistent with a 5-point stencil discretization, so that $\alpha = 500$. The unfiltered 2D restriction still corresponds to $\beta = 1$, while 6 additional operations are required to apply the 2D pre-filtering operator, consistent with two sweeps (one horizontal and one vertical) of 1D filters, each represented by a 3-point stencil, hence $\beta = 7$ for F-MLMC. Table 2 summarizes the normalized costs $C_\ell/C_L$ and cost reduction factors $C_\ell/C_{\ell-1}$ of the simulators $f_\ell$ (MLMC) and $\tilde{f}_\ell$ (F-MLMC) for this 2D application. The baseline for the following experiments is a crude, single-level MC ($L = 0$) computed with a sample of size 100; hence the total budget is $C = 100C_L$.

As discussed in section 4.1, the $k$-th entry $(\boldsymbol{\theta}_L)_k$ of the exact discretized field $\boldsymbol{\theta}_L$ can be computed explicitly as $(\boldsymbol{\theta}_L)_k = (\mathbf{L}_L\mathbf{e}_k)_k$, for $k = 1, \dots, n_L$, where $\mathbf{L}_L$ denotes the discrete diffusion operator
<figure><img src="image_12.png" /><figcaption>Figure 8: Illustration of the 2D discretization of the differential operator $\nabla \cdot \mathbf{K} \nabla$ in eq. (59) on a coarse uniform grid, represented by solid lines, whose sub-tessellation into a finer, geometrically nested, uniform grid is represented by dotted lines. The discretized solution points on the fine grid (respectively, on the coarse grid) are located at the center of the fine (respectively, coarse) cells and represented by ● (respectively, ○). The fine (respectively, coarse) discretized scalar field $D_{11}$ is defined at points located on the fine (respectively, coarse) vertical edges and represented by ■ (respectively, □). The fine (respectively, coarse) discretized scalar field $D_{22}$ is defined at points located on the fine (respectively, coarse) horizontal edges and represented by ◆ (respectively, ◦).</figcaption></figure>

<table><caption>Table 2: Normalized costs $C_ℓ/C_L$ and cost reduction factors $C_ℓ/C_{ℓ-1}$ of the simulators $f_ℓ$ (MLMC) and $\bar{f}_ℓ$ (F-MLMC) corresponding to the 2D example of section 4.3.</caption><thead><tr><th rowspan="2">ℓ</th><th colspan="2">MLMC<br>(γ = 1/375)</th><th colspan="2">F-MLMC<br>(γ = 7/375)</th></tr><tr><th>C<sub>ℓ</sub>/C<sub>L</sub></th><th>C<sub>ℓ</sub>/C<sub>ℓ-1</sub></th><th>C<sub>ℓ</sub>/C<sub>L</sub></th><th>C<sub>ℓ</sub>/C<sub>ℓ-1</sub></th></tr></thead><tbody><tr><td>L</td><td>1.00</td><td>3.97</td><td>1.00</td><td>3.79</td></tr><tr><td>L - 1</td><td>0.25</td><td>3.88</td><td>0.26</td><td>3.30</td></tr><tr><td>L - 2</td><td>0.07</td><td>3.56</td><td>0.08</td><td>2.35</td></tr><tr><td>L - 3</td><td>0.02</td><td>N/A</td><td>0.03</td><td>N/A</td></tr></tbody></table>

defined in eq. (60) on the finest level $L$ and $e_k$ denotes the $k$-th canonical basis vector of $\mathbb{R}^{n_L}$ (i.e., the $k$-th column of $\mathbf{I}_{n_L}$). The top-left sub-figure of fig. 9 shows the variance field $\theta_L$ for the tensor field $\mathbf{K}$ obtained by eq. (68) from the length-scale fields $D_{11}$ and $D_{22}$ with parameters $\mu_\zeta = 0.12$ and $D_\zeta = 0.2$ depicted in fig. 7(a). 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 $C = 100C_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{\theta}] = \theta$. Consequently, the MSE solely consists of the variance of the estimators.

Figure 10 shows the variance of the considered estimators. Although the variance of the MLMC estimators (middle row) visually is 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 $\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. 10 demonstrates, at least visually, that the addition of filtering effectively damps these high-frequency components and reduced the variance. These observations are
<figure><img src="image_13.png" /><figcaption>Figure 9: The top-left sub-figure shows the exact discretized field $\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. 7(a)). The expectation is approximated from 500 estimators, each constructed with a computational budget $C = 100C_L$.</figcaption></figure>

confirmed by fig. 11, 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. (32). The 2D spectral variance of an unbiased estimator $\hat{\theta}_L$ of $\theta_L$, representing either the MC estimator on the finest level $L$, or the MLMC or F-MLMC estimator, is defined as $\nu = (\nu_k)_{k=0}^{n_L-1} \in \mathbb{R}^{n_L}$, with

$$ \forall k = 0, \dots, n_L - 1, \quad \nu_k := 4n_L^{-2}\mathbb{E}[(\mathbf{h}_k^L)^\mathrm{T}(\hat{\theta}_L - \theta_L))^2] = 4n_L^{-2}\mathbb{V}[(\mathbf{h}_k^L)^\mathrm{T}\hat{\theta}_L], \quad (69) $$

where $\mathbf{h}_k^L$ denotes the $k$-th column of $\mathbf{H}_L$. Again, the columns of $\mathbf{H}_L$ are re-ordered so that, in fig. 11, 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). Specifically, the 2D Hartley matrix with reordered columns is constructed as the Kronecker product of two 1D Hartley matrices, each with reordered columns. 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. 10, but it may also deteriorate the overall variance, as was the case in the 1D illustration (see section 4.2).
<figure><img src="image_14.png" /><figcaption>Figure 10: 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 K corresponds to $\mu_\zeta = 0.12$ and $D_\zeta = 0.2$ (fig. 7(a)). The variance is approximated from 500 estimators, each constructed with a computational budget $C = 100C_L$.</figcaption></figure>

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

$$ \nu_k^{\text{cml}} = \sum_{j=0}^{k} \sum_{i=0}^{2k+1} \nu_{jn_L^x+i}, \quad k = 0, \dots, n_L^y - 1. \qquad (70) $$

As such, the first components of $\nu^{\text{cml}}$ represent the cumulative variance associated with the low frequencies (large scales), while the last component $\nu_{n_L^y-1}^{\text{cml}}$ coincides with the total variance of the estimator. Figure 12 presents the cumulative variance of the MLMC and F-MLMC estimators for 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. 12(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. 11), but it also reduces the error on the lower frequencies, thus leading to lower total variance. In the second case (fig. 12(b)), corresponding to a tensor field K with smaller scales (see fig. 7(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
<figure><img src="image_15.png" /><figcaption>Figure 11: Spectral variance $\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. 7(a)). The variance is approximated from 500 estimators, each constructed with a computational budget $C = 100C_L$.</figcaption></figure>

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. 13.

# 5 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 prompted us to introduce 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
<figure><img src="image_16.png" /><figcaption>Figure 12: 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 = 100C_L$.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>Figure 13: 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 $C = 100C_L$.</figcaption></figure>

experiments, the MLMC estimator was still able to 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}: \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 (but actually deteriorated) 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) [9]. A similar issue exists for the estimation of covariance matrices, for which multilevel estimators are not guaranteed to be almost surely SPD [12], although advanced (but computationally expensive) approaches have been proposed to design multilevel estimators that are SPD by construction [14, 15]. 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 [4, 65] for multigrid methods. Furthermore, extensions of the MLBLUE techniques [66, 67] 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, sections 4.3–4.4], possibly in combination with the spectral pre-restriction smoothing technique of [30, 32]. Another direction for future work could be to extend the filtering ideas developed in this paper to multilevel quasi Monte Carlo integration (see, e.g., [68, 69]).

As demonstrated numerically in section 4, the introduction of pre-restriction filtering operators significantly improves the multilevel estimator in terms of variance reduction. However, in our formulation, this comes at the expense of sampling the input fields on the finest level and resorting to finest-to-coarse restriction operators, thereby increasing the sampling cost on the coarse levels (see our discussion on the cost model in section 4). It may thus be interesting to study the conditions under which the improvements in variance reduction outweigh the deterioration of the sampling cost. Furthermore, it is sometimes possible to sample coupled pairs of input fields directly on their discretization level (as opposed to our formulation where the input fields are always sampled on the finest level). Such alternative approaches would result in a more favourable cost model, since the use of restriction operators could then be avoided, but would require alternative smoothing/filtering strategies, e.g., the spectral strategies employed in [30, 32]. For our target problem of estimating the normalization coefficients of a diffusion-based covariance operator described in section 4.1, future work could focus on
alternative (F-)MLMC estimators inspired by [33], i.e., based on the joint distribution of the coupled input fields given by eq. (64), thus avoiding the use of restriction operators. In this case, spectral truncation is not readily available but could be performed, for example, after projecting the input fields onto a spectral basis. A non-trivial challenge could then lie in maintaining the telescoping sum property of the resulting multilevel estimator to ensure its unbiasedness. These ideas may be investigated in follow-up research projects.

Finally, 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 involved in the normalization problem of section 4, using the same grid hierarchy, possibly along the lines of [70], [71–74], or [75, 76].

## 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).

## A Orthogonality of the Hartley matrix

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

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

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

$$ (\mathbf{H}\mathbf{H}^\mathrm{T})_{i,j} = \sum_{k=0}^{n-1} (\cos \alpha_{ik} + \sin \alpha_{ik})(\cos \alpha_{jk} + \sin \alpha_{jk}) = \sum_{k=0}^{n-1} [\cos(\alpha_{ik} - \alpha_{jk}) + \sin(\alpha_{ik} + \alpha_{jk})] \quad (72) $$

$$ = \sum_{k=0}^{n-1} \cos \frac{2(i-j)k\pi}{n} + \sum_{k=0}^{n-1} \sin \frac{2(i+j+1)k\pi}{n}. \quad (73) $$

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

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

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

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

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

## B Prolongation of Hartley vectors

For $j, k = 0, \dots, n_0 - 1$,

$$ (\mathbf{P}\mathbf{h}_k^0)_{2j} = (\mathbf{P}\mathbf{h}_k^0)_{2j+1} = (\mathbf{h}_k^0)_j = \cos \frac{(4j+2)k\pi}{n_1} + \sin \frac{(4j+2)k\pi}{n_1}. \quad (76) $$
Upon writing $(4j+2)k\pi = (4j+1)k\pi + k\pi$ and applying elementary trigonometric identities,

$$ (\mathbf{Ph}_k^0)_{2j} = c_k(\mathbf{h}_k^1)_{2j} + \sin \frac{k\pi}{n_1} \left( \cos \frac{(4j+1)k\pi}{n_1} - \sin \frac{(4j+1)k\pi}{n_1} \right). \quad (77) $$

Noticing that

$$ \sin \frac{k\pi}{n_1} = \sin \left( \frac{(n_0 + k)\pi}{2n_0} - \frac{\pi}{2} \right) = -c_{n_0 + k}, \quad (78) $$

$$ \cos \frac{(4j+1)k\pi}{n_1} = \cos \left( \frac{(4j+1)(n_0+k)\pi}{n_1} - \frac{\pi}{2} \right) = \sin \frac{(4j+1)(n_0+k)\pi}{n_1}, \quad (79) $$

$$ \sin \frac{(4j+1)k\pi}{n_1} = \sin \left( \frac{(4j+1)(n_0+k)\pi}{n_1} - \frac{\pi}{2} \right) = -\cos \frac{(4j+1)(n_0+k)\pi}{n_1}, \quad (80) $$

we conclude that $(\mathbf{Ph}_k^0)_{2j} = c_k(\mathbf{h}_k^1)_{2j} - c_{n_0+k}(\mathbf{h}_{n_0+k}^1)_{2j}$. Then, upon writing $(4j+2)k\pi = (4j+3)k\pi - k\pi$, a similar derivation leads to $(\mathbf{Ph}_k^0)_{2j+1} = c_k(\mathbf{h}_k^1)_{2j+1} - c_{n_0+k}(\mathbf{h}_{n_0+k}^1)_{2j+1}$, thus proving eq. (33).

### C Restriction of Hartley vectors

For $j = 0, \dots, n_0 - 1$ and $k = 0, \dots, 2n_0 - 1$,

$$ (\mathbf{Rh}_k^1)_j = (\mathbf{h}_k^1)_{2j} + (\mathbf{h}_k^1)_{2j+1} \quad (81) $$

$$ = \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} \quad (82) $$

$$ = 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} \quad (83) $$

$$ = 2c_k \left[ \cos \frac{(2j+1)k\pi}{n_0} + \sin \frac{(2j+1)k\pi}{n_0} \right]. \quad (84) $$

Hence, for $j, k = 0, \dots, n_0 - 1$, it follows immediately that $\mathbf{Rh}_k^1 = 2c_k \mathbf{h}_k^0$. Furthermore,

$$ (\mathbf{Rh}_{n_0+k}^1)_j = 2c_{n_0+k} \left[ \cos\left(\frac{(2j+1)k\pi}{n_0} + \pi\right) + \sin\left(\frac{(2j+1)k\pi}{n_0} + \pi\right) \right] = -2c_{n_0+k} \mathbf{h}_k^0, \quad (85) $$

thus proving eq. (35).

### D 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. (71). To do so, we start by recalling or proving results on the node-centered Fourier basis $\tilde{\mathbf{F}} \in \mathbb{C}^{n \times n}$ and Hartley bases $\tilde{\mathbf{H}}^{\pm} \in \mathbb{R}^{n \times n}$ defined by

$$ \tilde{\mathbf{F}} := \tilde{\mathbf{H}}_c + i\tilde{\mathbf{H}}_s, \quad \tilde{\mathbf{H}}^{\pm} := \tilde{\mathbf{H}}_c \pm \tilde{\mathbf{H}}_s, \quad (\tilde{\mathbf{H}}_c)_{j,k} := \frac{1}{\sqrt{n}} \cos \frac{2jk\pi}{n}, \quad (\tilde{\mathbf{H}}_s)_{j,k} := \frac{1}{\sqrt{n}} \sin \frac{2jk\pi}{n}, \quad (86) $$

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

$$ \tilde{\mathbf{H}}_c^2 + \tilde{\mathbf{H}}_s^2 = \mathbf{I}_n, \qquad \tilde{\mathbf{H}}_c \tilde{\mathbf{H}}_s = \tilde{\mathbf{H}}_s \tilde{\mathbf{H}}_c = \mathbf{0}_n, \qquad \tilde{\mathbf{H}}_c \mathbf{1}_n = \sqrt{n} \mathbf{e}_1, \qquad \tilde{\mathbf{H}}_s \mathbf{1}_n = \mathbf{0}_n, \qquad (87) $$

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

*Proof.* The first two identities follow from [50, Lemma 1], while the last two identities follow from [77, 1.342 (1&2)]. □

**Corollary 1.** $\tilde{\mathbf{H}}^{\pm} \mathbf{1}_n = \sqrt{n} \mathbf{e}_1$.

**Corollary 2.** $\tilde{\mathbf{H}}^{\pm}$ are orthogonal, i.e., $(\tilde{\mathbf{H}}^{\pm})^2 = \mathbf{I}_n$. Furthermore, $\tilde{\mathbf{H}}^{+} \tilde{\mathbf{H}}^{-} = \tilde{\mathbf{H}}^{-} \tilde{\mathbf{H}}^{+} = \tilde{\mathbf{F}}^2$.

*Proof.* By the definitions eq. (86) and lemma 1, we have

$$ \tilde{\mathbf{H}}^{\pm} \mathbf{1}_n = \tilde{\mathbf{H}}_c \mathbf{1}_n \pm \tilde{\mathbf{H}}_s \mathbf{1}_n = \tilde{\mathbf{H}}_c \mathbf{1}_n = \sqrt{n} \mathbf{e}_1, \qquad (88) $$

$$ (\tilde{\mathbf{H}}^{\pm})^2 = (\tilde{\mathbf{H}}_c^2 + \tilde{\mathbf{H}}_s^2) \pm (\tilde{\mathbf{H}}_c \tilde{\mathbf{H}}_s + \tilde{\mathbf{H}}_s \tilde{\mathbf{H}}_c) = \tilde{\mathbf{H}}_c^2 + \tilde{\mathbf{H}}_s^2 = \mathbf{I}_n, \qquad (89) $$

$$ \tilde{\mathbf{F}}^2 = \tilde{\mathbf{H}}_c^2 - \tilde{\mathbf{H}}_s^2 + i(\tilde{\mathbf{H}}_c \tilde{\mathbf{H}}_s + \tilde{\mathbf{H}}_s \tilde{\mathbf{H}}_c) = \tilde{\mathbf{H}}_c^2 - \tilde{\mathbf{H}}_s^2, \qquad (90) $$

$$ \tilde{\mathbf{H}}^{-} \tilde{\mathbf{H}}^{+} = \tilde{\mathbf{H}}_{c}^{2} - \tilde{\mathbf{H}}_{s}^{2} + (\tilde{\mathbf{H}}_{c} \tilde{\mathbf{H}}_{s} - \tilde{\mathbf{H}}_{s} \tilde{\mathbf{H}}_{c}) = \tilde{\mathbf{F}}^{2} = \tilde{\mathbf{H}}^{+} \tilde{\mathbf{H}}^{-}. \qquad (91) $$

**Lemma 2.** Let $\mathbf{A} = \text{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, \dots, n-1$. Then $\tilde{\mathbf{H}}_c \mathbf{A} \tilde{\mathbf{H}}_s + \tilde{\mathbf{H}}_s \mathbf{A} \tilde{\mathbf{H}}_c = \tilde{\mathbf{H}}_c \mathbf{A} \tilde{\mathbf{H}}_s - \tilde{\mathbf{H}}_s \mathbf{A} \tilde{\mathbf{H}}_c = \mathbf{0}_n$.

*Proof.* Proven in the intermediary steps of the proof of [50, Theorem 1]. □

**Theorem 1.** Let $\mathbf{A} = \text{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, \dots, n-1$. Then $\tilde{\mathbf{H}}^{+} \mathbf{A} \tilde{\mathbf{H}}^{+} = \tilde{\mathbf{H}}^{-} \mathbf{A} \tilde{\mathbf{H}}^{-} = \mathbf{\Lambda}$, where $\mathbf{\Lambda} = \sqrt{n} \text{Diag}(\tilde{\mathbf{H}}^{+} \mathbf{a}) = \sqrt{n} \text{Diag}(\tilde{\mathbf{H}}^{-} \mathbf{a}) = \sqrt{n} \text{Diag}(\tilde{\mathbf{H}}_c \mathbf{a})$.

*Proof.* The fact that $\tilde{\mathbf{H}}^{+} \mathbf{A} \tilde{\mathbf{H}}^{+} = \sqrt{n} \text{Diag}(\tilde{\mathbf{H}}^{+} \mathbf{a}) = \mathbf{\Lambda} = \tilde{\mathbf{H}}_c \mathbf{A} \tilde{\mathbf{H}}_c + \tilde{\mathbf{H}}_s \mathbf{A} \tilde{\mathbf{H}}_s$ follows from [50, Theorem 1], since symmetric, circulant matrices belong to the larger class of matrices considered in [50]. Then, from the definition of $\tilde{\mathbf{H}}^{-}$ and lemma 2, we have

$$ \tilde{\mathbf{H}}^{-} \mathbf{A} \tilde{\mathbf{H}}^{-} = (\tilde{\mathbf{H}}_c \mathbf{A} \tilde{\mathbf{H}}_c + \tilde{\mathbf{H}}_s \mathbf{A} \tilde{\mathbf{H}}_s) - (\tilde{\mathbf{H}}_c \mathbf{A} \tilde{\mathbf{H}}_s + \tilde{\mathbf{H}}_s \mathbf{A} \tilde{\mathbf{H}}_c) = \tilde{\mathbf{H}}_c \mathbf{A} \tilde{\mathbf{H}}_c + \tilde{\mathbf{H}}_s \mathbf{A} \tilde{\mathbf{H}}_s = \mathbf{\Lambda}. \qquad (92) $$

Finally, from lemma 1 and corollary 1,

$$ \text{diag}(\mathbf{\Lambda}) = \mathbf{\Lambda} \mathbf{1}_n = \tilde{\mathbf{H}}^{-} \mathbf{A} \tilde{\mathbf{H}}^{-} \mathbf{1}_n = \tilde{\mathbf{H}}_c \mathbf{A} \tilde{\mathbf{H}}_c \mathbf{1}_n + \tilde{\mathbf{H}}_s \mathbf{A} \tilde{\mathbf{H}}_s \mathbf{1}_n = \sqrt{n} \tilde{\mathbf{H}}^{-} \mathbf{A} \mathbf{e}_1 = \sqrt{n} \tilde{\mathbf{H}}_c \mathbf{A} \mathbf{e}_1, \qquad (93) $$

which concludes the proof, since $\mathbf{A} \mathbf{e}_1 = \mathbf{a}$. □

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

*Proof.* From theorem 1 and corollary 2, we have

$$ \tilde{\mathbf{H}}^{+} \mathbf{A} \tilde{\mathbf{H}}^{-} = \tilde{\mathbf{H}}^{+} \mathbf{A} \tilde{\mathbf{H}}^{+} \tilde{\mathbf{H}}^{+} \tilde{\mathbf{H}}^{-} = \tilde{\mathbf{H}}^{+} \mathbf{A} \tilde{\mathbf{H}}^{+} \tilde{\mathbf{F}}^{2} = \mathbf{\Lambda} \tilde{\mathbf{F}}^{2}, \qquad (94) $$

$$ \tilde{\mathbf{H}}^{-} \mathbf{A} \tilde{\mathbf{H}}^{+} = \tilde{\mathbf{H}}^{-} \mathbf{A} \tilde{\mathbf{H}}^{-} \tilde{\mathbf{H}}^{-} \tilde{\mathbf{H}}^{+} = \tilde{\mathbf{H}}^{-} \mathbf{A} \tilde{\mathbf{H}}^{-} \tilde{\mathbf{F}}^{2} = \mathbf{\Lambda} \tilde{\mathbf{F}}^{2}, \qquad (95) $$

$$ (\tilde{\mathbf{H}}^{+} \mathbf{A} \tilde{\mathbf{H}}^{-})^{\mathrm{T}} = (\mathbf{\Lambda} \tilde{\mathbf{F}}^{2})^{\mathrm{T}} = \tilde{\mathbf{F}}^{2} \mathbf{\Lambda} = \tilde{\mathbf{H}}^{-} \mathbf{A} \tilde{\mathbf{H}}^{+} = \mathbf{\Lambda} \tilde{\mathbf{F}}^{2}, \qquad (96) $$

where the last row follows from the symmetry of $\mathbf{A}$, $\tilde{\mathbf{H}}^{\pm}$ and $\tilde{\mathbf{F}}^{2}$. □
**Lemma 3.** $H = \sqrt{n}(\check{H}^+C + \check{H}^-S)$, with $C := \text{Diag}(\{\cos \frac{k\pi}{n}\}_{k=0}^{n-1})$ and $S := \text{Diag}(\{\sin \frac{k\pi}{n}\}_{k=0}^{n-1})$.

*Proof.* From elementary trigonometric identities, we have

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

$$ \sin \frac{(2j+1)k\pi}{n} = \sin \frac{2jk\pi}{n} \cos \frac{k\pi}{n} + \cos \frac{2jk\pi}{n} \sin \frac{k\pi}{n}, \qquad (98) $$

then lemma 3 follows. □

**Lemma 4.** $C\check{F}^2S + S\check{F}^2C = 0_n$, where $C$ and $S$ are defined as in lemma 3.

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

$$ (C\check{F}^2S)_{j,k} = \cos \frac{j\pi}{n} \sin \frac{(n-j)\pi}{n} = \cos \frac{j\pi}{n} \sin \left(\pi - \frac{j\pi}{n}\right) = \cos \frac{j\pi}{n} \sin \frac{j\pi}{n}, \qquad (99) $$

$$ (C\check{F}^2S)_{k,j} = \sin \frac{j\pi}{n} \cos \frac{(n-j)\pi}{n} = \sin \frac{j\pi}{n} \cos \left(\pi - \frac{j\pi}{n}\right) = -\cos \frac{j\pi}{n} \sin \frac{j\pi}{n}, \qquad (100) $$

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

We are now ready to state and prove theorem 2.

**Theorem 2.** Let $A = \text{Circ}(a) \in \mathbb{R}^{n \times n}$ be a symmetric, circulant matrix and let $\Lambda$ be defined as in theorem 1. Then $H^T AH = n\Lambda$.

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

$$ n^{-1}H^T AH = C\check{H}^+A\check{H}^+C + S\check{H}^-A\check{H}^-S + C\check{H}^+A\check{H}^-S + S\check{H}^-A\check{H}^+C \qquad (101) $$

$$ = C\Lambda C + S\Lambda S + C\Lambda\check{F}^2S + S\Lambda\check{F}^2C = \Lambda(C^2 + S^2 + C\check{F}^2S + S\check{F}^2C), \qquad (102) $$

so that $n^{-1}H^T AH = \Lambda(C^2 + S^2) = \Lambda$ since $C^2 + S^2 = I_n$ by elementary trigonometry. □

# E Optimality of the Galerkin operator

Let $n_1, n_0 \in \mathbb{N}$ be such that $n_1 > n_0 > 0$, and let $A_1 \in \mathbb{R}^{n_1 \times n_1}$ and $A_0 \in \mathbb{R}^{n_0 \times n_0}$ be two linear operators. Let $R \in \mathbb{R}^{n_0 \times n_1}$ and $P \in \mathbb{R}^{n_1 \times n_0}$ be a restriction operator and a prolongation operator, respectively, such that $R = P^T$ and $P^T W_1 P = W_0$ (see eq. (20)). We seek the operator $A_0^*$ that minimizes $\|A_1 - PA_0R\|_{F,W_1}^2$ among the linear operators $A_0 \in \mathbb{R}^{n_0 \times n_0}$. We restrict ourselves to the case where, for $\ell \in \{0, 1\}$, $W_\ell = n_\ell^{-1}I_{n_\ell}$, so that $P^T P = RR^T = (n_1/n_0)I_{n_0}$ and

$$ A_0^* := \arg\min_{A_0 \in \mathbb{R}^{n_0 \times n_0}} \|A_1 - PA_0R\|_{F,W_1}^2 = \arg\min_{A_0 \in \mathbb{R}^{n_0 \times n_0}} \|A_1 - PA_0R\|_{F}^2, \qquad (103) $$

where $\| \cdot \|_F$ denotes the classical (unweighted) Frobenius norm. We start by writing

$$ \|A_1 - PA_0R\|_F^2 = \|A_1 - PA_0R\|_F^2 = \|\text{vec}(A_1 - PA_0R)\|_2^2 \qquad (104) $$

$$ = \|\text{vec}(A_1) - \text{vec}(PA_0R)\|_2^2 = \|\text{vec}(A_1) - (R^T \otimes P) \text{vec}(A_0)\|_2^2, \qquad (105) $$
so that eq. (103) is recast as an ordinary linear least squares problem whose solution is

$$ \mathrm{vec}(\mathbf{A}_0^*) = ((\mathbf{R}^\mathrm{T} \otimes \mathbf{P})^\mathrm{T} (\mathbf{R}^\mathrm{T} \otimes \mathbf{P}))^{-1} (\mathbf{R}^\mathrm{T} \otimes \mathbf{P})^\mathrm{T} \mathrm{vec}(\mathbf{A}_1) \quad (106) $$

$$ = ((\mathbf{R} \otimes \mathbf{P}^\mathrm{T})(\mathbf{R}^\mathrm{T} \otimes \mathbf{P}))^{-1} (\mathbf{R} \otimes \mathbf{P}^\mathrm{T}) \mathrm{vec}(\mathbf{A}_1) \quad (107) $$

$$ = (\mathbf{R}\mathbf{R}^\mathrm{T} \otimes \mathbf{P}^\mathrm{T}\mathbf{P})^{-1} \mathrm{vec}(\mathbf{P}^\mathrm{T}\mathbf{A}_1\mathbf{R}^\mathrm{T}) = ((n_1/n_0)\mathbf{I}_{n_0} \otimes (n_1/n_0)\mathbf{I}_{n_0})^{-1} \mathrm{vec}(\mathbf{R}\mathbf{A}_1\mathbf{P}) \quad (108) $$

$$ = (n_0/n_1)^2 (\mathbf{I}_{n_0} \otimes \mathbf{I}_{n_0}) \mathrm{vec}(\mathbf{R}\mathbf{A}_1\mathbf{P}) = (n_0/n_1)^2 \mathrm{vec}(\mathbf{R}\mathbf{A}_1\mathbf{P}), \quad (109) $$

thus proving that $\mathbf{A}_0^* = (n_0/n_1)^2\mathbf{R}\mathbf{A}_1\mathbf{P}$.

## F Filtered restriction of Hartley vectors

We start by looking at the effect of the Shapiro filter on the fine Hartley basis vectors,

$$ \forall j, k = 0, \dots, n_1 - 1, \quad (\mathbf{S}_1 \mathbf{h}_k^1)_j = \frac{1}{4} [(\mathbf{h}_k^1)_{j-1} + 2(\mathbf{h}_k^1)_j + (\mathbf{h}_k^1)_{j+1}] . \quad (110) $$

Elementary trigonometric identities imply that, for $j, k = 0, \dots, n_1 - 1$,

$$ \cos \frac{(2j - 1)k\pi}{n_1} + \cos \frac{(2j + 1)k\pi}{n_1} = 2c_k \cos \frac{2jk\pi}{n_1}, \quad (111) $$

$$ \cos \frac{(2j + 1)k\pi}{n_1} + \cos \frac{(2j + 3)k\pi}{n_1} = 2c_k \cos \frac{(2j + 2)k\pi}{n_1}, \quad (112) $$

$$ \sin \frac{(2j - 1)k\pi}{n_1} + \sin \frac{(2j + 1)k\pi}{n_1} = 2c_k \sin \frac{2jk\pi}{n_1}, \quad (113) $$

$$ \sin \frac{(2j + 1)k\pi}{n_1} + \sin \frac{(2j + 3)k\pi}{n_1} = 2c_k \sin \frac{(2j + 2)k\pi}{n_1}, \quad (114) $$

leading to

$$ (\mathbf{S}_1 \mathbf{h}_k^1)_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}_k^1)_j, \quad (115) $$

where the last identity is obtained by applying the same trigonometric identities as for eqs. (111) to (114), and thus showing that $\mathbf{S}_1\mathbf{h}_k^1 = c_k^2\mathbf{h}_k^1$, for $k = 0, \dots, n_1 - 1$. Then, eq. (33) implies that, for $k = 0, \dots, n_0 - 1$, $\bar{\mathbf{P}}\mathbf{h}_k^0 = \mathbf{S}_1\mathbf{P}\mathbf{h}_k^0 = c_k\mathbf{S}_1\mathbf{h}_k^1 - c_{n_0+k}\mathbf{S}_1\mathbf{h}_{n_0+k}^1$, from which eq. (49) follows. Furthermore, eq. (50) follows immediately from the direct application of eq. (35) to $\bar{\mathbf{R}}\mathbf{h}_k^1 = \mathbf{R}\mathbf{S}_1\mathbf{h}_k^1 = c_k^2\mathbf{R}\mathbf{h}_k^1$.

## G Factorization of the diffusion-based covariance matrix

Let $\mathbf{\Delta}$ and $\mathbf{W}$ as defined in section 4. Because $\mathbf{\Delta}$ is self-adjoint with respect to $\langle \cdot, \cdot \rangle_{\mathbf{W}}$, we have

$$ \mathbf{W}(\mathbf{I} - \mathbf{\Delta}) = \mathbf{W} - \mathbf{W}\mathbf{\Delta} = \mathbf{W} - \mathbf{\Delta}^\mathrm{T}\mathbf{W} = (\mathbf{I} - \mathbf{\Delta})^\mathrm{T}\mathbf{W}, \quad (116) $$

and thus, for $q \in \mathbb{N}$,

$$ \mathbf{W}(\mathbf{I} - \mathbf{\Delta})^q = \mathbf{W}(\mathbf{I} - \mathbf{\Delta})(\mathbf{I} - \mathbf{\Delta})^{q-1} = (\mathbf{I} - \mathbf{\Delta})^\mathrm{T}\mathbf{W}(\mathbf{I} - \mathbf{\Delta})^{q-1} = \dots = [(\mathbf{I} - \mathbf{\Delta})^\mathrm{T}]^q \mathbf{W}. \quad (117) $$

It follows that

$$ \mathbf{A} := (\mathbf{I} - \mathbf{\Delta})^{-q} \mathbf{W}^{-1} = (\mathbf{W}(\mathbf{I} - \mathbf{\Delta})^q)^{-1} = (([\mathbf{I} - \mathbf{\Delta}]^q)^\mathrm{T} \mathbf{W})^{-1} = \mathbf{W}^{-1} [(\mathbf{I} - \mathbf{\Delta})^{-q}]^\mathrm{T}, \quad (118) $$

and thus

$$ \mathbf{L} := (\mathbf{I} - \mathbf{\Delta})^{-2q} \mathbf{W}^{-1} = (\mathbf{I} - \mathbf{\Delta})^{-q} \mathbf{A} = (\mathbf{I} - \mathbf{\Delta})^{-q} \mathbf{W}^{-1} [(\mathbf{I} - \mathbf{\Delta})^{-q}]^\mathrm{T} = \mathbf{A}\mathbf{W}\mathbf{A}^\mathrm{T}. \quad (119) $$
References

[1] M. B. Giles. "Multilevel Monte Carlo Methods." In: *Acta Numerica* 24 (2015), pp. 259-328. DOI: 10.1017/S096249291500001X.

[2] B. Peherstorfer, K. Willcox, and M. Gunzburger. "Survey of Multifidelity Methods in Uncertainty Propagation, Inference, and Optimization." In: *SIAM Review* 60.3 (2018), pp. 550-591. ISSN: 0036-1445. DOI: 10.1137/16M1082469.

[3] A. Brandt. "Guide to Multigrid Development." In: *Multigrid Methods*. Ed. by W. Hackbusch and U. Trottenberg. Lecture Notes in Mathematics. Berlin, Heidelberg: Springer, 1982, pp. 220-312. ISBN: 978-3-540-39544-7. DOI: 10.1007/BFb0069930.

[4] W. Hackbusch. *Multi-Grid Methods and Applications*. Vol. 4. Springer Series in Computational Mathematics. Berlin, Heidelberg: Springer, 1985. ISBN: 978-3-642-05722-9 978-3-662-02427-0. DOI: 10.1007/978-3-662-02427-0.

[5] U. Trottenberg, C. W. Oosterlee, and A. Schuller. *Multigrid*. Elsevier, 2000. ISBN: 978-0-08-047956-9. URL: https://www.elsevier.com/books/multigrid/trottenberg/978-0-08-047956-9.

[6] W. L. Briggs, V. E. Henson, and S. F. McCormick. *A Multigrid Tutorial*. 2nd. Other Titles in Applied Mathematics. Society for Industrial and Applied Mathematics, 2000. ISBN: 978-0-89871-462-3. DOI: 10.1137/1.9780898719505.

[7] S. Heinrich. "Multilevel Monte Carlo Methods." In: *Large-Scale Scientific Computing*. Ed. by S. Margenov, J. Waśniewski, and P. Yalamov. Lecture Notes in Computer Science. Berlin, Heidelberg: Springer, 2001, pp. 58-67. ISBN: 978-3-540-45346-8. DOI: 10.1007/3-540-45346-6_5.

[8] M. B. Giles. "Multilevel Monte Carlo Path Simulation." In: *Operations Research* 56.3 (2008), pp. 607-617. ISSN: 0030-364X, 1526-5463. DOI: 10.1287/opre.1070.0496.

[9] C. Bierig and A. Chernov. "Convergence Analysis of Multilevel Monte Carlo Variance Estimators and Application for Random Obstacle Problems." In: *Numerische Mathematik* 130.4 (2015), pp. 579-613. ISSN: 0029-599X, 0945-3245. DOI: 10.1007/s00211-014-0676-3.

[10] C. Bierig and A. 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. ISSN: 2194-0401, 2194-041X. DOI: 10.1007/s40072-015-0063-9.

[11] P. Mycek and M. De Lozzo. "Multilevel Monte Carlo Covariance Estimation for the Computation of Sobol' Indices." In: *SIAM/ASA Journal on Uncertainty Quantification* 7.4 (2019), pp. 1323-1348. ISSN: 2166-2525. DOI: 10.1137/18M1216389.

[12] H. Hoel, K. J. H. Law, and R. Tempone. "Multilevel Ensemble Kalman Filtering." In: *SIAM Journal on Numerical Analysis* 54.3 (2016), pp. 1813-1839. ISSN: 0036-1429. DOI: 10.1137/15M100955X.

[13] M. Destouches, P. Mycek, and S. Gürol. *Multivariate Extensions of the Multilevel Best Linear Unbiased Estimator for Ensemble-Variational Data Assimilation*. Technical Report TR-PA-23-67. Cerfacs, 2023. DOI: 10.48550/arXiv.2306.07017. arXiv: 2306.07017.

[14] A. Maurais, T. Alsup, B. Peherstorfer, and Y. Marzouk. "Multi-Fidelity Covariance Estimation in the Log-Euclidean Geometry." In: *Proceedings of the 40th International Conference on Machine Learning*. PMLR, 2023, pp. 24214-24235. URL: https://proceedings.mlr.press/v202/maurais23a.html.
[15] A. Maurais, T. Alsup, B. Peherstorfer, and Y. M. Marzouk. "Multifidelity Covariance Estimation via Regression on the Manifold of Symmetric Positive Definite Matrices." In: *SIAM Journal on Mathematics of Data Science* (2025), pp. 189-223. DOI: 10.1137/23M159247X.

[16] M. Destouches, P. Mycek, S. Gürol, A. T. Weaver, S. Gratton, and E. Simon. "Multilevel Monte Carlo Methods for Ensemble Variational Data Assimilation." In: *EGUsphere* (2024), pp. 1-33. DOI: 10.5194/egusphere-2024-3628.

[17] L. M. Stein. *Interpolation of Spatial Data*. Springer-Verlag, New York, NY, 1999.

[18] A. Weaver and P. Courtier. "Correlation Modelling on the Sphere Using a Generalized Diffusion Equation." In: *Quarterly Journal of the Royal Meteorological Society* 127.575 (2001), pp. 1815-1846. ISSN: 1477-870X. DOI: 10.1002/qj.49712757518.

[19] R. J. Purser, W.-S. Wu, D. F. Parrish, and N. M. Roberts. "Numerical Aspects of the Application of Recursive Filters to Variational Statistical Analysis. Part I: Spatially Homogeneous and Isotropic Gaussian Covariances." In: *Monthly Weather Review* 131.8 (2003), pp. 1524-1535. ISSN: 1520-0493, 0027-0644. DOI: 10.1175/1520-0493(2003)131<1524:NAOTAO>2.0.CO;2.

[20] F. Lindgren, H. Rue, and J. 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. ISSN: 1467-9868. DOI: 10.1111/j.1467-9868.2011.00777.x.

[21] T. Bui-Thanh, O. Ghattas, J. Martin, and G. 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), A2494-A2523. ISSN: 1064-8275, 1095-7197. DOI: 10.1137/12089586X.

[22] 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 (2017), S873-S897. ISSN: 1064-8275. DOI: 10.1137/16M1083591.

[23] A. T. Weaver, M. Chrust, B. Ménétrier, and A. Piacentini. "An Evaluation of Methods for Normalizing Diffusion-based Covariance Operators in Variational Data Assimilation." In: *Quarterly Journal of the Royal Meteorological Society* 147.734 (2021), pp. 289-320. ISSN: 0035-9009, 1477-870X. DOI: 10.1002/qj.3918.

[24] K. A. Cliffe, M. B. Giles, R. Scheichl, and A. L. 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. ISSN: 1432-9360, 1433-0369. DOI: 10.1007/s00791-011-0160-x.

[25] A. Barth, C. Schwab, and N. Zollinger. "Multi-Level Monte Carlo Finite Element Method for Elliptic PDEs with Stochastic Coefficients." In: *Numerische Mathematik* 119.1 (2011), pp. 123-161. ISSN: 0029-599X, 0945-3245. DOI: 10.1007/s00211-011-0377-0.

[26] S. Mishra and C. Schwab. "Sparse Tensor Multi-Level Monte Carlo Finite Volume Methods for Hyperbolic Conservation Laws with Random Initial Data." In: *Mathematics of Computation* 81.280 (2012), pp. 1979-2018. ISSN: 0025-5718, 1088-6842. DOI: 10.1090/S0025-5718-2012-02574-9.

[27] S. Mishra, C. Schwab, and J. Šukys. "Multilevel Monte Carlo Finite Volume Methods for Shallow Water Equations with Uncertain Topography in Multi-dimensions." In: *SIAM Journal on Scientific Computing* 34.6 (2012), B761-B784. ISSN: 1064-8275, 1095-7197. DOI: 10.1137/110857295.

[28] A. Abdulle, A. Barth, and C. Schwab. "Multilevel Monte Carlo Methods for Stochastic Elliptic Multiscale PDEs." In: *Multiscale Modeling & Simulation* 11.4 (2013), pp. 1033-1070. ISSN: 1540-3459, 1540-3467. DOI: 10.1137/120894725.
[29] C. J. Gittelson, J. Könnö, C. Schwab, and R. Stenberg. "The Multi-Level Monte Carlo Finite Element Method for a Stochastic Brinkman Problem." In: *Numerische Mathematik* 125.2 (2013), pp. 347–386. ISSN: 0945-3245. DOI: 10.1007/s00211-013-0537-5.

[30] A. L. Teckentrup, R. Scheichl, M. B. Giles, and E. Ullmann. "Further Analysis of Multilevel Monte Carlo Methods for Elliptic PDEs with Random Coefficients." In: *Numerische Mathematik* 125.3 (2013), pp. 569–600. ISSN: 0945-3245. DOI: 10.1007/s00211-013-0546-4.

[31] J. Charrier, R. Scheichl, and A. L. Teckentrup. "Finite Element Error Analysis of Elliptic PDEs with Random Coefficients and Its Application to Multilevel Monte Carlo Methods." In: *SIAM Journal on Numerical Analysis* 51.1 (2013), pp. 322–352. ISSN: 0036-1429, 1095-7170. DOI: 10.1137/110853054.

[32] A. Istratuca and A. L. Teckentrup. "Smoothed Circulant Embedding with Applications to Multi-level Monte Carlo Methods for PDEs with Random Coefficients." In: *IMA Journal of Numerical Analysis* (2025), drae102. ISSN: 0272-4979. DOI: 10.1093/imanum/drae102.

[33] M. Croci, M. B. Giles, M. E. Rognes, and P. E. Farrell. "Efficient White Noise Sampling and Coupling for Multilevel Monte Carlo with Nonnested Meshes." In: *SIAM/ASA Journal on Uncertainty Quantification* 6.4 (2018), pp. 1630–1655. DOI: 10.1137/18M1175239.

[34] R. Wienands and W. Joppich. *Practical Fourier Analysis for Multigrid Methods*. New York: Chapman and Hall/CRC, 2004. ISBN: 978-0-429-14328-1. DOI: 10.1201/9781420034998.

[35] G. Da Prato and J. Zabczyk. *Stochastic Equations in Infinite Dimensions*. 2nd ed. Encyclopedia of Mathematics and Its Applications. Cambridge: Cambridge University Press, 2014. ISBN: 978-1-107-05584-1. DOI: 10.1017/CBO9781107295513.

[36] T. Hytönen, J. Van Neerven, M. Veraar, and L. Weis. *Analysis in Banach Spaces*. Cham: Springer International Publishing, 2016. ISBN: 978-3-319-48519-5 978-3-319-48520-1. DOI: 10.1007/978-3-319-48520-1.

[37] G. Madec, M. Bell, A. Blaker, C. Bricaud, D. Bruciaferri, M. Castrillo, D. Calvert, J. Chanut, E. Clementi, A. Coward, I. Epicoco, C. Éthé, J. Ganderton, J. Harle, K. Hutchinson, D. Iovino, D. Lea, T. Lovato, M. Martin, N. Martin, F. Mele, D. Martins, S. Masson, P. Mathiot, F. Mele, S. Mocavero, S. Müller, A. J. G. Nurser, S. Paronuzzi, M. Peltier, R. Person, C. Rousset, S. Rynders, G. Samson, S. Téchené, M. Vancoppenolle, and C. Wilson. *NEMO Ocean Engine Reference Manual*. 27. Reference Manual. Zenodo, 2023. DOI: 10.5281/zenodo.8167700.

[38] P. Wesseling. "Cell-Centered Multigrid for Interface Problems." In: *Journal of Computational Physics* 79.1 (1988), pp. 85–91. ISSN: 0021-9991. DOI: 10.1016/0021-9991(88)90005-8.

[39] P. Hemker. "On the Order of Prolongations and Restrictions in Multigrid Procedures." In: *Journal of Computational and Applied Mathematics* 32.3 (1990), pp. 423–429. ISSN: 03770427. DOI: 10.1016/0377-0427(90)90047-4.

[40] M. Khalil and P. Wesseling. "Vertex-Centered and Cell-Centered Multigrid for Interface Problems." In: *Journal of Computational Physics* 98.1 (1992), pp. 1–10. ISSN: 0021-9991. DOI: 10.1016/0021-9991(92)90168-X.

[41] M. Mohr and R. Wienands. "Cell-Centred Multigrid Revisited." In: *Computing and Visualization in Science* 7.3 (2004), pp. 129–140. ISSN: 1433-0369. DOI: 10.1007/s00791-004-0137-0.

[42] M. Kac and A. J. F. Siegert. "An Explicit Representation of a Stationary Gaussian Process." In: *The Annals of Mathematical Statistics* 18.3 (1947), pp. 438–442. ISSN: 0003-4851, 2168-8990. DOI: 10.1214/aoms/1177730391.
[43] O. P. Le Maître and O. M. Knio. *Spectral Methods for Uncertainty Quantification*. Scientific Computation. Dordrecht: Springer Netherlands, 2010. ISBN: 978-90-481-3519-6. DOI: 10.1007/978-90-481-3520-2.

[44] W. Betz, I. Papaioannou, and D. Straub. "Numerical Methods for the Discretization of Random Fields by Means of the Karhunen-Loève Expansion." In: *Computer Methods in Applied Mechanics and Engineering* 271 (2014), pp. 109-129. ISSN: 0045-7825. DOI: 10.1016/j.cma.2013.12.010.

[45] I. Graham, F. Kuo, D. Nuyens, R. Scheichl, and I. Sloan. "Quasi-Monte Carlo Methods for Elliptic PDEs with Random Coefficients and Applications." In: *Journal of Computational Physics* 230.10 (2011), pp. 3668-3694. ISSN: 00219991. DOI: 10.1016/j.jcp.2011.01.023.

[46] I. G. Graham, F. Y. Kuo, D. Nuyens, R. Scheichl, and I. H. Sloan. "Analysis of Circulant Embedding Methods for Sampling Stationary Random Fields." In: *SIAM Journal on Numerical Analysis* 56.3 (2018), pp. 1871-1895. ISSN: 0036-1429, 1095-7170. DOI: 10.1137/17M1149730.

[47] M. R. El Amri, P. Mycek, S. Ricci, and M. De Lozzo. *Multilevel Surrogate-based Control Variates*. Tech. rep. arXiv:2306.10800. arXiv, 2024. DOI: 10.48550/arXiv.2306.10800. arXiv: 2306.10800.

[48] R. Hartley. "A More Symmetrical Fourier Analysis Applied to Transmission Problems." In: *Proceedings of the IRE* 30.3 (1942), pp. 144-150. ISSN: 2162-6634. DOI: 10.1109/JRPROC.1942.234333.

[49] R. N. Bracewell. "Discrete Hartley Transform." In: *JOSA* 73.12 (1983), pp. 1832-1835. DOI: 10.1364/JOSA.73.001832.

[50] D. Bini and P. Favati. "On a Matrix Algebra Related to the Discrete Hartley Transform." In: *SIAM Journal on Matrix Analysis and Applications* 14.2 (1993), pp. 500-507. ISSN: 0895-4798. DOI: 10.1137/0614035.

[51] R. Shapiro. "Smoothing, Filtering, and Boundary Effects." In: *Reviews of Geophysics* 8.2 (1970), p. 359. ISSN: 8755-1209. DOI: 10.1029/RG008i002p00359.

[52] 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. ISSN: 0021-9991. DOI: 10.1016/j.jcp.2013.07.001.

[53] F. Falissard. "Uneven-Order Decentered Shapiro Filters for Boundary Filtering." In: *Journal of Computational Physics* 292 (2015), pp. 168-175. ISSN: 0021-9991. DOI: 10.1016/j.jcp.2015.03.003.

[54] P. Whittle. "Stochastic processes in several dimensions." In: *Bull. Inst. Internat. Statist.* 40 (1963), pp. 974-994.

[55] P. Guttorp and T. Gneiting. "Studies in the History of Probability and Statistics XLIX: On the Matérn Correlation Family." In: *Biometrika* 93.4 (2006), pp. 989-995. ISSN: 0006-3444. DOI: 10.1093/biomet/93.4.989.

[56] A. T. Weaver and I. Mirouze. "On the Diffusion Equation and Its Application to Isotropic and Anisotropic Correlation Modelling in Variational Assimilation." In: *Quarterly Journal of the Royal Meteorological Society* 139.670 (2013), pp. 242-260. ISSN: 00359009. DOI: 10.1002/qj.1955.

[57] 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. ISSN: 01689274. DOI: 10.1016/j.apnum.2007.01.003.

[58] E. Hallman, I. C. F. Ipsen, and A. K. Saibaba. "Monte Carlo Methods for Estimating the Diagonal of a Real Symmetric Matrix." In: *SIAM Journal on Matrix Analysis and Applications* 44.1 (2023), pp. 240-269. ISSN: 0895-4798. DOI: 10.1137/22M1476277.
[59] A. T. Weaver, J. Tshimanga, and A. Piacentini. "Correlation Operators Based on an Implicitly Formulated Diffusion Equation Solved with the Chebyshev Iteration." In: *Quarterly Journal of the Royal Meteorological Society* 142.694 (2016), pp. 455-471. ISSN: 1477-870X. DOI: 10.1002/qj.2664.

[60] M. Chrust, A. T. Weaver, P. Browne, H. Zuo, and M. A. Balmaseda. "Impact of Ensemble-Based Hybrid Background-Error Covariances in ECMWF's next-Generation Ocean Reanalysis System." In: *Quarterly Journal of the Royal Meteorological Society* 151.767 (2025), e4914. ISSN: 1477-870X. DOI: 10.1002/qj.4914.

[61] O. Goux, S. Gürol, A. T. Weaver, Y. Diouane, and O. Guillet. "Impact of Correlated Observation Errors on the Conditioning of Variational Data Assimilation Problems." In: *Numerical Linear Algebra with Applications* 31.1 (2024), e2529. ISSN: 1099-1506. DOI: 10.1002/nla.2529.

[62] G. H. Golub and R. S. Varga. "Chebyshev Semi-Iterative Methods, Successive Overrelaxation Iterative Methods, and Second Order Richardson Iterative Methods: Part I." In: *Numerische Mathematik* 3.1 (1961), pp. 147-156. ISSN: 0945-3245. DOI: 10.1007/BF01386013.

[63] G. H. Golub and R. S. Varga. "Chebyshev Semi-Iterative Methods, Successive Overrelaxation Iterative Methods, and Second Order Richardson Iterative Methods: Part II." In: *Numerische Mathematik* 3.1 (1961), pp. 157-168. ISSN: 0945-3245. DOI: 10.1007/BF01386014.

[64] G. H. Golub and C. F. Van Loan. *Matrix Computations*. Fourth edition. Johns Hopkins Studies in the Mathematical Sciences. Baltimore: The Johns Hopkins University Press, 2013. ISBN: 978-1-4214-0794-4. URL: https://epubs.siam.org/doi/book/10.1137/1.9781421407944.

[65] L. Debreu, E. Neveu, E. Simon, F.-X. Le Dimet, and A. Vidard. "Multigrid Solvers and Multi-grid Preconditioners for the Solution of Variational Data Assimilation Problems." In: *Quarterly Journal of the Royal Meteorological Society* 142.694 (2016), pp. 515-528. ISSN: 1477-870X. DOI: 10.1002/qj.2676.

[66] D. Schaden and E. Ullmann. "On Multilevel Best Linear Unbiased Estimators." In: *SIAM/ASA Journal on Uncertainty Quantification* 8.2 (2020), pp. 601-635. DOI: 10.1137/19M1263534.

[67] D. Schaden and E. Ullmann. "Asymptotic Analysis of Multilevel Best Linear Unbiased Estimators." In: *SIAM/ASA Journal on Uncertainty Quantification* 9.3 (2021), pp. 953-978. DOI: 10.1137/20M1321607.

[68] L. Herrmann and C. Schwab. "Multilevel Quasi-Monte Carlo Integration with Product Weights for Elliptic PDEs with Lognormal Coefficients." In: *ESAIM: Mathematical Modelling and Numerical Analysis* 53.5 (2019), pp. 1507-1552. ISSN: 0764-583X, 1290-3841. DOI: 10.1051/m2an/2019016.

[69] M. Croci, M. Giles, and P. E. Farrell. "Multilevel Quasi Monte Carlo Methods for Elliptic PDEs with Random Field Coefficients via Fast White Noise Sampling." In: *SIAM Journal on Scientific Computing* 43.4 (2021), A2840-A2868. ISSN: 1064-8275. DOI: 10.1137/20M1329044.

[70] S. Osborn, P. S. Vassilevski, and U. Villa. "A Multilevel, Hierarchical Sampling Technique for Spatially Correlated Random Fields." In: *SIAM Journal on Scientific Computing* 39.5 (2017), S543-S562. ISSN: 1064-8275. DOI: 10.1137/16M1082688.

[71] P. Kumar, C. W. Oosterlee, and R. 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 (2017). ISSN: 2152-5080, 2152-5099. DOI: 10.1615/Int.J.UncertaintyQuantification.2016018677.

[72] P. Kumar, C. Rodrigo, F. J. Gaspar, and C. W. Oosterlee. "On Local Fourier Analysis of Multigrid Methods for PDEs with Jumping and Random Coefficients." In: *SIAM Journal on Scientific Computing* 41.3 (2019), A1385-A1413. ISSN: 1064-8275. DOI: 10.1137/18M1173769.
[73] P. Robbe, D. Nuyens, and S. Vandewalle. "Recycling Samples in the Multigrid Multilevel (Quasi-)Monte Carlo Method." In: *SIAM Journal on Scientific Computing* 41.5 (2019), S37-S60. ISSN: 1064-8275. DOI: 10.1137/18M1194031.

[74] P. Robbe, D. Nuyens, and S. Vandewalle. "Enhanced Multi-Index Monte Carlo by Means of Multiple Semicoarsened Multigrid for Anisotropic Diffusion Problems." In: *Numerical Linear Algebra with Applications* 28.3 (2021), e2281. ISSN: 1099-1506. DOI: 10.1002/nla.2281.

[75] J. Goodman and A. D. Sokal. "Multigrid Monte Carlo Method. Conceptual Foundations." In: *Physical Review D* 40.6 (1989), pp. 2035-2071. DOI: 10.1103/PhysRevD.40.2035.

[76] Y. Kazashi, E. H. Müller, and R. Scheichl. *Multigrid Monte Carlo Revisited: Theory and Bayesian Inference*. 2024. DOI: 10.48550/arXiv.2407.12149. arXiv: 2407.12149.

[77] I. S. Gradshteyn and I. M. Ryshik. *Table of Integrals, Series, and Products*. 8th. Elsevier, 2015. ISBN: 978-0-12-384933-5. URL: https://linkinghub.elsevier.com/retrieve/pii/C20100648395.

[78] P. J. Davis. *Circulant Matrices: Second Edition*. 2nd. American Mathematical Society, 2012. ISBN: 978-0-8218-9165-0. URL: https://bookstore.ams.org/chel-338.