Randomized Hamiltonian Monte Carlo

• Randomized Hamiltonian Monte Carlo is a family of methods that replaces deterministic trajectory lengths with random integration times, reducing resonance effects and enhancing ergodicity.
• The approach incorporates randomizations in step sizes, momentum refresh mechanisms, and kinetic parameters to achieve dimension-free convergence under conditions like strong log-concavity.
• Numerical implementations use techniques such as randomized leapfrog and PDMP formulations to reduce discretization bias and improve sampling performance in high-dimensional and constrained settings.

Searching arXiv for core papers on randomized Hamiltonian Monte Carlo and closely related formulations. Randomized Hamiltonian Monte Carlo denotes a family of Hamiltonian Monte Carlo variants in which the deterministic trajectory length of standard HMC is replaced or supplemented by randomization in the integration time, Poisson refresh times, step sizes, kinetic parameters, or the Hamiltonian itself. In the broader MCMC literature, “Randomized HMC” can mean random integration time, random step sizes, or other randomization layers in non-reversible samplers to avoid periodicity (Apers et al., 2022). Foundational formulations include a continuous-time RHMC method in which the durations between momentum randomizations are i.i.d. exponential random variables whose mean is a free parameter (Bou-Rabee et al., 2015), and a PDMP formulation in which exact Hamiltonian flow is altered at the arrival times of a homogeneous Poisson process by randomly perturbing the momentum component (Deligiannidis et al., 2018).

1. Foundational formulations

The basic HMC construction targets a density of the form

$\bar\pi(x)\propto e^{-U(x)}$

or, in equivalent notation,

$\Pi(dq)=C_0^{-1}\exp(-\Phi(q))\,dq.$

It augments the position variable with momentum and uses the Hamiltonian

$$H(x,v)=U(x)+\tfrac12\|v\|^2 \qquad\text{or}\qquad H(q,p)=\frac{|p|^2}{2}+\Phi(q),$$

so that the extended target is

$\pi(x,v)\propto e^{-U(x)-\|v\|^2/2}.$

Standard HMC chooses a trajectory length and integrator step size, integrates Hamilton’s equations for that deterministic duration, and then applies a Metropolis correction to remove discretization bias (Bou-Rabee et al., 2015).

The 2015 RHMC construction replaces the fixed duration by a random one. At time $t_0$, with state $Z_{t_0}=(Q_{t_0},P_{t_0})$, one draws

$\delta t_0 \sim \mathrm{Exp}(1/\lambda),$

evolves the exact Hamiltonian dynamics

$\dot q = p,\qquad \dot p = -\nabla \Phi(q)$

up to time $t_1=t_0+\delta t_0$, and then performs a momentum randomization

$$P_{t_1}=\cos\phi\,p(t_1)+\sin\phi\,\xi,\qquad \xi\sim\mathcal N(0,I_D),$$

where $\Pi(dq)=C_0^{-1}\exp(-\Phi(q))\,dq.$0 is the Horowitz angle (Bou-Rabee et al., 2015). In the small time step size limit, where the algorithm is rejection-free and the computational cost is proportional to the mean duration, the method is geometrically ergodic under the same conditions that imply geometric ergodicity of the solution to underdamped Langevin equations (Bou-Rabee et al., 2015).

A later formulation presents RHMC directly as a nonreversible PDMP. On $\Pi(dq)=C_0^{-1}\exp(-\Phi(q))\,dq.$1, with target

$\Pi(dq)=C_0^{-1}\exp(-\Phi(q))\,dq.$2

the generator is

$\Pi(dq)=C_0^{-1}\exp(-\Phi(q))\,dq.$3

where $\Pi(dq)=C_0^{-1}\exp(-\Phi(q))\,dq.$4 is the AR(1) momentum randomization kernel

$\Pi(dq)=C_0^{-1}\exp(-\Phi(q))\,dq.$5

This representation emphasizes that RHMC is a continuous-time Markov process, not a discrete-time chain (Deligiannidis et al., 2018).

2. Nonreversibility, scaling limits, and dimension-free convergence

