Randomized Hamiltonian Monte Carlo (RHMC)

• Randomized Hamiltonian Monte Carlo (RHMC) is a sampling method that uses exponentially-distributed durations for Hamiltonian evolution to replace fixed trajectory lengths.
• It employs stochastic momentum refreshment at random Poisson event times, yielding a piecewise deterministic Markov process that naturally preserves the target Gibbs distribution.
• RHMC offers robust convergence, dimension-free contraction properties, and extends to complex settings such as Riemannian manifolds and piecewise smooth targets.

Randomized Hamiltonian Monte Carlo (RHMC) is a family of Hamiltonian Markov chain Monte Carlo methods in which deterministic Hamiltonian motion is interrupted at random times, typically with momentum refreshment at Poisson event times or, in discrete-time formulations, with randomized integration durations. In its canonical form, RHMC targets the Boltzmann–Gibbs extension

$$\mu(dq,dp)=C_0^{-1}(2\pi)^{-D/2} e^{-H(q,p)}\,dq\,dp, \qquad H(q,p)=\frac{|p|^2}{2}+\Phi(q),$$

whose $q$-marginal is the desired target distribution. The central motivation is to replace the fragile choice of a deterministic trajectory length in standard HMC by a randomized timing mechanism that suppresses resonance, regularizes sampling efficiency as a function of the duration parameter, and yields a natural piecewise deterministic Markov process (PDMP) description (Bou-Rabee et al., 2015).

1. Canonical definition and algorithmic structure

In the canonical Euclidean formulation, the target distribution is

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

and RHMC augments $q\in\mathbb R^D$ with momentum $p\in\mathbb R^D$, using Hamilton’s equations

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

The basic modification relative to standard HMC is that the durations of Hamiltonian evolution are randomized rather than fixed. In the original RHMC construction, the durations between momentum refreshments are i.i.d. exponential random variables,

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

with mean $\lambda$. Between jump times, the process follows exact Hamiltonian flow; at each jump time, the momentum is randomized through the Horowitz refresh map

$$\Gamma(q,p)=\bigl(q,\cos(\phi)\,p+\sin(\phi)\,\xi\bigr),\qquad \xi\sim\mathcal N(0,I_D),$$

where $\phi=\pi/2$ gives full refreshment and smaller $q$0 gives partial refreshment (Bou-Rabee et al., 2015).

This idealized RHMC is naturally a PDMP. Its infinitesimal generator is

$q$1

which decomposes into Hamiltonian transport and a Poisson jump operator (Bou-Rabee et al., 2015). A closely related continuous-time formulation writes the process as Hamiltonian dynamics interrupted by the jump SDE

$q$2

where $q$3 is a homogeneous Poisson process of rate $q$4, $q$5, and $q$6 (Riou-Durand et al., 2022).

A broader class, Generalized Randomized Hamiltonian Monte Carlo (GRHMC), allows state-dependent event rates

$q$7

and jump kernels that preserve the event-rate-weighted momentum law

$q$8

Classical RHMC is recovered as the special case of constant-rate refreshment (Kleppe, 2020).

2. PDMP interpretation, invariance, and relation to HMC

RHMC differs from standard HMC in both timing and stochastic structure. Standard HMC usually chooses a deterministic integration time, applies numerical integration, and then uses a Metropolis correction. RHMC, by contrast, is continuous-time in its ideal form, with deterministic Hamiltonian motion between random Poisson event times and no accept/reject step in the exact-flow model. Randomness enters through event times and momentum updates rather than through repeated deterministic path lengths (Deligiannidis et al., 2018).

The invariant distribution is the canonical Gibbs law

$q$9

or, in the Euclidean mass-matrix notation of GRHMC,

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

For constant-rate RHMC, the invariant measure follows from two facts: Hamiltonian flow preserves energy and volume, and the refreshment kernel preserves the Gaussian momentum marginal. In the GRHMC extension with state-dependent event rates, the invariance mechanism is modified but remains explicit: the jump kernel must preserve the tilted momentum density $\Pi(dq)=C_0^{-1} e^{-\Phi(q)}\,dq,$1, which compensates exactly for the state-dependent event bias (Kleppe, 2020).

