Source-linked AI summary

A Hierarchical Multilevel Markov Chain Monte Carlo Algorithm with Applications to Uncertainty Quantification in Subsurface Flow

T. J. Dodwell, C. Ketelsen, R. Scheichl, A. L. Teckentrup

arXiv:1303.7343v2math.NA

TL;DR

The paper addresses the high cost of Metropolis-Hastings inference for large-scale, high-dimensional PDE models. It introduces a multilevel Metropolis-Hastings estimator with a problem-dependent complexity theorem and evaluates it for subsurface flow. The multilevel method achieves one lower power of ε in cost asymptotically and substantially lower observed costs near ε = 10^−3.

  • Problem

    Large-scale, high-dimensional Bayesian inverse problems are costly because likelihood evaluations are expensive and standard MCMC may require many samples with low acceptance rates.

  • Method

    The paper uses a hierarchy of computational models, level-dependent posteriors, parallel chains, and a new fine-level proposal to estimate differences between successive models.

  • Results

    O(ε^−(d+1)−δ) versus O(ε^−(d+2)−δ) gives the multilevel estimator one lower power of ε in cost; near ε = 10^−3, observed costs are about O(10–50) times lower.

  • Takeaways & Limitations

    The method provides substantial computational-cost reductions for subsurface-flow inference and is presented as generally applicable across computational models.

  • Takeaways & Limitations

    The analysis does not prove the MCMC mean-square-error assumption or fully address the per-sample likelihood-cost assumption; the latter holds with γ = 1 in the best case of an optimal linear solver.

Abstract

from arXiv · show

In this paper we address the problem of the prohibitively large computational cost of existing Markov chain Monte Carlo methods for large--scale applications with high dimensional parameter spaces, e.g. in uncertainty quantification in porous media flow. We propose a new multilevel Metropolis-Hastings algorithm, and give an abstract, problem dependent theorem on the cost of the new multilevel estimator based on a set of simple, verifiable assumptions. For a typical model problem in subsurface flow, we then provide a detailed analysis of these assumptions and show significant gains over the standard Metropolis-Hastings estimator. Numerical experiments confirm the analysis and demonstrate the effectiveness of the method with consistent reductions of more than an order of magnitude in the cost of the multilevel estimator over the standard Metropolis-Hastings algorithm for tolerances $\varepsilon < 10^{-2}$.

1 Introduction

The paper targets the prohibitive cost of Metropolis-Hastings inference for large-scale, high-dimensional PDE applications and introduces a multilevel approach with a complexity analysis and subsurface-flow validation.

  • Motivation: High-dimensional subsurface-flow inference is expensive because likelihood evaluations require fine-grid PDE solves for many accepted and rejected MCMC proposals.Slow Monte Carlo convergence also requires many samples, while low acceptance rates can further increase overall complexity.
  • Motivation: Multilevel Monte Carlo reduces cost by using a hierarchy of models and estimating differences between successive model outputs rather than each output directly.The approach exploits linearity of expectation and cheaper coarse-grid computations.
  • Results: The paper reports significant computational-cost reductions for standard Metropolis-Hastings through the multilevel approach.The method builds on multilevel ideas previously successful for independent-sample Monte Carlo estimators in subsurface-flow problems.
  • Method: The proposed multilevel Metropolis-Hastings method uses level-dependent posteriors and parallel chains, with a new proposal distribution for the finer chain.This design avoids finest-level likelihood evaluations dominating every level and preserves the correct posterior at each level.
  • Analysis: The paper develops a problem-dependent complexity theorem and verifies its assumptions for a stationary single-phase subsurface-flow model with log-normal permeability.The theorem quantifies cost to achieve root mean square error below ε.
  • Results: O(ε^−(d+1)−δ) replaces O(ε^−(d+2)−δ) for the multilevel estimator, yielding one lower power of ε than the single-level method.For d = 2, numerical experiments confirm the theory and report absolute costs about O(10–50) times lower near ε = 10^−3.

2 Standard Markov chain Monte Carlo

