Multilevel Monte Carlo Algorithm Research Articles

Analysis of a Class of Multilevel Markov Chain Monte Carlo Algorithms Based on Independent Metropolis–Hastings

In this work, we present, analyze, and implement a class of multilevel Markov chain Monte Carlo (ML-MCMC) algorithms based on independent Metropolis–Hastings proposals for Bayesian inverse problems. In this context, the likelihood function involves solving a complex differential model, which is then approximated on a sequence of increasingly accurate discretizations. The key point of this algorithm is to construct highly coupled Markov chains together with the standard multilevel Monte Carlo argument to obtain a better cost-tolerance complexity than a single-level MCMC algorithm. Our method extends the ideas of Dodwell et al., [SIAM/ASA J. Uncertain. Quantif., 3 (2015), pp. 1075–1108] to a wider range of proposal distributions. We present a thorough convergence analysis of the ML-MCMC method proposed, and show, in particular, that (i) under some mild conditions on the (independent) proposals and the family of posteriors, there exists a unique invariant probability measure for the coupled chains generated by our method, and (ii) that such coupled chains are uniformly ergodic. We also generalize the cost-tolerance theorem of Dodwell et al. to our wider class of ML-MCMC algorithms. Finally, we propose a self-tuning continuation-type ML-MCMC algorithm. The presented method is tested on an array of academic examples, where some of our theoretical results are numerically verified. These numerical experiments evidence how our extended ML-MCMC method is robust when targeting some pathological posteriors, for which some of the previously proposed ML-MCMC algorithms fail.

Goal-oriented adaptive finite element multilevel Monte Carlo with convergence rates

In this study, we present an adaptive multilevel Monte Carlo (AMLMC) algorithm for approximating deterministic, real-valued, bounded linear functionals that depend on the solution of a linear elliptic PDE with a lognormal diffusivity coefficient and geometric singularities in bounded domains of Rd. Our AMLMC algorithm is built on the results of the weak convergence rates in the work (Moon et al., 2006) for an adaptive algorithm using isoparametric d-linear quadrilateral finite element approximations and the dual weighted residual error representation in a deterministic setting. Designed to suit the geometric nature of the singularities in the solution, our AMLMC algorithm uses a sequence of deterministic, non-uniform auxiliary meshes as a building block. The above-mentioned deterministic adaptive algorithm generates these meshes, corresponding to a geometrically decreasing sequence of tolerances. In particular, for a given realization of the diffusivity coefficient and accuracy level, AMLMC constructs its approximate sample using the first mesh in the hierarchy that satisfies the corresponding bias accuracy constraint. This adaptive approach is particularly useful for the lognormal case treated here, which lacks uniform coercivity and thus produces functional outputs that vary over orders of magnitude when sampled. Furthermore, we discuss iterative solvers and compare their efficiency with direct ones. To reduce computational work, we propose a stopping criterion for the iterative solver with respect to the quantity of interest, the realization of the diffusivity coefficient, and the desired level of AMLMC approximation.From the numerical experiments, based on a Fourier expansion of the diffusivity coefficient field, we observe improvements in efficiency compared with both standard Monte Carlo (MC) and standard MLMC (SMLMC) for a problem with a singularity similar to that at the tip of a slit modeling a crack.

Multilevel Markov Chain Monte Carlo for Bayesian Inversion of Parabolic Partial Differential Equations under Gaussian Prior

We analyze the convergence of a multilevel Markov chain Monte Carlo (MLMCMC) algorithm for the Bayesian estimation of solution functionals for linear, parabolic partial differential equations subje...

Multilevel Monte Carlo algorithm for quantum mechanics on a lattice

