Pith. sign in

REVIEW 8 minor 22 references

The Physicist's Guide to the HMC

T0 review · 0 major / 8 minor · reviewed 2026-08-09 · deepseek-v4-flash

Challenge this review Re-run · record.json Download PDF Read on arXiv ↗

Pith's one-line read The paper claims that hybrid Monte Carlo with exact Fourier acceleration becomes equivalent to direct sampling for quadratic actions, and remains viable for generic actions when augmented with long randomised trajectories and radial…

desk verdict A clear, honest recipe for HMC with exact Fourier acceleration; no new science, but the algorithms are correct and the guide is genuinely useful. read the letter →

arxiv 2501.19130 v1 pith:POSHSVPW submitted 2025-01-31 hep-lat cond-mat.str-elphysics.comp-phstat.CO

Johann Ostmeyer This is my paper · ORCID
classification hep-latcond-mat.str-elphysics.comp-phstat.CO PACS 02.70.Uu11.15.Ha
keywords hybridMonteCarloexactFourieraccelerationautocorrelationtimelatticefieldtheoryradialupdatesHamiltonianGaussiansamplingMCMC
verification ladder T0 review T1 audit T2 compute T3 formal

The pith

A machine-rendered reading of the paper's core claim, the machinery that carries it, and where it could break.

The reading

The paper aims to establish a practical recipe for hybrid Monte Carlo (HMC) sampling that reaches near-optimal efficiency for a wide class of probability distributions. Its central claim is that coupling HMC with exact Fourier acceleration (EFA) makes the sampler equivalent to direct sampling for quadratic actions $S(x)=\frac12 x^\mathsf{T} M x$, i.e. normal distributions, and only perturbatively worse when a small anharmonic perturbation $V(x)$ is added. The author argues that EFA, together with randomised sufficiently long trajectories and radial updates on non-compact spaces, keeps the algorithm viable for arbitrary continuous distributions. The argument is supported by the analytical results of the companion reference and by numerical benchmarks on the SSH model, the Ising model, and lattice gauge theories, all of which show integrated autocorrelation times near the theoretical floor.

What carries the argument

The central object is the split Hamiltonian $H=\frac12 p^\mathsf{T} M^{-1} p + S(x)$, with the exact Fourier acceleration (EFA) step of Algorithm 3 as the load-bearing mechanism. EFA uses the eigendecomposition $M=\Omega\, \mathrm{diag}(\omega^2)\, \Omega^\dagger$ (often a Fourier transform for translationally invariant systems) to evolve harmonic eigenmodes exactly: after rotating to the eigenbasis, each mode is advanced by the exact rotation $y_i(h)=\cos(h)\, y0_i + (1/\omega_i^2)\, \sin(h)\, q0_i$ and $q_i(h)=\cos(h)\, q0_i - \omega_i^2\, \sin(h)\, y0_i$. This eliminates the sampling error coming from the quadratic part, leaving only the anharmonic perturbation to be handled by the integrator; the leap-frog step and the Metropolis acceptance rule then convert the near-conservation of $H$ into a high acceptance rate and minimal autocorrelation.

What would settle it

Simulate a translationally invariant Gaussian target on a $32\times32$ lattice using EFA and measure the integrated autocorrelation time $\tau_{\mathrm{int}}$ of a nonlinear observable such as $\langle x^4\rangle$; if $\tau_{\mathrm{int}}$ grows with lattice size instead of staying near its minimum of 0.5, the equivalence to direct sampling fails. Alternatively, for a strongly anharmonic action, compare EFA with plain HMC at equal compute time; if EFA's acceptance or effective sample size is not better, the claimed advantage has not materialised.

Watch

Extended reading notes

Core claim

For a distribution $P(x)\propto \mathrm{e}^{-S(x)}$ with harmonically dominated action $S=\frac12 x^\mathsf{T} M x + V(x)$, the paper's claim is that HMC with EFA samples the harmonic part exactly: the molecular dynamics solves the equation of motion for each eigenmode of $M$ with a closed cosine/sine step, so $\Delta H$ comes only from the anharmonic perturbation. With trajectory length $T=\pi/2$ and acceptance near one, the integrated autocorrelation time stays close to its minimum and the algorithm behaves like direct sampling. The claim extends perturbatively, so small deviations from the Gaussian case are only mildly worse, and it remains a viable sampler for arbitrary actions when three adjustments are made: randomised trajectory length, sufficiently long trajectories, and radial updates for non-compact spaces. Algorithms 1–4 implement this as a turn-key numerical recipe.