RHMC is intrinsically nonreversible. In the PDMP formulation, the Hamiltonian drift part is anti-symmetric in $\Pi(dq)=C_0^{-1}\exp(-\Phi(q))\,dq.$6, while the refreshment part is symmetric and positive. This structure underlies both coupling-based and hypocoercive analyses (Deligiannidis et al., 2018).

Under the strong log-concavity and bounded Hessian assumption

$\Pi(dq)=C_0^{-1}\exp(-\Phi(q))\,dq.$7

dimension-free convergence rates can be proved. In particular, RHMC admits dimension-free exponential convergence in Wasserstein–2 and a dimension-free $\Pi(dq)=C_0^{-1}\exp(-\Phi(q))\,dq.$8 spectral gap under these assumptions (Deligiannidis et al., 2018). The same work establishes RHMC as the high-dimensional scaling limit of the first-coordinate dynamics of the Bouncy Particle Sampler: if BPS is initialized at stationarity, then the projected process

$\Pi(dq)=C_0^{-1}\exp(-\Phi(q))\,dq.$9

converges weakly to a one-dimensional RHMC process in Skorokhod space as the ambient dimension $$H(x,v)=U(x)+\tfrac12\|v\|^2 \qquad\text{or}\qquad H(q,p)=\frac{|p|^2}{2}+\Phi(q),$$0 (Deligiannidis et al., 2018). This places RHMC at the intersection of Hamiltonian Monte Carlo and the theory of nonreversible PDMP samplers.

The continuous-time viewpoint also clarifies the distinction from standard Metropolized HMC. Standard HMC requires a Metropolis step to correct discretization error, whereas ideal RHMC uses exact continuous-time dynamics with Poisson refreshes and no Metropolis step (Deligiannidis et al., 2018). That distinction is important: statements about exact invariance and rejection-free behavior in RHMC apply to the exact-flow PDMP, whereas discrete implementations inherit the usual numerical integration issues.

3. Long and randomized integration times for Gaussian targets

A particularly sharp complexity theory is available for Gaussian targets. For a $$H(x,v)=U(x)+\tfrac12\|v\|^2 \qquad\text{or}\qquad H(q,p)=\frac{|p|^2}{2}+\Phi(q),$$1-dimensional Gaussian with covariance matrix $$H(x,v)=U(x)+\tfrac12\|v\|^2 \qquad\text{or}\qquad H(q,p)=\frac{|p|^2}{2}+\Phi(q),$$2, written using the precision matrix $$H(x,v)=U(x)+\tfrac12\|v\|^2 \qquad\text{or}\qquad H(q,p)=\frac{|p|^2}{2}+\Phi(q),$$3,

$$H(x,v)=U(x)+\tfrac12\|v\|^2 \qquad\text{or}\qquad H(q,p)=\frac{|p|^2}{2}+\Phi(q),$$4

the Hamiltonian dynamics reduce, in the eigenbasis of $$H(x,v)=U(x)+\tfrac12\|v\|^2 \qquad\text{or}\qquad H(q,p)=\frac{|p|^2}{2}+\Phi(q),$$5, to independent harmonic oscillators with frequencies $$H(x,v)=U(x)+\tfrac12\|v\|^2 \qquad\text{or}\qquad H(q,p)=\frac{|p|^2}{2}+\Phi(q),$$6 when

$$H(x,v)=U(x)+\tfrac12\|v\|^2 \qquad\text{or}\qquad H(q,p)=\frac{|p|^2}{2}+\Phi(q),$$7

For fixed integration times, resonance effects create a fundamental lower bound: the paper cites that HMC with fixed integration times and leapfrog integrator requires at least

$$H(x,v)=U(x)+\tfrac12\|v\|^2 \qquad\text{or}\qquad H(q,p)=\frac{|p|^2}{2}+\Phi(q),$$8

gradient queries even for Gaussians (Apers et al., 2022).

The randomized alternative chooses a discrete set of long times