The paper reviews standard Metropolis–Hastings MCMC for Bayesian inference, including its sampling mechanism, convergence conditions, and error decomposition. It emphasizes that large parameter dimensions, fine discretizations, expensive likelihood evaluations, and many samples make the estimator prohibitively costly, motivating a multilevel alternative.

  • Posterior sampling: Bayesian inference uses Metropolis–Hastings chains because posterior distributions are generally intractable and cannot be sampled exactly.The likelihood is evaluated through an approximate model response at a chosen discretization level.
  • Algorithm: Each Metropolis–Hastings iteration proposes θ′ from q(θ′|θn), accepts it with probability αM,R, or retains θn otherwise.The resulting states form a Markov chain whose samples are used for inference.
  • Convergence: The pCN proposal is important because random-walk convergence constants can grow linearly with the input dimension R, whereas judicious proposals can avoid this dependence.The proposal distribution determines the chain’s convergence behavior and is also used in the multilevel method.
  • Convergence: Under conditions including E ⊂ D and irreducibility, the chain has the desired stationary distribution and sample averages converge to EνM,R[QM,R].Convergence holds for arbitrary initial values and burn-in n0 ≥ 0 under the stated theorem conditions.
  • Error and cost: The MCMC mean-square error combines variance and squared bias with finite-sample and non-independent-sampling errors, whose bounds are difficult to establish.A central-limit result gives asymptotic decay comparable to independent sampling, but constants depend on chain covariance structure and burn-in assumptions.
  • Error and cost: Large discretization parameters M and R, expensive likelihood evaluations, and large sample counts make the standard estimator prohibitively costly.The paper introduces a multilevel estimator to overcome this computational cost.

3 Multilevel Markov chain Monte Carlo algorithm

The hierarchical multilevel MCMC estimator combines level-dependent posterior sampling with correlated two-level chains to reduce fine-level computational effort while preserving unbiasedness. An abstract complexity theorem quantifies the resulting ε-cost under verifiable assumptions, including MCMC, model, interlevel, and per-step cost conditions.

  • Multilevel estimator: MLMCMC estimates the finest-level expectation through a telescoping sum of expectations across posterior distributions that vary by level.When level distributions coincide, this reduces to the standard MLMC telescoping construction.
  • Convergence and complexity: The multilevel estimator is unbiased for EνL[QL] when the level chains have the required stationary marginals and irreducibility conditions hold.The construction preserves unbiasedness even though coarser posterior distributions need not model the data as faithfully as the finest posterior.
  • Multilevel estimator: Coarser levels use cheaper models, while decreasing correction variance permits fewer samples on fine levels.These two mechanisms jointly reduce computational effort relative to estimating the finest-level quantity directly.
  • Coupled chains: The estimator couples fine and coarse chains through correlated states while retaining νℓ and νℓ−1 as their stationary marginal distributions, avoiding bias in the telescoping sum.The fine chain uses a new proposal based on coarse-level sampling, whereas the coarse chain uses standard Metropolis-Hastings.
  • Coupled chains: Fine-level parameters are partitioned into coarse modes already present on level ℓ−1 and newly introduced fine modes.The fine proposal combines independently generated coarse- and fine-mode proposals before the usual Metropolis-Hastings accept/reject step.
  • Coupled chains: 18?
  • Convergence and complexity: The complexity theorem bounds ε-cost under assumptions on MCMC error, model error, interlevel strong error, and the cost of advancing a level.Compared with the i.i.d. multilevel case, MCMC introduces an additional factor of L + 1 in the sampling error and cost.
  • Convergence and complexity: Because L + 1 ≲ |log ε|, the ε-cost bounds inherit an additional |log ε| factor from the multilevel MCMC sampling analysis.This logarithmic factor follows from choosing the finest level so that the estimator bias is below ε^2/2.

4 Model Problem

