then the cost of their construction must be taken into account for fair comparison with other
approaches. In such a case, it may not be optimal, in terms of cost/accuracy tradeoff, to construct
surrogate models on all levels, especially when only a limited budget is available. In particular, using
a subset of surrogate models based on the coarser levels may already lead to considerable variance
reduction, provided that the outputs of the coarse surrogate models are sufficiently correlated with
that of the high-fidelity simulator. Then, considering additional surrogate models on finer levels
might only result in marginal improvement, at the expense of a significant computational cost. On
the contrary, if the surrogate models have already been constructed for other purposes, and are, in
some sense, available “for free,” then their construction cost need not be considered, and the entire
set of surrogate models may then be used.

From the former perspective, that is when the cost of the surrogate construction is considered
as part of the estimation cost, one may devise more involved strategies seeking to optimize the
tradeoff between the construction cost and the model error, directly impacting the projected variance
reduction. For instance, in the context of polynomial chaos surrogate models, the truncation strategy
may be controlled to this end, as proposed in a stochastic Galerkin framework in [56], where the total
polynomial degree is optimized alongside the sample size and the CV parameter to minimize the
PC-based CV estimator's variance under a cost constraint. Another avenue to improve the proposed
approach would be to replace the MLMC part of the MLMC-MLCV strategy by a more efficient
multilevel approach, such as the multilevel best linear unbiased estimator (MLBLUE) [48, 49, 47].
In particular, an MLBLUE-MLCV strategy should be more efficient when a collection of low-fidelity
simulators (e.g., with degraded physics) with no clear cost/accuracy hierarchy is available. Finally,
although the proposed approaches apply to the estimation of arbitrary statistics, they were only
tested here on the estimation of expected values. The theoretical and algorithmic ingredients for
the estimation of variances are, however, described in this paper and may be tested in follow-up
investigations and numerical experiments. Specifically, the multifidelity estimation of variance-based
sensitivity indices is of particular interest to our team.

# 6 Acknowledgment

We wish to acknowledge 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 for their support, financial funding and own knowledge.
then the cost of their construction must be taken into account for fair comparison with other approaches. In such a case, it may not be optimal, in terms of cost/accuracy tradeoff, to construct surrogate models on all levels, especially when only a limited budget is available. In particular, using a subset of surrogate models based on the coarser levels may already lead to considerable variance reduction, provided that the outputs of the coarse surrogate models are sufficiently correlated with that of the high-fidelity simulator. Then, considering additional surrogate models on finer levels might only result in marginal improvement, at the expense of a significant computational cost. On the contrary, if the surrogate models have already been constructed for other purposes, and are, in some sense, available “for free,” then their construction cost need not be considered, and the entire set of surrogate models may then be used.

From the former perspective, that is when the cost of the surrogate construction is considered as part of the estimation cost, one may devise more involved strategies seeking to optimize the tradeoff between the construction cost and the model error, directly impacting the projected variance reduction. For instance, in the context of polynomial chaos surrogate models, the truncation strategy may be controlled to this end, as proposed in a stochastic Galerkin framework in [56], where the total polynomial degree is optimized alongside the sample size and the CV parameter to minimize the PC-based CV estimator’s variance under a cost constraint. Another avenue to improve the proposed approach would be to replace the MLMC part of the MLMC-MLCV strategy by a more efficient multilevel approach, such as the multilevel best linear unbiased estimator (MLBLUE) [48, 49, 47]. In particular, an MLBLUE-MLCV strategy should be more efficient when a collection of low-fidelity simulators (e.g., with degraded physics) with no clear cost/accuracy hierarchy is available. Finally, although the proposed approaches apply to the estimation of arbitrary statistics, they were only tested here on the estimation of expected values. The theoretical and algorithmic ingredients for the estimation of variances are, however, described in this paper and may be tested in follow-up investigations and numerical experiments. Specifically, the multifidelity estimation of variance-based sensitivity indices is of particular interest to our team.

# 6 Acknowledgment

We wish to acknowledge 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 for their support, financial funding and own knowledge.
# A Optimal parameter for expectation and variance CV estimators

We derive here the expression of the optimal CV parameter $\boldsymbol{\alpha}^*$ for the CV estimators of the expectation and of the variance. We define $Y = f(\mathbf{X})$, $Z_m = g_m(\mathbf{X})$ and $\mathbf{Z} = (Z_m)_{m=1}^M$. Then, given an input $n$-sample $\{\mathbf{X}^{(i)}\}_{i=1}^n$, we define $Y^{(i)} = f(\mathbf{X}^{(i)})$, $Z_m^{(i)} = g_m(\mathbf{X}^{(i)})$ and $\mathbf{Z}^{(i)} = (Z_m^{(i)})_{m=1}^M$.

## A.1 CV estimator of the expectation

For the expectation, we have $\boldsymbol{\Sigma} = \mathbb{C}[\hat{E}[\mathbf{Z}]]$ and $\mathbf{c} = \mathbb{C}[\hat{E}[Y], \hat{E}[\mathbf{Z}]]$, i.e.

$$[\boldsymbol{\Sigma}]_{m,m'} = \mathbb{C}[\hat{E}[Z_m], \hat{E}[Z_{m'}]] = n^{-2} \sum_{i,j=1}^{n} \mathbb{C}[Z_m^{(i)}, Z_{m'}^{(j)}] = n^{-1} \mathbb{C}[Z_m, Z_{m'}], \quad (A.1)$$

$$[\mathbf{c}]_m = \mathbb{C}[\hat{E}[Y], \hat{E}[Z_m]] = n^{-2} \sum_{i,j=1}^{n} \mathbb{C}[Y^{(i)}, Z_m^{(j)}] = n^{-1} \mathbb{C}[Y, Z_m], \quad (A.2)$$

so that $\boldsymbol{\Sigma} = n^{-1}\mathbb{C}[\mathbf{Z}]$, $\mathbf{c} = n^{-1}\mathbb{C}[Y, \mathbf{Z}]$, and, eventually,

$$\boldsymbol{\alpha}^* = \mathbb{C}[\mathbf{Z}]^{-1} \mathbb{C}[Y, \mathbf{Z}]. \quad (A.3)$$

Furthermore, we have $\mathbb{V}[\hat{E}[Y]] = n^{-1}\mathbb{V}[Y]$, so that

$$R^2 = \frac{\mathbb{C}[Y, \mathbf{Z}]^\top \mathbb{C}[\mathbf{Z}]^{-1} \mathbb{C}[Y, \mathbf{Z}]}{\mathbb{V}[Y]} = \mathbf{r}_{Y,Z}^\top \mathbf{R}_{\mathbf{Z}}^{-1} \mathbf{r}_{Y,Z}, \quad (A.4)$$

where

$$\mathbf{r}_{Y,Z} = (\mathbb{V}[Y]\mathbf{D}_{\mathbf{Z}})^{-1/2}\mathbb{C}[Y, \mathbf{Z}], \quad \mathbf{R}_{\mathbf{Z}} = \mathbf{D}_{\mathbf{Z}}^{-1/2}\mathbb{C}[\mathbf{Z}]\mathbf{D}_{\mathbf{Z}}^{-1/2}, \quad \mathbf{D}_{\mathbf{Z}} = \mathrm{Diag}(\mathbb{C}[\mathbf{Z}]). \quad (A.5)$$

Thus, $R^2 \in [0, 1]$ corresponds to the squared coefficient of multiple correlation between $Y$ and the control variates $Z_1, \dots, Z_M$.

## A.2 CV estimator of the variance

Similarly, for the variance, we have $\boldsymbol{\Sigma} = \mathbb{C}[\hat{V}[\mathbf{Z}]]$ and $\mathbf{c} = \mathbb{C}[\hat{V}[Y], \hat{V}[\mathbf{Z}]]$. We start by deriving useful identities. In what follows, for any random variable $A$, we denote the corresponding centered variable by $\bar{A} := A - \mathbb{E}[A]$. First, we remark that $\hat{V}[A] = \hat{V}[\bar{A}]$, so that, for any two random variables $Y$ and $Z$,

$$\mathbb{C}[\hat{V}[Y], \hat{V}[Z]] = \mathbb{E}[\hat{V}[\bar{Y}]\hat{V}[\bar{Z}]] - \mathbb{V}[Y]\mathbb{V}[Z]. \quad (A.6)$$

Furthermore, it can be shown that

$$\mathbb{E}[\hat{V}[\bar{Y}]\hat{V}[\bar{Z}]] = \left(\frac{n}{n-1}\right)^2 (a_n(\bar{Y}, \bar{Z}) + b_n(\bar{Y}, \bar{Z}) - c_n(\bar{Y}, \bar{Z}) - c_n(\bar{Z}, \bar{Y})), \quad (A.7)$$

with (see proof below)

$$a_n(\bar{Y}, \bar{Z}) := \mathbb{E}[\hat{E}[\bar{Y}^2]\hat{E}[\bar{Z}^2]] = \frac{1}{n}\mathbb{C}[\bar{Y}^2, \bar{Z}^2] + \mathbb{V}[Y]\mathbb{V}[Z] = a_n(\bar{Z}, \bar{Y}), \quad (A.8)$$

$$b_n(\bar{Y}, \bar{Z}) := \mathbb{E}[\hat{E}[\bar{Y}^2]\hat{E}[\bar{Z}^2]] = \frac{a_n(\bar{Y}, \bar{Z})}{n^2} + 2\frac{n-1}{n^3}\mathbb{C}[Y, Z]^2 = b_n(\bar{Z}, \bar{Y}), \quad (A.9)$$

$$c_n(\bar{Y}, \bar{Z}) := \mathbb{E}[\hat{E}[\bar{Y}^2]\hat{E}[\bar{Z}^2]] = \frac{a_n(\bar{Y}, \bar{Z})}{n} = c_n(\bar{Z}, \bar{Y}). \quad (A.10)$$
# A Optimal parameter for expectation and variance CV estimators

We derive here the expression of the optimal CV parameter $\boldsymbol{\alpha}^*$ for the CV estimators of the expectation and of the variance. We define $Y = f(\mathbf{X})$, $Z_m = g_m(\mathbf{X})$ and $\mathbf{Z} = (Z_m)_{m=1}^M$. Then, given an input $n$-sample $\{\mathbf{X}^{(i)}\}_{i=1}^n$, we define $Y^{(i)} = f(\mathbf{X}^{(i)})$, $Z_m^{(i)} = g_m(\mathbf{X}^{(i)})$ and $\mathbf{Z}^{(i)} = (Z_m^{(i)})_{m=1}^M$.

## A.1 CV estimator of the expectation

For the expectation, we have $\boldsymbol{\Sigma} = \mathbb{C}[\hat{E}[\mathbf{Z}]]$ and $\mathbf{c} = \mathbb{C}[\hat{E}[Y], \hat{E}[\mathbf{Z}]]$, i.e.

$$[\boldsymbol{\Sigma}]_{m,m'} = \mathbb{C}[\hat{E}[Z_m], \hat{E}[Z_{m'}]] = n^{-2} \sum_{i,j=1}^{n} \mathbb{C}[Z_m^{(i)}, Z_{m'}^{(j)}] = n^{-1} \mathbb{C}[Z_m, Z_{m'}], \quad (A.1)$$

$$[\mathbf{c}]_m = \mathbb{C}[\hat{E}[Y], \hat{E}[Z_m]] = n^{-2} \sum_{i,j=1}^{n} \mathbb{C}[Y^{(i)}, Z_m^{(j)}] = n^{-1} \mathbb{C}[Y, Z_m], \quad (A.2)$$

so that $\boldsymbol{\Sigma} = n^{-1}\mathbb{C}[\mathbf{Z}]$, $\mathbf{c} = n^{-1}\mathbb{C}[Y, \mathbf{Z}]$, and, eventually,

$$\boldsymbol{\alpha}^* = \mathbb{C}[\mathbf{Z}]^{-1} \mathbb{C}[Y, \mathbf{Z}]. \quad (A.3)$$

Furthermore, we have $\mathbb{V}[\hat{E}[Y]] = n^{-1}\mathbb{V}[Y]$, so that

$$R^2 = \frac{\mathbb{C}[Y, \mathbf{Z}]^\top \mathbb{C}[\mathbf{Z}]^{-1} \mathbb{C}[Y, \mathbf{Z}]}{\mathbb{V}[Y]} = \mathbf{r}_{Y,Z}^\top \mathbf{R}_{\mathbf{Z}}^{-1} \mathbf{r}_{Y,Z}, \quad (A.4)$$

where

$$\mathbf{r}_{Y,Z} = (\mathbb{V}[Y]\mathbf{D}_{\mathbf{Z}})^{-1/2}\mathbb{C}[Y, \mathbf{Z}], \quad \mathbf{R}_{\mathbf{Z}} = \mathbf{D}_{\mathbf{Z}}^{-1/2}\mathbb{C}[\mathbf{Z}]\mathbf{D}_{\mathbf{Z}}^{-1/2}, \quad \mathbf{D}_{\mathbf{Z}} = \mathrm{Diag}(\mathbb{C}[\mathbf{Z}]). \quad (A.5)$$

Thus, $R^2 \in [0, 1]$ corresponds to the squared coefficient of multiple correlation between $Y$ and the control variates $Z_1, \dots, Z_M$.

## A.2 CV estimator of the variance

Similarly, for the variance, we have $\boldsymbol{\Sigma} = \mathbb{C}[\hat{V}[\mathbf{Z}]]$ and $\mathbf{c} = \mathbb{C}[\hat{V}[Y], \hat{V}[\mathbf{Z}]]$. We start by deriving useful identities. In what follows, for any random variable $A$, we denote the corresponding centered variable by $\bar{A} := A - \mathbb{E}[A]$. First, we remark that $\hat{V}[A] = \hat{V}[\bar{A}]$, so that, for any two random variables $Y$ and $Z$,

$$\mathbb{C}[\hat{V}[Y], \hat{V}[Z]] = \mathbb{E}[\hat{V}[\bar{Y}]\hat{V}[\bar{Z}]] - \mathbb{V}[Y]\mathbb{V}[Z]. \quad (A.6)$$

Furthermore, it can be shown that

$$\mathbb{E}[\hat{V}[\bar{Y}]\hat{V}[\bar{Z}]] = \left(\frac{n}{n-1}\right)^2 (a_n(\bar{Y}, \bar{Z}) + b_n(\bar{Y}, \bar{Z}) - c_n(\bar{Y}, \bar{Z}) - c_n(\bar{Z}, \bar{Y})), \quad (A.7)$$

with (see proof below)

$$a_n(\bar{Y}, \bar{Z}) := \mathbb{E}[\hat{E}[\bar{Y}^2]\hat{E}[\bar{Z}^2]] = \frac{1}{n}\mathbb{C}[\bar{Y}^2, \bar{Z}^2] + \mathbb{V}[Y]\mathbb{V}[Z] = a_n(\bar{Z}, \bar{Y}), \quad (A.8)$$

$$b_n(\bar{Y}, \bar{Z}) := \mathbb{E}[\hat{E}[\bar{Y}]^2 \hat{E}[\bar{Z}]^2] = \frac{a_n(\bar{Y}, \bar{Z})}{n^2} + 2 \frac{n-1}{n^3} \mathbb{C}[Y, Z]^2 = b_n(\bar{Z}, \bar{Y}), \quad (A.9)$$

$$c_n(\bar{Y}, \bar{Z}) := \mathbb{E}[\hat{E}[\bar{Y}^2]\hat{E}[\bar{Z}]^2] = \frac{a_n(\bar{Y}, \bar{Z})}{n} = c_n(\bar{Z}, \bar{Y}). \quad (A.10)$$
We thus have

