PaperlayerView exact revision on arXiv

Source-linked AI summary

MCMC Methods for Functions: Modifying Old Algorithms to Make Them Faster

S. L. Cotter, G. O. Roberts, A. M. Stuart, D. White

arXiv:1202.0709v3stat.COstat.ME

TL;DR

Standard MCMC methods become increasingly slow when function-valued posteriors are approximated on finer meshes. The paper modifies several methods using proposals from stochastic dynamics that preserve Gaussian reference measures, obtaining robust mesh-refinement behavior and substantial speed-ups across applications.

  • Problem

    Finite-dimensional MCMC approximations to infinite-dimensional function-valued posteriors suffer increasing convergence costs as the approximation dimension grows.

  • Method

    The paper constructs function-space MCMC proposals from carefully chosen stochastic dynamical-system discretizations that exactly preserve a Gaussian reference measure.

  • Results

    Small modifications of standard MCMC methods produce substantial speed-ups for Bayesian function estimation, with pCN mixing robustly as discretization dimension increases.

  • Takeaways & Limitations

    The Gaussian-reference framework applies across Bayesian nonparametric models and related settings, including conditioned diffusions and random-truncation priors.

  • Takeaways & Limitations

    Vanilla random-walk MCMC is not well-defined for functions and has mixing time tending to infinity as d_u increases; the target is assumed absolutely continuous with respect to a Gaussian reference measure.

Abstract

from arXiv · show

Many problems arising in applications result in the need to probe a probability distribution for functions. Examples include Bayesian nonparametric statistics and conditioned diffusion processes. Standard MCMC algorithms typically become arbitrarily slow under the mesh refinement dictated by nonparametric description of the unknown function. We describe an approach to modifying a whole range of MCMC methods, applicable whenever the target measure has density with respect to a Gaussian process or Gaussian random field reference measure, which ensures that their speed of convergence is robust under mesh refinement. Gaussian processes or random fields are fields whose marginal distributions, when evaluated at any finite set of $N$ points, are $\mathbb{R}^N$-valued Gaussians. The algorithmic approach that we describe is applicable not only when the desired probability measure has density with respect to a Gaussian process or Gaussian random field reference measure, but also to some useful non-Gaussian reference measures constructed through random truncation. In the applications of interest the data is often sparse and the prior specification is an essential part of the overall modelling strategy. These Gaussian-based reference measures are a very flexible modelling tool, finding wide-ranging application. Examples are shown in density estimation, data assimilation in fluid mechanics, subsurface geophysics and image registration. The key design principle is to formulate the MCMC method so that it is, in principle, applicable for functions; this may be achieved by use of proposals based on carefully chosen time-discretizations of stochastic dynamical systems which exactly preserve the Gaussian reference measure. Taking this approach leads to many new algorithms which can be implemented via minor modification of existing algorithms, yet which show enormous speed-up on a wide range of applied problems.

1. INTRODUCTION

The paper addresses the computational difficulty of sampling function-valued posteriors as finite-dimensional approximations are refined. It develops function-space MCMC methods based on Gaussian-reference-preserving proposals, yielding convergence robust to discretization dimension and substantial speed-ups.

  • 1. INTRODUCTION: Gaussian process priors provide flexible models for functions and support differentiability, physical-model structure, and computational specifications suited to differential equations.Their covariance or precision operators can be chosen according to statistical, local-structure, and numerical needs.
  • 1. INTRODUCTION: The framework covers target measures with densities relative to Gaussian random fields, conditioned diffusion processes, and random-truncation priors that can help prevent overfitting.Random truncation allows the data to determine the scales about which the model is informative.
  • 1. INTRODUCTION: Nonparametric inference concerns functions in infinite-dimensional spaces, so finite-dimensional approximations improve with dimension d_u but standard MCMC convergence slows as d_u grows.The paper seeks algorithms well-defined in the infinite-dimensional limit, whose discretizations retain robust convergence.
  • 1.1 Illustration of the Key Idea: The central design principle uses discretized stochastic dynamical systems whose proposals exactly preserve the Gaussian reference measure when the potential is zero.Such proposals do not reject in the purely Gaussian case and require only minor modifications to familiar MCMC algorithms.
  • 1.1 Illustration of the Key Idea: The modified pCN method remains robust as d_u increases, whereas standard random walk mixing worsens with dimension and mesh refinement.The disparity in mixing rates becomes greater as the mesh is refined.
  • 1.2 Overview of the Paper: The paper develops a common framework and function-space versions of multiple MCMC methods, then illustrates them on nontrivial Bayesian nonparametric applications.The stated goals include desirable d_u-independent mixing, applied demonstrations, and theoretical understanding of the methods’ benefits.

