# Multi-fidelity Monte Carlo¶

This tutorial builds on from Multi-level Monte Carlo and Approximate Control Variate Monte Carlo and introduces an approximate control variate estimator called Multi-fidelity Monte Carlo (MFMC). Unlike MLMC this method does not assume a strict ordering of models.

## Many Model MFMC¶

To derive the MFMC estimator first recall the two model ACV estimator

$Q_{0,\mathcal{Z}}^\mathrm{MF}=Q_{0,\mathcal{Z}_{0}} + \eta\left(Q_{1,\mathcal{Z}_{0}}-\mu_{1,\mathcal{Z}_{1}}\right)$

The MFMC estimator can be derived with the following recursive argument. Partition the samples assigned to each model such that $$\mathcal{Z}_\alpha=\mathcal{Z}_{\alpha,1}\cup\mathcal{Z}_{\alpha,2}$$ and $$\mathcal{Z}_{\alpha,1}\cap\mathcal{Z}_{\alpha,2}=\emptyset$$. That is the samples at the next lowest fidelity model are the samples used at all previous levels plus an additional independent set, i.e. $$\mathcal{Z}_{\alpha,1}=\mathcal{Z}_{\alpha-1}$$. See MFMC sampling strategy. Note the differences between this scheme and the MLMC scheme.

 MFMC sampling strategy¶ MLMC sampling strategy¶

Starting from two models we introduce the next low fidelity model in a way that reduces the variance of the estimate $$\mu_{\alpha}$$, i.e.

$\begin{split}Q_{0,\mathcal{Z}}^\mathrm{MF}&=Q_{0,\mathcal{Z}_{0}} + \eta_1\left(Q_{1,\mathcal{Z}_{1}}-\left(\mu_{1,\mathcal{Z}_{1}}+\eta_2\left(Q_{2,\mathcal{Z}_1}-\mu_{2,\mathcal{Z}_2}\right)\right)\right)\\ &=Q_{0,\mathcal{Z}_{0}} + \eta_1\left(Q_{1,\mathcal{Z}_{1}}-\mu_{1,\mathcal{Z}_{1}}\right)+\eta_1\eta_2\left(Q_{2,\mathcal{Z}_1}-\mu_{2,\mathcal{Z}_2}\right)\end{split}$

We repeat this process for all low fidelity models to obtain

$Q_{0,\mathcal{Z}}^\mathrm{MF}=Q_{0,\mathcal{Z}_{0}} + \sum_{\alpha=1}^M\eta_\alpha\left(Q_{\alpha,\mathcal{Z}_{\alpha,1}}-\mu_{\alpha,\mathcal{Z}_{\alpha}}\right)$

The optimal control variate weights for the MFMC estimator, which minimize the variance of the estimator, are $$\eta=(\eta_1,\ldots,\eta_M)^T$$, where for $$\alpha=1\ldots,M$$

$\eta_\alpha = -\frac{\covar{Q_0}{Q_\alpha}}{\var{Q_\alpha}}$

With this choice of weights the variance reduction obtained is given by

$\gamma = 1-\rho_1^2\left(\frac{r_1-1}{r_1}+\sum_{\alpha=2}^M \frac{r_\alpha-r_{\alpha-1}}{r_\alpha r_{\alpha-1}}\frac{\rho_\alpha^2}{\rho_1^2}\right)$

Let us use MFMC to estimate the mean of our high-fidelity model.

import numpy as np
import matplotlib.pyplot as plt
import pyapprox as pya
from functools import partial
from pyapprox.tests.test_control_variate_monte_carlo import \
TunableModelEnsemble, ShortColumnModelEnsemble, PolynomialModelEnsemble
np.random.seed(1)

short_column_model = ShortColumnModelEnsemble()
model_ensemble = pya.ModelEnsemble(
[short_column_model.m0,short_column_model.m1,short_column_model.m2])

