Add code to compute beta vectors for MLBLUE

Author: jdjakem

pyapprox/multifidelity/multilevelblue.py

@@ -98,6 +98,20 @@ def BLUE_Psi(Sigma, costs, reg_blue, nsamples_per_subset):
     return mat, submats
 
 
+def BLUE_betas(Sigma, asketch, reg_blue, nsamples_per_subset):
+    nmodels = Sigma.shape[0]
+    subsets = get_model_subsets(nmodels)
+    Psi = BLUE_Psi(Sigma, None, reg_blue, nsamples_per_subset)[0]
+    Psi_inv = np.linalg.inv(Psi)
+    betas = np.empty((len(subsets), nmodels))
+    for ii, subset in enumerate(subsets):
+        Sigma_inv = np.linalg.inv(Sigma[np.ix_(subset, subset)])
+        R = _restriction_matrix(nmodels, subset)
+        betas[ii] = np.linalg.multi_dot(
+            (R.T, Sigma_inv, R, Psi_inv, asketch))[:, 0]*nsamples_per_subset[ii]
+    return betas
+
+
 def BLUE_RHS(Sigma, values):
     """
     Parameters

pyapprox/multifidelity/tests/test_control_variate_monte_carlo.py

@@ -44,6 +44,7 @@
 from pyapprox.benchmarks.multifidelity_benchmarks import (
     ShortColumnModelEnsemble, TunableModelEnsemble, PolynomialModelEnsemble
 )
+from pyapprox.multifidelity.multilevelblue import BLUE_betas, BLUE_variance
 
 
 skiptest = unittest.skipIf(
@@ -911,6 +912,20 @@ def test_MLBLUE(self):
             estimator.nsamples_per_subset).astype(int)
         # rounded_target_cost = np.sum(
         #     estimator.nsamples_per_subset*subset_costs)
+
+        # test equivalent formulations for variance
+        betas = BLUE_betas(cov, asketch, 1e-10, estimator.nsamples_per_subset)
+        tmp = np.array(
+            [np.linalg.multi_dot((betas[kk].T, cov, betas[kk]))
+             for kk in range(betas.shape[0])])
+        tmp[estimator.nsamples_per_subset > 0] /= (
+            estimator.nsamples_per_subset[estimator.nsamples_per_subset > 0])
+        variance = tmp.sum()
+        assert np.allclose(
+            variance,
+            BLUE_variance(asketch, cov, None, 1e-10,
+                          estimator.nsamples_per_subset))
+
         values = estimator.generate_data(model.functions, variable)
         for ii in range(model.nmodels):
             asketch = np.zeros((costs.shape[0], 1))

tutorials/multi_fidelity/plot_many_model_approximate_control_variate_monte_carlo.py

@@ -233,14 +233,15 @@
 #
 #where each model is the function of a single uniform random variable defined on the unit interval :math:`[0,1]`.
 
+from pyapprox.util.configure_plots import mathrm_labels, mathrm_label
 plt.figure()
 benchmark = setup_benchmark("polynomial_ensemble")
 poly_model = benchmark.fun
 cov = poly_model.get_covariance_matrix()
 model_costs = np.asarray([10**-ii for ii in range(cov.shape[0])])
 nhf_samples = 10
 nsample_ratios_base = np.array([2, 4, 8, 16])
-cv_labels = [r'$\mathrm{OCV-1}$', r'$\mathrm{OCV-2}$', r'$\mathrm{OCV-4}$']
+cv_labels = mathrm_labels(["OCV-1", "OCV-2", "OCV-4"])
 cv_rsquared_funcs = [
     lambda cov: get_control_variate_rsquared(cov[:2, :2]),
     lambda cov: get_control_variate_rsquared(cov[:3, :3]),
@@ -253,7 +254,6 @@
 plt.axhline(y=1, linestyle='--', c='k')
 plt.text(xloc, 1, r'$\mathrm{MC}$', fontsize=16)
 
-from pyapprox.util.configure_plots import mathrm_labels
 acv_labels = mathrm_labels(["MLMC", "MFMC", "ACVMF"])
 estimators = [
     multifidelity.get_estimator("mlmc", cov, model_costs, poly_model.variable),
@@ -272,7 +272,7 @@
                  label=acv_labels[ii])
 plt.legend()
 plt.xlabel(r'$\log_2(r_i)-i$')
-_ = plt.ylabel(r'$\mathrm{Variance}$ $\mathrm{reduction}$ $\mathrm{ratio}$ $\gamma$')
+_ = plt.ylabel(mathrm_label("Variance reduction ratio ")+r"$\gamma$")
 
 #%%
 #As the theory suggests MLMC and MFMC use multiple models to increase the speed to which we converge to the optimal 2 model CV estimator OCV-2. These two approaches reduce the variance of the estimator more quickly than the ACV estimator, but cannot obtain the optimal variance reduction.

tutorials/multi_fidelity/plot_monte_carlo.py

@@ -14,17 +14,17 @@
 The mean squared error (MSE) of this estimator can be expressed as
 
 .. math::
-   
+
    \mean{\left(Q_{\alpha,N}-\mean{Q}\right)^2}&=\mean{\left(Q_{\alpha,N}-\mean{Q_{\alpha,N}}+\mean{Q_{\alpha,N}}-\mean{Q}\right)^2}\\
    &=\mean{\left(Q_{\alpha,N}-\mean{Q_{\alpha,N}}\right)^2}+\mean{\left(\mean{Q_{\alpha,N}}-\mean{Q}\right)^2}\\
    &\qquad\qquad+\mean{2\left(Q_{\alpha,N}-\mean{Q_{\alpha,N}}\right)\left(\mean{Q_{\alpha,N}}-\mean{Q}\right)}\\
    &=\var{Q_{\alpha,N}}+\left(\mean{Q_{\alpha,N}}-\mean{Q}\right)^2
-   
+
 Here we used that :math:`\mean{\left(Q_{\alpha,N}-\mean{Q_{\alpha,N}}\right)}=0` so the third term on the second line is zero. Now using
 
 .. math::
 
-   \var{Q_{\alpha,N}}=\var{N^{-1}\sum_{n=1}^N f^{(n)}_\alpha}=N^{-1}\sum_{n=1}^N \var{f^{(n)}_\alpha}=N^{-1}\var{Q_\alpha}
+   \var{Q_{\alpha,N}}=\var{N^{-1}\sum_{n=1}^N f^{(n)}_\alpha}=N^{-2}\sum_{n=1}^N \var{f^{(n)}_\alpha}=N^{-1}\var{Q_\alpha}
 
 yields
 
@@ -33,7 +33,7 @@
    \mean{\left(Q_{\alpha}-\mean{Q}\right)^2}=\underbrace{N^{-1}\var{Q_\alpha}}_{I}+\underbrace{\left(\mean{Q_{\alpha}}-\mean{Q}\right)^2}_{II}
 
 From this expression we can see that the MSE can be decomposed into two terms;
-a so called stochastic error (I) and a deterministic bias (II). The first term is the variance of the Monte Carlo estimator which comes from using a finite number of samples. The second term is due to using an approximation of :math:`f`. These two errors should be balanced, however in the vast majority of all MC analyses a single model :math:`f_\alpha` is used and the choice of :math:`\alpha`, e.g. mesh resolution, is made a priori without much concern for the balancing bias and variance. 
+a so called stochastic error (I) and a deterministic bias (II). The first term is the variance of the Monte Carlo estimator which comes from using a finite number of samples. The second term is due to using an approximation of :math:`f`. These two errors should be balanced, however in the vast majority of all MC analyses a single model :math:`f_\alpha` is used and the choice of :math:`\alpha`, e.g. mesh resolution, is made a priori without much concern for the balancing bias and variance.
 
 Given a fixed :math:`\alpha` the modelers only recourse to reducing the MSE is to reduce the variance of the estimator. In the following we plot the variance of the MC estimate of a simple algebraic function :math:`f_1` which belongs to an ensemble of models