Estimation of distributions via multilevel Monte Carlo with stratified sampling

TL;DR

This paper introduces a stratified multilevel Monte Carlo algorithm that improves efficiency in estimating distribution functions by reducing variance through stratification and smoothing techniques.

Contribution

The paper presents a novel stratified MLMC method combined with kernel density smoothing, enhancing efficiency over standard MLMC in distribution estimation tasks.

Findings

Stratified MLMC reduces variance compared to standard MLMC.

Kernel density smoothing improves estimation efficiency.

The method's computational cost depends on error tolerance.

Abstract

We design and implement a novel algorithm for computing a multilevel Monte Carlo (MLMC) estimator of the cumulative distribution function of a quantity of interest in problems with random input parameters or initial conditions. Our approach combines a standard MLMC method with stratified sampling by replacing standard Monte Carlo at each level with stratified Monte Carlo with proportional allocation. We show that the resulting stratified MLMC algorithm is more efficient than its standard MLMC counterpart, due to the reduction in variance at each level provided by the stratification of the random parameter's domain. A smoothing approximation for the indicator function based on kernel density estimation yields a more efficient algorithm compared to the typically used polynomial smoothing. The difference in computational cost between the smoothing methods depends on the required error tolerance.

Equations

$$F(q)=\mathbb{E}[{\mathcal{I}}_{(-\infty ,q]}(Q)]={\int }_{-\infty }^{\infty }{\mathcal{I}}_{(-\infty ,q]}(s)f(s)\mathrm{d}s={\int }_{-\infty }^{q}f(s)\mathrm{d}s,$$

