Pith. sign in

REVIEW 2 major objections 5 minor 53 references

Unadjusted HMC and underdamped Langevin have near-dimension-free marginal bias under weak or sparse interactions.

Reviewed by Pith at T0; open to challenge. T0 means a machine referee read the full paper against a public rubric. the ladder, T0–T4 →

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

Unadjusted HMC and BAOAB Langevin exhibit delocalization of bias: W2,ℓ∞ bias scales as O(h√log d) under weak/sparse interactions, so O(√K) integration steps control K-marginal bias.

T0 review reviewed 2026-08-01 challenge →

load-bearing objection Solid extension of the delocalization-of-bias framework to leapfrog HMC and BAOAB, with a genuinely useful matrix-polynomial toolkit and credible proofs; the abstract's finite-time 'O(√K) steps' claim is plausible but not actually proved. the 2 major comments →

arxiv 2607.15208 v1 pith:QP6NCGDQ submitted 2026-07-16 stat.CO cs.LGmath.PRstat.ML

Delocalization of bias in unadjusted Hamiltonian Monte Carlo and underdamped Langevin

Yifan Chen , Xiaoou Cheng , Jonathan Niles-Weed , Jonathan Weare This is my paper
classification stat.CO cs.LGmath.PRstat.ML MSC 65C0560J2265C40
keywords delocalization of biasunadjusted Hamiltonian Monte Carlounderdamped LangevinBAOAB integratorWasserstein distancemarginal biasmatrix polynomial propagatorsstrongly log-concave sampling
verification ladder T0 review T1 audit T2 compute T3 formal T4 reserved

The pith

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

The reading

Unadjusted Hamiltonian Monte Carlo (HMC) and the BAOAB discretization of underdamped Langevin are biased at any finite step size, but this paper proves the bias is delocalized: the W2,ℓ∞ distance between the chain's invariant measure and the target grows with the step size h times only a logarithm of the full dimension d, provided the target's interactions are weak or sparse. More concretely, any K-dimensional marginal has W2 bias at most O(√K h √log d) up to condition-number factors, so O(√K) integration steps suffice to control the marginal bias. This matters because the standard fix—Metropolis–Hastings adjustment—forces the step size to shrink with d to keep acceptance rates positive, while these results say unadjusted samplers can avoid that bottleneck when only low-dimensional marginals are needed. The proof handles the position–momentum coupling of phase-space integrators through a new matrix-polynomial representation of their propagators.

Core claim

The central claim is that the delocalization-of-bias phenomenon, previously established for overdamped Langevin, extends to two widely used phase-space samplers. Under strong convexity and smoothness plus weak interactions, unadjusted HMC with total integration time mh=1/√(20β) and BAOAB with h≤(1−η)/(2√β) have invariant measures satisfying W2,ℓ∞(π_h,π)=O((β/α) h √log(2d)). Under sparse interactions with polynomial neighborhood growth, the bound is h√log(2d) [O((β/α) log(2d))]^{n/2+1}. Since any K-marginal W2 error is at most √K W2,ℓ∞, both bounds imply that O(√K) integration steps control the bias of every K-dimensional marginal up to log d factors. The BAOAB result holds for all sufficient

What carries the argument