$$H(x,v)=U(x)+\tfrac12\|v\|^2 \qquad\text{or}\qquad H(q,p)=\frac{|p|^2}{2}+\Phi(q),$$9

and draws $\pi(x,v)\propto e^{-U(x)-\|v\|^2/2}.$0 uniformly from $\pi(x,v)\propto e^{-U(x)-\|v\|^2/2}.$1 at each step. The key anti-resonance property is

$\pi(x,v)\propto e^{-U(x)-\|v\|^2/2}.$2

which prevents the chain from locking into unfavorable periodic orbits (Apers et al., 2022). For the Metropolis-adjusted randomized leapfrog method, the resulting complexity is

$\pi(x,v)\propto e^{-U(x)-\|v\|^2/2}.$3

gradient queries to obtain an $\pi(x,v)\propto e^{-U(x)-\|v\|^2/2}.$4-accurate sample in total variation distance (Apers et al., 2022).

A complementary analysis reports an intriguing connection between variable integration time and partial velocity refreshment of ideal HMC samplers. On quadratic potentials, efficiency can be improved through these means by a $\pi(x,v)\propto e^{-U(x)-\|v\|^2/2}.$5 factor in Wasserstein-2 distance, compared to classical constant integration time, fully refreshed HMC (Jiang, 2022). In that account, randomized integration times, Chebyshev-designed deterministic variable times, and partial velocity refreshment are different ways of reducing the dissipative behavior of the dynamics. This suggests a unifying interpretation: randomized HMC is not merely trajectory-length jitter, but a mechanism for breaking the fixed rotation structure that produces poor conditioning dependence.

4. Numerical integration, Metropolization, and randomized integrators

Randomization enters not only through refresh times but also through numerical approximation of the Hamiltonian flow. In the Gaussian randomized-time analysis, leapfrog has a particularly simple structure: for each eigenvalue $\pi(x,v)\propto e^{-U(x)-\|v\|^2/2}.$6,

$\pi(x,v)\propto e^{-U(x)-\|v\|^2/2}.$7

and leapfrog exactly follows the Hamiltonian flow of a modified Hamiltonian

$\pi(x,v)\propto e^{-U(x)-\|v\|^2/2}.$8

The difference between the original and modified Hamiltonians is

$\pi(x,v)\propto e^{-U(x)-\|v\|^2/2}.$9

This identity supports both an unadjusted randomized-time leapfrog HMC kernel and a Metropolis-adjusted version whose acceptance probability can be written as a function of positions only (Apers et al., 2022).

A different line of work randomizes the integrator itself. Randomized Runge–Kutta–Nyström methods replace a deterministic symplectic integrator by a high-order randomized RKN scheme tailored for approximating Hamiltonian flows within unadjusted HMC and unadjusted kinetic Langevin MCMC. The paper introduces $t_0$0- and $t_0$1-order $t_0$2-accurate randomized Runge–Kutta–Nyström methods and establishes quantitative $t_0$3-order $t_0$4-accuracy upper bounds under gradient and Hessian Lipschitz assumptions on the potential energy function (Bou-Rabee et al., 2023). In this setting, the randomization is in the time quadrature inside each step, rather than in the Hamiltonian itself or in the refresh schedule.

These two strands address different numerical issues. Randomized trajectory length primarily targets resonance and conditioning, whereas randomized RKN integration targets discretization bias and strong approximation error in unadjusted samplers. They are therefore complementary rather than competing interpretations of randomized HMC.

5. Geometric, constrained, and generalized extensions

Randomized trajectory length extends naturally to Riemannian and constrained settings. On a $t_0$5-dimensional Riemannian manifold $t_0$6, with metric matrix $t_0$7, the manifold Hamiltonian is

$t_0$8

and the corresponding extended target is

$t_0$9

Randomized Time Riemannian Manifold HMC evolves the exact Hamiltonian flow on the tangent bundle between exponentially distributed event times and refreshes the velocity from the Riemannian Gaussian on $Z_{t_0}=(Q_{t_0},P_{t_0})$0. Its generator is