$$ \mathbb{C}[\hat{V}[Y], \hat{V}[Z]] = \frac{1}{n} \mathbb{C}[\bar{Y}^2, \bar{Z}^2] + \frac{2}{n(n-1)} \mathbb{C}[Y, Z]^2, \quad (A.11) $$

$$ \mathbb{V}[\hat{V}[Y]] = \mathbb{C}[\hat{V}[Y], \hat{V}[Z]] = \frac{1}{n} \mathbb{V}[\bar{Y}^2] + \frac{2}{n(n-1)} \mathbb{V}[Y]^2, \quad (A.12) $$

eventually leading to

$$ \boldsymbol{\Sigma} = \mathbb{C}[\hat{V}[\mathbf{Z}]] = \frac{1}{n} \left( \mathbb{C}[\bar{\mathbf{Z}}^{\odot 2}] + \frac{2}{n-1} \mathbb{C}[\mathbf{Z}]^{\odot 2} \right), \quad (A.13) $$

$$ \mathbf{c} = \mathbb{C}[\hat{V}[Y], \hat{V}[\mathbf{Z}]] = \frac{1}{n} \left( \mathbb{C}[\bar{Y}^2, \bar{\mathbf{Z}}^{\odot 2}] + \frac{2}{n-1} \mathbb{C}[Y, \mathbf{Z}]^{\odot 2} \right), \quad (A.14) $$

so that

$$ \boldsymbol{\alpha}^* = \left[ \mathbb{C}[\bar{\mathbf{Z}}^{\odot 2}] + \frac{2}{n-1} \mathbb{C}[\mathbf{Z}]^{\odot 2} \right]^{-1} \left[ \mathbb{C}[\bar{Y}^2, \bar{\mathbf{Z}}^{\odot 2}] + \frac{2}{n-1} \mathbb{C}[Y, \mathbf{Z}]^{\odot 2} \right], \quad (A.15) $$

$$ R^2 = \frac{\left[ \mathbb{C}[\bar{Y}^2, \bar{\mathbf{Z}}^{\odot 2}] + \frac{2}{n-1} \mathbb{C}[Y, \mathbf{Z}]^{\odot 2} \right]^\top \boldsymbol{\alpha}^*}{\mathbb{V}[\bar{Y}^2] + \frac{2}{n-1} \mathbb{V}[Y]^2}. \quad (A.16) $$

As $n \to \infty$, we see that

$$ \boldsymbol{\alpha}^* \to \mathbb{C}[\bar{\mathbf{Z}}^{\odot 2}]^{-1} \mathbb{C}[\bar{Y}^2, \bar{\mathbf{Z}}^{\odot 2}] \quad (A.17) $$

$$ R^2 \to \mathbb{V}[\bar{Y}^2]^{-1} \mathbb{C}[\bar{Y}^2, \bar{\mathbf{Z}}^{\odot 2}]^\top \mathbb{C}[\bar{\mathbf{Z}}^{\odot 2}]^{-1} \mathbb{C}[\bar{Y}^2, \bar{\mathbf{Z}}^{\odot 2}] = \mathbf{r}_{Y,Z}^\top \mathbf{R}_{\bar{\mathbf{Z}}}^{-1} \mathbf{r}_{Y,Z} =: R_{\lim}^2, \quad (A.18) $$

where

$$ \mathbf{r}_{Y,Z} = (\mathbb{V}[\bar{Y}^2]\mathbf{D}_Z)^{-1/2}\mathbb{C}[\bar{Y}^2, \bar{\mathbf{Z}}^{\odot 2}], \quad \mathbf{R}_Z = \mathbf{D}_Z^{-1/2}\mathbb{C}[\bar{\mathbf{Z}}^{\odot 2}]\mathbf{D}_Z^{-1/2}, \quad (A.19) $$

with $\mathbf{D}_Z = \text{Diag}(\mathbb{C}[\bar{\mathbf{Z}}^{\odot 2}])$. Thus, $R_{\lim}^2 \in [0, 1]$ corresponds to the squared coefficient of multiple correlation between $\bar{Y}^2$ and $\bar{Z}_1^2, \dots, \bar{Z}_M^2$.

We now proceed to the proof of identities eqs. (A.8) to (A.10). First, for eq. (A.8), by definition

$$ a_n(\bar{Y}, \bar{Z}) = \frac{1}{n^2} \mathbb{E} \left[ \sum_{i=1}^n (\bar{Y}^{(i)})^2 \sum_{i=1}^n (\bar{Z}^{(i)})^2 \right] = \frac{1}{n^2} \sum_{i,j=1}^n \mathbb{E} \left[ (\bar{Y}^{(i)})^2 (\bar{Z}^{(j)})^2 \right]. \quad (A.20) $$

We distinguish two (disjoint) cases:

1. $i = j$: $\mathbb{E}[(\bar{Y}^{(i)})^2(\bar{Z}^{(j)})^2] = \mathbb{E}[(\bar{Y}^{(i)})^2(\bar{Z}^{(i)})^2] = \mathbb{E}[\bar{Y}^2\bar{Z}^2] = \mathbb{C}[\bar{Y}^2, \bar{Z}^2] + \mathbb{V}[Y]\mathbb{V}[Z]$. There are $n$ such terms in the sum.

2. $i \neq j$: $\mathbb{E}[(\bar{Y}^{(i)})^2(\bar{Z}^{(j)})^2] = \mathbb{E}[(\bar{Y}^{(i)})^2]\mathbb{E}[(\bar{Z}^{(j)})^2] = \mathbb{E}[\bar{Y}^2]\mathbb{E}[\bar{Z}^2] = \mathbb{V}[Y]\mathbb{V}[Z]$. There are $n(n-1)$ such terms in the sum.

Then eq. (A.8) follows. For eq. (A.9),

$$ b_n(\bar{Y}, \bar{Z}) = \mathbb{E} \left[ \left( \sum_{i=1}^n \bar{Y}^{(i)} \right)^2 \left( \sum_{i=1}^n \bar{Z}^{(i)} \right)^2 \right] = \frac{1}{n^4} \sum_{i,j,k,l=1}^n \mathbb{E} [\bar{Y}^{(i)} \bar{Y}^{(j)} \bar{Z}^{(k)} \bar{Z}^{(l)}]. \quad (A.21) $$

We distinguish five (disjoint) cases:
We thus have

$$ \mathbb{C}[\hat{V}[Y], \hat{V}[Z]] = \frac{1}{n} \mathbb{C}[\bar{Y}^2, \bar{Z}^2] + \frac{2}{n(n-1)} \mathbb{C}[Y, Z]^2, \quad (A.11) $$

$$ \mathbb{V}[\hat{V}[Y]] = \mathbb{C}[\hat{V}[Y], \hat{V}[Z]] = \frac{1}{n} \mathbb{V}[\bar{Y}^2] + \frac{2}{n(n-1)} \mathbb{V}[Y]^2, \quad (A.12) $$

eventually leading to

$$ \boldsymbol{\Sigma} = \mathbb{C}[\hat{V}[\mathbf{Z}]] = \frac{1}{n} \left( \mathbb{C}[\bar{\mathbf{Z}}^{\odot 2}] + \frac{2}{n-1} \mathbb{C}[\mathbf{Z}]^{\odot 2} \right), \quad (A.13) $$

$$ \mathbf{c} = \mathbb{C}[\hat{V}[Y], \hat{V}[\mathbf{Z}]] = \frac{1}{n} \left( \mathbb{C}[\bar{Y}^2, \bar{\mathbf{Z}}^{\odot 2}] + \frac{2}{n-1} \mathbb{C}[Y, \mathbf{Z}]^{\odot 2} \right), \quad (A.14) $$

so that

$$ \boldsymbol{\alpha}^* = \left[ \mathbb{C}[\bar{\mathbf{Z}}^{\odot 2}] + \frac{2}{n-1} \mathbb{C}[\mathbf{Z}]^{\odot 2} \right]^{-1} \left[ \mathbb{C}[\bar{Y}^2, \bar{\mathbf{Z}}^{\odot 2}] + \frac{2}{n-1} \mathbb{C}[Y, \mathbf{Z}]^{\odot 2} \right], \quad (A.15) $$

$$ R^2 = \frac{\left[ \mathbb{C}[\bar{Y}^2, \bar{\mathbf{Z}}^{\odot 2}] + \frac{2}{n-1} \mathbb{C}[Y, \mathbf{Z}]^{\odot 2} \right]^\top \boldsymbol{\alpha}^*}{\mathbb{V}[\bar{Y}^2] + \frac{2}{n-1} \mathbb{V}[Y]^2}. \quad (A.16) $$

As $n \to \infty$, we see that

$$ \boldsymbol{\alpha}^* \to \mathbb{C}[\bar{\mathbf{Z}}^{\odot 2}]^{-1} \mathbb{C}[\bar{Y}^2, \bar{\mathbf{Z}}^{\odot 2}] \quad (A.17) $$

$$ R^2 \to \mathbb{V}[\bar{Y}^2]^{-1} \mathbb{C}[\bar{Y}^2, \bar{\mathbf{Z}}^{\odot 2}]^\top \mathbb{C}[\bar{\mathbf{Z}}^{\odot 2}]^{-1} \mathbb{C}[\bar{Y}^2, \bar{\mathbf{Z}}^{\odot 2}] = \mathbf{r}_{Y,Z}^\top \mathbf{R}_{\bar{\mathbf{Z}}}^{-1} \mathbf{r}_{Y,Z} =: R_{\text{lim}}^2, \quad (A.18) $$

where

$$ \mathbf{r}_{Y,Z} = (\mathbb{V}[\bar{Y}^2]\mathbf{D}_Z)^{-1/2}\mathbb{C}[\bar{Y}^2, \bar{\mathbf{Z}}^{\odot 2}], \quad \mathbf{R}_Z = \mathbf{D}_Z^{-1/2}\mathbb{C}[\bar{\mathbf{Z}}^{\odot 2}]\mathbf{D}_Z^{-1/2}, \quad (A.19) $$

with $\mathbf{D}_Z = \text{Diag}(\mathbb{C}[\bar{\mathbf{Z}}^{\odot 2}])$. Thus, $R_{\text{lim}}^2 \in [0, 1]$ corresponds to the squared coefficient of multiple correlation between $\bar{Y}^2$ and $\bar{Z}_1^2, \dots, \bar{Z}_M^2$.

We now proceed to the proof of identities eqs. (A.8) to (A.10). First, for eq. (A.8), by definition

$$ a_n(\bar{Y}, \bar{Z}) = \frac{1}{n^2} \mathbb{E} \left[ \sum_{i=1}^n (\bar{Y}^{(i)})^2 \sum_{i=1}^n (\bar{Z}^{(i)})^2 \right] = \frac{1}{n^2} \sum_{i,j=1}^n \mathbb{E} \left[ (\bar{Y}^{(i)})^2 (\bar{Z}^{(j)})^2 \right]. \quad (A.20) $$

We distinguish two (disjoint) cases:

1. $i = j$: $\mathbb{E}[(\bar{Y}^{(i)})^2(\bar{Z}^{(j)})^2] = \mathbb{E}[(\bar{Y}^{(i)})^2(\bar{Z}^{(i)})^2] = \mathbb{E}[\bar{Y}^2\bar{Z}^2] = \mathbb{C}[\bar{Y}^2, \bar{Z}^2] + \mathbb{V}[Y]\mathbb{V}[Z]$. There are $n$ such terms in the sum.

2. $i \neq j$: $\mathbb{E}[(\bar{Y}^{(i)})^2(\bar{Z}^{(j)})^2] = \mathbb{E}[(\bar{Y}^{(i)})^2]\mathbb{E}[(\bar{Z}^{(j)})^2] = \mathbb{E}[\bar{Y}^2]\mathbb{E}[\bar{Z}^2] = \mathbb{V}[Y]\mathbb{V}[Z]$. There are $n(n-1)$ such terms in the sum.

Then eq. (A.8) follows. For eq. (A.9),

$$ b_n(\bar{Y}, \bar{Z}) = \mathbb{E} \left[ \left( \sum_{i=1}^n \bar{Y}^{(i)} \right)^2 \left( \sum_{i=1}^n \bar{Z}^{(i)} \right)^2 \right] = \frac{1}{n^4} \sum_{i,j,k,l=1}^n \mathbb{E} [\bar{Y}^{(i)} \bar{Y}^{(j)} \bar{Z}^{(k)} \bar{Z}^{(l)}]. \quad (A.21) $$

We distinguish five (disjoint) cases:
1. $i = j = k = \ell$: $\mathbb{E}[\bar{Y}^{(i)}\bar{Y}^{(j)}\bar{Z}^{(k)}\bar{Z}^{(\ell)}] = \mathbb{E}[(\bar{Y}^{(i)})^2(\bar{Z}^{(i)})^2] = \mathbb{E}[\bar{Y}^2\bar{Z}^2] = \mathbb{C}[\bar{Y}^2, \bar{Z}^2] + \mathbb{V}[Y]\mathbb{V}[Z]$. There are $n$ such terms in the sum.

2. $i = j \neq k = \ell$: $\mathbb{E}[\bar{Y}^{(i)}\bar{Y}^{(j)}\bar{Z}^{(k)}\bar{Z}^{(\ell)}] = \mathbb{E}[(\bar{Y}^{(i)})^2]\mathbb{E}[(\bar{Z}^{(k)})^2] = \mathbb{E}[\bar{Y}^2]\mathbb{E}[\bar{Z}^2] = \mathbb{V}[Y]\mathbb{V}[Z]$. There are $n(n-1)$ such terms in the sum.

3. $i = k \neq j = \ell$: $\mathbb{E}[\bar{Y}^{(i)}\bar{Y}^{(j)}\bar{Z}^{(k)}\bar{Z}^{(\ell)}] = \mathbb{E}[\bar{Y}^{(i)}\bar{Z}^{(i)}]\mathbb{E}[\bar{Y}^{(j)}\bar{Z}^{(j)}] = \mathbb{E}[\bar{Y}\bar{Z}]^2 = \mathbb{C}[Y, Z]^2$. There are $n(n-1)$ such terms in the sum.

4. $i = \ell \neq j = k$: $\mathbb{E}[\bar{Y}^{(i)}\bar{Y}^{(j)}\bar{Z}^{(k)}\bar{Z}^{(\ell)}] = \mathbb{E}[\bar{Y}^{(i)}\bar{Z}^{(i)}]\mathbb{E}[\bar{Y}^{(j)}\bar{Z}^{(j)}] = \mathbb{E}[\bar{Y}\bar{Z}]^2 = \mathbb{C}[Y, Z]^2$. There are $n(n-1)$ such terms in the sum.

5. All remaining cases (at least one of the indices $i, j, k, \ell$ is different from all the others): $\mathbb{E}[\bar{Y}^{(i)}\bar{Y}^{(j)}\bar{Z}^{(k)}\bar{Z}^{(\ell)}] = 0$.

Then eq. (A.9) follows. Finally, for eq. (A.10),

$$c_n(\bar{Y}, \bar{Z})\mathbb{E}[(\sum_{i=1}^n (\bar{Y}^{(i)})^2)(\sum_{i=1}^n \bar{Z}^{(i)})^2] = \frac{1}{n^3} \sum_{i,j,k=1}^n \mathbb{E}[(\bar{Y}^{(i)})^2 \bar{Z}^{(j)} \bar{Z}^{(k)}]. \quad (A.22)$$

We distinguish three (disjoint) cases:

1. $i = j = k$: $\mathbb{E}[(\bar{Y}^{(i)})^2\bar{Z}^{(j)}\bar{Z}^{(k)}] = \mathbb{E}[(\bar{Y}^{(i)})^2(\bar{Z}^{(i)})^2] = \mathbb{E}[\bar{Y}^2\bar{Z}^2] = \mathbb{M}^4[Y, Z]$. There are $n$ such terms in the sum.

2. $i \neq j = k$: $\mathbb{E}[(\bar{Y}^{(i)})^2\bar{Z}^{(j)}\bar{Z}^{(k)}] = \mathbb{E}[(\bar{Y}^{(i)})^2]\mathbb{E}[(\bar{Z}^{(j)})^2] = \mathbb{E}[\bar{Y}^2]\mathbb{E}[\bar{Z}^2] = \mathbb{V}[Y]\mathbb{V}[Z]$. There are $n(n-1)$ such terms in the sum.

3. All remaining cases (at least one of the indices $j, k$ is different from all the others): $\mathbb{E}[(\bar{Y}^{(i)})^2\bar{Z}^{(j)}\bar{Z}^{(k)}] = 0$.

Then eq. (A.10) follows.

**Remark A.1.** When the expected value $\mu_Z$ of $Z$ is known, as is the case when defining $Z$ from the prediction of certain surrogate models, such as PC expansion and Taylor polynomials (and, in some instances, GPs, see section 2.4.1), it is possible to replace $\hat{V}[Z]$ with $\hat{E}[\bar{Z}^{\odot 2}]$. The derivation of the optimal CV parameter is somewhat easier and leads to similar results as in the unknown expectation case. Specifically,