The PDMP perspective is not merely formal. It places RHMC alongside Bouncy Particle Sampler and Zig-Zag, but with a different deterministic component: RHMC follows nonlinear Hamiltonian ODEs between jumps, whereas BPS uses linear motion and state-dependent reflections. A notable consequence is that RHMC can be obtained as a high-dimensional weak limit of the first coordinate of BPS. Under strong convexity and smoothness assumptions,

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

the first location-velocity component of BPS converges weakly to a one-dimensional RHMC process with generator

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

as the ambient dimension tends to infinity (Deligiannidis et al., 2018).

This suggests a broader interpretation of RHMC as a canonical Hamiltonian PDMP: exact Hamiltonian transport provides conservative phase-space motion, while random refreshments induce the irreversibility and dissipation needed for sampling.

3. Geometric ergodicity, contraction, and dimension-free theory

The original RHMC paper proves geometric ergodicity in the exact-integration or small-step-size regime. The key tool is the Lyapunov function

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

for which

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

under coercivity assumptions on $\Pi(dq)=C_0^{-1} e^{-\Phi(q)}\,dq,$6. Combined with a minorization argument, this yields a Harris-type theorem: $\Pi(dq)=C_0^{-1} e^{-\Phi(q)}\,dq,$7 The paper emphasizes that these assumptions are essentially the same as those used for geometric ergodicity of underdamped Langevin dynamics (Bou-Rabee et al., 2015).

A later RHMC analysis establishes explicit dimension-free convergence-rate bounds under strong log-concavity and bounded Hessians,

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

Using synchronous coupling, the difference process satisfies

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

with $q\in\mathbb R^D$0, and contraction is proved in a weighted quadratic metric

$q\in\mathbb R^D$1

For a suitable refreshment rate

$q\in\mathbb R^D$2

the generator satisfies

$q\in\mathbb R^D$3

which yields exponential contraction in $q\in\mathbb R^D$4 and corresponding $q\in\mathbb R^D$5 semigroup bounds, with constants independent of the ambient dimension (Deligiannidis et al., 2018).

Another continuous-time RHMC result sharpens contraction bounds for strongly log-concave targets in a form directly comparable to kinetic Langevin diffusion. Under

$q\in\mathbb R^D$6

and with refreshment intensity

$q\in\mathbb R^D$7

the transition kernel $q\in\mathbb R^D$8 satisfies

$q\in\mathbb R^D$9

with rate

$p\in\mathbb R^D$0

The same work proves that, under the scaling

$p\in\mathbb R^D$1

RHMC converges at generator level to kinetic Langevin diffusion as $p\in\mathbb R^D$2, which identifies Langevin as a limiting high-persistence, high-frequency-refreshment regime of RHMC (Riou-Durand et al., 2022).

4. Randomized durations, resonance suppression, and Gaussian or quadratic models

The clearest analytic contrast between RHMC and fixed-duration HMC arises for Gaussian or quadratic targets. For the Gaussian target

$p\in\mathbb R^D$3

the original RHMC paper shows that, at refresh times,

$p\in\mathbb R^D$4

with $p\in\mathbb R^D$5 and $p\in\mathbb R^D$6. The integrated autocorrelation time for the mean of coordinate $p\in\mathbb R^D$7 is then

$p\in\mathbb R^D$8

whereas fixed-duration HMC gives

$p\in\mathbb R^D$9

The RHMC expression is smooth and monotone in $\dot q = p,\qquad \dot p = -\nabla \Phi(q).$0; the fixed-duration expression is oscillatory and singular at resonant values (Bou-Rabee et al., 2015).

This same resonance mechanism is central to nonasymptotic Gaussian-sampling theory for randomized integration lengths. For a Gaussian target

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