2. COMMON STRUCTURE

The paper formulates diverse function-valued Bayesian inverse and statistical problems as probability measures with densities relative to Gaussian reference fields. This framework also includes conditioned diffusions and random-truncation Gaussian priors.

  • Measure formulation: The common framework represents a target measure on a Hilbert space through a potential Φ relative to a random-field reference measure µ0.The potential is assumed numerically evaluable to the desired accuracy, with mesh refinement tied to increasing the finite-dimensional representation.
  • Applications: Nonparametric density estimation uses a Gaussian-process prior for the unnormalized log density, with the posterior expressed in the common form.The density is modeled through a transformation ensuring positivity and normalization.
  • Applications: Fluid data assimilation infers an initial velocity field from later velocity observations generated by Stokes or Navier–Stokes dynamics.The inverse problem may use Eulerian observations of the velocity field or Lagrangian particle trajectories.
  • Applications: Subsurface geophysics infers log-permeability from hydraulic-head measurements, while image registration infers momentum and reparameterization functions from noisy curve observations.Writing k(x)=exp(u(x)) enforces physically required positivity of permeability.
  • Conditioned diffusions: The same absolutely-continuous-to-Gaussian framework applies to conditioned diffusion problems with endpoint, continuous-path, or discrete path observations.The Girsanov formula expresses nonzero-drift path measures relative to zero-drift measures.

3. SPECIFICATION OF THE REFERENCE MEASURE

The reference measure is primarily a centered Gaussian µ0=N(0,C), represented through covariance eigenpairs, mesh-based exact sampling, or the precision operator. Random truncation extends the construction to useful non-Gaussian priors.

  • Gaussian reference measure: The principal reference law is the centered Gaussian measure µ0=N(0,C), with covariance operator C and, when available, precision operator L=C^-1.Efficient implementation requires accessible information about the Gaussian reference measure.
  • Gaussian reference measure: Gaussian draws can be generated from known covariance eigenpairs using Karhunen–Loève truncation, directly on a mesh, or through efficient precision-operator computations.These alternatives are not mutually exclusive and connect naturally to FFT, finite-element, and finite-difference methods.
  • Karhunen–Loève representation: The Karhunen–Loève representation expresses Gaussian functions through independent normal coefficients in an eigenfunction basis, with traceclass covariance ensuring L2 convergence.Finite-dimensional prior approximation uses projection onto the first d modes.
  • Karhunen–Loève representation: Mesh refinement for the prior means increasing the number du of Karhunen–Loève terms used to represent the target function.If the truncated series can be summed quickly on a grid, it provides efficient exact Gaussian samples.
  • Non-Gaussian extensions: Random truncation makes du random, switching expansion coefficients on or off and producing non-Gaussian, almost surely C∞ functions.Sieve priors generalize this by giving individual basis functions separate Bernoulli switches.

4. MCMC METHODS FOR FUNCTIONS

The paper redesigns MCMC proposals at function-space level using stochastic dynamical-system discretizations that preserve the Gaussian reference measure. The resulting methods remain well-defined and mix more robustly as discretizations are refined.

  • Design principle: Function-space MCMC design yields discretized algorithms whose convergence is robust as the representation dimension du tends to infinity.The approach is applied to random-walk, Langevin, independence, Gibbs, and HMC methods.
  • Design principle: The proposals arise from SPDE discretizations invariant for the reference or target measure, with γ=0 preserving µ0 and γ=1 preserving µ.The target behaves like the reference measure in high-frequency components.
  • Random walk limitation: Standard random walk proposals become singular in function space, causing all moves to be rejected and finite-dimensional mixing time to diverge as du increases.This explains why algorithms designed only after discretization deteriorate under mesh refinement.
  • Crank–Nicolson proposals: Crank–Nicolson proposals are designed to remain well-defined and irreducible on the Hilbert space, using either precision-operator inversion or Gaussian draws.Their proposal variance is aligned with the prior covariance to control the large-du limit.
  • Crank–Nicolson proposals: The pCN proposal significantly improves on naive random walk, while its acceptance behavior and mixing remain robust under mesh refinement.Theoretical results establish equivalence of forward and reverse measures for the CN and pCN proposals, unlike standard random walk.
  • Crank–Nicolson proposals: CN and pCN generalize finite-dimensional random walks to targets defined by density relative to a Gaussian, with acceptance based on differences in log density.When Φ is zero, the proposals exactly preserve the Gaussian reference measure and therefore do not reject in the trivial case.