Load-bearing premise

The recipe's efficiency rests on knowing the harmonic matrix $M$ and being able to diagonalise it (ideally by a Fourier transform), and on the harmonic part dominating the action so that the residual $V(x)$ is a small perturbation.

Editorial extensions

If this is right

  • With EFA, users of HMC on Gaussian and near-Gaussian targets can set the trajectory length to $T=\pi/2$ and use a single leap-frog step, obtaining acceptance rates that stay in the 65–80% range and integrated autocorrelation times below about 0.7.
  • For translationally invariant systems, the eigendecomposition of $M$ is a Fourier transform, so EFA adds almost no overhead and applies directly to lattice problems with minimal extra cost.
  • Randomising the trajectory length and choosing it long enough on average prevents the slow-decorrelation modes that appear when the eigenfrequencies of the HMC dynamics are unknown.
  • On non-compact spaces, radial updates restore exponential convergence and can overcome barrier-induced ergodicity problems, completing the recipe for generic actions.

Reading between the lines

Editorial extensions of the paper, not claims the author makes directly.

  • For problems without translational symmetry, the eigendecomposition of $M$ is a genuine one-time cost that the paper does not budget for; the method's advantage is most certain when $M$ is diagonalisable by a fast transform.
  • A quantitative follow-up would be an analytic bound on acceptance-rate or autocorrelation loss as a function of the anharmonic fraction $|V(x)|/(x^\mathsf{T} M x)$, giving users a threshold for when EFA is worth the extra implementation effort.
  • Folding EFA into adaptive HMC samplers that already randomise path length (such as the No-U-Turn sampler) could yield a fully automatic near-optimal sampler for generic actions, though the paper stops short of testing that combination.
Share X Bluesky LinkedIn Reddit HN

Editorial analysis

A structured set of objections, weighed in public.

Desk editor's note, referee report, and a circularity audit.

Referee Report

0 major / 8 minor

Summary. The paper is a compact procedural guide to hybrid Monte Carlo with exact Fourier acceleration (EFA). It states that for actions S(x)=½x^T M x + V(x) with V(x) small, HMC with EFA solves the harmonic part exactly in the eigenbasis of M (Algorithms 1–3), and that with trajectory length T=π/2 each trajectory produces a sample from the quadratic target with acceptance probability one. It then collects adjustments needed for generic actions: randomising trajectory lengths, choosing long trajectories, including radial updates on non-compact spaces, and regularising the kinetic term for strongly anharmonic actions. Numerical illustrations are given for the SSH model, the 2D Ising model, and pure gauge theory, and the paper closes with a summary of the recommended practical recipe.

Significance. If the central claims hold, the recipe gives a concrete and simple way to achieve near-optimal sampling for near-Gaussian high-dimensional targets, which is directly relevant for lattice field theory and condensed-matter simulations. The algorithmic core is correct, and the self-contained presentation of Algorithms 1–4, together with the pointer to open-source code in Ref. [2], is a useful service to practitioners. The main caveats are that the 'only perturbatively worse' statement is not derived here but is delegated to Ref. [2], that the efficiency claim for arbitrary actions is explicitly limited in Sec. 3.1.4, and that the numerical evidence is presented without error bars or full parameter reporting. These caveats affect the strength of the advertised claims but not the correctness of the quadratic-action equivalence.

minor comments (8)
  1. [Abstract and Sec. 2] The phrase 'only perturbatively worse for perturbative deviations' is stronger than anything demonstrated in this manuscript; since all derivations are deliberately delegated to Ref. [2], please either state the precise perturbation result with its assumptions or soften the wording to 'as shown in Ref. [2]'.
  2. [Figs. 1 and 2] The dashed 'prediction' lines are normalized to the measured autocorrelation time at T=π/2, so they test only the functional scaling, not the absolute value of the autocorrelation time; please label them as scaling predictions rather than parameter-free predictions.
  3. [Figs. 1–3] The integrated autocorrelation times are reported without error bars, and the integrator step size h is not specified, so the statistical significance of the differences between EFA and no-FA, and of the agreement with the dashed lines, cannot be assessed; please add error bars (for example using the method of Ref. [9]) and report h and acceptance rates precisely.
  4. [Algorithm 3] The algorithm assumes M is positive definite with eigenvalues ω_i^2>0; for clarity, please state this assumption explicitly and indicate what should be done if M has zero or negative modes.
  5. [Algorithm 4] Please clarify that the radial update multiplies the whole vector by e^γ and that the acceptance probability follows from changing variables in R^d; the current description is correct but could be clearer for readers who are not already familiar with radial updates.
  6. [Sec. 2.1] The left panel uses physical parameters with ω0≈6 meV while the right panel uses ω0=1 eV; please add a sentence explaining this difference in units so that the reader does not mistake the two panels for the same physical regime.
  7. [References] Reference [12] lists the No-U-Turn sampler as 'Homan and Gelman'; the correct spelling is 'Hoffman and Gelman'.
  8. [Acknowledgements] The first sentence reads 'The authors thanks'; this should be 'The author thanks'.