$$[\boldsymbol{\Sigma}]_{m,m'} = \mathbb{C}[\hat{E}[\bar{Z}_m^2], \hat{E}[\bar{Z}_{m'}^2]] = \mathbb{C}[n^{-1}\sum_{i=1}^n (\bar{Z}_m^{(i)})^2, n^{-1}\sum_{i=1}^n (\bar{Z}_{m'}^{(i)})^2] \quad (A.23)$$

$$= n^{-2} \sum_{i,j=1}^{n} \mathbb{C}[(\bar{Z}_m^{(i)})^2, (\bar{Z}_{m'}^{(j)})^2] = n^{-2} \sum_{i=1}^{n} \mathbb{C}[(\bar{Z}_m^{(i)})^2, (\bar{Z}_{m'}^{(i)})^2] \quad (A.24)$$

$$= n^{-1} \mathbb{C}[\bar{Z}_m^2, \bar{Z}_{m'}^2], \quad (A.25)$$

i.e. $\boldsymbol{\Sigma} = n^{-1}\mathbb{C}[\bar{Z}^{\odot 2}]$. Regarding the vector of covariances $\mathbf{c}$,

$$[\mathbf{c}]_m = \mathbb{C}[\hat{V}[\bar{Y}], \hat{E}[\bar{Z}_m^2]] = \mathbb{E}[\hat{V}[\bar{Y}]\hat{E}[\bar{Z}_m^2]] - \mathbb{V}[Y]\mathbb{V}[Z_m] \quad (A.26)$$

$$= \frac{n}{n-1}(\mathbb{E}[\hat{E}[\bar{Y}^2]\hat{E}[\bar{Z}_m^2]] - \mathbb{E}[\hat{E}[\bar{Y}^2]\hat{E}[\bar{Z}_m^2]]) - \mathbb{V}[Y]\mathbb{V}[Z_m] \quad (A.27)$$

$$= \frac{n}{n-1}(a_n(\bar{Y}, \bar{Z}_m) - c_n(\bar{Z}_m, \bar{Y})) - \mathbb{V}[Y]\mathbb{V}[Z_m] = a_n(\bar{Y}, \bar{Z}_m) - \mathbb{V}[Y]\mathbb{V}[Z_m] \quad (A.28)$$

$$= n^{-1}\mathbb{C}[\bar{Y}^2, \bar{Z}_m^2], \quad (A.29)$$
1. $i = j = k = \ell$: $\mathbb{E}[\bar{Y}^{(i)}\bar{Y}^{(j)}\bar{Z}^{(k)}\bar{Z}^{(\ell)}] = \mathbb{E}[(\bar{Y}^{(i)})^2(\bar{Z}^{(i)})^2] = \mathbb{E}[\bar{Y}^2\bar{Z}^2] = \mathbb{C}[\bar{Y}^2, \bar{Z}^2] + \mathbb{V}[Y]\mathbb{V}[Z]$. There are $n$ such terms in the sum.

2. $i = j \neq k = \ell$: $\mathbb{E}[\bar{Y}^{(i)}\bar{Y}^{(j)}\bar{Z}^{(k)}\bar{Z}^{(\ell)}] = \mathbb{E}[(\bar{Y}^{(i)})^2]\mathbb{E}[(\bar{Z}^{(k)})^2] = \mathbb{E}[\bar{Y}^2]\mathbb{E}[\bar{Z}^2] = \mathbb{V}[Y]\mathbb{V}[Z]$. There are $n(n-1)$ such terms in the sum.

3. $i = k \neq j = \ell$: $\mathbb{E}[\bar{Y}^{(i)}\bar{Y}^{(j)}\bar{Z}^{(k)}\bar{Z}^{(\ell)}] = \mathbb{E}[\bar{Y}^{(i)}\bar{Z}^{(i)}]\mathbb{E}[\bar{Y}^{(j)}\bar{Z}^{(j)}] = \mathbb{E}[\bar{Y}\bar{Z}]^2 = \mathbb{C}[Y, Z]^2$. There are $n(n-1)$ such terms in the sum.

4. $i = \ell \neq j = k$: $\mathbb{E}[\bar{Y}^{(i)}\bar{Y}^{(j)}\bar{Z}^{(k)}\bar{Z}^{(\ell)}] = \mathbb{E}[\bar{Y}^{(i)}\bar{Z}^{(i)}]\mathbb{E}[\bar{Y}^{(j)}\bar{Z}^{(j)}] = \mathbb{E}[\bar{Y}\bar{Z}]^2 = \mathbb{C}[Y, Z]^2$. There are $n(n-1)$ such terms in the sum.

5. All remaining cases (at least one of the indices $i, j, k, \ell$ is different from all the others): $\mathbb{E}[\bar{Y}^{(i)}\bar{Y}^{(j)}\bar{Z}^{(k)}\bar{Z}^{(\ell)}] = 0$.

Then eq. (A.9) follows. Finally, for eq. (A.10),

$$c_n(\bar{Y}, \bar{Z})\mathbb{E}[(\sum_{i=1}^n (\bar{Y}^{(i)})^2)(\sum_{i=1}^n \bar{Z}^{(i)})^2] = \frac{1}{n^3} \sum_{i,j,k=1}^n \mathbb{E}[(\bar{Y}^{(i)})^2 \bar{Z}^{(j)} \bar{Z}^{(k)}]. \quad (A.22)$$

We distinguish three (disjoint) cases:

1. $i = j = k$: $\mathbb{E}[(\bar{Y}^{(i)})^2\bar{Z}^{(j)}\bar{Z}^{(k)}] = \mathbb{E}[(\bar{Y}^{(i)})^2(\bar{Z}^{(i)})^2] = \mathbb{E}[\bar{Y}^2\bar{Z}^2] = \mathbb{M}^4[Y, Z]$. There are $n$ such terms in the sum.

2. $i \neq j = k$: $\mathbb{E}[(\bar{Y}^{(i)})^2\bar{Z}^{(j)}\bar{Z}^{(k)}] = \mathbb{E}[(\bar{Y}^{(i)})^2]\mathbb{E}[(\bar{Z}^{(j)})^2] = \mathbb{E}[\bar{Y}^2]\mathbb{E}[\bar{Z}^2] = \mathbb{V}[Y]\mathbb{V}[Z]$. There are $n(n-1)$ such terms in the sum.

3. All remaining cases (at least one of the indices $j, k$ is different from all the others): $\mathbb{E}[(\bar{Y}^{(i)})^2\bar{Z}^{(j)}\bar{Z}^{(k)}] = 0$.

Then eq. (A.10) follows.

**Remark A.1.** When the expected value $\mu_Z$ of $Z$ is known, as is the case when defining $Z$ from the prediction of certain surrogate models, such as PC expansion and Taylor polynomials (and, in some instances, GPs, see section 2.4.1), it is possible to replace $\hat{V}[Z]$ with $\hat{E}[\bar{Z}^{\odot 2}]$. The derivation of the optimal CV parameter is somewhat easier and leads to similar results as in the unknown expectation case. Specifically,

$$[\boldsymbol{\Sigma}]_{m,m'} = \mathbb{C}[\hat{E}[\bar{Z}_m^2], \hat{E}[\bar{Z}_{m'}^2]] = \mathbb{C}[n^{-1} \sum_{i=1}^n (\bar{Z}_m^{(i)})^2, n^{-1} \sum_{i=1}^n (\bar{Z}_{m'}^{(i)})^2] \quad (A.23)$$

$$= n^{-2} \sum_{i,j=1}^{n} \mathbb{C}[(\bar{Z}_m^{(i)})^2, (\bar{Z}_{m'}^{(j)})^2] = n^{-2} \sum_{i=1}^{n} \mathbb{C}[(\bar{Z}_m^{(i)})^2, (\bar{Z}_{m'}^{(i)})^2] \quad (A.24)$$

$$= n^{-1} \mathbb{C}[\bar{Z}_m^2, \bar{Z}_{m'}^2], \quad (A.25)$$

i.e. $\boldsymbol{\Sigma} = n^{-1}\mathbb{C}[\bar{\boldsymbol{Z}}^{\odot 2}]$. Regarding the vector of covariances $\mathbf{c}$,

$$[\mathbf{c}]_m = \mathbb{C}[\hat{V}[\bar{Y}], \hat{E}[\bar{Z}_m^2]] = \mathbb{E}[\hat{V}[\bar{Y}]\hat{E}[\bar{Z}_m^2]] - \mathbb{V}[Y]\mathbb{V}[Z_m] \quad (A.26)$$

$$= \frac{n}{n-1} (\mathbb{E}[\hat{E}[\bar{Y}^2]\hat{E}[\bar{Z}_m^2]] - \mathbb{E}[\hat{E}[\bar{Y}^2]\hat{E}[\bar{Z}_m^2]]) - \mathbb{V}[Y]\mathbb{V}[Z_m] \quad (A.27)$$

$$= \frac{n}{n-1}(a_n(\bar{Y}, \bar{Z}_m) - c_n(\bar{Z}_m, \bar{Y})) - \mathbb{V}[Y]\mathbb{V}[Z_m] = a_n(\bar{Y}, \bar{Z}_m) - \mathbb{V}[Y]\mathbb{V}[Z_m] \quad (A.28)$$

$$= n^{-1}\mathbb{C}[\bar{Y}^2, \bar{Z}_m^2], \quad (A.29)$$
i.e. $\mathbf{c} = n^{-1}\mathbb{C}[\bar{Y}^2, \bar{Z}^{\odot 2}]$, so that

$$ \boldsymbol{\alpha}^* = \mathbb{C}[\bar{\mathbf{Z}}^{\odot 2}]^{-1} \mathbb{C}[\bar{Y}^2, \bar{\mathbf{Z}}^{\odot 2}], \qquad (A.30) $$

$$ R^2 = \frac{\mathbb{C}[\bar{Y}^2, \bar{\mathbf{Z}}^{\odot 2}]^\top \mathbb{C}[\bar{\mathbf{Z}}^{\odot 2}]^{-1} \mathbb{C}[\bar{Y}^2, \bar{\mathbf{Z}}^{\odot 2}]}{\mathbb{V}[\bar{Y}^2] + \frac{2}{n-1} \mathbb{V}[Y]^2} \xrightarrow{n \to \infty} R_{\text{lim}}^2, \qquad (A.31) $$

with the same definitions as in eq. (A.18).

# B Variance reduction of bi-fidelity (A)CV estimators of the expectation

In this appendix, we illustrate the superiority of the (exact) CV estimator of the expectation over the bi-fidelity ACV estimator introduced in [39, section 3] (see also [19, section IV]), hereafter referred to as the MFMC estimator, in terms of variance reduction. We consider here a high-fidelity simulator $f$ and a low-fidelity version $g$. We denote by $\rho$ the Pearson correlation coefficient between $f(\mathbf{X})$ and $g(\mathbf{X})$, and we define $w := c_g/c_f$ as the ratio of the expected cost $c_g$ of evaluating $g(\mathbf{X})$ to the expected cost $c_f$ of evaluating $f(\mathbf{X})$. Note that our definition of $w$ is the inverse of the definition of $w$ in [39, 19], allowing our $w$ to remain finite when the expected cost of evaluating $g(\mathbf{X})$ vanishes. We shall assume hereafter that $w, \rho^2 \in (0, 1)$.

## B.1 The MFMC estimator

For a given (expected) computational budget $\mathcal{C}$, the optimal MFMC estimator [39, 19] reads

$$ \hat{\mu}_{\mathcal{C}}^{\text{MFMC}} = \frac{1}{N} \sum_{i=1}^{N} f(\mathbf{X}^{(i)}) - \alpha \left( \frac{1}{N} \sum_{i=1}^{N} g(\mathbf{X}^{(i)}) - \frac{1}{N'} \sum_{i=1}^{N'} g(\mathbf{X}^{(i)}) \right), \qquad (B.1) $$

where

$$ \alpha = \frac{\mathbb{C}[Y, Z]}{\mathbb{V}[Z]}, \quad N = \frac{\mathcal{C}}{c_f(1 + \eta w)}, \quad \eta = \sqrt{\frac{\rho^2}{w(1 - \rho^2)}}, \quad N' = \eta N, \qquad (B.2) $$

and where $\mathcal{X} := \{\mathbf{X}^{(1)}, \dots, \mathbf{X}^{(N')}\}$ is an $N'$-sample of $\mathbf{X}$ (we assume here that $N' > N$; see [39] and corollary B.6 for a discussion on this assumption). Note that $N$ and $N'$ need to be appropriately rounded to integers for the estimator eq. (B.1) to make sense. For the same (expected) budget $\mathcal{C}$, one can afford $p := \mathcal{C}/c_f$ evaluations of $f(\mathbf{X})$, i.e., the expected cost of $p$ evaluations of $f(\mathbf{X})$ is $\mathcal{C}$. Let $\hat{\mu}_p := p^{-1} \sum_{i=1}^p f(\mathbf{X}^{(i)})$ denote the standard MC estimator of the expected value of $f(\mathbf{X})$ using a $p$-sample of $\mathbf{X}$. We now proceed to proving some useful properties of the MFMC estimator eq. (B.1).

**Proposition B.1.** The variance of the MFMC estimator defined by eqs. (B.1) and (B.2) is such that $\mathbb{V}[\hat{\mu}_{\mathcal{C}}^{\text{MFMC}}] = \beta^{\text{MFMC}}\mathbb{V}[\hat{\mu}_p]$, with $\beta^{\text{MFMC}} = (1 + \eta w)^2(1 - \rho^2)$.

*Proof.* From [39, Eq. (6)] (see also [19, Eq. (9)]), we have

$$ \beta^{\text{MFMC}} = (1 + \eta w) \left[ 1 - \left( 1 - \frac{1}{\eta} \right) \rho^2 \right]. \qquad (B.3) $$
i.e. $\mathbf{c} = n^{-1}\mathbb{C}[\bar{Y}^2, \bar{Z}^{\odot 2}]$, so that

$$ \alpha^* = \mathbb{C}[\bar{Z}^{\odot 2}]^{-1}\mathbb{C}[\bar{Y}^2, \bar{Z}^{\odot 2}], \qquad (A.30) $$

$$ R^2 = \frac{\mathbb{C}[\bar{Y}^2, \bar{Z}^{\odot 2}]^\top \mathbb{C}[\bar{Z}^{\odot 2}]^{-1} \mathbb{C}[\bar{Y}^2, \bar{Z}^{\odot 2}]}{\mathbb{V}[\bar{Y}^2] + \frac{2}{n-1} \mathbb{V}[Y]^2} \xrightarrow{n \to \infty} R_{\text{lim}}^2, \qquad (A.31) $$

with the same definitions as in eq. (A.18).

# B Variance reduction of bi-fidelity (A)CV estimators of the expectation

In this appendix, we illustrate the superiority of the (exact) CV estimator of the expectation over the bi-fidelity ACV estimator introduced in [39, section 3] (see also [19, section IV]), hereafter referred to as the MFMC estimator, in terms of variance reduction. We consider here a high-fidelity simulator $f$ and a low-fidelity version $g$. We denote by $\rho$ the Pearson correlation coefficient between $f(\mathbf{X})$ and $g(\mathbf{X})$, and we define $w := c_g/c_f$ as the ratio of the expected cost $c_g$ of evaluating $g(\mathbf{X})$ to the expected cost $c_f$ of evaluating $f(\mathbf{X})$. Note that our definition of $w$ is the inverse of the definition of $w$ in [39, 19], allowing our $w$ to remain finite when the expected cost of evaluating $g(\mathbf{X})$ vanishes. We shall assume hereafter that $w, \rho^2 \in (0, 1)$.

## B.1 The MFMC estimator

For a given (expected) computational budget $\mathcal{C}$, the optimal MFMC estimator [39, 19] reads

$$ \hat{\mu}_{\mathcal{C}}^{\text{MFMC}} = \frac{1}{N} \sum_{i=1}^{N} f(\mathbf{X}^{(i)}) - \alpha \left( \frac{1}{N} \sum_{i=1}^{N} g(\mathbf{X}^{(i)}) - \frac{1}{N'} \sum_{i=1}^{N'} g(\mathbf{X}^{(i)}) \right), \qquad (B.1) $$

where

$$ \alpha = \frac{\mathbb{C}[Y, Z]}{\mathbb{V}[Z]}, \quad N = \frac{\mathcal{C}}{c_f(1 + \eta w)}, \quad \eta = \sqrt{\frac{\rho^2}{w(1 - \rho^2)}}, \quad N' = \eta N, \qquad (B.2) $$

and where $\mathcal{X} := \{\mathbf{X}^{(1)}, \dots, \mathbf{X}^{(N')}\}$ is an $N'$-sample of $\mathbf{X}$ (we assume here that $N' > N$; see [39] and corollary B.6 for a discussion on this assumption). Note that $N$ and $N'$ need to be appropriately rounded to integers for the estimator eq. (B.1) to make sense. For the same (expected) budget $\mathcal{C}$, one can afford $p := \mathcal{C}/c_f$ evaluations of $f(\mathbf{X})$, i.e., the expected cost of $p$ evaluations of $f(\mathbf{X})$ is $\mathcal{C}$. Let $\hat{\mu}_p := p^{-1} \sum_{i=1}^p f(\mathbf{X}^{(i)})$ denote the standard MC estimator of the expected value of $f(\mathbf{X})$ using a $p$-sample of $\mathbf{X}$. We now proceed to proving some useful properties of the MFMC estimator eq. (B.1).

**Proposition B.1.** The variance of the MFMC estimator defined by eqs. (B.1) and (B.2) is such that $\mathbb{V}[\hat{\mu}_{\mathcal{C}}^{\text{MFMC}}] = \beta^{\text{MFMC}}\mathbb{V}[\hat{\mu}_p]$, with $\beta^{\text{MFMC}} = (1 + \eta w)^2(1 - \rho^2)$.

*Proof.* From [39, Eq. (6)] (see also [19, Eq. (9)]), we have

$$ \beta^{\text{MFMC}} = (1 + \eta w) \left[ 1 - \left( 1 - \frac{1}{\eta} \right) \rho^2 \right]. \qquad (B.3) $$
Furthermore,

$$1 - \left(1 - \frac{1}{\eta}\right) \rho^2 = 1 - \rho^2 + \frac{\rho^2}{\eta} = (1 - \rho^2) \left(1 + \frac{1}{\eta} \frac{\rho^2}{1 - \rho^2}\right). \qquad (B.4)$$

Finally, injecting the expression of $\eta$ defined in eq. (B.2), we notice that

$$\frac{1}{\eta} \frac{\rho^2}{1 - \rho^2} = \sqrt{\frac{w(1 - \rho^2)}{\rho^2}} \frac{\rho^2}{1 - \rho^2} = \sqrt{\frac{w\rho^2}{1 - \rho^2}} = \eta w, \qquad (B.5)$$

which concludes the proof. □

We note that $\beta^{\text{MFMC}}$ does not depend on $C$, and that the smaller $\beta^{\text{MFMC}}$, the greater the variance reduction with respect to the standard, high-fidelity MC estimator for the same (expected) budget. We also note that variance reduction is achieved if and only if $\beta^{\text{MFMC}} < 1$.

**Proposition B.2.** $\eta > 1$ (i.e., $N' > N$) if and only if $\rho^2 > w/(1+w)$.

$$\text{Proof.} \quad \eta > 1 \iff \eta^2 > 1 \iff \rho^2 > w(1 - \rho^2) \iff (1 + w)\rho^2 > w, \text{ hence the result follows.} \quad \square$$

**Proposition B.3.** $\beta^{\text{MFMC}} < 1$ if and only if $\rho^2 > 4w/(1+w)^2$.

*Proof.* We start by noticing that $1/(1-\rho^2) = 1 + \rho^2/(1-\rho^2) = 1 + \eta^2 w$. Then, it follows that $\beta^{\text{MFMC}} < 1 \iff (1+\eta w)^2 < 1/(1-\rho^2) = 1 + \eta^2 w \iff [2-(1-w)\eta]\eta w < 0$. Because $w, \eta > 0$ and $(1-w) \in (0,1)$, we deduce that $\beta^{\text{MFMC}} < 1 \iff \eta > 2/(1-w) \iff \eta^2 > 4/(1-w)^2$. Injecting the expression of $\eta$ defined in eq. (B.2), we get $\beta^{\text{MFMC}} < 1 \iff \rho^2(1-w)^2 > 4w(1-\rho^2)$, eventually leading to the desired result. □

**Proposition B.4.** For any $w \in (0,1)$, $4w/(1+w)^2 < 1$.

*Proof.* The function $a: x \mapsto 4x/(1+x)^2$ is continuous on $[0,1]$, and its derivative is defined on $[0,1]$ by $a'(x) = 4(1-x)/(1+x)^3$, which is positive for $x \in (0,1)$, so $a$ is strictly increasing on $(0,1)$. Furthermore, for $x \in [0,1]$, $a(x) = 1 \iff x = 1$, implying that, for any $x \in (0,1)$, $a(x) < 1$, which concludes the proof. □

proposition B.4 implies that, for any $w \in (0,1)$, variance reduction can always be achieved provided that $f(\mathbf{X})$ and $g(\mathbf{X})$ are sufficiently correlated.

**Proposition B.5.** For any $w \in (0,1)$, $4w/(1+w)^2 > 2w/(1+w)$.

*Proof.* Taking the ratio of the left-hand side to the right-hand side, we have

$$\frac{4w/(1+w)^2}{2w/(1+w)} = \frac{2}{1+w} > 1, \qquad (B.6)$$

since $w < 1$. □

**Corollary B.6.** If $\beta^{\text{MFMC}} < 1$, then $N' > N$.

*Proof.* From proposition B.3, if $\beta^{\text{MFMC}} < 1$, then $\rho^2 > 4w/(1+w)^2$. Then, from proposition B.5, we deduce that $\rho^2 > 2w/(1+w) > w/(1+w)$. Finally, we conclude from proposition B.2 that $\eta > 1$, i.e., $N' > N$. □
Furthermore,

$$1 - \left(1 - \frac{1}{\eta}\right) \rho^2 = 1 - \rho^2 + \frac{\rho^2}{\eta} = (1 - \rho^2) \left(1 + \frac{\rho^2}{\eta(1 - \rho^2)}\right). \quad (B.4)$$

Finally, injecting the expression of $\eta$ defined in eq. (B.2), we notice that

$$\frac{1}{\eta} \frac{\rho^2}{1 - \rho^2} = \sqrt{\frac{w(1 - \rho^2)}{\rho^2}} \frac{\rho^2}{1 - \rho^2} = \sqrt{\frac{w\rho^2}{1 - \rho^2}} = \eta w, \quad (B.5)$$

which concludes the proof. □

We note that $\beta^{\text{MFMC}}$ does not depend on $C$, and that the smaller $\beta^{\text{MFMC}}$, the greater the variance reduction with respect to the standard, high-fidelity MC estimator for the same (expected) budget. We also note that variance reduction is achieved if and only if $\beta^{\text{MFMC}} < 1$.

**Proposition B.2.** $\eta > 1$ (i.e., $N' > N$) if and only if $\rho^2 > w/(1+w)$.

$$\text{Proof.} \quad \eta > 1 \iff \eta^2 > 1 \iff \rho^2 > w(1 - \rho^2) \iff (1 + w)\rho^2 > w, \text{ hence the result follows.} \quad \square$$

**Proposition B.3.** $\beta^{\text{MFMC}} < 1$ if and only if $\rho^2 > 4w/(1+w)^2$.

*Proof.* We start by noticing that $1/(1-\rho^2) = 1 + \rho^2/(1-\rho^2) = 1 + \eta^2 w$. Then, it follows that $\beta^{\text{MFMC}} < 1 \iff (1+\eta w)^2 < 1/(1-\rho^2) = 1 + \eta^2 w \iff [2-(1-w)\eta]\eta w < 0$. Because $w, \eta > 0$ and $(1-w) \in (0,1)$, we deduce that $\beta^{\text{MFMC}} < 1 \iff \eta > 2/(1-w) \iff \eta^2 > 4/(1-w)^2$. Injecting the expression of $\eta$ defined in eq. (B.2), we get $\beta^{\text{MFMC}} < 1 \iff \rho^2(1-w)^2 > 4w(1-\rho^2)$, eventually leading to the desired result. □

**Proposition B.4.** For any $w \in (0,1)$, $4w/(1+w)^2 < 1$.

*Proof.* The function $a: x \mapsto 4x/(1+x)^2$ is continuous on $[0,1]$, and its derivative is defined on $[0,1]$ by $a'(x) = 4(1-x)/(1+x)^3$, which is positive for $x \in (0,1)$, so $a$ is strictly increasing on $(0,1)$. Furthermore, for $x \in [0,1]$, $a(x) = 1 \iff x = 1$, implying that, for any $x \in (0,1)$, $a(x) < 1$, which concludes the proof. □

proposition B.4 implies that, for any $w \in (0,1)$, variance reduction can always be achieved provided that $f(\mathbf{X})$ and $g(\mathbf{X})$ are sufficiently correlated.

**Proposition B.5.** For any $w \in (0,1)$, $4w/(1+w)^2 > 2w/(1+w)$.

*Proof.* Taking the ratio of the left-hand side to the right-hand side, we have

$$\frac{4w/(1+w)^2}{2w/(1+w)} = \frac{2}{1+w} > 1, \quad (B.6)$$

since $w < 1$. □

**Corollary B.6.** If $\beta^{\text{MFMC}} < 1$, then $N' > N$.

*Proof.* From proposition B.3, if $\beta^{\text{MFMC}} < 1$, then $\rho^2 > 4w/(1+w)^2$. Then, from proposition B.5, we deduce that $\rho^2 > 2w/(1+w) > w/(1+w)$. Finally, we conclude from proposition B.2 that $\eta > 1$, i.e., $N' > N$. □
## B.2 The exact CV estimator

We now turn to the exact CV estimator of the expected value of $f(\mathbf{X})$, i.e., considering that $\tau := \mathbb{E}[g(\mathbf{X})]$ is known. For a given (expected) computational budget $\mathcal{C}$, the (optimal) CV estimator reads

$$ \hat{\mu}_{\mathcal{C}}^{\text{CV}} = \frac{1}{n} \sum_{i=1}^{n} f(\mathbf{X}^{(i)}) - \alpha \left( \frac{1}{n} \sum_{i=1}^{n} g(\mathbf{X}^{(i)}) - \tau \right), \quad \text{with } \alpha = \frac{\mathbb{C}[Y, Z]}{\mathbb{V}[Z]}, \quad n = \frac{\mathcal{C}}{c_f(1+w)}, \quad (\text{B.7}) $$

and where $\mathcal{X} := \{\mathbf{X}^{(1)}, \dots, \mathbf{X}^{(n)}\}$ is an $n$-sample of $\mathbf{X}$. We have the following properties.

**Proposition B.7.** The variance of the CV estimator eq. (B.7) is such that $\mathbb{V}[\hat{\mu}_{\mathcal{C}}^{\text{CV}}] = \beta^{\text{CV}}\mathbb{V}[\hat{\mu}_{p}]$, with $\beta^{\text{CV}} = (1+w)(1-\rho^2)$.

*Proof.* From eq. (2.5) and appendix A.1, we know that $\mathbb{V}[\hat{\mu}_{\mathcal{C}}^{\text{CV}}] = (1-\rho^2)\mathbb{V}[\hat{\mu}_{n}]$. The result follows from the fact that $\mathbb{V}[\hat{\mu}_{n}] = \mathbb{V}[f(\mathbf{X})]/n$ and $\mathbb{V}[\hat{\mu}_{p}] = \mathbb{V}[f(\mathbf{X})]/p$, with $n = p/(1+w)$. $\square$

**Proposition B.8.** $\beta^{\text{CV}} < 1$ if and only if $\rho^2 > w/(1+w)$.

*Proof.* $\beta^{\text{CV}} < 1 \iff (1+w)(\rho^2-1) > -1 \iff \rho^2 > 1 - 1/(1+w) = w/(1+w)$. $\square$

Combined with propositions B.3 and B.8, proposition B.5 indicates that the minimum value of $\rho^2$ that guarantees variance reduction is always (i.e., for any $w \in (0,1)$) at least twice as large for the MFMC estimator than for the CV estimator. For instance, for $w = 1/3$, the MFMC estimator requires $\rho^2 > 3/4$ to have a reduced variance, while the CV estimator only requires $\rho^2 > 1/4$.

**Proposition B.9.** For $w, \rho^2 \in (0,1)$, $\beta^{\text{CV}} < 1$ if and only if $(1+w)\beta^{\text{CV}} < \beta^{\text{MFMC}}$.

*Proof.* $\beta^{\text{CV}} < 1 \iff \rho^2 > w/(1+w) \iff 1-\rho^2 < 1/(1+w) \iff 1/(1-\rho^2) > 1+w$. Because $1/(1-\rho^2) - 1 = \rho^2/(1-\rho^2)$, $\beta^{\text{CV}} < 1 \iff \rho^2/(1-\rho^2) > w \iff w\rho^2/(1-\rho^2) > w^2$. We note that $w\rho^2/(1-\rho^2) = (\eta w)^2$, so that $\beta^{\text{CV}} < 1 \iff \eta w > w$. Therefore, because $\beta^{\text{MFMC}}/\beta^{\text{CV}} = (1+\eta w)^2/(1+w)$, $\beta^{\text{CV}} < 1 \iff \beta^{\text{MFMC}}/\beta^{\text{CV}} > (1+w)^2/(1+w) = 1+w$, which concludes the proof. $\square$

**Corollary B.10.** For $w, \rho^2 \in (0,1)$, if $\beta^{\text{CV}} < 1$, then $\beta^{\text{CV}} < \beta^{\text{MFMC}}$.

*Proof.* Because $w > 0$, $\beta^{\text{MFMC}}/\beta^{\text{CV}} > 1+w \implies \beta^{\text{MFMC}}/\beta^{\text{CV}} > w$. Hence the result follows from proposition B.9. $\square$

Corollary B.10 implies that, in the $(w, \rho^2)$ variance reduction region of the CV estimator, which is larger than that of the MFMC estimator by propositions B.3, B.5, and B.8, the variance reduction of the CV estimator is always greater than that of the MFMC estimator (i.e., $\beta^{\text{CV}} < \beta^{\text{MFMC}}$).

Figure 5 depicts the $(w, \rho^2)$ variance reduction regions for the CV (fig. 5a) and MFMC (fig. 5b) estimators. We clearly see that for any $w \in (0,1)$, a higher value of $\rho^2$ is required for the MFMC estimator to achieve variance reduction than for the CV estimator (red solid line). Figure 5 also illustrates an immediate consequence of propositions B.1 and B.7, namely that $\lim_{w \to 0} \beta^{\text{CV}} = \lim_{w \to 0} \beta^{\text{MFMC}} = 1-\rho^2$. We also observe that $\beta^{\text{CV}}$ converges to $1-\rho^2$ faster than $\beta^{\text{MFMC}}$, consistent with corollary B.10.

These properties illustrate the advantages of using surrogate-based control variates. Indeed, when the low-fidelity model $g$ corresponds to a surrogate model, not only $w \ll 1$, but the exact expectation $\mathbb{E}[g(\mathbf{X})]$ may be known, so that the exact CV approach can be used, resulting in a variance reduction factor close to $\rho^2$.
## B.2 The exact CV estimator

We now turn to the exact CV estimator of the expected value of $f(\mathbf{X})$, i.e., considering that $\tau := \mathbb{E}[g(\mathbf{X})]$ is known. For a given (expected) computational budget $\mathcal{C}$, the (optimal) CV estimator reads

$$ \hat{\mu}_{\mathcal{C}}^{\text{CV}} = \frac{1}{n} \sum_{i=1}^{n} f(\mathbf{X}^{(i)}) - \alpha \left( \frac{1}{n} \sum_{i=1}^{n} g(\mathbf{X}^{(i)}) - \tau \right), \quad \text{with } \alpha = \frac{\mathbb{C}[Y, Z]}{\mathbb{V}[Z]}, \quad n = \frac{\mathcal{C}}{c_f(1+w)}, \quad (\text{B.7}) $$

and where $\mathcal{X} := \{\mathbf{X}^{(1)}, \dots, \mathbf{X}^{(n)}\}$ is an $n$-sample of $\mathbf{X}$. We have the following properties.

**Proposition B.7.** The variance of the CV estimator eq. (B.7) is such that $\mathbb{V}[\hat{\mu}_{\mathcal{C}}^{\text{CV}}] = \beta^{\text{CV}}\mathbb{V}[\hat{\mu}_{p}]$, with $\beta^{\text{CV}} = (1+w)(1-\rho^2)$.

*Proof.* From eq. (2.5) and appendix A.1, we know that $\mathbb{V}[\hat{\mu}_{\mathcal{C}}^{\text{CV}}] = (1-\rho^2)\mathbb{V}[\hat{\mu}_{n}]$. The result follows from the fact that $\mathbb{V}[\hat{\mu}_{n}] = \mathbb{V}[f(\mathbf{X})]/n$ and $\mathbb{V}[\hat{\mu}_{p}] = \mathbb{V}[f(\mathbf{X})]/p$, with $n = p/(1+w)$. $\square$

**Proposition B.8.** $\beta^{\text{CV}} < 1$ if and only if $\rho^2 > w/(1+w)$.

$$ \text{Proof.} \quad \beta^{\text{CV}} < 1 \iff (1+w)(\rho^2 - 1) > -1 \iff \rho^2 > 1 - 1/(1+w) = w/(1+w). \quad \square $$

Combined with propositions B.3 and B.8, proposition B.5 indicates that the minimum value of $\rho^2$ that guarantees variance reduction is always (i.e., for any $w \in (0, 1)$) at least twice as large for the MFMC estimator than for the CV estimator. For instance, for $w = 1/3$, the MFMC estimator requires $\rho^2 > 3/4$ to have a reduced variance, while the CV estimator only requires $\rho^2 > 1/4$.

**Proposition B.9.** For $w, \rho^2 \in (0, 1)$, $\beta^{\text{CV}} < 1$ if and only if $(1+w)\beta^{\text{CV}} < \beta^{\text{MFMC}}$.

*Proof.* $\beta^{\text{CV}} < 1 \iff \rho^2 > w/(1+w) \iff 1-\rho^2 < 1/(1+w) \iff 1/(1-\rho^2) > 1+w$. Because $1/(1-\rho^2) - 1 = \rho^2/(1-\rho^2)$, $\beta^{\text{CV}} < 1 \iff \rho^2/(1-\rho^2) > w \iff w\rho^2/(1-\rho^2) > w^2$. We note that $w\rho^2/(1-\rho^2) = (\eta w)^2$, so that $\beta^{\text{CV}} < 1 \iff \eta w > w$. Therefore, because $\beta^{\text{MFMC}}/\beta^{\text{CV}} = (1+\eta w)^2/(1+w)$, $\beta^{\text{CV}} < 1 \iff \beta^{\text{MFMC}}/\beta^{\text{CV}} > (1+w)^2/(1+w) = 1+w$, which concludes the proof. $\square$

**Corollary B.10.** For $w, \rho^2 \in (0, 1)$, if $\beta^{\text{CV}} < 1$, then $\beta^{\text{CV}} < \beta^{\text{MFMC}}$.

*Proof.* Because $w > 0$, $\beta^{\text{MFMC}}/\beta^{\text{CV}} > 1+w \implies \beta^{\text{MFMC}}/\beta^{\text{CV}} > w$. Hence the result follows from proposition B.9. $\square$

Corollary B.10 implies that, in the $(w, \rho^2)$ variance reduction region of the CV estimator, which is larger than that of the MFMC estimator by propositions B.3, B.5, and B.8, the variance reduction of the CV estimator is always greater than that of the MFMC estimator (i.e., $\beta^{\text{CV}} < \beta^{\text{MFMC}}$).

Figure 5 depicts the $(w, \rho^2)$ variance reduction regions for the CV (fig. 5a) and MFMC (fig. 5b) estimators. We clearly see that for any $w \in (0, 1)$, a higher value of $\rho^2$ is required for the MFMC estimator to achieve variance reduction than for the CV estimator (red solid line). Figure 5 also illustrates an immediate consequence of propositions B.1 and B.7, namely that $\lim_{w \to 0} \beta^{\text{CV}} = \lim_{w \to 0} \beta^{\text{MFMC}} = 1-\rho^2$. We also observe that $\beta^{\text{CV}}$ converges to $1-\rho^2$ faster than $\beta^{\text{MFMC}}$, consistent with corollary B.10.

These properties illustrate the advantages of using surrogate-based control variates. Indeed, when the low-fidelity model $g$ corresponds to a surrogate model, not only $w \ll 1$, but the exact expectation $\mathbb{E}[g(\mathbf{X})]$ may be known, so that the exact CV approach can be used, resulting in a variance reduction factor close to $\rho^2$.
<figure><img src="image_7.png" /><figcaption>Figure 5: Contour lines of $\beta^{CV}$ (fig. 5a) and $\beta^{MFMC}$ (fig. 5b) as functions of $w$ and $\rho^2$. In each plot, the green region corresponds to pairs $(w, \rho^2)$ for which the variance is reduced, while the red region corresponds to pairs for which the variance is increased. The boundary between these two regions is represented by a red line.</figcaption></figure>

## C Optimal MLMC-MLCV variance estimator

The MLMC-MLCV variance estimator is given by eq. (3.20), with

$$ \hat{T}_{0}^{\text{MLCV}}(\alpha_{0}) = \hat{V}^{(0)}[Y_{0}] - \alpha_{0}^{\text{T}}(\hat{V}^{(0)}[\mathbf{Z}] - \sigma_{\mathbf{Z}}^{2}), \qquad (C.1) $$

$$ \hat{T}_{\ell}^{\text{MLCV}}(\alpha_{\ell}) = \hat{V}^{(\ell)}[Y_{\ell}] - \hat{V}^{(\ell)}[Y_{\ell-1}] - \alpha_{\ell}^{\text{T}}[\hat{V}^{(\ell)}[\mathbf{Z}_{1:\ell}] - \hat{V}^{(\ell)}[\tilde{\mathbf{Z}}] - (\sigma_{\tilde{\mathbf{Z}}_{1:\ell}}^{2} - \sigma_{\tilde{\mathbf{Z}}}^{2})], \quad \text{for } \ell > 1, \qquad (C.2) $$

where

$$ \mathbf{Z} = (Z_0, \dots, Z_L) = (g_0(\mathbf{X}), \dots, g_L(\mathbf{X})), \qquad \sigma_{\mathbf{Z}}^2 = \mathbb{V}[\mathbf{Z}], \qquad (C.3) $$

$$ \mathbf{Z}_{1:} = (Z_1, \dots, Z_L) = (g_1(\mathbf{X}), \dots, g_L(\mathbf{X})), \qquad \sigma_{\tilde{\mathbf{Z}}_{1:}}^2 = \mathbb{V}[\tilde{\mathbf{Z}}_{1:}], \qquad (C.4) $$

$$ \tilde{\mathbf{Z}} = (\tilde{Z}_0, \dots, \tilde{Z}_{L-1}) = (\tilde{g}_0(\mathbf{X}), \dots, \tilde{g}_{L-1}(\mathbf{X})), \qquad \sigma_{\tilde{\mathbf{Z}}}^2 = \mathbb{V}[\tilde{\mathbf{Z}}], \qquad (C.5) $$

$$ \tilde{g}_{\ell-1} = g_{\ell} - h_{\ell}, \quad \text{for } \ell > 0. \qquad (C.6) $$

Furthermore, we assume that the expected values of the control variates are known, so that we use the following variance estimators:

$$ \hat{V}^{(0)}[\mathbf{Z}] = \hat{E}^{(0)}[\tilde{\mathbf{Z}}^2], \quad \hat{V}^{(\ell)}[\mathbf{Z}_{1:\ell}] - \hat{V}^{(\ell)}[\tilde{\mathbf{Z}}] = \hat{E}^{(\ell)}[\tilde{\mathbf{Z}}_{1:\ell}^{\odot 2}] - \hat{E}^{(\ell)}[\tilde{\tilde{\mathbf{Z}}}^{\odot 2}] = \hat{E}^{(\ell)}[\tilde{\mathbf{Z}}_{1:\ell}^{\odot 2} - \tilde{\tilde{\mathbf{Z}}}^{\odot 2}]. \qquad (C.7) $$

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

$$ \alpha_{0}^{*} = \Sigma_{0}^{-1} c_{0}, \qquad \Sigma_{0} = \mathbb{C}[\tilde{\mathbf{Z}}^{\odot 2}], \qquad c_{0} = \mathbb{C}[\tilde{\mathbf{Z}}^{\odot 2}, \bar{Y}_{0}^{2}], \qquad (C.8) $$

$$ \alpha_{\ell}^{*} = \Sigma_{\ell}^{-1} c_{\ell}, \qquad \Sigma_{\ell} = \mathbb{C}[\tilde{\mathbf{Z}}_{1:\ell}^{\odot 2} - \tilde{\tilde{\mathbf{Z}}}^{\odot 2}], \qquad c_{\ell} = \mathbb{C}[\tilde{\mathbf{Z}}_{1:\ell}^{\odot 2} - \tilde{\tilde{\mathbf{Z}}}^{\odot 2}, \bar{Y}_{\ell}^{2} - \bar{Y}_{\ell-1}^{2}], \quad \text{for } \ell > 0. \qquad (C.9) $$
<figure><img src="image_7.png" /><figcaption>(a) Standard (exact) CV.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>(b) MFMC.</figcaption></figure>

<caption>Figure 5: Contour lines of $\beta^{CV}$ (fig. 5a) and $\beta^{MFMC}$ (fig. 5b) as functions of $w$ and $\rho^2$. In each plot, the green region corresponds to pairs $(w, \rho^2)$ for which the variance is reduced, while the red region corresponds to pairs for which the variance is increased. The boundary between these two regions is represented by a red line.</caption>

## C Optimal MLMC-MLCV variance estimator

The MLMC-MLCV variance estimator is given by eq. (3.20), with

$$ \hat{T}_{0}^{\text{MLCV}}(\alpha_{0}) = \hat{V}^{(0)}[Y_{0}] - \alpha_{0}^{\text{T}}(\hat{V}^{(0)}[\mathbf{Z}] - \sigma_{\mathbf{Z}}^{2}), \quad (C.1) $$

$$ \hat{T}_{\ell}^{\text{MLCV}}(\alpha_{\ell}) = \hat{V}^{(\ell)}[Y_{\ell}] - \hat{V}^{(\ell)}[Y_{\ell-1}] - \alpha_{\ell}^{\text{T}}[\hat{V}^{(\ell)}[\mathbf{Z}_{1\cdot}] - \hat{V}^{(\ell)}[\tilde{\mathbf{Z}}] - (\sigma_{\mathbf{Z}_{1\cdot}}^{2} - \sigma_{\tilde{\mathbf{Z}}}^{2})], \quad \text{for } \ell > 1, \quad (C.2) $$

where

$$ \mathbf{Z} = (Z_0, \dots, Z_L) = (g_0(\mathbf{X}), \dots, g_L(\mathbf{X})), \qquad \sigma_{\mathbf{Z}}^2 = \mathbb{V}[\mathbf{Z}], \quad (C.3) $$

$$ \mathbf{Z}_{1\cdot} = (Z_1, \dots, Z_L) = (g_1(\mathbf{X}), \dots, g_L(\mathbf{X})), \qquad \sigma_{\mathbf{Z}_{1\cdot}}^2 = \mathbb{V}[\mathbf{Z}_{1\cdot}], \quad (C.4) $$

$$ \tilde{\mathbf{Z}} = (\tilde{Z}_0, \dots, \tilde{Z}_{L-1}) = (\tilde{g}_0(\mathbf{X}), \dots, \tilde{g}_{L-1}(\mathbf{X})), \qquad \sigma_{\tilde{\mathbf{Z}}}^2 = \mathbb{V}[\tilde{\mathbf{Z}}], \quad (C.5) $$

$$ \tilde{g}_{\ell-1} = g_{\ell} - h_{\ell}, \quad \text{for } \ell > 0. \quad (C.6) $$

Furthermore, we assume that the expected values of the control variates are known, so that we use the following variance estimators:

$$ \hat{V}^{(0)}[\mathbf{Z}] = \hat{E}^{(0)}[\tilde{\mathbf{Z}}^2], \quad \hat{V}^{(\ell)}[\mathbf{Z}_{1\cdot}] - \hat{V}^{(\ell)}[\tilde{\mathbf{Z}}] = \hat{E}^{(\ell)}[\tilde{\mathbf{Z}}^{\odot 2}] - \hat{E}^{(\ell)}[\tilde{\tilde{\mathbf{Z}}}^{\odot 2}] = \hat{E}^{(\ell)}[\tilde{\mathbf{Z}}_{1\cdot}^{\odot 2} - \tilde{\tilde{\mathbf{Z}}}^{\odot 2}]. \quad (C.7) $$

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

$$ \alpha_{0}^{*} = \Sigma_{0}^{-1} c_{0}, \qquad \Sigma_{0} = \mathbb{C}[\tilde{\mathbf{Z}}^{\odot 2}], \qquad c_{0} = \mathbb{C}[\tilde{\mathbf{Z}}^{\odot 2}, \bar{Y}_{0}^{2}], \quad (C.8) $$

$$ \alpha_{\ell}^{*} = \Sigma_{\ell}^{-1} c_{\ell}, \qquad \Sigma_{\ell} = \mathbb{C}[\tilde{\mathbf{Z}}_{1\cdot}^{\odot 2} - \tilde{\tilde{\mathbf{Z}}}^{\odot 2}], \qquad c_{\ell} = \mathbb{C}[\tilde{\mathbf{Z}}_{1\cdot}^{\odot 2} - \tilde{\tilde{\mathbf{Z}}}^{\odot 2}, \bar{Y}_{\ell}^{2} - \bar{Y}_{\ell-1}^{2}], \quad \text{for } \ell > 0. \quad (C.9) $$
Note that $\Sigma_{\ell}$ is the same for all $\ell > 0$. For PC-based control variates, $g_{\ell}(\mathbf{X}) = \sum_{k=0}^{P_g^{\ell}} \mathfrak{g}_{\ell,k} \Psi_k(\mathbf{X})$ and $h_{\ell}(\mathbf{X}) = \sum_{k=0}^{P_h^{\ell}} \mathfrak{h}_{\ell,k} \Psi_k(\mathbf{X})$, letting $P_{\tilde{g}}^{\ell} := \max(P_g^{\ell}, P_h^{\ell})$, we have

$$ \tilde{Z}_{\ell-1} = \sum_{k=0}^{P_{\tilde{g}}^{\ell}} \tilde{\mathfrak{g}}_{\ell-1,k} \Psi_k(\mathbf{X}), \quad \text{with } \tilde{\mathfrak{g}}_{\ell-1,k} := \mathfrak{g}_{\ell,k} - \mathfrak{h}_{\ell,k}, \quad \text{for } \ell = 1, \dots, L, \qquad (C.10) $$

where $\mathfrak{g}_{\ell,k} := 0$ for $k > P_g^{\ell}$ and $\mathfrak{h}_{\ell,k} := 0$ for $k > P_h^{\ell}$. Consequently,

$$ \mathbb{E}[Z_{\ell}] = \mathfrak{g}_{\ell,0}, \qquad \mathbb{E}[\tilde{Z}_{\ell-1}] = \tilde{\mathfrak{g}}_{\ell-1,0}, \qquad (C.11) $$

$$ \mathbb{V}[Z_{\ell}] = \sum_{k=1}^{P_g^{\ell}} \mathfrak{g}_{\ell,k}^2, \qquad \mathbb{V}[\tilde{Z}_{\ell-1}] = \sum_{k=1}^{P_g^{\ell}} \tilde{\mathfrak{g}}_{\ell-1,k}^2. \qquad (C.12) $$

Furthermore,

$$ [\Sigma_0]_{m,m'} = \sum_{i,j=1}^{P_g^m} \sum_{q,r=1}^{P_g^{m'}} \mathfrak{g}_{m,i} \mathfrak{g}_{m,j} \mathfrak{g}_{m',q} \mathfrak{g}_{m',r} (\Phi_{ijqr} - \delta_{ij}\delta_{qr}), \quad \text{for } m, m' = 0, \dots, L, \qquad (C.13) $$

where $\Phi_{ijqr} := \mathbb{E}[\Psi_i(\mathbf{X})\Psi_j(\mathbf{X})\Psi_q(\mathbf{X})\Psi_r(\mathbf{X})]$ are the entries of the fourth-order Galerkin product tensor, which is a well-known object in stochastic Galerkin methods [20, 32], and $\delta_{ij}$ denotes the Kronecker delta. Besides, noticing that $\bar{\mathfrak{Z}}_{1:}^{\circ 2} - \tilde{\mathfrak{Z}}^{\circ 2} = (\bar{\mathfrak{Z}}_1: - \tilde{\mathfrak{Z}}) \odot (\bar{\mathfrak{Z}}_1: + \tilde{\mathfrak{Z}}) = \bar{\mathfrak{W}} \odot (\bar{\mathfrak{Z}}_1: + \tilde{\mathfrak{Z}})$, with $\bar{\mathfrak{W}} = \mathbf{W} - \mu_{\mathbf{W}}$ and the definitions in eq. (3.30), we have

$$ [\Sigma_{\ell}]_{m,m'} = A_{m,m'} + B_{m,m'} + C_{m,m'} + C_{m',m}, \quad \text{for } \ell, m, m' = 1, \dots, L, \qquad (C.14) $$

where

$$ A_{m,m'} = \mathbb{C}[\bar{W}_m \bar{Z}_m, \bar{W}_{m'} \bar{Z}_{m'}] \qquad (C.15) $$

$$ = \sum_{i=1}^{P_h^m} \sum_{j=1}^{P_g^m} \sum_{q=1}^{P_h^{m'}} \sum_{r=1}^{P_g^{m'}} \mathfrak{h}_{m,i} \mathfrak{g}_{m,j} \mathfrak{h}_{m',q} \mathfrak{g}_{m',r} (\Phi_{ijqr} - \delta_{ij}\delta_{qr}), \qquad (C.16) $$

$$ B_{m,m'} = \mathbb{C}[\bar{W}_m \bar{\tilde{Z}}_{m-1}, \bar{W}_{m'} \bar{\tilde{Z}}_{m'-1}] \qquad (C.17) $$

$$ = \sum_{i=1}^{P_h^m} \sum_{j=1}^{P_g^m} \sum_{q=1}^{P_h^{m'}} \sum_{r=1}^{P_g^{m'}} \mathfrak{h}_{m,i} \tilde{\mathfrak{g}}_{m-1,j} \mathfrak{h}_{m',q} \tilde{\mathfrak{g}}_{m'-1,r} (\Phi_{ijqr} - \delta_{ij}\delta_{qr}), \qquad (C.18) $$

$$ C_{m,m'} = \mathbb{C}[\bar{W}_m \bar{\tilde{Z}}_{m-1}, \bar{W}_{m'} \bar{Z}_{m'}] \qquad (C.19) $$

$$ = \sum_{i=1}^{P_h^m} \sum_{j=1}^{P_g^m} \sum_{q=1}^{P_h^{m'}} \sum_{r=1}^{P_g^{m'}} \mathfrak{h}_{m,i} \tilde{\mathfrak{g}}_{m-1,j} \mathfrak{h}_{m',q} \tilde{\mathfrak{g}}_{m',r} (\Phi_{ijqr} - \delta_{ij}\delta_{qr}). \qquad (C.20) $$

## D Taylor surrogate for the numerical test case

In our example, $f_{\ell}$ is only differentiable in $[-\pi, \pi]^3 \times [\nu_{\min}, \nu_{\max}] \times ((-1,0) \cup (0,1))^3$. Consequently, the Taylor polynomial surrogate cannot be used directly as defined in (2.17) and (2.18), because
Note that $\Sigma_\ell$ is the same for all $\ell > 0$. For PC-based control variates, $g_\ell(\mathbf{X}) = \sum_{k=0}^{P_g^\ell} \mathfrak{g}_{\ell,k} \Psi_k(\mathbf{X})$ and $h_\ell(\mathbf{X}) = \sum_{k=0}^{P_h^\ell} \mathfrak{h}_{\ell,k} \Psi_k(\mathbf{X})$, letting $P_{\tilde{g}}^\ell := \max(P_g^\ell, P_h^\ell)$, we have

$$ \tilde{Z}_{\ell-1} = \sum_{k=0}^{P_{\tilde{g}}^{\ell}} \tilde{\mathfrak{g}}_{\ell-1,k} \Psi_k(\mathbf{X}), \quad \text{with } \tilde{\mathfrak{g}}_{\ell-1,k} := \mathfrak{g}_{\ell,k} - \mathfrak{h}_{\ell,k}, \quad \text{for } \ell = 1, \dots, L, \qquad (C.10) $$

where $\mathfrak{g}_{\ell,k} := 0$ for $k > P_g^\ell$ and $\mathfrak{h}_{\ell,k} := 0$ for $k > P_h^\ell$. Consequently,

$$ \mathbb{E}[Z_{\ell}] = \mathfrak{g}_{\ell,0}, \qquad \mathbb{E}[\tilde{Z}_{\ell-1}] = \tilde{\mathfrak{g}}_{\ell-1,0}, \qquad (C.11) $$

$$ \mathbb{V}[Z_{\ell}] = \sum_{k=1}^{P_g^\ell} \mathfrak{g}_{\ell,k}^2, \qquad \mathbb{V}[\tilde{Z}_{\ell-1}] = \sum_{k=1}^{P_g^\ell} \tilde{\mathfrak{g}}_{\ell-1,k}^2. \qquad (C.12) $$

Furthermore,

$$ [\Sigma_0]_{m,m'} = \sum_{i,j=1}^{P_g^m} \sum_{q,r=1}^{P_g^{m'}} \mathfrak{g}_{m,i} \mathfrak{g}_{m,j} \mathfrak{g}_{m',q} \mathfrak{g}_{m',r} (\Phi_{ijqr} - \delta_{ij}\delta_{qr}), \quad \text{for } m, m' = 0, \dots, L, \qquad (C.13) $$

where $\Phi_{ijqr} := \mathbb{E}[\Psi_i(\mathbf{X})\Psi_j(\mathbf{X})\Psi_q(\mathbf{X})\Psi_r(\mathbf{X})]$ are the entries of the fourth-order Galerkin product tensor, which is a well-known object in stochastic Galerkin methods [20, 32], and $\delta_{ij}$ denotes the Kronecker delta. Besides, noticing that $\bar{\mathfrak{Z}}_{1:}^{\circ 2} - \tilde{\mathfrak{Z}}^{\circ 2} = (\bar{\mathfrak{Z}}_1: - \tilde{\mathfrak{Z}}) \odot (\bar{\mathfrak{Z}}_1: + \tilde{\mathfrak{Z}}) = \bar{\mathfrak{W}} \odot (\bar{\mathfrak{Z}}_1: + \tilde{\mathfrak{Z}})$, with $\bar{\mathfrak{W}} = \mathbf{W} - \mu_{\mathbf{W}}$ and the definitions in eq. (3.30), we have

$$ [\Sigma_{\ell}]_{m,m'} = A_{m,m'} + B_{m,m'} + C_{m,m'} + C_{m',m}, \quad \text{for } \ell, m, m' = 1, \dots, L, \qquad (C.14) $$

where

$$ A_{m,m'} = \mathbb{C}[\bar{W}_m \bar{Z}_m, \bar{W}_{m'} \bar{Z}_{m'}] \qquad (C.15) $$

$$ = \sum_{i=1}^{P_h^m} \sum_{j=1}^{P_g^m} \sum_{q=1}^{P_h^{m'}} \sum_{r=1}^{P_g^{m'}} \mathfrak{h}_{m,i} \mathfrak{g}_{m,j} \mathfrak{h}_{m',q} \mathfrak{g}_{m',r} (\Phi_{ijqr} - \delta_{ij}\delta_{qr}), \qquad (C.16) $$

$$ B_{m,m'} = \mathbb{C}[\bar{W}_m \bar{\tilde{Z}}_{m-1}, \bar{W}_{m'} \bar{\tilde{Z}}_{m'-1}] \qquad (C.17) $$

$$ = \sum_{i=1}^{P_h^m} \sum_{j=1}^{P_g^m} \sum_{q=1}^{P_h^{m'}} \sum_{r=1}^{P_g^{m'}} \mathfrak{h}_{m,i} \tilde{\mathfrak{g}}_{m-1,j} \mathfrak{h}_{m',q} \tilde{\mathfrak{g}}_{m'-1,r} (\Phi_{ijqr} - \delta_{ij}\delta_{qr}), \qquad (C.18) $$

$$ C_{m,m'} = \mathbb{C}[\bar{W}_m \bar{\tilde{Z}}_{m-1}, \bar{W}_{m'} \bar{Z}_{m'}] \qquad (C.19) $$

$$ = \sum_{i=1}^{P_h^m} \sum_{j=1}^{P_g^m} \sum_{q=1}^{P_h^{m'}} \sum_{r=1}^{P_g^{m'}} \mathfrak{h}_{m,i} \tilde{\mathfrak{g}}_{m-1,j} \mathfrak{h}_{m',q} \tilde{\mathfrak{g}}_{m',r} (\Phi_{ijqr} - \delta_{ij}\delta_{qr}). \qquad (C.20) $$

## D Taylor surrogate for the numerical test case

In our example, $f_\ell$ is only differentiable in $[-\pi, \pi]^3 \times [\nu_{\min}, \nu_{\max}] \times ((-1,0) \cup (0,1))^3$. Consequently, the Taylor polynomial surrogate cannot be used directly as defined in (2.17) and (2.18), because
the Jacobian and Hessian matrices are not defined at $\boldsymbol{\mu_X}$. Instead, we define the first-order Taylor surrogate as

$$f_{\ell}(\mathbf{X}) \approx g_{\ell}^{\mathrm{T}_1}(\mathbf{X}) = f_{\ell}(\boldsymbol{\mu_X}) + \sum_{i=1}^{7} g_{\ell,i}^{\mathrm{T}_1}(\boldsymbol{\mu_X}; \mathbf{X}), \qquad (\mathrm{D.1})$$

where, for $i = 1, \dots, 4$, $g_{\ell,i}^{\mathrm{T}_1}(\boldsymbol{\mu_X}; \mathbf{X}) = (X_i - \mu_{X_i}) \frac{\partial f_{\ell}}{\partial X_i}(\boldsymbol{\mu_X})$, and, for $i = 5, \dots, 7$,

$$g_{\ell,i}^{\mathrm{T}_1}(\boldsymbol{\mu_X}; \mathbf{X}) = (X_i - \mu_{X_i}) \times \begin{cases} \frac{\partial f_{\ell}}{\partial X_i}(\boldsymbol{\mu_X}), & \mu_{X_i} \neq 0, \\ \lim_{\boldsymbol{\mu}'_X \to \boldsymbol{\mu}_X^{i,0^{-}}} \frac{\partial f_{\ell}}{\partial X_i}(\boldsymbol{\mu}'_X), & \mu_{X_i} = 0, X_i < 0, \\ \lim_{\boldsymbol{\mu}'_X \to \boldsymbol{\mu}_X^{i,0^{+}}} \frac{\partial f_{\ell}}{\partial X_i}(\boldsymbol{\mu}'_X), & \mu_{X_i} = 0, X_i > 0, \\ 0 & \mu_{X_i} = X_i = 0, \end{cases} \qquad (\mathrm{D.2})$$

where $\boldsymbol{\mu}_X^{i,0^{\pm}} = (\mu_{X_1}, \dots, \mu_{X_{i-1}}, 0^{\pm}, \mu_{X_{i+1}}, \dots, \mu_{X_7})$, which is now well-defined.

With our choice of distributions for $\mathbf{X}$ given in (4.11), we have $\boldsymbol{\mu_X} = (0, 0, 0, 0.005, 0, 0, 0)$, and the first-order Taylor surrogate is defined by (D.1), with

$$g_{\ell,1}^{\mathrm{T}_1}(\boldsymbol{\mu_X}; \mathbf{X}) = 7X_1 \sum_{k=1}^{K} B_k^{\ell}(T; \boldsymbol{\mu_X}) \sum_{j=1}^{N_{\ell}} w_j \sin(k\pi x_j) \mathcal{F}_2(x_j), \qquad (\mathrm{D.3})$$

$$g_{\ell,2}^{\mathrm{T}_1}(\boldsymbol{\mu_X}; \mathbf{X}) = 0, \qquad (\mathrm{D.4})$$

$$g_{\ell,3}^{\mathrm{T}_1}(\boldsymbol{\mu_X}; \mathbf{X}) = 0, \qquad (\mathrm{D.5})$$

$$g_{\ell,4}^{\mathrm{T}_1}(\boldsymbol{\mu_X}; \mathbf{X}) = -(X_4 - 0.005)\pi^2 T \sum_{k=1}^{K} k^2 \exp(-0.005k^2\pi^2 T) A_k^{\ell}(\boldsymbol{\mu_X}) \sum_{j=1}^{N_{\ell}} w_j \sin(k\pi x_j), \qquad (\mathrm{D.6})$$

$$g_{\ell,5}^{\mathrm{T}_1}(\boldsymbol{\mu_X}; \mathbf{X}) = 400 |X_5| \sum_{k=1}^{K} B_k^{\ell}(T; \boldsymbol{\mu_X}) \sum_{j=1}^{N_{\ell}} w_j \sin(k\pi x_j) \mathcal{F}_1(x_j), \qquad (\mathrm{D.7})$$

$$g_{\ell,6}^{\mathrm{T}_1}(\boldsymbol{\mu_X}; \mathbf{X}) = 400 |X_6| \sum_{k=1}^{K} B_k^{\ell}(T; \boldsymbol{\mu_X}) \sum_{j=1}^{N_{\ell}} w_j \sin(k\pi x_j) \mathcal{F}_1(x_j), \qquad (\mathrm{D.8})$$

$$g_{\ell,7}^{\mathrm{T}_1}(\boldsymbol{\mu_X}; \mathbf{X}) = 400 |X_7| \sum_{k=1}^{K} B_k^{\ell}(T; \boldsymbol{\mu_X}) \sum_{j=1}^{N_{\ell}} w_j \sin(k\pi x_j) \mathcal{F}_1(x_j). \qquad (\mathrm{D.9})$$

Because of the piecewise definition of $g_{\ell}^{\mathrm{T}_1}$, the identity $\mathbb{E}[g_{\ell}^{\mathrm{T}_1}(\mathbf{X})] = f_{\ell}(\boldsymbol{\mu_X})$ for a regular first-order Taylor surrogate no longer holds. Instead, we have

$$\mathbb{E}[g_{\ell}^{\mathrm{T}_1}(\mathbf{X})] = f_{\ell}(\boldsymbol{\mu_X}) + 600 \sum_{k=1}^{K} B_k^{\ell}(T; \boldsymbol{\mu_X}) \sum_{j=1}^{N_{\ell}} w_j \sin(k\pi x_j) \mathcal{F}_1(x_j). \qquad (\mathrm{D.10})$$

## References

[1] T. W. Anderson. *An Introduction to Multivariate Statistical Analysis*. Wiley Series in Probability and Statistics. Wiley-Interscience, Hoboken, N.J, 3rd ed edition, 2003.
the Jacobian and Hessian matrices are not defined at $\boldsymbol{\mu_X}$. Instead, we define the first-order Taylor surrogate as

$$f_{\ell}(\mathbf{X}) \approx g_{\ell}^{\mathrm{T}_1}(\mathbf{X}) = f_{\ell}(\boldsymbol{\mu_X}) + \sum_{i=1}^{7} g_{\ell,i}^{\mathrm{T}_1}(\boldsymbol{\mu_X}; \mathbf{X}), \qquad (\mathrm{D.1})$$

where, for $i = 1, \dots, 4$, $g_{\ell,i}^{\mathrm{T}_1}(\boldsymbol{\mu_X}; \mathbf{X}) = (X_i - \mu_{X_i}) \frac{\partial f_{\ell}}{\partial X_i}(\boldsymbol{\mu_X})$, and, for $i = 5, \dots, 7$,

$$g_{\ell,i}^{\mathrm{T}_1}(\boldsymbol{\mu_X}; \mathbf{X}) = (X_i - \mu_{X_i}) \times \left\{ \begin{array}{ll} \frac{\partial f_{\ell}}{\partial X_i}(\boldsymbol{\mu_X}), & \mu_{X_i} \neq 0, \\ \lim_{\boldsymbol{\mu}'_X \to \boldsymbol{\mu}_X^{i,0^{-}}} \frac{\partial f_{\ell}}{\partial X_i}(\boldsymbol{\mu}'_X), & \mu_{X_i} = 0, X_i < 0, \\ \lim_{\boldsymbol{\mu}'_X \to \boldsymbol{\mu}_X^{i,0^{+}}} \frac{\partial f_{\ell}}{\partial X_i}(\boldsymbol{\mu}'_X), & \mu_{X_i} = 0, X_i > 0, \\ 0 & \mu_{X_i} = X_i = 0, \end{array} \right. \qquad (\mathrm{D.2})$$

where $\boldsymbol{\mu}_X^{i,0^{\pm}} = (\mu_{X_1}, \dots, \mu_{X_{i-1}}, 0^{\pm}, \mu_{X_{i+1}}, \dots, \mu_{X_7})$, which is now well-defined.

With our choice of distributions for $\mathbf{X}$ given in (4.11), we have $\boldsymbol{\mu_X} = (0, 0, 0, 0.005, 0, 0, 0)$, and the first-order Taylor surrogate is defined by (D.1), with

$$g_{\ell,1}^{\mathrm{T}_1}(\boldsymbol{\mu_X}; \mathbf{X}) = 7X_1 \sum_{k=1}^{K} B_k^{\ell}(T; \boldsymbol{\mu_X}) \sum_{j=1}^{N_{\ell}} w_j \sin(k\pi x_j) \mathcal{F}_2(x_j), \qquad (\mathrm{D.3})$$

$$g_{\ell,2}^{\mathrm{T}_1}(\boldsymbol{\mu_X}; \mathbf{X}) = 0, \qquad (\mathrm{D.4})$$

$$g_{\ell,3}^{\mathrm{T}_1}(\boldsymbol{\mu_X}; \mathbf{X}) = 0, \qquad (\mathrm{D.5})$$

$$g_{\ell,4}^{\mathrm{T}_1}(\boldsymbol{\mu_X}; \mathbf{X}) = -(X_4 - 0.005)\pi^2 T \sum_{k=1}^{K} k^2 \exp(-0.005k^2\pi^2 T) A_k^{\ell}(\boldsymbol{\mu_X}) \sum_{j=1}^{N_{\ell}} w_j \sin(k\pi x_j), \qquad (\mathrm{D.6})$$

$$g_{\ell,5}^{\mathrm{T}_1}(\boldsymbol{\mu_X}; \mathbf{X}) = 400 |X_5| \sum_{k=1}^{K} B_k^{\ell}(T; \boldsymbol{\mu_X}) \sum_{j=1}^{N_{\ell}} w_j \sin(k\pi x_j) \mathcal{F}_1(x_j), \qquad (\mathrm{D.7})$$

$$g_{\ell,6}^{\mathrm{T}_1}(\boldsymbol{\mu_X}; \mathbf{X}) = 400 |X_6| \sum_{k=1}^{K} B_k^{\ell}(T; \boldsymbol{\mu_X}) \sum_{j=1}^{N_{\ell}} w_j \sin(k\pi x_j) \mathcal{F}_1(x_j), \qquad (\mathrm{D.8})$$

$$g_{\ell,7}^{\mathrm{T}_1}(\boldsymbol{\mu_X}; \mathbf{X}) = 400 |X_7| \sum_{k=1}^{K} B_k^{\ell}(T; \boldsymbol{\mu_X}) \sum_{j=1}^{N_{\ell}} w_j \sin(k\pi x_j) \mathcal{F}_1(x_j). \qquad (\mathrm{D.9})$$

Because of the piecewise definition of $g_{\ell}^{\mathrm{T}_1}$, the identity $\mathbb{E}[g_{\ell}^{\mathrm{T}_1}(\mathbf{X})] = f_{\ell}(\boldsymbol{\mu_X})$ for a regular first-order Taylor surrogate no longer holds. Instead, we have

$$\mathbb{E}[g_{\ell}^{\mathrm{T}_1}(\mathbf{X})] = f_{\ell}(\boldsymbol{\mu_X}) + 600 \sum_{k=1}^{K} B_k^{\ell}(T; \boldsymbol{\mu_X}) \sum_{j=1}^{N_{\ell}} w_j \sin(k\pi x_j) \mathcal{F}_1(x_j). \qquad (\mathrm{D.10})$$

## References

[1] T. W. Anderson. *An Introduction to Multivariate Statistical Analysis*. Wiley Series in Probability and Statistics. Wiley-Interscience, Hoboken, N.J, 3rd ed edition, 2003.
[2] Ivo Babuška, Fabio Nobile, and Raúl Tempone. A Stochastic Collocation Method for Elliptic Partial Differential Equations with Random Input Data. *SIAM Journal on Numerical Analysis*, 45(3):1005-1034, 2007.

[3] Marc Berveiller, Bruno Sudret, and Maurice Lemaire. Stochastic finite element: A non intrusive approach by regression. *European Journal of Computational Mechanics*, 15(1-3):81-92, 2006.

[4] Claudio Bierig and Alexey Chernov. Convergence analysis of multilevel Monte Carlo variance estimators and application for random obstacle problems. *Numerische Mathematik*, 130(4):579-613, 2015.

[5] Claudio Bierig and Alexey Chernov. Estimation of arbitrary order central statistical moments by the multilevel Monte Carlo method. *Stochastics and Partial Differential Equations Analysis and Computations*, 4(1):3-40, 2016.

[6] Géraud Blatman and Bruno Sudret. Efficient computation of global sensitivity indices using sparse polynomial chaos expansions. *Reliability Engineering & System Safety*, 95(11):1216-1229, 2010.

[7] Géraud Blatman and Bruno Sudret. Adaptive sparse polynomial chaos expansion based on least angle regression. *Journal of Computational Physics*, 230(6):2345-2367, 2011.

[8] François-Xavier Briol, Chris J. Oates, Mark Girolami, Michael A. Osborne, and Dino Sejdi-novic. Probabilistic Integration: A Role in Statistical Computation? *Statistical Science*, 34(1):1-22, 2019.

[9] Robert H. Cameron and William T. Martin. The orthogonal development of non-linear functionals in series of Fourier-Hermite functionals. *Annals of Mathematics*, 48(2):385-392, 1947.

[10] K. A. Cliffe, M. B. Giles, R. Scheichl, and A. L. Teckentrup. Multilevel Monte Carlo methods and applications to elliptic PDEs with random coefficients. *Computing and Visualization in Science*, 14(1):3-15, 2011.

[11] Patrick R. Conrad and Youssef M. Marzouk. Adaptive Smolyak Pseudospectral Approximations. *SIAM Journal on Scientific Computing*, 35(6):A2643-A2670, 2013.

[12] Paul G. Constantine, Michael S. Eldred, and Eric T. Phipps. Sparse pseudospectral approximation method. *Computer Methods in Applied Mechanics and Engineering*, 229-232:1-12, 2012.

[13] Josef Dick, Frances Y. Kuo, and Ian H. Sloan. High-dimensional integration: The quasi-monte carlo way. *Acta Numerica*, 22:133-288, 2013.

[14] Bradley Efron, Trevor Hastie, Iain Johnstone, and Robert Tibshirani. Least Angle Regression. *The Annals of Statistics*, 32(2):407-451, 2004.

[15] Reda EL Amri, Rodolphe Le Riche, Céline Helbert, Christophette Blanchet-Scalliet, and Sébastien Da Veiga. A sampling criterion for constrained bayesian optimization with uncertainties. *arXiv preprint arXiv:2103.05706*, 2021.

[16] Jamie Fox and Giray Ökten. Polynomial chaos as a control variate method. *SIAM Journal on Scientific Computing*, 43(3):A2268-A2294, 2021.
[2] Ivo Babuška, Fabio Nobile, and Raúl Tempone. A Stochastic Collocation Method for Elliptic Partial Differential Equations with Random Input Data. *SIAM Journal on Numerical Analysis*, 45(3):1005-1034, 2007.

[3] Marc Berveiller, Bruno Sudret, and Maurice Lemaire. Stochastic finite element: A non intrusive approach by regression. *European Journal of Computational Mechanics*, 15(1-3):81-92, 2006.

[4] Claudio Bierig and Alexey Chernov. Convergence analysis of multilevel Monte Carlo variance estimators and application for random obstacle problems. *Numerische Mathematik*, 130(4):579-613, 2015.

[5] Claudio Bierig and Alexey Chernov. Estimation of arbitrary order central statistical moments by the multilevel Monte Carlo method. *Stochastics and Partial Differential Equations Analysis and Computations*, 4(1):3-40, 2016.

[6] Géraud Blatman and Bruno Sudret. Efficient computation of global sensitivity indices using sparse polynomial chaos expansions. *Reliability Engineering & System Safety*, 95(11):1216-1229, 2010.

[7] Géraud Blatman and Bruno Sudret. Adaptive sparse polynomial chaos expansion based on least angle regression. *Journal of Computational Physics*, 230(6):2345-2367, 2011.

[8] François-Xavier Briol, Chris J. Oates, Mark Girolami, Michael A. Osborne, and Dino Sejdi-novic. Probabilistic Integration: A Role in Statistical Computation? *Statistical Science*, 34(1):1-22, 2019.

[9] Robert H. Cameron and William T. Martin. The orthogonal development of non-linear functionals in series of Fourier-Hermite functionals. *Annals of Mathematics*, 48(2):385-392, 1947.

[10] K. A. Cliffe, M. B. Giles, R. Scheichl, and A. L. Teckentrup. Multilevel Monte Carlo methods and applications to elliptic PDEs with random coefficients. *Computing and Visualization in Science*, 14(1):3-15, 2011.

[11] Patrick R. Conrad and Youssef M. Marzouk. Adaptive Smolyak Pseudospectral Approximations. *SIAM Journal on Scientific Computing*, 35(6):A2643-A2670, 2013.

[12] Paul G. Constantine, Michael S. Eldred, and Eric T. Phipps. Sparse pseudospectral approximation method. *Computer Methods in Applied Mechanics and Engineering*, 229-232:1-12, 2012.

[13] Josef Dick, Frances Y. Kuo, and Ian H. Sloan. High-dimensional integration: The quasi-monte carlo way. *Acta Numerica*, 22:133-288, 2013.

[14] Bradley Efron, Trevor Hastie, Iain Johnstone, and Robert Tibshirani. Least Angle Regression. *The Annals of Statistics*, 32(2):407-451, 2004.

[15] Reda EL Amri, Rodolphe Le Riche, Céline Helbert, Christophette Blanchet-Scalliet, and Sébastien Da Veiga. A sampling criterion for constrained bayesian optimization with uncertainties. *arXiv preprint arXiv:2103.05706*, 2021.

[16] Jamie Fox and Giray Ökten. Polynomial chaos as a control variate method. *SIAM Journal on Scientific Computing*, 43(3):A2268-A2294, 2021.
[17] Sarthak Garg, Neeraj Sood, and Costas D. Sarris. Uncertainty quantification of ray-tracing based wireless propagation models with a control variate-polynomial chaos expansion method. In *2015 IEEE International Symposium on Antennas and Propagation & USNC/URSI National Radio Science Meeting*, pages 1776–1777, 2015.

[18] G Geraci, M Eldred, and G Iaccarino. A multifidelity control variate approach for the multilevel Monte Carlo technique. Technical report, Stanford University, 2015.

[19] Gianluca Geraci, Michael S. Eldred, and Gianluca Iaccarino. A multifidelity multilevel Monte Carlo method for uncertainty propagation in aerospace applications. In *19th AIAA Non-Deterministic Approaches Conference, Grapevine, Texas, 2017*. American Institute of Aeronautics and Astronautics.

[20] Roger G. Ghanem and Pol D. Spanos. *Stochastic Finite Elements: A Spectral Approach*. Courier Corporation, 2003.

[21] Michael B. Giles. Multilevel Monte Carlo Path Simulation. *Operations Research*, 56(3):607–617, 2008.

[22] Michael B. Giles. Multilevel Monte Carlo methods. *Acta Numerica*, 24:259–328, 2015.

[23] Michael B. Giles and Benjamin J. Waterhouse. Multilevel quasi-Monte Carlo path simulation. *Advanced Financial Modelling, Radon Series on Computational and Applied Mathematics*, 8:165–181, 2009.

[24] Alex A. Gorodetsky, Gianluca Geraci, Michael S. Eldred, and John D. Jakeman. A generalized approximate control variate framework for multifidelity uncertainty quantification. *Journal of Computational Physics*, 408:109257, 2020.

[25] Zixi Gu and Costas D. Sarris. Multi-parametric uncertainty quantification with a hybrid monte-carlo / polynomial chaos expansion fdtd method. In *2015 IEEE MTT-S International Microwave Symposium*, pages 1–3, 2015.

[26] Abdul-Lateef Haji-Ali, Fabio Nobile, and Raúl Tempone. Multi-index Monte Carlo: When sparsity meets sampling. *Numerische Mathematik*, 132(4):767–806, 2016.

[27] Tim Hesterberg. Weighted average importance sampling and defensive mixture distributions. *Technometrics*, 37(2):185–194, 1995.

[28] Janis Janusevskis and Rodolphe Le Riche. Simultaneous kriging-based estimation and optimization of mean response. *Journal of Global Optimization*, 55(2):313–336, 2013.

[29] S. S Lavenberg, T. L. Moeller, and P. D. Welch. The application of control variables to the simulation of closed queueing networks. In *Proceedings of the 9th Conference on Winter Simulation - Volume 1*, WSC '77, pages 152–154, Gaitherburg, Maryland, USA, 1977. Winter Simulation Conference.

[30] S. S. Lavenberg and P. D. Welch. A Perspective on the Use of Control Variables to Increase the Efficiency of Monte Carlo Simulations. *Management Science*, 27(3):322–335, 1981.

[31] Stephen S. Lavenberg, Thomas L. Moeller, and Peter D. Welch. Statistical Results on Control Variables with Application to Queueing Network Simulation. *Operations Research*, 30(1):182–202, 1982.
[17] Sarthak Garg, Neeraj Sood, and Costas D. Sarris. Uncertainty quantification of ray-tracing based wireless propagation models with a control variate-polynomial chaos expansion method. In *2015 IEEE International Symposium on Antennas and Propagation & USNC/URSI National Radio Science Meeting*, pages 1776–1777, 2015.

[18] G Geraci, M Eldred, and G Iaccarino. A multifidelity control variate approach for the multilevel Monte Carlo technique. Technical report, Stanford University, 2015.

[19] Gianluca Geraci, Michael S. Eldred, and Gianluca Iaccarino. A multifidelity multilevel Monte Carlo method for uncertainty propagation in aerospace applications. In *19th AIAA Non-Deterministic Approaches Conference, Grapevine, Texas, 2017*. American Institute of Aeronautics and Astronautics.

[20] Roger G. Ghanem and Pol D. Spanos. *Stochastic Finite Elements: A Spectral Approach*. Courier Corporation, 2003.

[21] Michael B. Giles. Multilevel Monte Carlo Path Simulation. *Operations Research*, 56(3):607–617, 2008.

[22] Michael B. Giles. Multilevel Monte Carlo methods. *Acta Numerica*, 24:259–328, 2015.

[23] Michael B. Giles and Benjamin J. Waterhouse. Multilevel quasi-Monte Carlo path simulation. *Advanced Financial Modelling, Radon Series on Computational and Applied Mathematics*, 8:165–181, 2009.

[24] Alex A. Gorodetsky, Gianluca Geraci, Michael S. Eldred, and John D. Jakeman. A generalized approximate control variate framework for multifidelity uncertainty quantification. *Journal of Computational Physics*, 408:109257, 2020.

[25] Zixi Gu and Costas D. Sarris. Multi-parametric uncertainty quantification with a hybrid monte-carlo / polynomial chaos expansion fdtd method. In *2015 IEEE MTT-S International Microwave Symposium*, pages 1–3, 2015.

[26] Abdul-Lateef Haji-Ali, Fabio Nobile, and Raúl Tempone. Multi-index Monte Carlo: When sparsity meets sampling. *Numerische Mathematik*, 132(4):767–806, 2016.

[27] Tim Hesterberg. Weighted average importance sampling and defensive mixture distributions. *Technometrics*, 37(2):185–194, 1995.

[28] Janis Janusevskis and Rodolphe Le Riche. Simultaneous kriging-based estimation and optimization of mean response. *Journal of Global Optimization*, 55(2):313–336, 2013.

[29] S. S Lavenberg, T. L. Moeller, and P. D. Welch. The application of control variables to the simulation of closed queueing networks. In *Proceedings of the 9th Conference on Winter Simulation - Volume 1*, WSC '77, pages 152–154, Gaitherburg, Maryland, USA, 1977. Winter Simulation Conference.

[30] S. S. Lavenberg and P. D. Welch. A Perspective on the Use of Control Variables to Increase the Efficiency of Monte Carlo Simulations. *Management Science*, 27(3):322–335, 1981.

[31] Stephen S. Lavenberg, Thomas L. Moeller, and Peter D. Welch. Statistical Results on Control Variables with Application to Queueing Network Simulation. *Operations Research*, 30(1):182–202, 1982.
[32] O. P. Le Maître and Omar M. Knio. *Spectral Methods for Uncertainty Quantification*. Scientific Computation. Springer Netherlands, Dordrecht, 2010.

[33] Christiane Lemieux. Control Variates. In *Wiley StatsRef: Statistics Reference Online*, pages 1–8. American Cancer Society, 2017.

[34] Nora Lüthen, Stefano Marelli, and Bruno Sudret. Sparse polynomial chaos expansions: Literature survey and benchmark. *SIAM/ASA Journal on Uncertainty Quantification*, 9(2):593–649, 2021.

[35] MD McKay, RJ Beckman, and WJ Conover. A comparison of three methods for selecting values of input variables in the analysis of output from a computer code. *Technometrics*, 21(2):239–245, 1979.

[36] Paul Mycek and Matthias De Lozzo. Multilevel Monte Carlo Covariance Estimation for the Computation of Sobol' Indices. *SIAM/ASA Journal on Uncertainty Quantification*, 7(4):1323–1348, 2019.

[37] Habib N. Najm. Uncertainty quantification and polynomial chaos techniques in computational fluid dynamics. *Annual Review of Fluid Mechanics*, 41(1):35–52, 2009.

[38] Barry L. Nelson. Control Variate Remedies. *Operations Research*, 38(6):974–992, 1990.

[39] Leo W. T. Ng and Karen E. Willcox. Multifidelity approaches for optimization under uncertainty. *International Journal for Numerical Methods in Engineering*, 100(10):746–772, 2014.

[40] F. Nobile, R. Tempone, and C. G. Webster. An Anisotropic Sparse Grid Stochastic Collocation Method for Partial Differential Equations with Random Input Data. *SIAM Journal on Numerical Analysis*, 46(5):2411–2442, 2008.

[41] J. Park and I. W. Sandberg. Universal Approximation Using Radial-Basis-Function Networks. *Neural Computation*, 3(2):246–257, 06 1991.

[42] Raghu Pasupathy, Bruce W. Schmeiser, Michael R. Taafe, and Jin Wang. Control-variate estimation using estimated control means. *IIE Transactions*, 44(5):381–385, 2012.

[43] Benjamin Peherstorfer, Karen Willcox, and Max Gunzburger. Optimal Model Management for Multifidelity Monte Carlo Estimation. *SIAM Journal on Scientific Computing*, 38(5):A3163–A3194, 2016.

[44] Benjamin Peherstorfer, Karen Willcox, and Max Gunzburger. Survey of Multifidelity Methods in Uncertainty Propagation, Inference, and Optimization. *SIAM Review*, 60(3):550–591, 2018.

[45] Carl Edward Rasmussen. *Gaussian Processes in Machine Learning*, pages 63–71. Springer Berlin Heidelberg, Berlin, Heidelberg, 2004.

[46] Matthew T. Reagan, Habib N. Najm, Roger G. Ghanem, and Omar M. Knio. Uncertainty quantification in reacting-flow simulations through non-intrusive spectral projection. *Combustion and Flame*, 132(3):545–555, 2003.

[47] Daniel Schaden. *Variance Reduction with Multilevel Estimators*. PhD thesis, Technische Universität München, Munich, 2021.
[32] O. P. Le Maître and Omar M. Knio. *Spectral Methods for Uncertainty Quantification*. Scientific Computation. Springer Netherlands, Dordrecht, 2010.

[33] Christiane Lemieux. Control Variates. In *Wiley StatsRef: Statistics Reference Online*, pages 1–8. American Cancer Society, 2017.

[34] Nora Lüthen, Stefano Marelli, and Bruno Sudret. Sparse polynomial chaos expansions: Literature survey and benchmark. *SIAM/ASA Journal on Uncertainty Quantification*, 9(2):593–649, 2021.

[35] MD McKay, RJ Beckman, and WJ Conover. A comparison of three methods for selecting values of input variables in the analysis of output from a computer code. *Technometrics*, 21(2):239–245, 1979.

[36] Paul Mycek and Matthias De Lozzo. Multilevel Monte Carlo Covariance Estimation for the Computation of Sobol' Indices. *SIAM/ASA Journal on Uncertainty Quantification*, 7(4):1323–1348, 2019.

[37] Habib N. Najm. Uncertainty quantification and polynomial chaos techniques in computational fluid dynamics. *Annual Review of Fluid Mechanics*, 41(1):35–52, 2009.

[38] Barry L. Nelson. Control Variate Remedies. *Operations Research*, 38(6):974–992, 1990.

[39] Leo W. T. Ng and Karen E. Willcox. Multifidelity approaches for optimization under uncertainty. *International Journal for Numerical Methods in Engineering*, 100(10):746–772, 2014.

[40] F. Nobile, R. Tempone, and C. G. Webster. An Anisotropic Sparse Grid Stochastic Collocation Method for Partial Differential Equations with Random Input Data. *SIAM Journal on Numerical Analysis*, 46(5):2411–2442, 2008.

[41] J. Park and I. W. Sandberg. Universal Approximation Using Radial-Basis-Function Networks. *Neural Computation*, 3(2):246–257, 06 1991.

[42] Raghu Pasupathy, Bruce W. Schmeiser, Michael R. Taafe, and Jin Wang. Control-variate estimation using estimated control means. *IIE Transactions*, 44(5):381–385, 2012.

[43] Benjamin Peherstorfer, Karen Willcox, and Max Gunzburger. Optimal Model Management for Multifidelity Monte Carlo Estimation. *SIAM Journal on Scientific Computing*, 38(5):A3163–A3194, 2016.

[44] Benjamin Peherstorfer, Karen Willcox, and Max Gunzburger. Survey of Multifidelity Methods in Uncertainty Propagation, Inference, and Optimization. *SIAM Review*, 60(3):550–591, 2018.

[45] Carl Edward Rasmussen. *Gaussian Processes in Machine Learning*, pages 63–71. Springer Berlin Heidelberg, Berlin, Heidelberg, 2004.

[46] Matthew T. Reagan, Habib N. Najm, Roger G. Ghanem, and Omar M. Knio. Uncertainty quantification in reacting-flow simulations through non-intrusive spectral projection. *Combustion and Flame*, 132(3):545–555, 2003.

[47] Daniel Schaden. *Variance Reduction with Multilevel Estimators*. PhD thesis, Technische Universität München, Munich, 2021.
[48] Daniel Schaden and Elisabeth Ullmann. On Multilevel Best Linear Unbiased Estimators. *SIAM/ASA Journal on Uncertainty Quantification*, 8(2):601–635, 2020.

[49] Daniel Schaden and Elisabeth Ullmann. Asymptotic Analysis of Multilevel Best Linear Unbiased Estimators. *SIAM/ASA Journal on Uncertainty Quantification*, 9(3):953–978, 2021.

[50] Ralph C. Smith. *Uncertainty Quantification: Theory, Implementation, and Applications*. Society for Industrial and Applied Mathematics, USA, 2013.

[51] Bruno Sudret. Global sensitivity analysis using polynomial chaos expansions. *Reliability Engineering & System Safety*, 93(7):964–979, 2008.

[52] T. J. Sullivan. *Introduction to Uncertainty Quantification*, volume 63. Springer Cham, 2015.

[53] Aretha Leonore Teckentrup. *Multilevel Monte Carlo Methods and Uncertainty Quantification*. PhD thesis, University of Bath, 2013.

[54] Rohit K. Tripathy and Ilias Bilionis. Deep uq: Learning deep neural network surrogate models for high dimensional uncertainty quantification. *Journal of Computational Physics*, 375:565–588, 2018.

[55] Norbert Wiener. The Homogeneous Chaos. *American Journal of Mathematics*, 60(4):897–936, 1938.

[56] Hang Yang, Yuji Fujii, K. W. Wang, and Alex A. Gorodetsky. Control Variate Polynomial Chaos: Optimal Fusion of Sampling and Surrogates for Multifidelity Uncertainty Quantification. *arXiv:2201.10745 [stat]*, 2022.
[48] Daniel Schaden and Elisabeth Ullmann. On Multilevel Best Linear Unbiased Estimators. *SIAM/ASA Journal on Uncertainty Quantification*, 8(2):601–635, 2020.

[49] Daniel Schaden and Elisabeth Ullmann. Asymptotic Analysis of Multilevel Best Linear Unbiased Estimators. *SIAM/ASA Journal on Uncertainty Quantification*, 9(3):953–978, 2021.

[50] Ralph C. Smith. *Uncertainty Quantification: Theory, Implementation, and Applications*. Society for Industrial and Applied Mathematics, USA, 2013.

[51] Bruno Sudret. Global sensitivity analysis using polynomial chaos expansions. *Reliability Engineering & System Safety*, 93(7):964–979, 2008.

[52] T. J. Sullivan. *Introduction to Uncertainty Quantification*, volume 63. Springer Cham, 2015.

[53] Aretha Leonore Teckentrup. *Multilevel Monte Carlo Methods and Uncertainty Quantification*. PhD thesis, University of Bath, 2013.

[54] Rohit K. Tripathy and Ilias Bilionis. Deep uq: Learning deep neural network surrogate models for high dimensional uncertainty quantification. *Journal of Computational Physics*, 375:565–588, 2018.

[55] Norbert Wiener. The Homogeneous Chaos. *American Journal of Mathematics*, 60(4):897–936, 1938.

[56] Hang Yang, Yuji Fujii, K. W. Wang, and Alex A. Gorodetsky. Control Variate Polynomial Chaos: Optimal Fusion of Sampling and Surrogates for Multifidelity Uncertainty Quantification. *arXiv:2201.10745 [stat]*, 2022.