Sampling Methods

This chapter introduces several fundamental concepts related to sampling methods. In particular, the statistical properties of the Monte Carlo estimator are discussed (Monte Carlo (MC)) and strategies for multilevel and multifidelity sampling are introduced within this context. Hereafter, multilevel refers to the possibility of exploiting distinct discretization levels (i.e. space/time resolution) within a single model form, whereas multifidelity involves the use of more than one model form. In Multifidelity Monte Carlo, we describe the multifidelity Monte Carlo and its single fidelity model version, the control variate Monte Carlo, that we align with multifidelity sampling, and in Multilevel Monte Carlo, we describe the multilevel Monte Carlo algorithm that we align with multilevel sampling. In A multilevel-multifidelity approach, we show that these two approaches can be combined to create multilevel-multifidelity sampling approaches.

Monte Carlo (MC)

Monte Carlo is a popular algorithm for stochastic simulations due to its simplicity, flexibility, and the provably convergent behavior that is independent of the number of input uncertainties. A quantity of interest $Q:\mathrm{\Xi }\to \mathbb{R}$, represented as a random variable (RV), can be introduced as a function of a random vector $\mathit{\xi }\in \mathrm{\Xi }\subset {\mathbb{R}}^d$. The goal of any MC simulation is computing statistics for $Q$, e.g. the expected value $\mathbb{E}\left[Q\right]$. The MC estimator ${\hat{Q}}_N^{MC}$ for $\mathbb{E}\left[Q\right]$ is defined as follows

$${\hat{Q}}_N^{MC}={\displaystyle \frac{1}{N}}\sum _{i=1}^N{Q}^{(i)},$$

where $Q^{(i)}=Q({\mathit{\xi }}^{(i)})$ and $N$ is used to indicate the number of realizations of the model.

The MC estimator is unbiased, i.e., its bias is zero and its convergence to the true statistics is $\mathcal{O}(N^{-1/2})$. Moreover, since each set of realizations for $Q$ is different, another crucial property of any estimator is its own variance:

(36)$$\mathbb{V}ar\left({\hat{Q}}_N^{MC}\right)={\displaystyle \frac{\mathbb{V}ar\left(Q\right)}{N}}.$$

Furthermore, it is possible to show, in the limit $N\to \mathrm{\infty }$, that the error $(\mathbb{E}\left[Q\right]-{\hat{Q}}_N^{MC})\sim \sqrt{{\displaystyle \frac{\mathbb{V}ar\left(Q\right)}{N}}}\mathcal{N}(0,1)$, where $\mathcal{N}(0,1)$ represents a standard normal RV. As a consequence, it is possible to define a confidence interval for the MC estimator which has an amplitude proportional to the standard deviation of the estimator. Indeed, the variance of the estimator plays a fundamental role in the quality of the numerical results: the reduction of the estimator variance correspond to an error reduction in the statistics.

Multifidelity Monte Carlo

A closer inspection of Eq. (36) indicates that only an increase in the number of simulations $N$ might reduce the overall variance, since $\mathbb{V}ar\left(Q\right)$ is an intrinsic property of the model under analysis. However, more sophisticated techniques have been proposed to accelerate the convergence of a MC simulation. For instance, an incomplete list of these techniques can include stratified sampling, importance sampling, Latin hypercube, deterministic Sobol’ sequences and control variates (see [Owe13]). In particular, the control variate approach, is based on the idea of replacing the RV $Q$ with one that has the same expected value, but with a smaller variance. The goal is to reduce the numerator in Eq. (36), and hence the value of the estimator variance without requiring a larger number of simulations. In a practical setting, the control variate makes use of an auxiliary RV $G=G(\mathit{\xi })$ for which the expected value $\mathbb{E}\left[G\right]$ is known. Indeed, the alternative estimator can be defined as

(37)$${\hat{Q}}_N^{MCCV}={\hat{Q}}_N^{MC}-\beta ({\hat{G}}_N^{MC}-\mathbb{E}\left[G\right]),\mathrm{where}\beta \in \mathbb{R}.$$

The MC control variate estimator ${\hat{Q}}_N^{MCCV}$ is unbiased, but its variance now has a more complex dependence not only on the $\mathbb{V}ar\left(Q\right)$, but also on $\mathbb{V}ar\left(G\right)$ and the covariance between $Q$ and $G$ since

$$\mathbb{V}ar\left({\hat{Q}}_N^{MCCV}\right)={\displaystyle \frac{1}{N}}(\mathbb{V}ar\left({\hat{Q}}_N^{MC}\right)+{\beta }^2\mathbb{V}ar\left({\hat{G}}_N^{MC}\right)-2\beta \mathrm{Cov}(Q,G)).$$

The parameter $\beta$ can be used to minimize the overall variance leading to

$$\beta ={\displaystyle \frac{\mathrm{Cov}(Q,G)}{\mathbb{V}ar\left(G\right)}},$$

for which the estimator variance follows as

$$\mathbb{V}ar\left({\hat{Q}}_N^{MCCV}\right)=\mathbb{V}ar\left({\hat{Q}}_N^{MC}\right)(1-{\rho }^2).$$

Therefore, the overall variance of the estimator ${\hat{Q}}_N^{MCCV}$ is proportional to the variance of the standard MC estimator ${\hat{Q}}_N^{MC}$ through a factor $1-{\rho }^2$ where $\rho ={\displaystyle \frac{\mathrm{Cov}(Q,G)}{\sqrt{\mathbb{V}ar\left(Q\right)\mathbb{V}ar\left(G\right)}}}$ is the Pearson correlation coefficient between $Q$ and $G$. Since $0<{\rho }^2<1$, the variance $\mathbb{V}ar\left({\hat{Q}}_N^{MCCV}\right)$ is always less than the corresponding $\mathbb{V}ar\left({\hat{Q}}_N^{MC}\right)$. The control variate technique can be seen as a very general approach to accelerate a MC simulation. The main step is to define a convenient control variate function which is cheap to evaluate and well correlated to the target function. For instance, function evaluations obtained through a different (coarse) resolution may be employed or even coming from a more crude physical/engineering approximation of the problem. A viable way of building a well correlated control variate is to rely on a low-fidelity model (i.e. a crude approximation of the model of interest) to estimate the control variate using estimated control means (see [NgW14, PTSW14] for more details). In this latter case, clearly the expected value of the low-fidelity model is not known and needs to be computed.