4.3 MALA Proposal Distributions

The MALA extensions incorporate gradient information by discretizing stochastic dynamics invariant for the target measure rather than only the Gaussian reference. CNL and pCNL provide precision-based and prior-sampling implementations.

  • Target-invariant proposals: MALA proposals discretize an equation invariant for the target measure µ, incorporating steepest-descent information from the potential Φ.This contrasts with CN proposals, which use only the Gaussian reference measure µ0.
  • CNL: The Crank–Nicolson Langevin proposal CNL is accepted using the general Metropolis–Hastings acceptance rule and can be implemented by inverting I+ζL.The standard CNL method corresponds to θ=1/2.
  • pCNL: The preconditioned Crank–Nicolson Langevin proposal pCNL uses a draw w from the Gaussian reference measure µ0.Its implementation therefore requires efficient sampling from µ0, and the standard pCNL method also uses θ=1/2.

4.4 Independence Sampler

Setting δ = 2 in the pCN proposal produces an independence sampler that draws directly from the prior while retaining the stated acceptance probability.

  • δ = 2 turns the pCN proposal into an independence sampler with v = w, a draw from the prior.
  • The resulting independence sampler retains the acceptance probability given in (4.11).
  • An analogous choice, δ = 2 in the MALA proposal, provides a generalized independence-sampler construction.

4.5 Random Proposal Variance

The proposal variance δ can be randomized independently of the current state and auxiliary draw, while preserving proposal symmetry and the fixed-δ acceptance probability.

  • Randomizing δ according to any probability distribution ν on [0,∞), independently from w, defines a randomized pCN proposal.
  • For every δ, the measure η0(du,dv) = q(u,dv;δ)µ0(du) is well-defined and symmetric in u,v.
  • Both CN and pCN proposals can therefore use independently randomized δ values without changing the acceptance probability from (4.11).

4.6 Metropolis-Within-Gibbs: Blocking in Karhunen–Lo´eve Coordinates

The paper represents functions through Karhunen–Loève coefficients and uses coefficient blocks to construct Metropolis-within-Gibbs samplers defined at the function level. This supports robustness under mesh refinement, unlike standard singleton-coordinate partitioning.

  • Functions are expanded in a Karhunen–Loève basis, allowing the probability measure to be viewed as a measure on coefficients u = {ξi}∞.
  • Under the Gaussian prior, coefficient subsets ξI and ξI− are independent Gaussian variables with diagonal covariance operators.
  • Metropolis-within-Gibbs updates use partitions {Ij} covering N, with Φ expressed as a function of the coefficient expansion.
  • Function-defined blocking yields methods robust under mesh refinement when implemented, whereas standard singleton-coordinate samplers do not behave well.

4.7 Metropolis-Within-Gibbs: Random Truncation and Sieve Priors

Metropolis-within-Gibbs methods are extended to coefficient updates alongside variables controlling random truncation or active modes. The constructions include conditional sampling and reversible moves for the truncation variables.

  • The samplers alternate between updating function-expansion coefficients and parameters determining which modes are active.
  • For random-truncation priors, Metropolis-within-Gibbs updates the conditional distributions ξ|du and du|ξ, with ξ Gaussian under the prior and du distributed according to p(i).
  • A biased random walk for du|ξ can satisfy detailed balance with respect to p(i), after which the move uses the stated acceptance probability.
  • When p(i) decreases monotonically, local moves du → dv = du ± 1 provide a straightforward truncation proposal, though other stencils may improve mixing.
  • For sieve priors, the same framework samples ξ|χ and χ|ξ, with reversible χ|ξ proposals having the corresponding acceptance probability.

4.8 Hybrid Monte Carlo Methods

The paper develops a new HMC integrator whose Gaussian-case flow is exact on function space, preserving performance as discretizations grow arbitrarily large.

  • HMC introduces momentum or velocity variables alongside the state variable to construct a Hamiltonian flow.
  • The key innovation is an integrator for this Hamiltonian flow that is exact when Φ ≡ 0, even on function space.
  • Because the integrator remains exact in the Gaussian case, it behaves well for nonparametric problems with arbitrarily large or infinite-dimensional discretizations.

5. COMPUTATIONAL ILLUSTRATIONS

