Metropolized Integrators for Exact Sampling

• Metropolized integrators are numerical methods that patch deterministic or stochastic proposals with a Metropolis test to exactly preserve the invariant Gibbs–Boltzmann measure.
• They ensure long-term stability by counteracting energy drift in stiff, high-dimensional, or hard-core systems using local moves and multiple-timestep techniques.
• Advanced schemes, including force-gradient corrections and optimized splitting, improve acceptance rates and computational efficiency in applications like molecular dynamics and lattice QCD.

Metropolized integrators are a class of numerical methods that combine explicit integrators for deterministic or stochastic differential equations (SDEs), such as those found in molecular dynamics (MD) or Markov chain Monte Carlo (MCMC), with a Metropolis–Hastings accept/reject criterion to ensure exact preservation of the desired equilibrium or invariant measure. Unlike traditional explicit integrators, these methods guarantee long-term stability and exact sampling of the target (often Gibbs–Boltzmann) measure, even when applied with large time steps or to stiff/high-dimensional systems. The key principle is to patch or modify an otherwise approximate deterministic or stochastic step with a global or local Metropolis test, yielding a discrete-time Markov process with rigorous equilibrium and ergodicity properties (Bou-Rabee, 2013, Bou-Rabee et al., 2010, Tupper et al., 2014).

1. Construction and Principles of Metropolized Integrators

The archetype of Metropolized integrators begins with a deterministic or stochastic dynamical law, e.g., Langevin SDEs,

$$dq = m^{-1}p\,dt, \quad dp = -\nabla U(q)\,dt - \gamma p\,dt + \sqrt{2\gamma k_B T m}\,dW_t,$$

with Hamiltonian $H(q,p) = \frac{1}{2} p^T m^{-1} p + U(q)$ and invariant Gibbs–Boltzmann measure $\nu(q,p)\propto \exp(-H(q,p)/k_B T)$. The integrator employs a splitting strategy:

• Ornstein–Uhlenbeck (OU) “bath” step on $p$ (solved exactly),
• Hamiltonian flow (local time-stepping, e.g., velocity-Verlet, for $(q,p)$),

followed by a Metropolis–Hastings acceptance step: $$\alpha = \min\left(1, \exp\left[-\beta (H(q_*,p_*) - H(q_0,p_0))\right]\right),$$ where $(q_*,p_*)$ is the proposal and $(q_0,p_0)$ the current state. Upon rejection, momentum is flipped to maintain detailed balance. The OU step is then applied to the resulting momenta. This procedure exactly preserves the invariant measure and is stable for arbitrary simulation lengths (Bou-Rabee, 2013, Bou-Rabee et al., 2010).

For SDEs with more general structure, including variable diffusion, proposals use the local structure, and the Metropolis test enforces exact detailed balance with respect to a prescribed equilibrium density (Tupper et al., 2014).

2. Stability, Accuracy, and Invariant Measure Preservation

Metropolized integrators are unconditionally stable in the sense that no matter the time step, the Markov chain they define leaves the target distribution invariant. For MD systems, they are provably stable under infinitely long simulation and circumvent the energy drift issues inherent to explicit integrators such as Verlet (Bou-Rabee, 2013). This is due to exact detailed balance, symplectic proposal mechanisms, and measure invariance.

On finite time intervals, such integrators are weakly convergent (e.g., order $h^{1/2}$ in the case of variable-diffusion SDEs (Tupper et al., 2014)), and, for proposals based on second-order integrators (e.g., velocity-Verlet), the rejection probability vanishes as $O(h^3)$, ensuring negligible bias for sufficiently small $h$ (Bou-Rabee, 2013, Bou-Rabee et al., 2010).

Even for stiff or hard-core potentials, where standard explicit integrators may catastrophically fail, the Metropolis test automatically rejects any proposals drifting into forbidden or singular regimes, preventing unphysical trajectories (Bou-Rabee, 2013, Bou-Rabee et al., 2010).

3. Scaling and High-Dimensional Performance

In high dimensions or large $N$-particle systems, global Metropolis moves cause the energy increment $\Delta H$ to grow extensively ($O(N)$), leading to exponentially vanishing acceptance rates unless the step size is reduced as $h \sim N^{-1/2}$ (Bou-Rabee, 2013, Bou-Rabee et al., 2010). To mitigate this scaling issue, several techniques are used:

• Local proposals: Updating only a subset (e.g., single atom/block) per Metropolis step keeps $\Delta H=O(1)$, allowing fixed $h$ independent of system size.
• Fast force evaluation: Efficient neighbor list or multipole algorithms are essential to maintain linear or near-linear computational cost in system size.
• Multiple-timestep algorithms and higher-order symplectic splittings: Decompose forces into stiff and soft contributions and use nested integration or locally adaptive step sizes to mitigate stiffness and improve acceptance rates (Kamleh, 2011).

These strategies yield scalable Metropolized approaches for large biomolecular systems and lattice field theory simulative frameworks.