With a slight change in notation, it is possible to write

$${\hat{Q}}^{CVMC}=\hat{Q}+{\alpha }_1({\hat{Q}}_1-{\hat{\mu }}_1),$$

where $\hat{Q}$ represents the MC estimator for the high-fidelity model, ${\hat{Q}}_1$ the MC estimator for the low-fidelity model and ${\hat{\mu }}_1$ a different approximation for $\mathbb{E}[Q_1]$. If $N$ samples are used for approximating $\hat{Q}$ and ${\hat{Q}}_1$ and a total of $r_{1}N$ samples for the low-fidelity models are available, an optimal solution, which guarantees the best use of the low-fidelity resources, can be obtained following [NgW14] as

$${\alpha }_1=-{\rho }_1\sqrt{\frac{\mathbb{V}ar[Q]}{\mathbb{V}ar[Q_1]}}$$

$$r_1=\sqrt{\frac{\mathcal{C}}{{\mathcal{C}}_1}\frac{{\rho }_1^2}{1-{\rho }_1^2}},$$

where $\mathcal{C}$ and ${\mathcal{C}}_1$ represent the cost of evaluating the high- and low-fidelity models, respectively and ${\rho }_1$ is the correlation between the two models. This solution leads to the following expression for the estimator variance

$$\mathbb{V}ar[{\hat{Q}}^{CVMC}]=\mathbb{V}ar[\hat{Q}](1-\frac{r_1-1}{r_1}{\rho }_1^2),$$

which shows similarities with the variance of a control variate estimator with the only difference being the term $\frac{r_1-1}{r_1}$ that, by multiplying the correlation ${\rho }_1$, effectively penalizes the estimator due to the need for estimating the low-fidelity mean.

Another common case encountered in practice is the availability of more than a low-fidelity model. In this case, the multifidelity Monte Carlo can be extended following [PWG16, PWG18] as

$${\hat{Q}}^{MFMC}=\hat{Q}+\sum _{i=1}^M{\alpha }_i({\hat{Q}}_i-{\hat{\mu }}_i),$$

where ${\hat{Q}}_i$ represents the generic ith low-fidelity model.

The MFMC estimator is still unbiased (similarly to MC) and share similarities with CVMC; indeed one can recover CVMC directly from it. For each low-fidelity model we use $N_i{r}_i$ samples, as in the CVMC case, however for $i\ge 2$, the term $\widehat{Q_i}$ is approximated with exactly the same samples of the previous model, while each ${\hat{\mu }}_i$ is obtained by adding to this set a number of $(r_i-r_{i-1})N_i$ additional independent samples. Following [PWG16] the weights can be obtained as

(38)$${\alpha }_i=-{\rho }_i\sqrt{\frac{\mathbb{V}ar[Q]}{\mathbb{V}ar[Q_i]}}.$$

The optimal resource allocation problem is also obtainable in closed-form if, as demonstrated in [PWG16] the following conditions, for the models’ correlations and costs, hold

$$|{\rho }_1|>|{\rho }_2|>\cdots >|{\rho }_M|$$

$$\frac{{\mathcal{C}}_{i-1}}{{\mathcal{C}}_i}>\frac{{\rho }_{i-1}^2-{\rho }_i^2}{{\rho }_i^2-{\rho }_{i+1}^2},$$

leading to

$$r_i=\sqrt{\frac{\mathcal{C}}{{\mathcal{C}}_i}\frac{{\rho }_i^2-{\rho }_{i+1}^2}{1-{\rho }_1^2}}.$$

Multilevel Monte Carlo

In general engineering applications, the quantity of interest $Q$ is obtained as the result of the numerical solution of a partial partial differential equation (possibly a system of them). Therefore, the dependence on the physical $\mathbf{x}\in \mathrm{\Omega }\subset {\mathbb{R}}^n$ and/or temporal $t\in T\subset {\mathbb{R}}^{\mathbb{+}}$ coordinates should be included, hence $Q=Q(\mathbf{x},\mathit{\xi },t)$. A finite spatial/temporal resolution is always employed to numerically solve a PDE, implying the presence of a discretization error in addition to the stochastic error. The term discretization is applied generically with reference to either the spatial tessellation, the temporal resolution, or both (commonly, they are linked). For a generic tessellation with $M$ degrees-of-freedom (DOFs), the PDE solution of $Q$ is referred to as $Q_M$. Since $Q_M\to Q$ for $M\to \mathrm{\infty }$, then $\mathbb{E}\left[Q_M\right]\to \mathbb{E}\left[Q\right]$ for $M\to \mathrm{\infty }$ with a prescribed order of convergence. A MC estimator in presence of a finite spatial resolution and finite sampling is

$${\hat{Q}}_{M,N}^{MC}=\frac{1}{N}\sum _{i=1}^N{Q}_M^{(i)}$$

for which the mean square error (MSE) is

$$\mathbb{E}[({\hat{Q}}_{M,N}^{MC}-\mathbb{E}\left[Q\right]{)}^2]=N^{-1}\mathbb{V}ar\left(Q_M\right)+{\left(\mathbb{E}\left[Q_M-Q\right]\right)}^2,$$

where the first term represents the variance of the estimator, and the second term ${\left(\mathbb{E}[Q_M-Q]\right)}^2$ reflects the bias introduced by the (finite) spatial discretization. The two contributions appear to be independent of each other; accurate MC estimates can only be obtained by drawing the required $N$ number of simulations of $Q_M(\mathit{\xi })$ at a sufficiently fine resolution $M$. Since the numerical cost of a PDE is related to the number of DOFs of the tessellation, the total cost of a MC simulation for a PDE can easily become intractable for complex multi-physics applications that are computationally intensive.

