<figure><img src="image_1.png" /></figure>

HAL
open science

# Multilevel Surrogate-based Control Variates

Mohamed Reda El Amri, Paul Mycek, Sophie Ricci, Matthias De Lozzo

► To cite this version:

Mohamed Reda El Amri, Paul Mycek, Sophie Ricci, Matthias De Lozzo. Multilevel Surrogate-based Control Variates. 2024. hal-04129561v2

HAL Id: hal-04129561

https://hal.science/hal-04129561v2

Preprint submitted on 18 Jun 2024

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.
# Multilevel Surrogate-based Control Variates

Mohamed Reda El Amri, Paul Mycek, Sophie Ricci, Matthias De Lozzo

► To cite this version:

Mohamed Reda El Amri, Paul Mycek, Sophie Ricci, Matthias De Lozzo. Multilevel Surrogate-based Control Variates. IRT Saint-Exupéry; Cerfacs. 2024. hal-04129561v2

HAL Id: hal-04129561

https://hal.science/hal-04129561v2

Submitted on 18 Jun 2024

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>

Distributed under a Creative Commons Attribution 4.0 International License
# Multilevel Surrogate-based Control Variates*

Reda El Amri†

Paul Mycek‡§

Sophie Ricci‡§

Matthias De Lozzo§

## Technical report

### Abstract

Estimating statistical parameters of the output of a numerical simulator with random inputs by crude Monte Carlo sampling can be inaccurate when simulations are computationally expensive, so variance reduction techniques must be considered, such as the method of control variates (CVs) or the multilevel Monte Carlo (MLMC) method. In this paper, we propose three multifidelity variance reduction strategies relying on surrogate-based CVs, possibly combined with MLMC. MLCV extends the CV method by using multiple auxiliary random variables devised from surrogate models for simulators of different levels. MLMC-CV improves the MLMC estimator by using a CVs based on a surrogate of the correction term at each level. Further variance reduction is achieved by using the surrogate-based CVs of all the levels in the MLMC-MLCV strategy. Alternative solutions that reduce the subset of surrogates used for the multilevel estimation are also introduced. These strategies are validated on an uncertain 1D heat equation test case from the MLMC literature, where the statistic is the mean integrated temperature at a given time. The results are assessed in terms of the accuracy and computational cost of the multilevel estimators. It is shown that when the lower fidelity outputs are strongly correlated with the high-fidelity outputs, a significant variance reduction is obtained when using surrogate models for the coarser levels only. Taking advantage of pre-existing surrogate models proves to be an even more efficient strategy.

**Key words:** Multifidelity, multilevel Monte Carlo, control variates, surrogate models, polynomial chaos, variance reduction.

## 1 Introduction

In recent years, the propagation of uncertainties in numerical simulators has become an essential step in the study of physical phenomena. Therefore, uncertainty quantification (UQ) has emerged as an important element in scientific computing [50, 52]. For complex nonlinear systems, the task of quantifying the effect of uncertainties on the simulator behaviour can pose major challenges as closed-form solutions often do not exist. Sampling-based algorithms are considered the default approach when it comes to UQ for complex nonlinear simulators. Monte Carlo (MC) sampling is arguably the most popular and flexible method. It consists of extracting statistical information from a sample of simulator responses. Due to its non-intrusive nature, MC is straightforward to implement and provides unbiased estimators for classical statistical parameters, e.g., the expectation and the variance. Furthermore, the convergence of MC estimators in terms of their root-mean-square error

<footnote>*This work was funded by the PIA framework (CGI, ANR) and the industrial members of the IRT Saint Exupéry project R-Evol: Airbus, Liebherr, Altran Technologies, Capgemini DEMS France, CENAERO and Cerfacs.</footnote>

<footnote>†IFP Énergies Nouvelles, Solaize, France (mohamed-reda.el-amri@ifpen.fr).</footnote>

<footnote>‡CECI, CNRS - Cerfacs, Toulouse, France (mycek@cerfacs.fr, ricci@cerfacs.fr).</footnote>

<footnote>§IRT Saint Exupéry, Toulouse, France (matthias.delozzo@irt-saintexupery.com).</footnote>
# Multilevel Surrogate-based Control Variates*

Reda El Amri†

Paul Mycek‡§

Sophie Ricci‡§

Matthias De Lozzo§

## Technical report

### Abstract

Estimating statistical parameters of the output of a numerical simulator with random inputs by crude Monte Carlo sampling can be inaccurate when simulations are computationally expensive, so variance reduction techniques must be considered, such as the method of control variates (CVs) or the multilevel Monte Carlo (MLMC) method. In this paper, we propose three multifidelity variance reduction strategies relying on surrogate-based CVs, possibly combined with MLMC. MLCV extends the CV method by using multiple auxiliary random variables devised from surrogate models for simulators of different levels. MLMC-CV improves the MLMC estimator by using a CVs based on a surrogate of the correction term at each level. Further variance reduction is achieved by using the surrogate-based CVs of all the levels in the MLMC-MLCV strategy. Alternative solutions that reduce the subset of surrogates used for the multilevel estimation are also introduced. These strategies are validated on an uncertain 1D heat equation test case from the MLMC literature, where the statistic is the mean integrated temperature at a given time. The results are assessed in terms of the accuracy and computational cost of the multilevel estimators. It is shown that when the lower fidelity outputs are strongly correlated with the high-fidelity outputs, a significant variance reduction is obtained when using surrogate models for the coarser levels only. Taking advantage of pre-existing surrogate models proves to be an even more efficient strategy.

**Key words:** Multifidelity, multilevel Monte Carlo, control variates, surrogate models, polynomial chaos, variance reduction.

## 1 Introduction

In recent years, the propagation of uncertainties in numerical simulators has become an essential step in the study of physical phenomena. Therefore, uncertainty quantification (UQ) has emerged as an important element in scientific computing [50, 52]. For complex nonlinear systems, the task of quantifying the effect of uncertainties on the simulator behaviour can pose major challenges as closed-form solutions often do not exist. Sampling-based algorithms are considered the default approach when it comes to UQ for complex nonlinear simulators. Monte Carlo (MC) sampling is arguably the most popular and flexible method. It consists of extracting statistical information from a sample of simulator responses. Due to its non-intrusive nature, MC is straightforward to implement and provides unbiased estimators for classical statistical parameters, e.g., the expectation and the variance. Furthermore, the convergence of MC estimators in terms of their root-mean-square error

<footnote>*This work was funded by the PIA framework (CGI, ANR) and the industrial members of the IRT Saint Exupéry project R-Evol: Airbus, Liebherr, Altran Technologies, Capgemini DEMS France, CENAERO and Cerfacs.</footnote>

<footnote>†IFP Énergies Nouvelles, Solaize, France (mohamed-reda.el-amri@ifpen.fr).</footnote>

<footnote>‡CECI, CNRS - Cerfacs, Toulouse, France (mycek@cerfacs.fr, ricci@cerfacs.fr).</footnote>

<footnote>§IRT Saint Exupéry, Toulouse, France (matthias.delozzo@irt-saintexupery.com).</footnote>
(RMSE) as a function of the sample size *n* is independent of the stochastic input dimension, but, on the other hand, it is very slow, namely at a rate of *n*⁻⁰.⁵. For computationally expensive simulators, the cost of MC is often considered impractically high. In some cases, slight improvements can be obtained through the use of importance sampling [27], latin hypercube sampling [35] or quasi-Monte Carlo (QMC) [13].

Another well-known approach in UQ consists in replacing the simulator by a surrogate model. Fast-to-evaluate surrogates can be used to approximate the high-fidelity model response and therefore reduce drastically the computational cost of estimating statistics. Polynomial chaos (PC) [34, 37], Gaussian process models [45], radial basis functions [41] and (deep) neural networks [54] are commonly used surrogate models. The downside of the surrogate-based approach is that it introduces approximation error, which causes biases in the statistics estimators. Besides, this approximation error tends to increase in high dimension.

One recent MC sampling framework based on control variates (CV) [29, 30, 31, 33] has been extensively developed and used. In a sampling-based CV strategy, one seeks to reduce the variance of the MC estimator of a random variable, arising from the high-fidelity model, by exploiting its correlation with an auxiliary random variable that arises from low-fidelity models approximating the same input-output relationship. In classic CV theory, the mean of the auxiliary random variable is assumed to be known. Unfortunately, in many cases such an assumption is not valid. This creates the need to use another estimator for the auxiliary random variable [42, 39, 43, 44, 24, 48, 49, 47], which involves an additional computational cost. Recently, the adoption of surrogate models as such auxiliary random variables, in particular PC models, has been explored in [16, 56, 17, 25]. The benefits are twofold: (1) the prediction of the surrogate models may be highly correlated with the output high-fidelity model, leading to a reduced variance of the CV estimator as compared to the MC estimation; (2) certain surrogate models provide exact statistics that are needed by the CV approach and for the others, these statistics can be approximated at negligible cost with arbitrary precision.

Another variance reduction technique, called the multilevel MC (MLMC) method [21, 10, 22], uses a hierarchy of models. Originally devised for the estimation of expected values, MLMC has since been extended to the estimation of other statistics, see, e.g., [4, 5] for the estimation of variance and higher order central moments and [36] for the computation of Sobol' indices. Typically, multilevel methods are based on a sequence of levels which correspond to a hierarchy of simulators with increasing accuracy and cost. From a practical standpoint, the different levels often correspond to simulators with increasing mesh resolutions. This translates into a lower accuracy of the so-called *coarse* levels, whereas the *finer* levels correspond to accurate simulators. By construction, MLMC results in an unbiased estimator. It relies on a telescopic sum of terms based on the differences between the successive simulators. It debiases the MC estimator associated with the lowest-level simulator. Many other unbiased multilevel estimators have been devised in recent years. The multi-index MC estimator [26] is an extension of the MLMC estimator such that the telescoping sum idea is used in multiple directions. In [23], the authors developed a variant of MLMC which uses QMC samples instead of independent MC samples for each level.

In this work, we propose new estimation techniques to quantify efficiently the output uncertainties of simulators with limited computational budgets. These approaches build on existing statistical techniques to combine the best of them: MC sampling to guarantee the unbiasedness of the estimators and the independence of their accuracy with respect to the input dimension, multilevel estimation to reduce their variance by exploiting some hierarchical fidelity structure, and surrogate-based CVs to further reduce the variance, especially on the coarser levels, as efficiently as possible. It is worth noting that we do not seek to estimate statistics from accurate surrogate models, but rather from sampling techniques assisted by multifidelity simulators and coarse-to-fine
(RMSE) as a function of the sample size *n* is independent of the stochastic input dimension, but, on the other hand, it is very slow, namely at a rate of *n*⁻⁰.⁵. For computationally expensive simulators, the cost of MC is often considered impractically high. In some cases, slight improvements can be obtained through the use of importance sampling [27], latin hypercube sampling [35] or quasi-Monte Carlo (QMC) [13].

Another well-known approach in UQ consists in replacing the simulator by a surrogate model. Fast-to-evaluate surrogates can be used to approximate the high-fidelity model response and therefore reduce drastically the computational cost of estimating statistics. Polynomial chaos (PC) [34, 37], Gaussian process models [45], radial basis functions [41] and (deep) neural networks [54] are commonly used surrogate models. The downside of the surrogate-based approach is that it introduces approximation error, which causes biases in the statistics estimators. Besides, this approximation error tends to increase in high dimension.

One recent MC sampling framework based on control variates (CV) [29, 30, 31, 33] has been extensively developed and used. In a sampling-based CV strategy, one seeks to reduce the variance of the MC estimator of a random variable, arising from the high-fidelity model, by exploiting its correlation with an auxiliary random variable that arises from low-fidelity models approximating the same input-output relationship. In classic CV theory, the mean of the auxiliary random variable is assumed to be known. Unfortunately, in many cases such an assumption is not valid. This creates the need to use another estimator for the auxiliary random variable [42, 39, 43, 44, 24, 48, 49, 47], which involves an additional computational cost. Recently, the adoption of surrogate models as such auxiliary random variables, in particular PC models, has been explored in [16, 56, 17, 25]. The benefits are twofold: (1) the prediction of the surrogate models may be highly correlated with the output high-fidelity model, leading to a reduced variance of the CV estimator as compared to the MC estimation; (2) certain surrogate models provide exact statistics that are needed by the CV approach and for the others, these statistics can be approximated at negligible cost with arbitrary precision.

Another variance reduction technique, called the multilevel MC (MLMC) method [21, 10, 22], uses a hierarchy of models. Originally devised for the estimation of expected values, MLMC has since been extended to the estimation of other statistics, see, e.g., [4, 5] for the estimation of variance and higher order central moments and [36] for the computation of Sobol' indices. Typically, multilevel methods are based on a sequence of levels which correspond to a hierarchy of simulators with increasing accuracy and cost. From a practical standpoint, the different levels often correspond to simulators with increasing mesh resolutions. This translates into a lower accuracy of the so-called *coarse* levels, whereas the *finer* levels correspond to accurate simulators. By construction, MLMC results in an unbiased estimator. It relies on a telescopic sum of terms based on the differences between the successive simulators. It debiases the MC estimator associated with the lowest-level simulator. Many other unbiased multilevel estimators have been devised in recent years. The multi-index MC estimator [26] is an extension of the MLMC estimator such that the telescoping sum idea is used in multiple directions. In [23], the authors developed a variant of MLMC which uses QMC samples instead of independent MC samples for each level.

In this work, we propose new estimation techniques to quantify efficiently the output uncertainties of simulators with limited computational budgets. These approaches build on existing statistical techniques to combine the best of them: MC sampling to guarantee the unbiasedness of the estimators and the independence of their accuracy with respect to the input dimension, multilevel estimation to reduce their variance by exploiting some hierarchical fidelity structure, and surrogate-based CVs to further reduce the variance, especially on the coarser levels, as efficiently as possible. It is worth noting that we do not seek to estimate statistics from accurate surrogate models, but rather from sampling techniques assisted by multifidelity simulators and coarse-to-fine
surrogate models. A first strategy, named *multilevel control variates* (MLCV), combines multiple CVs based on surrogate models of simulators with different fidelity levels. Although this strategy does not build on the MLMC approach specifically, it still exploits multilevel information through the surrogates constructed at different levels. The second strategy, named *multilevel Monte Carlo with control variates* (MLMC-MLCV), utilizes CVs in an information fusion framework to exploit synergies between the flexible MLMC sampling and the correlation shared between the high- and low-fidelity components. The numerical results show that the unbiased MLMC-MLCV estimators can converge faster than the existing estimators.

This technical report is organized as follows. We introduce notations and the necessary mathematical background in section 2. In section 3, we briefly discuss the MLMC estimator and then define our proposed MLMC-MLCV estimator that combines multilevel sampling and surrogate-based CV. In section 4, we conduct numerical experiments to support the theoretical results. Section 5 proposes concluding remarks.

# 2 Background

In this section, we summarize the important results of statistical estimation using control variates and we show how surrogate models can be leveraged in this setting.

## 2.1 Notation

We first introduce a few notations. Throughout the rest of this paper, the high-fidelity numerical model is abstractly represented by the deterministic mapping

$$f: \Xi \to \mathbb{R} \\ \mathbf{x} \mapsto f(\mathbf{x}), \qquad (2.1)$$

where $\mathbf{x} := (x_1 \cdots x_d)^\top$ is a vector of $d$ uncertain input parameters evolving in a measurable space denoted by $\Xi = \Xi_1 \times \cdots \times \Xi_d$, with $\Xi_i \subset \mathbb{R}$. Following a probabilistic approach, the $d$-dimensional input is modelled as a random vector $\mathbf{X}: \Omega \to \mathbb{R}^d$ with known joint probability distribution. In this work, we seek an accurate estimator of some statistic $\theta$ of the output random variable $Y := f(\mathbf{X}): \Omega \to \mathbb{R}$, e.g., its expected value ($\theta = \mathbb{E}[Y]$) or variance ($\theta = \mathbb{V}[Y]$), at a reasonable computational cost $\mathcal{C}$.

## 2.2 Crude Monte Carlo

In practice, the MC estimator $\hat{\theta}$ of the statistic $\theta$ is based on $n$ observations $Y^{(1)}, \dots, Y^{(n)}$ defined as $Y^{(i)} = f(\mathbf{X}^{(i)})$ where $\mathcal{X} = \{\mathbf{X}^{(1)}, \dots, \mathbf{X}^{(n)}\}$ is a $n$-sample of $\mathbf{X}$, that is, a collection of independent and identically distributed random variables with the same distribution as $\mathbf{X}$. For instance, the *sample mean* and (unbiased) *sample variance* estimators,

$$\hat{E}[Y] = n^{-1} \sum_{i=1}^{n} Y^{(i)} \quad \text{and} \quad \hat{V}[Y] = (n - 1)^{-1} \sum_{i=1}^{n} (Y^{(i)} - \hat{E}[Y])^2 \qquad (2.2)$$

are MC estimators of the expectation and variance of $Y$, respectively. It is well known that these estimators are *unbiased*, that is, $\mathbb{E}[\hat{\theta}] = \theta$, so that the mean square error (MSE) of $\hat{\theta}$ reduces to the variance of the estimator:

$$\text{MSE}(\hat{\theta}, \theta) := \mathbb{E}[(\hat{\theta} - \theta)^2] = \mathbb{V}[\hat{\theta}] + (\mathbb{E}[\hat{\theta}] - \theta)^2 = \mathbb{V}[\hat{\theta}]. \qquad (2.3)$$
surrogate models. A first strategy, named *multilevel control variates* (MLCV), combines multiple CVs based on surrogate models of simulators with different fidelity levels. Although this strategy does not build on the MLMC approach specifically, it still exploits multilevel information through the surrogates constructed at different levels. The second strategy, named *multilevel Monte Carlo with control variates* (MLMC-MLCV), utilizes CVs in an information fusion framework to exploit synergies between the flexible MLMC sampling and the correlation shared between the high- and low-fidelity components. The numerical results show that the unbiased MLMC-MLCV estimators can converge faster than the existing estimators.

This technical report is organized as follows. We introduce notations and the necessary mathematical background in section 2. In section 3, we briefly discuss the MLMC estimator and then define our proposed MLMC-MLCV estimator that combines multilevel sampling and surrogate-based CV. In section 4, we conduct numerical experiments to support the theoretical results. Section 5 proposes concluding remarks.

# 2 Background

In this section, we summarize the important results of statistical estimation using control variates and we show how surrogate models can be leveraged in this setting.

## 2.1 Notation

We first introduce a few notations. Throughout the rest of this paper, the high-fidelity numerical model is abstractly represented by the deterministic mapping

$$f: \Xi \to \mathbb{R} \\ \mathbf{x} \mapsto f(\mathbf{x}), \qquad (2.1)$$

where $\mathbf{x} := (x_1 \cdots x_d)^T$ is a vector of $d$ uncertain input parameters evolving in a measurable space denoted by $\Xi = \Xi_1 \times \cdots \times \Xi_d$, with $\Xi_i \subset \mathbb{R}$. Following a probabilistic approach, the $d$-dimensional input is modelled as a random vector $\mathbf{X}: \Omega \to \mathbb{R}^d$ with known joint probability distribution. In this work, we seek an accurate estimator of some statistic $\theta$ of the output random variable $Y := f(\mathbf{X}): \Omega \to \mathbb{R}$, e.g., its expected value ($\theta = \mathbb{E}[Y]$) or variance ($\theta = \mathbb{V}[Y]$), at a reasonable computational cost $\mathcal{C}$.

## 2.2 Crude Monte Carlo

In practice, the MC estimator $\hat{\theta}$ of the statistic $\theta$ is based on $n$ observations $Y^{(1)}, \dots, Y^{(n)}$ defined as $Y^{(i)} = f(\mathbf{X}^{(i)})$ where $\mathcal{X} = \{\mathbf{X}^{(1)}, \dots, \mathbf{X}^{(n)}\}$ is a $n$-sample of $\mathbf{X}$, that is, a collection of independent and identically distributed random variables with the same distribution as $\mathbf{X}$. For instance, the *sample mean* and (unbiased) *sample variance* estimators,

$$\hat{E}[Y] = n^{-1} \sum_{i=1}^{n} Y^{(i)} \quad \text{and} \quad \hat{V}[Y] = (n-1)^{-1} \sum_{i=1}^{n} (Y^{(i)} - \hat{E}[Y])^2 \qquad (2.2)$$

are MC estimators of the expectation and variance of $Y$, respectively. It is well known that these estimators are *unbiased*, that is, $\mathbb{E}[\hat{\theta}] = \theta$, so that the mean square error (MSE) of $\hat{\theta}$ reduces to the variance of the estimator:

$$\text{MSE}(\hat{\theta}, \theta) := \mathbb{E}[(\hat{\theta} - \theta)^2] = \mathbb{V}[\hat{\theta}] + (\mathbb{E}[\hat{\theta}] - \theta)^2 = \mathbb{V}[\hat{\theta}]. \qquad (2.3)$$
The convergence of these MC estimators is known to be slow, so that variance reduction techniques are needed when dealing with computationally expensive simulators.

## 2.3 Control variates

In this section, we present a well-known variance reduction technique using auxiliary random variables $Z_1, \dots, Z_M$ as control variates. We denote by $\tau_1, \dots, \tau_M$ statistical parameters related to the control variates $Z_1, \dots, Z_M$, and assume they are known exactly. Then, the CV estimator is defined as

$$ \hat{\theta}^{\text{CV}}(\alpha) = \hat{\theta} - \alpha^{\top}(\hat{\tau} - \tau), \quad (2.4) $$

where $\hat{\theta}$ and $\hat{\tau} = (\hat{\tau}_1 \dots \hat{\tau}_M)^{\top}$ are unbiased MC estimators of $\theta$ and $\tau = (\tau_1 \dots \tau_M)^{\top}$, respectively, based on a common input $n$-sample $\mathcal{X}$, and where $\alpha \in \mathbb{R}^M$ is the control parameter. We note that the CV estimator is unbiased by construction, regardless of the value of the parameter $\alpha$, and that its variance reads

$$ \mathbb{V}[\hat{\theta}^{\text{CV}}(\alpha)] = \mathbb{V}[\hat{\theta}] + \alpha^{\top}\Sigma\alpha - 2\alpha^{\top}\mathbf{c}. \quad (2.5) $$

with $\mathbf{c} := \mathbb{C}[\hat{\tau}, \hat{\theta}] \in \mathbb{R}^M$ and $\Sigma := \mathbb{C}[\hat{\tau}] \in \mathbb{R}^{M \times M}$. To fully take advantage of the control variates, the parameter $\alpha$ is selected so as to minimize the variance (2.5) of the CV estimator. Assuming that the covariance matrix $\Sigma$ is non-singular (and thus symmetric positive definite, SPD), the first- and second-order optimality conditions of the minimization problem imply that there exists a unique optimal solution,

$$ \alpha^* = \Sigma^{-1}\mathbf{c}. \quad (2.6) $$

The optimal CV estimator is thus $\hat{\theta}^{\text{CV}}(\alpha^*)$, and its variance is given by

$$ \mathbb{V}[\hat{\theta}^{\text{CV}}(\alpha^*)] = (1 - R^2)\mathbb{V}[\hat{\theta}], \quad (2.7) $$

where $R^2 = \mathbb{V}[\hat{\theta}]^{-1}\mathbf{c}^{\top}\Sigma^{-1}\mathbf{c}$ is nonnegative, as $\Sigma^{-1}$ is SPD. Denoting by $\mathbf{D} = \text{Diag}(\Sigma) \in \mathbb{R}^{M \times M}$ the diagonal matrix consisting of the diagonal of $\Sigma$, it can be shown further that $R^2 = \mathbf{r}^{\top}\mathbf{R}^{-1}\mathbf{r}$, where $\mathbf{r} = (\mathbb{V}[\hat{\theta}]\mathbf{D})^{-1/2}\mathbf{c}$ is the vector of the Pearson correlation coefficients between $\hat{\theta}$ and $\hat{\tau}$, and $\mathbf{R} = \mathbf{D}^{-1/2}\Sigma\mathbf{D}^{-1/2}$ is the correlation matrix of $\hat{\tau}$. Thus, $R^2 \in [0, 1]$ corresponds to the squared coefficient of multiple correlation between $\hat{\theta}$ and the elements of $\hat{\tau}$ [1, section 2.5.2]. Consequently, $0 \le \mathbb{V}[\hat{\theta}^{\text{CV}}(\alpha^*)] \le \mathbb{V}[\hat{\theta}]$, and we refer to $R^2$ as the variance reduction factor of the CV estimator. This shows that the variance of the optimal CV estimator $\hat{\theta}^{\text{CV}}(\alpha^*)$ is always reduced (or, rigorously speaking, not increased) compared to the MC estimator $\hat{\theta}$. Furthermore, the higher $R^2$, the greater the reduction in variance.

**Remark 2.1.** The requirement that $\Sigma$ be non-singular is a reasonable one. Indeed, let us suppose that $\Sigma$ is singular. Then, because $\Sigma$ is positive semi-definite by construction, this implies that there exists a nonzero vector $\eta \in \mathbb{R}^M \setminus \{0\}$ such that $\eta^{\top}\Sigma\eta = \mathbb{V}[\eta^{\top}\hat{\tau}] = 0$, indicating that any one element of $\hat{\tau}$ can be expressed as an affine function of the others. As such, it does not bring any additional information to the CV estimator, so that at least one of the control variates can simply be discarded.

A desirable property of the CV estimator is that increasing the number of control variates improves the CV estimator. Specifically, under mild assumptions, proposition 2.1 states that the variance of the CV estimator is reduced (or, rigorously speaking, not increased) when adding a new control variate.

**Proposition 2.1.** Let $\hat{\tau}_+ := [\hat{\tau}^{\top} \ \hat{\tau}_{M+1}]^{\top}$, $\tau_+ := [\tau^{\top} \ \tau_{M+1}]^{\top}$, and define $\mathbf{c}_+ = \mathbb{C}[\hat{\tau}_+, \hat{\theta}]$, and $\Sigma_+ = \mathbb{C}[\hat{\tau}_+]$. We further assume that $\Sigma$ and $\Sigma_+$ are non-singular and that $\mathbb{V}[\hat{\tau}_{M+1}] > 0$. Let $\hat{\theta}^{\text{CV}}(\alpha_+^*) := \hat{\theta} - \alpha_+^{\top}(\hat{\tau}_+ - \tau_+)$, with $\alpha_+^* = \Sigma_+^{-1}\mathbf{c}_+$, be the optimal CV estimator based on $M+1$ control variates, and let $R_+^2$ denote its variance reduction factor. Then $R_+^2 \ge R^2$.
The convergence of these MC estimators is known to be slow, so that variance reduction techniques are needed when dealing with computationally expensive simulators.

## 2.3 Control variates

In this section, we present a well-known variance reduction technique using auxiliary random variables $Z_1, \dots, Z_M$ as control variates. We denote by $\tau_1, \dots, \tau_M$ statistical parameters related to the control variates $Z_1, \dots, Z_M$, and assume they are known exactly. Then, the CV estimator is defined as

$$ \hat{\theta}^{\text{CV}}(\alpha) = \hat{\theta} - \alpha^{\top}(\hat{\tau} - \tau), \qquad (2.4) $$

where $\hat{\theta}$ and $\hat{\tau} = (\hat{\tau}_1 \dots \hat{\tau}_M)^{\top}$ are unbiased MC estimators of $\theta$ and $\tau = (\tau_1 \dots \tau_M)^{\top}$, respectively, based on a common input $n$-sample $\mathcal{X}$, and where $\alpha \in \mathbb{R}^M$ is the control parameter. We note that the CV estimator is unbiased by construction, regardless of the value of the parameter $\alpha$, and that its variance reads

$$ \mathbb{V}[\hat{\theta}^{\text{CV}}(\alpha)] = \mathbb{V}[\hat{\theta}] + \alpha^{\top}\Sigma\alpha - 2\alpha^{\top}\mathbf{c}. \qquad (2.5) $$

with $\mathbf{c} := \mathbb{C}[\hat{\tau}, \hat{\theta}] \in \mathbb{R}^M$ and $\Sigma := \mathbb{C}[\hat{\tau}] \in \mathbb{R}^{M \times M}$. To fully take advantage of the control variates, the parameter $\alpha$ is selected so as to minimize the variance (2.5) of the CV estimator. Assuming that the covariance matrix $\Sigma$ is non-singular (and thus symmetric positive definite, SPD), the first- and second-order optimality conditions of the minimization problem imply that there exists a unique optimal solution,

$$ \alpha^* = \Sigma^{-1}\mathbf{c}. \qquad (2.6) $$

The optimal CV estimator is thus $\hat{\theta}^{\text{CV}}(\alpha^*)$, and its variance is given by

$$ \mathbb{V}[\hat{\theta}^{\text{CV}}(\alpha^*)] = (1 - R^2)\mathbb{V}[\hat{\theta}], \qquad (2.7) $$