exact HMC in mode $\dot q = p,\qquad \dot p = -\nabla \Phi(q).$2 evolves as

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

If $\dot q = p,\qquad \dot p = -\nabla \Phi(q).$4 is fixed, some modes may satisfy $\dot q = p,\qquad \dot p = -\nabla \Phi(q).$5, producing periodic or near-periodic behavior. If $\dot q = p,\qquad \dot p = -\nabla \Phi(q).$6 is randomized from

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

then for every $\dot q = p,\qquad \dot p = -\nabla \Phi(q).$8,

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

and the product of mode-wise memory coefficients decays exponentially: $\delta t \sim \mathrm{Exp}(1/\lambda),$0 The resulting leapfrog HMC with randomized integration times and Metropolis correction attains

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

gradient evaluations to sample a $\delta t \sim \mathrm{Exp}(1/\lambda),$2-dimensional Gaussian to total variation error $\delta t \sim \mathrm{Exp}(1/\lambda),$3, improving on the $\delta t \sim \mathrm{Exp}(1/\lambda),$4 lower bound for HMC with fixed integration times (Apers et al., 2022).

A complementary quadratic analysis interprets the same phenomenon in terms of dissipation. For ideal HMC with randomized integration times

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

one has

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

Choosing

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

for quadratic targets with eigenvalues $\delta t \sim \mathrm{Exp}(1/\lambda),$8 yields an expected total integration time

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

which improves by a $\lambda$0 factor over classical constant-time, fully refreshed HMC. The same study argues that randomized integration time and partial refreshment should be viewed as two ways of controlling the effective dissipation level $\lambda$1 (Jiang, 2022).

5. Numerical methods, generalized event mechanisms, and geometric extensions

Exact Hamiltonian flow is rarely available outside special models, so practical RHMC implementations require numerical integration. One route is to approximate the continuous-time GRHMC process directly by adaptive ODE solvers. In Numerical GRHMC (NGRHMC), the Hamiltonian equations

$\lambda$2

are integrated with an explicit embedded sixth-order RKN pair, RKN6(4)6FD, together with adaptive error control

$\lambda$3

for recommended tolerances $\lambda$4. Event times are located by integrating the hazard

$\lambda$5

along the numerical trajectory until it hits an $\lambda$6 draw. The resulting process is approximate rather than exactly target-preserving, because there is no Metropolis correction, but the reported numerical bias is negligible relative to Monte Carlo error at these tolerances (Kleppe, 2020).

GRHMC also extends naturally to state-dependent event rates. A central example chooses

$\lambda$7

for which the integrated hazard equals the Mahalanobis arc length of the position trajectory. This specification makes event times correspond to exponentially distributed arc lengths rather than exponentially distributed elapsed times (Kleppe, 2020).

RHMC ideas have also been extended to piecewise smooth targets. For continuous densities with discontinuous gradients, the proposed numerical rule is to truncate each ODE step at a boundary crossing, switch to the new region-specific gradient, and continue integrating, thereby preserving the order of the adaptive Runge–Kutta solver. For piecewise smooth densities, GRHMC is augmented with reflection/refraction momentum updates at boundaries, including a randomized reflection rule

$\lambda$8

which preserves the required boundary momentum law (Tran et al., 25 Apr 2025).

A geometric analogue, Randomized Time Riemannian Manifold Hamiltonian Monte Carlo (RT-RMHMC), carries the same randomized-duration mechanism to manifolds and constrained submanifolds. On a Riemannian manifold $\lambda$9, the lifted target is

$$\Gamma(q,p)=\bigl(q,\cos(\phi)\,p+\sin(\phi)\,\xi\bigr),\qquad \xi\sim\mathcal N(0,I_D),$$0

with Hamiltonian

$$\Gamma(q,p)=\bigl(q,\cos(\phi)\,p+\sin(\phi)\,\xi\bigr),\qquad \xi\sim\mathcal N(0,I_D),$$1

The generator is