The model problem applies multilevel MCMC to uncertainty quantification for single-phase subsurface flow with log-normal permeability, combining spatial and stochastic discretizations. The analysis verifies complexity assumptions and shows that multilevel MCMC saves one power of ε^-1 over standard MCMC under stated solver assumptions.

  • Model setting: The application concerns stationary single-phase subsurface flow, where uncertainty in pressure and velocity arises from a random permeability field.The governing model couples Darcy’s law with incompressibility, while permeability is inferred from limited geological information.
  • Model setting: The permeability is modeled as a positive log-normal random field whose logarithm has zero mean and exponential covariance with r = 1.This model permits broad permeability variation and is the setting used for the paper’s theoretical analysis.
  • Discretization: The computational model combines finite-element spatial discretization with a Karhunen–Loève expansion of log k using independent standard Gaussian variables.The spatial level controls degrees of freedom, while truncating the expansion controls the number of random variables.
  • Discretization: For the covariance model with r = 1, the KL eigenvalues satisfy μ_n ≲ n^-2, giving an intrinsic ordering of parameter importance and motivating truncation.Truncating after R terms approximates permeability using R standard normal random variables.
  • Complexity analysis: The analysis verifies assumptions M1 and M2 with α = β = 1/d − δ and α′ = β′ = 1/2 − δ for any δ > 0.These rates are established under the stated Gaussian-field and domain assumptions.
  • Complexity analysis: Assuming multigrid-like sample costs C_ℓ ≲ M_ℓ log(M_ℓ), the multilevel estimator saves one power of ε^-1 over standard MCMC.The paper also notes that the corresponding saving for i.i.d.-sample multilevel estimators is two powers of ε^-1 in dimensions d = 2, 3.

5 Numerics

The numerical study implements hierarchical MLMCMC for a finite-element subsurface-flow model and evaluates autocorrelation, sampling cost, multilevel variance, and convergence. The experiments report reduced fine-level sampling costs, negligible dependence bias under moderate subsampling, and agreement with predicted convergence behavior.

  • 5.1 Implementation Details: The model uses piecewise linear finite elements on uniformly refined triangular meshes, a truncated KL expansion for log-normal permeability, and UMFPACK linear solves.The mesh has Mℓ=(mℓ+1)^2 grid points, with mℓ+1=2^ℓm0+1.
  • 5.1 Implementation Details: The effective sample size Neff=N/τQ measures how many independent samples would yield the same variance as the correlated MCMC estimator.This quantity is used to compare sampling efficiency across levels.
  • 5.1 Implementation Details: Recursive independence sampling builds level-wise chains from coarser-level independent samples, while pCN proposals generate newly added fine modes.Subsampling rates tk and burn-in lengths are selected using estimated integrated autocorrelation times.
  • 5.3 Comparison of MLMCMC with a standard single-level MCMC estimator: The estimated bias decayed approximately as O(h), while the multilevel variance converged almost as O(h^2) and the rejection probability slightly faster than O(h).The authors caution that the observed rates may reflect pessimistic bounds or pre-asymptotic behavior.
  • 5.3 Comparison of MLMCMC with a standard single-level MCMC estimator: For the four-level estimator, about 562 fine-level PDE solves matched the sampling error of standard MCMC requiring about 4×10^6 fine-level solves.Although roughly 4×10^7 coarse auxiliary-chain solves were used, each was about 250 times cheaper than a level-4 solve.
  • 5.3 Comparison of MLMCMC with a standard single-level MCMC estimator: Recursive independence sampling also reduced the cost per independent sample for standard MCMC and provided an efficient burn-in strategy through shorter fine-level autocorrelation times.Most burn-in occurs on the coarsest level.

6 Conclusion

The conclusion presents MLMCMC as a hierarchy-based alternative for costly Bayesian inverse problems and reports its analysis and implementation for groundwater flow. Its applicability and gains depend on model-specific convergence and likelihood-cost rates, while multilevel acceptance must be verified for each proposal distribution.

  • 6 Conclusion: The method reduces cost by estimating differences in the quantity of interest between successive computational models rather than estimating each quantity independently.It was analyzed and implemented for a single-phase Darcy flow problem in groundwater modelling.
  • 6 Conclusion: The algorithm is generally applicable because it does not rely on a specific computational model, although model-dependent rates govern standard and multilevel costs.The relevant rates include discretisation-error exponents and the growth rate of likelihood-computation cost.
  • 6 Conclusion: The expected gain is larger for more challenging models with smaller convergence-rate parameters and larger likelihood-cost growth.This conclusion follows from the theorem’s dependence on α, α′, β, β′, and γ.
  • 6 Conclusion: The method supports various proposal distributions, but convergence of the multilevel acceptance probability to 1 must generally be verified for each proposal.The paper expects this property to hold for most proposal distributions.
Loading 1303.7343v2…