Circularity Check

1 steps flagged · score 2.0 of 10

Central HMC-EFA equivalence is self-contained; only a peripheral 'prediction' line in Fig. 2 is normalized to a measured autocorrelation value.

  1. fitted input called prediction [Figure 2 caption, Section 3.1.2]
    "The dashed line shows the prediction for τint(𝑇)=𝜏int(𝑇 = 𝜋/2)· (1−cos(𝑇)/(𝜋/2𝑇))−1 from corollary 2 in Ref. [2]."

    The dashed line is labelled a 'prediction' but its multiplicative factor is τint(T=π/2), which is the integrated autocorrelation measured from the same EFA simulations at trajectory length T=π/2. Using this measured value as the normalization forces the curve to pass through the data point at T=π/2 by construction, so agreement at that point carries no confirmatory weight. The only genuinely predictive content is the shape of the T-dependence; the absolute scale is an input from the data being compared. This is a mild, peripheral instance of a fitted input being presented as a prediction; it does not affect the central equivalence claim, which is demonstrated by Algorithm 3.

full rationale

The paper's central claim — that HMC with EFA is equivalent to direct sampling for quadratic actions and nearly optimal for near-quadratic actions — is supported self-containedly in Algorithms 1-3: with V=0 and T=π/2, the exact harmonic integrator maps p0 ~ N(0,M) to x_f = M^{-1} p0 ~ N(0,M^{-1}) and preserves H, so acceptance is 1 and each trajectory is an exact target sample. This does not reduce to a fit or to a definition. The numerical demonstrations against the SSH model (external first-principles data) and the Ising model (external benchmark) provide independent evidence for the efficiency claims. The only concrete reduction I can exhibit is the dashed 'prediction' line in Fig. 2, whose normalization uses the measured τint(T=π/2); that one point is reproduced by construction, though the T-dependence remains a genuine corollary prediction. This is a minor presentational circularity in a peripheral figure, not a load-bearing defect in the derivation. Accordingly the circularity score is low.

Assumptions & free parameters 5 free parameters · 5 assumptions · 1 invented entities

The guide introduces no new entities. Its recommendations rely on eigendecompositions of the harmonic matrix and on asymptotic properties of the action for radial updates; these are domain assumptions rather than invented objects.

free parameters (5)
  • No-FA trajectory length T (Fig. 3 left) = 0.6
    Tuned by hand to minimize τ_int on a 15x15 lattice, which may bias the EFA comparison.
  • No-FA trajectory length T (Fig. 3 right) = 1/√β
    Chosen for 'correct scaling' rather than optimized; affects the comparison of autocorrelation times.
  • Prediction normalization τ_int(T=π/2) = measured value
    The dashed prediction lines in Figs. 1 and 2 use the measured autocorrelation time at a reference point as normalization, so the 'prediction' is a fitted curve shape.
  • Regulator μ in Eq. (6) = unspecified
    Any μ large enough to make M+μ positive definite is valid; the paper admits the optimal kinetic term cannot be derived without knowledge of V.
  • Radial update standard deviation σ = sqrt(2/(a d))
    Default choice based on asymptotic exponent a of the action; if a is unknown this becomes a hand-specified parameter.
assumptions (5)
  • domain assumption Target distributions have the form P(x) ∝ e^{-S(x)} with continuous variables x and a known action S.
    Stated in Sec. 2; the entire recipe is built on this representation.
  • domain assumption The harmonic matrix M is symmetric positive definite and diagonalizable as M = Ω diag(ω^2) Ω†.
    Algorithm 3 assumes this eigendecomposition is available; for translationally invariant problems Ω is a Fourier transform.
  • standard math The leap-frog integrator is symplectic and, with a Metropolis accept/reject step, samples the target distribution exactly.
    Standard HMC theory, invoked via Refs. [3,4,7]; not proven in this paper.
  • domain assumption For radial updates, the action behaves as S(x) ≈ c |x|^a for large |x|, with known exponent a.
    Algorithm 4 requires this asymptotic form to set the proposal scale σ.
  • domain assumption The anharmonic perturbation V is small enough that the exact integration of the harmonic part leaves a small ΔH.
    Used to justify high acceptance with a single leap-frog step and the 'perturbatively worse' claim; not quantified in this paper.
