1 Introduction and Motivation
Simulation of many complex systems, such as subsurface flows in porous media [1, 2] or reaction initiation in heterogeneous explosives [3], is complicated by a lack of information about key properties such as permeability or initial porosity. Uncertainty in the medium’s properties or initial state propagates into uncertainty in predicted quantities of interest (QoIs), such as mass flow rate or the material’s temperature.
Probabilistic methods, which treat an uncertain input or initial state of the system as a random variable, provide a natural venue to quantify predictive uncertainty in a QoI. These techniques render the QoI random as well, i.e., it takes on values that are distributed according to some probability density function (PDF). Such approaches include stochastic finite element methods (FEMs), which characterize the random parameter fields in terms of a finite set of random variables, e.g., via a spectral representation or a KarhunenLoève expansion. This finite set of random variables defines a finitedimensional outcome space on which the solution to the resulting stochastic partial differential equation (PDE) is defined. Examples of stochastic FEMs include stochastic Galerkin, which expands the solution of a stochastic PDE in terms of orthogonal basis functions, and stochastic collocation, which samples the random parameters at predetermined values or “nodes”
[4]. While such methods perform well when the number of stochastic parameters (aka “stochastic dimension”) is low and these parameters exhibit a long correlation length, for many stochastic degrees of freedom and short correlation lengths their performance decreases dramatically
[5]. In addition, since nonlinearity degrades the solution’s regularity in the outcome space, stochastic FEMs also struggle with solving highly nonlinear problems [6]. Another class of probabilistic techniques involves the derivation of deterministic equations for the statistical moments
[7, 8] or PDF [9, 10]of the QoI. While these methods do not suffer from the “curse of dimensionality”, they require a closure approximation for the derived moment or PDF equations. Such closures often involve perturbation expansions of relevant quantities into series in the powers of the random parameters’ variances, which limits the applicability of such methods to parameters with low coefficients of variation.
Monte Carlo simulations (MC) [11]
remain the most robust and straightforward way to solve PDEs with random parameters/initial conditions. The method samples the random variables from their distribution, solves the deterministic PDE for each realization, and computes the resulting statistics of the QoI. While combining a nonintrusive character with a convergence that is independent of the stochastic dimension, MC converges very slowly: the standard deviation of the MC estimator for the QoI’s expectation value is inversely proportional to
where is the number of realizations [11]. While this drawback spurred the development of alternative probabilistic methods such as those listed above, efforts to combine MC with the multigrid concept by Heinrich [12, 13] and later by Giles [14] have sparked renewed interest in MC under the form of the multilevel Monte Carlo (MLMC) method. MLMC aims to achieve the same solution error as MC but at a lower computational cost by correcting realizations on a coarse spatial grid with sampling at finer levels of discretization. By sampling predominantly at the coarsest levels where samples are cheaper to compute, it aims to outperform MC which only performs realizations on the finest spatial grid. The related technique of multifidelity Monte Carlo (also referred to as solverbased MLMC as a opposed to traditional gridbased MLMC) generalizes this approach by using models of varying fidelities and speeds on the different levels [15, 16, 17].Most work on MLMC has been aimed at computing estimators for expectation values and variances of QoIs (see, e.g., [18, 19, 20, 21]). However, the information contained in these first two moments is insufficient to understand, e.g., the probability of rare events. This task requires knowledge of the QoI’s full PDF, or equivalently, its cumulative distribution function (CDF) [22, 23, 24]. Assuming the QoI to be a function of a continuous random input variable , i.e., , where is a measurable function with the sample space, the CDF of at a point is given by
(1) 
where the indicator function is defined as
(2) 
Application of MLMC to the estimation of distributions only started a few years ago. Giles et al. [25] developed an algorithm to estimate PDFs and CDFs using the indicator function approach, while Bierig et al. [26] approximated PDFs via a truncated moment sequence and the method of maximum entropy. Both approaches were used to approximate CDFs of molecular species in biochemical reaction networks [27]. Elfverson et al. [28] used MLMC to estimate failure probabilities, which are singlepoint evaluations of the CDF. Lu et al. [29] obtained a more efficient MLMC algorithm by calibrating the polynomial smoothing of the indicator function in [25] to optimize the smoothing bandwidth for a given value of the error tolerance, thereby achieving a faster variance decay with level, and by enabling an a posteriori switch to MC if the latter turned out to have a lower computational cost. Krumscheid et al. [30]
developed an algorithm for approximating general parametric expectations, including CDFs and characteristic functions, and simultaneously deriving robustness indicators such as quantiles and conditional valuesatrisk.
To reduce the computational cost of MLMC further, the multigrid approach can be combined with a sampling strategy at each level that is more efficient than standard MC. For example, quasiMonte Carlo uses quasirandom, rather than random or pseudorandom, sequences to achieve faster convergence than MC, and has been used to speed up the MLMC computation of the mean system state [31]. Furthermore, a number of socalled “variance reduction” techniques have been developed to obtain estimators with a lower variance than MC for the same number of realizations. These include stratification, antithetic sampling, importance sampling, and control variates [11]. Recently, importance sampling was incorporated into MLMC for a more efficient computation of expectation values of QoIs [32]. Ullman et al. [33]
estimated failure probabilities for rare events by combining subset simulation using Markov chain Monte Carlo
[34] with multilevel failure domains defined on a hierarchy of discrete spatial grids with decreasing mesh sizes. We propose to divide the domain of a random input parameter or initial condition using stratified Monte Carlo to achieve variance reduction at each discretization level, and to estimate the CDF of an output QoI via the resulting “stratified” Multilevel Monte Carlo approach.Section 2 provides a mathematical formulation of MLMC and sMLMC estimators of CDFs and the cost and error associated with computing them. Section 3 discusses two testbed problems for assessing the performance of MLMC and sMLMC compared to standard MC. Conclusions and future research directions are reserved for Section 4.
2 Multilevel Monte Carlo for CDFs
Usually the QoI cannot be simulated directly. Instead, one discretizes on a spatial grid and introduces a sequence of random variables, where is the number of cells in , which converges to as increases. We assume that converges both in the mean and in the sense of distribution to as , i.e.,
(3) 
as for independent of and . We approximate the statistics of rather than of , i.e., our goal is to estimate the CDF of on some compact interval . At each point , is given by
(4) 
We assume the PDF of to be at least times continuously differentiable on for some and . Then, given the values of at a set of equidistant points with separation distance , we estimate at any point
via a piecewise polynomial interpolation of degree
. The resulting approximation of is given by(5) 
where (
) are, e.g., Lagrange basis polynomials or, in the case of cubic spline interpolation, thirddegree polynomials. The goal is to find an unbiased estimator
of , so that an estimator of is computed as(6) 
The estimator converges to the CDF of the original QoI as (to reduce the discretization error) and its variance decreases (to reduce the sampling error).
2.1 Standard Monte Carlo (MC)
The MC estimator for based on independent samples of is defined by
(7) 
where is the sample of . Let denote the MC estimator of , and be a piecewise polynomial interpolation of given its value at a set of points with . The error between the two, , consists of a sampling part related to estimating and a discretization part related to approximating by . Specifically, the mean squared error (MSE) of is bounded by
(8)  
where denotes the norm and refers to the variance operator. Here and are, respectively, the sampling and discretization error, in the mean squared sense. For the root mean squared error (RMSE) to be lower than a prescribed tolerance , it is sufficient to limit both and to, at most, .
For and all other CDF estimators considered in this work, we assume that the number of interpolation points is large enough for the interpolation error to be negligible and for to represent the error between the estimator and the true CDF of .
2.2 Standard Multilevel Monte Carlo (MLMC) without smoothing
Rather than sampling on a single spatial mesh, one considers a sequence of approximations () of associated with discrete meshes {}. Here the number of cells in mesh and , where is the spatial dimension. The idea behind this approach is to start by performing cheaptocompute samples on a coarse mesh, and then gradually correct the resulting estimate of by sampling on finer grids, where generating a realization is more computationally expensive. One rewrites as a telescopic sum
(9a)  
where () are given by  
(9b) 
The MLMC estimator for is defined as
(10) 
While remains approximately constant with , decreases with , allowing the estimator to have the same overall sampling error as its MC counterpart using a decreasing number of samples as increases.
The MSE of the MLMC estimator for is bounded by
(11)  
where and are, respectively, the sampling and discretization error, in the mean squared sense. From the triangle inequality it follows that
(12) 
Hence, to determine the maximum level of an MLMC simulation with given tolerance , we check if is satisfied for the current level . Once the value of is found, we can compare the performance of MLMC to MC by performing the latter on this finest level, i.e., for in Section 2.1 equal to , reusing the samples already computed with MLMC at this level. Since , this strategy ensures that the bound for in (11) is the same as the bound for in (8). To achieve an RSME error of at most , it is sufficient that and .
2.3 Standard Multilevel Monte Carlo (MLMC) with smoothing
The jump discontinuity in the indicator function may lead to a slow decay of
and make MLMC slower than MC for sufficiently large values of the error tolerance [29]. To accelerate the variance decay and to improve the computational efficiency of MLMC, a sigmoidtype smoothing function can be used to remove the singularity in the indicator function. We consider two different smoothing functions, namely a polynomial proposed by Giles et al. [25] (which we will refer to as ) and the CDF of a Gaussian kernel in the context of kernel density estimation [35] (denoted by ).2.3.1 Smoothing via Giles’ polynomial
Under the assumption, made in Section 2, that the PDF of is at least times continuously differentiable on for some and , we define a smoothing function that satisfies [25]