Monte Carlo simulations of quantum field theories on a lattice become increasingly expensive as the continuum limit is approached since the cost per independent sample grows with a high power of the inverse lattice spacing. Simulations on fine lattices suffer from critical slowdown, the rapid growth of autocorrelations in the Markov chain. This causes a strong increase in the number of lattice configurations that have to be generated to obtain statistically significant results. This paper discusses hierarchical sampling methods to tame the growth in autocorrelations. Combined with multilevel variance reduction, this significantly reduces the computational cost of simulations for given tolerances $\epsilon_{\text{disc}}$ on the discretisation error and $\epsilon_{\text{stat}}$ on the statistical error. For observables with lattice errors of order $\alpha$ and integrated autocorrelation times that grow like $\tau_{\mathrm{int}}\propto a^{-z}$, multilevel Monte Carlo (MLMC) reduces the cost from $\mathcal{O}(\epsilon_{\text{stat}}^{-2}\epsilon_{\text{disc}}^{-(1+z)/\alpha})$ to $\mathcal{O}(\epsilon_{\text{stat}}^{-2}\vert\log \epsilon_{\text{disc}} \vert^2+\epsilon_{\text{disc}}^{-1/\alpha})$ or $\mathcal{O}(\epsilon_{\text{stat}}^{-2}+\epsilon_{\text{disc}}^{-1/\alpha})$. Higher gains are expected for simulations of quantum field theories in $D$ dimensions. The efficiency of the approach is demonstrated on two model systems, including a topological oscillator that is badly affected by critical slowdown from topological charge freezing. On fine lattices, the new methods are orders of magnitude faster than standard Hybrid Monte Carlo sampling. For high resolutions, MLMC can be used to accelerate even the cluster algorithm for the topological oscillator. Performance is further improved through perturbative matching which guarantees efficient coupling of theories on the multilevel hierarchy.

Accelerated Multilevel Monte Carlo With Kernel‐Based Smoothing and Latinized Stratification

AbstractHeterogeneity and a paucity of measurements of key material properties undermine the veracity of quantitative predictions of subsurface flow and transport. For such model forecasts to be useful as a management tool, they must be accompanied by computationally expensive uncertainty quantification, which yields confidence intervals, probability of exceedance, and so forth. We design and implement novel multilevel Monte Carlo (MLMC) algorithms that accelerate estimation of the cumulative distribution functions (CDFs) of quantities of interest, for example, water breakthrough time or oil production rate. Compared to standard non‐smoothed MLMC, the new estimators achieve a significant variance reduction at each discretization level by smoothing the indicator function with a Gaussian kernel or replacing standard Monte Carlo (MC) with the recently developed hierarchical Latinized stratified sampling (HLSS). After validating the kernel‐smoothed MLMC and HLSS‐enhanced MLMC methods on a single‐phase flow test bed, we demonstrate that they are orders of magnitude faster than standard MC for estimating the CDF of breakthrough times in multiphase flow problems.

Estimation of distributions via multilevel Monte Carlo with stratified sampling

We design and implement a novel algorithm for computing a multilevel Monte Carlo (MLMC) estimator of the joint cumulative distribution function (CDF) of a vector-valued quantity of interest in problems with random input parameters and initial conditions. Our approach combines MLMC with stratified sampling of the input sample space by replacing standard Monte Carlo at each level with stratified Monte Carlo initialized with proportionally allocated samples. We show that the resulting stratified MLMC (sMLMC) algorithm is more efficient than its standard MLMC counterpart due to the additional variance reduction provided by the stratification of the random parameter's domain, especially at the coarsest levels. Additional computational cost savings are obtained by smoothing the indicator function with a Gaussian kernel, which proves to be an efficient and robust alternative to recently developed polynomial-based techniques.

An adaptive multilevel Monte Carlo algorithm for the stochastic drift–diffusion–Poisson system

We present an adaptive multilevel Monte Carlo algorithm for solving the stochastic drift–diffusion–Poisson system with non-zero recombination rate. The a-posteriori error is estimated to enable goal-oriented adaptive mesh refinement for the spatial dimensions, while the a-priori error is estimated to guarantee linear convergence of the H1 error. In the adaptive mesh refinement, efficient estimation of the error indicator gives rise to better error control. For the stochastic dimensions, we use the multilevel Monte Carlo method to solve this system of stochastic partial differential equations. Finally, the advantage of the technique developed here compared to uniform mesh refinement is discussed using a realistic numerical example.

Quantifying uncertain system outputs via the multilevel Monte Carlo method — Part I: Central moment estimation

