HMC Algorithm

WarningDeprecated

The fixed-length HMC sampler (update_method = "hamiltonian-mc") is deprecated and will be removed in a future release. Use update_method = "nuts" instead — NUTS adapts trajectory length automatically and avoids the need to tune the number of leapfrog steps.

Hamiltonian Monte Carlo (HMC) simulates Hamiltonian dynamics to propose distant moves in parameter space while maintaining a high acceptance rate. The bgms implementation uses fixed-length leapfrog integration with a Metropolis accept/reject correction.

The implementation is in src/mcmc/algorithms/hmc.h and hmc.cpp. The sampler wrapper is HMCSampler in src/mcmc/samplers/hmc_sampler.h, which inherits warmup adaptation from GradientSamplerBase.

Algorithm

Each HMC iteration proceeds as follows:

  1. Sample momentum — Draw \(p \sim \mathcal{N}(0,\, M)\), where \(M\) is the diagonal mass matrix.
  2. Leapfrog integration — Simulate \(L\) steps of the Störmer–Verlet integrator with step size \(\epsilon\):
    • Half-step momentum: \(p \leftarrow p + \tfrac{\epsilon}{2}\,\nabla\log \pi(q)\)
    • Full-step position: \(q \leftarrow q + \epsilon\,M^{-1} p\)
    • Half-step momentum: \(p \leftarrow p + \tfrac{\epsilon}{2}\,\nabla\log \pi(q)\)
  3. Metropolis correction — Accept the proposal \((q^*, p^*)\) with probability \[ \alpha = \min\!\bigl(1,\; \exp\bigl(H(q, p) - H(q^*, p^*)\bigr)\bigr), \] where \(H(q, p) = -\log \pi(q) + \tfrac{1}{2}\,p^\top M^{-1} p\) is the Hamiltonian.

The number of leapfrog steps \(L\) is controlled by hmc_num_leapfrogs (default 100 from the R interface, configurable via SamplerConfig::num_leapfrogs).

Why NUTS is preferred

HMC requires \(L\) to be set in advance. A value that is too small wastes gradient evaluations on short trajectories; a value that is too large wastes them on trajectories that loop back. NUTS eliminates this tuning parameter by growing a binary tree until a U-turn is detected, adapting the trajectory length to the local geometry of the posterior. See NUTS Algorithm.

Constrained variant

For models with edge selection on continuous variables (GGM, Mixed MRF), HMCSampler dispatches to a constrained integration path that uses RATTLE projections at each leapfrog step. This keeps excluded precision entries exactly zero throughout the trajectory. The constrained path includes an optional reversibility check identical to the one used by NUTS (see Constrained Leapfrog (RATTLE)).

A fixed trajectory length interacts poorly with the constraint manifold: the integrator cannot shorten the trajectory when it encounters high-curvature regions or non-reversible projections. NUTS handles this by terminating the tree early. When update_method = "hamiltonian-mc" is used with edge selection, bgms issues a warning:

hamiltonian-mc with edge selection on a GGM uses constrained integration (RATTLE), which can be numerically fragile with a fixed trajectory length. Consider using ‘nuts’ instead, which adapts trajectory length and avoids degenerate regions.

Step-size heuristic

Before warmup begins, heuristic_initial_step_size() finds a reasonable starting step size by iteratively doubling or halving a candidate \(\epsilon\) until a single leapfrog step yields an acceptance probability near the target (default 0.8). The algorithm:

  1. Evaluate the Hamiltonian at the current position with a fresh momentum draw.
  2. Run one leapfrog step and evaluate the proposed Hamiltonian.
  3. If \(\alpha > \alpha_{\text{target}}\), double \(\epsilon\); otherwise halve it.
  4. Repeat until the acceptance probability crosses the target or max_attempts (default 20) is reached.

Three overloads exist:

Variant Mass matrix Integration
Identity \(M = I\) Unconstrained
Diagonal User-supplied \(M^{-1}\) Unconstrained
Constrained User-supplied \(M^{-1}\) RATTLE projections

The constrained overload (heuristic_initial_step_size_constrained) projects momentum onto the cotangent space before evaluating kinetic energy and uses leapfrog_constrained() for the trial step.

Warmup adaptation

HMCSampler inherits from GradientSamplerBase, which owns an HMCAdaptationController. The adaptation schedule is identical to NUTS (see Warmup Schedule):

  • Stage 1 — Fast step-size adaptation via dual averaging.
  • Stage 2 — Diagonal mass matrix estimation in doubling windows, with step-size re-initialization after each window.
  • Stage 3a — Final step-size tuning with the learned mass matrix.
  • Stage 3b — Proposal-SD tuning for edge selection (model-driven).
  • Stage 3c — Step-size re-adaptation after edge-selection structure changes.

The default target acceptance rate for HMC is 0.65 (vs. 0.80 for NUTS).

Sampler dispatch

The chain runner factory in chain_runner.cpp maps two string identifiers to HMCSampler: "hmc" and "hamiltonian-mc". The R interface uses "hamiltonian-mc" as the update_method value. On the R side, validate_sampler() issues a .Deprecated() warning when this method is selected.

Source files

File Contents
src/mcmc/algorithms/hmc.h kinetic_energy(), heuristic_initial_step_size() (3 overloads), hmc_step() (2 overloads)
src/mcmc/algorithms/hmc.cpp Implementations of all HMC algorithm functions
src/mcmc/samplers/hmc_sampler.h HMCSampler class (header-only)
src/mcmc/samplers/hmc_adaptation.h DualAveraging, DiagMassMatrixAccumulator, HMCAdaptationController