costs = np.asarray([100, 50, 5])
target_cost = int(1e4)
idx = [0,1,2]
cov = short_column_model.get_covariance_matrix()[np.ix_(idx,idx)]

# define the sample allocation
nhf_samples,nsample_ratios = pya.allocate_samples_mfmc(
cov, costs, target_cost)[:2]
# generate sample sets
samples,values =pya.generate_samples_and_values_mfmc(
nhf_samples,nsample_ratios,model_ensemble,
short_column_model.generate_samples)
# compute mean using only hf data
hf_mean = values[0][0].mean()
# compute mlmc control variate weights
eta = pya.get_mfmc_control_variate_weights(cov)
# compute MLMC mean
mfmc_mean = pya.compute_approximate_control_variate_mean_estimate(eta,values)

# get the true mean of the high-fidelity model
true_mean = short_column_model.get_means()[0]
print('MLMC error',abs(mfmc_mean-true_mean))
print('MC error',abs(hf_mean-true_mean))


Out:

MLMC error 0.017508096391419592
MC error 0.009640551732770564


## Optimal Sample Allocation¶

Similarly to MLMC, the optimal number of samples that minimize the variance of the MFMC estimator can be determined analytically (see [PWGSIAM2016]). Recalling that $$C_\mathrm{tot}$$ is the total budget then the optimal number of high fidelity samples is

$N_0 = \frac{C_\mathrm{tot}}{\V{w}^T\V{r}}$

where $$\V{r}=[r_0,\ldots,r_M]^T$$ are the sample ratios defining the number of samples assigned to each level, i.e. $$N_\alpha=r_\alpha N_0$$. The sample ratios are

$r_\alpha=\left(\frac{w_0(\rho^2_{0,\alpha}-\rho^2_{0,\alpha+1})}{w_\alpha(1-\rho^2_{0,1})}\right)^{\frac{1}{2}}$

where $$\V{w}=[w_0,w_M]^T$$ are the relative costs of each model, and $$\rho_{j,k}$$ is the correlation between models $$j$$ and $$k$$.

Now lets us compare MC with MFMC using optimal sample allocations

poly_model = PolynomialModelEnsemble()
model_ensemble = pya.ModelEnsemble(poly_model.models)
cov = poly_model.get_covariance_matrix()
target_costs = np.array([1e1,1e2,1e3,1e4],dtype=int)
costs = np.asarray([10**-ii for ii in range(cov.shape[0])])
model_labels=[r'$f_0$',r'$f_1$',r'$f_2$',r'$f_3$',r'$f_4$']
variances, nsamples_history = [],[]
npilot_samples = 5
estimators = [pya.MC,pya.MFMC]
for target_cost in target_costs:
for estimator in estimators:
est = estimator(cov,costs)
nhf_samples,nsample_ratios = est.allocate_samples(target_cost)[:2]
variances.append(est.get_variance(nhf_samples,nsample_ratios))
nsamples_history.append(est.get_nsamples(nhf_samples,nsample_ratios))
variances = np.asarray(variances)
nsamples_history = np.asarray(nsamples_history)

fig,axs=plt.subplots(1,2,figsize=(2*8,6))
# plot sample allocation
pya.plot_acv_sample_allocation(nsamples_history[1::2],costs,model_labels,axs[1])
mfmc_total_costs = np.array(nsamples_history[1::2]).dot(costs)
mfmc_variances = variances[1::2]
axs[0].loglog(mfmc_total_costs,mfmc_variances,':',label=r'$\mathrm{MFMC}$')
mc_total_costs = np.array(nsamples_history[::2]).dot(costs)
mc_variances = variances[::2]
axs[0].loglog(mc_total_costs,mc_variances,label=r'$\mathrm{MC}$')
axs[0].set_ylim(axs[0].get_ylim()[0],1e-2)
_ = axs[0].legend()
#plt.show()
#fig # necessary for jupyter notebook to reshow plot in new cell


Total running time of the script: ( 0 minutes 1.487 seconds)

Gallery generated by Sphinx-Gallery