Multilevel Monte Carlo for the mean

The multilevel Monte Carlo (MLMC) algorithm has been introduced, starting from the control variate idea, for situation in which additional discretization levels can be defined. The basic idea, borrowed from the multigrid approach, is to replace the evaluation of the statistics of $Q_M$ with a sequence of evaluations at coarser levels. If it is possible to define a sequence of discretization levels $\{M_{\ell }:\ell =0,\dots ,L\}$ with $M_0<M_1<\cdots <M_L\stackrel{\mathrm{def}}{=}M$, the expected value $\mathbb{E}\left[Q_M\right]$ can be decomposed, exploiting the linearity of the expected value operator as

$$\mathbb{E}\left[Q_M\right]=\mathbb{E}\left[Q_{M_0}\right]+\sum _{\ell =1}^L\mathbb{E}[Q_{M_{\ell }}-Q_{M_{\ell -1}}].$$

If the difference function $Y_{\ell }$ is defined according to

$$\begin{array}{r}{Y}_{\ell }=\{\begin{array}{rl}{Q}_{M_0}& \mathrm{if}\ell =0\\ Q_{M_{\ell }}-Q_{M_{\ell -1}}& \mathrm{if}0<\ell \le L,\end{array}\end{array}$$

the expected value $\mathbb{E}\left[Q_M\right]=\sum _{\ell =0}^L\mathbb{E}\left[Y_{\ell }\right]$. A multilevel MC estimator is obtained when a MC estimator is adopted independently for the evaluation of the expected value of $Y_{\ell }$ on each level. The resulting multilevel estimator ${\hat{Q}}_M^{\mathrm{ML}}$ is

$${\hat{Q}}_M^{\mathrm{ML}}=\sum _{\ell =0}^L{\hat{Y}}_{\ell ,N_{\ell }}^{\mathrm{MC}}=\sum _{\ell =0}^L\frac{1}{N_{\ell }}\sum _{i=1}^{N_{\ell }}{Y}_{\ell }^{(i)}.$$

Since the multilevel estimator is unbiased, the advantage of using this formulation is in its reduced estimator variance $\sum _{\ell =0}^L{N}_{\ell }^{-1}\mathbb{V}ar\left(Y_{\ell }\right)$: since $Q_M\to Q$, the difference function $Y_{\ell }\to 0$ as the level $\ell$ increases. Indeed, the corresponding number of samples $N_{\ell }$ required to resolve the variance associated with the $\ell$th level is expected to decrease with $\ell$.

The MLMC algorithm can be interpreted as a strategy to optimally allocate resources. If the total cost of the MLMC algorithm is written as

$$\mathcal{C}({\hat{Q}}_M^{ML})=\sum _{\ell =0}^L{N}_{\ell }{\mathcal{C}}_{\ell },$$

with ${\mathcal{C}}_{\ell }$ being the cost of the evaluation of $Y_{\ell }$ (involving either one or two discretization evaluations), then the following constrained minimization problem can be formulated where an equality constraint enforces a stochastic error (from MLMC estimator variance) equal to the residual bias error (${\epsilon }^2/2$)

(39)$$f(N_{\ell },\lambda )=\sum _{\ell =0}^L{N}_{\ell }{\mathcal{C}}_{\ell }+\lambda (\sum _{\ell =0}^L{N}_{\ell }^{-1}\mathbb{V}ar\left(Y_{\ell }\right)-{\epsilon }^2/2).$$

using a Lagrange multiplier $\lambda$. This equality constraint reflects a balance between the two contributions to MSE, reflecting the goal to not over-resolve one or the other. The result of the minimization is

$$N_{\ell }=\frac{2}{{\epsilon }^2}\left[\sum _{k=0}^L{\left(\mathbb{V}ar\left(Y_k\right){\mathcal{C}}_k\right)}^{1/2}\right]\sqrt{\frac{\mathbb{V}ar\left(Y_{\ell }\right)}{{\mathcal{C}}_{\ell }}},$$

defining an optimal sample allocation per discretization level.

MLMC extension to the variance

Despite the original introduction of the MLMC approach for the computation of the mean estimator in [Gil08, Gil15], it is possible to estimate higher-order moments with a MLMC sampling strategy, as for instance the variance.

A single level unbiased estimator for the variance of a generic QoI at the highest level $M_L$ of the hierarchy can be written as

(40)$$\mathbb{V}ar\left[Q_{M_L}\right]\approx \frac{1}{N_{M_L}-1}\sum _{i=1}^{N_{M_L}}{(Q_{M_L}^{(i)}-\mathbb{E}\left[Q_L\right])}^2.$$

The multilevel version of Eq. (40) can be obtained via a telescopic expansion in term of difference of estimators over subsequent levels. To simplify the notation and for simplicity of exposure from now on we only indicate the level, i.e. $M_{\ell }=\ell$.

The expansion is obtained by re-writing Eq. (40) as

$$\begin{array}{r}\begin{array}{rl}\mathbb{V}ar\left[Q_L\right]& \approx \frac{1}{N_L-1}\sum _{i=1}^{N_L}{(Q_L^{(i)}-\mathbb{E}\left[Q_L\right])}^2\\ & \approx \sum _{\ell =0}^L\frac{1}{N_{\ell }-1}({(Q_{\ell }^{(i)}-\mathbb{E}\left[Q_{\ell }\right])}^2-{(Q_{\ell -1}^{(i)}-\mathbb{E}\left[Q_{\ell -1}\right])}^2).\end{array}\end{array}$$

It is important here to note that since the estimators at the levels $\ell$ and $\ell -1$ are computed with the same number of samples both estimators use the factor $1/(N_{\ell }-1)$ to obtain their unbiased version. Moreover, each estimator is indeed written with respect to its own mean value, i.e. the mean value on its level, either $\ell$ or $\ell -1$. This last requirement leads to the computation of a local expected value estimator with respect to the same samples employed for the difference estimator. If we now denote with ${\hat{Q}}_{\ell ,2}$ the sampling estimator for the second order moment of the QoI $Q_{\ell }$ we can write

(41)$$\begin{array}{r}\mathbb{V}ar\left[Q_L\right]\approx {\hat{Q}}_{L,2}^{\mathrm{ML}}=\sum _{\ell =0}^L{\hat{Q}}_{\ell ,2}-{\hat{Q}}_{\ell -1,2},\end{array}$$

where

(42)$${\hat{Q}}_{\ell ,2}=\frac{1}{N_{\ell }-1}\sum _{i=1}^{N_{\ell }}{(Q_{\ell }^{(i)}-{\hat{Q}}_{\ell })}^2\mathrm{and}{\hat{Q}}_{\ell -1,2}=\frac{1}{N_{\ell }-1}\sum _{i=1}^{N_{\ell }}{(Q_{\ell -1}^{(i)}-{\hat{Q}}_{\ell -1})}^2.$$

Note that ${\hat{Q}}_{\ell ,2}$ and ${\hat{Q}}_{\ell -1,2}$ are explicitly sharing the same samples $N_{\ell }$.

For this estimator we are interested in minimizing its cost while also prescribing its variance as done for the expected value. This is accomplished by evaluating the variance of the multilevel variance estimator ${\hat{Q}}_{L,2}^{ML}$

$$\mathbb{V}ar\left[{\hat{Q}}_{L,2}^{\mathrm{ML}}\right]=\sum _{\ell =0}^L\mathbb{V}ar[{\hat{Q}}_{\ell ,2}-{\hat{Q}}_{\ell -1,2}]=\sum _{\ell =0}^L\mathbb{V}ar\left[{\hat{Q}}_{\ell ,2}\right]+\mathbb{V}ar\left[{\hat{Q}}_{\ell -1,2}\right]-2\mathbb{C}ov({\hat{Q}}_{\ell ,2},{\hat{Q}}_{\ell -1,2}),$$

where the covariance term is a result of the dependence described in (42).

The previous expression can be evaluated once the variance for the sample estimator of the second order order moment $\mathbb{V}ar\left[{\hat{Q}}_{\ell ,2}\right]$ and the covariance term $\mathbb{C}ov({\hat{Q}}_{\ell ,2},{\hat{Q}}_{\ell -1,2})$ are known. These terms can be evaluated as:

$$\mathbb{V}ar\left[{\hat{Q}}_{\ell ,2}\right]\approx \frac{1}{N_{\ell }}({\hat{Q}}_{\ell ,4}-\frac{N_{\ell }-3}{N_{\ell }-1}{\left({\hat{Q}}_{\ell ,2}\right)}^2),$$

where ${\hat{Q}}_{\ell ,4}$ denotes the sampling estimator for the fourth order central moment.

The expression for the covariance term is more involved and can be written as

$$\begin{array}{r}\begin{array}{rl}\mathbb{C}ov({\hat{Q}}_{\ell ,2},{\hat{Q}}_{\ell -1,2})& \approx \frac{1}{N_{\ell }}\mathbb{E}[{\hat{Q}}_{\ell ,2},{\hat{Q}}_{\ell -1,2}]\\ & +\frac{1}{N_{\ell }{N}_{\ell -1}}(\mathbb{E}{\left[Q_{\ell }{Q}_{\ell -1}\right]}^2-2\mathbb{E}\left[Q_{\ell }{Q}_{\ell -1}\right]\mathbb{E}\left[Q_{\ell }\right]\mathbb{E}\left[Q_{\ell -1}\right]+{\left(\mathbb{E}\left[Q_{\ell }\right]\mathbb{E}\left[Q_{\ell -1}\right]\right)}^2).\end{array}\end{array}$$

The first term of the previous expression is evaluated by estimating and combining several sampling moments as

$$\begin{array}{r}\begin{array}{rl}\mathbb{E}[{\hat{Q}}_{\ell ,2},{\hat{Q}}_{\ell -1,2}]& =\frac{1}{N_{\ell }}\left(\mathbb{E}\left[Q_{\ell }^2{Q}_{\ell -1}^2\right]\right)-\mathbb{E}\left[Q_{\ell }^2\right]\mathbb{E}\left[Q_{\ell -1}^2\right]-2\mathbb{E}\left[Q_{\ell -1}\right]\mathbb{E}\left[Q_{\ell }^2{Q}_{\ell -1}\right]\\ & +2\mathbb{E}\left[Q_{\ell -1}^2\right]\mathbb{E}\left[Q_{\ell }^2\right]-2\mathbb{E}\left[Q_{\ell }\right]\mathbb{E}\left[Q_{\ell }{Q}_{\ell -1}^2\right]+2\mathbb{E}{\left[Q_{\ell }\right]}^2\mathbb{E}\left[Q_{\ell -1}^2\right]\\ & +4\mathbb{E}\left[Q_{\ell }\right]\mathbb{E}\left[Q_{\ell -1}\right]\mathbb{E}\left[Q_{\ell }{Q}_{\ell -1}\right]-4\mathbb{E}{\left[Q_{\ell }\right]}^2\mathbb{E}{\left[Q_{\ell -1}\right]}^2.\end{array}\end{array}$$

It is important to note here that the previous expression can be computed only if several sampling estimators for product of the QoIs at levels $\ell$ and $\ell -1$ are available. These quantities are not required in the standard MLMC implementation for the mean and therefore for the estimation of the variance more data need to be stored to assemble the quantities on each level.

An optimization problem, similar to the one formulated for the mean in the previous section, can be written in the case of variance

(43)$$\begin{array}{r}\underset{N_{\ell }}{min}\sum _{\ell =0}^L{\mathcal{C}}_{\ell }{N}_{\ell }\mathrm{s}.\mathrm{t}.\mathbb{V}ar\left[{\hat{Q}}_{L,2}^{\mathrm{ML}}\right]={\epsilon }^2/2.\end{array}$$

This optimization problem can be solved in two different ways, namely an analytical approximation and by solving a non-linear optimization problem. The analytical approximation follows the approach described in [PKN17] and introduces a helper variable

$${\hat{V}}_{2,\ell }:=\mathbb{V}ar\left[{\hat{Q}}_{\ell ,2}\right]\cdot N_{\ell }.$$

Next, the following constrained minimization problem is formulated

(44)$$f(N_{\ell },\lambda )=\sum _{\ell =0}^L{N}_{\ell }{\mathcal{C}}_{\ell }+\lambda (\sum _{\ell =0}^L{N}_{\ell }^{-1}{\hat{V}}_{2,\ell }-{\epsilon }^2/2),$$

and a closed form solution is obtained

(45)$$N_{\ell }=\frac{2}{{\epsilon }^2}\left[\sum _{k=0}^L{\left({\hat{V}}_{2,k}{\mathcal{C}}_k\right)}^{1/2}\right]\sqrt{\frac{{\hat{V}}_{2,\ell }}{{\mathcal{C}}_{\ell }}},$$

similarly as for the expected value in (39).

The second approach uses numerical optimization directly on the non-linear optimization problem (43) to find an optimal sample allocation. Dakota uses OPTPP as the default optimizer and switches to NPSOL if it is available.

Both approaches for finding the optimal sample allocation when allocating for the variance are currently implemented in Dakota. The analytical solution is employed by default while the optimization is enabled using a keyword. We refer to the reference manual for a discussion of the keywords to select these different options.

MLMC extension to the standard deviation

The extension of MLMC for the standard deviation is slightly more complicated by the presence of the square root, which prevents a straightforward expansion over levels.

One possible way of obtaining a biased estimator for the standard deviation is

$${\hat{\sigma }}_L^{ML}=\sqrt{\sum _{\ell =0}^L{\hat{Q}}_{\ell ,2}-{\hat{Q}}_{\ell -1,2}}.$$

To estimate the variance of the standard deviation estimator, it is possible to leverage the result, derived in the previous section for the variance, and write the variance of the standard deviation as a function of the variance and its estimator variance. If we can estimate the variance ${\hat{Q}}_{L,2}$ and its estimator variance $\mathbb{V}ar\left[{\hat{Q}}_{L,2}\right]$, the variance for the standard deviation ${\hat{\sigma }}_L^{ML}$ can be approximated as

$$\mathbb{V}ar\left[{\hat{\sigma }}_L^{ML}\right]\approx \frac{1}{4{\hat{Q}}_{L,2}}\mathbb{V}ar\left[{\hat{Q}}_{L,2}\right].$$

Similarly to the variance case, a numerical optimization problem can be solved to obtain the sample allocation for the estimator of the standard deviation given a prescribed accuracy target.

MLMC extension to the scalarization function

Often, especially in the context of optimization, it is necessary to estimate statistics of a metric defined as a linear combination of mean and standard deviation of a QoI. A classical reliability measure $c^{ML}[Q]$ can be defined, for the quantity $Q$, starting from multilevel (ML) statistics, as

$$c_L^{ML}[Q]={\hat{Q}}_L^{ML}+\alpha {\hat{\sigma }}_L^{ML}.$$

To obtain the sample allocation, in the MLMC context, it is necessary to evaluate the variance of $c_L^{ML}[Q]$, which can be written as

$$\mathbb{V}ar[c_L^{ML}[Q]]=\mathbb{V}ar\left[{\hat{Q}}_L^{ML}\right]+{\alpha }^2\mathbb{V}ar\left[{\hat{\sigma }}_L^{ML}\right]+2\alpha \mathbb{C}ov[{\hat{Q}}_L^{ML},{\hat{\sigma }}_L^{ML}].$$

This expression requires, in addition to the already available terms $\mathbb{V}ar\left[{\hat{Q}}_L^{ML}\right]$ and $\mathbb{V}ar\left[{\hat{\sigma }}_L^{ML}\right]$, also the covariance term $\mathbb{C}ov[{\hat{Q}}_L^{ML},{\hat{\sigma }}_L^{ML}]$. This latter term can be written knowing that shared samples are only present on the same level

$$\begin{array}{r}\begin{array}{rl}\mathbb{C}ov[{\hat{Q}}_L^{ML},{\hat{\sigma }}_L^{ML}]& =\mathbb{C}ov[\sum _{\ell =0}^L{\hat{Q}}_{\ell }-{\hat{Q}}_{\ell -1},\sum _{\ell =0}^L{\hat{\sigma }}_{\ell }-{\hat{\sigma }}_{\ell -1}]\\ & =\sum _{\ell =0}^L\mathbb{C}ov[{\hat{Q}}_{\ell }-{\hat{Q}}_{\ell -1},{\hat{\sigma }}_{\ell }-{\hat{\sigma }}_{\ell -1}],\end{array}\end{array}$$

which leads to the need for evaluating the following four contributions

$$\mathbb{C}ov[{\hat{Q}}_{\ell }-{\hat{Q}}_{\ell -1},{\hat{\sigma }}_{\ell }-{\hat{\sigma }}_{\ell -1}]=\mathbb{C}ov[{\hat{Q}}_{\ell },{\hat{\sigma }}_{\ell }]-\mathbb{C}ov[{\hat{Q}}_{\ell },{\hat{\sigma }}_{\ell -1}]-\mathbb{C}ov[{\hat{Q}}_{\ell -1},{\hat{\sigma }}_{\ell }]+\mathbb{C}ov[{\hat{Q}}_{\ell -1},{\hat{\sigma }}_{\ell -1}].$$

In Dakota, we adopt the following approximation, for two arbitrary levels $\ell$ and $\kappa \in \{\ell -1,\ell ,\ell +1\}$

$$\rho [{\hat{Q}}_{\ell },{\hat{\sigma }}_{\kappa }]\approx \rho [{\hat{Q}}_{\ell },{\hat{Q}}_{\kappa ,2}]$$

(we indicate with ${\hat{Q}}_{\kappa ,2}$ the second central moment for $Q$ at the level $\kappa$), which corresponds to assuming that the correlation between expected value and variance is a good approximation of the correlation between the expected value and the standard deviation. This assumption is particularly convenient because it is possible to obtain in closed form the covariance between expected value and variance and, therefore, we can adopt the following approximation

$$\begin{array}{r}\begin{array}{c}\frac{\mathbb{C}ov[{\hat{Q}}_{\ell },{\hat{\sigma }}_{\kappa }]}{\sqrt{\mathbb{V}ar\left[{\hat{Q}}_{\ell }\right]\mathbb{V}ar\left[{\hat{\sigma }}_{\kappa }\right]}}\approx \frac{\mathbb{C}ov[{\hat{Q}}_{\ell },{\hat{Q}}_{\kappa ,2}]}{\sqrt{\mathbb{V}ar\left[{\hat{Q}}_{\ell }\right]\mathbb{V}ar\left[{\hat{Q}}_{\kappa ,2}\right]}}\\ \mathbb{C}ov[{\hat{Q}}_{\ell },{\hat{\sigma }}_{\kappa }]\approx \mathbb{C}ov[{\hat{Q}}_{\ell },{\hat{Q}}_{\kappa ,2}]\frac{\sqrt{\mathbb{V}ar\left[{\hat{\sigma }}_{\kappa }\right]}}{\sqrt{\mathbb{V}ar\left[{\hat{Q}}_{\kappa ,2}\right]}}.\end{array}\end{array}$$

Finally, we can derive the term $\mathbb{C}ov[{\hat{Q}}_{\ell },{\hat{Q}}_{\kappa ,2}]$ for all possible cases

$$\begin{array}{r}\mathbb{C}ov[{\hat{Q}}_{\ell },{\hat{Q}}_{\kappa ,2}]=\{\begin{array}{ll}\frac{1}{N_{\ell }}(\mathbb{E}\left[Q_{\ell }{Q}_{\kappa }^2\right]-\mathbb{E}\left[Q_{\ell }\right]\mathbb{E}\left[Q_{\kappa }^2\right]-2\mathbb{E}\left[Q_{\kappa }\right]\mathbb{E}\left[Q_{\ell }{Q}_{\kappa }\right]+2\mathbb{E}\left[Q_{\ell }\right]\mathbb{E}\left[Q_{\kappa }^2\right]),& \text{if }\kappa \ne \ell \\ \frac{{\hat{Q}}_{\ell ,3}}{N_{\ell }},& \text{if }\kappa =\ell .\end{array}\end{array}$$

Even for this case, the sample allocation problem can be solved by resorting to a numerical optimization given a prescribed target.

A multilevel-multifidelity approach

The MLMC approach described in Multilevel Monte Carlo can be related to a recursive control variate technique in that it seeks to reduce the variance of the target function in order to limit the sampling at high resolution. In addition, the difference function $Y_{\ell }$ for each level can itself be the target of an additional control variate (refer to Multifidelity Monte Carlo). A practical scenario is when not only different resolution levels are available (multilevel part), but also a cheaper computational model can be used (multifidelity part). The combined approach is a multilevel-multifidelity algorithm [FDKI17, GIE15, NT15], and in particular, a multilevel-control variate Monte Carlo sampling approach.

$Y_l$ correlations

If the target QoI can be generated from both a high-fidelity (HF) model and a cheaper, possibly biased low-fidelity (LF) model, it is possible to write the following estimator

(46)$$\mathbb{E}\left[Q_M^{\mathrm{HF}}\right]=\sum _{l=0}^{L_{\mathrm{HF}}}\mathbb{E}\left[Y_{\ell }^{\mathrm{HF}}\right]\approx \sum _{l=0}^{L_{\mathrm{HF}}}{\hat{Y}}_{\ell }^{\mathrm{HF}}=\sum _{l=0}^{L_{\mathrm{HF}}}{Y}_{\ell }^{\mathrm{HF},\star },$$

where

$$Y_{\ell }^{\mathrm{HF},\star }=Y_{\ell }^{\mathrm{HF}}+{\alpha }_{\ell }({\hat{Y}}_{\ell }^{\mathrm{LF}}-\mathbb{E}\left[Y_{\ell }^{\mathrm{LF}}\right]).$$

The estimator $Y_{\ell }^{\mathrm{HF},\star }$ is unbiased with respect to ${\hat{Y}}_{\ell }^{\mathrm{HF}}$, hence with respect to the true value $\mathbb{E}\left[Y_{\ell }^{\mathrm{HF}}\right]$. The control variate is obtained by means of the LF model realizations for which the expected value can be computed in two different ways: ${\hat{Y}}_{\ell }^{\mathrm{LF}}$ and $\mathbb{E}\left[Y_{\ell }^{\mathrm{LF}}\right]$. A MC estimator is employed for each term but the estimation of $\mathbb{E}\left[Y_{\ell }^{\mathrm{LF}}\right]$ is more resolved than ${\hat{Y}}_{\ell }^{\mathrm{LF}}$. For ${\hat{Y}}_{\ell }^{\mathrm{LF}}$, we choose the number of LF realizations to be equal to the number of HF realizations, $N_{\ell }^{\mathrm{HF}}$. For the more resolved $\mathbb{E}\left[Y_{\ell }^{\mathrm{LF}}\right]$, we augment with an additional and independent set of realizations ${\mathrm{\Delta }}_{\ell }^{\mathrm{LF}}$, hence $N_{\ell }^{\mathrm{LF}}=N_{\ell }^{\mathrm{HF}}+{\mathrm{\Delta }}_{\ell }^{\mathrm{LF}}$. The set ${\mathrm{\Delta }}_{\ell }^{\mathrm{LF}}$ is written, for convenience, as proportional to $N_{\ell }^{\mathrm{HF}}$ by means of a parameter $r_{\ell }\in {\mathbb{R}}_0^{+}$

$$N_{\ell }^{\mathrm{LF}}=N_{\ell }^{\mathrm{HF}}+{\mathrm{\Delta }}_{\ell }^{\mathrm{LF}}=N_{\ell }^{\mathrm{HF}}+r_{\ell }{N}_{\ell }^{\mathrm{HF}}=N_{\ell }^{\mathrm{HF}}(1+r_{\ell }).$$

The set of samples ${\mathrm{\Delta }}_{\ell }^{\mathrm{LF}}$ is independent of $N_{\ell }^{\mathrm{HF}}$, therefore the variance of the estimator can be written as (for further details see [GIE15])

(47)$$\begin{array}{r}\begin{array}{rl}\mathbb{V}ar\left({\hat{Q}}_M^{MLMF}\right)& =\sum _{l=0}^{L_{\mathrm{HF}}}({\displaystyle \frac{1}{N_{\ell }^{\mathrm{HF}}}}\mathbb{V}ar\left(Y_{\ell }^{\mathrm{HF}}\right)+{\displaystyle \frac{{\alpha }_{\ell }^2{r}_{\ell }}{(1+r_{\ell })N_{\ell }^{\mathrm{HF}}}}\mathbb{V}ar\left(Y_{\ell }^{\mathrm{HF}}\right)\\ & +2{\displaystyle \frac{{\alpha }_{\ell }{r}_{\ell }^2}{(1+r_{\ell })N_{\ell }^{\mathrm{HF}}}}{\rho }_{\ell }\sqrt{\mathbb{V}ar\left(Y_{\ell }^{\mathrm{HF}}\right)\mathbb{V}ar\left(Y_{\ell }^{\mathrm{LF}}\right)}),\end{array}\end{array}$$

The Pearson’s correlation coefficient between the HF and LF models is indicated by ${\rho }_{\ell }$ in the previous equations. Assuming the vector $r_{\ell }$ as a parameter, the variance is minimized per level, mimicking the standard control variate approach, and thus obtaining the optimal coefficient as ${\alpha }_{\ell }=-{\rho }_{\ell }\sqrt{{\displaystyle \frac{\mathbb{V}ar\left(Y_{\ell }^{\mathrm{HF}}\right)}{\mathbb{V}ar\left(Y_{\ell }^{\mathrm{LF}}\right)}}}$. By making use of the optimal coefficient ${\alpha }_{\ell }$, it is possible to show that the variance $\mathbb{V}ar\left(Y_{\ell }^{\mathrm{HF},\star }\right)$ is proportional to the variance $\mathbb{V}ar\left(Y_{\ell }^{\mathrm{HF}}\right)$ through a factor ${\mathrm{\Lambda }}_{\ell }(r_{\ell })$, which is an explicit function of the ratio $r_{\ell }$:

(48)$$\begin{array}{r}\begin{array}{rl}\mathbb{V}ar\left({\hat{Q}}_M^{MLMF}\right)& =\sum _{l=0}^{L_{\mathrm{HF}}}{\displaystyle \frac{1}{N_{\ell }^{\mathrm{HF}}}}\mathbb{V}ar\left(Y_{\ell }^{\mathrm{HF}}\right){\mathrm{\Lambda }}_{\ell }(r_{\ell })\mathrm{where}\\ {\mathrm{\Lambda }}_{\ell }(r_{\ell })& =(1-{\displaystyle \frac{r_{\ell }}{1+r_{\ell }}}{\rho }_{\ell }^2).\end{array}\end{array}$$

Note that ${\mathrm{\Lambda }}_{\ell }(r_{\ell })$ represents a penalty with respect to the classical control variate approach presented in Multifidelity Monte Carlo, which stems from the need to evaluate the unknown function $\mathbb{E}\left[Y_{\ell }^{\mathrm{LF}}\right]$. However, the ratio $r_{\ell }/(r_{\ell }+1)$ is dependent on the additional number of LF evaluations ${\mathrm{\Delta }}_{\ell }^{\mathrm{LF}}$, hence it is fair to assume that it can be made very close to unity by choosing an affordably large $r_{\ell }$, i.e., ${\mathrm{\Delta }}_{\ell }^{\mathrm{LF}}>>N_{\ell }^{\mathrm{HF}}$.

The optimal sample allocation is determined taking into account the relative cost between the HF and LF models and their correlation (per level). In particular the optimization problem introduced in Eq. (39) is replaced by

$${\mathrm{argmin}}_{N_{\ell }^{\mathrm{HF}},r_{\ell }}(\mathcal{L}),\mathrm{where}\mathcal{L}=\sum _{\ell =0}^{L_{\mathrm{HF}}}{N}_{\ell }^{\mathrm{HF}}{\mathcal{C}}_{\ell }^{\mathrm{eq}}+\lambda (\sum _{\ell =0}^{L_{\mathrm{HF}}}{\displaystyle \frac{1}{N_{\ell }^{\mathrm{HF}}}}\mathbb{V}ar\left(Y_{\ell }^{\mathrm{HF}}\right){\mathrm{\Lambda }}_{\ell }(r_{\ell })-{\epsilon }^2/2),$$

where the optimal allocation is obtained as well as the optimal ratio $r_{\ell }$. The cost per level includes now the sum of the HF and LF realization cost, therefore it can be expressed as ${\mathcal{C}}_{\ell }^{\mathrm{eq}}={\mathcal{C}}_{\ell }^{\mathrm{HF}}+{\mathcal{C}}_{\ell }^{\mathrm{LF}}(1+r_{\ell })$.

If the cost ratio between the HF and LF model is $w_{\ell }={\mathcal{C}}_{\ell }^{\mathrm{HF}}/{\mathcal{C}}_{\ell }^{\mathrm{LF}}$ then the optimal ratio is

$$r_{\ell }^{\star }=-1+\sqrt{{\displaystyle \frac{{\rho }_{\ell }^2}{1-{\rho }_{\ell }^2}}{w}_{\ell }},$$

and the optimal allocation is

$$\begin{array}{rl}{N}_{\ell }^{\mathrm{HF},\star }& =\frac{2}{{\epsilon }^2}[\sum _{k=0}^{L_{\mathrm{HF}}}{\left({\displaystyle \frac{\mathbb{V}ar\left(Y_k^{\mathrm{HF}}\right){\mathcal{C}}_k^{\mathrm{HF}}}{1-{\rho }_{\ell }^2}}\right)}^{1/2}{\mathrm{\Lambda }}_k(r_k^{\star })]\sqrt{(1-{\rho }_{\ell }^2)\frac{\mathbb{V}ar\left(Y_{\ell }^{\mathrm{HF}}\right)}{{\mathcal{C}}_{\ell }^{\mathrm{HF}}}}.\end{array}$$

It is clear that the efficiency of the algorithm is related not only to the efficiency of the LF model, i.e. how fast a simulation runs with respect to the HF model, but also to the correlation between the LF and HF model.

$Q_l$ correlations

A potential refinement of the previous approach [GEI17] consists in exploiting the QoI on each pair of levels, $\ell$ and $\ell -1$, to build a more correlated LF function. For instance, it is possible to use

$${\mathring{Y}}_{\ell }^{\mathrm{LF}}={\gamma }_{\ell }{Q}_{\ell }^{\mathrm{LF}}-Q_{\ell -1}^{\mathrm{LF}}$$

and maximize the correlation between $Y_{\ell }^{\mathrm{HF}}$ and ${\mathring{Y}}_{\ell }^{\mathrm{LF}}$ through the coefficient ${\gamma }_{\ell }$.

Formally the two formulations are completely equivalent if $Y_{\ell }^{\mathrm{LF}}$ is replaced with ${\mathring{Y}}_{\ell }^{\mathrm{LF}}$ in Eq. (46) and they can be linked through the two ratios

$$\begin{array}{r}\begin{array}{rl}{\theta }_{\ell }& ={\displaystyle \frac{\mathrm{Cov}(Y_{\ell }^{\mathrm{HF}},{\mathring{Y}}_{\ell }^{\mathrm{LF}})}{\mathrm{Cov}(Y_{\ell }^{\mathrm{HF}},Y_{\ell }^{\mathrm{LF}})}}\\ {\tau }_{\ell }& ={\displaystyle \frac{\mathbb{V}ar\left({\mathring{Y}}_{\ell }^{\mathrm{LF}}\right)}{\mathbb{V}ar\left(Y_{\ell }^{\mathrm{LF}}\right)}},\end{array}\end{array}$$

obtaining the following variance for the estimator

$$\mathbb{V}ar\left({\hat{Q}}_M^{MLMF}\right)={\displaystyle \frac{1}{N_{\ell }^{\mathrm{HF}}}}\mathbb{V}ar\left(Y_{\ell }^{\mathrm{HF}}\right)(1-{\displaystyle \frac{r_{\ell }}{1+r_{\ell }}}{\rho }_{\ell }^2{\displaystyle \frac{{\theta }_{\ell }^2}{{\tau }_{\ell }}}).$$

Therefore, a way to increase the variance reduction is to maximize the ratio ${\displaystyle \frac{{\theta }_{\ell }^2}{{\tau }_{\ell }}}$ with respect to the parameter ${\gamma }_{\ell }$. It is possible to solve analytically this maximization problem obtaining

$${\gamma }_{\ell }^{\star }={\displaystyle \frac{\mathrm{Cov}(Y_{\ell }^{\mathrm{HF}},Q_{\ell -1}^{\mathrm{LF}})\mathrm{Cov}(Q_{\ell }^{\mathrm{LF}},Q_{\ell -1}^{\mathrm{LF}})-\mathbb{V}ar\left(Q_{\ell -1}^{\mathrm{LF}}\right)\mathrm{Cov}(Y_{\ell }^{\mathrm{HF}},Q_{\ell }^{\mathrm{LF}})}{\mathbb{V}ar\left(Q_{\ell }^{\mathrm{LF}}\right)\mathrm{Cov}(Y_{\ell }^{\mathrm{HF}},Q_{\ell -1}^{\mathrm{LF}})-\mathrm{Cov}(Y_{\ell }^{\mathrm{HF}},Q_{\ell }^{\mathrm{LF}})\mathrm{Cov}(Q_{\ell }^{\mathrm{LF}},Q_{\ell -1}^{\mathrm{LF}})}}.$$

The resulting optimal allocation of samples across levels and model forms is given by

$$\begin{array}{r}\begin{array}{rl}{r}_{\ell }^{\star }& =-1+\sqrt{{\displaystyle \frac{{\rho }_l^2{\displaystyle \frac{{\theta }_{\ell }^2}{{\tau }_{\ell }}}}{1-{\rho }_{\ell }^2{\displaystyle \frac{{\theta }_{\ell }^2}{{\tau }_{\ell }}}}}{w}_{\ell }},\mathrm{where}{w}_{\ell }={\mathcal{C}}_{\ell }^{\mathrm{HF}}/{\mathcal{C}}_{\ell }^{\mathrm{LF}}\\ {\mathrm{\Lambda }}_{\ell }& =1-{\rho }_{\ell }^2{\displaystyle \frac{{\theta }_{\ell }^2}{{\tau }_{\ell }}}{\displaystyle \frac{r_{\ell }^{\star }}{1+r_{\ell }^{\star }}}\\ N_{\ell }^{\mathrm{HF},\star }& =\frac{2}{{\epsilon }^2}[\sum _{k=0}^{L_{\mathrm{HF}}}{\left({\displaystyle \frac{\mathbb{V}ar\left(Y_k^{\mathrm{HF}}\right){\mathcal{C}}_k^{\mathrm{HF}}}{1-{\rho }_{\ell }^2{\displaystyle \frac{{\theta }_{\ell }^2}{{\tau }_{\ell }}}}}\right)}^{1/2}{\mathrm{\Lambda }}_k(r_k^{\star })]\sqrt{(1-{\rho }_{\ell }^2{\displaystyle \frac{{\theta }_{\ell }^2}{{\tau }_{\ell }}})\frac{\mathbb{V}ar\left(Y_{\ell }^{\mathrm{HF}}\right)}{{\mathcal{C}}_{\ell }^{\mathrm{HF}}}}\end{array}\end{array}$$