In this work we introduce and analyze a novel multilevel Monte Carlo (MLMC) estimator for the accurate approximation of central moments of system outputs affected by uncertainties. Central moments play a central role in many disciplines to characterize a random system output's distribution and are of primary importance in many prediction, optimization, and decision making processes under uncertainties. We detail how to effectively tune the MLMC algorithm for central moments of any order and present a complete practical algorithm that is implemented as part of a Python library [1]. In fact, we validate the methodology on selected reference problems and apply it to an aerodynamic relevant test case, namely the transonic RAE 2822 airfoil affected by operating and geometric uncertainties.

Numerical simulation of the stochastic Burgers’ equation using MLMC and CBC algorithm

Burgers turbulence is a model for developing tools to study the Navier-Stokes turbulence, many hydrodynamical problemsin the theory related to random Lagrangian systems. Many questions that are generally asked in Navier-Stokes turbulence can be answered using Burgers turbulence. The aim of the present paper is to apply Multi-level Monte Carlo (MLMC) on stochastic Burgers’ equation by using component-by-component algorithm (CBC). CBC algorithm is developed by the concept of circulant matrix that reduces the cost as a Quasi Monte Carlo technique from O(s n) to O(s n log n) where s is the dimension of integral with equi-distributed points. In this paper, Burgers’ equation is discretized using the finite-volume technique, the MLMC with different random samples is applied and the stability is tested. The results show that MLMC is suitable only for some cases of stochastic differential equations (SDEs) when using pseudo random generator, which is Monte Carlo with high cost than using CBC.

MULTILEVEL MONTE CARLO ON A HIGH-DIMENSIONAL PARAMETER SPACE FOR TRANSMISSION PROBLEMS WITH GEOMETRIC UNCERTAINTIES

In the framework of uncertainty quantification, we consider a quantity of interest which depends non-smoothly on the high-dimensional parameter representing the uncertainty. We show that, in this situation, the multilevel Monte Carlo algorithm is a valid option to compute moments of the quantity of interest (here we focus on the expectation), as it allows to bypass the precise location of discontinuities in the parameter space. We illustrate how such lack of smoothness occurs for the point evaluation of the solution to a (Helmholtz) transmission problem with uncertain interface, if the point can be crossed by the interface for some realizations. For this case, we provide a space regularity analysis for the solution, in order to state converge results in the L1-norm for the finite element discretization. The latter are then used to determine the optimal distribution of samples among the Monte Carlo levels. Particular emphasis is given on the robustness of our estimates with respect to the dimension of the parameter space.

Multilevel Monte Carlo Method for Stochastic Analysis of Fluid‐Structure Interaction

AbstractAbstract: Uncertainty quantification for multi‐physics problems like fluid‐structure interaction (FSI) is challenging in terms of the required computational effort and complexity of physics involved. A Multilevel Monte Carlo (MLMC) approach for stochastic analysis of FSI problems is presented in the current study. An overall reduction in computational cost is achieved by variance reduction in MLMC method. The efficacy of MLMC algorithm for expensive problems like FSI is compared with classical Monte Carlo(MC) method and MLMC is found to achieve the required accuracy at a significantly reduced computational cost. An FSI benchmark test case with uncertain inputs is presented here. The possibility of uncertainty quantification for FSI systems with a plausible computational cost using MLMC is demonstrated.

An FEM-MLMC algorithm for a moving shutter diffraction in time stochastic model

We consider a moving shutter and non-deterministic generalization of the diffraction in time model introduced by Moshinsky several decades ago to study a class of quantum transients. We first develop a moving-mesh finite element method (FEM) to simulate the determisitic version of the model. We then apply the FEM and multilevel Monte Carlo (MLMC) algorithm to the stochastic moving-domain model for simulation of approximate statistical moments of the density profile of the stochastic transients.

Computational Complexity Analysis for Monte Carlo Approximations of Classically Scaled Population Processes