invented entities (1)
  • None
    purpose: No new particles, forces, dimensions, or conserved quantities are introduced.
    The paper is a recipe combining existing algorithms; it postulates no new physical or mathematical entities.

how reviews work

0 comments
Cite this review

Pith. "Pith review of The Physicist's Guide to the HMC." pith.science (2026). https://pith.science/paper/POSHSVPW

@misc{pith2026250119130,
  author       = {Pith},
  title        = {Pith review of: The Physicist's Guide to the HMC},
  year         = {2026},
  howpublished = {\url{https://pith.science/paper/POSHSVPW}},
  note         = {Machine review of arXiv:2501.19130}
}
abstract

The hybrid Monte Carlo (HMC) algorithm is arguably the most efficient sampling method for general probability distributions of continuous variables. Together with exact Fourier acceleration (EFA) the HMC becomes equivalent to direct sampling for quadratic actions $S(x)=\frac12 x^\mathsf{T} M x$ (i.e. normal distributions $x\sim \mathrm{e}^{-S(x)}$), only perturbatively worse for perturbative deviations of the action from the quadratic case, and it remains viable for arbitrary actions. In this work the most recent improvements of the HMC including EFA and radial updates are collected into a numerical recipe.

Figures

Figures reproduced from arXiv: 2501.19130 by the authors.

Figure 1
Figure 1. Integrated autocorrelation time 𝜏int of the phonon [PITH_FULL_IMAGE:figures/full_fig_p005_1.png] view at source ↗
Figure 2
Figure 2. Integrated autocorrelation time 𝜏int of the absolute value of the magnetisation |𝑚| in the 2D Ising model using HMC simulations [13] on a 15 × 15 lattice. Left: different coupling strengths (𝐽 ≈ 0.44 is the critical coupling [14]); Right: weak coupling 𝐽 = 0.2 and different trajectory lengths 𝑇. All simulations used EFA and have similar acceptance ≳ 80%. The measurement frequency has been adjusted so that the HMC ti… view at source ↗
Figure 3
Figure 3. Integrated autocorrelation time 𝜏int of the plaquette expectation value ⟨𝑃(𝛽)⟩ in pure gauge theory HMC simulations. Left: 2D, U(1) weak coupling 𝛽 = 10 and different lattice sizes 𝐿; Right: 4D, 𝐿 = 10 lattice, SU(3) and different coupling strengths 𝛽. All simulations required the same compute time per trajectory and volume. Measurements were performed every trajectory. Trajectory lengths without FA were chosen as f… view at source ↗

Discussion (0). Continue with ORCID to comment.

Reference graph

Works this paper leans on

22 extracted references · 16 canonical work pages

  1. [2]

    J.OstmeyerandP.Buividovich, MinimalAutocorrelationinHybridMonteCarlosimulations using Exact Fourier Acceleration, 2404.09723

    arXiv

  2. [1]

    Duane, A.D

    S. Duane, A.D. Kennedy, B.J. Pendleton and D. Roweth,Hybrid Monte Carlo,Phys. Lett. B 195 (1987) 216

    1987

  3. [3]

    Omelyan, I

    I. Omelyan, I. Mryglod and R. Folk,Symplectic analytically integrable decomposition algorithms: classification, derivation, and application to molecular dynamics, quantum and celestial mechanics simulations, Computer Physics Communications151 (2003) 272

    work page 2003

  4. [4]

    Optimised Trotter Decompositions for Classical and Quantum Computing

    J. Ostmeyer,Optimised Trotter decompositions for classical and quantum computing,J. Phys. A56(2023) 285303 [2211.02691]

    work page Pith review arXiv 2023

  5. [5]

    Batrouni, G.R

    G.G. Batrouni, G.R. Katz, A.S. Kronfeld, G.P. Lepage, B. Svetitsky and K.G. Wilson, Langevin simulations of lattice field theories, Phys. Rev. D32(1985) 2736

    work page 1985

  6. [6]

    Choice of Integrator in the Hybrid Monte Carlo Algorithm

    T. Takaishi,Choice of integrator in the hybrid Monte Carlo algorithm, Comput. Phys. Commun. 133(2000) 6 [hep-lat/9909134]

    work page Pith review arXiv 2000

  7. [7]

    Neal,Handbook of Markov Chain Monte Carlo, Chapman and Hall/CRC (5, 2011), 10.1201/b10905, [1206.1901]

    R.M. Neal,Handbook of Markov Chain Monte Carlo, Chapman and Hall/CRC (5, 2011), 10.1201/b10905, [1206.1901]

    arXiv 2011

  8. [8]

    First-principle quantum Monte-Carlo study of charge carrier mobility in organic molecular semiconductors

    J. Ostmeyer, T. Nematiaram, A. Troisi and P. Buividovich,First-principles quantum Monte Carlo study of charge-carrier mobility in organic molecular semiconductors,Phys. Rev. Applied 22 (2024) L031004 [2312.14914]

    work page Pith review arXiv 2024

Show all 22 references
  1. [9]

    Wolff,Monte Carlo errors with less errors,Computer Physics Communications156 (2004) 143

    U. Wolff,Monte Carlo errors with less errors,Computer Physics Communications156 (2004) 143

    2004

  2. [10]

    Mackenzie,An Improved Hybrid Monte Carlo Method, Phys

    P.B. Mackenzie,An Improved Hybrid Monte Carlo Method, Phys. Lett. B226(1989) 369

    1989

  3. [11]

    Apers, S

    S. Apers, S. Gribling and D. Szilágyi,Hamiltonian Monte Carlo for efficient Gaussian sampling: long and random steps, 2022. 8 The Physicist’s Guide to the HMC Johann Ostmeyer

    2022

  4. [12]

    Homan and A

    M.D. Homan and A. Gelman,The No-U-turn sampler: adaptively setting path lengths in Hamiltonian Monte Carlo,J. Mach. Learn. Res.15(2014) 1593 [1111.4246]

    2014 arXiv

  5. [13]

    Ostmeyer, E

    J. Ostmeyer, E. Berkowitz, T. Luu, M. Petschlies and F. Pittler,The Ising Model with Hybrid Monte Carlo, Comput. Phys. Commun.265 (2021) 107978

    2021

  6. [14]

    Onsager,Crystal Statistics

    L. Onsager,Crystal Statistics. I. A Two-Dimensional Model with an Order-Disorder Transition, Phys. Rev.65 (1944) 117

    1944

  7. [15]

    Kennedy and X

    A. Kennedy and X. Yu,The convergence of HMC on non-compact Riemannian Manifolds, forthcoming

  8. [16]

    Ostmeyer,Exponential speed up in Monte Carlo sampling through Radial Updates, 2411.18218

    J. Ostmeyer,Exponential speed up in Monte Carlo sampling through Radial Updates, 2411.18218

    arXiv

  9. [17]

    Temmen, E

    F. Temmen, E. Berkowitz, A. Kennedy, T. Luu, J. Ostmeyer and X. Yu,Overcoming Ergodicity Problems of the Hybrid Monte Carlo Method using Radial Updates, in41st International Symposium on Lattice Field Theory, 10, 2024 [2410.19148]

    2024 arXiv

  10. [18]

    Temmen, E

    F. Temmen, E. Berkowitz, A. Kennedy, T. Luu, J. Ostmeyer and X. Yu,Fully ergodic simulations of the Hubbard model using Radial Updates, forthcoming

  11. [19]

    Wilson,Confinement of Quarks,Phys

    K.G. Wilson,Confinement of Quarks,Phys. Rev. D10 (1974) 2445

    1974

  12. [20]

    Gattringer and C.B

    C. Gattringer and C.B. Lang,Quantum chromodynamics on the lattice, vol. 788, Springer, Berlin (2010), 10.1007/978-3-642-01850-3

    2010 doi

  13. [21]

    Duane and B.J

    S. Duane and B.J. Pendleton,Gauge invariant fourier acceleration,Physics Letters B206 (1988) 101

    1988

  14. [22]

    Ostmeyer,comp-avg: Compare Averages of time series and more, 10.5281/zenodo.10794620 (2024)

    J. Ostmeyer,comp-avg: Compare Averages of time series and more, 10.5281/zenodo.10794620 (2024). 9

    2024 doi

Pith tools

Reviewed August 9, 2026 · model on record in the stance chip above.