where $R^2 = \mathbb{V}[\hat{\theta}]^{-1}\mathbf{c}^{\top}\Sigma^{-1}\mathbf{c}$ is nonnegative, as $\Sigma^{-1}$ is SPD. Denoting by $\mathbf{D} = \text{Diag}(\Sigma) \in \mathbb{R}^{M \times M}$ the diagonal matrix consisting of the diagonal of $\Sigma$, it can be shown further that $R^2 = \mathbf{r}^{\top}\mathbf{R}^{-1}\mathbf{r}$, where $\mathbf{r} = (\mathbb{V}[\hat{\theta}]\mathbf{D})^{-1/2}\mathbf{c}$ is the vector of the Pearson correlation coefficients between $\hat{\theta}$ and $\hat{\tau}$, and $\mathbf{R} = \mathbf{D}^{-1/2}\Sigma\mathbf{D}^{-1/2}$ is the correlation matrix of $\hat{\tau}$. Thus, $R^2 \in [0, 1]$ corresponds to the squared coefficient of multiple correlation between $\hat{\theta}$ and the elements of $\hat{\tau}$ [1, section 2.5.2]. Consequently, $0 \le \mathbb{V}[\hat{\theta}^{\text{CV}}(\alpha^*)] \le \mathbb{V}[\hat{\theta}]$, and we refer to $R^2$ as the variance reduction factor of the CV estimator. This shows that the variance of the optimal CV estimator $\hat{\theta}^{\text{CV}}(\alpha^*)$ is always reduced (or, rigorously speaking, not increased) compared to the MC estimator $\hat{\theta}$. Furthermore, the higher $R^2$, the greater the reduction in variance.

**Remark 2.1.** The requirement that $\Sigma$ be non-singular is a reasonable one. Indeed, let us suppose that $\Sigma$ is singular. Then, because $\Sigma$ is positive semi-definite by construction, this implies that there exists a nonzero vector $\eta \in \mathbb{R}^M \setminus \{0\}$ such that $\eta^{\top}\Sigma\eta = \mathbb{V}[\eta^{\top}\hat{\tau}] = 0$, indicating that any one element of $\hat{\tau}$ can be expressed as an affine function of the others. As such, it does not bring any additional information to the CV estimator, so that at least one of the control variates can simply be discarded.

A desirable property of the CV estimator is that increasing the number of control variates improves the CV estimator. Specifically, under mild assumptions, proposition 2.1 states that the variance of the CV estimator is reduced (or, rigorously speaking, not increased) when adding a new control variate.

**Proposition 2.1.** Let $\hat{\tau}_+ := [\hat{\tau}^{\top} \ \hat{\tau}_{M+1}]^{\top}$, $\tau_+ := [\tau^{\top} \ \tau_{M+1}]^{\top}$, and define $\mathbf{c}_+ = \mathbb{C}[\hat{\tau}_+, \hat{\theta}]$, and $\Sigma_+ = \mathbb{C}[\hat{\tau}_+]$. We further assume that $\Sigma$ and $\Sigma_+$ are non-singular and that $\mathbb{V}[\hat{\tau}_{M+1}] > 0$. Let $\hat{\theta}^{\text{CV}}(\alpha_+^*) := \hat{\theta} - \alpha_+^{\top}(\hat{\tau}_+ - \tau_+)$, with $\alpha_+^* = \Sigma_+^{-1}\mathbf{c}_+$, be the optimal CV estimator based on $M+1$ control variates, and let $R_+^2$ denote its variance reduction factor. Then $R_+^2 \ge R^2$.
*Proof.* We have $R_+^2 = \mathbf{r}_+^\top \mathbf{R}_+^{-1} \mathbf{r}_+$, where $\mathbf{r}_+ = [\mathbf{r}^\top \ \gamma]^\top$ and $\mathbf{R}_+ = \begin{bmatrix} \mathbf{R} & \mathbf{u} \\ \mathbf{u}^\top & 1 \end{bmatrix}$, and where $\mathbf{u} = (\mathbb{V}[\hat{\tau}_{M+1}]\mathbf{D})^{-1/2}\mathbb{C}[\hat{\tau}, \hat{\tau}_{M+1}]$ and $\gamma = (\mathbb{V}[\hat{\theta}]\mathbb{V}[\hat{\tau}_{M+1}])^{-1/2}\mathbb{C}[\hat{\theta}, \hat{\tau}_{M+1}]$. Note that $\mathbf{R}$ and $\mathbf{R}_+$ are both non-singular, because so are $\mathbf{\Sigma}$ and $\mathbf{\Sigma}_+$. The augmented matrix $\mathbf{R}_+$ may be inverted by block using the Schur complement $s = 1 - \mathbf{u}^\top \mathbf{R}^{-1} \mathbf{u} \neq 0$ of the (1,1)-block $\mathbf{R}$. It follows that $R_+^2 = R^2 + s^{-1}(\gamma - \mathbf{u}^\top \mathbf{R}^{-1} \mathbf{r})^2$. Now, because $\mathbf{R}_+$ is SPD (it is positive semidefinite by construction and non-singular by assumption), $\mathbf{x}^\top \mathbf{R}_+ \mathbf{x} > 0$ for any choice of $\mathbf{x} \neq \mathbf{0}$. The particular choice of $\mathbf{x} = [-\mathbf{u}^\top \mathbf{R}^{-1} \ \mathbf{1}]^\top$ implies that $s > 0$, which in turn implies that $R_+^2 \geq R^2$. $\square$

**Corollary 2.2.** Under the assumptions of proposition 2.1, the equality $R_+^2 = R^2$ holds if and only if $\gamma = \mathbf{u}^\top \mathbf{R}^{-1} \mathbf{r}$.

*Proof.* Straightforward from $R_+^2 = R^2 + s^{-1}(\gamma - \mathbf{u}^\top \mathbf{R}^{-1} \mathbf{r})^2$ with $s > 0$. $\square$

For the CV estimation of specific statistics, the expressions of $\mathbf{c}$ and $\mathbf{\Sigma}$ may be further reduced to involve statistics on $Y$ and $\mathbf{Z} := (Z_1 \cdots Z_M)^\top$ directly. These expressions, as well as the consequences on eqs. (2.6) and (2.7), are given in appendix A for the CV estimators of $\mathbb{E}[Y]$ and $\mathbb{V}[Y]$.

It should be noted at this point that, in practice, the optimal parameter $\alpha^*$ needs to be approximated. Specifically, the statistics involved in $\mathbf{c}$ and $\mathbf{\Sigma}$ need to be estimated. This can be done either using the same sample $\mathcal{X}$ as for the CV estimator itself, or using an independent, pilot sample $\mathcal{X}'$. The former strategy introduces a bias in the resulting CV estimator, while the latter guarantees unbiasedness but requires additional high-fidelity simulator evaluations, so that the former strategy is generally preferred. Furthermore, statistical remedies such as jackknifing, splitting or bootstrapping have been proposed to reduce or eliminate the bias introduced by the former strategy [38]. In both strategies, however, neither the theoretical variance reduction given by eq. (2.7) nor proposition 2.1 hold anymore.

**Remark 2.2.** In some instances, because it only involves the control variates, $\mathbf{\Sigma}$ may be known exactly (see remark 3.1 for an example in a multilevel setting involving PC-based control variates) or estimated accurately using many independent samples at negligible cost. Unfortunately, this is not the case for $\mathbf{c}$, which involves the output $Y = f(\mathbf{X})$ of the high-fidelity simulator $f$.

## 2.4 Surrogate-based control variates

The efficiency of the CV approach relies on the strong correlation between $Y = f(\mathbf{X})$ and the control variates $\mathbf{Z}$. In a multifidelity framework, the control variates correspond to the output of low-fidelity versions of $f$, typically in the form of simulators with simplified physics and/or coarser discretization. In many applications, the exact statistical measures $\tau$ of such control variates $\mathbf{Z}$ may not be available. One way to circumvent this limitation is to use additional samples to also estimate $\tau$. This type of estimators led to an approximate CV class of methods [24, 39, 43], where different model management and sample allocation strategies may be used to find a suitable trade-off between the additional cost of sampling $\mathbf{Z}$ and the resulting variance reduction. Recently, an optimal strategy was proposed [48, 49, 47] leading to the so-called multilevel best linear unbiased estimator (MLBLUE), as well as a unifying framework for a large class of multilevel and multifidelity estimators. However, ACV estimators typically yield lower variance reduction than their exact CV counterparts, as illustrated in appendix B on bi-fidelity (A)CV estimators of the expectation.

In this paper, we focus on the case where the low-fidelity models correspond to surrogate models of the high-fidelity simulator $f$, which allows us to use the original CV approach described in
*Proof.* We have $R_+^2 = \mathbf{r}_+^\top \mathbf{R}_+^{-1} \mathbf{r}_+$, where $\mathbf{r}_+ = [\mathbf{r}^\top \ \gamma]^\top$ and $\mathbf{R}_+ = \begin{bmatrix} \mathbf{R} & \mathbf{u} \\ \mathbf{u}^\top & 1 \end{bmatrix}$, and where $\mathbf{u} = (\mathbb{V}[\hat{\tau}_{M+1}]\mathbf{D})^{-1/2}\mathbb{C}[\hat{\tau}, \hat{\tau}_{M+1}]$ and $\gamma = (\mathbb{V}[\hat{\theta}]\mathbb{V}[\hat{\tau}_{M+1}])^{-1/2}\mathbb{C}[\hat{\theta}, \hat{\tau}_{M+1}]$. Note that $\mathbf{R}$ and $\mathbf{R}_+$ are both non-singular, because so are $\mathbf{\Sigma}$ and $\mathbf{\Sigma}_+$. The augmented matrix $\mathbf{R}_+$ may be inverted by block using the Schur complement $s = 1 - \mathbf{u}^\top \mathbf{R}^{-1} \mathbf{u} \neq 0$ of the (1,1)-block $\mathbf{R}$. It follows that $R_+^2 = R^2 + s^{-1}(\gamma - \mathbf{u}^\top \mathbf{R}^{-1} \mathbf{r})^2$. Now, because $\mathbf{R}_+$ is SPD (it is positive semidefinite by construction and non-singular by assumption), $\mathbf{x}^\top \mathbf{R}_+ \mathbf{x} > 0$ for any choice of $\mathbf{x} \neq \mathbf{0}$. The particular choice of $\mathbf{x} = [-\mathbf{u}^\top \mathbf{R}^{-1} \ \mathbf{1}]^\top$ implies that $s > 0$, which in turn implies that $R_+^2 \geq R^2$. $\square$

**Corollary 2.2.** Under the assumptions of proposition 2.1, the equality $R_+^2 = R^2$ holds if and only if $\gamma = \mathbf{u}^\top \mathbf{R}^{-1} \mathbf{r}$.

*Proof.* Straightforward from $R_+^2 = R^2 + s^{-1}(\gamma - \mathbf{u}^\top \mathbf{R}^{-1} \mathbf{r})^2$ with $s > 0$. $\square$

For the CV estimation of specific statistics, the expressions of $\mathbf{c}$ and $\mathbf{\Sigma}$ may be further reduced to involve statistics on $Y$ and $\mathbf{Z} := (Z_1 \cdots Z_M)^\top$ directly. These expressions, as well as the consequences on eqs. (2.6) and (2.7), are given in appendix A for the CV estimators of $\mathbb{E}[Y]$ and $\mathbb{V}[Y]$.

It should be noted at this point that, in practice, the optimal parameter $\alpha^*$ needs to be approximated. Specifically, the statistics involved in $\mathbf{c}$ and $\mathbf{\Sigma}$ need to be estimated. This can be done either using the same sample $\mathcal{X}$ as for the CV estimator itself, or using an independent, pilot sample $\mathcal{X}'$. The former strategy introduces a bias in the resulting CV estimator, while the latter guarantees unbiasedness but requires additional high-fidelity simulator evaluations, so that the former strategy is generally preferred. Furthermore, statistical remedies such as jackknifing, splitting or bootstrapping have been proposed to reduce or eliminate the bias introduced by the former strategy [38]. In both strategies, however, neither the theoretical variance reduction given by eq. (2.7) nor proposition 2.1 hold anymore.

**Remark 2.2.** In some instances, because it only involves the control variates, $\mathbf{\Sigma}$ may be known exactly (see remark 3.1 for an example in a multilevel setting involving PC-based control variates) or estimated accurately using many independent samples at negligible cost. Unfortunately, this is not the case for $\mathbf{c}$, which involves the output $Y = f(\mathbf{X})$ of the high-fidelity simulator $f$.

## 2.4 Surrogate-based control variates

The efficiency of the CV approach relies on the strong correlation between $Y = f(\mathbf{X})$ and the control variates $\mathbf{Z}$. In a multifidelity framework, the control variates correspond to the output of low-fidelity versions of $f$, typically in the form of simulators with simplified physics and/or coarser discretization. In many applications, the exact statistical measures $\tau$ of such control variates $\mathbf{Z}$ may not be available. One way to circumvent this limitation is to use additional samples to also estimate $\tau$. This type of estimators led to an approximate CV class of methods [24, 39, 43], where different model management and sample allocation strategies may be used to find a suitable trade-off between the additional cost of sampling $\mathbf{Z}$ and the resulting variance reduction. Recently, an optimal strategy was proposed [48, 49, 47] leading to the so-called multilevel best linear unbiased estimator (MLBLUE), as well as a unifying framework for a large class of multilevel and multifidelity estimators. However, ACV estimators typically yield lower variance reduction than their exact CV counterparts, as illustrated in appendix B on bi-fidelity (A)CV estimators of the expectation.

In this paper, we focus on the case where the low-fidelity models correspond to surrogate models of the high-fidelity simulator $f$, which allows us to use the original CV approach described in
section 2.3. Indeed, compared to the cost of evaluating the high-fidelity simulator $f$, the cost of evaluating surrogate models can be neglected, so that the statistical parameters $\tau$ of the output of the surrogate models can be estimated arbitrarily accurately at negligible cost. In fact, in some instances, the exact statistical parameters $\tau$ may even be directly available, as detailed hereafter. As a direct consequence, as discussed above, the question of sample allocation need not be addressed, and greater variance reduction is expected. Using surrogate models in a CV strategy has been explored previously in [16, 56], where the available computational budget is allocated both to the construction of the surrogates and to the actual CV estimations. In [56], the authors introduce an approach that optimally balances the computational effort needed to select the optimal degree of the polynomial chaos (PC) expansion used in a stochastic Galerkin CV approach. In our work, we focus on the situation in which the surrogate models are available, so that we do not consider any optimization strategy for the construction of the surrogates. Nevertheless, for the fairness of comparison, we still report the pre-processing cost of constructing the surrogate models.

We now briefly describe three different surrogate models commonly used in a UQ framework, and discuss the availability of statistics of their output.

### 2.4.1 Gaussian process modelling

We assume that the simulator $f$ is a realization of a Gaussian process (GP) $F$ indexed by $\mathbf{x}$ and defined by its mean function $m_F$ and covariance kernel $k_F$,

$$ \mathbb{E}[F(\mathbf{x})] = m_F(\mathbf{x}), \qquad (2.8) $$

$$ \mathbb{C}[F(\mathbf{x}), F(\mathbf{x}')] = k_F(\mathbf{x}, \mathbf{x}'), \quad \forall \mathbf{x}, \mathbf{x}' \in \Xi. \qquad (2.9) $$

In practice, one can parametrize the forms of $m_F$ and $k_F$. For example, in the widely used ordinary GP method, a stationary GP is assumed. In this case, $m_F$ is set as a constant $m_F(\mathbf{x}) = m$. More importantly, it is assumed that $k_F(\mathbf{x}, \mathbf{x}') = \bar{k}_F(\mathbf{x} - \mathbf{x}')$, and $k_F(\mathbf{x}, \mathbf{x}) = \bar{k}_F(0) = \sigma^2$ is a constant. Popular forms of kernels include polynomial, exponential, Gaussian, and Matérn functions. For example, the Gaussian kernel can be written as $k_F(\mathbf{x}, \mathbf{x}') = \sigma^2 \exp(-\frac{1}{2} \|\mathbf{x} - \mathbf{x}'\|_h^2)$, where the weighted norm is defined as $\|\mathbf{x} - \mathbf{x}'\|_h = (\sum_{i=1}^d \frac{(x_i - x_i')^2}{h_i^2})^{1/2}$ where $h_1, \dots, h_d$ are correlation lengths. The hyperparameters $\sigma$ and $h_i$ can be obtained by maximum likelihood. Then, given $n$ observations $\mathcal{F} = (f(\mathbf{x}^{(1)}), \dots, f(\mathbf{x}^{(n)}))^{\top}$ of $F$ at $\mathcal{X} = (\mathbf{x}^{(1)}, \dots, \mathbf{x}^{(n)})^{\top}$, the posterior $\tilde{F}$ of $F$ can be defined as

$$ \tilde{F} = [F \mid F(\mathcal{X}) = \mathcal{F}], \qquad (2.10) $$

whose expectation and covariance are given by

$$ m_{\tilde{F}}(\mathbf{x}) = \mathbb{E}[\tilde{F}(\mathbf{x})] = m_F(\mathbf{x}) + k_F(\mathbf{x}, \mathcal{X})^\top k_F(\mathcal{X}, \mathcal{X})^{-1} (\mathcal{F} - m_F(\mathcal{X})), \qquad (2.11) $$

$$ k_{\tilde{F}}(\mathbf{x}, \mathbf{x}') = \mathbb{C}[\tilde{F}(\mathbf{x}), \tilde{F}(\mathbf{x}')] = k_F(\mathbf{x}, \mathbf{x}') - k_F(\mathbf{x}, \mathcal{X})^\top k_F(\mathcal{X}, \mathcal{X})^{-1} k_F(\mathbf{x}, \mathbf{x}'). \qquad (2.12) $$

Thereafter, the high-fidelity model $f$ at $\mathbf{x}$ will be approximated by the conditional expectation,

$$ g^{\text{GP}}(\mathbf{x}) = m_{\tilde{F}}(\mathbf{x}). \qquad (2.13) $$

Thus, the expectation and variance of $g^{\text{GP}}(\mathbf{X})$ are defined as

$$ \mathbb{E}[g^{\text{GP}}(\mathbf{X})] = \int_{\Xi} m_{\tilde{F}}(\mathbf{x}) p_{\mathbf{X}}(\mathbf{x}) \mathrm{d}\mathbf{x}, \qquad (2.14) $$

$$ \mathbb{V}[g^{\text{GP}}(\mathbf{X})] = \int_{\Xi} (m_{\tilde{F}}(\mathbf{x}) - \mathbb{E}[g^{\text{GP}}(\mathbf{X})])^2 p_{\mathbf{X}}(\mathbf{x}) \mathrm{d}\mathbf{x}. \qquad (2.15) $$
section 2.3. Indeed, compared to the cost of evaluating the high-fidelity simulator $f$, the cost of evaluating surrogate models can be neglected, so that the statistical parameters $\tau$ of the output of the surrogate models can be estimated arbitrarily accurately at negligible cost. In fact, in some instances, the exact statistical parameters $\tau$ may even be directly available, as detailed hereafter. As a direct consequence, as discussed above, the question of sample allocation need not be addressed, and greater variance reduction is expected. Using surrogate models in a CV strategy has been explored previously in [16, 56], where the available computational budget is allocated both to the construction of the surrogates and to the actual CV estimations. In [56], the authors introduce an approach that optimally balances the computational effort needed to select the optimal degree of the polynomial chaos (PC) expansion used in a stochastic Galerkin CV approach. In our work, we focus on the situation in which the surrogate models are available, so that we do not consider any optimization strategy for the construction of the surrogates. Nevertheless, for the fairness of comparison, we still report the pre-processing cost of constructing the surrogate models.

We now briefly describe three different surrogate models commonly used in a UQ framework, and discuss the availability of statistics of their output.

### 2.4.1 Gaussian process modelling

We assume that the simulator $f$ is a realization of a Gaussian process (GP) $F$ indexed by $\mathbf{x}$ and defined by its mean function $m_F$ and covariance kernel $k_F$,

$$ \mathbb{E}[F(\mathbf{x})] = m_F(\mathbf{x}), \qquad (2.8) $$

$$ \mathbb{C}[F(\mathbf{x}), F(\mathbf{x}')] = k_F(\mathbf{x}, \mathbf{x}'), \quad \forall \mathbf{x}, \mathbf{x}' \in \Xi. \qquad (2.9) $$

In practice, one can parametrize the forms of $m_F$ and $k_F$. For example, in the widely used ordinary GP method, a stationary GP is assumed. In this case, $m_F$ is set as a constant $m_F(\mathbf{x}) = m$. More importantly, it is assumed that $k_F(\mathbf{x}, \mathbf{x}') = \bar{k}_F(\mathbf{x} - \mathbf{x}')$, and $k_F(\mathbf{x}, \mathbf{x}) = \bar{k}_F(0) = \sigma^2$ is a constant. Popular forms of kernels include polynomial, exponential, Gaussian, and Matérn functions. For example, the Gaussian kernel can be written as $k_F(\mathbf{x}, \mathbf{x}') = \sigma^2 \exp(-\frac{1}{2}\|\mathbf{x} - \mathbf{x}'\|_h^2)$, where the weighted norm is defined as $\|\mathbf{x} - \mathbf{x}'\|_h = \left(\sum_{i=1}^d \frac{(x_i - x_i')^2}{h_i^2}\right)^{1/2}$ where $h_1, \dots, h_d$ are correlation lengths. The hyperparameters $\sigma$ and $h_i$ can be obtained by maximum likelihood. Then, given $n$ observations $\mathcal{F} = (f(\mathbf{x}^{(1)}), \dots, f(\mathbf{x}^{(n)}))^{\top}$ of $F$ at $\mathcal{X} = (\mathbf{x}^{(1)}, \dots, \mathbf{x}^{(n)})^{\top}$, the posterior $\tilde{F}$ of $F$ can be defined as

$$ \tilde{F} = [F \mid F(\mathcal{X}) = \mathcal{F}], \qquad (2.10) $$

whose expectation and covariance are given by

$$ m_{\tilde{F}}(\mathbf{x}) = \mathbb{E}[\tilde{F}(\mathbf{x})] = m_F(\mathbf{x}) + k_F(\mathbf{x}, \mathcal{X})^\top k_F(\mathcal{X}, \mathcal{X})^{-1} (\mathcal{F} - m_F(\mathcal{X})), \qquad (2.11) $$

$$ k_{\tilde{F}}(\mathbf{x}, \mathbf{x}') = \mathbb{C}[\tilde{F}(\mathbf{x}), \tilde{F}(\mathbf{x}')] = k_F(\mathbf{x}, \mathbf{x}') - k_F(\mathbf{x}, \mathcal{X})^\top k_F(\mathcal{X}, \mathcal{X})^{-1} k_F(\mathbf{x}, \mathbf{x}'). \qquad (2.12) $$

Thereafter, the high-fidelity model $f$ at $\mathbf{x}$ will be approximated by the conditional expectation,

$$ g^{\text{GP}}(\mathbf{x}) = m_{\tilde{F}}(\mathbf{x}). \qquad (2.13) $$

Thus, the expectation and variance of $g^{\text{GP}}(\mathbf{X})$ are defined as

$$ \mathbb{E}[g^{\text{GP}}(\mathbf{X})] = \int_{\Xi} m_{\tilde{F}}(\mathbf{x}) p_{\mathbf{X}}(\mathbf{x}) \mathrm{d}\mathbf{x}, \qquad (2.14) $$

$$ \mathbb{V}[g^{\text{GP}}(\mathbf{X})] = \int_{\Xi} (m_{\tilde{F}}(\mathbf{x}) - \mathbb{E}[g^{\text{GP}}(\mathbf{X})])^2 p_{\mathbf{X}}(\mathbf{x}) \mathrm{d}\mathbf{x}. \qquad (2.15) $$
These two statistics can be approximated empirically by taking a large sample of $\mathbf{X}$, as the GP model is inexpensive. Analytical formulas exist for some pairs of distributions and covariance functions (see [8, Table 1] and [28, 15]).

### 2.4.2 Taylor polynomials

We assume that the numerical simulator $f$ is infinitely differentiable and that the moments of $\mathbf{X}$ are finite. Under this assumption, it is possible to expand the original simulator $f$ around the input's expected value $\mu_{\mathbf{X}} := \mathbb{E}[\mathbf{X}]$ as the infinite polynomial series according to Taylor's theorem,

$$f(\mathbf{X}) = g^{\mathrm{T}\infty}(\mathbf{X}) = \sum_{|\beta| \le p} \frac{(\mathbf{X} - \mu_{\mathbf{X}})^{\beta}}{\beta!} D^{\beta} f(\mu_{\mathbf{X}}), \quad (2.16)$$

where $\beta \in \mathbb{N}^d$, $|\beta| = \sum_{i=1}^d \beta_i$, $\beta! = \prod_{i=1}^d \beta_i!$, $x^\beta = \prod_{i=1}^d x_i^{\beta_i}$ and $D^\beta f = \frac{\partial^{|\beta|} f}{\partial \beta_1 x_1 \cdots \partial \beta_d x_d}$. Thus, the function $f$ may be approximated by the first- and second-order Taylor polynomials,

$$g^{\mathrm{T}_1}(\mathbf{X}) = f(\mu_{\mathbf{X}}) + \mathbf{J}_f(\mu_{\mathbf{X}})(\mathbf{X} - \mu_{\mathbf{X}}), \quad (2.17)$$

$$g^{\mathrm{T}_2}(\mathbf{X}) = f(\mu_{\mathbf{X}}) + \mathbf{J}_f(\mu_{\mathbf{X}})(\mathbf{X} - \mu_{\mathbf{X}}) + \frac{1}{2}(\mathbf{X} - \mu_{\mathbf{X}})^{\mathrm{T}} \mathbf{H}_f(\mu_{\mathbf{X}})(\mathbf{X} - \mu_{\mathbf{X}}), \quad (2.18)$$

where $\mathbf{J}_f(\mu_{\mathbf{X}}) = \nabla_f(\mu_{\mathbf{X}})^\mathrm{T} \in \mathbb{R}^{1 \times d}$ and $\mathbf{H}_f(\mu_{\mathbf{X}}) \in \mathbb{R}^{d \times d}$ are the Jacobian and Hessian matrices of $f$ at $\mu_{\mathbf{X}}$, respectively. In practical computations, the derivatives may be approximated by numerical differentiation if they are not provided. The expectation and variance of $g^{\mathrm{T}_1}(\mathbf{X})$ and $g^{\mathrm{T}_2}(\mathbf{X})$ are then defined as

$$\mathbb{E}[g^{\mathrm{T}_1}(\mathbf{X})] = f(\mu_{\mathbf{X}}), \quad (2.19)$$

$$\mathbb{V}[g^{\mathrm{T}_1}(\mathbf{X})] = \mathbf{J}_f(\mu_{\mathbf{X}})^{\odot 2} \sigma_{\mathbf{X}}^2, \quad (2.20)$$

$$\mathbb{E}[g^{\mathrm{T}_2}(\mathbf{X})] = f(\mu_{\mathbf{X}}) + \frac{1}{2} \mathrm{Tr}(\mathbf{H}_f(\mu_{\mathbf{X}}) \Sigma_{\mathbf{X}}^2), \quad (2.21)$$

$$\mathbb{V}[g^{\mathrm{T}_2}(\mathbf{X})] = \mathrm{Tr}(\mathbf{J}_f(\mu_{\mathbf{X}})^{\mathrm{T}} \mathbf{J}_f(\mu_{\mathbf{X}}) \Sigma_{\mathbf{X}}^2) + \frac{1}{2} \mathrm{Tr}(\mathbf{H}_f(\mu_{\mathbf{X}})^{\odot 2} \sigma_{\mathbf{X}}^2 (\sigma_{\mathbf{X}}^2)^{\mathrm{T}}), \quad (2.22)$$

where for any matrix $\mathbf{A}$, $\mathbf{A}^{\odot 2}$ denotes the element-wise square of $\mathbf{A}$, $\mathrm{Tr}(\mathbf{A})$ is the trace of $\mathbf{A}$, $\sigma_{\mathbf{X}}^2 = [\sigma_{X_1}^2 \cdots \sigma_{X_d}^2]^{\mathrm{T}} \in \mathbb{R}^d$ is the element-wise variance of $\mathbf{X}$, and $\Sigma_{\mathbf{X}}^2 = \mathrm{Diag}(\sigma_{\mathbf{X}}^2) \in \mathbb{R}^{d \times d}$.

### 2.4.3 Polynomial chaos expansion

We consider the truncated polynomial chaos (PC) expansion of $f(\mathbf{X})$ [55, 9, 20, 32],

$$f(\mathbf{X}) \approx g^{\mathrm{PC}_P}(\mathbf{X}) = \sum_{k=0}^{P} g_k \Psi_k(\mathbf{X}), \quad (2.23)$$

where the coefficients $g_k$ are real scalars and $\Psi_k$ are orthonormal multivariate polynomials:

$$\forall i, j \ge 0, \quad \langle \Psi_i, \Psi_j \rangle_{p_X} := \int_{\Xi} \Psi_i(\mathbf{x}) \Psi_j(\mathbf{x}) p_X(\mathbf{x}) \mathrm{d}\mathbf{x} = \mathbb{E}[\Psi_i(\mathbf{X}) \Psi_j(\mathbf{X})] = \delta_{ij}, \quad (2.24)$$

where $\delta_{ij}$ denotes the Kronecker delta. The coefficients $g_k$ may be approximated using non-intrusive techniques such as non-intrusive (pseudo)spectral projection [46, 12, 11], regression [3, 51, 6, 7],
These two statistics can be approximated empirically by taking a large sample of $\mathbf{X}$, as the GP model is inexpensive. Analytical formulas exist for some pairs of distributions and covariance functions (see [8, Table 1] and [28, 15]).

### 2.4.2 Taylor polynomials

We assume that the numerical simulator $f$ is infinitely differentiable and that the moments of $\mathbf{X}$ are finite. Under this assumption, it is possible to expand the original simulator $f$ around the input's expected value $\mu_{\mathbf{X}} := \mathbb{E}[\mathbf{X}]$ as the infinite polynomial series according to Taylor's theorem,