$${\mathcal{I}}_{(-\infty ,q]}(s)=\{\begin{array}{ll}1& \text{for }s\in (-\infty ,q]\\ 0& \text{for }s\in (q,+\infty ).\end{array}$$

$$\mathbb{E}[Q_M-Q]=\mathcal{O}(M^{-{\alpha }_1}),\mathbb{E}[{\mathcal{I}}_{(-\infty ,q]}(Q_M)-{\mathcal{I}}_{(-\infty ,q]}(Q)]=\mathcal{O}(M^{-{\alpha }_2}),$$

$$F_M(q)=\mathbb{E}[{\mathcal{I}}_{(-\infty ,q]}(Q_M)].$$

$$F_{h,M}(q)=\sum _{n=0}^S\mathbb{E}[{\mathcal{I}}_{(-\infty ,q_n]}(Q_M)]{\varphi }_n(q)\equiv \sum _{n=0}^S\mathbb{E}[{\mathcal{I}}_n(Q_M)]{\varphi }_n(q),$$

$${\hat{F}}_{h,M}(q)=\sum _{n=0}^S{\hat{\mathcal{I}}}_{n,M}{\varphi }_n(q).$$

$${\hat{\mathcal{I}}}_{n,M}^{\text{MC}}=\frac{1}{N_{\text{MC}}}\sum _{j=1}^{N_{\text{MC}}}{\mathcal{I}}_n(Q_M^{(j)}),$$

$$\mathbb{E}[\Vert e({\hat{F}}_{h,M}^{\text{MC}}){\Vert }_{\infty }^2]$$

$$\le \underset{0\le n\le S}{\max }{N}_{\text{MC}}^{-1}\mathbb{V}[{\mathcal{I}}_n(Q_M)]+\underset{0\le n\le S}{\max }|\mathbb{E}[{\mathcal{I}}_n(Q_M)-{\mathcal{I}}_n(Q)]{|}^2,$$

$$\mathbb{E}[{\mathcal{I}}_n(Q_M)]=\mathbb{E}[{\mathcal{I}}_n(Q_0)]+\sum _{l=0}^{L_{\text{max}}}\mathbb{E}[{\mathcal{I}}_n(Q_{M_l})-{\mathcal{I}}_n(Q_{M_{l-1}})]\equiv \sum _{l=0}^{L_{\text{max}}}\mathbb{E}[{\mathcal{I}}_n(Y_l)],$$

$${\mathcal{I}}_n(Y_l)=\{\begin{array}{ll}{\mathcal{I}}_n(Q_{M_l})-{\mathcal{I}}_n(Q_{M_{l-1}})& 1\le l\le L_{\text{max}}\\ {\mathcal{I}}_n(Q_{M_l})& l=0.\end{array}$$

$${\hat{\mathcal{I}}}_{n,M}^{\text{MLMC}}=\sum _{l=0}^{L_{\text{max}}}{\hat{\mathcal{I}}}_n^{\text{MC}}(Y_l)=\sum _{l=0}^{L_{\text{max}}}\frac{1}{N_l}\sum _{j=1}^{N_l}{\mathcal{I}}_n(Y_l^{(j)}).$$

$$\mathbb{E}[\Vert e({\hat{F}}_{h,M}^{\text{MLMC}}){\Vert }_{\infty }^2]$$

$$\le \underset{0\le n\le S}{\max }\sum _{l=0}^{L_{\text{max}}}{N}_l^{-1}\mathbb{V}[{\mathcal{I}}_n(Y_l)]+\underset{0\le n\le S}{\max }|\mathbb{E}[{\mathcal{I}}_n(Q_{M_{L_{\text{max}}}})-{\mathcal{I}}_n(Q)]{|}^2,$$

$$\underset{0\le n\le S}{\max }|\mathbb{E}[{\mathcal{I}}_n(Y_{L_{\text{max}}})|\approx \underset{0\le n\le S}{\max }|\mathbb{E}[{\mathcal{I}}_n(Q_{M_{L_{\text{max}}}})-{\mathcal{I}}_n(Q)]|.$$

$$g_{\text{G}}(s)=\{\begin{array}{ll}1& s<-1\\ 0& s>1.\end{array}$$

$${\hat{g}}_{n,M}^{\text{MLMC}}=\sum _{l=0}^{L_{\text{max}}}{\hat{g}}_n^{\text{MC}}(Y_l),$$

$$g_n(Y_l^{(j)})=\{\begin{array}{ll}{g}_{\text{G}}\left(\frac{Q_{M_l}^{(j)}-q_n}{{\delta }_{\text{G},l}}\right)-g_{\text{G}}\left(\frac{Q_{M_{l-1}}^{(j)}-q_n}{{\delta }_{\text{G},l}}\right)& 1\le l\le L_{\text{max}}\\ g_{\text{G}}\left(\frac{Q_{M_l}^{(j)}-q_n}{{\delta }_{\text{G},l}}\right)& l=0.\end{array}$$

$${\hat{F}}_{h,{\delta }_{\text{G}},M}^{\text{MLMC}}(q)=\sum _{n=0}^S{\hat{g}}_{n,M}^{\text{MLMC}}{\varphi }_n(q)=\sum _{n=0}^S\sum _{l=0}^{L_{\text{max}}}{\hat{g}}_n^{\text{MC}}(Y_l){\varphi }_n(q).$$

$$\mathbb{E}[\Vert e({\hat{F}}_{h,{\delta }_{\text{G}},M}^{\text{MLMC}}){\Vert }_{\infty }^2]\le$$

$$+\underset{e_3^{\text{ML}}}{\underbrace{\Vert \mathbb{E}[{\hat{F}}_{h,{\delta }_{\text{G}},M}^{\text{MLMC}}]-F_{h,M}{\Vert }_{\infty }^2}}.$$

$$\frac{1}{N_l^0}\left|\sum _{j=1}^{N_l^0}\left[g_{\text{G}}\left(\frac{Q_{M_l}^{(j)}-q_n}{{\delta }_{\text{G},l,n}}\right)-{\mathcal{I}}_n(Q_{M_l}^{(j)})\right]\right|=\frac{ϵ}{2}$$

$${\delta }_{\text{G},l}=\underset{0\le n\le S}{\max }{\delta }_{\text{G},l,n}.$$

$${\hat{f}}_{{\delta }_{\text{K}}}(q)=\frac{1}{n}\sum _{i=1}^n{\mathcal{K}}_{{\delta }_{\text{K}}}(q-q_i)=\frac{1}{n{\delta }_{\text{K}}}\sum _{i=1}^n\mathcal{K}\left(\frac{q-q_i}{{\delta }_{\text{K}}}\right)$$

$${\hat{f}}_{{\delta }_{\text{K}}}(q)=\frac{1}{n{\delta }_{\text{K}}}\sum _{i=1}^n\frac{1}{\sqrt{2\pi }}\exp \left[-\frac{1}{2}{\left(\frac{q-q_i}{{\delta }_{\text{K}}}\right)}^2\right].$$

$${\hat{F}}_{{\delta }_{\text{K}}}(q)=\frac{1}{n}\sum _{i=1}^n\Phi \left(\frac{q-q_i}{{\delta }_{\text{K}}}\right).$$

$$\mathbb{E}[Q_M(W)]=\sum _{i=1}^r{p}_i{\zeta }_i,$$

$${\zeta }_i={\int }_{{\mathcal{D}}_i}{Q}_M(w)\mathrm{d}{F}_W^{(i)}(w)=\frac{1}{p_i}{\int }_{{\mathcal{D}}_i}{Q}_M(w)\mathrm{d}{F}_W(w).$$

$${\hat{Q}}_M^{\text{sMC}}=\sum _{i=1}^r\frac{p_i}{n_i}\sum _{j=1}^{n_i}{Q}_M^{(j,i)},$$

$$\mathbb{V}[{\hat{Q}}_M^{\text{sMC}}]=\sum _{i=1}^r\frac{{\sigma }_i^2{p}_i^2}{n_i},{\sigma }_i^2={\int }_{{\mathcal{D}}_i}(Q_M(w)-{\zeta }_i{)}^2\mathrm{d}{F}_W^{(i)}(w)$$

$${\hat{Q}}_M^{\text{sMC, p}}$$

$$\mathbb{V}[{\hat{Q}}_M^{\text{sMC, p}}]$$

$$\mathbb{V}[{\hat{Q}}_M^{\text{sMC, o}}]$$

Full text

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 Karhunen-Loève expansion. This finite set of random variables defines a finite-dimensional 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 $\sqrt{N}$ where $N$ 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 solver-based MLMC as a opposed to traditional grid-based 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 $Q$ to be a function of a continuous random input variable $W$, i.e., $Q=Q(W)$, where $W:\Omega\rightarrow\mathbb{R}$ is a measurable function with $\Omega$ the sample space, the CDF $F$ of $Q$ at a point $q\in\mathbb{R}$ is given by

[TABLE]

where the indicator function $\mathcal{I}_{(-\infty,q]}(s)$ is defined as

[TABLE]

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 single-point 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 values-at-risk.

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, quasi-Monte Carlo uses quasi-random, rather than random or pseudo-random, 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 so-called “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 $Q$ cannot be simulated directly. Instead, one discretizes $Q$ on a spatial grid $\mathcal{T}_{M}$ and introduces a sequence of random variables, $Q_{M}$ where $M$ is the number of cells in $\mathcal{T}_{M}$, which converges to $Q$ as $M$ increases. We assume that $Q_{M}$ converges both in the mean and in the sense of distribution to $Q$ as $M\rightarrow\infty$, i.e.,

[TABLE]

as $M\rightarrow\infty$ for $\alpha_{1},\alpha_{2}\in\mathbb{R}$ independent of $M$ and $q$. We approximate the statistics of $Q_{M}$ rather than of $Q$, i.e., our goal is to estimate the CDF $F_{M}$ of $Q_{M}$ on some compact interval $[a,b]\subset\mathbb{R}$. At each point $q\in[a,b]$, $F_{M}(q)$ is given by

[TABLE]

We assume the PDF $f_{M}$ of $Q_{M}$ to be at least $d$-times continuously differentiable on $[a-\xi_{0},b+\xi_{0}]$ for some $d\in\mathbb{N}_{0}$ and $\xi_{0}>0$. Then, given the values of $F_{M}$ at a set of $S+1$ equidistant points $\mathcal{S}_{h}=\{a=q_{0}<q_{1}<\dots<q_{S}=b\}$ with separation distance $h$, we estimate $F_{M}(q)$ at any point $q\in[a,b]$ via a piecewise polynomial interpolation of degree $\max(d,1)$. The resulting approximation $F_{h,M}(q)$ of $F_{M}(q)$ is given by

[TABLE]

where $\phi_{n}$ ($n=0,\dots,S$) are, e.g., Lagrange basis polynomials or, in the case of cubic spline interpolation, third-degree polynomials. The goal is to find an unbiased estimator $\hat{\mathcal{I}}_{n,M}$ of $\mathbb{E}[\mathcal{I}_{n}(Q_{M})]$, so that an estimator $\hat{F}_{h,M}(q)$ of $F_{h,M}(q)$ is computed as

[TABLE]

The estimator $\hat{F}_{h,M}$ converges to the CDF $F$ of the original QoI $Q$ as $M\rightarrow\infty$ (to reduce the discretization error) and its variance decreases (to reduce the sampling error).

2.1 Standard Monte Carlo (MC)

The MC estimator for $\mathbb{E}[\mathcal{I}_{n}(Q_{M})]$ based on $N_{\text{MC}}$ independent samples of $Q_{M}$ is defined by

[TABLE]

where $Q_{M}^{(j)}$ is the $j^{\text{th}}$ sample of $Q_{M}$. Let $\hat{F}_{h,M}^{\text{MC}}$ denote the MC estimator of $F_{h,M}$, and $F_{h}(q)=\sum_{n=0}^{S}\mathbb{E}[\mathcal{I}_{n}(Q)]\ \phi_{n}$ be a piecewise polynomial interpolation of $F(q)$ given its value at a set of points $\{q_{n}\}=\{a=q_{0}<q_{1}<\dots<q_{S}=b\}$ with $n=0,\dots,S$. The error between the two, $e(\hat{F}_{h,M}^{\text{MC}})\equiv\hat{F}_{h,M}^{\text{MC}}-F_{h}$, consists of a sampling part related to estimating $\mathbb{E}[\mathcal{I}_{n}(Q)]$ and a discretization part related to approximating $Q$ by $Q_{M}$. Specifically, the mean squared error (MSE) of $\hat{F}_{h,M}^{\text{MC}}$ is bounded by

[TABLE]

where $\|\cdot\|_{\infty}$ denotes the $L^{\infty}$ norm and $\mathbb{V}[\cdot]$ refers to the variance operator. Here $e_{1}^{\text{MC}}$ and $e_{2}^{\text{MC}}$ 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 $\epsilon$, it is sufficient to limit both $e_{1}^{\text{MC}}$ and $e_{2}^{\text{MC}}$ to, at most, $\epsilon^{2}/2$.

For $\hat{F}_{h,M}^{\text{MC}}$ and all other CDF estimators considered in this work, we assume that the number of interpolation points $S+1$ is large enough for the interpolation error to be negligible and for $e$ to represent the error between the estimator and the true CDF $F$ of $Q$.

2.2 Standard Multilevel Monte Carlo (MLMC) without smoothing

Rather than sampling $Q_{M}$ on a single spatial mesh, one considers a sequence of approximations $Q_{M_{l}}$ ($l=0,\dots,L_{\text{max}}$) of $Q$ associated with discrete meshes {$\mathcal{T}_{M_{l}},l=0,\dots,L_{\text{max}}$}. Here $M_{l}$ the number of cells in mesh $\mathcal{T}_{M_{l}}$ and $M_{l-1}=2^{-d}M_{l}$, where $d$ is the spatial dimension. The idea behind this approach is to start by performing cheap-to-compute samples on a coarse mesh, and then gradually correct the resulting estimate of $F_{h,M}$ by sampling on finer grids, where generating a realization is more computationally expensive. One rewrites $\mathbb{E}[\mathcal{I}_{n}(Q_{M})]$ as a telescopic sum

[TABLE]

The MLMC estimator for $\mathbb{E}[\mathcal{I}_{n}(Q_{M})]$ is defined as

[TABLE]

While $\mathbb{V}[\mathcal{I}_{n}(Q_{M_{l}})]$ remains approximately constant with $l$, $\mathbb{V}[\mathcal{I}_{n}(Y_{l})]$ decreases with $l$, allowing the estimator $\hat{\mathcal{I}}_{n,M}^{\text{MLMC}}$ to have the same overall sampling error as its MC counterpart $\hat{\mathcal{I}}_{n,M}^{\text{MC}}$ using a decreasing number of samples $N_{l}$ as $l$ increases.

The MSE of the MLMC estimator $\hat{F}_{h,M}^{\text{MLMC}}$ for $F_{h,M}$ is bounded by

[TABLE]

where $e_{1}^{\text{ML}}$ and $e_{2}^{\text{ML}}$ are, respectively, the sampling and discretization error, in the mean squared sense. From the triangle inequality it follows that

[TABLE]

Hence, to determine the maximum level $L_{\text{max}}$ of an MLMC simulation with given tolerance $\epsilon$, we check if $\max_{0\leq n\leq S}|\mathbb{E}[\mathcal{I}_{n}(Y_{L})|\leq\epsilon/\sqrt{2}$ is satisfied for the current level $L$. Once the value of $L_{\text{max}}$ is found, we can compare the performance of MLMC to MC by performing the latter on this finest level, i.e., for $M$ in Section 2.1 equal to $M_{L_{\text{max}}}$, re-using the samples already computed with MLMC at this level. Since $M=M_{L_{\text{max}}}$, this strategy ensures that the bound for $e_{2}^{\text{ML}}$ in (11) is the same as the bound for $e_{2}^{\text{MC}}$ in (8). To achieve an RSME error of at most $\epsilon$, it is sufficient that $e_{1}^{\text{ML}}\leq\epsilon^{2}/2$ and $e_{2}^{\text{ML}}\leq\epsilon^{2}/2$.

2.3 Standard Multilevel Monte Carlo (MLMC) with smoothing

The jump discontinuity in the indicator function may lead to a slow decay of $\mathbb{V}[\mathcal{I}_{n}(Y_{l})]$ and make MLMC slower than MC for sufficiently large values of the error tolerance $\epsilon$ [29]. To accelerate the variance decay and to improve the computational efficiency of MLMC, a sigmoid-type smoothing function $g$ 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 $g_{\text{G}}$) and the CDF of a Gaussian kernel in the context of kernel density estimation [35] (denoted by $g_{\text{K}}$).

2.3.1 Smoothing via Giles’ polynomial

Under the assumption, made in Section 2, that the PDF $f_{M}$ of $Q_{M}$ is at least $d$-times continuously differentiable on $[a-\xi_{0},b+\xi_{0}]$ for some $d\in\mathbb{N}_{0}$ and $\xi_{0}>0$, we define a smoothing function $g_{\text{G}}:\mathbb{R}\rightarrow\mathbb{R}$ that satisfies [25]

cost of computing $g_{\text{G}}(s)\leq C$ $\forall s\in\mathbb{R}$ for some constant $C$ 2. 2.

$g_{\text{G}}$ is Lipschitz continuous 3. 3.

$g_{\text{G}}(s)=1$ for $s<-1$ and $g_{\text{G}}(s)=0$ for $s>1$ 4. 4.

$\int_{-1}^{1}s^{k}\mathcal{I}_{]-\infty,0]}(s-g_{\text{G}}(s))\text{d}s=0$ for $k=0,\dots,d-1$.

The function $g_{\text{G}}$ can be constructed as the uniquely determined polynomial of degree (at most) $d+1$, such that $\int_{-1}^{1}s^{k}g_{\text{G}}(s)\text{d}s=(-1)^{k}/(k+1)$ for $k=0,\dots,d-1$, $g_{\text{G}}(1)=0$ and $g_{\text{G}}(-1)=1$. To obtain an appropriate smoothing function, $g_{\text{G}}$ is extended with

[TABLE]

The indicator function $\mathcal{I}_{n}(Q_{M_{l}})$ at each level $l$ ($l=0,\dots,L_{\text{max}}$) is replaced by $g_{\text{G}}((Q_{M_{l}}-q_{n})/\delta_{\text{G},l})$, where the bandwidth $\delta_{\text{G},l}$ is a measure of the width over which the discontinuity in $\mathcal{I}_{n}(Q_{M_{l}})$ is smoothed out. The MLMC estimator with smoothing for $\mathbb{E}[\mathcal{I}_{n}(Q_{M})]$ is given by

[TABLE]

where $\hat{g}_{n}^{\text{MC}}(Y_{l})=N_{l}^{-1}\sum_{j=1}^{N_{l}}g_{n}(Y_{l}^{(j)})$ with

[TABLE]

The MLMC estimator with smoothing $\hat{F}_{h,\delta_{\text{G}},M}^{\text{MLMC}}$ for $F_{h,M}$ is given by a piecewise polynomial interpolation with degree $\max(d,1)$

[TABLE]

The MSE of $\hat{F}_{h,\delta_{\text{G}},M}^{\text{MLMC}}$ is bounded by

[TABLE]

Compared to (11), (17) contains an additional term $e_{3}^{\text{ML}}$, which is the (mean squared) smoothing error. To achieve an RSME error of at most $\epsilon$, it is sufficient that $e_{1}^{\text{ML}}\leq\epsilon^{2}/4$, $e_{2}^{\text{ML}}\leq\epsilon^{2}/2$ and $e_{3}^{\text{ML}}\leq\epsilon^{2}/4$.

Finding the optimal value for the bandwidth $\delta_{\text{G},l}$ at each level $l$ ($l=0,\dots,L_{\text{max}}$) is crucial to the performance of MLMC. The error $e_{3}^{\text{ML}}$ should be as close as possible to $\epsilon^{2}/4$ 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 $\epsilon$, at level $l=0$ estimate $\delta_{\text{G},l,n}$ for each interpolation point $q_{n}$ in $S_{h}=\{q_{n},n=0,\dots,S\}$ by solving

[TABLE]

based on a set of initial samples $\{Q_{M_{l}}^{(j)}\}_{j=1}^{N_{l}^{0}}$. 2. 2.

Define the smoothing parameter for level $l=0$, $\delta_{\text{G},l}$, as

[TABLE] 3. 3.

Repeat steps 1 and 2 for each new level.

We follow the numerical algorithm in A to compute $\hat{F}_{h,\delta_{\text{G}},M}^{\text{MLMC}}$ 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 $L_{\text{max}}$ (which ensures that $e_{2}^{\text{MC}}\leq\epsilon^{2}/2$), we also compute the number of MC samples of $Q_{L_{\text{max}}}$ required to satisfy $e_{1}^{\text{MC}}\leq\epsilon^{2}/2$ and the resulting computational cost of the MC estimator $\hat{F}_{h,M}^{\text{MC}}$ (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 $q_{1},\dots,q_{n}$ be independent and identically distributed samples drawn from the distribution of a QoI $Q$ with unknown PDF $f$. Then the KDE of $f$ is given by

[TABLE]

where $\mathcal{K}$ is a kernel and $\delta_{\text{K}}$ is the bandwidth. We refer to $K_{\delta_{\text{K}}}(z)=(1/\delta_{\text{K}})\mathcal{K}(z/\delta_{\text{K}})$ for $z\in\mathbb{R}$ as a scaled kernel. For a Gaussian kernel, the KDE of $f$ is

[TABLE]

A corresponding estimate of the CDF $F$ is then obtained by considering the CDF of the Gaussian kernel, which yields

[TABLE]

Here $\Phi$ is the CDF of the standard normal distribution and plays the role of indicator smoothing function. The bandwidth $\delta_{\text{K}}$ is the counterpart of $\delta_{\text{G}}$ defined in Section 2.3.1. Based on (22), at each level $l$ ($l=0,\dots,L_{\text{max}}$) we replace $\mathcal{I}_{n}(Q_{M_{l}})$ with $\Phi[(q_{n}-Q_{M_{l}})/\delta_{\text{K},l}]\equiv g_{\text{K}}[(q_{n}-Q_{M_{l}})/\delta_{\text{K},l}]$ and define an MLMC estimator with smoothing for $F_{h,M}$ similar to the one given by (16) but with the smoothing function $g_{\text{K}}$. The methodology of Section 2.3.1 is then used to bound its MSE and find the optimal value of the bandwidth $\delta_{\text{K},l}$ at each level $l$ ($l=0,\dots,L_{\text{max}}$). The algorithm in A.1 is deployed to compute $\hat{F}_{h,\delta_{\text{K}},M}^{\text{MLMC}}$ and to measure its computational cost.

2.4 Stratified Multilevel Monte Carlo (sMLMC)

In stratified Monte Carlo (sMC), one divides the domain $\mathcal{D}$ of the uncertain input parameter $W$ into $r$ mutually exclusive and exhaustive regions $\mathcal{D}_{i}$ ($i=1,\dots,r$) called strata. Let $p_{i}=\mathbb{P}(W\in\mathcal{D}_{i})$ be the probability of $W$ being in stratum $\mathcal{D}_{i}$. The expectation value of $Q_{M}(W)$ is

[TABLE]

where $\zeta_{i}$ ($i=1,\dots,r$) is the $i^{\text{th}}$ stratum mean given by

[TABLE]

Here $F_{W}^{(i)}$ is the conditional CDF of $W$ given that $W\in\mathcal{D}_{i}$. The sMC estimator for $\mathbb{E}[Q_{M}(W)]$ is given by

[TABLE]

where $n_{i}$ is the number of independent samples generated in the $i$th stratum for each $i=1,\dots,r$ with $\sum_{i=1}^{r}n_{i}\equiv N$, and $Q_{M}^{(j,i)}$ is the $j$th sample of $Q_{M}$ that has a corresponding input parameter ($W$) in stratum $i$. The variance of this estimator is

[TABLE]

with $\sigma_{i}^{2}$ being the variance of $Q_{M}$ within stratum $i$.

Two common choices for $n_{i}$ are proportional and optimal allocations [11]. For proportional allocation, $n_{i}=Np_{i}$ and

[TABLE]

which shows that stratification produces an estimator with a lower variance than its MC counterpart. For optimal allocation, $n_{i}=\alpha_{i}N$ with $\alpha_{i}=\sigma_{i}p_{i}/(\sum_{k=1}^{r}\sigma_{k}p_{k})$ for $i=1,\dots,r$, yielding

[TABLE]

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 $Q_{M}$ with $\mathcal{I}_{n}(Q_{M})$ in (25) leads to

[TABLE]

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 $\mathbb{E}[\mathcal{I}_{n}(Q_{M})]$ is defined as

[TABLE]

where $\mathcal{I}_{n}(Y_{l})$ with $l=0,\dots,L_{\text{max}}$ are defined in (9b). The MSE of the sMLMC estimator $\hat{F}_{h,M}^{\text{sMLMC}}$ for $F_{h,M}$ is bounded by

[TABLE]

where the bound for $e_{2}^{\text{sML}}$ is the same as the bound for $e_{2}^{\text{MC}}$ in (8), provided $M_{L_{\text{max}}}=M$. To reduce the computational cost further, we again smooth the indicator function, yielding an additional error term in (32) similar to $e_{3}^{\text{ML}}$ in (17). The resulting sMLMC estimator with smoothing for $\mathbb{E}[\mathcal{I}_{n}(Q_{M})]$ is defined as

[TABLE]

where $\hat{g}_{n}^{\text{sMC}}(Y_{l})=\sum_{i=1}^{r}(p_{i}/n_{i,l})\sum_{j=1}^{n_{i,l}}g_{n}(Y_{l}^{(j,i)})$ and

[TABLE]

for polynomial-based smoothing, or

[TABLE]

for KDE-based smoothing. The sMLMC estimator with smoothing for $F_{h,M}$ is then given by

[TABLE]

with $\delta$ given by $\delta_{\text{G}}$ or $\delta_{\text{K}}$. To compute $\hat{F}_{h,\delta_{\text{G}},M}^{\text{sMLMC}}$ and $\hat{F}_{h,\delta_{\text{K}},M}^{\text{sMLMC}}$ 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 $F_{h,M}$ as an average over $N_{\text{real}}$ independent realizations of the corresponding algorithm,

[TABLE]

where $\bar{w}_{l}^{(k)}$ is the average computational cost of computing a sample of $Q_{M_{l}}$ on level $l$ for realization $k$, and $L_{\text{max},k}^{\text{MLMC}}$ denotes the finest discretization level at which the sampling is performed for this realization.

For sMLMC without smoothing,

[TABLE]

where $w_{i,l}^{(k)}$ is the average computational cost for computing a sample of $Q_{M_{l}}$ in stratum $i$ on level $l$ for realization $k$, and $L_{\text{max},k}^{\text{sMLMC}}$ 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 $N_{\text{real}}$ realizations of the MC algorithm and perform the $k^{\text{th}}$ realization on level $L_{\text{max},k}^{\text{MLMC}}$ of the corresponding realization of MLMC without smoothing. This provides a single average cost,

[TABLE]

against which to compare $\mathcal{C}(\hat{F}_{h,M}^{\text{MLMC}})$, $\mathcal{C}(\hat{F}_{h,\delta,M}^{\text{MLMC}})$, $\mathcal{C}(\hat{F}_{h,M}^{\text{sMLMC}})$ and $\mathcal{C}(\hat{F}_{h,\delta,M}^{\text{sMLMC}})$ (with $\delta=\delta_{\text{G}}$ for polynomial-based smoothing or $\delta=\delta_{\text{K}}$ for KDE-based smoothing). Here $N_{\text{MC}}^{(k)}$ is the number of samples computed in the $k^{\text{th}}$ 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 $W$ is drawn from a truncated lognormal PDF $f_{W}(w;\mu,\sigma,w_{l},w_{u})$ defined on $[w_{l},w_{u}]$ with mean $\mu$ and variance $\sigma^{2}$, i.e.,

[TABLE]

We assume the PDF $f$ of the QoI $Q(W)$ to be at least 3 times continuously differentiable such that we can use a cubic-spline interpolation to approximate its CDF $F$. To compare the performance of MC, MLMC and sMLMC, we measure their average computational cost over $N_{\text{real}}=50$ independent runs for error tolerances $\epsilon=0.01$, 0.008 and 0.005. The highest possible maximum level $L_{\text{max}}^{\star}$ for MLMC/sMLMC we consider is $L_{\text{max}}^{\star}=7$. At each level we use proportional allocation for the sMC algorithm. We choose an identical number of warmup samples $N_{l}^{0}$ ($l=0,\dots,L_{\text{max}}^{\star}$) 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 ($N_{l}^{0}$ is lower when smoothing is applied). For sMLMC, we define the number of warmup samples $n_{i,l}^{0}$ in each stratum $i=1,\dots,r$, based on $N_{l}^{0}$ with $l=0,\dots,L_{\text{max}}^{\star}$ and the chosen allocation strategy (identical for all independent sMLMC runs). The values of $N_{l}^{0}$ 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 $\epsilon$, the $k^{\text{th}}$ run involves the following steps.

Perform non-smoothed MLMC, yielding a maximum level $L_{\text{max},k}^{\text{MLMC}}$. 2. 2.

Perform MC with $M=M_{L_{\text{max},k}^{\text{MLMC}}}$, re-using already computed samples from the MLMC run. 3. 3.

Perform smoothed MLMC using Giles’ polynomial with a computed smoothing parameter $\delta_{\text{G},l}$ at each level $l$. 4. 4.

Perform smoothed MLMC using KDE with a computed smoothing parameter $\delta_{\text{K},l}$ at each level $l$. 5. 5.

Perform non-smoothed sMLMC. 6. 6.

Perform smoothed sMLMC using KDE.

3.1 Linear diffusion equation

Consider

[TABLE]

Here $D\in\Omega_{D}=[1,4]$ is a lognormal random variable with PDF $f_{D}(w;3,3,1,4)$. Our goal is to compute the CDF of a QoI

[TABLE]

We discretize (41) in space using a central finite difference scheme, and then apply the Crank-Nicholson method to the resulting system of initial value problems. The matrix associated with this linear system is tridiagonal Toeplitz, hence we apply the Thomas algorithm to solve it at reduced computational complexity. This numerical scheme is second-order accurate in both space and time. We approximate the CDF $F$ of Q on the interval [14,28], which contains the support of $F$, and use $S+1=29$ interpolation points, i.e., set $h=0.5$.

We first compare the computational cost of MC, $\mathcal{C}(\hat{F}_{h,M}^{\text{MC}})$, to that of non-smoothed MLMC, $\mathcal{C}(\hat{F}_{h,M}^{\text{MLMC}})$, and MLMC with the polynomial-based, $\mathcal{C}(\hat{F}_{h,\delta_{\text{G}},M}^{\text{MLMC}})$, and KDE-based, $\mathcal{C}(\hat{F}_{h,\delta_{\text{K}},M}^{\text{MLMC}})$, smoothing of the indicator function. Figure 1 (left) shows that even without smoothing the cost of MLMC is about an order of magnitude lower than that of MC, at the lowest error tolerance considered ($\epsilon=0.005$). Applying the smoothing further reduces the computational cost by almost an order of magnitude at this tolerance (Fig. 1, right). We find that KDE-based smoothing yields a lower computational cost than its polynomial-based counterpart, and hence only consider KDE for smoothing the indicator function in the sMLMC algorithm.

Figure 2 (left) illustrates the evolution of $\tilde{V}[\mathcal{I}_{n}(Q_{M_{l}})]$, $\tilde{V}[\mathcal{I}_{n}(Y_{l})]$ and $\tilde{V}[g_{n}(Y_{l})]$ (using KDE-based smoothing) with level $l$ for a single run and $\epsilon=0.005$. Here $\tilde{V}$ denotes a sample estimate of $\mathbb{V}$. While $\tilde{V}[\mathcal{I}_{n}(Q_{M_{l}})]$ remains approximately constant as the spatial resolution increases, $\tilde{V}[\mathcal{I}_{n}(Y_{l})]$ and $\tilde{V}[g_{n}(Y_{l})]$ decay as the spatial mesh is refined, so that fewer samples are needed at higher levels of discretization (Fig. 2, right).

Next, we stratify the sample space of $D$, $\Omega_{D}$, into 8 or 16 strata of equal width and compare the computational cost of the (non-stratified) MLMC algorithm without smoothing, $\mathcal{C}(\hat{F}_{h,M}^{\text{MLMC}})$, to that of sMLMC without smoothing, $\mathcal{C}(\hat{F}_{h,M}^{\text{sMLMC}})$ and with KDE-based smoothing, $\mathcal{C}(\hat{F}_{h,\delta_{\text{K}},M}^{\text{sMLMC}})$, for 8 and 16 strata and the error tolerances considered. Figure 1 (right) demonstrates that stratifying the non-smoothed MLMC algorithm yields similar computational cost savings as applying smoothing to the indicator function. Combining KDE-based smoothing with stratification further reduces the computational cost by almost an order of magnitude for $\epsilon=0.005$, yielding between one and two orders of magnitude cost savings compared to the non-smoothed MLMC method at this tolerance. We note that for 16 strata and $\epsilon=0.01$, smoothing increases, rather than decreases, the overall computational cost of the sMLMC algorithm. In this case, the number of samples on each level is already quite low for non-smoothed sMLMC and smoothing the indicator function is not beneficial. As predicted by (28), Figure 2 (left) demonstrates that stratification yields an additional variance reduction. Here and for the smoothed case, we compare the quantity $N_{l}\cdot\max_{0\leq n\leq S}\tilde{V}[\hat{\mathcal{I}}_{n}^{\text{sMC, p}}(Y_{l})]$ with $\tilde{V}[\mathcal{I}_{n}(Y_{l})]$. Combining stratification with smoothing yields further variance reduction. This is translated into a decrease in the required number of samples (Fig. 2, right).

3.2 Inviscid Burgers’ equation

Consider

[TABLE]

The initial state $U_{1}\in\Omega_{U_{1}}=[0,2]$ is a lognormal random variable with PDF $f_{U_{1}}(w;1.5,1,0,2)$. Our goal is to compute the CDF of a QoI

[TABLE]

We discretize (43) using the Godunov method, which is first-order accurate in both space and time. This is a conservative finite volume scheme which solves a Riemann problem at each inter-cell boundary. We approximate the CDF $F$ of Q over the interval [15,65], which contains the support of $F$, and use $S+1=101$ interpolation points, i.e., set $h=0.5$.

We compare the computational cost of MC, $\mathcal{C}(\hat{F}_{h,M}^{\text{MC}})$, to that of non-smoothed MLMC, $\mathcal{C}(\hat{F}_{h,M}^{\text{MLMC}})$, and MLMC with polynomial-based, $\mathcal{C}(\hat{F}_{h,\delta_{\text{G}},M}^{\text{MLMC}})$, and KDE-based, $\mathcal{C}(\hat{F}_{h,\delta_{\text{K}},M}^{\text{MLMC}})$, smoothing of the indicator function. Figure 3 (left) illustrates that, as in the linear diffusion problem, MLMC is more efficient than MC even without smoothing. However, the computational cost savings are much lower in this case. The KDE-based smoothing is more efficient than the polynomial-based smoothing, especially at the lowest error tolerance $\epsilon=0.005$. However, the gain in computational efficiency from smoothing is relatively modest.

Next, we stratify the sample space of $U_{1}$, $\Omega_{U_{1}}$, into 8 or 16 strata of equal width and compare the computational cost of the (non-stratified) MLMC algorithm without smoothing, $\mathcal{C}(\hat{F}_{h,M}^{\text{MLMC}})$, to that of sMLMC without smoothing, $\mathcal{C}(\hat{F}_{h,M}^{\text{sMLMC}})$, and with KDE-based smoothing, $\mathcal{C}(\hat{F}_{h,\delta_{\text{K}},M}^{\text{sMLMC}})$, for 8 and 16 strata and the error tolerances considered. Figure 3 (right) shows that the computational saving from applying stratification are much more substantial than those originating from smoothing the indicator function, unlike for the linear diffusion problem where the non-smoothed sMLMC and smoothed MLMC algorithm had a similar computational cost. Combining KDE-based smoothing with stratification again yields the most efficient result, saving about an order of magnitude in computational cost compared to standard MLMC with or without smoothing.

4 Conclusions

We constructed a stratified Multilevel Monte Carlo (sMLMC) algorithm for estimating the cumulative distribution function (CDF) of a quantify of interest (QoI) in a problem with random input parameters or a random initial state. Our method combines the benefits of multigrid from standard MLMC with variance reduction from stratified sampling in each of the discretization levels. We also explored the use of Gaussian kernel density estimator (KDE) in lieu of the currently used polynomial-based smothers.

Our study yields the following major conclusions:

For all test problems and error tolerances considered, the computational cost of non-smoothed sMLMC is smaller than that of non-smoothed MLMC. This can be attributed to the fact that, as spatial resolution increases, the variance decays faster for sMLMC than for MLMC due to the additional variance reduction from the stratification at each level. 2. 2.

Stratifying MLMC can yield either similar or substantially higher computational cost savings compared to smoothing the indicator function, depending on the considered problem. 3. 3.

Smoothing the indicator function in the sMLMC algorithm yields noticeable additional computational cost savings, with a few exceptions at high tolerances where smoothing either produces a negligible improvement or even a reduction in the algorithm’s efficiency. Like in the non-smoothed case, smoothed sMLMC is more efficient than its smoothed MLMC counterpart. 4. 4.

For the same level of smoothing error, KDE-based smoothing yields a more efficient algorithm than its polynomial-based counterpart. The gain in efficiency depends on the required error tolerance and the problem considered.

5 Acknowledgments

This work was supported in part by Defense Advanced Research Project Agency under award number 101513612, by Air Force Office of Scientific Research under award number FA9550-17-1-0417, and by U.S. Department of Energy under award number DE-SC0019130.

Appendix A Numerical algorithms for computing $\hat{F}_{h,\delta,M}^{\text{MLMC}}$ and $\hat{F}_{h,\delta,M}^{\text{sMLMC}}$

A.1 MLMC with smoothing

A.2 sMLMC with smoothing