Improved MCMC Algorithms for Function Space Sampling

• Improved MCMC algorithms are advanced techniques that modify traditional methods to achieve robust, mesh-free mixing in infinite-dimensional and high-dimensional settings.
• They employ function-space matched proposals such as the pCN method to preserve Gaussian reference measures, ensuring consistent performance during mesh refinement.
• These approaches enhance Bayesian inverse problem solving and data assimilation by reducing autocorrelation times and mitigating challenges like the curse of dimensionality.

Improved Markov Chain Monte Carlo (MCMC) algorithms refer to a spectrum of approaches that modify, extend, or redesign traditional MCMC methods—typically Metropolis–Hastings or Gibbs sampling—to accelerate mixing, ensure robust performance in challenging regimes (such as high dimensions, function-space targets, or multimodal posteriors), and improve computational efficiency per unit sample. In contemporary literature, “improved” often means addressing bottlenecks like mesh-refinement degeneracy, curse of dimensionality, poor scaling in the presence of hierarchical or spatial structure, and inefficiency due to naive proposal mechanisms or adaptation rules. An important subclass of such improvements involves algorithmic innovations that preserve target measures under mesh refinement, construct proposals directly in function space or infinite-dimensional Hilbert space, and adopt SPDE-discretization-based proposal mechanisms.

1. Function-Space-Matched Proposal Mechanisms

A foundational innovation in improved MCMC algorithms for function-valued unknowns is the design of proposal distributions that are “compatible” with the Gaussian process or Gaussian random field prior. When the target is absolutely continuous with respect to a centered Gaussian reference measure μ₀ (with covariance operator 𝒞), i.e.,

$d\mu(u) \propto \exp(-\Phi(u))\, d\mu_0(u),$

the corresponding MCMC algorithm should, when $\Phi \equiv 0$, preserve μ₀ exactly. Classical finite-dimensional random walks or Langevin proposals, e.g., $v = u + \sqrt{2\delta}\,\xi_0$ with $\xi_0 \sim N(0, I)$, become ineffective in the high-dimensional limit because the forward and reverse proposal distributions become mutually singular as the discretization is refined.

The preconditioned Crank–Nicolson (pCN) proposal addresses this by constructing moves as

$$v = \sqrt{1 - \beta^2}\,u + \beta w,\qquad w \sim \mathcal{N}(0, \mathcal{C}),$$

where $\beta \in [0, 1]$ is related to the time discretization step via a Crank–Nicolson discretization of an infinite-dimensional Ornstein–Uhlenbeck (OU) process. Such proposals exactly preserve the reference measure μ₀, ensuring that acceptance rates and mixing do not deteriorate upon mesh refinement—an essential property for algorithms designed to operate on function space.

2. SPDE Discretization and Acceptance Probability Structure

The link between SPDEs and proposal design is formalized by recognizing that many infinite-dimensional Bayesian inverse problems or nonparametric models can be expressed as sampling from the invariant measure of an SPDE such as

$$du/ds = -\mathcal{K}(\mathcal{L}u + \gamma D\Phi(u)) + \sqrt{2\mathcal{K}}\, db/ds,$$

where $\mathcal{L}$ is the precision operator (inverse of the covariance), $D\Phi(u)$ is the Fréchet derivative of the potential, and $\mathcal{K}$ prescribes the direction of exploration (often set to $I$ or $\mathcal{C}$). Discretizing this dynamic with, e.g., a Crank–Nicolson (θ = 1/2) scheme ensures that for $\Phi \equiv 0$ the transition kernel is reversible with respect to μ₀ and the forward/reverse proposal measures are equivalent.

Under nontrivial $\Phi$, the Metropolis–Hastings acceptance probability simplifies to

$$a(u, v) = \min\left\{1, \exp[\Phi(u) - \Phi(v)]\right\},$$

which is dimension-independent due to the structure of the proposal. This contrasts with standard random walk or Langevin proposals, for which the acceptance probability degenerates to zero as the number of function degrees of freedom increases.

3. Robustness Under Mesh Refinement

A critical motivation for these improved algorithms is the mesh refinement necessary for increasingly accurate function representation, e.g., via truncated Karhunen–Loève expansions:

$u(x) = \sum_{i=1}^{d_u} \xi_i \varphi_i(x).$

As $d_u$ increases, classical Metropolis-within-Gibbs or random-coefficient updating schemes succumb to the curse of dimensionality, with integrated autocorrelation times (IACT) escalating rapidly. By contrast, pCN and related SPDE-discretization-based proposals preserve global balance with respect to the infinite-dimensional reference measure and maintain stable IACT and acceptance rates as $d_u \to \infty$.

4. Applications to Bayesian Nonparametrics and Inverse Problems

Improved MCMC schemes have direct and substantial implications for applied domains:

a) Density estimation: When modeling an unknown density via $\rho(x) = \exp(u(x)) / \int \exp(u(s))ds$ with a Gaussian process prior on $u$, classical coefficient-wise updates are inefficient as basis size grows. Using pCN, the entire function is updated coherently, respecting the geometry of μ₀ and dramatically reducing the IACT.

b) Data assimilation in fluid mechanics: Inverse problems for PDE initial conditions (e.g., recovering $u(x)$ in Navier–Stokes from velocity field measurements) yield a posterior over infinite-dimensional function spaces. The pCN proposal leads to acceptance rates and mixing times that do not degrade as the discretization is refined, evidenced by numerical studies in both Eulerian and Lagrangian observation schemes.

5. Algorithmic Generalization and Extensions

The proposal framework admits multiple generalizations:

• Gradient-informed pCN (pCN-Langevin or pCNL): Incorporating gradient information via $\gamma = 1$ in the SPDE discretization.
• Independence samplers: Obtained from the pCN family by taking $\delta = 2$ in the Crank–Nicolson scheme.
• Randomized step sizes or more complex updating schedules: Possible within the invariance-preserving discretization framework.
• Gibbs sampling and reversible-jump extensions: Within basis coefficient representations, allowing further adaptation to models with variable dimension or random truncation priors.

6. Numerical Performance and Empirical Evidence

Empirical studies reported in density estimation and data assimilation confirm order-of-magnitude improvements in integrated autocorrelation time—mixing times remain nearly constant even as function discretization is refined. Acceptance rates are near unity for the prior (Φ ≡ 0) and remain high in the presence of data, while classical algorithms degrade sharply as dimension grows.

7. Synthesis and Impact

Improved MCMC algorithms in the sense articulated here are characterized by proposals constructed to be invariant with respect to the function-space reference measure (typically Gaussian). This invariance enforces mesh-free mixing and robust acceptance probabilities, which are vital for the practical feasibility of Bayesian inference and uncertainty quantification in infinite-dimensional settings. The approach is widely applicable, allowing integration of advanced MCMC schemes (Gibbs, Langevin, reversible-jump) into the function-space framework and supporting the increasing scientific demand for scalable MCMC in high- and infinite-dimensional models encountered in modern statistical modeling and inverse problems (Cotter et al., 2012).