4. Algorithmic Variants and Advanced Schemes

The core Metropolized framework has been extended in several significant directions:

• Multiple-timescale integrators (generalized nested leapfrog): In hybrid Monte Carlo (HMC) simulations, especially in lattice gauge theory, novel time-splitting arrangements allow tighly tuned, reversible, and symplectic composition of integrators, only constraining larger time steps to be integer multiples of the smallest one, thus enabling higher acceptance at fixed cost (Kamleh, 2011).
• Force-gradient and processed splitting integrators: By augmenting traditional symplectic integrator steps with explicit commutator or processing corrections, it is possible to suppress the leading-order shadow Hamiltonian error, thus dramatically raising the acceptance probability for a fixed number of force evaluations (Clark et al., 2010, Blanes et al., 2020, Blanes et al., 2014). For instance, symmetrically processed kernels can yield 3–5× greater acceptance per computation budget over Verlet/leapfrog in high-dimensions.
• Interaction with stochastic particle systems: Metropolization has been adapted to synergize with mean-field, ensemble, or interacting particle samplers (ALDI, CBS, stochastic SVGD), ensuring exact invariance and dramatically reducing discretization bias for complex Bayesian inverse problems (Sprungk et al., 2023). Particle- and block-wise Metropolis updates provide flexibility for parallel and scalable implementations.

5. Mixing Time Theory and Complexity Bounds

Rigorous mixing-time analyses show that, in log-concave or isoperimetrically regular settings, Metropolized Hamiltonian Monte Carlo (HMC) equipped with leapfrog/splitting integrators achieves faster convergence than Metropolis-adjusted Langevin (MALA) or random walk (MRW) samplers. Specifically, the total gradient-evaluation complexity for Metropolized HMC scales as $O(d^{11/12}\log(1/\epsilon))$, outperforming MALA ($O(d \log(1/\epsilon))$) and MRW ($O(d^2 \log(1/\epsilon))$) in high dimension (Chen et al., 2019). When only minimal isoperimetry and Hessian regularity are assumed, optimally tuning the number of leapfrog steps ($K \sim d^{1/4}$) and step size ($h \sim d^{-1/4}$) yields a mixing time of $\tilde O(d^{1/4})$, a strict improvement over the $d^{1/2}$ dependence of the best-known bounds for MALA in the same regime (Chen et al., 2023). This advantage is attributed to the ability of multi-step proposals to properly decorrelate, and to new KL divergence controls on proposal overlaps.

6. Practical Guidelines and Implementation Considerations

Efficient implementation depends on several factors:

• Step size and acceptance rate: One should maximize the step size subject to maintaining a moderate acceptance (e.g., 65–90% depending on dimensionality), often by leveraging per-coordinate or block-wise moves and employing optimized or processed integrator families over vanilla Verlet (Bou-Rabee, 2013, Blanes et al., 2014).
• Measurement-based tuning: Poisson-bracket estimators are used to tune free parameters in splitting schemes to maximize acceptance probability per CPU cost (Clark et al., 2010).
• Handling constraints and multiple time scales: Constraints (e.g., RATTLE for holonomic bonds) and multiple time steps (e.g., RESPA splits) can be incorporated at the proposal level, with post-proposal Metropolization ensuring equilibrium preservation even in the presence of high stiffness or bond constraints (Bou-Rabee et al., 2010, Kamleh, 2011).
• Scaling with system size: For extensive systems, local move strategies and parallelization (block or particle-wise Metropolis) are advised to avoid collapse of global step acceptance (Bou-Rabee, 2013, Sprungk et al., 2023).
• Processing and higher order: For very high dimension, symmetrically processed or multi-stage splitting integrators optimized for low energy error guarantee higher acceptance and practical efficiency gains (Blanes et al., 2020, Blanes et al., 2014).

7. Applications and Extension Scenarios

Metropolized integrators are utilized in diverse domains:

• Molecular dynamics for biomolecular and fluid simulations under thermal equilibrium, where unconditional long-term stability and exact sampling of the canonical ensemble are required (Bou-Rabee, 2013, Bou-Rabee et al., 2010).
• Hybrid Monte Carlo and Lattice QCD: Efficient sampling of high-dimensional field configurations using multi-timescale, force-gradient, and Metropolis-corrected integrators (Kamleh, 2011, Clark et al., 2010).
• Stochastic particle-based inference and Bayesian inversion: Ensemble samplers for complex posteriors metropolized to eliminate discretization bias and yield scalable and variance-reducing algorithms (Sprungk et al., 2023).
• General SDEs with variable diffusion: Integration techniques guaranteeing detailed balance and exact equilibrium for non-constant-diffusion SDEs (Tupper et al., 2014).

The above points collectively establish Metropolized integrators as a mature, scalable, and theoretically grounded approach for exact, long-time and high-dimensional simulation and sampling tasks across computational physics, statistics, and applied mathematics.