We analyze and compare the computational complexity of different simulation strategies for Monte Carlo in the setting of classically scaled population processes. This allows a range of widely used competing strategies to be judged systematically. Our setting includes stochastically modeled biochemical systems. We consider the task of approximating the expected value of some path functional of the state of the system at a fixed time point. We study the use of standard Monte Carlo when samples are produced by exact simulation and by approximation with tau-leaping or an Euler-Maruyama discretization of a diffusion approximation. Appropriate modifications of recently proposed multilevel Monte Carlo algorithms are also studied for the tau-leaping and Euler-Maruyama approaches. In order to quantify computational complexity in a tractable yet meaningful manner, we consider a parameterization that, in the mass action chemical kinetics setting, corresponds to the classical system size scaling. We base the analysis on a novel asymptotic regime where the required accuracy is a function of the model scaling parameter. Our new analysis shows that, under the specific assumptions made in the manuscript, if the bias inherent in the diffusion approximation is smaller than the required accuracy, then multilevel Monte Carlo for the diffusion approximation is most efficient, besting multilevel Monte Carlo with tau-leaping by a factor of a logarithm of the scaling parameter. However, if the bias of the diffusion model is greater than the error tolerance, or if the bias can not be bounded analytically, multilevel versions of tau-leaping are often the optimal choice.

A multilevel Monte Carlo finite difference method for random scalar degenerate convection–diffusion equations

This paper proposes a finite difference multilevel Monte Carlo algorithm for degenerate parabolic convection–diffusion equations where the convective and diffusive fluxes are allowed to be random. We establish a notion of stochastic entropy solutions to these equations. Our chief goal is to efficiently compute approximations to statistical moments of these stochastic entropy solutions. To this end, we design a multilevel Monte Carlo method based on a finite volume scheme for each sample. We present a novel convergence rate analysis of the combined multilevel Monte Carlo finite volume method, allowing in particular for low [Formula: see text]-integrability of the random solution with [Formula: see text], and low deterministic convergence rates (here, the theoretical rate is [Formula: see text]). We analyze the design and error versus work of the multilevel estimators. We obtain that the maximal rate (based on optimizing possibly the pessimistic upper bounds on the discretization error) is obtained for [Formula: see text], for finite volume convergence rate of [Formula: see text]. We conclude with numerical experiments.

A Continuation Multi Level Monte Carlo (C-MLMC) method for uncertainty quantification in compressible inviscid aerodynamics

In this work we apply the Continuation Multi-Level Monte Carlo (C-MLMC) algorithm proposed in Collier et al. (2014) to efficiently propagate operating and geometric uncertainties in inviscid compressible aerodynamics numerical simulations. The key idea of MLMC is that one can draw MC samples simultaneously and independently on several approximations of the problem under investigations on a hierarchy of nested computational grids (levels). The expectation of an output quantity is computed as a sample average using the coarsest solutions and corrected by averages of the differences of the solutions of two consecutive grids in the hierarchy. By this way, most of the computational effort is transported from the finest level (as in a standard Monte Carlo approach) to the coarsest one. The continuation algorithm (C-MLMC) is a robust and self-tuning version that estimates on the fly the optimal number of level and realizations per level. In this work we describe in detail how C-MLMC can be adapted to perform uncertainty quantification analysis in compressible aerodynamics and we apply it to the relevant test cases of a quasi 1D convergent–divergent Laval nozzle and the 2D transonic RAE-2822 airfoil.

A CLT for infinitely stratified estimators, with applications to debiased MLMC

This paper develops a general central limit theorem (CLT) for post-stratified Monte Carlo estimators with an associated infinite number of strata. In addition, consistency of the corresponding variance estimator is established in the same setting. With these results in hand, one can then construct asymptotically valid confidence interval procedures for such infinitely stratified estimators. We then illustrate our general theory, by applying it to the specific case of debiased multi-level Monte Carlo (MLMC) algorithms. This leads to the first asymptotically valid confidence interval procedure for such stratified debiased MLMC procedures.

Numerical approximation of statistical solutions of planar, incompressible flows

We present a finite difference-(multi-level) Monte Carlo algorithm to efficiently compute statistical solutions of the two-dimensional incompressible Navier–Stokes equations (NSE), with periodic boundary conditions and for arbitrary high Reynolds number. We propose a reformulation of statistical solutions in the sense of Foiaş and Prodi in the vorticity-stream function form. The vorticity-stream function formulation of the NSE in two-space dimensions is discretized with a finite difference scheme. We obtain a convergence rate error estimate for this approximation which is explicit in the viscosity parameter [Formula: see text], under realistic assumptions on the solution regularity. We also prove convergence and complexity estimates, for the (multi-level) Monte Carlo finite difference algorithm to compute statistical solutions. Numerical experiments illustrating the validity of our estimates are presented. They show that the multi-level Monte Carlo algorithm can significantly accelerate the computation of statistical solutions in the sense of Foiaş and Prodi, even for very high Reynolds numbers.