The proof is carried by two tools. The first is the W2,ℓ∞ metric, a Wasserstein distance built from the ℓ∞ norm, which upper-bounds the W2 error of every low-dimensional marginal without paying for the full dimension. The second is a family of multivariate damped Chebyshev matrix polynomials—defined by a three-term recurrence in ordered, non-commuting matrices—that exactly represent the multi-step propagators of the leap-frog and BAOAB integrators on phase space. These polynomials turn the intractable position–momentum mixing into coefficient bounds: ℓ2 contraction estimates are converted into ℓ∞ estimates using either near-diagonality (weak interactions) or neighborhood sparsity (sparse int

Load-bearing premise

The load-bearing premise is that the interaction structure is benign in a pointwise sense—every off-diagonal Hessian entry bounded uniformly by C^(O)α at every point (weak interactions), or the interaction graph has at most polynomial neighborhood growth (sparse interactions); if a single pair of variables interacts more strongly at any point, the ℓ∞ propagator contraction that the proof needs is no longer guaranteed.

What would settle it

Simulate unadjusted HMC on V = (α/2)‖q‖² + c q₁q₂ with c set just above the threshold (e.g., c = α/10 for HMC) and d large; measure the W2 distance of a one-dimensional marginal of the invariant measure. If the delocalized bound O((β/α)h√log d) still holds, the uniform off-diagonal condition is not necessary; if the error grows with d, the condition is load-bearing. The same test with c just below the threshold should confirm the predicted plateau.

Watch this falsifier. Get emailed when new claim-graph text bears on it.

If this is right

  • For weakly interacting targets, the invariant-measure bias of unadjusted HMC and BAOAB is O((β/α) h √log d), independent of d except for a log factor.
  • For sparse interaction graphs with polynomial neighborhood growth, the same bias is h√log d times a condition-number-log factor, so the marginal bias stays nearly dimension-free even when interactions are not globally weak.
  • Any K-dimensional marginal has W2 bias at most √K times the W2,ℓ∞ bound, so O(√K) integration steps give controlled marginal bias up to log d factors.
  • The BAOAB contraction is valid for all large friction values, so the delocalization effect carries over to the overdamped-limit integrator commonly used for molecular sampling.
  • These bounds give a theoretical rationale for running unadjusted HMC and underdamped Langevin without Metropolis adjustment when the estimand is a low-dimensional marginal.

Where Pith is reading between the lines

These are editorial extensions of the paper, not claims the author makes directly.

  • The uniform off-diagonal Hessian condition is stronger than a total-variation or row-sum notion of weak interaction; a natural testable relaxation is an average or ℓ2 row-sum bound on the off-diagonal part, which might preserve the delocalized bias for a wider class of dense graphs.
  • The matrix-polynomial propagator representation is scheme-agnostic, so the same coefficient machinery should yield delocalization bounds for other splittings (e.g., OBABO/ABOBA) and could extend to randomized integration times if the randomness is handled as an additional average over the polynomial degree.
  • A practical consequence the paper does not spell out: in Bayesian inverse problems where only marginals are reported, the step size can be chosen from the local curvature scale of the coordinates of interest rather than the global dimension, potentially reducing wall-clock time relative to adjusted HMC.
  • The Gaussian example suggests the delocalization threshold might be considerably sharper than the worst-case weak-interaction condition; an exact bias formula for quadratic perturbations of product measures would reveal how the off-diagonal strength enters the constant.
Share X Bluesky LinkedIn Reddit HN

Editorial analysis

A structured set of objections, weighed in public.

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

Referee Report

2 major / 5 minor

Summary. The paper extends the 'delocalization of bias' phenomenon from the overdamped Langevin algorithm to two unadjusted kinetic samplers: leapfrog-discretized Hamiltonian Monte Carlo (HMC) and the BAOAB discretization of underdamped Langevin (UL). Under a strong-log-concavity/log-smoothness assumption, the authors prove stationary-bias bounds in the W_{2,ℓ∞} metric. For weakly interacting targets (Assumption 1.3), Theorem 1.4 gives W_{2,ℓ∞}(π_h,π)=O((β/α) h √log(2d)), with step-size restrictions depending on the algorithm. For sparsely interacting targets (Assumption 1.5), Theorem 1.6 gives W_{2,ℓ∞}(π_h,π)=h √log(2d) [O((β/α) log(2d))]^{n/2+1}. The analysis introduces a matrix-polynomial representation of leapfrog and BAOAB propagators, with ℓ∞ norm bounds obtained from ℓ2 contraction plus weak or sparse structure. The abstract goes further, claiming that O(√K) integration steps suffice to control the W2 bias of any K-dimensional marginal.

Significance. If the main theorems are correct, this is a valuable extension of [17] to samplers more commonly used in high-dimensional problems. The matrix-polynomial framework is a genuine technical contribution and appears likely to be reusable for other splitting integrators. The paper is careful to state structural assumptions and to provide explicit propagation and discretization-error estimates in the appendices, with no fitted parameters entering the bounds. The principal caveat is that the advertised 'O(√K) integration steps suffice' result is not actually proved; the theorems establish only bounds on the invariant law π_h. In addition, one displayed discretization-error estimate in the UL argument has a scaling that appears numerically inconsistent with the claimed O(1) exponential factor. These issues are local to the presentation and proof, but they are load-bearing for the paper's central claims.

major comments (2)
  1. [Abstract; §6; Theorems 1.4 and 1.6] The abstract and Section 6 claim that O(√K) integration steps suffice to control the W2 bias of any K-dimensional marginal. The theorems proved are stationary-bias results: they bound W_{2,ℓ∞}(π_h,π), not the error of a finite-time chain. Combining the K-marginal projection bound W_2(π_h^{(K)},π^{(K)})≤√K W_{2,ℓ∞}(π_h,π) with h∝1/√K gives √K leapfrog steps per outer HMC loop, but no theorem controls W_2(ρ_T^{(K)},π^{(K)}) after T steps from a specified initial law. The ingredients for such a theorem are present (contraction in Proposition 2.1, discretization error in §5), and a finite-time corollary is plausible, but it is not stated. As written, the paper does not prove the advertised iteration-complexity statement. Please add a finite-time corollary (with explicit dependence on initial W_{2,ℓ∞} and ε) or revise the abstract and Section 6 to claim only 'bias of the stationary law.'
  2. [§5.2, Eq. (5.5); used in §F.2, Eq. (F.22)] The exponential in Proposition 5.2 is exp((l−1)h( (3/2) C_5 √β + (1/16) C_6/β^{3/2} + √β )). Under Assumption 1.1, C_5≍β, and with l=⌈1/(2√β h)⌉ one has (l−1)h ≈1/(2√β). The first term is then approximately (3/4)β, so the exponential is exp(O(β)), not O(1) as claimed in Appendix F. This would destroy the final O((β/α) h √log(2d)) bound in the UL case. Inspecting the recurrence (E.37)–(E.39) suggests the coefficient should be C_5/√β, not C_5√β: the term h|Δ_p| in the position recurrence becomes h√β times the weighted coordinate √a|Δ_p|. If the displayed formula is a typo, it must be corrected and the downstream constants re-checked; if the displayed formula is intended, the proof of the UL weak-interaction theorem fails.
minor comments (5)
  1. [Assumption 1.5; Theorem 1.6] The constants C and n in s_k ≤ C(k+1)^n should be explicitly required to be independent of d (or the theorem should allow them to depend on d and state the consequence). As written, the exponent n/2+1 in Theorem 1.6 could depend on d, which would undermine the dimension-free interpretation.
  2. [Propositions 2.2 and 4.5] The definition of r_i contains a stray artifact: '⌈ih√βe+ log√d/log(5/3)⌉' should presumably be '⌈i h √β + log√d/log(5/3)⌉'. Please check the original LaTeX and correct the displayed formulas.
  3. [Lemma C.2] The sentence 'The (i,j) entry of the product is nonzero only if j∈N_r(i), or if i∈N_r(k) and j∈N_r(k) for some k' is garbled and does not match the concluding statement. The intended statement is that the (i,j) entry is supported on j∈N_{2r}(i), as used below. Please rewrite.
  4. [Proposition 2.1 and Theorem 1.4] For HMC, m is required to satisfy mh=1/√(20β), but m is not explicitly declared to be an integer and no discussion is given of rounding when h does not divide this time. A sentence explaining the integer-rounding convention would be helpful.
  5. [Example 1.2] The proof of Example 1.2 uses the formula for π_h from [5], but it is stated only for 'when HMC converges' and explicitly excludes resonant values of mh. This is acceptable, but the set of excluded h should be described more precisely, even if measure-zero.

Circularity Check

0 steps flagged

No significant circularity; the bias bounds are derived from structural assumptions via explicit propagator and discretization estimates.

full rationale

The paper's central results (Theorems 1.4 and 1.6) are proved from Assumptions 1.1, 1.3, and 1.5 through explicit iterative inequalities (Propositions 2.1 and 2.2), ℓ∞ propagator bounds (Propositions 4.3–4.6), and discretization-error estimates (Propositions 5.1, 5.2, and Appendices E–F). The constants C^(O) = 1/50 and 1/20 are structural assumption constants, not fitted parameters, and no prediction is defined in terms of the target quantity. The self-citations to [17] supply the W_{2,ℓ∞} metric, the coupling framework, and two elementary lemmas (Lemma E.1 and Lemma E.2) used to obtain log(2d) factors; these are auxiliary, externally checkable facts and do not assume the target HMC/UL bias result. Citations to [34] and [11] are independent prior contraction results, not self-citations. The only notable issue is that the abstract's phrase 'O(√K) integration steps suffice' is stated as though it were a finite-time complexity guarantee, whereas the proved results bound the stationary bias; this is a correctness/overreach concern rather than circularity, and it does not involve fitting or definitional equivalence. No circular step can be exhibited by reducing an equation of the paper to its own input.

Axiom & Free-Parameter Ledger

5 free parameters · 5 axioms · 1 invented entities

No numerical fitting is involved; all listed constants are explicit structural or step-size parameters. The central claim depends on the assumed forms of V and on smallness constants C^(O), not on data. The only newly introduced object is a proof device (the matrix polynomials), which has no independent empirical content.

free parameters (5)
  • C^(O) (weak-interaction strength) = 1/50 (HMC), 1/20 (UL)
    Ad hoc small constant chosen so the ℓ∞ contraction inequalities in Propositions 4.3 and 4.4 hold; not fit to data.
  • HMC integration time per outer loop mh = 1/√(20β)
    Hand-chosen to optimize the contraction rate; appears as a step-size condition in Theorem 1.4.
  • UL step-size bound = h ≤ (1−η)/(2√β)
    Needed for weight-matrix equivalence (2.10) and BAOAB contraction; inherited from the prior BAOAB analysis [34].
  • Sparsity thresholds r_i = ⌈ih√βe + log√d/log(5/3)⌉ (HMC); ⌈ie²βh²/(1−η) + log√d⌉ (UL)
    Hand-chosen thresholds balancing sparse-neighborhood support against geometric decay in the ℓ∞ propagator bounds.
  • Contraction block lengths N and K_1 = N = ⌈(200β/α) log(2√d)⌉; K_1 = ⌈(16√β(1−η) log(14√2 d))/(hα) + 1⌉
    Chosen so the contraction factor is at most 1/2; pure proof constants, not fitted to data.
axioms (5)
  • domain assumption Assumption 1.1: V ∈ C²(R^d), α-strongly convex and β-smooth with 0 < α ≤ β < ∞.
    Underlies all contraction and discretization-error estimates; stated in Section 1.3.
  • domain assumption Assumption 1.3: the off-diagonal Hessian part has |∇²V^(O)|_{ℓ∞} ≤ C^(O)α uniformly.
    Weak-interaction structure used to preserve ℓ∞ contraction of the propagators.
  • domain assumption Assumption 1.5: V(x) = Σ_i V_i(X_i) with each V_i depending only on a local neighborhood, and the sparsity parameters s_k grow at most polynomially.
    Sparse-interaction structure used to control locality of Hessian products in Lemma C.2.
  • domain assumption The coupled exact chain Y starts at stationarity π⊗N(0,I) and the discrete chain converges to an invariant distribution π_h.
    The bias is defined as the limit of the coupled error as k→∞; convergence is assumed or inferred from the contraction bounds.
  • standard math Known contraction results [11, Lemma 19] and [34, Theorem 5.1] can be adapted from state-dependent Hessians to arbitrary matrix sequences satisfying αI ⪯ H_i ⪯ βI.
    The paper reproduces the adaptation, but the underlying techniques are prior art and are not machine-checked.
invented entities (1)
  • Multivariate damped Chebyshev polynomials no independent evidence
    purpose: Represent the multi-step position-difference propagators of leapfrog HMC and BAOAB UL as ordered, non-commuting matrix polynomials.
    A new mathematical device introduced in this paper; it is internally defined and used for proofs, with no falsifiable handle outside the paper.

reviewed 2026-08-01 · how reviews work

0 comments
Cite this review

Pith. "Pith review of Delocalization of bias in unadjusted Hamiltonian Monte Carlo and underdamped Langevin." pith.science (2026). https://pith.science/paper/QP6NCGDQ

@misc{pith2026260715208,
  author       = {Pith},
  title        = {Pith review of: Delocalization of bias in unadjusted Hamiltonian Monte Carlo and underdamped Langevin},
  year         = {2026},
  howpublished = {\url{https://pith.science/paper/QP6NCGDQ}},
  note         = {Machine review of arXiv:2607.15208}
}
Share X Bluesky LinkedIn Reddit HN
abstract

Unadjusted samplers such as unadjusted Hamiltonian Monte Carlo and underdamped Langevin are well-known to be biased. Metropolis--Hastings adjustment has been conventionally incorporated into Hamiltonian Monte Carlo to eliminate the bias. However, this adjustment can significantly increase the iteration complexity due to the small step size required for reasonable Metropolis acceptance rates. In this work, we extend the \emph{delocalization of bias} phenomenon, previously established for the overdamped Langevin algorithm, to these two unadjusted algorithms. We show that to control the $W_2$ bias of any $K$-dimensional marginal of a high-dimensional distribution, $O(\sqrt{K})$ integration steps suffice up to $\log d$ terms, assuming either weak or sparse interactions among variables. The discrete-time integrators here introduce technical difficulties beyond those of the overdamped setting, which we address through a broadly applicable matrix-polynomial framework that characterizes their propagators. Our result for the underdamped Langevin algorithm is valid for all large friction parameters, implying that the Leimkuhler-Matthews integrator for the overdamped Langevin dynamics also exhibits delocalization of bias.

Figures

Figures reproduced from arXiv: 2607.15208 by Jonathan Niles-Weed, Jonathan Weare, Xiaoou Cheng, Yifan Chen.

Figure 1
Figure 1. Figure 1: Illustration of a potential V (x) with sparse interactions. Adapted from [17]. Compared to our results for the unadjusted Langevin algorithm in [17], we have an improvement from p h log(2d) to h p log(2d). The new scaling aligns with [34], where the Wasserstein-2 bias of the BAOAB scheme is h √ d. It does not show a second order accuracy due to lack of higher order regularity conditions on V [39, 8]. We al… view at source ↗

discussion (0)

Sign in with ORCID, Apple, or X to comment. Anyone can read and Pith papers without signing in.

Reference graph

Works this paper leans on

53 extracted references · 5 linked inside Pith

  1. [1]

    GROMACS: High performance molecular simulations through multi-level paral- lelism from laptops to supercomputers.SoftwareX, 1:19–25, 2015

    Mark James Abraham, Teemu Murtola, Roland Schulz, Szil´ ard P´ all, Jeremy C Smith, Berk Hess, and Erik Lindahl. GROMACS: High performance molecular simulations through multi-level paral- lelism from laptops to supercomputers.SoftwareX, 1:19–25, 2015

    2015

  2. [2]

    Stegun.Handbook of Mathematical Functions: With Formulas, Graphs, and Mathematical Tables

    Milton Abramowitz and Irene A. Stegun.Handbook of Mathematical Functions: With Formulas, Graphs, and Mathematical Tables. Dover, 1964

    1964

  3. [3]

    Allen and Dominic J

    Michael P. Allen and Dominic J. Tildesley.Computer Simulation of Liquids. Oxford University Press, 2017

    2017

  4. [4]

    Molecular dynamics simulations at constant pressure and/or temperature.The Journal of chemical physics, 72(4):2384–2393, 1980

    Hans C Andersen. Molecular dynamics simulations at constant pressure and/or temperature.The Journal of chemical physics, 72(4):2384–2393, 1980

    1980

  5. [5]

    Hamiltonian Monte Carlo for efficient Gaussian sampling: Long and random steps.Journal of Machine Learning Research, 25(348):1–30, 2024

    Simon Apers, Sander Gribling, and D´ aniel Szil´ agyi. Hamiltonian Monte Carlo for efficient Gaussian sampling: Long and random steps.Journal of Machine Learning Research, 25(348):1–30, 2024

    2024

  6. [6]

    Optimal tuning of the Hybrid Monte-Carlo Algorithm.Bernoulli, 19(5A):1501–1534, 2013

    Alexandros Beskos, Natesh Pillai, Gareth Roberts, Jesus-Maria Sanz-Serna, and Andrew Stuart. Optimal tuning of the Hybrid Monte-Carlo Algorithm.Bernoulli, 19(5A):1501–1534, 2013

    2013

  7. [7]

    Mixing time guarantees for unadjusted Hamiltonian Monte Carlo.Bernoulli, 29(1):75–104, 2023

    Nawaf Bou-Rabee and Andreas Eberle. Mixing time guarantees for unadjusted Hamiltonian Monte Carlo.Bernoulli, 29(1):75–104, 2023

    2023

  8. [8]

    Unadjusted Hamiltonian MCMC with stratified Monte Carlo time integration.The Annals of Applied Probability, 35(1):360–392, 2025

    Nawaf Bou-Rabee and Milo Marsden. Unadjusted Hamiltonian MCMC with stratified Monte Carlo time integration.The Annals of Applied Probability, 35(1):360–392, 2025

    2025

  9. [9]

    Tail-sensitive KL and R´ enyi convergence of unadjusted Hamiltonian Monte Carlo via one-shot couplings.arXiv preprint arXiv:2601.09019, 2026

    Nawaf Bou-Rabee, Siddharth Mitra, and Andre Wibisono. Tail-sensitive KL and R´ enyi convergence of unadjusted Hamiltonian Monte Carlo via one-shot couplings.arXiv preprint arXiv:2601.09019, 2026

    arXiv 2026

  10. [10]

    Randomized Hamiltonian Monte Carlo.The Annals of Applied Probability, 27(4):2159–2194, 2017

    Nawaf Bou-Rabee and Jes´ us Mar ´ ıa Sanz-Serna. Randomized Hamiltonian Monte Carlo.The Annals of Applied Probability, 27(4):2159–2194, 2017

    2017

  11. [11]

    Convergence of unadjusted Hamiltonian Monte Carlo for mean-field models.Electronic Journal of Probability, 28:1–40, 2023

    Nawaf Bou-Rabee and Katharina Schuh. Convergence of unadjusted Hamiltonian Monte Carlo for mean-field models.Electronic Journal of Probability, 28:1–40, 2023

    2023

  12. [12]

    Second order quantitative bounds for unadjusted generalized Hamiltonian Monte Carlo.arXiv preprint arXiv:2306.09513, 2023

    Evan Camrud, Alain Durmus, Pierre Monmarch´ e, and Gabriel Stoltz. Second order quantitative bounds for unadjusted generalized Hamiltonian Monte Carlo.arXiv preprint arXiv:2306.09513, 2023

    Pith/arXiv arXiv 2023

  13. [13]

    Hoffman, Daniel Lee, Ben Goodrich, Michael Be- tancourt, Marcus A

    Bob Carpenter, Andrew Gelman, Matthew D. Hoffman, Daniel Lee, Ben Goodrich, Michael Be- tancourt, Marcus A. Brubaker, Jiqiang Guo, Peter Li, and Allen Riddell. Stan: A probabilistic programming language.Journal of Statistical Software, 76(1):1–32, 2017

    2017

  14. [14]

    Recent developments in Amber biomolecular simulations.Journal of Chemical Information and Modeling, 65(15):7835–7843, 2025

    David A Case, David S Cerutti, Vin ´ ıcius Wilian D Cruzeiro, Thomas A Darden, Robert E Duke, Mahdieh Ghazimirsaeed, George M Giambasu, Timothy J Giese, Andreas W Gotz, Julie A Harris, et al. Recent developments in Amber biomolecular simulations.Journal of Chemical Information and Modeling, 65(15):7835–7843, 2025

    2025

  15. [15]

    Chada, Benedict Leimkuhler, Daniel Paulin, and Peter A

    Neil K. Chada, Benedict Leimkuhler, Daniel Paulin, and Peter A. Whalley. Unbiased kinetic Langevin Monte Carlo with inexact gradients.arXiv preprint arXiv:2311.05025, 2023

    arXiv 2023

  16. [16]

    Reflection coupling for unadjusted generalized Hamiltonian Monte Carlo in the nonconvex stochastic gradient case.IMA Journal of Numerical Analysis, page draf045, 2025

    Martin Chak and Pierre Monmarch´ e. Reflection coupling for unadjusted generalized Hamiltonian Monte Carlo in the nonconvex stochastic gradient case.IMA Journal of Numerical Analysis, page draf045, 2025. 22 DELOCALIZATION OF BIAS IN UNADJUSTED HMC AND UNDERDAMPED LANGEVIN

    2025

  17. [17]

    Convergence of unad- justed langevin in high dimensions: Delocalization of bias.Communications on Pure and Applied Mathematics, 79(6):1467–1491, 2026

    Yifan Chen, Xiaoou Cheng, Jonathan Niles-Weed, and Jonathan Weare. Convergence of unad- justed langevin in high dimensions: Delocalization of bias.Communications on Pure and Applied Mathematics, 79(6):1467–1491, 2026

    2026

  18. [18]

    When does Metropolized Hamiltonian Monte Carlo provably outperform Metropolis-adjusted Langevin algorithm?arXiv preprint arXiv:2304.04724, 2023

    Yuansi Chen and Khashayar Gatmiry. When does Metropolized Hamiltonian Monte Carlo provably outperform Metropolis-adjusted Langevin algorithm?arXiv preprint arXiv:2304.04724, 2023

    arXiv 2023

  19. [19]

    Zongchen Chen and Santosh S. Vempala. Optimal convergence rate of Hamiltonian Monte Carlo for strongly logconcave distributions.Theory of Computing, 18(9):1–18, 2022

    2022

  20. [20]

    Chatterji, Peter L

    Xiang Cheng, Niladri S. Chatterji, Peter L. Bartlett, and Michael I. Jordan. Underdamped Langevin MCMC: A non-asymptotic analysis. InProceedings of the 31st Conference on Learning Theory (COLT), 2018

    2018

  21. [21]

    Tiangang Cui, Shuigen Liu, and Xin T. Tong. Stein’s method for marginals on large graphical models.arXiv preprint arXiv:2410.11771, 2024

    arXiv 2024

  22. [22]

    Dalalyan and Lionel Riou-Durand

    Arnak S. Dalalyan and Lionel Riou-Durand. On sampling from a log-concave density using kinetic Langevin diffusions.Bernoulli, 26(3):1956 – 1988, 2020

    1956

  23. [23]

    Kennedy, Brian J

    Simon Duane, Anthony D. Kennedy, Brian J. Pendleton, and Duncan Roweth. Hybrid Monte Carlo. Physics letters B, 195(2):216–222, 1987

    1987

  24. [24]

    Asymptotic bias of inexact Markov chain Monte Carlo methods in high dimension.The Annals of Applied Probability, 34(4):3435–3468, 2024

    Alain Durmus and Andreas Eberle. Asymptotic bias of inexact Markov chain Monte Carlo methods in high dimension.The Annals of Applied Probability, 34(4):3435–3468, 2024

    2024

  25. [25]

    Chodera, Robert T

    Peter Eastman, Jason Swails, John D. Chodera, Robert T. McGibbon, Yutong Zhao, Kyle A. Beauchamp, Lee-Ping Wang, Andrew C. Simmonett, Matthew P. Harrigan, Chaya D. Stern, Rafal P. Wiewiora, Bernard R. Brooks, and Vijay S. Pande. OpenMM 7: Rapid development of high performance algorithms for molecular dynamics.PLOS Computational Biology, 13(7):e1005659, 2017

    2017

  26. [26]

    Couplings and quantitative contraction rates for Langevin dynamics.The Annals of Probability, 47(4):1982–2010, 2019

    Andreas Eberle, Arnaud Guillin, and Raphael Zimmer. Couplings and quantitative contraction rates for Langevin dynamics.The Annals of Probability, 47(4):1982–2010, 2019

    1982

  27. [27]

    Elsevier, 2023

    Daan Frenkel and Berend Smit.Understanding Molecular Simulation: From Algorithms to Appli- cations. Elsevier, 2023

    2023

  28. [28]

    HMC and underdamped Langevin united in the unadjusted convex smooth case.SIAM/ASA Journal on Uncertainty Quan- tification, 13(1):278–303, 2025

    Nicola ¨ ı Gouraud, Pierre Le Bris, Adrien Majka, and Pierre Monmarch´ e. HMC and underdamped Langevin united in the unadjusted convex smooth case.SIAM/ASA Journal on Uncertainty Quan- tification, 13(1):278–303, 2025

    2025

  29. [29]

    GROMACS stochastic dynamics and BAOAB are equivalent configurational sampling algorithms.Journal of Chemical Theory and Computation, 18(10):5792–5798, 2022

    Stefanie Kieninger and Bettina G Keller. GROMACS stochastic dynamics and BAOAB are equivalent configurational sampling algorithms.Journal of Chemical Theory and Computation, 18(10):5792–5798, 2022

    2022

  30. [30]

    A hierarchical entropy method for the delocalization of bias in high-dimensional Langevin Monte Carlo.arXiv preprint arXiv:2509.08619, 2025

    Daniel Lacker and Fuzhong Zhou. A hierarchical entropy method for the delocalization of bias in high-dimensional Langevin Monte Carlo.arXiv preprint arXiv:2509.08619, 2025

    Pith/arXiv arXiv 2025

  31. [31]

    Springer, 2015

    Ben Leimkuhler and Charles Matthews.Molecular Dynamics: With Deterministic and Stochastic Numerical Methods. Springer, 2015

    2015

  32. [32]

    Rational construction of stochastic numerical methods for molecular sampling.Applied Mathematics Research eXpress, 2013(1):34–56, 2013

    Benedict Leimkuhler and Charles Matthews. Rational construction of stochastic numerical methods for molecular sampling.Applied Mathematics Research eXpress, 2013(1):34–56, 2013

    2013

  33. [33]

    Robust and efficient configurational molecular sam- pling via Langevin dynamics.The Journal of Chemical Physics, 138(17):174102, 2013

    Benedict Leimkuhler and Charles Matthews. Robust and efficient configurational molecular sam- pling via Langevin dynamics.The Journal of Chemical Physics, 138(17):174102, 2013

    2013

  34. [34]

    Leimkuhler, Daniel Paulin, and Peter A

    Benedict J. Leimkuhler, Daniel Paulin, and Peter A. Whalley. Contraction and convergence rates for discretized kinetic Langevin dynamics.SIAM Journal on Numerical Analysis, 62(3):1226–1258, 2024

    2024

  35. [35]

    Yi-An Ma, Niladri Chatterji, Xiang Cheng, Nicolas Flammarion, Peter Bartlett, and Michael I. Jordan. Is there an analog of Nesterov acceleration for MCMC?Bernoulli, 27(3):1942–1992, 2021

    1942

  36. [36]

    Mixing of Hamiltonian Monte Carlo on strongly log-concave distributions 2: Numerical integrators

    Oren Mangoubi and Aaron Smith. Mixing of Hamiltonian Monte Carlo on strongly log-concave distributions 2: Numerical integrators. InProceedings of the 22nd International Conference on Artificial Intelligence and Statistics (AISTATS), 2019

    2019

  37. [37]

    Mixing of Hamiltonian Monte Carlo on strongly log-concave distributions: Continuous dynamics.The Annals of Applied Probability, 31(5):2019–2045, 2021

    Oren Mangoubi and Aaron Smith. Mixing of Hamiltonian Monte Carlo on strongly log-concave distributions: Continuous dynamics.The Annals of Applied Probability, 31(5):2019–2045, 2021

    2019

  38. [38]

    Dimensionally tight bounds for second-order Hamiltonian Monte Carlo

    Oren Mangoubi and Nisheeth Vishnoi. Dimensionally tight bounds for second-order Hamiltonian Monte Carlo. InAdvances in Neural Information Processing Systems 31 (NeurIPS), 2018. DELOCALIZATION OF BIAS IN UNADJUSTED HMC AND UNDERDAMPED LANGEVIN 23

    2018

  39. [39]

    High-dimensional MCMC with a standard splitting scheme for the under- damped Langevin diffusion.Electronic Journal of Statistics, 15(2):4117–4166, 2021

    Pierre Monmarch´ e. High-dimensional MCMC with a standard splitting scheme for the under- damped Langevin diffusion.Electronic Journal of Statistics, 15(2):4117–4166, 2021

    2021

  40. [40]

    Almost sure contraction for diffusions onR d

    Pierre Monmarch´ e. Almost sure contraction for diffusions onR d. Application to generalised Langevin diffusions.Stochastic Processes and their Applications, 161:316–349, 2023

    2023

  41. [41]

    An entropic approach for Hamiltonian Monte Carlo: The idealized case.The Annals of Applied Probability, 34(2):2243–2293, 2024

    Pierre Monmarch´ e. An entropic approach for Hamiltonian Monte Carlo: The idealized case.The Annals of Applied Probability, 34(2):2243–2293, 2024

    2024

  42. [42]

    Radford M. Neal. MCMC using Hamiltonian dynamics. InHandbook of Markov Chain Monte Carlo, pages 113–162. Chapman and Hall/CRC, 2011

    2011

  43. [43]

    Scalable molecular dynamics on CPU and GPU architectures with NAMD.The Journal of Chemical Physics, 153(4):044130, 2020

    James C Phillips, David J Hardy, Julio DC Maia, John E Stone, Jo˜ ao V Ribeiro, Rafael C Bernardi, Ronak Buch, Giacomo Fiorin, J´ erˆ ome H´ enin, Wei Jiang, et al. Scalable molecular dynamics on CPU and GPU architectures with NAMD.The Journal of Chemical Physics, 153(4):044130, 2020

    2020

  44. [44]

    Tinker 8: Software tools for molecular design.Journal of Chemical Theory and Computation, 14(10):5273–5289, 2018

    Joshua A Rackers, Zhi Wang, Chao Lu, Marie L Laury, Louis Lagard` ere, Michael J Schnieders, Jean-Philip Piquemal, Pengyu Ren, and Jay W Ponder. Tinker 8: Software tools for molecular design.Journal of Chemical Theory and Computation, 14(10):5273–5289, 2018

    2018

  45. [45]

    Metropolis adjusted Langevin trajectories: a robust alter- native to Hamiltonian Monte Carlo.arXiv preprint arXiv:2202.13230, 2022

    Lionel Riou-Durand and Jure Vogrinc. Metropolis adjusted Langevin trajectories: a robust alter- native to Hamiltonian Monte Carlo.arXiv preprint arXiv:2202.13230, 2022

    Pith/arXiv arXiv 2022

  46. [46]

    Black-box unadjusted Hamiltonian Monte Carlo.arXiv preprint arXiv:2412.08876, 2024

    Jakob Robnik, Reuben Cohn-Gordon, and Uroˇ s Seljak. Black-box unadjusted Hamiltonian Monte Carlo.arXiv preprint arXiv:2412.08876, 2024

    Pith/arXiv arXiv 2024

  47. [47]

    Wiecki, and Christopher Fonnesbeck

    John Salvatier, Thomas V. Wiecki, and Christopher Fonnesbeck. Probabilistic programming in Python using PyMC3.PeerJ Computer Science, 2:e55, 2016

    2016

  48. [48]

    Katharina Schuh and Peter A. Whalley. Convergence of kinetic Langevin samplers for non-convex potentials.arXiv preprint arXiv:2405.09992, 2024

    Pith/arXiv arXiv 2024

  49. [49]

    The randomized midpoint method for log-concave sampling

    Ruoqi Shen and Yin Tat Lee. The randomized midpoint method for log-concave sampling. In Advances in Neural Information Processing Systems 32 (NeurIPS), 2019

    2019

  50. [50]

    LAMMPS - a flexible simulation tool for particle-based materials modeling at the atomic, meso, and continuum scales.Computer Physics Communications, 271:108171, 2022

    Aidan P Thompson, H Metin Aktulga, Richard Berger, Dan S Bolintineanu, W Michael Brown, Paul S Crozier, Pieter J In’t Veld, Axel Kohlmeyer, Stan G Moore, Trung Dac Nguyen, et al. LAMMPS - a flexible simulation tool for particle-based materials modeling at the atomic, meso, and continuum scales.Computer Physics Communications, 271:108171, 2022

    2022

  51. [51]

    Zhang, Jason M

    Matthew S. Zhang, Jason M. Altschuler, and Sinho Chewi. Algorithmic warm starts for Hamiltonian Monte Carlo.arXiv preprint arXiv:2603.22741, 2026. The appendices contain all technical details. Section A gives the proof for Gaussian distributions. Section B derives the matrix polynomials. Sections C and D provide bounds for the propagators of the leap-frog...

    arXiv 2026

  52. [52]

    and the matrices defined through the intermediate iterates z1 :=U ξ0 bao,h(X0, P0), z ′ 1 :=U ξ0 bao,h(Y0, P′ 0), z k+1 :=U ab,hzk , z ′ k+1 :=U ab,hz′ k , zi :=U ξi−1 abao,h · · ·Uξ1 abao,hz1 , z ′ i :=U ξi−1 abao,h · · ·Uξ1 abao,hz′ 1 ,2≤i≤k , satisfying W(z 1 −z ′

  53. [53]

    By Proposition 4.4 and| I0 W −1|ℓ∞ = 1 + b√ a−b2 , (F.19) (a)≤6 1 + b√ a−b2 (1− 1 2 c(h))k−1|∆Z0|2,ℓ∞w

    =M w bao(H0)W∆Z 0 , W(z k+1 −z ′ k+1) =M w ab(Hk)W(z k −z ′ k), W(z i −z ′ i) =M w abao(Hi−1)W(z i−1 −z ′ i−1),2≤i≤k . By Proposition 4.4 and| I0 W −1|ℓ∞ = 1 + b√ a−b2 , (F.19) (a)≤6 1 + b√ a−b2 (1− 1 2 c(h))k−1|∆Z0|2,ℓ∞w . Discretization error.Write (∆ l q,∆ l p) = (U l baoab,h −U l hoh,h)(q∗, p∗). The matrix representation ofΠ 1Uk−(i+1)˜l baoab,h (U˜l b...

This paper was first reviewed by deepseek-v4-flash on August 1, 2026.