cost of computing for some constant

is Lipschitz continuous

for and for

for .
The function can be constructed as the uniquely determined polynomial of degree (at most) , such that for , and . To obtain an appropriate smoothing function, is extended with
(13) 
The indicator function at each level () is replaced by , where the bandwidth is a measure of the width over which the discontinuity in is smoothed out. The MLMC estimator with smoothing for is given by
(14) 
where with
(15) 
The MLMC estimator with smoothing for is given by a piecewise polynomial interpolation with degree
(16) 
The MSE of is bounded by
(17)  
Compared to (11), (17) contains an additional term , which is the (mean squared) smoothing error. To achieve an RSME error of at most , it is sufficient that , and .
Finding the optimal value for the bandwidth at each level () is crucial to the performance of MLMC. The error should be as close as possible to in order to maximize the smoothing and, hence, the variance decay with increasing level. A possible strategy consists of the following steps.

Given the error tolerance , at level estimate for each interpolation point in by solving
(18) based on a set of initial samples .

Define the smoothing parameter for level , , as
(19) 
Repeat steps 1 and 2 for each new level.
We follow the numerical algorithm in A to compute and to measure the associated computational cost (see Section 2.5). This algorithm is inspired by the approaches [18, 29]. To compare the performance of the MLMC simulation with that of the corresponding MC simulation at the highest level (which ensures that ), we also compute the number of MC samples of required to satisfy and the resulting computational cost of the MC estimator (see A.1).
2.3.2 Smoothing based on kernel density estimation (KDE)
We propose an alternative way of smoothing the indicator function. It is grounded in Kernel Density Estimation (KDE), a nonparametric estimation of the PDF of a random variable. Let be independent and identically distributed samples drawn from the distribution of a QoI with unknown PDF . Then the KDE of is given by
(20) 
where is a kernel and is the bandwidth. We refer to for as a scaled kernel. For a Gaussian kernel, the KDE of is
(21) 
A corresponding estimate of the CDF is then obtained by considering the CDF of the Gaussian kernel, which yields
(22) 
Here
is the CDF of the standard normal distribution and plays the role of indicator smoothing function. The bandwidth
is the counterpart of defined in Section 2.3.1. Based on (22), at each level () we replace with and define an MLMC estimator with smoothing for similar to the one given by (16) but with the smoothing function . The methodology of Section 2.3.1 is then used to bound its MSE and find the optimal value of the bandwidth at each level (). The algorithm in A.1 is deployed to compute and to measure its computational cost.2.4 Stratified Multilevel Monte Carlo (sMLMC)
In stratified Monte Carlo (sMC), one divides the domain of the uncertain input parameter into mutually exclusive and exhaustive regions () called strata. Let be the probability of being in stratum . The expectation value of is
(23) 
where () is the stratum mean given by
(24) 
Here is the conditional CDF of given that . The sMC estimator for is given by
(25) 
where is the number of independent samples generated in the th stratum for each with , and is the th sample of that has a corresponding input parameter () in stratum . The variance of this estimator is
(26) 
with being the variance of within stratum .
Two common choices for are proportional and optimal allocations [11]. For proportional allocation, and
(27)  
(28) 
which shows that stratification produces an estimator with a lower variance than its MC counterpart. For optimal allocation, with for , yielding
(29) 
It can be shown [11] that this is the smallest variance possible for an sMC estimator, which, given (28), is also smaller than the variance of the corresponding MC estimator.
Replacing with in (25) leads to
(30) 
To combine the benefits of stratification with those of the multigrid approach, we replace MC at each level of the MLMC algorithm with sMC and refer to the resulting algorithm as “stratified Multilevel Monte Carlo” (sMLMC). In analogy to (10), the sMLMC estimator for is defined as
(31)  
where with are defined in (9b). The MSE of the sMLMC estimator for is bounded by
(32)  
where the bound for is the same as the bound for in (8), provided . To reduce the computational cost further, we again smooth the indicator function, yielding an additional error term in (32) similar to in (17). The resulting sMLMC estimator with smoothing for is defined as
(33) 
where and
(34) 
for polynomialbased smoothing, or
(35) 
for KDEbased smoothing. The sMLMC estimator with smoothing for is then given by
(36) 
with given by or . To compute and and measure the associated computational cost, we deploy the algorithm in A.2.
2.5 Relative cost of MLMC and sMLMC versus MC
We estimate the total cost of computing the MLMC estimator without smoothing of as an average over independent realizations of the corresponding algorithm,
(37) 
where is the average computational cost of computing a sample of on level for realization , and denotes the finest discretization level at which the sampling is performed for this realization.
For sMLMC without smoothing,
(38) 
where is the average computational cost for computing a sample of in stratum on level for realization , and refers to the finest discretization level for this realization.
To compare the performance of MLMC and sMLMC with or without smoothing with that of MC, we consider realizations of the MC algorithm and perform the realization on level of the corresponding realization of MLMC without smoothing. This provides a single average cost,
(39) 
against which to compare , , and (with for polynomialbased smoothing or for KDEbased smoothing). Here is the number of samples computed in the realization of the MC algorithm.
3 Numerical results
We consider two testbed problems: linear diffusion with a random diffusion coefficient, and inviscid Burgers’ equation with a random initial condition. In both cases, the random input variable is drawn from a truncated lognormal PDF defined on with mean and variance , i.e.,
(40) 
We assume the PDF of the QoI to be at least 3 times continuously differentiable such that we can use a cubicspline interpolation to approximate its CDF . To compare the performance of MC, MLMC and sMLMC, we measure their average computational cost over independent runs for error tolerances , 0.008 and 0.005. The highest possible maximum level for MLMC/sMLMC we consider is . At each level we use proportional allocation for the sMC algorithm. We choose an identical number of warmup samples () for all independent runs of the MLMC algorithm, but use different values for MLMC with and without smoothing to eliminate as much as possible any oversampling ( is lower when smoothing is applied). For sMLMC, we define the number of warmup samples in each stratum , based on with and the chosen allocation strategy (identical for all independent sMLMC runs). The values of for sMLMC is typically lower than those for MLMC and different between the cases with and without smoothing, again to avoid oversampling.
For each value of , the run involves the following steps.

Perform nonsmoothed MLMC, yielding a maximum level .

Perform MC with , reusing already computed samples from the MLMC run.

Perform smoothed MLMC using Giles’ polynomial with a computed smoothing parameter at each level .

Perform smoothed MLMC using KDE with a computed smoothing parameter at each level .

Perform nonsmoothed sMLMC.

Perform smoothed sMLMC using KDE.
3.1 Linear diffusion equation
Consider
(41a)  
(41b)  
(41c) 
Here
Comments
There are no comments yet.