Approximate Bayesian Computation by Subset Simulation using hierarchical state-space models

A new multi-level Markov Chain Monte Carlo algorithm for Approximate Bayesian Computation, ABC-SubSim, has recently appeared that exploits the Subset Simulation method for efficient rare-event simulation. ABC-SubSim adaptively creates a nested decreasing sequence of data-approximating regions in the output space that correspond to increasingly closer approximations of the observed output vector in this output space. At each level, multiple samples of the model parameter vector are generated by a component-wise Metropolis algorithm so that the predicted output corresponding to each parameter value falls in the current data-approximating region. Theoretically, if continued to the limit, the sequence of data-approximating regions would converge on to the observed output vector and the approximate posterior distributions, which are conditional on the data-approximation region, would become exact, but this is not practically feasible. In this paper we study the performance of the ABC-SubSim algorithm for Bayesian updating of the parameters of dynamical systems using a general hierarchical state-space model. We note that the ABC methodology gives an approximate posterior distribution that actually corresponds to an exact posterior where a uniformly distributed combined measurement and modeling error is added. We also note that ABC algorithms have a problem with learning the uncertain error variances in a stochastic state-space model and so we treat them as nuisance parameters and analytically integrate them out of the posterior distribution. In addition, the statistical efficiency of the original ABC-SubSim algorithm is improved by developing a novel strategy to regulate the proposal variance for the component-wise Metropolis algorithm at each level. We demonstrate that Self-regulated ABC-SubSim is well suited for Bayesian system identification by first applying it successfully to model updating of a two degree-of-freedom linear structure for three cases: globally, locally and un-identifiable model classes, and then to model updating of a two degree-of-freedom nonlinear structure with Duffing nonlinearities in its interstory force-deflection relationship.

A Multilevel Adaptive Reaction-splitting Simulation Method for Stochastic Reaction Networks

In this work, we present a novel multilevel Monte Carlo method for kinetic simulation of stochastic reaction networks characterized by having simultaneously fast and slow reaction channels. To produce efficient simulations, our method adaptively classifies the reactions channels into fast and slow channels. To this end, we first introduce a state-dependent quantity named level of activity of a reaction channel. Then, we propose a low-cost heuristic that allows us to adaptively split the set of reaction channels into two subsets characterized by either a high or a low level of activity. Based on a time-splitting technique, the increments associated with high-activity channels are simulated using the tau-leap method, while those associated with low-activity channels are simulated using an exact method. This path simulation technique is amenable for coupled path generation and a corresponding multilevel Monte Carlo algorithm. To estimate expected values of observables of the system at a prescribed final time...

Mixed finite element analysis of lognormal diffusion and multilevel Monte Carlo methods

This work is motivated by the need to develop efficient tools for uncertainty quantification in subsurface flows associated with radioactive waste disposal studies. We consider single phase flow problems in random porous media described by correlated lognormal distributions. We are interested in the error introduced by a finite element discretisation of these problems. In contrast to several recent works on the analysis of standard nodal finite element discretisations, we consider here mass-conservative lowest order Raviart–Thomas mixed finite elements. This is very important since local mass conservation is highly desirable in realistic groundwater flow problems. Due to the limited spatial regularity and the lack of uniform ellipticity and boundedness of the operator the analysis is non-trivial in the presence of lognormal random fields. We establish finite element error bounds for Darcy velocity and pressure, as well as for a more accurate recovered pressure approximation. We then apply the error bounds to prove convergence of the multilevel Monte Carlo algorithm for estimating statistics of these quantities. Moreover, we prove convergence for a class of bounded, linear functionals of the Darcy velocity. An important special case is the approximation of the effective permeability in a 2D flow cell. We perform numerical experiments to confirm the convergence results.