$$f(\mathbf{X}) = g^{\mathrm{T}\infty}(\mathbf{X}) = \sum_{|\beta| \le p} \frac{(\mathbf{X} - \mu_{\mathbf{X}})^{\beta}}{\beta!} D^{\beta} f(\mu_{\mathbf{X}}), \quad (2.16)$$

where $\beta \in \mathbb{N}^d$, $|\beta| = \sum_{i=1}^d \beta_i$, $\beta! = \prod_{i=1}^d \beta_i!$, $x^\beta = \prod_{i=1}^d x_i^{\beta_i}$ and $D^\beta f = \frac{\partial^{|\beta|} f}{\partial \beta_1 x_1 \cdots \partial \beta_d x_d}$. Thus, the function $f$ may be approximated by the first- and second-order Taylor polynomials,

$$g^{\mathrm{T}_1}(\mathbf{X}) = f(\mu_{\mathbf{X}}) + \mathbf{J}_f(\mu_{\mathbf{X}})(\mathbf{X} - \mu_{\mathbf{X}}), \quad (2.17)$$

$$g^{\mathrm{T}_2}(\mathbf{X}) = f(\mu_{\mathbf{X}}) + \mathbf{J}_f(\mu_{\mathbf{X}})(\mathbf{X} - \mu_{\mathbf{X}}) + \frac{1}{2}(\mathbf{X} - \mu_{\mathbf{X}})^{\mathrm{T}} \mathbf{H}_f(\mu_{\mathbf{X}})(\mathbf{X} - \mu_{\mathbf{X}}), \quad (2.18)$$

where $\mathbf{J}_f(\mu_{\mathbf{X}}) = \nabla_f(\mu_{\mathbf{X}})^{\mathrm{T}} \in \mathbb{R}^{1 \times d}$ and $\mathbf{H}_f(\mu_{\mathbf{X}}) \in \mathbb{R}^{d \times d}$ are the Jacobian and Hessian matrices of $f$ at $\mu_{\mathbf{X}}$, respectively. In practical computations, the derivatives may be approximated by numerical differentiation if they are not provided. The expectation and variance of $g^{\mathrm{T}_1}(\mathbf{X})$ and $g^{\mathrm{T}_2}(\mathbf{X})$ are then defined as

$$\mathbb{E}[g^{\mathrm{T}_1}(\mathbf{X})] = f(\mu_{\mathbf{X}}), \quad (2.19)$$

$$\mathbb{V}[g^{\mathrm{T}_1}(\mathbf{X})] = \mathbf{J}_f(\mu_{\mathbf{X}})^{\odot 2} \sigma_{\mathbf{X}}^2, \quad (2.20)$$

$$\mathbb{E}[g^{\mathrm{T}_2}(\mathbf{X})] = f(\mu_{\mathbf{X}}) + \frac{1}{2} \mathrm{Tr}(\mathbf{H}_f(\mu_{\mathbf{X}}) \Sigma_{\mathbf{X}}^2), \quad (2.21)$$

$$\mathbb{V}[g^{\mathrm{T}_2}(\mathbf{X})] = \mathrm{Tr}(\mathbf{J}_f(\mu_{\mathbf{X}})^{\mathrm{T}} \mathbf{J}_f(\mu_{\mathbf{X}}) \Sigma_{\mathbf{X}}^2) + \frac{1}{2} \mathrm{Tr}(\mathbf{H}_f(\mu_{\mathbf{X}})^{\odot 2} \sigma_{\mathbf{X}}^2 (\sigma_{\mathbf{X}}^2)^{\mathrm{T}}), \quad (2.22)$$

where for any matrix $\mathbf{A}$, $\mathbf{A}^{\odot 2}$ denotes the element-wise square of $\mathbf{A}$, $\mathrm{Tr}(\mathbf{A})$ is the trace of $\mathbf{A}$, $\sigma_{\mathbf{X}}^2 = [\sigma_{X_1}^2 \cdots \sigma_{X_d}^2]^{\mathrm{T}} \in \mathbb{R}^d$ is the element-wise variance of $\mathbf{X}$, and $\Sigma_{\mathbf{X}}^2 = \mathrm{Diag}(\sigma_{\mathbf{X}}^2) \in \mathbb{R}^{d \times d}$.

### 2.4.3 Polynomial chaos expansion

We consider the truncated polynomial chaos (PC) expansion of $f(\mathbf{X})$ [55, 9, 20, 32],

$$f(\mathbf{X}) \approx g^{\mathrm{PC}_P}(\mathbf{X}) = \sum_{k=0}^{P} g_k \Psi_k(\mathbf{X}), \quad (2.23)$$

where the coefficients $g_k$ are real scalars and $\Psi_k$ are orthonormal multivariate polynomials:

$$\forall i, j \ge 0, \quad \langle \Psi_i, \Psi_j \rangle_{p_X} := \int_{\Xi} \Psi_i(\mathbf{x}) \Psi_j(\mathbf{x}) p_X(\mathbf{x}) \mathrm{d}\mathbf{x} = \mathbb{E}[\Psi_i(\mathbf{X}) \Psi_j(\mathbf{X})] = \delta_{ij}, \quad (2.24)$$

where $\delta_{ij}$ denotes the Kronecker delta. The coefficients $g_k$ may be approximated using non-intrusive techniques such as non-intrusive (pseudo)spectral projection [46, 12, 11], regression [3, 51, 6, 7],
interpolation [2, 40], or from intrusive approaches such as the stochastic Galerkin method [20, 32]. Assuming $\Psi_0 \equiv 1$, the expectation and variance of $g^{\text{PCP}}(\mathbf{X})$ are then given by

$$ \mathbb{E}[g^{\text{PCP}}(\mathbf{X})] = \mathfrak{g}_0 \quad \text{and} \quad \mathbb{V}[g^{\text{PCP}}(\mathbf{X})] = \sum_{k=1}^{P} \mathfrak{g}_k^2. \qquad (2.25) $$

# 3 Multilevel estimators

In this section, we present so-called multilevel statistical estimation techniques based on a sequence of simulators ($f_\ell$)$_{\ell=0}^L$, with increasing accuracy and computational cost, where $f_L \equiv f$ denotes the high-fidelity numerical simulator. The levels are ordered from the coarsest ($\ell = 0$) to the finest ($\ell = L$). We denote by $Y_\ell$ the random variable $Y_\ell = f_\ell(\mathbf{X})$ and $\theta_\ell$ the statistic of $f_\ell(\mathbf{X})$ increasingly close to $\theta_L \equiv \theta$.

## 3.1 Multilevel Monte Carlo

The statistic $\theta_L$ can be expressed as the telescoping sum

$$ \theta_L = \sum_{\ell=0}^{L} T_\ell, \qquad (3.1) $$

where $T_\ell = \theta_\ell - \theta_{\ell-1}$, and, by convention, $\theta_{-1} := 0$. The MLMC estimator $\hat{\theta}_L^{\text{MLMC}}$ of $\theta_L$ is then defined as [21, 22, 53]

$$ \hat{\theta}^{\text{MLMC}} = \sum_{\ell=0}^{L} \hat{T}_{\ell}^{(\ell)}, \qquad (3.2) $$

