High-dimensional and higher-order multifidelity Monte Carlo estimators
Alessio Quaglino Simone Pezzuto Rolf Krause


An experiment in distributed Internet address management using blockchains
    [
    December 30, 2023
==========================================================================





Multifidelity Monte Carlo methods rely on a hierarchy of possibly less accurate but statistically correlated
simplified or reduced models, in order to accelerate the estimation of statistics of high-fidelity models without 
compromising the accuracy of the estimates. This approach has recently gained widespread attention in uncertainty quantification <cit.>. This is partly due to the availability of optimal strategies for the estimation of the expectation of scalar quantities-of-interest <cit.>. In practice, the optimal strategy for the expectation is also used for the estimation of variance and sensitivity indices <cit.>. However, a general strategy is still lacking for vector-valued problems, nonlinearly statistically-dependent models, and estimators for which a closed-form expression of the error is unavailable. The focus of the present work is to generalize the standard multifidelity estimators to the above cases. The proposed generalized estimators lead to an optimization problem that can be solved analytically and whose coefficients can be estimated numerically with few runs of the high- and low-fidelity models. We analyze the performance of the proposed approach on a selected number of experiments, with a particular focus on cardiac electrophysiology, where a hierarchy of physics-based low-fidelity models is readily available.






§ INTRODUCTION

Sampling methods for Uncertainty Quantification (UQ) require multiple evaluations of the problem at hand. 
Often, each sample requires computing the solution of partial differential equations (PDEs), which is generally
a computationally demanding task. For example, in cardiac electrophysiology, a single patient-tailored simulation 
based on the bidomain model can take thousands of node-hours
on a large cluster <cit.>. Therefore, UQ 
for such models at organ scale is currently unfeasible with plain Monte Carlo methods, although its importance has been highlighted in the literature <cit.>.
A more sophisticated approach is to perform most of the simulations on a
hierarchy of low-resolution models, yielding the so-called (geometric) multilevel Monte
Carlo (MLMC) method <cit.>, which extends the idea of control
variates methods <cit.>. The use of the
highest-resolution level guarantees convergence, while a significant portion of
the computational load is offset to the low-resolution hierarchy. While this
can offer a significant speedup of the estimation, it is not always possible to
coarsen complex geometries and to robustly transfer
information from one level to another, to an extent that the full
potential of this approach can be reached <cit.>.  This is
for example the case in the context of cardiac
modeling <cit.>. Moreover, the estimation of high-order moments requires
special care <cit.>.

An alternative approach, which does not rely on geometry coarsening, is to use model reduction strategies, 
such as projection-based or surrogate models <cit.>. 
Such an approximation is build upon the observation that in many applications it is clearly possible to
distinguish between offline and online phases of the workflow,
where the former typically can be very expensive and encompass everything that
can be precomputed in advance, e.g., training a surrogate model before any
patient-specific data becomes available, while the latter should be very cheap and
consist only of evaluations of the reduced model. Unfortunately, this approach has two
significant drawbacks. On the one hand, complex error estimates must be
provided to ensure that the approximation error is within acceptable bounds. On
the other hand, repeating the offline training is very expensive but may often be
necessary, e.g., in the case of patient-specific UQ. Alternative to that, physics-based 
reduction strategies are also possible: in cardiac electrophysiology, the activation map
can readily be computed via the eikonal equation, which provides a physiologically-meaningful solution <cit.>, 
strongly linked to the bidomain equation <cit.>. Error bounds, however,
are not trivially found, limiting the applicability of the methodology.

A more recent idea is to replace the above paradigm of model selection with that of model fusion, 
where low-fidelity models are allowed to be inaccurate, in the sense of not providing
an approximation within certain error bounds, provided that they exhibit some degree of
correlation <cit.> or even just statistical
dependence <cit.> to the high-fidelity one. This
approach has a twofold advantage. Firstly, it is the statistical
dependence, rather than the error bounds of coarse models, that is crucial to
ensure that propagating uncertainties via the low-fidelity models provides
useful information on the statistics of the high-fidelity quantity-of-interest.
Secondly, in the case of PDEs, low-fidelity models are not restricted to low-resolution
geometries, therefore solving them can be several orders of magnitude faster
than a coarse model. If a statistical dependence between the models is
present, in practice, most of the computational effort 
is taken by the estimation of the linear correlation coefficients or a nonlinear map between the outputs of the low- and the high-fidelity models.
For such reasons, multifidelity Monte Carlo (MFMC) methods
have become very popular over the last years and their applications span the
fields of UQ, inverse problems, and optimization <cit.>.

Several approaches can be used to combine the outputs of the model hierarchy. The use of a Bayesian regression,
mapping the low-fidelity output to the high-fidelity one, has been advocated in <cit.>. In <cit.>, a Gaussian Process Regression is considered in the context of cardiac electrophysiology. The Bayesian nature of this
approach automatically augments the estimate with full probability distributions and credible intervals,
obtained independently of the degree of statistical dependence between models and the number of samples
employed. This makes it a very suitable approach for scenarios where resources are scarce, such as 
clinical applications. In contrast to that, the multifidelity Monte Carlo method <cit.> 
provides point estimates, rather than probability distributions. However, these are equipped with convergence estimates 
for the asymptotic case, i.e., their error vanishes as the number of samples goes to infinity. This makes this approach 
particularly interesting when computational resources are large (but finite) and accurate
estimates are desired.

Current MFMC methods have been successfully used for estimating the expectation,
the variance, the Sobol indices, and rare events <cit.>. For the expectation
of a scalar quantity-of-interest, the optimal multifidelity estimator for a given computational budget can be found analytically <cit.>.
Similarly, <cit.> proposed an estimator for variance
and Sobol indices, with an analytical derivation only for the former.
Interestingly,
the authors also report that the optimal weights and sampling strategy for the
expectation multifidelity estimator
performs very well also for the variance and Sobol indices, although not being optimal.
This observation suggests that a multifidelity estimator for a general statistic
of interest may perform relatively well even when the weights and sampling strategy
are not optimal.

In this work, we explore the above idea, hence extending the aforementioned multifidelity estimators in several ways.
First, we generalize the scalar-valued multifidelity estimator to vector-valued quantity-of-interests,
proposing an optimal estimator in terms of mean squared error (MSE) that closely resembles the one
obtained in the scalar case. Second, we analyze a broad class of high-order statistics
of interest, for which a closed-form expression of the error is unavailable and therefore the optimal sampling strategy is unknown. To overcome this limitation, we suggest to approximate the MSE by estimating the correlations between the estimators, via sampling of the low- and high-fidelity models. Third, we suggest a strategy to leverage nonlinear statistical dependence between the models, which is closely connected to <cit.>.

The paper is organized as follows: Section 2 introduces the standard multifidelity Monte Carlo estimators of scalar quantities-of-interest. Section 3 describes the proposed generalized multifidelity estimators. Section 4 presents the results of the numerical experiments for the assessment of each estimator. Section 5 details a real-world application to cardiac electrophysiology. Section 6 illustrates an application to nonlinear incompressible elasticity.





§ STANDARD MULTIFIDELITY MONTE CARLO ESTIMATORS



Let upper-case letters denote random
variables (e.g. S, Ψ) and lower-case letters the values these take
(e.g. s, ψ), where ψ and ψ(s)
denote both, the value of the quantity-of-interest (QoI) as well as the function that provides the
QoI with respect to the random input s.
Let ψ^(i)(s) be the i-th model of the multifidelity hierarchy, with i=1,…,K. For ease of notations, we
denote by ψ_h(s) := ψ^(1)(s) and ψ_l(s) := ψ^(K)(s) the high- and lowest-fidelity models.

The main task in the context of uncertainty propagation is the computation
of integrals. In particular, the sought expectation of the QoI is obtained by:


    [Ψ_h] = ∫_Ωψ_h(s)  p(s) s,

where p(s) denotes the probability density of the random inputs.
Similarly, one can compute other statistical indicators
(e.g. moments, densities etc), by evaluating integrals of the form ∫
h(ψ_h(s))  p(s)ds. The method of choice in high-dimensional settings is
direct Monte Carlo, where the expectation of Equation (<ref>) is
estimated by:


    Ψ̂^(i)_m =1/m∑_j=1^m ψ^(i)(s_j),


where i=1 and the s_j are independent, identically distributed samples drawn from p(s). 
The convergence rate is independent of the dimension of 
S and the error decays as 𝒪(1/√(m)) 
<cit.>.
Therefore, in problems where each evaluation of ψ_h(s) poses a 
significant computational burden, the use of direct Monte Carlo can become 
impractical or even infeasible.



 §.§ Scalar multifidelity expectation estimator

At this stage, a multifidelity estimator is introduced


    Ψ̂_h = Ψ̂^(1)_m_1 + ∑_i=2^K α_i ( Ψ̂^(i)_m_i - Ψ̂^(i)_m_i-1).


The key aspect of the multifidelity estimate is that the magnitude of the discrepancy
|ψ_h(s) -  ψ^(i)(s)| does not affect its accuracy. More precisely, let 
the Pearson correlation coefficient of the i-th model be


    ρ_i,j := [ψ^(i),ψ^(j)]/σ_i σ_j


and its variance σ^2_i := [Ψ̂^(i)]. Then, it is possible to show <cit.> that the MSE of the multifidelity estimator is


    [Ψ̂_h] = σ^2_1/m_1 + ∑^K_i=2( 1/m_i-1 - 1/m_i) ( α_i^2 σ^2_i -2 α_i ρ_1,iσ_1 σ_i ).


Moreover, the cost of the estimator is 


    C(Ψ̂_h) = ∑^K_i=1 w_i m_i,


where w_i is the computational cost of the i-th model. Under some mild assumptions on ρ_i,j and w_i, it is possible to perform a constrained minimization of the variance (<ref>) with respect to α and m, for a prescribed computational cost C(Ψ̂_h) = B > 0, obtaining


    α^*_i := ρ_1,iσ_1/σ_i,    m^*_1=B/∑_i=1^K w_i r_i,    m^*_i=m^*_1 r_i,    i=2,…,K


where 


    r_i := √(w_1(ρ^2_1,i-ρ^2_1,i+1)/w_i(1-ρ^2_1,2)).


It follows that the ratio of errors (i.e. the variance reduction) of the standard and the multifidelity estimators, for the same computational budget, is 


    E(Ψ̂^*_h)/E(Ψ̂^(1)_n) = ( ∑_i=1^K √(w_i/w_1(ρ^2_1,i-ρ^2_1,i+1)))^2,


where n=B/w_1. Therefore, the variance is reduced if the costs w_i and the differences of the squared correlation coefficients are low.
This condition replaces the classical assumption on the deterministic pointwise errors |ψ_h(s) -  ψ^(i)(s)| of the low-fidelity models with respect to the high-fidelity one.




 §.§ Variance and Sobol indices estimators

The standard approach to develop multifidelity estimators for higher-order moments q(Ψ_h) is to substitute Ψ̂^(i)_m in (<ref>) with the unbiased single-level estimator q̂. For example, in the case of variance and Sobol indices, the use of the following estimators is proposed in <cit.>


    V̂   = 1/m-1∑_i=1^m ( ψ(s_i) - Ψ̂)^2, 
    V̂_̂ĵ   = 2/2m-1( ∑_i=1^m ψ(s_i) ψ(y^j_i) - m ( Ψ̂+ Ψ̂'/2)^2 + V̂+ V̂'/4), 
    T̂_̂ĵ   = 1/2m∑_i=1^m ( ψ(s'_i) - ψ(y^j_i) )^2,


where s' is a second set of samples and y^j_i is equal to s'_i, except for the j-th component which is equal to that of s_i.  However, obtaining error estimates in the form of (<ref>) is generally too difficult or impractical. Therefore, it is advocated in <cit.> to employ the optimal strategy (m^*,α^*) for the expectation also for the computation of other estimators. In the following section, we propose a different approach based on estimating the correlations between each pairs of summands (as computed via two different models) in the above estimators.






§ GENERALIZED MULTIFIDELITY MONTE CARLO ESTIMATORS




 §.§ Vector-valued multifidelity expectation estimator
 
We now consider the case of a vector-valued quantity-of-interest. In particular, we focus on a spatially-dependent ψ, i.e. ψ(x). It follows that the random variables and their estimators are also spatially dependent, i.e., Ψ(x) and Ψ̂(x). We consider the finite-dimensional case of computing them at a given set of points { x_j }_j=1^N. Given a computational budget B, the optimal scalar-valued estimator can therefore be computed at each point x_j. By prescribing a tolerance ϵ^2 for the integrated MSE of the optimal estimator, expressed as a function of the computational budget B, as given in <cit.>, we obtain


    ∫_Ωσ_1^2(x)/B( ∑_i=1^K √(w_i/w_1(ρ^2_1,i(x)-ρ^2_1,i+1(x))))^2 dx < ϵ^2.


By assuming a spatially-uniform budget B(x) = B, linearity of the integral can be exploited, so that the optimal budget B^* that satisfies the above equation is easily found. However, this will yield a spatially-dependent number of model evaluations m_i(x). Since it is not possible to perform a model evaluation only for a subset of the points, the following choices are made

    m_i = max_j=1, …, N m_i (x_j).


Unfortunately, in general this will lead to a large number of model evaluations. Instead, we define the integrated vector error as 

    ∫_Ω E(Ψ̂_h,x,α)  dx := ∑_j=1^N E(Ψ̂_h,x_j,α) | Ω_j |,


where | Ω_j | is the measure of the local integration area and 


    E(Ψ̂_h,x_j,α)    = σ^2_1(x_j)/m_1 + ∑^K_i=2( 1/m_i-1 - 1/m_i) ( α^2_i(x_j) σ^2_i(x_j) -2 α_i(x_j) ρ_1,i(x_j) σ_1(x_j) σ_i(x_j) ).


Therefore, the error of the vector-valued case can be seen as a weighted
average of the scalar-valued errors. Hence, the resulting estimator can be also
used for the generic spa­tially-in­de­pendent vector-valued case, given some
importance weights | Ω_j |. Setting the partial derivative of
E(Ψ̂_h,x_j,α) with respect to α_i(x_j) equal to zero, we
obtain that


    α^*_i(x_j) := ρ_1,i(x_j)σ_1(x_j)/σ_i(x_j),   ∀ j=1⋯ N.


It follows that


    E(Ψ̂_h,x_j,α^*)    = σ^2_1(x_j)/m_1 - ∑^K_i=2( 1/m_i-1 - 1/m_i) ρ^2_1,i(x_j) σ^2_1(x_j).


Therefore, the error sum can be written as


    ∑_j=1^N E(Ψ̂_h,x_j,α) = 1/m_1σ^2_1 - ∑^K_i=2( 1/m_i-1 - 1/m_i)  ρ_1,i^2 σ_1^2 ,


where


    σ^2_1 := ∑_j=1^N σ^2_1(x_j) | Ω_j |,   ρ^2_1,i := 1/σ^2_1( ∑_j=1^N ρ^2_1,i(x_j) σ^2_1(x_j) | Ω_j | )


Equation (<ref>) has the same form as for the scalar case, hence it leads to


    m^*_1=B/∑_i=1^K w_i r_i,    m^*_i=m^*_1 r_i,   r_i := √(w_1( ρ^2_1,i-ρ^2_1,i+1)/w_i(1-ρ^2_1,2)).





 §.§ Higher-order and sensitivity estimators

In general, given a single-level estimator q̂, its corresponding multifidelity estimator is defined as


    q̂_h = q̂^(1)_m_1 + ∑_i=2^K α_i ( q̂^(i)_m_i - q̂^(i)_m_i-1).


Regardless of the true expression of q̂, we implicitly define a function q via the relation


    q̂^(j)_m_i = 1/m_i∑_n=1^m_i q(Ψ^(j)(s_n)).


Under the assumption that the s_n are i.i.d. and m_i-1 < m_i, it can be shown that 


    [q̂^(j)_m_i,q̂^(k)_m_i-1-q̂^(k)_m_i]= 1/m_i( 1 - min( 1, m_i/m_i-1) ) c_jk =0,   ∀ (j,k),


where c_jk = [q(Ψ^(j)),q(Ψ^(k))]. Therefore, it possible to write the MSE of (<ref>) as


    [q̂_h] = [q̂^(1)_m_1] +    ∑_i=2^K α_i^2 ( [q̂^(i)_m_i-1] - [q̂^(i)_m_i]  ) 
    
    	- 2    ∑_i=2^K α_i ( [q̂^(1)_m_1,q̂^(i)_m_i-1] - [q̂^(1)_m_1,q̂^(i)_m_i] ).


While in general it is not possible to find analytic expressions for the above terms, they can be estimated numerically in a preprocessing step by the definition of q in equation (<ref>) as 


    [q̂^(j)_m_i,q̂^(k)_m_i-1] ≈1/m_i1/N-1∑_n=1^N ( q(Ψ^(j)(s_n)) - q̂^(j)_N ) (q(Ψ^(k)(s_n)) - q̂^(k)_N ),    m_i-1≤ m_i,


where N is a small number (i.e. N ≈ 10-20). By defining 


    ( σ^q_i )^2     := 1/N-1∑_n=1^N ( q(Ψ^(i)(s_n)) - q̂^(i)_N )^2, 
    ρ^q_1,i    := 1/σ^q_1 σ^q_i (N-1)∑_n=1^N ( q(Ψ^(1)(s_n)) - q̂^(1)_N ) (q(Ψ^(i)(s_n)) - q̂^(i)_N ),


we obtain an expression similar to (<ref>)


    [q̂_h] = (σ^q_1)^2/m_1 + ∑^K_i=2( 1/m_i-1 - 1/m_i) ( α_i^2 (σ^q_i)^2 -2 α_i ρ^q_1,iσ^q_1 σ^q_i ).


Therefore, the optimization procedure presented above leads to a similar result with σ^q and ρ^q in place of σ and ρ.




 §.§ Nonlinearly statistically-dependent estimators

We now consider the case where the correlation coefficient ρ_1,i between the models is low, but the two model exhibit a high degree of nonlinear statistical dependence. For ease of notation, we consider to estimate the expectation of a scalar Ψ, but the same procedure can be applied to estimate any statistics and to the vector-valued case, by combining it with the strategies from the previous sections. Let us consider an estimator of the form 


    ĝ_h = Ψ̂^(1)_m_1 + ∑_i=2^K α_i ( ĝ^(i)_m_i - ĝ^(i)_m_i-1),


where 


    ĝ^(j)_m_i = 1/m_i∑_n=1^m_i g(Ψ^(j)(s_n)).


The function g is an additional 1-dimensional model, for example a Gaussian Process regression, which is fitted using N_g samples. The idea is that for a strongly nonlinear statistical dependence, N_g+N_ρ≪ N, where N and N_ρ are the number of samples used to estimate the correlation for, respectively, the estimators (<ref>) and (<ref>). Following the same procedure as in the previous section, it follows that


    [ĝ_h] = (σ_1)^2/m_1 + ∑^K_i=2( 1/m_i-1 - 1/m_i) ( α_i^2 (σ^g_i)^2 -2 α_i ρ^g_1,iσ_1 σ^g_i ).


Therefore, the optimization procedure presented above leads to a similar result with σ^q and ρ^q in place of σ and ρ, where


    ( σ^g_i )^2     := 1/N-1∑_n=1^N ( g(Ψ^(i)(s_n)) - ĝ^(i)_N )^2, 
    ρ^g_1,i    := 1/σ_1 σ^g_i (N-1)∑_n=1^N ( Ψ^(1)(s_n) - Ψ̂^(1)_N ) (g(Ψ^(i)(s_n)) - ĝ^(i)_N ).







 §.§ Choosing the computational budget

In order to choose the budget, we express the integrated MSE of the optimal estimator as a function of the computational budget B, as given in <cit.>


    E(Ψ̂^*_h,B) = σ_1^2/B( ∑_i=1^K √(w_i/w_1(ρ^2_1,i-ρ^2_1,i+1)))^2 < ϵ^2,


from which it follows that 


    B^* = ( σ_1/ϵ∑_i=1^K √(w_i/w_1(ρ^2_1,i-ρ^2_1,i+1)))^2.










§ NUMERICAL ASSESSMENT

In this section, we test the proposed generalized estimators on two simple benchmark problems. The first test is a well-known problem in the context of multifidelity methods <cit.>. There, we consider the expectation, variance, and Sobol index estimators. In the second test, we focus on comparing the linear and nonlinearly statistically-dependent estimators. In Figures <ref> and <ref>, the two approaches are referred to as Linear and GPR, respectively. All of the multifidelity routines are implemented in a in-house, newly developed Python library called SLOTH. 







 §.§ Ishigami model

We consider the Ishigami model <cit.> and the low-fidelity models proposed in <cit.>


    f^(1)   =  sin(z_1) + 5 sin^2(z_2) +  1/10 z_3^4 sin(z_1),    z_i ∼𝒰(-π,π) 
    
    	f^(2)   =  sin(z_1) + 4.75 sin^2(z_2) + 1/10 z_3^4 sin(z_1), 
    
    	f^(3)   =  sin(z_1) + 3 sin^2(z_2) + 9/10 z_3^2 sin(z_1).


For this model, expectation and variance can be computed analytically as [f(z)]=2.5 and [f(z)] ≈ 10.845. The reference values of the Sobol indices can be found in <cit.>. We study numerically the errors in the proposed multifidelity estimators. For the expectation and variance estimators, we use both approaches of linear correlation and nonlinear statistical dependence. Moreover, we consider the Sobol indices as a vector-valued output, using the multidimensional estimator proposed in Section <ref>.



  §.§.§ Preprocessing budget

Figure <ref> shows the histograms of the estimated optimal m with 5000 replicates of the preprocessing step. In each replicate, m is obtained by analytic minimization (<ref>) of the error (<ref>), computed with either 20 (blue) or 100
(red) samples. Figure <ref> shows the histogram of the computed
estimators for 2000 replicates of the full analysis, again performed with
either 20 (blue) or 100 (red) preprocessing samples. From
Figure <ref>,
it is clear that the number of samples used for preprocessing has an
impact on the accuracy of the computed m (and α). This is due to
the accuracy of the variance and correlation estimates, which appear in the
formula for computing m and α.
Fortunately, the
second figure shows that this does not have an impact on the computed
estimators. In all of the examples, we have used 100 samples as preprocessing
budget.



  §.§.§ Expectation and variance estimation

Figure <ref> shows the computed error of the expectation and variance estimators using 100 replicates of the proposed estimators (<ref>) and (<ref>). We see that in this case, taking advantage of the nonlinear statistical dependence has no advantage over using only linear correlations. Table <ref> summarizes the optimal α and m for the case of linear correlations, both averaged over 100 runs. As discussed in <cit.>, the optimal budget allocations for expectation and variance are similar. Therefore, in principle there is no advantage in using a separate allocation for the variance. However, for this problem, we get a better approximation of the optimal α (0.8826 vs. 0.9455, as given in <cit.>). This improves the accuracy of the computed variance with no extra cost, since the samples used for estimating ρ and σ in the standard multifidelity estimators, which are required to estimate the expectation, can be re-used for the proposed generalized ones. Moreover, this is achieved without using any a priori error estimate and can therefore be used also in the cases where such bounds are not available.







  §.§.§ Sobol indices estimation

Here, we consider the main and total effect Sobol indices as vector-valued outputs of three components each. Therefore, we obtained a single optimal value m_k for each set of sensitivities, as discussed in Section <ref>, rather than a separate optimum for each index. Table <ref> summarizes the estimates for the optimal α and m, both averaged over 100 runs. We see that the optimal allocation differs between the two sets, and also with respect to the expectation and variance, shown in Table <ref>. Therefore, in this case the optimum for the expectation does not guarantee a nearly optimal allocation of the samples. Figure <ref> shows the computed error of the proposed multifidelity Sobol indices estimators using 100 runs of the multifidelity algorithm, using only linear correlations. In particular, the errors are similar for the main effect sensitivities, while they differ for the total effect ones. We have used uniform weights | Ω_j |=1 in (<ref>), however a different choice would have been possible if a higher precision had been required for one of the Sobol indices.








 §.§ Quintic model

We here consider the following model


    f^(1)   =  sin(z_1) +  sin^2(z_2) +  1/10 z_3^5,    z_i ∼𝒰(-π,π) 
    
    	f^(2)   =  sin(z_1) +  sin^2(z_2) +  2 z_3^3, 
    
    	f^(3)   =  sin(z_1) +  sin^2(z_2) + 20 z_3.


For this model, we estimated the exact expectation and variance using a standard Monte Carlo method. For the proposed multifidelity estimators of the expectation and variance, we use both approaches of linear correlation and statistical dependence, which are referred to in Figures <ref> as linear and GPR, respectively. For both strategies, we also give the optimal budget allocations in Tables <ref> and <ref>.








  §.§.§ Expectation and variance estimation

Tables <ref> and <ref> summarize the values ρ and the optimal α and m for the case of, respectively, linear correlations and nonlinear statistical dependence. Both are averaged over 100 runs and the nonlinear model is estimated using 100 preprocessing samples. We see that in particular, the use of the nonlinear model significantly increases the correlation, to the point that the second model becomes redundant for estimating the expectation. Similarly, for the variance the number of model evaluations required by the high-fidelity reduces from 24 to 14 when the nonlinear dependence is exploited. Figure <ref> shows the computed error of the expectation and variance estimators as a function of the budget, using 100 replicates of the multifidelity algorithm. From this chart, it is evident that taking advantage of the nonlinear statistical dependence, as done by the estimator (<ref>), provides a significant error reduction over the use of the more standard multifidelity estimator (<ref>), which is only based on linear correlations. 











§ APPLICATION TO CARDIAC ELECTROPHYSIOLOGY


We apply here the multifidelity methodology to a problem from
cardiac electrophysiology.  Briefly, the heart is an electrically
active organ, with an electric wave spreading throughout the tissue
thus to dictate the mechanical contraction.  The activation pattern,
defined as the first arrival time of the wave, is of clinical interest,
as an abnormal pattern leads to a mechanical dysfunction.
The propagation of the
traveling wave is anisotropic, being faster along the cardiac
fibers (usually 3- to 6-fold the transverse conduction velocity).
Fiber distribution is roughly known from histological studies <cit.>.
Patient-specific models incorporate the fiber distribution either
from diffusion-tensor MRI <cit.> or tailored to the anatomy
with a rule-based approach <cit.>. In both cases, epistemic
uncertainty is present.

In what follows, we consider the problem of quantifying the uncertainty
in the activation map, named (x), x∈Ω under the assumption
of randomly distributed cardiac fibers, (x).

The high-fidelity model is the monodomain equation <cit.>. This
is a reaction-diffusion equation, with nonlinear reaction, coupled
to a possibly large system of ODEs, encoding for the cellular
membrane model. The solution of the problem are the spatio-temporal
evolution of the transmembrane potential (x,t) and some auxiliary membrane
variables. The activation map is eventually obtained from the potential
as the time of first crossing of a given threshold:

    _h(x) = inf_t { t(x,t) ≥ V_m,thres}.

The monodomain equation reads as follows:

    β( C_m∂/∂ t
    + I_ion(,𝐳) - I_stim(x,t) ) = ∇·(∇),

where C_m is the membrane capacitance, β the surface-to-volume
ratio, 𝐳 the set of auxiliary variables, and  the
monodomain conductivity tensor, namely

    (x) = (x)  + ((x)-(x)) (x)⊗(x),

with  and  respectively the longitudinal and
transverse conductivity.


The low-fidelity model is the zeroth-order eikonal approximation of the
monodomain equation <cit.>.
The eikonal model is obtained by assuming that (x,t)=U(t-_l(x))
for some action potential shape U(ξ) and activation map _l,
and then neglecting the diffusion term. The model reads:

    {   θ/√(β)√(∇·∇) = 1,
       
       (x_i) = τ_i,
        i=1,…,K.
    .

The boundary conditions set the initial activation time to τ_i
at some given locations x_i∈Ω. The parameter θ is a
scaling factor linked to the membrane model. We remark that for the
purposes of the analysis the function U(ξ) is irrelevant.


In this study, we consider a patient-specific anatomy.
The electrophysiology parameters were already fitted to the patient data in
a previous study <cit.>, with good
correlation in terms of activation map and surface ECG. The numerical grid
is composed by approximately 20 millions voxels at full resolution (0.2 mm),
and 200 000 voxels at coarse resolution (1 mm).

For the monodomain simulation, we used Propag-5 software <cit.>.
Only the full resolution grid is considered, necessary to ensure
small numerical error <cit.>. 
The membrane model is from <cit.>. Activation is initiated with
a current stimulus in a 1 mm^3 region for 2 ms.
Regarding the eikonal model, we used our GPGPU implementation
described in <cit.>. We performed eikonal simulations at
both coarse and full resolution.






















In summary, the full hierarchy of models for the multifidelity framework includes:
the discretized monodomain model, f^(1)↦_h, the discretized
eikonal model with high spatial resolution, f^(2)↦_l,1,
and the eikonal model at coarse spatial resolution,
f^(3)↦_l,2.




 §.§ Random fiber distribution and sampling


The fiber distribution is modeled through a rule-based approach.
In details, the fiber orientation is locally determined by 3 angles, namely
α, ϕ and γ. These angles are related to the Euler's angles,
but conveniently adapted so that α is consistent with the definition
given by <cit.>: it is the angle between the local fiber
direction and the valves' plane of the heart.
Transmurally, the angle α varies cubically from
π/3 at the endocardium to -π/3 at the epicardium <cit.>

    α(r) = π/3(1-2r)^3,

where r measures the endocardium-to-epicardium distance, normalized to unity.
(r=0 corresponds to the endocardium, r=1 to the epicardium.)

The random fiber field is obtained with the above rule, but assuming that
the angle α is a Gaussian random (scalar field), named A(x), with
mean given by (<ref>) and correlation function as follows

    k(x,y) = σexp(-d(x,y)^2/2ℓ^2),

where σ is the standard deviation, ℓ the correlation length,
and d(x,y) the geodesic distance between x and y.

The sampling strategy is provided in details in <cit.>. In
summary, the random angle A(x) is approximated by a truncated Karhuhen-Loève
expansion:

    A(x) = A_0(x)
    + σ∑^K_i=1√(λ_i) Z_i φ_i(x),

with Z_i independent standard normal random variables, and
λ_i and φ_i being eigenvalues and eigenvectors
of the Hilbert-Schmidt operator associated to the kernel k(x,y).
The optimal truncation is obtained from a low-rank approximation
of the discretized correlation matrix
via a pivoted Cholesky decomposition <cit.>.
The distance function, required during the assembly of the
correlation matrix, is evaluated by solving an eikonal problem.







 §.§ Expectation estimation

Figure <ref> shows the correlations between the expectation estimators for the bidomain and eikonal models at the same mesh resolution h=0.02 cm (top) and at a coarser resolution h=0.1 cm (bottom). Figure <ref> shows the estimated error in the estimated expectation when only the bidomain model is used (top) and its reduction when the fine eikonal model is added (bottom). Figure <ref> (top) shows the computed expectation when all three models are used for the estimation. Table <ref> summarizes the optimal number of model evaluations, estimated RMSEs e_k with increasing number of models, and correlations ρ̅_k for a tolerance of ϵ=0.75 ms. As a comparison, standard MC would require 92 samples of the high-fidelity model in order to achieve the same accuracy. The estimated node time for standard and multifidelity MC is, respectively, 277h and 36h.





 §.§ Variance estimation


Figure <ref> shows the correlations between the variance estimators for the bidomain and eikonal models at the same mesh resolution h=0.02 cm (top) and at a coarser resolution h=0.1 cm (bottom). Figure <ref> shows the estimated error in the estimated standard deviation when only the bidomain model is used (top) and its reduction when the fine eikonal model is added (bottom). Figure <ref> (bottom) shows the computed standard deviation when all three models are used for the estimation. Table <ref> summarizes the optimal number of model evaluations, estimated RMSEs e_k with increasing number of models, and correlations ρ̅_k for a tolerance of ϵ=10.5 ms^2. As a comparison, standard MC would require 48 samples of the high-fidelity model in order to achieve the same accuracy. The estimated node time for standard and multifidelity MC is, respectively, 145h and 35h.






































§ APPLICATION TO NONLINEAR MECHANICS


The last application is inspired by to the indentation problem
of hyperelastic material, as described for instance in [REF]. Briefly, we
consider a 2×2×1 ^3 cube clamped at the bottom,
while top face is indented of 0.1 
in a 1 ^2 square region at its center.
We consider the finite deformation regime of a incompressible neo-Hookean
material with strain energy density function given as follows:

    𝒲(𝐅) = μ/2(𝐅:𝐅-3) - μln(J).

The deformation gradient tensor is defined as 𝐅=𝐈+∇𝐮,
where 𝐮Ω→ℝ^3 is the (unknown) displacement.
The shear modulus is set to vary in space following a spatially-correlated
random field, evaluated via Karhuhen-Loève expantion:

    μ(x) = μ_0 + σ∑_i=1^K √(λ_i) Z_i ϕ_i(x),     Z_i ∼𝒩(0,1).

Functions ϕ_i(x) and coefficients λ_i are, respectively,
eigenfunctions and eigenvalues of the Hilbert-Schmidt operator for
the squared-exponential correlation function, as in the previous section.
The correlation length is ℓ = 1/2√(2), while mean and
standard deviation of μ are set to 80 and 20.
Samples with non-positive μ(x) are discarded.

The aim of this experiment is to estimated the pointwise average of the
von Mises stress, which is related to the second invariant of the Cauchy
stress:

    S = μ/J√(32𝐁:𝐁),
        𝐁=𝐅𝐅^T,

given the uncertainty in the shear modulus.

The high-fidelity model is the Finite Element (FE) discretization of the problem
in mixed form with Taylor-Hood ℙ^2-ℙ^1 elements. We introduce
the pressure field p to enforce the incompressibility constraint J=1.
The discrete problem is then solved with the Newton's method and a direct solver
for the linearised system.

The low-fidelity problem is based on the compressible version of the
neo-Hookean material,

    𝒲_lf(𝐅) = μ/2(𝐅:𝐅-3) - μln(J)
    + λ/2(J-1)^2,

with λ = 10μ_0 = 800. We discretized the problem with linear
FE for the displacement. Since the Poisson's ratio is about 0.45, we do not expect
locking for this setup. In addition, the linear system within each Newton's iteration
is solved via algebraic multigrid.



 §.§ Expectation estimation of Von Mises stresses

Figure <ref> (left) shows the computed expectation when both the high- and low-fidelity models are used for the estimation. Figure <ref> (center) shows the correlations between the expectation estimators for the two models. Figure <ref> (right) shows the estimated error in the estimated expectation when both models are used. We note that the high concentration of the error along the displaced boundary would have significantly increased the number of required model evaluations, had each point been treated as a separate scalar-valued QoI. Table <ref> summarizes the optimal number of model evaluations, estimated RMSEs e_k with increasing number of models, and correlations ρ̅_k for a budget of B=2000 s. As a comparison, standard MC would require 293 samples of the high-fidelity model in order to achieve the same accuracy. The estimated CPU time for standard and multifidelity MC is, respectively, 27000 s and 2000 s.
























 §.§ Variance estimation of Von Mises stresses

Table <ref> (last three columns) summarizes the optimal number of model evaluations, estimated RMSEs e_k with increasing number of models, and correlations ρ̅_k for a budget of B=2000 s. As a comparison, standard MC would require 184 samples of the high-fidelity model in order to achieve the same accuracy. The estimated CPU time for standard and multifidelity MC is, respectively, 16000 s and 2000 s. The qualitative profiles of the estimated variance, the correlation, and the error, follow closely those of the expectation and are therefore not reported in the figures.
























§ CONCLUSIONS

In this work, we extended and generalized previously proposed multifidelity
Monte Carlo estimators for the expectation <cit.>, the variance, and sensitivity indices <cit.>. We tested our  methodology on selected benchmarks, showing the good performance and usability of the proposed estimators. Moreover, we detailed their application to a real-world example from cardiac electrophysiology, which is unapproachable by standard Monte Carlo techniques.

Rather than considering the low- and high-fidelity models as opposing, which should be selected according to time constraints or accuracy considerations, the multifidelity framework combines them and exploits their best features. By doing so, a full UQ analysis of the complex problems, such as the monodomain equation in cardiac electrophysiology, can be performed in a total time comparable to that of using only the low-fidelity model, without sacrificing the complexity of the high-fidelity solution.

Future works will focus on reduced-physics models for electrophysiology, incorporating the repolarization and the ECG into the multifidelity hierarchy. Moreover, we intend to study multifidelity estimators for inverse problems.




§ ACKNOWLEDGEMENTS

The authors acknowledge financial support by the Theo Rossi di Montelera Foundation, the Metis Foundation Sergio Mantegazza, the Fidinam Foundation, and the Horten Foundation to the Center for Computational Medicine in Cardiology. This work was also supported by grants from the Swiss National Supercomputing Centre (CSCS) under project ID s778.