The computational illustrations compare standard and function-space MCMC methods across density estimation, data assimilation, subsurface geophysics, and image registration. They report faster mixing, convergence behavior dependent on proposal tuning and prior choice, and recovery of unknown functions and noise parameters.

  • 5.1 Density Estimation: pCN and RTM-pCN outperform MwG by an order of magnitude in density-estimation sampling.Their integrated autocorrelation times, trace plots, and correlation functions show faster performance; MwG is heavily mesh-dependent because it updates one Fourier component at a time.
  • 5.1 Density Estimation: RTM-pCN has approximately double pCN’s asymptotic variance but can reduce runtime per unit error through adaptive basis size.The reported improvement is primarily attributed to fewer random-number generations.
  • 5.2 Data Assimilation in Fluid Mechanics: Different pCN initial states converge to the same distribution in the Lagrangian data-assimilation experiment.The proposal variance was selected for an average acceptance probability of approximately 25%.
  • 5.2 Data Assimilation in Fluid Mechanics: Randomizing β produces roughly comparable convergence to static βopt in the Eulerian data-assimilation experiment.The randomized proposal mixes values from 0.1 × βopt to 1.9 × βopt and may help explore multimodal posteriors through both large and small steps.
  • 5.3 Subsurface Geophysics: Only the sieve prior/Sieve-pCN combination converges within the available computational time in the under-determined subsurface inverse problem.The MwG and Gaussian-prior pCN chains fail to converge under these test conditions, highlighting the importance of prior modelling assumptions.
  • 5.4 Image Registration: Increasing the number of observations makes the posterior increasingly peaked near the true observational variance 0.01 while recovering the true function and precision parameter.This is presented as a form of posterior consistency in the image-registration problem.

6. THEORETICAL ANALYSIS

The theoretical analysis explains why Crank–Nicolson-based proposals are well-defined on function space and retain favorable acceptance behavior, while standard random-walk proposals do not. It also connects these properties to mesh-robust MCMC convergence.

  • 6. THEORETICAL ANALYSIS: Crank–Nicolson-type function-space algorithms are supported by theory explaining their behavior and acceptance properties in nonparametric settings.The analysis complements numerical demonstrations with function-space results.
  • 6. THEORETICAL ANALYSIS: θ = 1 is the unique choice making the forward and reverse proposal measures equivalent under the stated assumptions.This equivalence is required for a well-defined function-space MCMC method.
  • 6. THEORETICAL ANALYSIS: The Crank–Nicolson proposal preserves the Gaussian prior measure, since u ∼ N(0,C) implies v ∼ N(0,C).This property mirrors the invariant measure of the underlying stochastic differential equation.
  • 6. THEORETICAL ANALYSIS: Standard random-walk proposals are not defined on function space because the forward and reverse joint measures are not absolutely continuous.The result applies for both K = I and K = C.
  • 6. THEORETICAL ANALYSIS: Under the stated assumptions, both pCN and CN algorithms are defined on the Hilbert space for fixed δ, with acceptance probability behaving continuously as δ → 0.This extends finite-dimensional intuition to the function-space setting.
  • 6. THEORETICAL ANALYSIS: Scaling-limit analyses show that standard MCMC proposal variance must shrink with discretization dimension, increasing the number of steps under mesh refinement.The cited literature relates this scaling to optimal proposal variance and acceptance probability.

7. CONCLUSIONS

The conclusions emphasize designing MCMC methods directly on function space and then discretizing them, using Gaussian-measure-preserving proposals. The resulting methods require small changes to existing algorithms, apply across diverse problems, and show efficacy against standard methods.

  • 7. CONCLUSIONS: Many applications naturally produce posterior problems defined relative to Gaussian random field reference measures or related structures.This includes the broad application setting motivating the framework.
  • 7. CONCLUSIONS: Designing MCMC methods on function space before discretization provides better insight into nonparametric algorithm design than discretizing first and applying standard methods.
  • 7. CONCLUSIONS: Gaussian-reference proposals should accept with probability one when sampling only the reference measure, motivating stochastic-dynamical-system discretizations that preserve it.
  • 7. CONCLUSIONS: The framework yields random-walk, Langevin, and Hybrid Monte Carlo Metropolis methods that can be implemented through small modifications of existing codes.
  • 7. CONCLUSIONS: Applications demonstrate efficacy relative to standard methods and compatibility with Gibbs sampling and conjugate Gamma-prior noise-precision estimation.
  • 7. CONCLUSIONS: Gaussian-prior ubiquity makes the technology immediately applicable across many applications, while the underlying philosophy suggests extensions to non-Gaussian priors.
Loading 1202.0709v3…