$$\Gamma(q,p)=\bigl(q,\cos(\phi)\,p+\sin(\phi)\,\xi\bigr),\qquad \xi\sim\mathcal N(0,I_D),$$2

where $$\Gamma(q,p)=\bigl(q,\cos(\phi)\,p+\sin(\phi)\,\xi\bigr),\qquad \xi\sim\mathcal N(0,I_D),$$3 fully refreshes tangent velocities from the metric Gaussian. For embedded manifolds, the practical method uses constrained RATTLE integration, a random total trajectory time $$\Gamma(q,p)=\bigl(q,\cos(\phi)\,p+\sin(\phi)\,\xi\bigr),\qquad \xi\sim\mathcal N(0,I_D),$$4, and Metropolis acceptance

$$\Gamma(q,p)=\bigl(q,\cos(\phi)\,p+\sin(\phi)\,\xi\bigr),\qquad \xi\sim\mathcal N(0,I_D),$$5

The paper proves invariance in continuous time and ergodicity for the discretized constrained chain under stated assumptions (Whalley et al., 2022).

6. Terminology, adjacent meanings of “RHMC,” and related methods

The acronym “RHMC” is not unique. In current literature, it may denote Randomized Hamiltonian Monte Carlo, Riemannian Hamiltonian Monte Carlo, or Relativistic Hamiltonian Monte Carlo, depending on context.

In computational statistics, “RHMC” frequently means Riemannian Hamiltonian Monte Carlo, i.e. HMC with a position-dependent metric

$$\Gamma(q,p)=\bigl(q,\cos(\phi)\,p+\sin(\phi)\,\xi\bigr),\qquad \xi\sim\mathcal N(0,I_D),$$6

A recent line of work uses this acronym to denote Hessian-informed Riemannian methods and studies how to reduce diagonal-metric fixed-point costs from $$\Gamma(q,p)=\bigl(q,\cos(\phi)\,p+\sin(\phi)\,\xi\bigr),\qquad \xi\sim\mathcal N(0,I_D),$$7 to $$\Gamma(q,p)=\bigl(q,\cos(\phi)\,p+\sin(\phi)\,\xi\bigr),\qquad \xi\sim\mathcal N(0,I_D),$$8 on coordinate-friendly targets. This usage is conceptually distinct from randomized-duration RHMC, even though both belong to the broader Hamiltonian MCMC family (Luu et al., 4 Jun 2026).

Another distinct usage is Relativistic Hamiltonian Monte Carlo, which replaces the quadratic kinetic energy by

$$\Gamma(q,p)=\bigl(q,\cos(\phi)\,p+\sin(\phi)\,\xi\bigr),\qquad \xi\sim\mathcal N(0,I_D),$$9

so that the induced velocity

$\phi=\pi/2$0

is bounded by $\phi=\pi/2$1. This “RHMC” is not randomized-duration HMC; its relevance is indirect, through the broader theme that modifying Hamiltonian structure can improve robustness (Lu et al., 2016).

Within the randomized-duration family itself, RHMC is closely connected to GHMC and MALT. One comparison paper treats RHMC as a robust alternative to standard HMC on anisotropic targets, using the exact Gaussian formula

$\phi=\pi/2$2

and derives lag-$\phi=\pi/2$3 correlations

$\phi=\pi/2$4

That work emphasizes both the robustness of RHMC to anisotropy and a limitation for even functions, since

$\phi=\pi/2$5

controls square-function correlations and does not vanish as $\phi=\pi/2$6 (Riou-Durand et al., 2022).

Taken together, these results establish Randomized Hamiltonian Monte Carlo as a distinct Hamiltonian-PDMP methodology defined by randomized flight times, typically exponential, with momentum refreshment at event times. Its most stable theoretical advantages arise in settings where deterministic path lengths create resonance or periodicity, and its modern extensions show that the same timing principle can be carried to state-dependent event mechanisms, manifold constraints, and piecewise smooth targets (Bou-Rabee et al., 2015).