$Z_{t_0}=(Q_{t_0},P_{t_0})$1

where $Z_{t_0}=(Q_{t_0},P_{t_0})$2 refreshes momentum while holding position fixed (Whalley et al., 2022). In this setting, the conservation of the stationary distribution is proved in continuous time, and the Metropolized RATTLE-based discrete implementation is shown to be ergodic under geometric assumptions on the constraint manifold (Whalley et al., 2022).

For isometrically embedded constraint manifolds, the embedded Hamiltonian simplifies to

$Z_{t_0}=(Q_{t_0},P_{t_0})$3

with $Z_{t_0}=(Q_{t_0},P_{t_0})$4 restricted to the tangent space. The corresponding Randomized Time Constrained HMC samples tangent-space Gaussian velocities by orthogonal projection and then integrates constrained Hamiltonian dynamics with a RATTLE integrator (Whalley et al., 2022). Numerical studies in that work include spheres, Stiefel manifolds, and a high-dimensional covariance estimation problem.

Recent generalization targets nonsmooth structure directly. Generalized Randomized Hamiltonian Monte Carlo processes for sampling continuous densities with discontinuous gradient and piecewise smooth targets are proposed in order to combine the advantages of Hamiltonian Monte Carlo methods with the nature of continuous time processes in the form of piecewise deterministic Markov processes (Tran et al., 25 Apr 2025). This broadens the scope of randomized HMC beyond the differentiable setting required by traditional gradient-based samplers.

6. Other randomization mechanisms, interpretation, and limitations

Randomization can also be applied to kinetic parameters rather than to integration times. Quantum-Inspired Hamiltonian Monte Carlo allows a particle to have a random mass matrix with a probability distribution rather than a fixed mass (Liu et al., 2019). In that framework, each trajectory uses a draw

$Z_{t_0}=(Q_{t_0},P_{t_0})$5

and then performs a standard leapfrog trajectory with that mass. The stationary joint density over $Z_{t_0}=(Q_{t_0},P_{t_0})$6 is

$Z_{t_0}=(Q_{t_0},P_{t_0})$7

so the marginal in $Z_{t_0}=(Q_{t_0},P_{t_0})$8 remains

$Z_{t_0}=(Q_{t_0},P_{t_0})$9

provided the mass distribution is independent of the state (Liu et al., 2019). The same paper shows that, in the special case of implicit adaptation $\delta t_0 \sim \mathrm{Exp}(1/\lambda),$0, the time-averaged distribution equals the target, whereas in the general state-dependent case one generally has $\delta t_0 \sim \mathrm{Exp}(1/\lambda),$1 (Liu et al., 2019). This distinction separates admissible randomization from explicit adaptive schemes that alter the invariant law.

A still more abstract Gaussian theory studies HMC with random Hamiltonians. In that setting, the family of Gaussian distributions and their mixtures are invariant under Gaussian HMC operators, each such operator is a contraction on the space of parameters, and random sequences of such operators can be analyzed through stochastic fixed-point equations and Lyapunov exponents (Lu et al., 27 Feb 2026). The centered Gaussian update takes the form

$\delta t_0 \sim \mathrm{Exp}(1/\lambda),$2

making the contraction mechanism explicit (Lu et al., 27 Feb 2026).

Several persistent misconceptions can therefore be stated precisely. Randomized HMC is not a single algorithmic object; the term covers random integration time, random step sizes, and other randomization layers in non-reversible samplers (Apers et al., 2022). It is also not uniformly “Metropolis-free”: exact continuous-time RHMC dispenses with rejection because it uses exact Hamiltonian flow, whereas discrete-time implementations on Euclidean space or manifolds typically use leapfrog or RATTLE together with a Metropolis correction (Bou-Rabee et al., 2015). Finally, the strongest complexity improvements currently available are specific to Gaussian targets or to strongly log-concave targets with bounded Hessian, so claims of universal $\delta t_0 \sim \mathrm{Exp}(1/\lambda),$3-type acceleration beyond those settings remain interpretive rather than general theorems (Apers et al., 2022).