where, at each level $\ell$, $\hat{T}_\ell^{(\ell)}$ is an unbiased MC estimator of $T_\ell$, based on an input sample $\mathcal{X}^{(\ell)} = \{\mathbf{X}^{(\ell,i)}\}_{i=1}^{n_\ell}$ such that the members of $\mathcal{X}^{(\ell)}$ and $\mathcal{X}^{(\ell')}$ are mutually independent for $\ell \neq \ell'$. In many instances, $\hat{T}_\ell^{(\ell)}$ is actually defined as $\hat{T}_\ell^{(\ell)} = \hat{\theta}_\ell^{(\ell)} - \hat{\theta}_{\ell-1}^{(\ell)}$, where $\hat{\theta}_k^{(\ell)}$ denotes the unbiased MC estimator of $\theta_k$ based on the simulator $f_k$ using the $n_\ell$-sample $\mathcal{X}^{(\ell)}$.

For example, the MLMC estimator $\hat{E}^{\text{MLMC}}[Y]$ of the expectation $\mathbb{E}[Y]$ is given by

$$ \hat{E}^{\text{MLMC}}[Y] = \hat{E}^{(0)}[Y_0] + \sum_{\ell=1}^{L} \hat{E}^{(\ell)}[Y_\ell] - \hat{E}^{(\ell)}[Y_{\ell-1}], \qquad (3.3) $$

where $\hat{E}^{(\ell)}[Y_k] = n_\ell^{-1} \sum_{i=1}^{n_\ell} f_k(\mathbf{X}^{(\ell,i)}),$ with $k \in \{\ell, \ell-1\}$. We stress that the correction terms at each level $\ell$ are computed from the same input sample $\mathcal{X}^{(\ell)}$, but using two successive simulators, $f_\ell$ and $f_{\ell-1}$. Similarly, the MLMC estimator $\hat{V}^{\text{MLMC}}[Y]$ of the variance $\mathbb{V}[Y]$ is defined as [4]

$$ \hat{V}^{\text{MLMC}}[Y] = \hat{V}^{(0)}[Y_0] + \sum_{\ell=1}^{L} \hat{V}^{(\ell)}[Y_\ell] - \hat{V}^{(\ell)}[Y_{\ell-1}], \qquad (3.4) $$

where $\hat{V}^{(\ell)}[Y_k] = \frac{n_\ell}{n_\ell-1}(\hat{E}^{(\ell)}[Y_k^2] - \hat{E}^{(\ell)}[Y_k]^2)$ is the single-level unbiased MC variance estimator. Owing to the independence of the estimators $\hat{T}_\ell^{(\ell)}$, the variance of the MLMC estimator is

$$ \mathbb{V}[\hat{\theta}^{\text{MLMC}}] = \sum_{\ell=0}^{L} \mathbb{V}[\hat{T}_{\ell}^{(\ell)}]. \qquad (3.5) $$
interpolation [2, 40], or from intrusive approaches such as the stochastic Galerkin method [20, 32]. Assuming $\Psi_0 \equiv 1$, the expectation and variance of $g^{PC_P}(\mathbf{X})$ are then given by

$$ \mathbb{E}[g^{PC_P}(\mathbf{X})] = \mathfrak{g}_0 \quad \text{and} \quad \mathbb{V}[g^{PC_P}(\mathbf{X})] = \sum_{k=1}^{P} \mathfrak{g}_k^2. \qquad (2.25) $$

# 3 Multilevel estimators

In this section, we present so-called multilevel statistical estimation techniques based on a sequence of simulators ($f_\ell$)$_{\ell=0}^L$, with increasing accuracy and computational cost, where $f_L \equiv f$ denotes the high-fidelity numerical simulator. The levels are ordered from the coarsest ($\ell = 0$) to the finest ($\ell = L$). We denote by $Y_\ell$ the random variable $Y_\ell = f_\ell(\mathbf{X})$ and $\theta_\ell$ the statistic of $f_\ell(\mathbf{X})$ increasingly close to $\theta_L \equiv \theta$.

## 3.1 Multilevel Monte Carlo

The statistic $\theta_L$ can be expressed as the telescoping sum

$$ \theta_L = \sum_{\ell=0}^{L} T_{\ell}, \qquad (3.1) $$

where $T_\ell = \theta_\ell - \theta_{\ell-1}$, and, by convention, $\theta_{-1} := 0$. The MLMC estimator $\hat{\theta}_L^{\text{MLMC}}$ of $\theta_L$ is then defined as [21, 22, 53]

$$ \hat{\theta}^{\text{MLMC}} = \sum_{\ell=0}^{L} \hat{T}_{\ell}^{(\ell)}, \qquad (3.2) $$

where, at each level $\ell$, $\hat{T}_\ell^{(\ell)}$ is an unbiased MC estimator of $T_\ell$, based on an input sample $\mathcal{X}^{(\ell)} = \{\mathbf{X}^{(\ell,i)}\}_{i=1}^{n_\ell}$ such that the members of $\mathcal{X}^{(\ell)}$ and $\mathcal{X}^{(\ell')}$ are mutually independent for $\ell \neq \ell'$. In many instances, $\hat{T}_\ell^{(\ell)}$ is actually defined as $\hat{T}_\ell^{(\ell)} = \hat{\theta}_\ell^{(\ell)} - \hat{\theta}_{\ell-1}^{(\ell)}$, where $\hat{\theta}_k^{(\ell)}$ denotes the unbiased MC estimator of $\theta_k$ based on the simulator $f_k$ using the $n_\ell$-sample $\mathcal{X}^{(\ell)}$.

For example, the MLMC estimator $\hat{E}^{\text{MLMC}}[Y]$ of the expectation $\mathbb{E}[Y]$ is given by

$$ \hat{E}^{\text{MLMC}}[Y] = \hat{E}^{(0)}[Y_0] + \sum_{\ell=1}^{L} \hat{E}^{(\ell)}[Y_\ell] - \hat{E}^{(\ell)}[Y_{\ell-1}], \qquad (3.3) $$

where $\hat{E}^{(\ell)}[Y_k] = n_\ell^{-1} \sum_{i=1}^{n_\ell} f_k(\mathbf{X}^{(\ell,i)})$, with $k \in \{\ell, \ell-1\}$. We stress that the correction terms at each level $\ell$ are computed from the same input sample $\mathcal{X}^{(\ell)}$, but using two successive simulators, $f_\ell$ and $f_{\ell-1}$. Similarly, the MLMC estimator $\hat{V}^{\text{MLMC}}[Y]$ of the variance $\mathbb{V}[Y]$ is defined as [4]

$$ \hat{V}^{\text{MLMC}}[Y] = \hat{V}^{(0)}[Y_0] + \sum_{\ell=1}^{L} \hat{V}^{(\ell)}[Y_\ell] - \hat{V}^{(\ell)}[Y_{\ell-1}], \qquad (3.4) $$

where $\hat{V}^{(\ell)}[Y_k] = \frac{n_\ell}{n_\ell-1}(\hat{E}^{(\ell)}[Y_k^2] - \hat{E}^{(\ell)}[Y_k]^2)$ is the single-level unbiased MC variance estimator. Owing to the independence of the estimators $\hat{T}_\ell^{(\ell)}$, the variance of the MLMC estimator is

$$ \mathbb{V}[\hat{\theta}^{\text{MLMC}}] = \sum_{\ell=0}^{L} \mathbb{V}[\hat{T}_{\ell}^{(\ell)}]. \qquad (3.5) $$
In practice, the MLMC method relies on the allocation of the total (expected) computational cost of the MLMC estimator across the different levels, with

$$ \mathrm{cost}(\hat{\theta}^{\mathrm{MLMC}}) = \sum_{\ell=0}^{L} n_{\ell}(C_{\ell} + C_{\ell-1}), \qquad (3.6) $$

where $C_\ell$ is the (expected) cost of one evaluation of the simulator $f_\ell$, with $C_{-1} := 0$ by convention. Thus, a key aspect is played by the choice of the number of samples $n_\ell$ allocated to each level $\ell$. The goal is to find the sample sizes $n_0, \dots, n_L$ that minimize the variance of the estimator $\mathbb{V}[\hat{\theta}^{\mathrm{MLMC}}]$ given a computational budget $C$. Thus, the sample allocation problem reads

$$ \begin{array}{ll} \underset{n_0, \dots, n_L \in \mathbb{N}^*}{\mathrm{minimize}} & \mathbb{V}[\hat{\theta}^{\mathrm{MLMC}}] \\ \text{subject to} & \mathrm{cost}(\hat{\theta}^{\mathrm{MLMC}}) = C. \end{array} \qquad (3.7) $$

In practice, $\mathbb{V}[\hat{\theta}^{\mathrm{MLMC}}]$ is not known, and we instead rely on the assumption that $\mathbb{V}[\hat{T}_\ell^{(\ell)}] \lesssim n_\ell^{-1} \mathcal{V}_\ell$, with $\mathcal{V}_\ell$ independent of $n_\ell$ [53, 36]. This is a reasonable assumption that holds for the MLMC estimation of the expectation, variance and covariance [36, Table 1]. Note that, for the estimation of the expectation, we have $\mathbb{V}[\hat{T}_\ell^{(\ell)}] = n_\ell^{-1} \mathcal{V}_\ell$, with $\mathcal{V}_\ell = \mathbb{V}[Y_\ell - Y_{\ell-1}]$ and $Y_{-1} \equiv 0$. The sample allocation problem eq. (3.7) is then replaced with

$$ \begin{array}{ll} \underset{n_0, \dots, n_L \in \mathbb{N}^*}{\mathrm{minimize}} & \displaystyle\sum_{\ell=0}^{L} n_{\ell}^{-1} \mathcal{V}_{\ell} \\ \text{subject to} & \mathrm{cost}(\hat{\theta}^{\mathrm{MLMC}}) = C, \end{array} \qquad (3.8) $$

which is equivalent for the expectation, and an approximation for other statistics. This minimization problem has a unique solution which can be computed analytically (see, e.g., [21, 10, 36]).

## 3.2 Multilevel surrogate-based control variate strategies

In this section, we introduce various surrogate-based control variate strategies in a multilevel framework where a hierarchy of simulators ($f_\ell$)$_{\ell=0}^L$ is available. These strategies rely on using the random variables ($g_\ell(\mathbf{X})$)$_{\ell=0}^L$ as control variates, where $g_\ell$ is a surrogate model of the simulator $f_\ell$.

### 3.2.1 Multilevel control variates (MLCV)

The first strategy, referred to as multilevel control variates and hereafter abbreviated MLCV, consists of using the surrogate models at all levels to build the control variates in eq. (2.4). Thus, $\tau_\ell$ corresponds to the statistical measure of the random variable $Z_\ell = g_\ell(\mathbf{X})$, and $\hat{\tau}_\ell$ to its unbiased MC estimator. For instance, the MLCV estimator of the expectation based on an $n_L$-sample $\mathcal{X}^{(L)} = \{\mathbf{X}^{(1)}, \dots, \mathbf{X}^{(n_L)}\}$ reads

$$ \hat{E}^{\mathrm{MLCV}}[Y](\boldsymbol{\alpha}) = \hat{E}^{(L)}[Y_L] - \sum_{\ell=0}^{L} \alpha_{\ell} \left( \hat{E}^{(L)}[Z_{\ell}] - \mu_{Z_{\ell}} \right), \qquad (3.9) $$

where $\hat{E}^{(L)}[Y_L] = n_L^{-1} \sum_{i=1}^{n_L} f_L(\mathbf{X}^{(i)}), \hat{E}^{(L)}[Z_\ell] = n_L^{-1} \sum_{i=1}^{n_L} g_\ell(\mathbf{X}^{(i)}),$ and $\mu_{Z_\ell} = \mathbb{E}[g_\ell(\mathbf{X})]$. Note that this approach does not build on the MLMC methodology described in section 3.1, but still exploits multilevel information through the surrogates constructed at different levels.
In practice, the MLMC method relies on the allocation of the total (expected) computational cost of the MLMC estimator across the different levels, with

$$ \text{cost}(\hat{\theta}^{\text{MLMC}}) = \sum_{\ell=0}^{L} n_{\ell}(C_{\ell} + C_{\ell-1}), \qquad (3.6) $$

where $C_\ell$ is the (expected) cost of one evaluation of the simulator $f_\ell$, with $C_{-1} := 0$ by convention. Thus, a key aspect is played by the choice of the number of samples $n_\ell$ allocated to each level $\ell$. The goal is to find the sample sizes $n_0, \dots, n_L$ that minimize the variance of the estimator $\mathbb{V}[\hat{\theta}^{\text{MLMC}}]$ given a computational budget $C$. Thus, the sample allocation problem reads

$$ \begin{array}{ll} \underset{n_0, \dots, n_L \in \mathbb{N}^*}{\text{minimize}} & \mathbb{V}[\hat{\theta}^{\text{MLMC}}] \\ \text{subject to} & \text{cost}(\hat{\theta}^{\text{MLMC}}) = C. \end{array} \qquad (3.7) $$

In practice, $\mathbb{V}[\hat{\theta}^{\text{MLMC}}]$ is not known, and we instead rely on the assumption that $\mathbb{V}[\hat{T}_\ell^{(\ell)}] \lesssim n_\ell^{-1}V_\ell$, with $V_\ell$ independent of $n_\ell$ [53, 36]. This is a reasonable assumption that holds for the MLMC estimation of the expectation, variance and covariance [36, Table 1]. Note that, for the estimation of the expectation, we have $\mathbb{V}[\hat{T}_\ell^{(\ell)}] = n_\ell^{-1}V_\ell$, with $V_\ell = \mathbb{V}[Y_\ell - Y_{\ell-1}]$ and $Y_{-1} \equiv 0$. The sample allocation problem eq. (3.7) is then replaced with

$$ \begin{array}{ll} \underset{n_0, \dots, n_L \in \mathbb{N}^*}{\text{minimize}} & \displaystyle\sum_{\ell=0}^{L} n_{\ell}^{-1} V_{\ell} \\ \text{subject to} & \text{cost}(\hat{\theta}^{\text{MLMC}}) = C, \end{array} \qquad (3.8) $$

which is equivalent for the expectation, and an approximation for other statistics. This minimization problem has a unique solution which can be computed analytically (see, e.g., [21, 10, 36]).

## 3.2 Multilevel surrogate-based control variate strategies

In this section, we introduce various surrogate-based control variate strategies in a multilevel framework where a hierarchy of simulators ($f_\ell$)$_{\ell=0}^L$ is available. These strategies rely on using the random variables ($g_\ell(\mathbf{X})$)$_{\ell=0}^L$ as control variates, where $g_\ell$ is a surrogate model of the simulator $f_\ell$.

### 3.2.1 Multilevel control variates (MLCV)

The first strategy, referred to as multilevel control variates and hereafter abbreviated MLCV, consists of using the surrogate models at all levels to build the control variates in eq. (2.4). Thus, $\tau_\ell$ corresponds to the statistical measure of the random variable $Z_\ell = g_\ell(\mathbf{X})$, and $\hat{\tau}_\ell$ to its unbiased MC estimator. For instance, the MLCV estimator of the expectation based on an $n_L$-sample $\mathcal{X}^{(L)} = \{\mathbf{X}^{(1)}, \dots, \mathbf{X}^{(n_L)}\}$ reads

$$ \hat{E}^{\text{MLCV}}[Y](\boldsymbol{\alpha}) = \hat{E}^{(L)}[Y_L] - \sum_{\ell=0}^{L} \alpha_{\ell} \left( \hat{E}^{(L)}[Z_{\ell}] - \mu_{Z_{\ell}} \right), \qquad (3.9) $$

where $\hat{E}^{(L)}[Y_L] = n_L^{-1} \sum_{i=1}^{n_L} f_L(\mathbf{X}^{(i)}), \hat{E}^{(L)}[Z_\ell] = n_L^{-1} \sum_{i=1}^{n_L} g_\ell(\mathbf{X}^{(i)}),$ and $\mu_{Z_\ell} = \mathbb{E}[g_\ell(\mathbf{X})]$. Note that this approach does not build on the MLMC methodology described in section 3.1, but still exploits multilevel information through the surrogates constructed at different levels.
### 3.2.2 Multilevel Monte Carlo with control variates (MLMC-CV)

The second strategy consists of improving the MLMC estimator eq. (3.2) using the surrogate-based control variates $Z_0, \dots, Z_L$. Specifically, the MLMC-CV estimator improves the MC estimation of each of the correction terms of the MLMC estimator by using a surrogate model of the corresponding correction term as a control variate. Note that the MLMF approach proposed in [18, 19] is based on a similar strategy, using arbitrary low-fidelity models at level in an approximate CV setting. The MLMC-CV estimator reads

$$ \hat{\theta}^{\text{MLMC-CV}}(\alpha_1, \dots, \alpha_L) = \sum_{\ell=0}^{L} \hat{T}_{\ell}^{\text{CV}}(\alpha_{\ell}), \qquad (3.10) $$

where $\hat{T}_{\ell}^{\text{CV}}(\alpha_{\ell})$ is the CV estimator of the multilevel correction term $T_{\ell}$ (cf. eq. (3.11)),

$$ \hat{T}_{\ell}^{\text{CV}}(\alpha_{\ell}) = \hat{T}_{\ell}^{(\ell)} - \alpha_{\ell}(\hat{U}_{\ell}^{(\ell)} - U_{\ell}), \qquad (3.11) $$

with $\hat{U}_{\ell}^{(\ell)}$ an unbiased MC estimator of the control variate statistic $U_k = \tau_k - \tau_{k-1}$, based on the same input sample $\mathcal{X}^{(\ell)}$ as $\hat{T}_{\ell}^{(\ell)}$, again with members of $\mathcal{X}^{(\ell)}$ and $\mathcal{X}^{(\ell')}$ being independent for $\ell \neq \ell'$. Note that, in eq. (3.11), because a single control variate is used at each level, the definition of $\hat{U}_{\ell}^{(\ell)}$ is used in the specific case where $k = \ell$. The more general definition when $k$ and $\ell$ are not necessarily equal will be useful later in section 3.2.3, where we consider multiple control variates per level. In practice, $\hat{U}_{\ell}^{(\ell)}$ may be defined as $\hat{U}_{\ell}^{(\ell)} = \hat{\tau}_{k}^{(\ell)} - \hat{\tau}_{k-1}^{(\ell)}$, where $\hat{\tau}_{k}^{(\ell)}$ is an unbiased estimator of $\tau_k$ using the $n_{\ell}$-sample $\mathcal{X}^{(\ell)}$.

The optimal value $\alpha_{\ell}^*$ for $\alpha_{\ell}$ is obtained individually for each $\ell = 0, \dots, L$ as the optimal (single) CV parameter for $\hat{T}_{\ell}^{\text{CV}}(\alpha_{\ell})$,

$$ \alpha_{\ell}^* = \frac{\mathbb{C}[\hat{T}_{\ell}^{(\ell)}, \hat{U}_{\ell}^{(\ell)}]}{\mathbb{V}[\hat{U}_{\ell}^{(\ell)}]}, \qquad (3.12) $$

and the resulting variance of the control variate estimator of the correction is

$$ \mathbb{V}[\hat{T}_{\ell}^{\text{CV}}(\alpha_{\ell}^*)] = (1 - \rho_{\ell}^2)\mathbb{V}[\hat{T}_{\ell}^{(\ell)}], \quad \text{with } \rho_{\ell} = \frac{\mathbb{C}[\hat{T}_{\ell}^{(\ell)}, \hat{U}_{\ell}^{(\ell)}]}{\mathbb{V}[\hat{T}_{\ell}^{(\ell)}]^{1/2}\mathbb{V}[\hat{U}_{\ell}^{(\ell)}]^{1/2}} \in [-1, 1]. \qquad (3.13) $$

The correction estimators $(\hat{T}_{\ell}^{\text{CV}})_{\ell=0}^L$ being mutually independent, the variance of the optimal MLMC-CV estimator is

$$ \mathbb{V}[\hat{\theta}^{\text{MLMC-CV}}(\alpha_0^*, \dots, \alpha_L^*)] = \sum_{\ell=0}^{L} \mathbb{V}[\hat{T}_{\ell}^{\text{CV}}(\alpha_{\ell}^*)] = \sum_{\ell=0}^{L} (1 - \rho_{\ell}^2) \mathbb{V}[\hat{T}_{\ell}^{(\ell)}] \le \sum_{\ell=0}^{L} \mathbb{V}[\hat{T}_{\ell}^{(\ell)}], \qquad (3.14) $$

indicating that the variance of the MLMC-CV estimator is smaller (as long as $\rho_{\ell}^2 > 0$) than that of the MLMC estimator; see eq. (3.5). We remark that the variance reduction depends on the (squared) correlation between $\hat{T}_{\ell}^{(\ell)}$ and $\hat{U}_{\ell}^{(\ell)}$, which, in turn, typically relates to some measure of similarity between high-fidelity corrections $Y_{\ell} - Y_{\ell-1}$ and the corresponding CV corrections (see appendix A for the expectation and variance estimators in a single-level setting).

In our surrogate-based approach, we may define control variates as $Z_{\ell} = g_{\ell}(\mathbf{X})$, where $g_{\ell}$ is a surrogate of $f_{\ell}$, so that $\tau_{\ell}$ and $\tau_{\ell-1}$ could be estimated using samples of $Z_{\ell}$ and $Z_{\ell-1}$, respectively. It is then crucial to construct the surrogates such that their successive differences $g_{\ell} - g_{\ell-1}$ are good approximations of $f_{\ell} - f_{\ell-1}$, to ensure a high similarity between $Y_{\ell} - Y_{\ell-1}$ and $Z_{\ell} - Z_{\ell-1}$. However,
### 3.2.2 Multilevel Monte Carlo with control variates (MLMC-CV)

The second strategy consists of improving the MLMC estimator eq. (3.2) using the surrogate-based control variates $Z_0, \dots, Z_L$. Specifically, the MLMC-CV estimator improves the MC estimation of each of the correction terms of the MLMC estimator by using a surrogate model of the corresponding correction term as a control variate. Note that the MLMF approach proposed in [18, 19] is based on a similar strategy, using arbitrary low-fidelity models at level in an approximate CV setting. The MLMC-CV estimator reads

$$ \hat{\theta}^{\text{MLMC-CV}}(\alpha_1, \dots, \alpha_L) = \sum_{\ell=0}^{L} \hat{T}_{\ell}^{\text{CV}}(\alpha_{\ell}), \qquad (3.10) $$

where $\hat{T}_{\ell}^{\text{CV}}(\alpha_{\ell})$ is the CV estimator of the multilevel correction term $T_{\ell}$ (cf. eq. (3.11)),

$$ \hat{T}_{\ell}^{\text{CV}}(\alpha_{\ell}) = \hat{T}_{\ell}^{(\ell)} - \alpha_{\ell}(\hat{U}_{\ell}^{(\ell)} - U_{\ell}), \qquad (3.11) $$

with $\hat{U}_{\ell}^{(\ell)}$ an unbiased MC estimator of the control variate statistic $U_k = \tau_k - \tau_{k-1}$, based on the same input sample $\mathcal{X}^{(\ell)}$ as $\hat{T}_{\ell}^{(\ell)}$, again with members of $\mathcal{X}^{(\ell)}$ and $\mathcal{X}^{(\ell')}$ being independent for $\ell \neq \ell'$. Note that, in eq. (3.11), because a single control variate is used at each level, the definition of $\hat{U}_{\ell}^{(\ell)}$ is used in the specific case where $k = \ell$. The more general definition when $k$ and $\ell$ are not necessarily equal will be useful later in section 3.2.3, where we consider multiple control variates per level. In practice, $\hat{U}_{\ell}^{(\ell)}$ may be defined as $\hat{U}_{\ell}^{(\ell)} = \hat{\tau}_{k}^{(\ell)} - \hat{\tau}_{k-1}^{(\ell)}$, where $\hat{\tau}_{k}^{(\ell)}$ is an unbiased estimator of $\tau_k$ using the $n_{\ell}$-sample $\mathcal{X}^{(\ell)}$.

The optimal value $\alpha_{\ell}^*$ for $\alpha_{\ell}$ is obtained individually for each $\ell = 0, \dots, L$ as the optimal (single) CV parameter for $\hat{T}_{\ell}^{\text{CV}}(\alpha_{\ell})$,

$$ \alpha_{\ell}^* = \frac{\mathbb{C}[\hat{T}_{\ell}^{(\ell)}, \hat{U}_{\ell}^{(\ell)}]}{\mathbb{V}[\hat{U}_{\ell}^{(\ell)}]}, \qquad (3.12) $$

and the resulting variance of the control variate estimator of the correction is

$$ \mathbb{V}[\hat{T}_{\ell}^{\text{CV}}(\alpha_{\ell}^*)] = (1 - \rho_{\ell}^2)\mathbb{V}[\hat{T}_{\ell}^{(\ell)}], \quad \text{with } \rho_{\ell} = \frac{\mathbb{C}[\hat{T}_{\ell}^{(\ell)}, \hat{U}_{\ell}^{(\ell)}]}{\mathbb{V}[\hat{T}_{\ell}^{(\ell)}]^{1/2}\mathbb{V}[\hat{U}_{\ell}^{(\ell)}]^{1/2}} \in [-1, 1]. \qquad (3.13) $$

The correction estimators $(\hat{T}_{\ell}^{\text{CV}})_{\ell=0}^L$ being mutually independent, the variance of the optimal MLMC-CV estimator is

$$ \mathbb{V}[\hat{\theta}^{\text{MLMC-CV}}(\alpha_0^*, \dots, \alpha_L^*)] = \sum_{\ell=0}^{L} \mathbb{V}[\hat{T}_{\ell}^{\text{CV}}(\alpha_{\ell}^*)] = \sum_{\ell=0}^{L} (1 - \rho_{\ell}^2) \mathbb{V}[\hat{T}_{\ell}^{(\ell)}] \le \sum_{\ell=0}^{L} \mathbb{V}[\hat{T}_{\ell}^{(\ell)}], \qquad (3.14) $$

indicating that the variance of the MLMC-CV estimator is smaller (as long as $\rho_{\ell}^2 > 0$) than that of the MLMC estimator; see eq. (3.5). We remark that the variance reduction depends on the (squared) correlation between $\hat{T}_{\ell}^{(\ell)}$ and $\hat{U}_{\ell}^{(\ell)}$, which, in turn, typically relates to some measure of similarity between high-fidelity corrections $Y_{\ell} - Y_{\ell-1}$ and the corresponding CV corrections (see appendix A for the expectation and variance estimators in a single-level setting).

In our surrogate-based approach, we may define control variates as $Z_{\ell} = g_{\ell}(\mathbf{X})$, where $g_{\ell}$ is a surrogate of $f_{\ell}$, so that $\tau_{\ell}$ and $\tau_{\ell-1}$ could be estimated using samples of $Z_{\ell}$ and $Z_{\ell-1}$, respectively. It is then crucial to construct the surrogates such that their successive differences $g_{\ell} - g_{\ell-1}$ are good approximations of $f_{\ell} - f_{\ell-1}$, to ensure a high similarity between $Y_{\ell} - Y_{\ell-1}$ and $Z_{\ell} - Z_{\ell-1}$. However,
constructing $g_\ell$ as a surrogate of $f_\ell$ does not give any guarantee on the quality of $g_\ell - g_{\ell-1}$ as a surrogate of $f_\ell - f_{\ell-1}$. Instead, in addition to the surrogates models $g_\ell$ of $f_\ell$, we construct surrogate models $h_\ell$ of the differences $f_\ell - f_{\ell-1}$, and we define auxiliary surrogate models $\tilde{g}_{\ell-1} = g_\ell - h_\ell$, for $\ell = 1, \dots, L$. On each level $\ell$, we then use samples of $Z_\ell = g_\ell(\mathbf{X})$ and $\tilde{Z}_{\ell-1} = \tilde{g}_{\ell-1}(\mathbf{X})$ for the estimation of $\tau_\ell$ and $\tau_{\ell-1}$, respectively. As a result, the variance reduction now depends on the similarity between $Y_\ell - Y_{\ell-1}$ and $W_\ell := Z_\ell - \tilde{Z}_{\ell-1} = g_\ell(\mathbf{X}) - \tilde{g}_{\ell-1}(\mathbf{X}) = h_\ell(\mathbf{X})$, where $h_\ell$ has been constructed to ensure such a similarity. Specifically, for the expectation, the MLMC-CV estimator reads

$$ \hat{E}^{\text{MLMC-CV}}[Y](\alpha_0, \dots, \alpha_L) = \sum_{\ell=0}^{L} \left( \hat{E}^{(\ell)}[Y_\ell] - \hat{E}^{(\ell)}[Y_{\ell-1}] \right) - \alpha_\ell \left( \hat{E}^{(\ell)}[Z_\ell] - \hat{E}^{(\ell)}[\tilde{Z}_{\ell-1}] - (\mu_{Z_\ell} - \mu_{\tilde{Z}_{\ell-1}}) \right), \quad (3.15) $$

with optimal values of $\alpha_\ell$ given by (see appendix A.1, with $M = 1$)

$$ \alpha_{\ell}^{*} = \frac{\mathbb{C}[Y_{\ell} - Y_{\ell-1}, W_{\ell}]}{\mathbb{V}[Z_{\ell} - Z_{\ell-1}]}, \quad \forall \ell = 0, \dots, L, \qquad (3.16) $$

resulting in level-dependent reduction factors $1 - \rho_\ell^2$, where

$$ \rho_{\ell}^{2} = \frac{\mathbb{C}[Y_{\ell} - Y_{\ell-1}, W_{\ell}]^{2}}{\mathbb{V}[Y_{\ell} - Y_{\ell-1}]\mathbb{V}[W_{\ell}]} \qquad (3.17) $$

is the squared correlation coefficient between $Y_\ell - Y_{\ell-1}$ and $W_\ell = h_\ell(\mathbf{X})$. It should be noted that, because of the linearity of the expectation operator and its MC estimator, the use of $\tilde{g}_{\ell-1}$ is superfluous. Indeed, we may directly define the MLMC-CV estimator of the expectation as

$$ \hat{E}^{\text{MLMC-CV}}[Y](\alpha_0, \dots, \alpha_L) = \sum_{\ell=0}^{L} (\hat{E}^{(\ell)}[Y_\ell] - \hat{E}^{(\ell)}[Y_{\ell-1}]) - \alpha_\ell (\hat{E}^{(\ell)}[W_\ell] - \mu_{W_\ell}), \quad (3.18) $$

with $\mu_{W_\ell} = \mathbb{E}[W_\ell]$. For the variance, the MLMC-CV estimator reads

$$ \hat{V}^{\text{MLMC-CV}}[Y](\alpha_0, \dots, \alpha_L) = \sum_{\ell=0}^{L} \left( \hat{V}^{(\ell)}[Y_\ell] - \hat{V}^{(\ell)}[Y_{\ell-1}] \right) - \alpha_\ell \left( \hat{V}^{(\ell)}[Z_\ell] - \hat{V}^{(\ell)}[\tilde{Z}_{\ell-1}] - (\sigma_{Z_\ell}^2 - \sigma_{\tilde{Z}_{\ell-1}}^2) \right), \quad (3.19) $$

with $\sigma_{Z_\ell}^2 = \mathbb{V}[Z_\ell]$ and $\sigma_{\tilde{Z}_{\ell-1}}^2 = \mathbb{V}[\tilde{Z}_{\ell-1}]$. The resulting level-dependent reduction factors are then related to the correlation between $(Y_\ell - Y_{\ell-1} - \mathbb{E}[Y_\ell - Y_{\ell-1}])^2$ and $(h_\ell(\mathbf{X}) - \mathbb{E}[h_\ell(\mathbf{X})])^2$ (see appendix A.2, with $M = 1$). Further details on the construction of $h_\ell$ will be given in section 4.

### 3.2.3 Multilevel Monte Carlo with multilevel control variates (MLMC-MLCV)

We now propose to further improve the MLMC-CV estimator eq. (3.10) by combining MLMC with the MLCV approach described in section 3.2.1, resulting in the MLMC-MLCV method. The approach consists in using, at each level $\ell$, the surrogate-based control variates of all the levels $\ell' = 0, \dots, L$, rather than only using those of level $\ell$, as was previously done with the MLMC-CV approach.
constructing $g_\ell$ as a surrogate of $f_\ell$ does not give any guarantee on the quality of $g_\ell - g_{\ell-1}$ as a surrogate of $f_\ell - f_{\ell-1}$. Instead, in addition to the surrogates models $g_\ell$ of $f_\ell$, we construct surrogate models $h_\ell$ of the differences $f_\ell - f_{\ell-1}$, and we define auxiliary surrogate models $\tilde{g}_{\ell-1} = g_\ell - h_\ell$, for $\ell = 1, \dots, L$. On each level $\ell$, we then use samples of $Z_\ell = g_\ell(\mathbf{X})$ and $\tilde{Z}_{\ell-1} = \tilde{g}_{\ell-1}(\mathbf{X})$ for the estimation of $\tau_\ell$ and $\tau_{\ell-1}$, respectively. As a result, the variance reduction now depends on the similarity between $Y_\ell - Y_{\ell-1}$ and $W_\ell := Z_\ell - \tilde{Z}_{\ell-1} = g_\ell(\mathbf{X}) - \tilde{g}_{\ell-1}(\mathbf{X}) = h_\ell(\mathbf{X})$, where $h_\ell$ has been constructed to ensure such a similarity. Specifically, for the expectation, the MLMC-CV estimator reads

$$ \hat{E}^{\text{MLMC-CV}}[Y](\alpha_0, \dots, \alpha_L) = \sum_{\ell=0}^{L} \left( \hat{E}^{(\ell)}[Y_\ell] - \hat{E}^{(\ell)}[Y_{\ell-1}] \right) - \alpha_\ell \left( \hat{E}^{(\ell)}[Z_\ell] - \hat{E}^{(\ell)}[\tilde{Z}_{\ell-1}] - (\mu_{Z_\ell} - \mu_{\tilde{Z}_{\ell-1}}) \right), \quad (3.15) $$

with optimal values of $\alpha_\ell$ given by (see appendix A.1, with $M = 1$)

$$ \alpha_{\ell}^{*} = \frac{\mathbb{C}[Y_{\ell} - Y_{\ell-1}, W_{\ell}]}{\mathbb{V}[Z_{\ell} - Z_{\ell-1}]}, \quad \forall \ell = 0, \dots, L, \qquad (3.16) $$

resulting in level-dependent reduction factors $1 - \rho_\ell^2$, where

$$ \rho_{\ell}^{2} = \frac{\mathbb{C}[Y_{\ell} - Y_{\ell-1}, W_{\ell}]^{2}}{\mathbb{V}[Y_{\ell} - Y_{\ell-1}]\mathbb{V}[W_{\ell}]} \qquad (3.17) $$

is the squared correlation coefficient between $Y_\ell - Y_{\ell-1}$ and $W_\ell = h_\ell(\mathbf{X})$. It should be noted that, because of the linearity of the expectation operator and its MC estimator, the use of $\tilde{g}_{\ell-1}$ is superfluous. Indeed, we may directly define the MLMC-CV estimator of the expectation as

$$ \hat{E}^{\text{MLMC-CV}}[Y](\alpha_0, \dots, \alpha_L) = \sum_{\ell=0}^{L} (\hat{E}^{(\ell)}[Y_\ell] - \hat{E}^{(\ell)}[Y_{\ell-1}]) - \alpha_\ell (\hat{E}^{(\ell)}[W_\ell] - \mu_{W_\ell}), \quad (3.18) $$

with $\mu_{W_\ell} = \mathbb{E}[W_\ell]$. For the variance, the MLMC-CV estimator reads

$$ \hat{V}^{\text{MLMC-CV}}[Y](\alpha_0, \dots, \alpha_L) = \sum_{\ell=0}^{L} \left( \hat{V}^{(\ell)}[Y_\ell] - \hat{V}^{(\ell)}[Y_{\ell-1}] \right) - \alpha_\ell \left( \hat{V}^{(\ell)}[Z_\ell] - \hat{V}^{(\ell)}[\tilde{Z}_{\ell-1}] - (\sigma_{Z_\ell}^2 - \sigma_{\tilde{Z}_{\ell-1}}^2) \right), \quad (3.19) $$

with $\sigma_{Z_\ell}^2 = \mathbb{V}[Z_\ell]$ and $\sigma_{\tilde{Z}_{\ell-1}}^2 = \mathbb{V}[\tilde{Z}_{\ell-1}]$. The resulting level-dependent reduction factors are then related to the correlation between $(Y_\ell - Y_{\ell-1} - \mathbb{E}[Y_\ell - Y_{\ell-1}])^2$ and $(h_\ell(\mathbf{X}) - \mathbb{E}[h_\ell(\mathbf{X})])^2$ (see appendix A.2, with $M = 1$). Further details on the construction of $h_\ell$ will be given in section 4.

### 3.2.3 Multilevel Monte Carlo with multilevel control variates (MLMC-MLCV)

We now propose to further improve the MLMC-CV estimator eq. (3.10) by combining MLMC with the MLCV approach described in section 3.2.1, resulting in the MLMC-MLCV method. The approach consists in using, at each level $\ell$, the surrogate-based control variates of all the levels $\ell' = 0, \dots, L$, rather than only using those of level $\ell$, as was previously done with the MLMC-CV approach.
The MLMC-MLCV estimator then reads

$$ \hat{\theta}^{\text{MLMC-MLCV}}(\boldsymbol{\alpha}_0, \dots, \boldsymbol{\alpha}_L) = \sum_{\ell=0}^{L} \hat{T}_{\ell}^{\text{MLCV}}(\boldsymbol{\alpha}_{\ell}), \qquad (3.20) $$

where $\boldsymbol{\alpha}_\ell$ denotes the CV parameter at level $\ell$, and $\hat{T}_\ell^{\text{MLCV}}(\boldsymbol{\alpha}_\ell)$ is the MLCV estimator of $T_\ell$,

$$ \hat{T}_{\ell}^{\text{MLCV}}(\boldsymbol{\alpha}_{\ell}) = \hat{T}_{\ell}^{(\ell)} - \boldsymbol{\alpha}_{\ell}^{\text{T}}(\hat{\mathbf{U}}_{\ell}^{(\ell)} - \mathbf{U}_{\ell}), \qquad (3.21) $$

with

$$ \mathbf{U}_0 = (\tau_k)_{k=0}^L \qquad \hat{\mathbf{U}}_0^{(0)} = (\hat{\tau}_k^{(0)})_{k=0}^L \qquad (3.22) $$

$$ \mathbf{U}_{\ell} = (U_k)_{k=1}^L, \quad \text{for } \ell > 0, \qquad \hat{\mathbf{U}}_{\ell}^{(\ell)} = (\hat{U}_k^{(\ell)})_{k=1}^L, \quad \text{for } \ell > 0, \qquad (3.23) $$

and with $U_k$ and $\hat{U}_k^{(\ell)}$ defined as in section 3.2.2.

Because each term eq. (3.21) is an unbiased (multiple) CV estimator of $T_\ell$, the resulting estimator eq. (3.20) is also unbiased, and the optimal (variance minimizing) value $\boldsymbol{\alpha}_\ell^*$ of $\boldsymbol{\alpha}_\ell$ is given individually for each $\ell = 0, \dots, L$ as the optimal (multiple) CV parameter for $\hat{T}_\ell^{(\ell)}(\boldsymbol{\alpha}_\ell)$,

$$ \boldsymbol{\alpha}_{\ell}^{*} = \mathbb{C}[\hat{\mathbf{U}}_{\ell}^{(\ell)}]^{-1} \mathbb{C}[\hat{\mathbf{U}}_{\ell}^{(\ell)}, \hat{T}_{\ell}^{(\ell)}]. \qquad (3.24) $$

The resulting variance is given by

$$ \mathbb{V}[\hat{\theta}^{\text{MLMC-MLCV}}(\boldsymbol{\alpha}_0^*, \dots, \boldsymbol{\alpha}_L^*)] = \sum_{\ell=0}^{L} \mathbb{V}[\hat{T}_{\ell}^{\text{MLCV}}(\boldsymbol{\alpha}_{\ell}^*)] = \sum_{\ell=0}^{L} (1 - R_{\ell}^2) \mathbb{V}[\hat{T}_{\ell}^{(\ell)}] \qquad (3.25) $$

with

$$ R_{\ell}^{2} = \mathbb{V}[\hat{T}_{\ell}^{(\ell)}]^{-1} \mathbb{C}[\hat{\mathbf{U}}^{(\ell)}, \hat{T}_{\ell}^{(\ell)}]^{\top} \boldsymbol{\alpha}_{\ell}^{*} \in [0, 1]. \qquad (3.26) $$

Again, owing to the fact that $R_\ell^2 \le 1$, the variance of the MLMC-MLCV estimator is always less than or equal to the variance of the MLMC estimator given by $\mathbb{V}[\hat{\theta}_L^{\text{MLMC}}] = \sum_{\ell=0}^L \mathbb{V}[\hat{T}_\ell^{(\ell)}]$ (see eq. (3.5)).

In our surrogate-based approach, $L + 1$ surrogate-based control variates can be used for the coarsest level $\ell = 0$, namely $g_0, \dots, g_L$, so that $\boldsymbol{\alpha}_0 \in \mathbb{R}^{L+1}$. At correction levels $\ell > 0$, we can use $L$ control variates based on $g_1, \dots, g_L$ and $\tilde{g}_0, \dots, \tilde{g}_{L-1}$, as described in section 3.2.2, so that $\boldsymbol{\alpha}_\ell \in \mathbb{R}^L$, for $\ell > 0$.

The MLMC-MLCV estimator of the expectation $\mathbb{E}[f(\mathbf{X})]$ is defined by eq. (3.20), with

$$ \hat{T}_0^{\text{MLCV}}(\boldsymbol{\alpha}_0) = \hat{E}^{(0)}[Y_0] - \boldsymbol{\alpha}_0^{\text{T}}(\hat{E}^{(0)}[\mathbf{Z}] - \boldsymbol{\mu}_{\mathbf{Z}}), \qquad (3.27) $$

$$ \hat{T}_{\ell}^{\text{MLCV}}(\boldsymbol{\alpha}_{\ell}) = \hat{E}^{(\ell)}[Y_{\ell}] - \hat{E}^{(\ell)}[Y_{\ell-1}] - \boldsymbol{\alpha}_{\ell}^{\text{T}}(\hat{E}^{(\ell)}[\mathbf{W}] - \boldsymbol{\mu}_{\mathbf{W}}), \quad \text{for } \ell > 0, \qquad (3.28) $$

where

$$ \mathbf{Z} = (Z_0, \dots, Z_L) = (g_0(\mathbf{X}), \dots, g_L(\mathbf{X})), \qquad \boldsymbol{\mu}_{\mathbf{Z}} = \mathbb{E}[\mathbf{Z}], \qquad (3.29) $$

$$ \mathbf{W} = (W_1, \dots, W_L) = (h_1(\mathbf{X}), \dots, h_L(\mathbf{X})), \qquad \boldsymbol{\mu}_{\mathbf{W}} = \mathbb{E}[\mathbf{W}]. \qquad (3.30) $$

The optimal values $\boldsymbol{\alpha}_\ell^*$ of the CV parameters are given by

$$ \boldsymbol{\alpha}_0^* = \boldsymbol{\Sigma}_0^{-1} \mathbf{c}_0, \qquad \boldsymbol{\Sigma}_0 = \mathbb{C}[\mathbf{Z}], \qquad \mathbf{c}_0 = \mathbb{C}[\mathbf{Z}, Y_0], \qquad (3.31) $$

$$ \boldsymbol{\alpha}_{\ell}^{*} = \boldsymbol{\Sigma}_{\ell}^{-1} \mathbf{c}_{\ell}, \qquad \boldsymbol{\Sigma}_{\ell} = \mathbb{C}[\mathbf{W}], \qquad \mathbf{c}_{\ell} = \mathbb{C}[\mathbf{W}, Y_{\ell} - Y_{\ell-1}], \quad \text{for } \ell > 0, \qquad (3.32) $$

resulting in $R_\ell^2 = \mathbb{V}[Y_\ell - Y_{\ell-1}]^{-1} \mathbf{c}_\ell^\top \boldsymbol{\Sigma}_\ell \mathbf{c}_\ell$, for $\ell = 0, \dots, L$. Note that $\boldsymbol{\Sigma}_\ell = \mathbb{C}[\mathbf{W}]$ is the same for all $\ell > 0$. The optimal MLMC-MLCV variance estimator is derived in appendix C.
The MLMC-MLCV estimator then reads

$$ \hat{\theta}^{\text{MLMC-MLCV}}(\boldsymbol{\alpha}_0, \dots, \boldsymbol{\alpha}_L) = \sum_{\ell=0}^{L} \hat{T}_{\ell}^{\text{MLCV}}(\boldsymbol{\alpha}_{\ell}), \qquad (3.20) $$

where $\boldsymbol{\alpha}_\ell$ denotes the CV parameter at level $\ell$, and $\hat{T}_\ell^{\text{MLCV}}(\boldsymbol{\alpha}_\ell)$ is the MLCV estimator of $T_\ell$,

$$ \hat{T}_{\ell}^{\text{MLCV}}(\boldsymbol{\alpha}_{\ell}) = \hat{T}_{\ell}^{(\ell)} - \boldsymbol{\alpha}_{\ell}^{\text{T}}(\hat{\mathbf{U}}_{\ell}^{(\ell)} - \mathbf{U}_{\ell}), \qquad (3.21) $$

with

$$ \mathbf{U}_0 = (\tau_k)_{k=0}^L \qquad \hat{\mathbf{U}}_0^{(0)} = (\hat{\tau}_k^{(0)})_{k=0}^L \qquad (3.22) $$

$$ \mathbf{U}_{\ell} = (U_k)_{k=1}^L, \quad \text{for } \ell > 0, \qquad \hat{\mathbf{U}}_{\ell}^{(\ell)} = (\hat{U}_k^{(\ell)})_{k=1}^L, \quad \text{for } \ell > 0, \qquad (3.23) $$

and with $U_k$ and $\hat{U}_k^{(\ell)}$ defined as in section 3.2.2.

Because each term eq. (3.21) is an unbiased (multiple) CV estimator of $T_\ell$, the resulting estimator eq. (3.20) is also unbiased, and the optimal (variance minimizing) value $\boldsymbol{\alpha}_\ell^*$ of $\boldsymbol{\alpha}_\ell$ is given individually for each $\ell = 0, \dots, L$ as the optimal (multiple) CV parameter for $\hat{T}_\ell^{(\ell)}(\boldsymbol{\alpha}_\ell)$,

$$ \boldsymbol{\alpha}_{\ell}^{*} = \mathbb{C}[\hat{\mathbf{U}}_{\ell}^{(\ell)}]^{-1} \mathbb{C}[\hat{\mathbf{U}}_{\ell}^{(\ell)}, \hat{T}_{\ell}^{(\ell)}]. \qquad (3.24) $$

The resulting variance is given by

$$ \mathbb{V}[\hat{\theta}^{\text{MLMC-MLCV}}(\boldsymbol{\alpha}_0^*, \dots, \boldsymbol{\alpha}_L^*)] = \sum_{\ell=0}^{L} \mathbb{V}[\hat{T}_{\ell}^{\text{MLCV}}(\boldsymbol{\alpha}_{\ell}^*)] = \sum_{\ell=0}^{L} (1 - R_{\ell}^2) \mathbb{V}[\hat{T}_{\ell}^{(\ell)}] \qquad (3.25) $$

with

$$ R_{\ell}^{2} = \mathbb{V}[\hat{T}_{\ell}^{(\ell)}]^{-1} \mathbb{C}[\hat{\mathbf{U}}^{(\ell)}, \hat{T}_{\ell}^{(\ell)}]^{\top} \boldsymbol{\alpha}_{\ell}^{*} \in [0, 1]. \qquad (3.26) $$

Again, owing to the fact that $R_\ell^2 \le 1$, the variance of the MLMC-MLCV estimator is always less than or equal to the variance of the MLMC estimator given by $\mathbb{V}[\hat{\theta}_L^{\text{MLMC}}] = \sum_{\ell=0}^L \mathbb{V}[\hat{T}_\ell^{(\ell)}]$ (see eq. (3.5)).

In our surrogate-based approach, $L + 1$ surrogate-based control variates can be used for the coarsest level $\ell = 0$, namely $g_0, \dots, g_L$, so that $\boldsymbol{\alpha}_0 \in \mathbb{R}^{L+1}$. At correction levels $\ell > 0$, we can use $L$ control variates based on $g_1, \dots, g_L$ and $\tilde{g}_0, \dots, \tilde{g}_{L-1}$, as described in section 3.2.2, so that $\boldsymbol{\alpha}_\ell \in \mathbb{R}^L$, for $\ell > 0$.

The MLMC-MLCV estimator of the expectation $\mathbb{E}[f(\mathbf{X})]$ is defined by eq. (3.20), with

$$ \hat{T}_{0}^{\text{MLCV}}(\boldsymbol{\alpha}_{0}) = \hat{E}^{(0)}[Y_{0}] - \boldsymbol{\alpha}_{0}^{\text{T}}(\hat{E}^{(0)}[\mathbf{Z}] - \boldsymbol{\mu}_{\mathbf{Z}}), \qquad (3.27) $$

$$ \hat{T}_{\ell}^{\text{MLCV}}(\boldsymbol{\alpha}_{\ell}) = \hat{E}^{(\ell)}[Y_{\ell}] - \hat{E}^{(\ell)}[Y_{\ell-1}] - \boldsymbol{\alpha}_{\ell}^{\text{T}}(\hat{E}^{(\ell)}[\mathbf{W}] - \boldsymbol{\mu}_{\mathbf{W}}), \quad \text{for } \ell > 0, \qquad (3.28) $$

where

$$ \mathbf{Z} = (Z_0, \dots, Z_L) = (g_0(\mathbf{X}), \dots, g_L(\mathbf{X})), \qquad \boldsymbol{\mu}_{\mathbf{Z}} = \mathbb{E}[\mathbf{Z}], \qquad (3.29) $$

$$ \mathbf{W} = (W_1, \dots, W_L) = (h_1(\mathbf{X}), \dots, h_L(\mathbf{X})), \qquad \boldsymbol{\mu}_{\mathbf{W}} = \mathbb{E}[\mathbf{W}]. \qquad (3.30) $$

The optimal values $\boldsymbol{\alpha}_\ell^*$ of the CV parameters are given by

$$ \boldsymbol{\alpha}_{0}^{*} = \boldsymbol{\Sigma}_{0}^{-1} \mathbf{c}_{0}, \qquad \boldsymbol{\Sigma}_{0} = \mathbb{C}[\mathbf{Z}], \qquad \mathbf{c}_{0} = \mathbb{C}[\mathbf{Z}, Y_{0}], \qquad (3.31) $$

$$ \boldsymbol{\alpha}_{\ell}^{*} = \boldsymbol{\Sigma}_{\ell}^{-1} \mathbf{c}_{\ell}, \qquad \boldsymbol{\Sigma}_{\ell} = \mathbb{C}[\mathbf{W}], \qquad \mathbf{c}_{\ell} = \mathbb{C}[\mathbf{W}, Y_{\ell} - Y_{\ell-1}], \quad \text{for } \ell > 0, \qquad (3.32) $$

resulting in $R_\ell^2 = \mathbb{V}[Y_\ell - Y_{\ell-1}]^{-1} \mathbf{c}_\ell^\top \boldsymbol{\Sigma}_\ell \mathbf{c}_\ell$, for $\ell = 0, \dots, L$. Note that $\boldsymbol{\Sigma}_\ell = \mathbb{C}[\mathbf{W}]$ is the same for all $\ell > 0$. The optimal MLMC-MLCV variance estimator is derived in appendix C.
**Remark 3.1.** In practice, $\Sigma$ and $c_\ell$ may be estimated using either a pilot sample or the same sample as for the estimation of $\hat{U}_\ell^{(\ell)}$ and $\hat{T}_\ell^{(\ell)}$. Alternatively, in the specific context of PC-based control variates, for the estimation of the expectation, a closed-form expression for $\Sigma$ can be obtained. Letting $Z_\ell = g_\ell(\mathbf{X}) = \sum_{k=0}^{P_\ell^g} g_{\ell,k} \Psi_k(\mathbf{X})$ and $W_\ell = h_\ell(\mathbf{X}) = \sum_{k=0}^{P_\ell^h} h_{\ell,k} \Psi_k(\mathbf{X})$, we have

$$ \forall m, m' = 0, \dots, L, \quad [\Sigma_0]_{m,m'} = \mathbb{C}[Z_m, Z_{m'}] = \sum_{k=1}^{\min(P_\ell^m, P_\ell^{m'})} \mathfrak{g}_{m,k} \mathfrak{g}_{m',k}, \quad (3.33) $$

$$ \forall m, m' = 1, \dots, L, \quad [\Sigma_\ell]_{m,m'} = \mathbb{C}[W_m, W_{m'}] = \sum_{k=1}^{\min(P_h^m, P_h^{m'})} \mathfrak{h}_{m,k} \mathfrak{h}_{m',k}, \quad \text{for } \ell > 0. \quad (3.34) $$

## 3.3 Practical details

A summary of the methods is presented in table 1, including two minor variations of the MLMC-CV and MLMC-MLCV estimators introduced previously, which consist of using fewer surrogate models for the multilevel estimation. We remark that all the MLMC-like estimators (including the MLMC estimator), hereafter abbreviated MLMC-, have variance $\sum_{\ell=0}^{L} (1 - R_\ell^2) \mathbb{V}[\hat{T}_\ell^{(\ell)}]$, with

* $R_\ell^2 = 0$ for plain MLMC;

* $R_\ell^2 = \rho_\ell^2$ as defined in eq. (3.13) for MLMC-CV; and

* $R_\ell^2$ defined by eq. (3.26) for MLMC-MLCV.

For all these methods, we will further assume that $\mathbb{V}[\hat{T}_\ell^{(\ell)}] \lesssim n_\ell^{-1} \mathcal{V}_\ell$ (see section 3.1), which implies that $\mathbb{V}[\hat{T}_\ell^{\text{MLMC-*}}(\alpha_\ell^*)] = (1 - R_\ell^2)\mathbb{V}[\hat{T}_\ell^{(\ell)}] \lesssim n_\ell^{-1} \mathcal{V}_\ell^{\text{CV}}$, with $\mathcal{V}_\ell^{\text{CV}} := (1 - R_\ell^2)\mathcal{V}_\ell$.

In the surrogate-based variants of MLMC, the cost of evaluating the surrogate models is assumed to be negligible compared to the costs of evaluating the simulators $f_0, \dots, f_L$. Therefore, the total computational cost of the MLMC-* estimator reduces to the cost of the MLMC estimator, $\text{cost}(\hat{\theta}^{\text{MLMC-MLCV}}) = \text{cost}(\hat{\theta}^{\text{MLMC}})$, given by eq. (3.6). Similarly to the case of the MLMC estimator (see, e.g., [36]), the optimal sample sizes $(n_\ell^*)_{\ell=0}^L$ such that $\sum_{\ell=0}^L \mathcal{V}_\ell^{\text{CV}}/n_\ell^*$ is minimal under a constrained computational budget of $\mathcal{C}$ are given by

$$ n_\ell^* = \frac{C}{S_L} \sqrt{\frac{\mathcal{V}_\ell^{\text{CV}}}{C_\ell + C_{\ell-1}}}, \quad \text{with } S_\ell := \sum_{\ell'=0}^{\ell} \sqrt{(C_{\ell'} + C_{\ell'-1})\mathcal{V}_{\ell'}^{\text{CV}}}, \quad (3.35) $$

so that $\sum_{\ell=0}^L \mathcal{V}_\ell^{\text{CV}}/n_\ell^* = S_L^2/\mathcal{C}$. As a consequence,

$$ \mathbb{V}[\hat{\theta}_L^{\text{MLMC-MLCV}}(\alpha_\ell^*)] \lesssim \frac{S_L^2}{C}, \quad (3.36) $$

with an equality between the left- and right-hand sides for the expectation estimators.

In practice, $\mathcal{V}_\ell^{\text{CV}}$ is not known and must be estimated for each level. In this work, we consider the sequential algorithm proposed in [36, Algorithm 2] for the MLMC. The algorithm starts from an initial, small number of samples $n_\ell^{\text{init}}$ on each level. Then, it selects the optimal level on which to increase the sample size by an inflation factor $r_\ell > 1$, i.e. the level $\ell^*$ on which the reduction in total variance relative to the additional computational effort achieved by inflating the sample size by $r_\ell^*$ is maximal.
**Remark 3.1.** In practice, $\Sigma$ and $c_\ell$ may be estimated using either a pilot sample or the same sample as for the estimation of $\hat{\mathbf{U}}_\ell^{(\ell)}$ and $\hat{T}_\ell^{(\ell)}$. Alternatively, in the specific context of PC-based control variates, for the estimation of the expectation, a closed-form expression for $\Sigma$ can be obtained. Letting $Z_\ell = g_\ell(\mathbf{X}) = \sum_{k=0}^{P_\ell^g} \mathfrak{g}_{\ell,k} \Psi_k(\mathbf{X})$ and $W_\ell = h_\ell(\mathbf{X}) = \sum_{k=0}^{P_\ell^h} \mathfrak{h}_{\ell,k} \Psi_k(\mathbf{X})$, we have

$$ \forall m, m' = 0, \dots, L, \quad [\Sigma_0]_{m,m'} = \mathbb{C}[Z_m, Z_{m'}] = \sum_{k=1}^{\min(P_h^m, P_h^{m'})} \mathfrak{g}_{m,k} \mathfrak{g}_{m',k}, \quad (3.33) $$

$$ \forall m, m' = 1, \dots, L, \quad [\Sigma_\ell]_{m,m'} = \mathbb{C}[W_m, W_{m'}] = \sum_{k=1}^{\min(P_h^m, P_h^{m'})} \mathfrak{h}_{m,k} \mathfrak{h}_{m',k}, \quad \text{for } \ell > 0. \quad (3.34) $$

## 3.3 Practical details

A summary of the methods is presented in table 1, including two minor variations of the MLMC-CV and MLMC-MLCV estimators introduced previously, which consist of using fewer surrogate models for the multilevel estimation. We remark that all the MLMC-like estimators (including the MLMC estimator), hereafter abbreviated MLMC-, have variance $\sum_{\ell=0}^L (1 - R_\ell^2) \mathbb{V}[\hat{T}_\ell^{(\ell)}]$, with

* $R_\ell^2 = 0$ for plain MLMC;

* $R_\ell^2 = \rho_\ell^2$ as defined in eq. (3.13) for MLMC-CV; and

* $R_\ell^2$ defined by eq. (3.26) for MLMC-MLCV.

For all these methods, we will further assume that $\mathbb{V}[\hat{T}_\ell^{(\ell)}] \lesssim n_\ell^{-1} \mathcal{V}_\ell$ (see section 3.1), which implies that $\mathbb{V}[\hat{T}_\ell^{\text{MLMC-}(\alpha_\ell^*)}] = (1 - R_\ell^2)\mathbb{V}[\hat{T}_\ell^{(\ell)}] \lesssim n_\ell^{-1} \mathcal{V}_\ell^{\text{CV}}$, with $\mathcal{V}_\ell^{\text{CV}} := (1 - R_\ell^2)\mathcal{V}_\ell$.

In the surrogate-based variants of MLMC, the cost of evaluating the surrogate models is assumed to be negligible compared to the costs of evaluating the simulators $f_0, \dots, f_L$. Therefore, the total computational cost of the MLMC-* estimator reduces to the cost of the MLMC estimator, $\text{cost}(\hat{\theta}^{\text{MLMC-MLCV}}) = \text{cost}(\hat{\theta}^{\text{MLMC}})$, given by eq. (3.6). Similarly to the case of the MLMC estimator (see, e.g., [36]), the optimal sample sizes $(n_\ell^*)_{\ell=0}^L$ such that $\sum_{\ell=0}^L \mathcal{V}_\ell^{\text{CV}}/n_\ell^*$ is minimal under a constrained computational budget of $C$ are given by

$$ n_\ell^* = \frac{C}{S_L} \sqrt{\frac{\mathcal{V}_\ell^{\text{CV}}}{C_\ell + C_{\ell-1}}}, \quad \text{with } S_\ell := \sum_{\ell'=0}^{\ell} \sqrt{(C_{\ell'} + C_{\ell'-1})\mathcal{V}_{\ell'}^{\text{CV}}}, \quad (3.35) $$

so that $\sum_{\ell=0}^L \mathcal{V}_\ell^{\text{CV}}/n_\ell^* = S_L^2/C$. As a consequence,

$$ \mathbb{V}[\hat{\theta}_L^{\text{MLMC-MLCV}}(\alpha_\ell^*)] \lesssim \frac{S_L^2}{C}, \quad (3.36) $$

with an equality between the left- and right-hand sides for the expectation estimators.

In practice, $\mathcal{V}_\ell^{\text{CV}}$ is not known and must be estimated for each level. In this work, we consider the sequential algorithm proposed in [36, Algorithm 2] for the MLMC. The algorithm starts from an initial, small number of samples $n_\ell^{\text{init}}$ on each level. Then, it selects the optimal level on which to increase the sample size by an inflation factor $r_\ell > 1$, i.e. the level $\ell^*$ on which the reduction in total variance relative to the additional computational effort achieved by inflating the sample size by $r_\ell^*$ is maximal.
<table><caption>Table 1: Summary of the methods. The first three are state-of-the art methods, the next three are the novel multilevel methods proposed in this paper, and the last two are variants of the proposed methods.</caption><thead><tr><th>Method</th><th>Form of the estimator</th><th>Eq.</th></tr></thead><tbody><tr><td>Monte Carlo (MC).</td><td>θ̂</td><td></td></tr><tr><td>Control Variates (CV) [29, 30, 31, 38, 33].</td><td>θ̂ = &sum;<sub>m=1</sub><sup>M</sup> α<sub>m</sub>(τ̂<sub>m</sub> - τ<sub>m</sub>).</td><td>2.4</td></tr><tr><td>Multilevel Monte Carlo (MLMC) [21, 22].</td><td>θ̂<sub>0</sub><sup>(0)</sup> + &sum;<sub>ℓ=1</sub><sup>L</sup> T̂<sub>ℓ</sub><sup>(ℓ)</sup>.</td><td>3.2</td></tr><tr><td>Multilevel Control Variates (MLCV). CV method where the CVs are based on surrogate models of simulators of different levels of fidelity.</td><td>θ̂ = &sum;<sub>ℓ=0</sub><sup>L</sup> α<sub>ℓ</sub>(τ̂<sub>ℓ</sub> - τ<sub>ℓ</sub>).</td><td>2.4</td></tr><tr><td>MLMC-CV. MLMC with one CV at each correction level based on surrogate models of the simulators on the corresponding level. Corresponds to a surrogate-based MLMF [18, 19] where the exact statistics of the CVs are known.</td><td>θ̂<sub>0</sub><sup>(0)</sup> - α<sub>0</sub><sup>(0)</sup> - τ<sub>0</sub> + &sum;<sub>ℓ=1</sub><sup>L</sup> (T̂<sub>ℓ</sub><sup>(ℓ)</sup> - α<sub>ℓ</sub>(Û<sub>ℓ</sub><sup>(ℓ)</sup> - U<sub>ℓ</sub>)).</td><td>3.10, 3.11</td></tr><tr><td>MLMC-MLCV. MLMC-CV with the CVs based on the surrogate models g<sub>0</sub>,..., g<sub>L</sub> on level 0 and on h<sub>1</sub>,..., h<sub>L</sub> on levels ℓ > 0.</td><td>θ̂<sub>0</sub><sup>(0)</sup> - α<sub>0</sub><sup>(0)</sup> + &sum;<sub>ℓ=1</sub><sup>L</sup> T̂<sub>ℓ</sub><sup>(ℓ)</sup>.</td><td>3.20, 3.21</td></tr><tr><td>MLMC-CV[0]. MLMC-CV using only one CV based on the surrogate g<sub>0</sub> on level 0 and no CVs on levels ℓ > 0</td><td>θ̂<sub>0</sub><sup>(0)</sup> - α<sub>0</sub><sup>(0)</sup> + &sum;<sub>ℓ=1</sub><sup>L</sup> T̂<sub>ℓ</sub><sup>(ℓ)</sup>.</td><td></td></tr><tr><td>MLMC-MLCV[0]. MLMC-MLCV using only CVs based on the surrogates g<sub>0</sub> and g<sub>1</sub> on level 0, and h<sub>1</sub> on levels ℓ > 0.</td><td>θ̂<sub>0</sub><sup>(0)</sup> - α<sub>0</sub><sup>(0)</sup> + &sum;<sub>ℓ=1</sub><sup>L</sup> (T̂<sub>ℓ</sub><sup>(ℓ)</sup> - α<sub>ℓ</sub>(Û<sub>1</sub><sup>(ℓ)</sup> - U<sub>1</sub>)).</td><td></td></tr></tbody></table>
<table><caption>Table 1: Summary of the methods. The first three are state-of-the art methods, the next three are the novel multilevel methods proposed in this paper, and the last two are variants of the proposed methods.</caption><thead><tr><th>Method</th><th>Form of the estimator</th><th>Eq.</th></tr></thead><tbody><tr><td>Monte Carlo (MC).</td><td>θ̂</td><td></td></tr><tr><td>Control Variates (CV) [29, 30, 31, 38, 33].</td><td>θ̂ = &sum;<sub>m=1</sub><sup>M</sup> α<sub>m</sub>(τ̂<sub>m</sub> - τ<sub>m</sub>).</td><td>2.4</td></tr><tr><td>Multilevel Monte Carlo (MLMC) [21, 22].</td><td>θ̂<sub>0</sub><sup>(0)</sup> + &sum;<sub>ℓ=1</sub><sup>L</sup> T̂<sub>ℓ</sub><sup>(ℓ)</sup>.</td><td>3.2</td></tr><tr><td>Multilevel Control Variates (MLCV). CV method where the CVs are based on surrogate models of simulators of different levels of fidelity.</td><td>θ̂ = &sum;<sub>ℓ=0</sub><sup>L</sup> α<sub>ℓ</sub>(τ̂<sub>ℓ</sub> - τ<sub>ℓ</sub>).</td><td>2.4</td></tr><tr><td>MLMC-CV. MLMC with one CV at each correction level based on surrogate models of the simulators on the corresponding level. Corresponds to a surrogate-based MLMF [18, 19] where the exact statistics of the CVs are known.</td><td>θ̂<sub>0</sub><sup>(0)</sup> - α<sub>0</sub><sup>(0)</sup> - τ<sub>0</sub> + &sum;<sub>ℓ=1</sub><sup>L</sup> (T̂<sub>ℓ</sub><sup>(ℓ)</sup> - α<sub>ℓ</sub>(Û<sub>ℓ</sub><sup>(ℓ)</sup> - U<sub>ℓ</sub>)).</td><td>3.10, 3.11</td></tr><tr><td>MLMC-MLCV. MLMC-CV with the CVs based on the surrogate models g<sub>0</sub>,..., g<sub>L</sub> on level 0 and on h<sub>1</sub>,..., h<sub>L</sub> on levels ℓ > 0.</td><td>θ̂<sub>0</sub><sup>(0)</sup> - α<sub>0</sub><sup>(0)</sup> + &sum;<sub>ℓ=1</sub><sup>L</sup> T̂<sub>ℓ</sub><sup>(ℓ)</sup>.</td><td>3.20, 3.21</td></tr><tr><td>MLMC-CV[0]. MLMC-CV using only one CV based on the surrogate g<sub>0</sub> on level 0 and no CVs on levels ℓ > 0</td><td>θ̂<sub>0</sub><sup>(0)</sup> - α<sub>0</sub><sup>(0)</sup> + &sum;<sub>ℓ=1</sub><sup>L</sup> T̂<sub>ℓ</sub><sup>(ℓ)</sup>.</td><td></td></tr><tr><td>MLMC-MLCV[0]. MLMC-MLCV using only CVs based on the surrogates g<sub>0</sub> and g<sub>1</sub> on level 0, and h<sub>1</sub> on levels ℓ > 0.</td><td>θ̂<sub>0</sub><sup>(0)</sup> - α<sub>0</sub><sup>(0)</sup> + &sum;<sub>ℓ=1</sub><sup>L</sup> (T̂<sub>ℓ</sub><sup>(ℓ)</sup> - α<sub>ℓ</sub>(Û<sub>1</sub><sup>(ℓ)</sup> - U<sub>1</sub>)).</td><td></td></tr></tbody></table>
**Algorithm 1** Simplified MLMC-* algorithm inspired by [36].

**Require:** $n_{\ell}^{\text{init}} > 1$, $r_{\ell} > 1$, surrogate models (depending on the method), and budget $C$.

1. Set consumed budget to $\tilde{C} = 0$ and $\delta n_{\ell} = n_{\ell}^{\text{init}}$ samples on levels $\ell \le L$;

2. **while** $\tilde{C} \le C$ **do**

3. compute $\delta n_{\ell}$ samples on each level by evaluating $f_{\ell}$ and the appropriate surrogates;

4. update sample size on each level: $n_{\ell} \leftarrow n_{\ell} + \delta n_{\ell}$;

5. update consumed budget: $\tilde{C} \leftarrow \tilde{C} + \sum_{\ell=0}^{L} \delta n_{\ell}(C_{\ell} + C_{\ell-1})$;

6. estimate the optimal CV parameter(s) on each level;

7. compute/update CV estimates for $\hat{T}_{\ell}^{(\ell)}$ and $V_{\ell}^{\text{CV}}$ from samples on levels $\ell \le L$;

8. select level $\ell^* = \arg \max_{0 \le \ell \le L} \frac{V_{\ell}^{\text{CV}}}{r_{\ell} n_{\ell}^2 (C_{\ell} + C_{\ell-1})}$;

9. $\delta n_{\ell^*} \leftarrow \lfloor (r_{\ell^*} - 1) n_{\ell^*} \rfloor$, $\delta n_{\ell \ne \ell^*} \leftarrow 0$;

10. **end while**

11. **return** $\hat{\theta}_{L}^{\text{MLMC-*}}$, the MLMC-* estimate of $\theta_{L}$.

# 4 Numerical experiments

We demonstrate the value of our MLMC-MLCV method on the uncertain heat equation problem proposed in [18], with seven random input variables modeling the uncertainty in the diffusion coefficient and the initial condition. Although academic, this problem is nonetheless challenging because of its moderate stochastic dimensionality, the non-linearity of the solution with respect to the inputs, and the high variance of the random output. The test problem is summarized in section 4.1. The surrogate models used for the control variates are described in section 4.2, and the results from numerical experiments are presented and discussed in section 4.3.

## 4.1 Problem description

We consider the partial differential equation describing the time-evolution of the temperature $u(x, t; \mathbf{X})$ in a 1D rod of unit length over the time interval $[0, T]$, with uncertain (random) initial data $u_0$ and thermal diffusivity $\nu$,

$$ \left\{ \begin{array}{ll} \displaystyle \frac{\partial u(x, t; \mathbf{X})}{\partial t} = \nu(\mathbf{X}) \frac{\partial^2 u(x, t; \mathbf{X})}{\partial x^2}, & x \in \mathcal{D} := (0, 1), \quad t \in [0, T], \\ \\ u(x, 0; \mathbf{X}) = u_0(x; \mathbf{X}), & x \in \mathcal{D}, \\ \\ u(0, t; \mathbf{X}) = u(1, t; \mathbf{X}) = 0, & t \in [0, T], \end{array} \right. \quad (4.1) $$

where $\mathbf{X}: \Omega \to \Xi$ is a random vector modelling the uncertainty in the input parameters, and where $\nu(\mathbf{X}) > 0$ almost surely. The solution of eq. (4.1) may be expressed as

$$ u(x, t; \mathbf{X}) = \sum_{k=1}^{\infty} a_k(\mathbf{X}) \exp(-\nu(\mathbf{X})k^2\pi^2t) \sin(k\pi x) \quad (4.2) $$

with

$$ a_k(\mathbf{X}) = 2 \int_{\mathcal{D}} u_0(x; \mathbf{X}) \sin(k\pi x) \, dx. \quad (4.3) $$
**Algorithm 1** Simplified MLMC-* algorithm inspired by [36].

**Require:** $n_{\ell}^{\text{init}} > 1$, $r_{\ell} > 1$, surrogate models (depending on the method), and budget $\mathcal{C}$.

1. Set consumed budget to $\tilde{\mathcal{C}} = 0$ and $\delta n_{\ell} = n_{\ell}^{\text{init}}$ samples on levels $\ell \le L$;

2. **while** $\tilde{\mathcal{C}} \le \mathcal{C}$ **do**

3. compute $\delta n_{\ell}$ samples on each level by evaluating $f_{\ell}$ and the appropriate surrogates;

4. update sample size on each level: $n_{\ell} \leftarrow n_{\ell} + \delta n_{\ell}$;

5. update consumed budget: $\tilde{\mathcal{C}} \leftarrow \tilde{\mathcal{C}} + \sum_{\ell=0}^{L} \delta n_{\ell}(\mathcal{C}_{\ell} + \mathcal{C}_{\ell-1})$;

6. estimate the optimal CV parameter(s) on each level;

7. compute/update CV estimates for $\hat{T}_{\ell}^{(\ell)}$ and $\mathcal{V}_{\ell}^{\text{CV}}$ from samples on levels $\ell \le L$;

8. select level $\ell^* = \arg\max_{0 \le \ell \le L} \frac{\mathcal{V}_{\ell}^{\text{CV}}}{r_{\ell} n_{\ell}^2 (\mathcal{C}_{\ell} + \mathcal{C}_{\ell-1})}$;

9. $\delta n_{\ell^*} \leftarrow \lfloor (r_{\ell^*} - 1) n_{\ell^*} \rfloor$, $\delta n_{\ell \ne \ell^*} \leftarrow 0$;

10. **end while**

11. **return** $\hat{\theta}_{L}^{\text{MLMC-*}}$, the MLMC-* estimate of $\theta_{L}$.

# 4 Numerical experiments

We demonstrate the value of our MLMC-MLCV method on the uncertain heat equation problem proposed in [18], with seven random input variables modeling the uncertainty in the diffusion coefficient and the initial condition. Although academic, this problem is nonetheless challenging because of its moderate stochastic dimensionality, the non-linearity of the solution with respect to the inputs, and the high variance of the random output. The test problem is summarized in section 4.1. The surrogate models used for the control variates are described in section 4.2, and the results from numerical experiments are presented and discussed in section 4.3.

## 4.1 Problem description

We consider the partial differential equation describing the time-evolution of the temperature $u(x, t; \mathbf{X})$ in a 1D rod of unit length over the time interval $[0, T]$, with uncertain (random) initial data $u_0$ and thermal diffusivity $\nu$,

$$ \left\{ \begin{array}{ll} \displaystyle \frac{\partial u(x, t; \mathbf{X})}{\partial t} = \nu(\mathbf{X}) \frac{\partial^2 u(x, t; \mathbf{X})}{\partial x^2}, & x \in \mathcal{D} := (0, 1), \quad t \in [0, T], \\ \\ u(x, 0; \mathbf{X}) = u_0(x; \mathbf{X}), & x \in \mathcal{D}, \\ \\ u(0, t; \mathbf{X}) = u(1, t; \mathbf{X}) = 0, & t \in [0, T], \end{array} \right. \quad (4.1) $$

where $\mathbf{X}: \Omega \to \Xi$ is a random vector modelling the uncertainty in the input parameters, and where $\nu(\mathbf{X}) > 0$ almost surely. The solution of eq. (4.1) may be expressed as

$$ u(x, t; \mathbf{X}) = \sum_{k=1}^{\infty} a_k(\mathbf{X}) \exp(-\nu(\mathbf{X})k^2\pi^2t) \sin(k\pi x) \quad (4.2) $$

with

$$ a_k(\mathbf{X}) = 2 \int_{\mathcal{D}} u_0(x; \mathbf{X}) \sin(k\pi x) \, dx. \quad (4.3) $$
The initial condition is chosen to have the same prescribed form as in [18]. Specifically, we consider $u_0(x; \mathbf{X}) = \mathcal{G}(\mathbf{X})\mathcal{F}_1(x) + \mathcal{I}(\mathbf{X})\mathcal{F}_2(x)$ with

$$ \mathcal{F}_1(x) = \sin(\pi x), \qquad (4.4) $$

$$ \mathcal{F}_2(x) = \sin(2\pi x) + \sin(3\pi x) + 50(\sin(9\pi x) + \sin(21\pi x)), \qquad (4.5) $$

$$ \mathcal{I}(\mathbf{X}) = \frac{7}{2} [\sin(X_1) + 7 \sin(X_2)^2 + 0.1 X_3^4 \sin(X_1)], \qquad (4.6) $$

$$ \mathcal{G}(\mathbf{X}) = 50(4|X_5| - 1)(4|X_6| - 1)(4|X_7| - 1), \qquad (4.7) $$

which allows to control the spectral content of the solution $u$. Furthermore, as in [18], the diffusion coefficient is modelled by $\nu(\mathbf{X}) = X_4$. The random output variable of interest is defined as the integral of the temperature along the rod at final time $T$,

$$ \mathcal{M}(\mathbf{X}) = \int_{\mathcal{D}} u(x, T; \mathbf{X}) \, dx \qquad (4.8) $$

$$ = \sum_{k=1}^{\infty} a_k(\mathbf{X}) \int_{\mathcal{D}} \exp(-\nu(\mathbf{X})k^2\pi^2T) \sin(k\pi x) \, dx \qquad (4.9) $$

$$ = \mathcal{G}(\mathbf{X})\mathcal{H}_1(\mathbf{X}) + \mathcal{I}(\mathbf{X})[\mathcal{H}_3(\mathbf{X}) + 50\mathcal{H}_9(\mathbf{X}) + 50\mathcal{H}_{21}(\mathbf{X})], \qquad (4.10) $$

where $\mathcal{H}_k(\mathbf{X}) = \frac{2}{k\pi} \exp(-\nu(\mathbf{X})k^2\pi^2T)$. In this experiment, we seek to estimate the expectation $\mathbb{E}[\mathcal{M}(\mathbf{X})]$, for a given uncertain setting. Consistently with [18], we consider the random variables $X_1, \dots, X_7$ to be independent and distributed as

$$ X_1, X_2, X_3 \sim \mathcal{U}[-\pi, \pi], \quad X_4 \sim \mathcal{U}[\nu_{\min}, \nu_{\max}], \quad X_5, X_6, X_7 \sim \mathcal{U}[-1, 1]. \qquad (4.11) $$

The expected value $\mathbb{E}[\mathcal{M}(\mathbf{X})]$ is then given by

$$ \mathbb{E}[\mathcal{M}(\mathbf{X})] = 50H_1 + \frac{49}{4}(H_3 + 50H_9 + 50H_{21}), \qquad (4.12) $$

where

$$ H_k = \mathbb{E}[\mathcal{H}_k(\mathbf{X})] = \frac{2}{k^3\pi^3T} \frac{\exp(-\nu_{\min}k^2\pi^2T) - \exp(-\nu_{\max}k^2\pi^2T)}{\nu_{\max} - \nu_{\min}}. \qquad (4.13) $$

Finally, we set $T = 0.5$, $\nu_{\min} = 0.001$ and $\nu_{\max} = 0.009$, resulting in $\mathbb{E}[\mathcal{M}(\mathbf{X})] \approx 41.98$.

Numerically, $\mathcal{M}(\mathbf{X})$ is approximated by truncating the Fourier expansion in eq. (4.9) to $K < \infty$ modes and by approximating the integrals in eqs. (4.3) and (4.9) by a trapezoidal quadrature rule with equispaced nodes in $[0, 1]$. The multilevel hierarchy of simulators $\{f_l\}_{l=0}^L$ is then defined according to the number of quadrature nodes $N_\ell$ used for the approximation at level $\ell$. Specifically, $\mathcal{M}(\mathbf{X})$ is approximated at level $\ell$ by

$$ Y_{\ell} = f_{\ell}(\mathbf{X}) = \sum_{k=1}^{K} A_{k}^{\ell}(\mathbf{X}) B_{k}^{\ell}(\mathbf{X}), \qquad (4.14) $$

with

$$ A_k^\ell(\mathbf{X}) = 2 \sum_{i=1}^{N_\ell} w_i u_0(x_i; \mathbf{X}) \sin(k\pi x_i), \quad B_k^\ell(\mathbf{X}) = \exp(-\nu(\mathbf{X})k^2\pi^2T) \sum_{i=1}^{N_\ell} w_i \sin(k\pi x_i), \qquad (4.15) $$

where $\{(x_i, w_i)\}_{i=1}^{N_\ell}$ are the pairs of quadrature nodes and associated weights on level $\ell$. It is then natural to assume that the computational cost $C_\ell$ of an evaluation of $f_\ell$ is $\mathcal{O}(KN_\ell)$. The statistic
The initial condition is chosen to have the same prescribed form as in [18]. Specifically, we consider $u_0(x; \mathbf{X}) = \mathcal{G}(\mathbf{X})\mathcal{F}_1(x) + \mathcal{I}(\mathbf{X})\mathcal{F}_2(x)$ with

$$ \mathcal{F}_1(x) = \sin(\pi x), \qquad (4.4) $$

$$ \mathcal{F}_2(x) = \sin(2\pi x) + \sin(3\pi x) + 50(\sin(9\pi x) + \sin(21\pi x)), \qquad (4.5) $$

$$ \mathcal{I}(\mathbf{X}) = \frac{7}{2} [\sin(X_1) + 7 \sin(X_2)^2 + 0.1 X_3^4 \sin(X_1)], \qquad (4.6) $$

$$ \mathcal{G}(\mathbf{X}) = 50(4|X_5| - 1)(4|X_6| - 1)(4|X_7| - 1), \qquad (4.7) $$

which allows to control the spectral content of the solution $u$. Furthermore, as in [18], the diffusion coefficient is modelled by $\nu(\mathbf{X}) = X_4$. The random output variable of interest is defined as the integral of the temperature along the rod at final time $T$,

$$ \mathcal{M}(\mathbf{X}) = \int_{\mathcal{D}} u(x, T; \mathbf{X}) \, dx \qquad (4.8) $$

$$ = \sum_{k=1}^{\infty} a_k(\mathbf{X}) \int_{\mathcal{D}} \exp(-\nu(\mathbf{X})k^2\pi^2T) \sin(k\pi x) \, dx \qquad (4.9) $$

$$ = \mathcal{G}(\mathbf{X})\mathcal{H}_1(\mathbf{X}) + \mathcal{I}(\mathbf{X})[\mathcal{H}_3(\mathbf{X}) + 50\mathcal{H}_9(\mathbf{X}) + 50\mathcal{H}_{21}(\mathbf{X})], \qquad (4.10) $$

where $\mathcal{H}_k(\mathbf{X}) = \frac{2}{k\pi} \exp(-\nu(\mathbf{X})k^2\pi^2T)$. In this experiment, we seek to estimate the expectation $\mathbb{E}[\mathcal{M}(\mathbf{X})]$, for a given uncertain setting. Consistently with [18], we consider the random variables $X_1, \dots, X_7$ to be independent and distributed as

$$ X_1, X_2, X_3 \sim \mathcal{U}[-\pi, \pi], \quad X_4 \sim \mathcal{U}[\nu_{\min}, \nu_{\max}], \quad X_5, X_6, X_7 \sim \mathcal{U}[-1, 1]. \qquad (4.11) $$

The expected value $\mathbb{E}[\mathcal{M}(\mathbf{X})]$ is then given by

$$ \mathbb{E}[\mathcal{M}(\mathbf{X})] = 50H_1 + \frac{49}{4}(H_3 + 50H_9 + 50H_{21}), \qquad (4.12) $$

where

$$ H_k = \mathbb{E}[\mathcal{H}_k(\mathbf{X})] = \frac{2}{k^3\pi^3T} \frac{\exp(-\nu_{\min}k^2\pi^2T) - \exp(-\nu_{\max}k^2\pi^2T)}{\nu_{\max} - \nu_{\min}}. \qquad (4.13) $$

Finally, we set $T = 0.5$, $\nu_{\min} = 0.001$ and $\nu_{\max} = 0.009$, resulting in $\mathbb{E}[\mathcal{M}(\mathbf{X})] \approx 41.98$.

Numerically, $\mathcal{M}(\mathbf{X})$ is approximated by truncating the Fourier expansion in eq. (4.9) to $K < \infty$ modes and by approximating the integrals in eqs. (4.3) and (4.9) by a trapezoidal quadrature rule with equispaced nodes in $[0, 1]$. The multilevel hierarchy of simulators $\{f_l\}_{l=0}^L$ is then defined according to the number of quadrature nodes $N_\ell$ used for the approximation at level $\ell$. Specifically, $\mathcal{M}(\mathbf{X})$ is approximated at level $\ell$ by

$$ Y_{\ell} = f_{\ell}(\mathbf{X}) = \sum_{k=1}^{K} A_{k}^{\ell}(\mathbf{X}) B_{k}^{\ell}(\mathbf{X}), \qquad (4.14) $$

with

$$ A_k^\ell(\mathbf{X}) = 2 \sum_{i=1}^{N_\ell} w_i u_0(x_i; \mathbf{X}) \sin(k\pi x_i), \quad B_k^\ell(\mathbf{X}) = \exp(-\nu(\mathbf{X})k^2\pi^2T) \sum_{i=1}^{N_\ell} w_i \sin(k\pi x_i), \qquad (4.15) $$

where $\{(x_i, w_i)\}_{i=1}^{N_\ell}$ are the pairs of quadrature nodes and associated weights on level $\ell$. It is then natural to assume that the computational cost $C_\ell$ of an evaluation of $f_\ell$ is $\mathcal{O}(KN_\ell)$. The statistic
of interest is thus $\theta = \theta_L = \mathbb{E}[f_L(\mathbf{X})]$, whose MC estimator $\hat{\theta}_L = n_L^{-1} \sum_{i=1}^{n_L} f_L(\mathbf{X}^{(L,i)})$ will represent the baseline estimator for our experiments. Besides, the quality of all the presented estimators will be assessed in terms of their root mean square error (RMSE) w.r.t. the exact statistic $\mathbb{E}[\mathcal{M}(\mathbf{X})]$ given by eq. (4.12).

In the following experiments, we set the number of quadrature nodes $K = 21$ and the number of levels to 4 (i.e. $L = 3$). Furthermore, we set $N_\ell = 120 \times 2^{L-\ell}$ so that evaluating $f_\ell$ is twice as expensive as evaluating $f_{\ell-1}$. Table 2 summarizes the number of quadrature nodes and the evaluation cost per level. Note that the costs are normalized so that $C_3 = 1$.

<table><caption>Table 2: Number of quadrature nodes $N_\ell$, simulator evaluation cost $C_\ell$ and MLMC correction evaluation cost $C_\ell + C_{\ell-1}$ per level.</caption><thead><tr><th>ℓ</th><th>0</th><th>1</th><th>2</th><th>3</th></tr></thead><tbody><tr><th>N<sub>ℓ</sub></th><td>15</td><td>30</td><td>60</td><td>120</td></tr><tr><th>C<sub>ℓ</sub></th><td>0.125</td><td>0.25</td><td>0.5</td><td>1</td></tr><tr><th>C<sub>ℓ</sub> + C<sub>ℓ-1</sub></th><td>0.125</td><td>0.375</td><td>0.75</td><td>1.5</td></tr></tbody></table>

## 4.2 Surrogate models

We will mostly use PC models (see section 2.4.3) for the surrogate-based CV estimators. Constructing high-quality surrogates in 7 dimensions can be hard, especially in the presence of non-linearities and with a limited sample size. To avoid overfitting, we resort to the least angle regression (LARS) procedure [14], which is a model-selection regression method that promotes sparsity. More precisely, we employ the basis-adaptive hybrid LARS algorithm proposed by [7, Fig. 5] for the selection of sparse PC bases. For a given design of experiment (DoE) for the PC surrogate construction, this algorithm applies the LARS procedure on candidate PC bases $A_p$ of increasing total polynomial degree $p = 1, \dots, p_{max}$, resulting for each $p$ in the selection of a limited number $|\tilde{A}_p| < |A_p|$ of basis polynomial functions. For each $p$, a PC surrogate is constructed by classical least-squares regression on the corresponding reduced (or active) PC basis $\tilde{A}_p$, and its quality is estimated using a corrected leave-one-out cross-validation procedure [7]. The best surrogate according to this quality measure is eventually retained, and the associated reduced PC basis is denoted by $\tilde{A}_{p*}$.

In our experiments, we set $p_{max} = 16$ and the construction budget to 400 times the evaluation cost of $f_3$, i.e. $C^{DoE} = 400$. This budget is distributed equally among the different levels, so that the associated evaluation cost on each level corresponds to the cost of 100 $f_3$ evaluations, i.e. $n_\ell^{DoE}C_\ell = 100$, as reported in tables 3 and 4. Once constructed, the quality of a surrogate $g$ of $f$ is assessed in terms of the $Q^2$ measure,

$$Q^2(g, f) = 1 - \frac{\mathbb{E}[(g(\mathbf{X}) - f(\mathbf{X}))^2]}{\mathbb{V}[f(\mathbf{X})]}, \quad (4.16)$$

estimated using a test sample of size $n_{test} = 10000$, which is more robust than the corrected leave-one-out measure used for the model selection. The higher the $Q^2$ value, the higher-quality the associated surrogate.

For the MLCV estimator, we learn a PC model $g_\ell$ of $f_\ell$ from a training DoE of size $n_\ell^{DoE}$ generated by latin hypercube sampling (LHS) improved by simulated annealing. The training DoE samples for the surrogate construction are different for each level and independent. Table 3 summarizes the properties of the different PC models. Except for $g_3$, which is built with only 100 points, all the PC
of interest is thus $\theta = \theta_L = \mathbb{E}[f_L(\mathbf{X})]$, whose MC estimator $\hat{\theta}_L = n_L^{-1} \sum_{i=1}^{n_L} f_L(\mathbf{X}^{(L,i)})$ will represent the baseline estimator for our experiments. Besides, the quality of all the presented estimators will be assessed in terms of their root mean square error (RMSE) w.r.t. the exact statistic $\mathbb{E}[\mathcal{M}(\mathbf{X})]$ given by eq. (4.12).

In the following experiments, we set the number of quadrature nodes $K = 21$ and the number of levels to 4 (i.e. $L = 3$). Furthermore, we set $N_\ell = 120 \times 2^{L-\ell}$ so that evaluating $f_\ell$ is twice as expensive as evaluating $f_{\ell-1}$. Table 2 summarizes the number of quadrature nodes and the evaluation cost per level. Note that the costs are normalized so that $C_3 = 1$.

<table><caption>Table 2: Number of quadrature nodes $N_\ell$, simulator evaluation cost $C_\ell$ and MLMC correction evaluation cost $C_\ell + C_{\ell-1}$ per level.</caption><thead><tr><th>ℓ</th><th>0</th><th>1</th><th>2</th><th>3</th></tr></thead><tbody><tr><th>N<sub>ℓ</sub></th><td>15</td><td>30</td><td>60</td><td>120</td></tr><tr><th>C<sub>ℓ</sub></th><td>0.125</td><td>0.25</td><td>0.5</td><td>1</td></tr><tr><th>C<sub>ℓ</sub> + C<sub>ℓ-1</sub></th><td>0.125</td><td>0.375</td><td>0.75</td><td>1.5</td></tr></tbody></table>

## 4.2 Surrogate models

We will mostly use PC models (see section 2.4.3) for the surrogate-based CV estimators. Constructing high-quality surrogates in 7 dimensions can be hard, especially in the presence of non-linearities and with a limited sample size. To avoid overfitting, we resort to the least angle regression (LARS) procedure [14], which is a model-selection regression method that promotes sparsity. More precisely, we employ the basis-adaptive hybrid LARS algorithm proposed by [7, Fig. 5] for the selection of sparse PC bases. For a given design of experiment (DoE) for the PC surrogate construction, this algorithm applies the LARS procedure on candidate PC bases $A_p$ of increasing total polynomial degree $p = 1, \dots, p_{max}$, resulting for each $p$ in the selection of a limited number $|\tilde{A}_p| < |A_p|$ of basis polynomial functions. For each $p$, a PC surrogate is constructed by classical least-squares regression on the corresponding reduced (or active) PC basis $\tilde{A}_p$, and its quality is estimated using a corrected leave-one-out cross-validation procedure [7]. The best surrogate according to this quality measure is eventually retained, and the associated reduced PC basis is denoted by $\tilde{A}_{p*}$.

In our experiments, we set $p_{max} = 16$ and the construction budget to 400 times the evaluation cost of $f_3$, i.e. $C^{DoE} = 400$. This budget is distributed equally among the different levels, so that the associated evaluation cost on each level corresponds to the cost of 100 $f_3$ evaluations, i.e. $n_\ell^{DoE}C_\ell = 100$, as reported in tables 3 and 4. Once constructed, the quality of a surrogate $g$ of $f$ is assessed in terms of the $Q^2$ measure,

$$Q^2(g, f) = 1 - \frac{\mathbb{E}[(g(\mathbf{X}) - f(\mathbf{X}))^2]}{\mathbb{V}[f(\mathbf{X})]}, \quad (4.16)$$

estimated using a test sample of size $n_{test} = 10000$, which is more robust than the corrected leave-one-out measure used for the model selection. The higher the $Q^2$ value, the higher-quality the associated surrogate.

For the MLCV estimator, we learn a PC model $g_\ell$ of $f_\ell$ from a training DoE of size $n_\ell^{DoE}$ generated by latin hypercube sampling (LHS) improved by simulated annealing. The training DoE samples for the surrogate construction are different for each level and independent. Table 3 summarizes the properties of the different PC models. Except for $g_3$, which is built with only 100 points, all the PC
models have a good $Q^2$, greater than 0.8. We observe that the basis-adaptive LARS algorithm has selected a decreasing polynomial degree $p^*$ with level $\ell$, while retaining only a limited number $|\tilde{A}_{p^*}|$ of polynomials in these reduced bases. Thus, although $g_1$ and $g_2$ have distinct degrees, the sizes of the associated reduced bases are similar.

<table><caption>Table 3: PC models for CV and MLCV estimators, with their sample size, degree and quality measure. These models are built with the basis-adaptive LARS algorithm of [7, Fig. 5].</caption><thead><tr><th>PC models</th><th>g<sub>0</sub></th><th>g<sub>1</sub></th><th>g<sub>2</sub></th><th>g<sub>3</sub></th></tr></thead><tbody><tr><td>n<sub>ℓ</sub><sup>DoE</sup></td><td>800</td><td>400</td><td>200</td><td>100</td></tr><tr><td>p<sup>*</sup></td><td>10</td><td>8</td><td>6</td><td>4</td></tr><tr><td>|Ã<sub>p<sup>*</sup></sub>|</td><td>211</td><td>77</td><td>74</td><td>14</td></tr><tr><td>Q<sup>2</sup></td><td>0.98</td><td>0.94</td><td>0.86</td><td>0.59</td></tr></tbody></table>

For the MLMC-based estimators, we learn a PC model $g_ℓ$ of $f_ℓ$ for $ℓ = 0, \dots, L$ and a PC model $h_ℓ$ of $f_ℓ - f_{ℓ-1}$ for $ℓ = 1, \dots, L$. In practice, training points used for the construction of $g_0, \dots, g_L$ may be reused for the construction of $h_1, \dots, h_L$ using nested DoEs. However, it should be noted that the generation of nested LHS DoEs is not as straightforward as for purely random DoEs. While using nested random DoEs is a perfectly valid strategy, we opt for an alternative choice based on LHS. Specifically, we first generate a DoE $\mathcal{X}_{\text{DoE}}^{(0)}$ of size $n_0^{\text{DoE}}$ using LHS improved by simulated annealing. Then, for $ℓ = 1, \dots, L$, we sequentially extract a DoE $\mathcal{X}_{\text{DoE}}^{(ℓ)} \subset \mathcal{X}_{\text{DoE}}^{(ℓ-1)}$ of size $n_ℓ^{\text{DoE}}$. The subset $\mathcal{X}_{\text{DoE}}^{(ℓ)}$ is selected such that it has minimal centered $L^2$ discrepancy among a pool of 1 million random candidate subsets. Table 4 summarizes the properties of these different PC models. The PC models again have a quality measure higher than 0.8, except for $g_3$ and $h_3$, which were constructed from only 100 evaluations. Note that the values of $p^*$, $|\tilde{A}_{p^*}|$ and $Q^2$ are of the same order of magnitude as those of the PC surrogate models reported in table 3.

<table><caption>Table 4: PC models for the estimators combining MLMC and CV techniques, with their sample size, degree and quality measure. More precisely, $g_ℓ$ is the surrogate of $f_ℓ$, $h_ℓ$ is the surrogate of $f_ℓ - f_{ℓ-1}$ and the samples used to train $h_ℓ$ are also used to train $g_ℓ$ and $g_{ℓ-1}$. These models are built with the basis-adaptive LARS algorithm of [7, Fig. 5].</caption><thead><tr><th>PC models</th><th>g<sub>0</sub></th><th>g<sub>1</sub></th><th>g<sub>2</sub></th><th>g<sub>3</sub></th><th>h<sub>1</sub></th><th>h<sub>2</sub></th><th>h<sub>3</sub></th></tr></thead><tbody><tr><td>n<sub>ℓ</sub><sup>DoE</sup></td><td>800</td><td>400</td><td>200</td><td>100</td><td>400</td><td>200</td><td>100</td></tr><tr><td>p<sup>*</sup></td><td>10</td><td>8</td><td>8</td><td>2</td><td>9</td><td>7</td><td>5</td></tr><tr><td>|Ã<sub>p<sup>*</sup></sub>|</td><td>211</td><td>83</td><td>57</td><td>6</td><td>112</td><td>71</td><td>30</td></tr><tr><td>Q<sup>2</sup></td><td>0.98</td><td>0.95</td><td>0.86</td><td>0.21</td><td>0.99</td><td>0.82</td><td>0.63</td></tr></tbody></table>

Hereafter, the CV and MLCV estimators use the PC models of table 3, while the MLMC-based estimators use the PC models presented in table 4.

## 4.3 Results

First, in section 4.3.1, we illustrate the use of one or several control variates to reduce the variance of a single-level MC estimator. Then, we compare the MLCV and MLMC-MLCV approaches with the MC and MLMC estimators in section 4.3.2, and we discuss variants of the MLMC-CV and MLMC-MLCV approaches, considering only a limited subset of the surrogate models, in section 4.3.3. We
models have a good $Q^2$, greater than 0.8. We observe that the basis-adaptive LARS algorithm has selected a decreasing polynomial degree $p^*$ with level $\ell$, while retaining only a limited number $|\tilde{\mathcal{A}}_{p^*}|$ of polynomials in these reduced bases. Thus, although $g_1$ and $g_2$ have distinct degrees, the sizes of the associated reduced bases are similar.

<table><caption>Table 3: PC models for CV and MLCV estimators, with their sample size, degree and quality measure. These models are built with the basis-adaptive LARS algorithm of [7, Fig. 5].</caption><thead><tr><th>PC models</th><th>g<sub>0</sub></th><th>g<sub>1</sub></th><th>g<sub>2</sub></th><th>g<sub>3</sub></th></tr></thead><tbody><tr><td>n<sub>ℓ</sub><sup>DoE</sup></td><td>800</td><td>400</td><td>200</td><td>100</td></tr><tr><td>p<sup>*</sup></td><td>10</td><td>8</td><td>6</td><td>4</td></tr><tr><td>|Ã<sub>p<sup>*</sup></sub>|</td><td>211</td><td>77</td><td>74</td><td>14</td></tr><tr><td>Q<sup>2</sup></td><td>0.98</td><td>0.94</td><td>0.86</td><td>0.59</td></tr></tbody></table>

For the MLMC-based estimators, we learn a PC model $g_ℓ$ of $f_ℓ$ for $ℓ = 0, \dots, L$ and a PC model $h_ℓ$ of $f_ℓ - f_{ℓ-1}$ for $ℓ = 1, \dots, L$. In practice, training points used for the construction of $g_0, \dots, g_L$ may be reused for the construction of $h_1, \dots, h_L$ using nested DoEs. However, it should be noted that the generation of nested LHS DoEs is not as straightforward as for purely random DoEs. While using nested random DoEs is a perfectly valid strategy, we opt for an alternative choice based on LHS. Specifically, we first generate a DoE $\mathcal{X}_{\text{DoE}}^{(0)}$ of size $n_0^{\text{DoE}}$ using LHS improved by simulated annealing. Then, for $ℓ = 1, \dots, L$, we sequentially extract a DoE $\mathcal{X}_{\text{DoE}}^{(ℓ)} \subset \mathcal{X}_{\text{DoE}}^{(ℓ-1)}$ of size $n_ℓ^{\text{DoE}}$. The subset $\mathcal{X}_{\text{DoE}}^{(ℓ)}$ is selected such that it has minimal centered $L^2$ discrepancy among a pool of 1 million random candidate subsets. Table 4 summarizes the properties of these different PC models. The PC models again have a quality measure higher than 0.8, except for $g_3$ and $h_3$, which were constructed from only 100 evaluations. Note that the values of $p^*$, $|\tilde{\mathcal{A}}_{p^*}|$ and $Q^2$ are of the same order of magnitude as those of the PC surrogate models reported in table 3.

<table><caption>Table 4: PC models for the estimators combining MLMC and CV techniques, with their sample size, degree and quality measure. More precisely, $g_ℓ$ is the surrogate of $f_ℓ$, $h_ℓ$ is the surrogate of $f_ℓ - f_{ℓ-1}$ and the samples used to train $h_ℓ$ are also used to train $g_ℓ$ and $g_{ℓ-1}$. These models are built with the basis-adaptive LARS algorithm of [7, Fig. 5].</caption><thead><tr><th>PC models</th><th>g<sub>0</sub></th><th>g<sub>1</sub></th><th>g<sub>2</sub></th><th>g<sub>3</sub></th><th>h<sub>1</sub></th><th>h<sub>2</sub></th><th>h<sub>3</sub></th></tr></thead><tbody><tr><td>n<sub>ℓ</sub><sup>DoE</sup></td><td>800</td><td>400</td><td>200</td><td>100</td><td>400</td><td>200</td><td>100</td></tr><tr><td>p<sup>*</sup></td><td>10</td><td>8</td><td>8</td><td>2</td><td>9</td><td>7</td><td>5</td></tr><tr><td>|Ã<sub>p<sup>*</sup></sub>|</td><td>211</td><td>83</td><td>57</td><td>6</td><td>112</td><td>71</td><td>30</td></tr><tr><td>Q<sup>2</sup></td><td>0.98</td><td>0.95</td><td>0.86</td><td>0.21</td><td>0.99</td><td>0.82</td><td>0.63</td></tr></tbody></table>

Hereafter, the CV and MLCV estimators use the PC models of table 3, while the MLMC-based estimators use the PC models presented in table 4.

## 4.3 Results

First, in section 4.3.1, we illustrate the use of one or several control variates to reduce the variance of a single-level MC estimator. Then, we compare the MLCV and MLMC-MLCV approaches with the MC and MLMC estimators in section 4.3.2, and we discuss variants of the MLMC-CV and MLMC-MLCV approaches, considering only a limited subset of the surrogate models, in section 4.3.3. We
conclude the analysis by reporting the estimation budget allocation across levels resulting from the various MLMC-based methods in section 4.3.4. In practice, all the MLMC-based estimators are built using algorithm 1. Unless stated otherwise, the parameters are set to $n_{\ell}^{\text{init}} = 30$ and $r_{\ell} = 1.1$ for $\ell = 0, \dots, 3$. The quality of the various estimators will be assessed in terms of their RMSE w.r.t. $\mathbb{E}[\mathcal{M}(\mathbf{X})]$, estimated from 500 repetitions of the experiment.

### 4.3.1 Single-level MC and control variates

In this first part, we consider only the finest simulator $f_3$ and try to reduce the variance of the MC estimator of $\mathbb{E}[f_3(\mathbf{X})]$ by means of surrogate-based CVs. For that purpose, we consider a first-order Taylor polynomial expansion $g_3^{T_1}$ around $\mu_{\mathbf{X}} = \mathbb{E}[\mathbf{X}]$ (see section 2.4.2 and appendix D) and the PC model $g_3^{\text{PC}}$ described in table 3.

<table><caption>Table 5: Pearson correlation coefficients between the finest simulator $f_3$ and the corresponding PC model (see table 3) and $T_1$ model around $\mu_{\mathbf{X}} = \mathbb{E}[\mathbf{X}]$, estimated with a sample of size 1000.</caption><thead><tr><th></th><th>Y<sub>3</sub></th><th>g<sub>3</sub><sup>PC</sup>(X)</th><th>g<sub>3</sub><sup>T<sub>1</sub></sup>(X)</th></tr></thead><tbody><tr><th>Y<sub>3</sub></th><td>1.00</td><td>0.80</td><td>0.57</td></tr><tr><th>g<sub>3</sub><sup>PC</sup>(X)</th><td>0.80</td><td>1.00</td><td>0.64</td></tr><tr><th>g<sub>3</sub><sup>T<sub>1</sub></sup>(X)</th><td>0.57</td><td>0.64</td><td>1.00</td></tr></tbody></table>

Table 5 shows that $g_3^{\text{PC}}(\mathbf{X})$ is well correlated with $Y_3 = f_3(\mathbf{X})$, with a Pearson coefficient of 0.8, while $g_3^{T_1}$ is less so, with a coefficient of 0.57, which may be explained by the strong non-linearity of $f_3$. According to eq. (2.7), these correlation coefficients lead to a theoretical variance reduction factor $R^2$ of 64% when using $g_3^{\text{PC}}(\mathbf{X})$ as a single CV, and of about 32% when using $g_3^{T_1}(\mathbf{X})$ as a single CV. Using both CVs results in a reduction factor of 65%. This minor increase in $R^2$ can be explained by the modest correlation coefficient of 0.64 between $g_3^{T_1}(\mathbf{X})$ and $g_3^{\text{PC}}(\mathbf{X})$. Note that a reduction factor of $R^2$ in the variance corresponds to a reduction factor of $1 - \sqrt{1 - R^2}$ in the standard deviation.

These theoretical expectations are reflected in fig. 1a with an RMSE reduction of about 20% when using $g_3^{T_1}$ and 40% when using $g_3^{\text{PC}}$ alone or jointly with $g_3^{T_1}$. This figure confirms that $g_3^{\text{PC}}$ provides a better CV than $g_3^{T_1}$, reducing the RMSE of the MC estimator twice as much, regardless of the computational budget. However, the construction cost of the surrogate is not the same. While constructing $g_3^{T_1}$ requires only one evaluation of $f_3$ and its Jacobian matrix, namely at $\mu_{\mathbf{X}}$, the construction of $g_3^{\text{PC}}$ involved 100 $f_3$ evaluations. In the case where the surrogate is built specifically for the estimation of the statistic, the real estimation cost is higher as it includes this construction cost. Figure 1b illustrates this difference by including the surrogate construction cost in the total evaluation cost. As a result, a significant offset appears when using $g_3^{\text{PC}}$, since part of the computational budget (namely 100) is used for the surrogate construction and thus not for the estimation. Specifically, the RMSEs of the CV estimators using $g_3^{\text{PC}}$ get below that of the CV estimator using only $g_3^{T_1}$ for a budget of 300 $f_3$ evaluations. For budgets under 200, using only $g_3^{T_1}$ is preferable, even without the analytical gradient available, which would then incur a construction cost of 8 $f_3$ evaluations to approximate the gradient using finite differences. This figure also illustrates proposition 2.1, that is, increasing the number of control variates improves (rigorously speaking, does not deteriorate) the variance of the CV estimator.
conclude the analysis by reporting the estimation budget allocation across levels resulting from the various MLMC-based methods in section 4.3.4. In practice, all the MLMC-based estimators are built using algorithm 1. Unless stated otherwise, the parameters are set to $n_{\ell}^{\text{init}} = 30$ and $r_{\ell} = 1.1$ for $\ell = 0, \dots, 3$. The quality of the various estimators will be assessed in terms of their RMSE w.r.t. $\mathbb{E}[\mathcal{M}(\mathbf{X})]$, estimated from 500 repetitions of the experiment.

### 4.3.1 Single-level MC and control variates

In this first part, we consider only the finest simulator $f_3$ and try to reduce the variance of the MC estimator of $\mathbb{E}[f_3(\mathbf{X})]$ by means of surrogate-based CVs. For that purpose, we consider a first-order Taylor polynomial expansion $g_3^{T_1}$ around $\mu_{\mathbf{X}} = \mathbb{E}[\mathbf{X}]$ (see section 2.4.2 and appendix D) and the PC model $g_3^{\text{PC}}$ described in table 3.

<table><caption>Table 5: Pearson correlation coefficients between the finest simulator $f_3$ and the corresponding PC model (see table 3) and $T_1$ model around $\mu_{\mathbf{X}} = \mathbb{E}[\mathbf{X}]$, estimated with a sample of size 1000.</caption><thead><tr><th></th><th>Y<sub>3</sub></th><th>g<sub>3</sub><sup>PC</sup>(<b>X</b>)</th><th>g<sub>3</sub><sup>T<sub>1</sub></sup>(<b>X</b>)</th></tr></thead><tbody><tr><td>Y<sub>3</sub></td><td>1.00</td><td>0.80</td><td>0.57</td></tr><tr><td>g<sub>3</sub><sup>PC</sup>(<b>X</b>)</td><td>0.80</td><td>1.00</td><td>0.64</td></tr><tr><td>g<sub>3</sub><sup>T<sub>1</sub></sup>(<b>X</b>)</td><td>0.57</td><td>0.64</td><td>1.00</td></tr></tbody></table>

Table 5 shows that $g_3^{\text{PC}}(\mathbf{X})$ is well correlated with $Y_3 = f_3(\mathbf{X})$, with a Pearson coefficient of 0.8, while $g_3^{T_1}$ is less so, with a coefficient of 0.57, which may be explained by the strong non-linearity of $f_3$. According to eq. (2.7), these correlation coefficients lead to a theoretical variance reduction factor $R^2$ of 64% when using $g_3^{\text{PC}}(\mathbf{X})$ as a single CV, and of about 32% when using $g_3^{T_1}(\mathbf{X})$ as a single CV. Using both CVs results in a reduction factor of 65%. This minor increase in $R^2$ can be explained by the modest correlation coefficient of 0.64 between $g_3^{T_1}(\mathbf{X})$ and $g_3^{\text{PC}}(\mathbf{X})$. Note that a reduction factor of $R^2$ in the variance corresponds to a reduction factor of $1 - \sqrt{1 - R^2}$ in the standard deviation.

These theoretical expectations are reflected in fig. 1a with an RMSE reduction of about 20% when using $g_3^{T_1}$ and 40% when using $g_3^{\text{PC}}$ alone or jointly with $g_3^{T_1}$. This figure confirms that $g_3^{\text{PC}}$ provides a better CV than $g_3^{T_1}$, reducing the RMSE of the MC estimator twice as much, regardless of the computational budget. However, the construction cost of the surrogate is not the same. While constructing $g_3^{T_1}$ requires only one evaluation of $f_3$ and its Jacobian matrix, namely at $\mu_{\mathbf{X}}$, the construction of $g_3^{\text{PC}}$ involved 100 $f_3$ evaluations. In the case where the surrogate is built specifically for the estimation of the statistic, the real estimation cost is higher as it includes this construction cost. Figure 1b illustrates this difference by including the surrogate construction cost in the total evaluation cost. As a result, a significant offset appears when using $g_3^{\text{PC}}$, since part of the computational budget (namely 100) is used for the surrogate construction and thus not for the estimation. Specifically, the RMSEs of the CV estimators using $g_3^{\text{PC}}$ get below that of the CV estimator using only $g_3^{T_1}$ for a budget of 300 $f_3$ evaluations. For budgets under 200, using only $g_3^{T_1}$ is preferable, even without the analytical gradient available, which would then incur a construction cost of 8 $f_3$ evaluations to approximate the gradient using finite differences. This figure also illustrates proposition 2.1, that is, increasing the number of control variates improves (rigorously speaking, does not deteriorate) the variance of the CV estimator.
<figure><img src="image_2.png" /><figcaption>(a) The evaluation cost does not include the cost of constructing the surrogate models.</figcaption></figure>

<caption>(b) The evaluation cost includes the cost of constructing the surrogate models.</caption>

Figure 1: RMSE of the MC and CV estimators of $\theta = E[f_3(\mathbf{X})]$ (see eq. (4.14)) with respect to the number $C \in \{100; 300; 1000; 3000; 10000\}$ of $f_3$ evaluations. The CV uses either a first-order Taylor approximation (T1) of $f_3$ or a PC approximation of $f_3$ trained from 100 evaluations ($C^{\text{DoE}} = 100$), or both. The RMSEs are computed using 500 replicates.

### 4.3.2 MLCV and MLMC-MLCV

Figure 2 compares the MLCV and MLMC-MLCV estimators proposed in eqs. (3.9) and (3.20) with the classical MC and MLMC estimators. This comparison is repeated for different budgets $C$, expressed in terms of the equivalent number of $f_3$ evaluations. Note that the MC and MLCV estimators only use evaluations of the finest simulator $f_3$, so that $C = n_3$, while MLMC-based estimators use all the simulators, so that $C$ is given in terms of the MLMC cost eq. (3.6), indeed corresponding to the equivalent number of $f_3$ evaluations, since the costs are normalized such that $C_3 = 1$. From here on, only PC-based surrogates will be used, so that we omit the superscript “PC” in the notations of the surrogate models.

**MLMC vs. MC** Figure 2a shows that, for a given sampling budget $C$, the MC estimator is the least accurate. This can be explained by the fact that it only has access to the finest simulator, $f_3$, whose cost only allows a limited number of evaluations. On the contrary, the MLMC estimator spreads this sampling budget over the four simulators. Ideally, the optimal sample allocation of MLMC eq. (3.35) results in many coarse, cheap evaluations and few fine, expensive evaluations. This is typically the case when the outputs of the simulators are highly correlated and the associated computational cost grows exponentially. Table 7 shows that the first assumption holds. The second assumption also holds since $C_\ell = O(KN_\ell) = O(N_\ell)$, $K$ being fixed, i.e., the evaluation cost grows linearly with the number of quadrature nodes. As a consequence, the MLMC estimator has a lower RMSE than the standard MC estimator.

**MLCV vs. MLMC** For this experiment, the MLCV estimator is more accurate than the MLMC estimator. Based on one control variate per level $\ell$, based on the PC model $g_\ell$ of $f_\ell$ from table 3, this estimator dedicates all the sampling budget to the finest simulator $f_3$, and uses these control variates based on $g_0, g_1, g_2$ and $g_3$ to reduce the variance of the MC estimator at no extra cost, as the evaluation cost of a PC model $g_\ell$ is negligible compared to the evaluation cost of $f_3$. This MLCV
<figure><img src="image_2.png" /><figcaption>(a) The evaluation cost does not include the cost of constructing the surrogate models.</figcaption></figure>

<caption>(b) The evaluation cost includes the cost of constructing the surrogate models.</caption>

Figure 1: RMSE of the MC and CV estimators of $\theta = E[f_3(\mathbf{X})]$ (see eq. (4.14)) with respect to the number $C \in \{100; 300; 1000; 3000; 10000\}$ of $f_3$ evaluations. The CV uses either a first-order Taylor approximation (T1) of $f_3$ or a PC approximation of $f_3$ trained from 100 evaluations ($C^{\text{DoE}} = 100$), or both. The RMSEs are computed using 500 replicates.

### 4.3.2 MLCV and MLMC-MLCV

Figure 2 compares the MLCV and MLMC-MLCV estimators proposed in eqs. (3.9) and (3.20) with the classical MC and MLMC estimators. This comparison is repeated for different budgets $C$, expressed in terms of the equivalent number of $f_3$ evaluations. Note that the MC and MLCV estimators only use evaluations of the finest simulator $f_3$, so that $C = n_3$, while MLMC-based estimators use all the simulators, so that $C$ is given in terms of the MLMC cost eq. (3.6), indeed corresponding to the equivalent number of $f_3$ evaluations, since the costs are normalized such that $C_3 = 1$. From here on, only PC-based surrogates will be used, so that we omit the superscript “PC” in the notations of the surrogate models.

**MLMC vs. MC** Figure 2a shows that, for a given sampling budget $C$, the MC estimator is the least accurate. This can be explained by the fact that it only has access to the finest simulator, $f_3$, whose cost only allows a limited number of evaluations. On the contrary, the MLMC estimator spreads this sampling budget over the four simulators. Ideally, the optimal sample allocation of MLMC eq. (3.35) results in many coarse, cheap evaluations and few fine, expensive evaluations. This is typically the case when the outputs of the simulators are highly correlated and the associated computational cost grows exponentially. Table 7 shows that the first assumption holds. The second assumption also holds since $C_\ell = O(KN_\ell) = O(N_\ell)$, $K$ being fixed, i.e., the evaluation cost grows linearly with the number of quadrature nodes. As a consequence, the MLMC estimator has a lower RMSE than the standard MC estimator.

**MLCV vs. MLMC** For this experiment, the MLCV estimator is more accurate than the MLMC estimator. Based on one control variate per level $\ell$, based on the PC model $g_\ell$ of $f_\ell$ from table 3, this estimator dedicates all the sampling budget to the finest simulator $f_3$, and uses these control variates based on $g_0, g_1, g_2$ and $g_3$ to reduce the variance of the MC estimator at no extra cost, as the evaluation cost of a PC model $g_\ell$ is negligible compared to the evaluation cost of $f_3$. This MLCV
<figure><img src="image_3.png" /></figure>

(a) The evaluation cost does not include the cost of construction of the surrogate models.

(b) The evaluation cost includes the cost of construction of the surrogate models.

<caption>Figure 2: RMSE of the MC, MLMC, MLCV and MLMC-MLCV estimators of $\theta = E[f_3(\mathbf{X})]$ (see eq. (4.14)) with respect to the sampling budget $C \in \{100; 300; 1000; 3000; 10000\}$ of $f_3$ evaluations. The surrogates used for the MLCV and MLMC-MLCV estimators are described in tables 3 and 4, respectively, and their total construction cost is $C^{\text{DoE}} = 400$. The RMSEs are estimated using 500 replicates.</caption>

<table><caption>Table 6: Pearson correlation coefficients between the finest output $Y_3 = f_3(\mathbf{X})$ and the corresponding control variates based on the PC models of table 3 used for the MLCV method, estimated with a sample of size 1000.</caption><thead><tr><th></th><th>Y<sub>3</sub></th><th>g<sub>0</sub>(X)</th><th>g<sub>1</sub>(X)</th><th>g<sub>2</sub>(X)</th><th>g<sub>3</sub>(X)</th></tr></thead><tbody><tr><th>Y<sub>3</sub></th><td>1.00</td><td>0.96</td><td>0.97</td><td>0.93</td><td>0.80</td></tr><tr><th>g<sub>0</sub>(X)</th><td>0.96</td><td>1.00</td><td>0.96</td><td>0.92</td><td>0.80</td></tr><tr><th>g<sub>1</sub>(X)</th><td>0.97</td><td>0.96</td><td>1.00</td><td>0.95</td><td>0.81</td></tr><tr><th>g<sub>2</sub>(X)</th><td>0.93</td><td>0.92</td><td>0.95</td><td>1.00</td><td>0.80</td></tr><tr><th>g<sub>3</sub>(X)</th><td>0.80</td><td>0.96</td><td>0.81</td><td>0.80</td><td>1.00</td></tr></tbody></table>

technique works particularly well in this case because the control variates are highly correlated to $f_3$. Indeed, table 6 shows that their Pearson coefficients are at least 0.8, which guarantees a theoretical reduction of at least 94% in the variance of the MC estimator (i.e. a reduction of at least about 76% in standard deviation), corresponding to the variance reduction when using a single control variate based on $g_3$. In fact, the variance reduction factor $R^2$ when using all the surrogates is only slightly higher, namely $R^2 \approx 95\%$ corresponding to a standard reduction factor of around 78%, which is reflected in fig. 2a.

**MLMC-MLCV vs. MLCV** Combining the MLMC and MLCV techniques allows the resulting MLMC-MLCV estimator to reduce the variance even more significantly. This can be explained by the very high correlation between $Y_0$, $Y_1$, $Y_2$ and $Y_3$ on the one hand, which ensures the good performance of the MLMC approach, and by the strong correlation between the control variates on the other hand, ensuring their good performance in combination with the MLMC technique. In
<figure><img src="image_3.png" /></figure>

(a) The evaluation cost does not include the cost of construction of the surrogate models.

(b) The evaluation cost includes the cost of construction of the surrogate models.

<caption>Figure 2: RMSE of the MC, MLMC, MLCV and MLMC-MLCV estimators of $\theta = E[f_3(\mathbf{X})]$ (see eq. (4.14)) with respect to the sampling budget $C \in \{100; 300; 1000; 3000; 10000\}$ of $f_3$ evaluations. The surrogates used for the MLCV and MLMC-MLCV estimators are described in tables 3 and 4, respectively, and their total construction cost is $C^{\text{DoE}} = 400$. The RMSEs are estimated using 500 replicates.</caption>

<table><caption>Table 6: Pearson correlation coefficients between the finest output $Y_3 = f_3(\mathbf{X})$ and the corresponding control variates based on the PC models of table 3 used for the MLCV method, estimated with a sample of size 1000.</caption><thead><tr><th></th><th>Y<sub>3</sub></th><th>g<sub>0</sub>(X)</th><th>g<sub>1</sub>(X)</th><th>g<sub>2</sub>(X)</th><th>g<sub>3</sub>(X)</th></tr></thead><tbody><tr><th>Y<sub>3</sub></th><td>1.00</td><td>0.96</td><td>0.97</td><td>0.93</td><td>0.80</td></tr><tr><th>g<sub>0</sub>(X)</th><td>0.96</td><td>1.00</td><td>0.96</td><td>0.92</td><td>0.80</td></tr><tr><th>g<sub>1</sub>(X)</th><td>0.97</td><td>0.96</td><td>1.00</td><td>0.95</td><td>0.81</td></tr><tr><th>g<sub>2</sub>(X)</th><td>0.93</td><td>0.92</td><td>0.95</td><td>1.00</td><td>0.80</td></tr><tr><th>g<sub>3</sub>(X)</th><td>0.80</td><td>0.96</td><td>0.81</td><td>0.80</td><td>1.00</td></tr></tbody></table>

technique works particularly well in this case because the control variates are highly correlated to $f_3$. Indeed, table 6 shows that their Pearson coefficients are at least 0.8, which guarantees a theoretical reduction of at least 94% in the variance of the MC estimator (i.e. a reduction of at least about 76% in standard deviation), corresponding to the variance reduction when using a single control variate based on $g_3$. In fact, the variance reduction factor $R^2$ when using all the surrogates is only slightly higher, namely $R^2 \approx 95\%$ corresponding to a standard reduction factor of around 78%, which is reflected in fig. 2a.

**MLMC-MLCV vs. MLCV** Combining the MLMC and MLCV techniques allows the resulting MLMC-MLCV estimator to reduce the variance even more significantly. This can be explained by the very high correlation between $Y_0$, $Y_1$, $Y_2$ and $Y_3$ on the one hand, which ensures the good performance of the MLMC approach, and by the strong correlation between the control variates on the other hand, ensuring their good performance in combination with the MLMC technique. In
<table><caption>Table 7: Pearson correlation coefficients between the outputs $Y_\ell = f_\ell(\mathbf{X})$ of the simulators and $g_\ell(\mathbf{X})$ of the PC models, as well as the successive differences $Y_\ell - Y_{\ell-1}$ and the outputs $h_\ell(\mathbf{X})$ of the associated PC models, for the MLMC-based estimators. The coefficients are estimated from a sample of size 1000.</caption><thead><tr><th></th><th>Y<sub>0</sub></th><th>Y<sub>1</sub></th><th>Y<sub>2</sub></th><th>Y<sub>3</sub></th><th>Y<sub>1</sub> - Y<sub>0</sub></th><th>Y<sub>2</sub> - Y<sub>1</sub></th><th>Y<sub>3</sub> - Y<sub>2</sub></th><th>g<sub>0</sub>(X)</th><th>g<sub>1</sub>(X)</th><th>g<sub>2</sub>(X)</th><th>g<sub>3</sub>(X)</th><th>h<sub>1</sub>(X)</th><th>h<sub>2</sub>(X)</th><th>h<sub>3</sub>(X)</th></tr></thead><tbody><tr><th>Y<sub>0</sub></th><td>1.00</td><td>0.97</td><td>0.97</td><td>0.97</td><td>-0.11</td><td>-0.04</td><td>-0.04</td><td>0.99</td><td>0.96</td><td>0.93</td><td>0.47</td><td>-0.12</td><td>-0.13</td><td>-0.10</td></tr><tr><th>Y<sub>1</sub></th><td>0.97</td><td>1.00</td><td>1.00</td><td>1.00</td><td>0.13</td><td>0.19</td><td>0.19</td><td>0.97</td><td>0.97</td><td>0.93</td><td>0.48</td><td>0.12</td><td>0.09</td><td>0.09</td></tr><tr><th>Y<sub>2</sub></th><td>0.97</td><td>1.00</td><td>1.00</td><td>1.00</td><td>0.14</td><td>0.20</td><td>0.20</td><td>0.96</td><td>0.97</td><td>0.92</td><td>0.48</td><td>0.13</td><td>0.10</td><td>0.10</td></tr><tr><th>Y<sub>3</sub></th><td>0.97</td><td>1.00</td><td>1.00</td><td>1.00</td><td>0.14</td><td>0.20</td><td>0.20</td><td>0.96</td><td>0.97</td><td>0.92</td><td>0.48</td><td>0.13</td><td>0.10</td><td>0.10</td></tr><tr><th>Y<sub>1</sub> - Y<sub>0</sub></th><td>-0.11</td><td>0.13</td><td>0.14</td><td>0.14</td><td>1.00</td><td>0.96</td><td>0.96</td><td>-0.10</td><td>0.07</td><td>0.01</td><td>0.03</td><td>0.99</td><td>0.90</td><td>0.80</td></tr><tr><th>Y<sub>2</sub> - Y<sub>1</sub></th><td>-0.04</td><td>0.19</td><td>0.20</td><td>0.20</td><td>0.96</td><td>1.00</td><td>1.00</td><td>-0.03</td><td>0.12</td><td>0.06</td><td>0.06</td><td>0.94</td><td>0.90</td><td>0.81</td></tr><tr><th>Y<sub>3</sub> - Y<sub>2</sub></th><td>-0.04</td><td>0.19</td><td>0.20</td><td>0.20</td><td>0.96</td><td>1.00</td><td>1.00</td><td>-0.03</td><td>0.12</td><td>0.06</td><td>0.06</td><td>0.95</td><td>0.90</td><td>0.81</td></tr><tr><th>g<sub>0</sub>(X)</th><td>0.99</td><td>0.97</td><td>0.96</td><td>0.96</td><td>-0.10</td><td>-0.03</td><td>-0.03</td><td>1.00</td><td>0.96</td><td>0.93</td><td>0.48</td><td>-0.11</td><td>-0.12</td><td>-0.09</td></tr><tr><th>g<sub>1</sub>(X)</th><td>0.96</td><td>0.97</td><td>0.97</td><td>0.97</td><td>0.07</td><td>0.12</td><td>0.12</td><td>0.96</td><td>1.00</td><td>0.94</td><td>0.47</td><td>0.06</td><td>0.04</td><td>0.04</td></tr><tr><th>g<sub>2</sub>(X)</th><td>0.93</td><td>0.93</td><td>0.92</td><td>0.92</td><td>0.01</td><td>0.06</td><td>0.06</td><td>0.93</td><td>0.94</td><td>1.00</td><td>0.47</td><td>-0.00</td><td>-0.03</td><td>-0.02</td></tr><tr><th>g<sub>3</sub>(X)</th><td>0.47</td><td>0.48</td><td>0.48</td><td>0.48</td><td>0.03</td><td>0.06</td><td>0.06</td><td>0.48</td><td>0.47</td><td>0.47</td><td>1.00</td><td>0.03</td><td>0.01</td><td>0.02</td></tr><tr><th>h<sub>1</sub>(X)</th><td>-0.12</td><td>0.12</td><td>0.13</td><td>0.13</td><td>0.99</td><td>0.94</td><td>0.95</td><td>-0.11</td><td>0.06</td><td>-0.00</td><td>0.03</td><td>1.00</td><td>0.91</td><td>0.80</td></tr><tr><th>h<sub>2</sub>(X)</th><td>-0.13</td><td>0.09</td><td>0.10</td><td>0.10</td><td>0.90</td><td>0.90</td><td>0.90</td><td>-0.12</td><td>0.04</td><td>-0.03</td><td>0.01</td><td>0.91</td><td>1.00</td><td>0.80</td></tr><tr><th>h<sub>3</sub>(X)</th><td>-0.10</td><td>0.09</td><td>0.10</td><td>0.10</td><td>0.80</td><td>0.81</td><td>0.81</td><td>-0.09</td><td>0.04</td><td>-0.02</td><td>0.02</td><td>0.80</td><td>0.80</td><td>1.00</td></tr></tbody></table>

particular, table 7 shows that $g_0(\mathbf{X})$, $g_1(\mathbf{X})$ and $g_2(\mathbf{X})$ are highly correlated with $Y_3$, with Pearson correlation coefficients greater than 0.9, while $g_3(\mathbf{X})$ is poorly correlated with $Y_3$, with a correlation coefficient of 0.48. Besides, $h_1(\mathbf{X})$, $h_2(\mathbf{X})$ and $h_3(\mathbf{X})$ are well-correlated with $Y_1 - Y_0$, $Y_2 - Y_1$ and $Y_3 - Y_2$, respectively, with correlation coefficients greater than 0.8.

**A priori estimation of the variance reduction factor and sample allocation** Further insights regarding the expected variance reduction of the MLMC-* estimators can be drawn from table 8, which reports the variance reduction factor $R_\ell^2$ w.r.t. pure MLMC on each level defined in eqs. (3.25) and (3.26), as well as the quantity $S_\ell^2$ defined in eq. (3.35). In addition, the variance of the MLMC-* estimators of the expectation with optimal sample allocation is $S_L^2/C$ (see eq. (3.36)), so that the variance reduction factor of an MLMC-* estimator compared to the high-fidelity MC estimator is $1 - S_L^2/(C_LV[Y_L])$, which is reported in the last rows and last column of the table. Consequently, the variance of the different MLMC-based estimators per unit cost can be compared directly through the ratio of their respective $S_L^2$. For instance, we observe that the variance of the MLMC-MLCV estimator is significantly reduced compared to that of the MLMC estimator, by about 1 - 40.04/2797.65 ≈ 98.6%, resulting in a reduction in standard deviation of about 88%. Again, this is well reflected in fig. 2a. Table 8 also reports the share of computational work across correction levels, namely $n_\ell^*(C_\ell+C_{\ell-1})/C$, with $n_\ell^*$ as defined in eq. (3.35), for the correction on level $\ell$. These anticipated shares are consistent with those obtained by using algorithm 1, reported in fig. 4b. It is important to note that the bottom rows of table 8 are deduced from estimates of $V[Y_\ell]$, $V_\ell$ and $R_\ell^2$, reported in the top rows. This means that, in practice, these estimates can be obtained from a few evaluations of $f_0, \dots, f_L$ (and of the appropriate surrogate models), to anticipate the expected variance reduction factor of the MLMC-* estimators and apply algorithm 1 to the most promising.

**Remark 4.1.** Table 8 highlights the fast decay of $V_\ell$ with $\ell$, which is a favorable scenario for MLMC. In fact, table 8 reports that the variance of the plain MLMC estimator is reduced by about 75% compared to that of the crude, high-fidelity MC estimator, i.e., a reduction by a factor 4.
<table><caption>Table 7: Pearson correlation coefficients between the outputs $Y_\ell = f_\ell(\mathbf{X})$ of the simulators and $g_\ell(\mathbf{X})$ of the PC models, as well as the successive differences $Y_\ell - Y_{\ell-1}$ and the outputs $h_\ell(\mathbf{X})$ of the associated PC models, for the MLMC-based estimators. The coefficients are estimated from a sample of size 1000.</caption><thead><tr><th></th><th>Y<sub>0</sub></th><th>Y<sub>1</sub></th><th>Y<sub>2</sub></th><th>Y<sub>3</sub></th><th>Y<sub>1</sub> - Y<sub>0</sub></th><th>Y<sub>2</sub> - Y<sub>1</sub></th><th>Y<sub>3</sub> - Y<sub>2</sub></th><th>g<sub>0</sub>(X)</th><th>g<sub>1</sub>(X)</th><th>g<sub>2</sub>(X)</th><th>g<sub>3</sub>(X)</th><th>h<sub>1</sub>(X)</th><th>h<sub>2</sub>(X)</th><th>h<sub>3</sub>(X)</th></tr></thead><tbody><tr><td>Y<sub>0</sub></td><td>1.00</td><td>0.97</td><td>0.97</td><td>0.97</td><td>-0.11</td><td>-0.04</td><td>-0.04</td><td>0.99</td><td>0.96</td><td>0.93</td><td>0.47</td><td>-0.12</td><td>-0.13</td><td>-0.10</td></tr><tr><td>Y<sub>1</sub></td><td>0.97</td><td>1.00</td><td>1.00</td><td>1.00</td><td>0.13</td><td>0.19</td><td>0.19</td><td>0.97</td><td>0.97</td><td>0.93</td><td>0.48</td><td>0.12</td><td>0.09</td><td>0.09</td></tr><tr><td>Y<sub>2</sub></td><td>0.97</td><td>1.00</td><td>1.00</td><td>1.00</td><td>0.14</td><td>0.20</td><td>0.20</td><td>0.96</td><td>0.97</td><td>0.92</td><td>0.48</td><td>0.13</td><td>0.10</td><td>0.10</td></tr><tr><td>Y<sub>3</sub></td><td>0.97</td><td>1.00</td><td>1.00</td><td>1.00</td><td>0.14</td><td>0.20</td><td>0.20</td><td>0.96</td><td>0.97</td><td>0.92</td><td>0.48</td><td>0.13</td><td>0.10</td><td>0.10</td></tr><tr><td>Y<sub>1</sub> - Y<sub>0</sub></td><td>-0.11</td><td>0.13</td><td>0.14</td><td>0.14</td><td>1.00</td><td>0.96</td><td>0.96</td><td>-0.10</td><td>0.07</td><td>0.01</td><td>0.03</td><td>0.99</td><td>0.90</td><td>0.80</td></tr><tr><td>Y<sub>2</sub> - Y<sub>1</sub></td><td>-0.04</td><td>0.19</td><td>0.20</td><td>0.20</td><td>0.96</td><td>1.00</td><td>1.00</td><td>-0.03</td><td>0.12</td><td>0.06</td><td>0.06</td><td>0.94</td><td>0.90</td><td>0.81</td></tr><tr><td>Y<sub>3</sub> - Y<sub>2</sub></td><td>-0.04</td><td>0.19</td><td>0.20</td><td>0.20</td><td>0.96</td><td>1.00</td><td>1.00</td><td>-0.03</td><td>0.12</td><td>0.06</td><td>0.06</td><td>0.95</td><td>0.90</td><td>0.81</td></tr><tr><td>g<sub>0</sub>(X)</td><td>0.99</td><td>0.97</td><td>0.96</td><td>0.96</td><td>-0.10</td><td>-0.03</td><td>-0.03</td><td>1.00</td><td>0.96</td><td>0.93</td><td>0.48</td><td>-0.11</td><td>-0.12</td><td>-0.09</td></tr><tr><td>g<sub>1</sub>(X)</td><td>0.96</td><td>0.97</td><td>0.97</td><td>0.97</td><td>0.07</td><td>0.12</td><td>0.12</td><td>0.96</td><td>1.00</td><td>0.94</td><td>0.47</td><td>0.06</td><td>0.04</td><td>0.04</td></tr><tr><td>g<sub>2</sub>(X)</td><td>0.93</td><td>0.93</td><td>0.92</td><td>0.92</td><td>0.01</td><td>0.06</td><td>0.06</td><td>0.93</td><td>0.94</td><td>1.00</td><td>0.47</td><td>-0.00</td><td>-0.03</td><td>-0.02</td></tr><tr><td>g<sub>3</sub>(X)</td><td>0.47</td><td>0.48</td><td>0.48</td><td>0.48</td><td>0.03</td><td>0.06</td><td>0.06</td><td>0.48</td><td>0.47</td><td>0.47</td><td>1.00</td><td>0.03</td><td>0.01</td><td>0.02</td></tr><tr><td>h<sub>1</sub>(X)</td><td>-0.12</td><td>0.12</td><td>0.13</td><td>0.13</td><td>0.99</td><td>0.94</td><td>0.95</td><td>-0.11</td><td>0.06</td><td>-0.00</td><td>0.03</td><td>1.00</td><td>0.91</td><td>0.80</td></tr><tr><td>h<sub>2</sub>(X)</td><td>-0.13</td><td>0.09</td><td>0.10</td><td>0.10</td><td>0.90</td><td>0.90</td><td>0.90</td><td>-0.12</td><td>0.04</td><td>-0.03</td><td>0.01</td><td>0.91</td><td>1.00</td><td>0.80</td></tr><tr><td>h<sub>3</sub>(X)</td><td>-0.10</td><td>0.09</td><td>0.10</td><td>0.10</td><td>0.80</td><td>0.81</td><td>0.81</td><td>-0.09</td><td>0.04</td><td>-0.02</td><td>0.02</td><td>0.80</td><td>0.80</td><td>1.00</td></tr></tbody></table>

particular, table 7 shows that $g_0(\mathbf{X})$, $g_1(\mathbf{X})$ and $g_2(\mathbf{X})$ are highly correlated with $Y_3$, with Pearson correlation coefficients greater than 0.9, while $g_3(\mathbf{X})$ is poorly correlated with $Y_3$, with a correlation coefficient of 0.48. Besides, $h_1(\mathbf{X})$, $h_2(\mathbf{X})$ and $h_3(\mathbf{X})$ are well-correlated with $Y_1 - Y_0$, $Y_2 - Y_1$ and $Y_3 - Y_2$, respectively, with correlation coefficients greater than 0.8.

**A priori estimation of the variance reduction factor and sample allocation** Further insights regarding the expected variance reduction of the MLMC-* estimators can be drawn from table 8, which reports the variance reduction factor $R_\ell^2$ w.r.t. pure MLMC on each level defined in eqs. (3.25) and (3.26), as well as the quantity $S_\ell^2$ defined in eq. (3.35). In addition, the variance of the MLMC-* estimators of the expectation with optimal sample allocation is $S_L^2/C$ (see eq. (3.36)), so that the variance reduction factor of an MLMC-* estimator compared to the high-fidelity MC estimator is $1 - S_L^2/(C_LV[Y_L])$, which is reported in the last rows and last column of the table. Consequently, the variance of the different MLMC-based estimators per unit cost can be compared directly through the ratio of their respective $S_L^2$. For instance, we observe that the variance of the MLMC-MLCV estimator is significantly reduced compared to that of the MLMC estimator, by about 1 - 40.04/2797.65 ≈ 98.6%, resulting in a reduction in standard deviation of about 88%. Again, this is well reflected in fig. 2a. Table 8 also reports the share of computational work across correction levels, namely $n_\ell^*(C_\ell+C_{\ell-1})/C$, with $n_\ell^*$ as defined in eq. (3.35), for the correction on level $\ell$. These anticipated shares are consistent with those obtained by using algorithm 1, reported in fig. 4b. It is important to note that the bottom rows of table 8 are deduced from estimates of $V[Y_\ell]$, $V_\ell$ and $R_\ell^2$, reported in the top rows. This means that, in practice, these estimates can be obtained from a few evaluations of $f_0, \dots, f_L$ (and of the appropriate surrogate models), to anticipate the expected variance reduction factor of the MLMC-* estimators and apply algorithm 1 to the most promising.

**Remark 4.1.** Table 8 highlights the fast decay of $V_\ell$ with $\ell$, which is a favorable scenario for MLMC. In fact, table 8 reports that the variance of the plain MLMC estimator is reduced by about 75% compared to that of the crude, high-fidelity MC estimator, i.e., a reduction by a factor 4.
Equivalently, this means that a 4-fold increase of the computational budget is needed to achieve the same variance with the high-fidelity MC estimator as that of the plain MLMC. This favorable setting for MLMC makes it all the more challenging to further reduce the variance. Yet, we observe that the proposed MLMC-* estimators are able to achieve a much smaller variance, especially by reducing the MLMC variance on the coarsest levels.

<table><caption>Table 8: Relevant quantities for the MLMC-based estimators, with $n_{\ell}^*$ and $S_{\ell}^2$ as defined in eq. (3.35), $C_{\ell}$ as reported in table 2, and where $V_{\ell}$ and $R_{\ell}^2$ have been estimated using an independent sample of size 10 000. The bold percentages reported in the last rows and last column correspond to the variance reduction factor of the MLMC-* estimators compared to the high-fidelity MC estimator.</caption><thead><tr><th>ℓ</th><th></th><th>0</th><th>1</th><th>2</th><th>3</th></tr></thead><tbody><tr><td></td><td>V[Y<sub>ℓ</sub>]</td><td>1.0850 × 10<sup>4</sup></td><td>1.1153 × 10<sup>4</sup></td><td>1.1205 × 10<sup>4</sup></td><td>1.1218 × 10<sup>4</sup></td></tr><tr><td></td><td>V<sub>ℓ</sub> = V[Y<sub>ℓ</sub> - Y<sub>ℓ-1</sub>]</td><td>1.0850 × 10<sup>4</sup></td><td>5.9029 × 10<sup>2</sup></td><td>1.0590</td><td>5.8160 × 10<sup>-2</sup></td></tr><tr><td rowspan="5">R<sub>ℓ</sub><sup>2</sup></td><td>MLMC</td><td>0</td><td>0</td><td>0</td><td>0</td></tr><tr><td>MLMC-CV</td><td>0.9838</td><td>0.9916</td><td>0.8224</td><td>0.6469</td></tr><tr><td>MLMC-MLCV</td><td>0.9840</td><td>0.9920</td><td>0.9173</td><td>0.9202</td></tr><tr><td>MLMC-CV[0]</td><td>0.9838</td><td>0</td><td>0</td><td>0</td></tr><tr><td>MLMC-MLCV[0]</td><td>0.9840</td><td>0.9916</td><td>0.8992</td><td>0.9027</td></tr><tr><td rowspan="5">n<sub>ℓ</sub><sup>*</sup> (C<sub>ℓ</sub> + C<sub>ℓ-1</sub>) / C</td><td>MLMC</td><td>69.63%</td><td>28.13%</td><td>1.68%</td><td>0.56%</td></tr><tr><td>MLMC-CV</td><td>71.00%</td><td>20.66%</td><td>5.69%</td><td>2.66%</td></tr><tr><td>MLMC-MLCV</td><td>73.66%</td><td>20.97%</td><td>4.05%</td><td>1.32%</td></tr><tr><td>MLMC-CV[0]</td><td>22.59%</td><td>71.69%</td><td>4.29%</td><td>1.42%</td></tr><tr><td>MLMC-MLCV[0]</td><td>72.84%</td><td>21.30%</td><td>4.42%</td><td>1.44%</td></tr><tr><td rowspan="5">S<sub>ℓ</sub><sup>2</sup></td><td>MLMC</td><td>1356.31</td><td>2673.54</td><td>2766.49</td><td>2797.65</td></tr><tr><td>MLMC-CV</td><td>21.99</td><td>36.64</td><td>41.33</td><td>43.62</td></tr><tr><td>MLMC-MLCV</td><td>21.72</td><td>35.85</td><td>38.99</td><td>40.04</td></tr><tr><td>MLMC-CV[0]</td><td>21.99</td><td>382.87</td><td>418.54</td><td>430.71</td></tr><tr><td>MLMC-MLCV[0]</td><td>21.76</td><td>36.35</td><td>39.84</td><td>41.01</td></tr><tr><td rowspan="5">1 - S<sub>ℓ</sub><sup>2</sup> / C<sub>ℓ</sub>V[Y<sub>ℓ</sub>]</td><td>MLMC</td><td>0.00%</td><td>4.11%</td><td>50.62%</td><td>75.06%</td></tr><tr><td>MLMC-CV</td><td>98.38%</td><td>98.69%</td><td>99.26%</td><td>99.61%</td></tr><tr><td>MLMC-MLCV</td><td>98.40%</td><td>98.71%</td><td>99.30%</td><td>99.64%</td></tr><tr><td>MLMC-CV[0]</td><td>98.38%</td><td>86.27%</td><td>92.53%</td><td>96.16%</td></tr><tr><td>MLMC-MLCV[0]</td><td>98.40%</td><td>98.70%</td><td>99.29%</td><td>99.63%</td></tr></tbody></table>

**Accounting for the construction cost of the surrogate model** These first results from fig. 2a highlight the interest of multilevel control variates, be it with the MLMC-MLCV estimator or simply with the MLCV one. These results suppose that the PC models are not built specifically for the study, so that the budget does not include the number of $f_3$-equivalent simulations required for their construction. Figure 2b illustrates the alternative case where the cost of the surrogate based CV estimators includes the cost of constructing the surrogates. As a result, an additional budget of $C^{DoE} = 400$ (see section 4.2), is allocated to the construction of the PC models. As was the case for the single-level PC-based CV estimators, this results in an offset of 400 in the total cost of the MLCV and MLMC-MLCV estimators. The effect is especially noticeable when the construction
Equivalently, this means that a 4-fold increase of the computational budget is needed to achieve the same variance with the high-fidelity MC estimator as that of the plain MLMC. This favorable setting for MLMC makes it all the more challenging to further reduce the variance. Yet, we observe that the proposed MLMC-* estimators are able to achieve a much smaller variance, especially by reducing the MLMC variance on the coarsest levels.

<table><caption>Table 8: Relevant quantities for the MLMC-based estimators, with $n_{\ell}^*$ and $S_{\ell}^2$ as defined in eq. (3.35), $C_{\ell}$ as reported in table 2, and where $V_{\ell}$ and $R_{\ell}^2$ have been estimated using an independent sample of size 10 000. The bold percentages reported in the last rows and last column correspond to the variance reduction factor of the MLMC-* estimators compared to the high-fidelity MC estimator.</caption><thead><tr><th>ℓ</th><th></th><th>0</th><th>1</th><th>2</th><th>3</th></tr></thead><tbody><tr><td></td><td>V[Y<sub>ℓ</sub>]</td><td>1.0850 × 10<sup>4</sup></td><td>1.1153 × 10<sup>4</sup></td><td>1.1205 × 10<sup>4</sup></td><td>1.1218 × 10<sup>4</sup></td></tr><tr><td></td><td>V<sub>ℓ</sub> = V[Y<sub>ℓ</sub> - Y<sub>ℓ-1</sub>]</td><td>1.0850 × 10<sup>4</sup></td><td>5.9029 × 10<sup>2</sup></td><td>1.0590</td><td>5.8160 × 10<sup>-2</sup></td></tr><tr><td rowspan="5">R<sub>ℓ</sub><sup>2</sup></td><td>MLMC</td><td>0</td><td>0</td><td>0</td><td>0</td></tr><tr><td>MLMC-CV</td><td>0.9838</td><td>0.9916</td><td>0.8224</td><td>0.6469</td></tr><tr><td>MLMC-MLCV</td><td>0.9840</td><td>0.9920</td><td>0.9173</td><td>0.9202</td></tr><tr><td>MLMC-CV[0]</td><td>0.9838</td><td>0</td><td>0</td><td>0</td></tr><tr><td>MLMC-MLCV[0]</td><td>0.9840</td><td>0.9916</td><td>0.8992</td><td>0.9027</td></tr><tr><td rowspan="5">n<sub>ℓ</sub><sup>*</sup> (C<sub>ℓ</sub> + C<sub>ℓ-1</sub>) / C</td><td>MLMC</td><td>69.63 %</td><td>28.13 %</td><td>1.68 %</td><td>0.56 %</td></tr><tr><td>MLMC-CV</td><td>71.00 %</td><td>20.66 %</td><td>5.69 %</td><td>2.66 %</td></tr><tr><td>MLMC-MLCV</td><td>73.66 %</td><td>20.97 %</td><td>4.05 %</td><td>1.32 %</td></tr><tr><td>MLMC-CV[0]</td><td>22.59 %</td><td>71.69 %</td><td>4.29 %</td><td>1.42 %</td></tr><tr><td>MLMC-MLCV[0]</td><td>72.84 %</td><td>21.30 %</td><td>4.42 %</td><td>1.44 %</td></tr><tr><td rowspan="5">S<sub>ℓ</sub><sup>2</sup></td><td>MLMC</td><td>1356.31</td><td>2673.54</td><td>2766.49</td><td>2797.65</td></tr><tr><td>MLMC-CV</td><td>21.99</td><td>36.64</td><td>41.33</td><td>43.62</td></tr><tr><td>MLMC-MLCV</td><td>21.72</td><td>35.85</td><td>38.99</td><td>40.04</td></tr><tr><td>MLMC-CV[0]</td><td>21.99</td><td>382.87</td><td>418.54</td><td>430.71</td></tr><tr><td>MLMC-MLCV[0]</td><td>21.76</td><td>36.35</td><td>39.84</td><td>41.01</td></tr><tr><td rowspan="5">1 - S<sub>ℓ</sub><sup>2</sup> / (C<sub>ℓ</sub>V[Y<sub>ℓ</sub>])</td><td>MLMC</td><td>0.00 %</td><td>4.11 %</td><td>50.62 %</td><td><strong>75.06 %</strong></td></tr><tr><td>MLMC-CV</td><td>98.38 %</td><td>98.69 %</td><td>99.26 %</td><td><strong>99.61 %</strong></td></tr><tr><td>MLMC-MLCV</td><td>98.40 %</td><td>98.71 %</td><td>99.30 %</td><td><strong>99.64 %</strong></td></tr><tr><td>MLMC-CV[0]</td><td>98.38 %</td><td>86.27 %</td><td>92.53 %</td><td><strong>96.16 %</strong></td></tr><tr><td>MLMC-MLCV[0]</td><td>98.40 %</td><td>98.70 %</td><td>99.29 %</td><td><strong>99.63 %</strong></td></tr></tbody></table>

**Accounting for the construction cost of the surrogate model** These first results from fig. 2a highlight the interest of multilevel control variates, be it with the MLMC-MLCV estimator or simply with the MLCV one. These results suppose that the PC models are not built specifically for the study, so that the budget does not include the number of $f_3$-equivalent simulations required for their construction. Figure 2b illustrates the alternative case where the cost of the surrogate based CV estimators includes the cost of constructing the surrogates. As a result, an additional budget of $C^{DoE} = 400$ (see section 4.2), is allocated to the construction of the PC models. As was the case for the single-level PC-based CV estimators, this results in an offset of 400 in the total cost of the MLCV and MLMC-MLCV estimators. The effect is especially noticeable when the construction
budget larger than the estimation budget, i.e. for $C \in \{100; 300\}$. For a sampling cost of 100, i.e. a total evaluation cost of 500, the MLCV estimator is still slightly more accurate than the MLMC estimator, while the MLMC-MLCV estimator still largely outperforms both.

### 4.3.3 Variants of MLMC-based CV estimators

The discussion about the surrogate construction budget prompts us to investigate variants of MLMC-MLCV using fewer surrogates, in order to reduce the total evaluation cost for limited budgets. The MLMC-CV estimator is described in section 3.2.2, while the variants MLMC-CV[0] and MLMC-MLCV[0] of MLMC-CV and MLMC-MLCV are introduced in table 1. The MLMC-CV[0] estimator only uses a control variate based on $g_0$, so that the construction cost drops to $C^{\text{DoE}} = 100$, while the MLMC-MLCV[0] estimator uses control variates based on $g_0$, $g_1$ and $h_1$, so that the construction cost drops to $C^{\text{DoE}} = 200$. The MLMC-CV estimator uses control variates based on surrogates at all levels, so that the construction cost remains $C^{\text{DoE}} = 400$.

<figure><img src="image_4.png" /><figcaption>Figure 3: RMSE of the MLMC-MLCV, MLMC-CV, MLMC-MLCV[0], and MLMC-CV[0] estimators of $\theta = \mathbb{E}[f_3(\mathbf{X})]$ (see eq. (4.14)) with respect to the sampling budget $C \in \{100; 300; 1000; 3000; 10000\}$ of $f_3$ evaluations. The surrogates used for the MLMC-based estimators are described in table 4. The RMSEs of the pure PC estimators, summarized in table 9 are shown in black for comparison. The RMSEs for the MLMC-* estimators are computed using 500 replicates, while the RMSEs of the PC-only estimators are computed from 200 replicates.</figcaption></figure>

Figure 3a shows that MLMC-CV[0] has much higher RMSE than the other variants, resulting from the fact that it only reduces the variance associated with the coarsest level of the MLMC estimator. This behavior is consistent with the quantities of table 8. In particular, the value of $S_L^2$ is about 10 times higher than for the other variants, accounting for its RMSE being about 3 times higher than for the other variants. The remaining variants have similar performances regardless of the estimation budget, which is consistent with the $S_L^2$ values given in table 8. On the other hand, when considering the construction cost of the surrogates, fig. 3b shows that MLMC-MLCV[0] performs best, as it uses only surrogate models related to the two coarsest levels, namely $g_0$, $g_1$ and $h_1$, so that the construction cost is reduced. Furthermore, these surrogates have excellent $Q^2$, and they are such that $g_0(\mathbf{X})$ is highly correlated with $Y_0$, $g_1(\mathbf{X})$ is highly correlated with $Y_1$, and $h_1(\mathbf{X})$ is highly correlated with $Y_1 - Y_0$, $Y_2 - Y_1$ and $Y_3 - Y_2$. Namely, the associated Pearson
budget larger than the estimation budget, i.e. for $C \in \{100; 300\}$. For a sampling cost of 100, i.e. a total evaluation cost of 500, the MLCV estimator is still slightly more accurate than the MLMC estimator, while the MLMC-MLCV estimator still largely outperforms both.

### 4.3.3 Variants of MLMC-based CV estimators

The discussion about the surrogate construction budget prompts us to investigate variants of MLMC-MLCV using fewer surrogates, in order to reduce the total evaluation cost for limited budgets. The MLMC-CV estimator is described in section 3.2.2, while the variants MLMC-CV[0] and MLMC-MLCV[0] of MLMC-CV and MLMC-MLCV are introduced in table 1. The MLMC-CV[0] estimator only uses a control variate based on $g_0$, so that the construction cost drops to $C^{\text{DoE}} = 100$, while the MLMC-MLCV[0] estimator uses control variates based on $g_0$, $g_1$ and $h_1$, so that the construction cost drops to $C^{\text{DoE}} = 200$. The MLMC-CV estimator uses control variates based on surrogates at all levels, so that the construction cost remains $C^{\text{DoE}} = 400$.

<figure><img src="image_4.png" /><figcaption>Figure 3: RMSE of the MLMC-MLCV, MLMC-CV, MLMC-MLCV[0], and MLMC-CV[0] estimators of $\theta = \mathbb{E}[f_3(\mathbf{X})]$ (see eq. (4.14)) with respect to the sampling budget $C \in \{100; 300; 1000; 3000; 10000\}$ of $f_3$ evaluations. The surrogates used for the MLMC-based estimators are described in table 4. The RMSEs of the pure PC estimators, summarized in table 9 are shown in black for comparison. The RMSEs for the MLMC-* estimators are computed using 500 replicates, while the RMSEs of the PC-only estimators are computed from 200 replicates.</figcaption></figure>

Figure 3a shows that MLMC-CV[0] has much higher RMSE than the other variants, resulting from the fact that it only reduces the variance associated with the coarsest level of the MLMC estimator. This behavior is consistent with the quantities of table 8. In particular, the value of $S_L^2$ is about 10 times higher than for the other variants, accounting for its RMSE being about 3 times higher than for the other variants. The remaining variants have similar performances regardless of the estimation budget, which is consistent with the $S_L^2$ values given in table 8. On the other hand, when considering the construction cost of the surrogates, fig. 3b shows that MLMC-MLCV[0] performs best, as it uses only surrogate models related to the two coarsest levels, namely $g_0$, $g_1$ and $h_1$, so that the construction cost is reduced. Furthermore, these surrogates have excellent $Q^2$, and they are such that $g_0(\mathbf{X})$ is highly correlated with $Y_0$, $g_1(\mathbf{X})$ is highly correlated with $Y_1$, and $h_1(\mathbf{X})$ is highly correlated with $Y_1 - Y_0$, $Y_2 - Y_1$ and $Y_3 - Y_2$. Namely, the associated Pearson
correlation coefficients reported in table 7 are all at least 0.94. Therefore, should one have to build the surrogates specifically for the CV estimation of a statistic, it is more advantageous to adopt the MLMC-MLCV[0] variant over the others. In our case, the construction budget is divided by two compared to MLMC-MLCV and MLMC-CV, for a similar performance in terms of RMSE. For completeness, we also compare the proposed estimators with plain PC models of the high-fidelity simulator $f_3$, whose characteristics are summarized in table 9. The black curve in fig. 3 shows that, for the same budget, combining MLMC sampling and multifidelity surrogate models is more advantageous than using a surrogate model of the highest fidelity model alone. By extrapolating, the latter should eventually become better than the MLMC-* estimators around a budget of 3000, which is much higher than practically affordable for expensive simulators. Note also that the PC models used here for comparison are optimized using the entire budget, while the MLMC-* estimators use PC models that are constructed from an arbitrary and non-optimized distribution of the sampling budget between the different model levels as reported in table 4. Strategies for optimizing this tradeoff, which constitutes a promising avenue for further improvement, will be investigated in future work.

<table><caption>Table 9: PC models of the high-fidelity simulator $f_L$ used directly for the estimation of $E[M(X)]$, with their sample size, degree and quality measure. These models are built from an optimized LHS DoE, using the basis-adaptive LARS algorithm of [7, Fig. 5]. The average values of $p^*$, $|\tilde{A}_{p^*}|$ and $Q^2$ over 200 replications is reported, along with their respective standard deviation between brackets.</caption><thead><tr><th>n<sub>L</sub><sup>DoE</sup></th><th>100</th><th>200</th><th>500</th><th>1000</th></tr></thead><tbody><tr><th>p<sup>*</sup></th><td>4(2)</td><td>6(1)</td><td>9(1)</td><td>11(1)</td></tr><tr><th>|A&#x0303;<sub>p<sup>*</sup></sub>|</th><td>20(13)</td><td>48(16)</td><td>118(31)</td><td>219(56)</td></tr><tr><th>Q<sup>2</sup></th><td>0.35(21)</td><td>0.84(2)</td><td>0.95(1)</td><td>0.98</td></tr></tbody></table>

### 4.3.4 Budget allocation

Lastly, fig. 4a shows the number $n_\ell$ of evaluations for each of the $L + 1$ correction levels of the MLMC-* telescopic sum. Precisely, $n_0$ is the number of evaluations of $f_0$, while $n_\ell$ is the number of evaluations of $f_\ell - f_{\ell-1}$, for $\ell > 0$. We observe a typical sample allocation for MLMC-like estimators in ideal cases, that is, many coarse evaluations, and fewer and fewer fine evaluations. The MLMC-CV[0] estimator slightly deviates from this pattern, with $n_1 \approx n_0$, which can be explained by the fact that $V_1 \approx 590$ is of the same order of magnitude as $(1 - R_0^2)V_0 \approx 176$. Figure 4b depicts the share of overall sampling cost associated with the different correction levels. Specifically, $n_0C_0C^{-1}$ is the share of sampling budget dedicated to evaluating $f_0$, and $n_\ell(C_\ell+C_{\ell-1})C^{-1}$ is the share dedicated to evaluating $f_\ell - f_{\ell-1}$, for $\ell > 0$. We see that, except for MLMC-CV[0] for the reasons explained above, most of the sampling budget (around 70%) is allocated to the coarsest level, and most of the remaining budget is dedicated to correction level $\ell = 1$. We note that the CV-based MLMC estimators dedicate slightly more budget to levels 2 and 3 than for the standard MLMC estimator. The sampling budget shares of fig. 4b are consistent with the theoretical optimal shared reported in table 8, including those of the MLMC-CV[0] estimator. This suggests that the sample allocation resulting from algorithm 1 seems to converge to the theoretical optimal sample allocation eq. (3.35).
correlation coefficients reported in table 7 are all at least 0.94. Therefore, should one have to build the surrogates specifically for the CV estimation of a statistic, it is more advantageous to adopt the MLMC-MLCV[0] variant over the others. In our case, the construction budget is divided by two compared to MLMC-MLCV and MLMC-CV, for a similar performance in terms of RMSE. For completeness, we also compare the proposed estimators with plain PC models of the high-fidelity simulator $f_3$, whose characteristics are summarized in table 9. The black curve in fig. 3 shows that, for the same budget, combining MLMC sampling and multifidelity surrogate models is more advantageous than using a surrogate model of the highest fidelity model alone. By extrapolating, the latter should eventually become better than the MLMC-* estimators around a budget of 3000, which is much higher than practically affordable for expensive simulators. Note also that the PC models used here for comparison are optimized using the entire budget, while the MLMC-* estimators use PC models that are constructed from an arbitrary and non-optimized distribution of the sampling budget between the different model levels as reported in table 4. Strategies for optimizing this tradeoff, which constitutes a promising avenue for further improvement, will be investigated in future work.

<table><caption>Table 9: PC models of the high-fidelity simulator $f_L$ used directly for the estimation of $E[M(X)]$, with their sample size, degree and quality measure. These models are built from an optimized LHS DoE, using the basis-adaptive LARS algorithm of [7, Fig. 5]. The average values of $p^*$, $|\tilde{A}_{p^*}|$ and $Q^2$ over 200 replications is reported, along with their respective standard deviation between brackets.</caption><thead><tr><th>n<sub>L</sub><sup>DoE</sup></th><th>100</th><th>200</th><th>500</th><th>1000</th></tr></thead><tbody><tr><th>p<sup>*</sup></th><td>4(2)</td><td>6(1)</td><td>9(1)</td><td>11(1)</td></tr><tr><th>|A&#x0303;<sub>p<sup>*</sup></sub>|</th><td>20(13)</td><td>48(16)</td><td>118(31)</td><td>219(56)</td></tr><tr><th>Q<sup>2</sup></th><td>0.35(21)</td><td>0.84(2)</td><td>0.95(1)</td><td>0.98</td></tr></tbody></table>

### 4.3.4 Budget allocation

Lastly, fig. 4a shows the number $n_\ell$ of evaluations for each of the $L + 1$ correction levels of the MLMC-* telescopic sum. Precisely, $n_0$ is the number of evaluations of $f_0$, while $n_\ell$ is the number of evaluations of $f_\ell - f_{\ell-1}$, for $\ell > 0$. We observe a typical sample allocation for MLMC-like estimators in ideal cases, that is, many coarse evaluations, and fewer and fewer fine evaluations. The MLMC-CV[0] estimator slightly deviates from this pattern, with $n_1 \approx n_0$, which can be explained by the fact that $V_1 \approx 590$ is of the same order of magnitude as $(1 - R_0^2)V_0 \approx 176$. Figure 4b depicts the share of overall sampling cost associated with the different correction levels. Specifically, $n_0C_0C^{-1}$ is the share of sampling budget dedicated to evaluating $f_0$, and $n_\ell(C_\ell+C_{\ell-1})C^{-1}$ is the share dedicated to evaluating $f_\ell - f_{\ell-1}$, for $\ell > 0$. We see that, except for MLMC-CV[0] for the reasons explained above, most of the sampling budget (around 70%) is allocated to the coarsest level, and most of the remaining budget is dedicated to correction level $\ell = 1$. We note that the CV-based MLMC estimators dedicate slightly more budget to levels 2 and 3 than for the standard MLMC estimator. The sampling budget shares of fig. 4b are consistent with the theoretical optimal shared reported in table 8, including those of the MLMC-CV[0] estimator. This suggests that the sample allocation resulting from algorithm 1 seems to converge to the theoretical optimal sample allocation eq. (3.35).
<figure><img src="image_5.png" /><figcaption>(a) Sample allocation. The bars represent the median sample size $n_\ell$ associated with the corresponding correction level $\ell$ of the MLMC-* telescopic sum. The black vertical lines represent the 25% and 75% quantiles.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>(b) Budget allocation. The stacked bars represent the proportion of the sampling cost dedicated to each level of the MLMC-* telescopic sum.</figcaption></figure>

Figure 4: Sample allocation and associated computational cost across the correction levels for the MLMC-based estimators of $\theta = E[f_3(\mathbf{X})]$ (see eq. (4.14)), with an estimation budget of $C = 10000$.

# 5 Conclusions

In this paper, we proposed multilevel variance reduction strategies relying on surrogate-based control variates. On the one hand, using specific surrogate models, such as polynomial chaos expansions or Taylor polynomial expansions, allows to directly access exact statistics (mean, variance) of the control variates. Even if these exact statistics are not directly accessible (e.g., when using GPs), they can be estimated very accurately at negligible cost. This contrasts with typical control variates relying on lower-fidelity models based on models/simulators with degraded physics or coarser discretizations, for which approximate control variate strategies need to be devised, resulting in a lower variance reduction. On the other hand, when multiple levels of fidelities (e.g., based on the discretization) with a clear cost/accuracy hierarchy are available, the surrogate-based control variate approach can be efficiently combined with multilevel strategies. The first strategy, MLCV, simply consists of using multiple control variates based on surrogate models of the simulators corresponding to the different levels. The main advantage is that the surrogate models corresponding to coarse levels may be constructed using larger sample sizes than for the finest level, resulting in a more accurate surrogate model (i.e., with lower model error). This strategy thus leads to a greater variance reduction compared to only using one surrogate model based on the finest level. This is supported by the numerical experiments we conducted, as well as by the theoretical variance reduction provided by proposition 2.1. The second strategy, MLMC-MLCV, allows to further improve the variance reduction by combining the surrogate-based control variates with an MLMC strategy. The most appropriate way to construct and utilize the surrogate models is, however, not straightforward, and was discussed in detail in sections 3.2.2 and 3.2.3. The additional variance reduction as compared to plain MLMC is demonstrated in our numerical experiments and supported by the theoretical variance reduction factor eq. (3.26) derived in section 3.2.3.

The construction cost of the surrogate models was discussed from two perspectives. When the surrogate models are constructed for the sole purpose of serving for the control variate estimation,
<figure><img src="image_5.png" /><figcaption>(a) Sample allocation. The bars represent the median sample size $n_\ell$ associated with the corresponding correction level $\ell$ of the MLMC-* telescopic sum. The black vertical lines represent the 25% and 75% quantiles.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>(b) Budget allocation. The stacked bars represent the proportion of the sampling cost dedicated to each level of the MLMC-* telescopic sum.</figcaption></figure>

Figure 4: Sample allocation and associated computational cost across the correction levels for the MLMC-based estimators of $\theta = E[f_3(\mathbf{X})]$ (see eq. (4.14)), with an estimation budget of $C = 10000$.

# 5 Conclusions

In this paper, we proposed multilevel variance reduction strategies relying on surrogate-based control variates. On the one hand, using specific surrogate models, such as polynomial chaos expansions or Taylor polynomial expansions, allows to directly access exact statistics (mean, variance) of the control variates. Even if these exact statistics are not directly accessible (e.g., when using GPs), they can be estimated very accurately at negligible cost. This contrasts with typical control variates relying on lower-fidelity models based on models/simulators with degraded physics or coarser discretizations, for which approximate control variate strategies need to be devised, resulting in a lower variance reduction. On the other hand, when multiple levels of fidelities (e.g., based on the discretization) with a clear cost/accuracy hierarchy are available, the surrogate-based control variate approach can be efficiently combined with multilevel strategies. The first strategy, MLCV, simply consists of using multiple control variates based on surrogate models of the simulators corresponding to the different levels. The main advantage is that the surrogate models corresponding to coarse levels may be constructed using larger sample sizes than for the finest level, resulting in a more accurate surrogate model (i.e., with lower model error). This strategy thus leads to a greater variance reduction compared to only using one surrogate model based on the finest level. This is supported by the numerical experiments we conducted, as well as by the theoretical variance reduction provided by proposition 2.1. The second strategy, MLMC-MLCV, allows to further improve the variance reduction by combining the surrogate-based control variates with an MLMC strategy. The most appropriate way to construct and utilize the surrogate models is, however, not straightforward, and was discussed in detail in sections 3.2.2 and 3.2.3. The additional variance reduction as compared to plain MLMC is demonstrated in our numerical experiments and supported by the theoretical variance reduction factor eq. (3.26) derived in section 3.2.3.

The construction cost of the surrogate models was discussed from two perspectives. When the surrogate models are constructed for the sole purpose of serving for the control variate estimation,