Practical and Scalable Hamiltonian Monte Carlo Without the Metropolis Test
Abstract
Hamiltonian Monte Carlo and underdamped Langevin Monte Carlo are state-of-the-art methods for taking samples from high-dimensional distributions with a differentiable density function. To generate samples, they numerically integrate Hamiltonian or Langevin dynamics. This numerical integration introduces an asymptotic bias in Monte Carlo estimators of expectation values, which can be eliminated by adjusting the dynamics with a Metropolis-Hastings (MH) proposal step. Alternatively, one can trade bias for variance by avoiding MH, and select an integration step size that ensures sufficiently small asymptotic bias, relative to the variance inherent in a finite set of samples. Such \emph{unadjusted} methods often significantly outperform their adjusted counterparts in high-dimensional problems where sampling would otherwise be prohibitively expensive, yet are rarely used in statistical applications due to the absence of an automated way of choosing a step size. We propose just such an automatic tuning scheme that takes a user-provided asymptotic bias tolerance and selects a step size that ensures it. The key to the method is a relationship we establish between the energy change during integration and the asymptotic bias. We show that this procedure rigorously bounds the asymptotic bias for Gaussian target distributions. We then numerically show that the procedure works beyond Gaussians. To demonstrate the practicality of the proposed scheme, we provide a comprehensive comparison of adjusted and unadjusted samplers on Bayesian inference problems and on a statistical physics model in more than one million parameters. With our tuning scheme, the unadjusted methods achieve close to optimal performance, significantly and consistently outperforming their adjusted counterparts.
Lay Summary
Making decisions based on scientific, engineering and financial outputs critically depends on their uncertainty. A common way to estimate uncertainty is to use sampling algorithms, but these methods can become extremely slow for modern problems involving millions of variables. Many state-of-the-art algorithms rely on a correction step to guarantee mathematically exact results. Making this correction typically significantly reduces efficiency, especially in high-dimensional problems. We show that this correction step can often be removed without sacrificing meaningful accuracy. We introduce an automatic tuning method that controls numerical errors while the algorithm runs, making these faster sampling methods practical without requiring expert manual tuning. Across sampling benchmarks and large-scale physics simulations with more than two million variables, our approach consistently outperforms existing methods, sometimes by more than a factor of one hundred. Our results suggest that large-scale uncertainty estimation and Bayesian inference can be made substantially faster and more practical in many real-world applications