Black-box unadjusted Hamiltonian Monte Carlo

## Abstract

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 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 error in the integration and asymptotic bias. For Gaussians, we show that this procedure rigorously bounds the asymptotic bias. We then numerically show that the procedure works beyond Gaussians, on typical Bayesian problems. To demonstrate the practicality of the proposed scheme, we provide a comprehensive comparison of adjusted and unadjusted samplers, showing that with our tuning scheme, the unadjusted methods achieve close to optimal performance and significantly and consistently outperform their adjusted counterparts.

## 1 Introduction

Sampling offers a way to compute expectation values $\mathbb{E}_{p}[f]=\int p(\boldsymbol{x})f(\boldsymbol{x})d\boldsymbol{x}$,
where $f(\boldsymbol{x})$ is some function of parameters $\boldsymbol{x}\in\mathbb{R}^{d}$ and $p(\boldsymbol{x})=e^{-\mathcal{L}(\boldsymbol{x})}/Z$ is a given probability density, with a possibly unknown normalization constant $Z=\int e^{-\mathcal{L}(\boldsymbol{x})}dx$.
This is a key tool across a variety of disciplines, from economics and social science (Gelman et al., 2013), to high energy physics (Duane
et al., 1987a), computational chemistry Leimkuhler and
Matthews (2015b), statistical physics (Leimkuhler and
Matthews, 2015a), and machine learning (Neal, 2012).
Often a gradient $\nabla\mathcal{L}(\boldsymbol{x})$ is available, either analytically, or via automatic differentiation (Griewank and
Walther, 2008; Margossian, 2019), as is often the case in Bayesian statistics (Štrumbelj
et al., 2024; Carpenter
et al., 2017), machine learning (Baydin et al., 2018), lattice quantum problems (Gattringer and
Lang, 2010), and cosmology (Campagne
et al., 2023; Ruiz-Zapatero
et al., 2024; Horowitz and
Lukic, 2025).

Markov Chain Monte Carlo (MCMC) (Metropolis et al., 1953) is a commonly employed class of sampling methods in which a Markov chain $\{\boldsymbol{x}_{i}\}_{i=1}^{n}$ is designed such that its stationary distribution is $\tilde{p}(\boldsymbol{x})$ and the expectation value $\mathbb{E}_{p}[f]$ can be approximated by the time average over the chain $\bar{f}=\frac{1}{n}\sum_{i=1}^{n}f(\boldsymbol{x}_{i})$. When a (smooth) gradient is available, Hamiltonian Monte Carlo (HMC) (Duane
et al., 1987b; Neal et al., 2011a; Betancourt, 2017) and underdamped Langevin Monte Carlo (Horowitz, 1991; Leimkuhler and
Matthews, 2015a) are the gold standard algorithms for generating Markov chain transitions.
Despite their success (Štrumbelj
et al., 2024), MCMC often presents a computational bottleneck (Gattringer and
Lang, 2010; Leimkuhler and
Matthews, 2015a; Simon-Onfroy et al., 2025), especially in the high-dimensional applications, forcing the practitioners to resort to more approximate methods, such as Laplace approximation (Millea and
Seljak, 2022) or an ensemble Kalman filter (Houtekamer and
Mitchell, 2005; Houtekamer and
Zhang, 2016).

#### Bias-variance tradeoff

Figure: Figure 1: Graphical representation of the bias decomposition from Equation ([2](https://arxiv.org/html/2412.08876v2#S1.E2)). Both adjusted and unadjusted methods have initialization bias, which typically decays exponentially fast for HMC-like algorithms (Margossian and Gelman, 2024) (with a rate denoted by $c$ in the Figure). Both methods also have the variance associated with the finite number of samples taken. It decays inversely proportionally to the number of samples $t$, with a proportionality constant that is determined by the autocorrelation time. Unadjusted methods also have asymptotic bias which depends strongly on the step size $\epsilon$, for the example from Appendix [C](https://arxiv.org/html/2412.08876v2#A3), as $\mathcal{O}(\epsilon^{4})$.
Refer to caption: x1.png

To understand the challenge faced by MCMC, consider the root mean squared error (RMSE) of the Monte Carlo estimator:

$$ $\mathrm{RMSE}[\bar{f}]=\mathbb{E}_{MC}[(\bar{f}-\mathbb{E}_{p}[f])^{2}]^{1/2},$ (1) $$

where $\mathbb{E}_{MC}[\cdot]$ denotes the expectation with respect to the Monte Carlo randomness. We will sometimes also refer to the relative RMSE, which we define as $\mathrm{RMSE}[\bar{f}]/\mathbb{E}_{p}[f]$.
The RMSE can be decomposed as

$$ $\mathrm{RMSE}[\bar{f}]^{2}=\mathrm{Bias}[\bar{f}]^{2}+\mathrm{Var}[\bar{f}],$ (2) $$

where $\mathrm{Bias}[\bar{f}]=\mathbb{E}_{MC}[\bar{f}-\mathbb{E}_{p}[f]]$ and the variance $\mathrm{Var}[\bar{f}]=\mathbb{E}_{MC}[(\bar{f}-\mathbb{E}_{MC}[\bar{f}])^{2}]$.
The bias term arises either from the chain not yet reaching its stationary state and being influenced by its initial position (initialization bias), or because the stationary distribution of the chain $\tilde{p}(\boldsymbol{x})$ does not equal the target distribution $p(\boldsymbol{x})$ (asymptotic bias). Initialization bias is an important issue (Margossian et al., 2024), but we will not study it here, noting that for long enough chains it can be eliminated by discarding initial samples (Margossian and
Gelman, 2024).
The variance term comes from the finite chain length and the correlation between the samples. This decomposition is illustrated in Figure [1](https://arxiv.org/html/2412.08876v2#S1.F1).

One of the main practical goals in designing MCMC methods is to obtain a desired RMSE at the lowest possible computational budget. This amounts to balancing the cost of generating samples, correlations between the samples, and the asymptotic bias. For example, the use of Metropolis-Hastings adjustment (Chib and
Greenberg, 1995) ensures that a chain satisfies detailed balance, so that the asymptotic bias vanishes, but for finite-length chains does not remove the variance. Noting that chains are finite in practice, one could negotiate the tradeoff differently and reduce the variance at the cost of introducing some bias. This possibility is the focus of the present work.

Figure: Figure 2: Sampling cost scaling with the dimensionality for product targets. Cost is measured by the number of gradient calls needed to achieve low error. MH adjusted methods are shown in dashed lines, and unadjusted methods are shown in solid lines. Step size for the unadjusted schemes is selected by the scheme proposed in this paper. $d^{1/4}$ power law grey lines are shown in the background. A Standard Gaussian target is shown on the left ($D=1$), and a product of banana-shaped Rosenbrock distributions ($D=2$) on the right. In both cases, the cost of the MH adjusted methods scales as $d^{1/4}$, while the cost of the unadjusted method remains constant with dimension.
Refer to caption: x2.png

#### Illustrative example

The shortcoming of performing MH is the scaling with the dimension that it implies for the sampler. To develop an intuition, we consider a simple problem, where the target is a product of $K$ independent $D$-dimensional distributions $q$:

$$ $p(\boldsymbol{x})=\prod_{i=1}^{K}q(x_{iD},x_{iD+1}\ldots x_{iD+D-1}),$ $$

where $D$ is a small number.
Suppose we are interested in the expectation value of a function of only one of those parameters, such as $f(\boldsymbol{x})=x_{1}^{2}$.
We will measure the performance as the number of evaluations of $\nabla p(\boldsymbol{x})$ (which is typically the bottleneck of computation and thus a proxy for the wall-clock time) that the sampler uses to get the relative RMSE (averaged over 128 chains) of $\mathbb{E}[f(\boldsymbol{x})]$ below $10\%$.
Figure [2](https://arxiv.org/html/2412.08876v2#S1.F2) shows performance as a function of the dimension of the problem $d=KD$. As can be seen, the number of gradient calls scales as $d^{1/4}$ for Metropolis adjusted methods, in accordance with (Beskos et al., 2013). This is because even though the Hamiltonian and Langevin dynamics operate independently on each copy of the distribution, the energy change $\Delta_{H}$ (defined below), and therefore the MH acceptance probability, involves a sum over all parameters. Optimal performance requires a fixed acceptance rate (Beskos et al., 2013; Neal
et al., 2011b), so to compensate for the increasing number of parameters, the step size needs to decrease, so that the integrator uses more gradient evaluations for a fixed trajectory length.
Unadjusted methods, by contrast, do not suffer from this scaling since they operate on each copy independently. They need to ensure that the bias is small, but this is independent of the number of parameters in these examples.

The example in Figure [2](https://arxiv.org/html/2412.08876v2#S1.F2) is idealized, but a similar scaling has been observed for mean field models (Durmus and
Eberle, 2023) and real problems, like cosmological field level inference (Simon-Onfroy et al., 2025).
Thus, particularly for high-dimensional problems, unadjusted samplers may be more efficient, in some cases by a large factor. Crucially, however, they are only viable if the step size can be chosen so that the asymptotic bias is small compared to the variance of the finite set of samples taken.

#### Our contributions

In this work, we propose an automated step size adaptation scheme for unadjusted HMC and LMC, making them usable in a black-box manner, that is, *without manually tuning any hyperparameters*. We expect this to have a significant impact on scientific applications that require gradient-based sampling in high dimensions.
The central idea is that the energy error arising from integrating with step size $\epsilon$ provides a measure of the asymptotic bias, and step size can be adaptively varied in a short initial chain to target an energy error which results in low asymptotic bias relative to the variance.
In Section [4](https://arxiv.org/html/2412.08876v2#S4) we show analytically for Gaussians that the energy error variance can be used to upper bound the bias, with no other knowledge about the target distribution.
In Section [5](https://arxiv.org/html/2412.08876v2#S5), we numerically confirm the analytical results and show that the same upper bound generally extends to non-Gaussian targets.
In Section [6](https://arxiv.org/html/2412.08876v2#S6) we use the bound to construct our adaptation algorithm and demonstrate its effectiveness on a range of standard benchmarks in Section[7](https://arxiv.org/html/2412.08876v2#S7).

## 2 Related work

State-of-the-art Markov chain kernel for densities with smooth gradient are Hamiltonian Monte Carlo (HMC) (Duane
et al., 1987b; Neal et al., 2011a; Betancourt, 2017) and *underdamped* Langevin Monte Carlo (LMC) (Leimkuhler and
Matthews, 2015a), not to be confused with *overdamped* Langevin dynamics and the corresponding Metropolis Adjusted Langevin Algorithm (MALA) (Rossky
et al., 1978).
In HMC and LMC, each parameter $x_{i}$ has an associated momentum variable $u_{i}$. The current state of the Markov chain $\boldsymbol{x}_{n}$ is used as an initial condition $\boldsymbol{x}(t=0)=\boldsymbol{x}_{n}$ of the dynamics, along with randomly drawn initial conditions for the momenta $u_{i}(0)\sim\mathcal{N}(0,1)$. The new state of the chain in HMC is then generated by solving for $\boldsymbol{x}(t=T)$, where $T$ is a predefined trajectory length, which is a hyperparameter of the algorithm, and $\boldsymbol{z}(t)\equiv(\boldsymbol{x}(t),\boldsymbol{u}(t))$ are solutions of the Hamiltonian equations with Hamiltonian function
$\mathcal{H}(\boldsymbol{z})=\frac{1}{2}\left\|\boldsymbol{u}\right\|^{2}+
\mathcal{L}(\boldsymbol{x})$:

$$ $\displaystyle\frac{d}{dt}\boldsymbol{x}(t)=\boldsymbol{u}(t)\qquad\frac{d}{dt} \boldsymbol{u}(t)=-\nabla\mathcal{L}(\boldsymbol{x}).$ (3) $$

Note that the value of the Hamiltonian function, also called the energy, is a conserved quantity of the Hamiltonian equations, meaning that the energy difference $\Delta_{H}(\boldsymbol{z}^{\prime},\boldsymbol{z})=H(\boldsymbol{z}^{\prime})-
H(\boldsymbol{z})$ vanishes for exact solutions of the Hamiltonian dynamics:
$\Delta_{H}(\boldsymbol{z}(t),\boldsymbol{z}(0))=0$.
In LMC, after every step with Hamiltonian dynamics, the velocity is partially randomized, which corresponds to the discretization of the Langevin stochastic differential equation (Leimkuhler and
Matthews, 2015a).

In statistics, adjusted versions of these dynamics are used almost exclusively. However,
unadjusted versions are used in some fields with high-dimensional distributions, such as Lattice quantum chromodynamics (Lüscher, 2018; Clark and
Kennedy, 2007) and Molecular Dynamics (Leimkuhler and
Matthews, 2015a). Commonly, underdamped Langevin dynamics or Nosé-Hoover thermostat (Evans and
Holian, 1985) are employed.
Another promising option is Microcanonical Langevin Monte Carlo (MCLMC) (Robnik et al., 2024; Robnik and
Seljak, 2024; Tuckerman et al., 2001; Ver Steeg and
Galstyan, 2021), which makes use of velocity norm preserving dynamics.
Domain knowledge and trial runs are used in these fields to select an appropriate step size, which ensures a sufficiently small asymptotic bias.
Unadjusted methods have also been analyzed theoretically, establishing a bound on the asymptotic bias (Durmus and
Eberle, 2023) and the mixing time (Bou-Rabee and
Eberle, 2023; Camrud
et al., 2024). For mean-field models, the bound on both is dimension-free (Bou-Rabee and
Schuh, 2023) and thus unadjusted methods provably outperform adjusted methods in this case. However, these bounds assume some global knowledge of the target distribution that is not available in practice.
A notable gap in the above work is an algorithm for choosing the step size $\epsilon$. Without a principled way to do this, general practitioners cannot use the unadjusted algorithm in a black-box fashion, which gravely limits their applicability.

## 3 Measuring asymptotic bias

#### Origin of bias

Generically, Hamiltonian equations like ([3](https://arxiv.org/html/2412.08876v2#S2.E3)) cannot be solved exactly, so a numerical integrator like velocity Verlet (Leimkuhler and
Matthews, 2015a) is used to approximate it.
Velocity Verlet is an example of a splitting method, where to solve the dynamics numerically, one first analytically solves for $\boldsymbol{x}$ at fixed $\boldsymbol{u}$ and vice versa. The first solution is called the position update and is given by $\Phi^{T}_{\epsilon}(\boldsymbol{z})=(\boldsymbol{x}+\epsilon\boldsymbol{u},
\boldsymbol{u})$ for Equation ([3](https://arxiv.org/html/2412.08876v2#S2.E3)), the second is called the velocity update and is given by $\Phi^{V}_{\epsilon}(\boldsymbol{z})=(\boldsymbol{x},\boldsymbol{u}-\epsilon
\nabla\mathcal{L}(\boldsymbol{x}))$. The joint solution is then approximated by a composition of these maps, which for velocity Verlet is,

$$ $\boldsymbol{z}(t+\epsilon)\approx\Phi_{\epsilon}(\boldsymbol{z}(t))=(\Phi^{V}_ {\epsilon/2}\circ\Phi^{T}_{\epsilon}\circ\Phi^{V}_{\epsilon/2})(\boldsymbol{z} (t)).$ (4) $$

Due to this approximation, the stationary distribution $\tilde{p}(\boldsymbol{x})$ of the sampler no longer equals the desired target distribution $p(\boldsymbol{x})$ and expectation values acquire an asymptotic bias, which vanishes in the limit of step size going to zero (Durmus and
Eberle, 2023).

In the case of LMC, the update is additionally complemented by the partial refreshments of the velocity (Leimkuhler and
Matthews, 2015b):

$$ $\boldsymbol{z}(t+\epsilon)\approx(\Phi^{O}_{\epsilon/2}\Phi_{\epsilon}\circ \Phi^{O}_{\epsilon/2})(\boldsymbol{z}(t)),$ (5) $$

where $\Phi^{O}_{\epsilon}(\boldsymbol{z})=(\boldsymbol{x},\,e^{-\epsilon/L}
\boldsymbol{u}+(\sqrt{1-e^{-2\epsilon/L}})\textbf{n})$ and $\textbf{n}\sim\mathcal{N}(0,I)$. Parameter $L$ determines the amount of momentum refreshment and plays a similar role as the trajectory length in HMC.

#### Choice of expectation

When minimizing bias, it is important to determine the expectation with respect to which the bias is defined. For instance, an important expectation in Bayesian statistics (where quantification of uncertainty is key) is of the second moments, so that we are concerned with expectations of the form $\mathbb{E}_{p}[x_{i}x_{j}]$, or more generally the covariance matrix $\Sigma_{p}=\mathbb{E}_{q}[(\boldsymbol{x}-\mathbb{E}[\boldsymbol{x}])(
\boldsymbol{x}-\mathbb{E}[\boldsymbol{x}])^{T}]$.
For Gaussian distributions, the covariance matrix
contains all the information about the posterior, but even for non-Gaussian distributions, they are often used as a posterior summary.

From a theoretical perspective, it may also be interesting to consider the bias of a wider set of functions, such as any Lipschitz-continous functions.
In what follows, we consider those two cases.

#### Covariance matrix error

To quantify expectation value error of the second moments, we introduce a scalar measure of the covariance matrix error, which we define as

$$ $b_{\mathit{cov}}^{2}(\Sigma_{p},\Sigma_{q})\equiv\frac{1}{d}\Tr{(I-\Sigma^{-1} _{p}\Sigma_{q})^{2}},$ (6) $$

where $\Sigma_{p}$ is the true covariance matrix of the target distribution $p$ and $\Sigma_{q}$ is its estimate under some other distribution $q$. In the simple case where the covariance matrices are diagonal, $b_{\mathit{cov}}$ is the relative error of the variance estimate, averaged over the parameters, i.e., $b_{\mathit{cov}}^{2}(\Sigma_{p},\Sigma_{q})=\frac{1}{d}\sum_{i=1}^{d}([\Sigma_
{p}]_{ii}-[\Sigma_{q}]_{ii})^{2}/[\Sigma_{p}]_{ii}^{2}$.
We emphasize that this is the convergence metric that is often used in practice (Grumitt
et al., 2022; Robnik et al., 2022).
The diagonal form is preferable in high dimensions as it does not require storage of the full covariance matrix.
Nonetheless, we advocate measuring error with ([6](https://arxiv.org/html/2412.08876v2#S3.E6)) when feasible
because it additionally penalizes the off-diagonal terms and has a number of nice properties:

- •
It is a divergence on the space of positive-definite matrices, meaning that it is non-negative and zero if and only if the two matrices are the same.
- •
It can be related to the effective sample size in the following way: for an empirical covariance matrix $\bar{\Sigma}$ estimated from $n$ exact samples from the Gaussian target distribution $p=\mathcal{N}(0,\Sigma)$,
the expected value of $b_{\mathit{cov}}^{2}$ depends only on the number of samples $n$ and not on $\Sigma$:
$\mathbb{E}_{MC}[b_{\mathit{cov}}^{2}(\Sigma,\bar{\Sigma})]=(d+1)/n,$
(7)
So for Gaussian distributions, $b_{\mathit{cov}}^{2}$ defines the effective sample size of the estimate, that is, the number of exact samples that would yield the same error $b_{\mathit{cov}}^{2}$. Concretely, given a target distribution $p$ in $d=9$, suppose it takes $1000$ samples from a Markov chain to achieve $b_{\mathit{cov}}^{2}=0.1$. This would correspond to $(9+1)/0.1=100$ effective samples, so $0.1$ effective samples per step.
- •
It is invariant to linear change of basis of the configuration space.

Proofs of these properties are provided in Appendix [A](https://arxiv.org/html/2412.08876v2#A1).

#### Wasserstein distance

The Wasserstein distance $\mathcal{W}_{\nu}(p,q)$ between densities $p$ and $q$ is (Kantorovich, 1960):

$$ $\mathcal{W}_{\nu}(p,q)=\big{(}\mathrm{inf}_{\pi\in\Pi(p,q)}\int\left\| \boldsymbol{x}-\boldsymbol{x}^{\prime}\right\|^{\nu}\pi(\boldsymbol{x}, \boldsymbol{x}^{\prime})d\boldsymbol{x}d\boldsymbol{x}^{\prime}\big{)}^{1/\nu},$ (8) $$

where $\Pi(p,q)$ is the set of probability densities on $\mathbb{R}^{d}\times\mathbb{R}^{d}$ with marginals $p$ and $q$.
$\mathcal{W}_{\nu}(p,q)$ can be interpreted as the optimal cost of transporting one density to the other: $\left\|\boldsymbol{x}-\boldsymbol{x}^{\prime}\right\|^{\nu}$ is the cost of transporting the probability mass, and $\pi$ defines a strategy of the transport, i.e. $\pi(\boldsymbol{x},\boldsymbol{x}^{\prime})$ is the probability density for taking a piece at location $\boldsymbol{x}^{\prime}$ and moving it to the location $\boldsymbol{x}^{\prime}$.

Wasserstein distance has several nice properties: it is invariant to change of basis, is a metric on the space of distributions, and most importantly in the present context, it upper bounds the bias of Lipschitz-continuous functions. That is, for any Lipschitz-continuous function $f$, with Lipschitz constant $L$ (meaning that $|f(\boldsymbol{x})-f(\boldsymbol{x}^{\prime})|<L\left\|\boldsymbol{x}-
\boldsymbol{x}^{\prime}\right\|$ for any $\boldsymbol{x},\boldsymbol{x}^{\prime}\in\mathbb{R}^{d}$), the bias associated with this function is upper bounded by $\mathcal{W}_{1}(p,\tilde{p})$, which is in turn upper bounded by $\mathcal{W}_{2}(p,\tilde{p})$:

$$ $\mathbb{E}_{p}[f]-\mathbb{E}_{\tilde{p}}[f]\leq L\mathcal{W}_{1}(p,\tilde{p}) \leq L\mathcal{W}_{2}(p,\tilde{p}),$ (9) $$

The first inequality is Kantorovich-Rubinstein duality (Villani, 2003), and the second is a consequence of Jensen’s inequality. This work will focus on $\mathcal{W}_{2}$.

#### Energy error variance per dimension

Our goal is to control $b_{\mathit{cov}}^{2}$ and $\mathcal{W}_{2}$ in unadjusted samplers. The first quantity is of practical interest in various fields like Bayesian inference, and the second provides a bound on the bias of all Lipschitz continuous functions. We will do this by monitoring the energy error $\Delta_{H}(\Phi_{\epsilon}(\boldsymbol{z}(t)),\boldsymbol{z}(t))$ induced by numerical integration with step size $\epsilon$. We define Energy Error Variance Per Dimension (EEVPD),

$$ $\mathrm{EEVPD}=\mathrm{Var}_{\boldsymbol{x}\sim\tilde{p},\boldsymbol{u}\sim N( 0,I)}[\Delta_{H}(\Phi_{\epsilon}(\boldsymbol{x},\boldsymbol{u}),(\boldsymbol{x },\boldsymbol{u}))]/d.$ (10) $$

This is a quantity that can easily be estimated in practice: $\boldsymbol{x}\sim\bar{p}$, $\boldsymbol{u}\sim\mathcal{N}(0,I)$ is the stationary distribution of the chain, so computing $\mathrm{EEVPD}$ amounts to collecting the samples from the stationary chain, evaluating the one-step energy error for each of those samples and computing their variance. This can be done online, using a running average of the first and second moment, so that the step size $\epsilon$ can be adaptively varied to target a desired value of $\mathrm{EEVPD}$.

## 4 Analytic results for Gaussian distributions

We will start by showing that EEVPD can be used to control the asymptotic covariance matrix bias $b^{2}_{\mathit{cov}}(\Sigma,\tilde{\Sigma})$ or Wasserstein distance $\mathcal{W}_{2}(p,\tilde{p})$
for Gaussian target distributions, $p=\mathcal{N}(0,\Sigma)$.The key tool is that the stationary distribution $\tilde{p}$ of unadjusted HMC and LMC for Gaussians is known exactly (Monmarché, 2022) and is also Gaussian, with $\tilde{p}=\mathcal{N}(0,\tilde{\Sigma})$. $\tilde{\Sigma}$ has the same eigenvectors as $\Sigma$, but its eigenvalue associated to the $i$-th eigenvector is

$$ $\tilde{\sigma}_{i}^{2}=\frac{\sigma_{i}^{2}}{1-\epsilon^{2}/4\sigma_{i}^{2}},$ (11) $$

where $\sigma^{2}_{i}$ is the corresponding eigenvalue of $\Sigma$.
Therefore the $\mathrm{EEVPD}$ has a closed form:

###### Lemma 4.1 .

For a Gaussian distribution with covariance matrix eigenvalues $\{\sigma_{i}^{2}\}_{i=1}^{d}$ and an HMC or LMC sampler with a stable velocity Verlet integrator, meaning that step size $\epsilon<2\min_{i}{\sigma_{i}}$,

$$ $\mathrm{EEVPD}=\frac{1}{d}\sum_{i=1}^{d}E(\epsilon^{2}/\sigma_{i}^{2}),$ $$

where $E(y)=\frac{y^{3}}{16(1-y/4)}$.

Our key results are then:

###### Theorem 4.2 .

For a Gaussian distribution $\mathcal{N}(0,\Sigma)$ with covariance matrix eigenvalues $\{\sigma_{i}^{2}\}_{i=1}^{d}$ and HMC or LMC sampler with a stable velocity Verlet integrator, meaning that the step size $\epsilon<2\min_{i}{\sigma_{i}}$, the covariance matrix bias is upper bounded by

$$ $b^{2}_{\mathit{cov}}(\Sigma,\tilde{\Sigma})\leq\varphi^{-1}(\mathrm{EEVPD}),$ $$

as long as $\mathrm{EEVPD}<0.397$.
Here,

$$ $\varphi(x)=\frac{4x^{3/2}}{(1+x^{1/2})^{2}}.$ $$

Similarly, the Wasserstein distance between the target and stationary distributions is upper bounded by

$$ $\mathcal{W}_{2}(p,\tilde{p})^{2}/d\leq\epsilon^{2}\varphi_{W}^{-1}(\mathrm{ EEVPD}),$ $$

as long as $\mathrm{EEVPD}<6.75$. Here,
$\varphi_{W}\equiv E\circ W^{-1}$ with
$W(y)=\frac{2(1-y/8-\sqrt{1-y/4})}{y(1-y/4)}$.
Both bounds are sharp and realized if and only if the target is isotropic, i.e. $\Sigma_{p}\propto I$.

All the proofs are given in Appendix [B](https://arxiv.org/html/2412.08876v2#A2). Note that the conditions on EEVPD are not a severe limitation in practice, see for example Table [1](https://arxiv.org/html/2412.08876v2#S6.T1), where significantly lower values of EEVPD are used.

###### Lemma 4.1 .

###### Theorem 4.2 .

## 5 Numeric results for non-Gaussian distributions

Figure: Figure 3: The asymptotic covariance matrix error $b_{\mathit{cov}}(\Sigma,\tilde{\Sigma})$ as a function of $\mathrm{EEVPD}$. Unadjusted HMC is shown on the left, and unadjusted MCLMC on the right. The relation is shown for various inference problems from Section [5](https://arxiv.org/html/2412.08876v2#S5). The analytical equality for the Standard Gaussians from Theorem [4.2](https://arxiv.org/html/2412.08876v2#S4.Thmtheorem2) is shown in black and agrees perfectly with the numerical results for HMC (left). The inequality for arbitrary Gaussian distributions is shown as a shaded grey region, and also applies perfectly for Gaussians. As can be seen, most non-Gaussian targets also abide by this inequality.
Refer to caption: x3.png

We now examine the validity of Theorem [4.2](https://arxiv.org/html/2412.08876v2#S4.Thmtheorem2) for standard Bayesian inference benchmark problems. We do not numerically verify the part of the theorem about the Wasserstein distance, due to the difficulty of numerically computing the Wasserstein distance for the problems considered here.
Even though the exact asymptotic error-$\mathrm{EEVPD}$ relation of Theorem [4.2](https://arxiv.org/html/2412.08876v2#S4.Thmtheorem2) applies only for HMC and LMC, unadjusted MCLMC has a notion of energy (Robnik et al., 2024), so we will also test MCLMC here.

Benchmark problem descriptions are in Appendix [D.1](https://arxiv.org/html/2412.08876v2#A4.SS1). For each problem, we show the asymptotic value of $b_{\mathit{cov}}(\Sigma,\tilde{\Sigma})^{2}$ as a function of $\mathrm{EEVPD}$ in Figure [3](https://arxiv.org/html/2412.08876v2#S5.F3). Asymptotic $b_{\mathit{cov}}$ is computed by running unadjusted chains with different step sizes, each using $10^{8}$ gradient calls.
We eliminate the initial $10^{4}$ calls to eliminate the initialization bias and use the subsequent samples to compute the expectation values for the covariance matrix and $\mathrm{EEVPD}$. We monitor $b_{\mathit{cov}}(\Sigma,\tilde{\Sigma})^{2}$ from Equation ([6](https://arxiv.org/html/2412.08876v2#S3.E6)) and check that it has converged to the asymptotic value. If the convergence has not yet been achieved, i.e., the bias is still decaying, we do not show these measurements on the plots. This happens for some of the harder problems at small step sizes, where the chains are not long enough for the variance to become negligible. We have checked that the variation between the chains is negligible, but nonetheless average the results over 4 independent chains.

Numerical results for unadjusted HMC on Gaussians agree perfectly with Theorem [4.2](https://arxiv.org/html/2412.08876v2#S4.Thmtheorem2): for the standard Gaussian, the equality holds, while for the Ill-conditioned Gaussian, the inequality holds for $\mathrm{EEVPD}<0.397$ as per the theorem.
The inequality also applies to the majority of non-Gaussian benchmark problems, illustrating its broader applicability.
One exception is the Brownian motion example, where it is off by approximately a factor of $1.5$ at small $\epsilon$, meaning that one would think one has $<2\%$ asymptotic error, when in fact it was $3\%$.
The Rosenbrock and Funnel examples are only shown at small step sizes, because the problem becomes numerically unstable at higher step sizes, incurring divergences.
Very similar results also apply to MCLMC, except that the bias at a fixed $\mathrm{EEVPD}$ is usually lower for MCLMC. This suggests that the tuning scheme we develop for unadjusted HMC can also be applied to MCLMC, but a slightly larger $\mathrm{EEVPD}$ may be used.

## 6 Automatic tuning scheme

Figure: Figure 4: Posterior density for the funnel problem. The 2D marginal distribution in the $\theta-z_{1}$ plane and the corresponding 1D marginals are shown. The contours are obtained by kernel density estimation, and samples are shown as dots. The ground truth, obtained by a very long NUTS chain, is shown in orange, and unadjusted MCLMC is shown in blue. Both chains are $10^{7}$ samples long to eliminate the variance error. The two methods give practically indistinguishable posteriors, demonstrating that the discretization bias was successfully suppressed by the $\mathrm{EEVPD}$ control.
Refer to caption: extracted/6448359/funnel.png

The core result of sections [4](https://arxiv.org/html/2412.08876v2#S4) and [5](https://arxiv.org/html/2412.08876v2#S5) is that by controlling $\mathrm{EEVPD}$, we can in turn control the asymptotic bias of expectations of interest.
Our tuning scheme for step size is therefore straightforward: for any unadjusted sampler with an appropriate notion of energy (here we consider unadjusted HMC, LMC and MCLMC), we keep a running estimate of $\mathrm{EEVPD}$, and adaptively vary $\epsilon$ to target a desired value of $\mathrm{EEVPD}$. In practice, an optimization algorithm such as dual averaging (Nesterov, 2009; Hoffman and
Gelman, 2014) can be used here, but the choice of optimizer has little effect on the results, we use the scheme in Appendix [D.2](https://arxiv.org/html/2412.08876v2#A4.SS2).
It remains to select the desired bias tolerance and the corresponding $\mathrm{EEVPD}$. This depends on the application and the accuracy requirements; we, nonetheless, provide some guidance. Certainly, asymptotic bias squared should be smaller than the mean squared error tolerance (because mean squared error is composed of the bias and the variance, which are both non-negative), but by how much? Appendix [C](https://arxiv.org/html/2412.08876v2#A3) shows that for estimating $\mathbb{E}[x^{2}]$ of a standard normal distribution with finite chain length of unadjusted HMC, the optimal squared bias should be one fifth of the mean squared error tolerance. Adopting this prescription gives the conversion table:

**Table 1: Tabulated values of $\mathrm{EEVPD}$ (third column) that ensure a desired asymptotic covariance matrix bias (second column). A useful, quick to compute approximation for small $b^{2}_{\mathit{cov}}$ is $\mathrm{EEVPD}=\varphi(b^{2}_{\mathit{cov}})\approx 4b_{\mathit{cov}}^{3}$. Optimal asymptotic bias should be smaller than the given relative root mean square error tolerance (Equation ([1](https://arxiv.org/html/2412.08876v2#S1.E1))), for example, one can use the prescription $\mathrm{Bias}^{2}<\mathrm{RMSE}^{2}/5$ from Appendix [C](https://arxiv.org/html/2412.08876v2#A3). This error is given in the first column.**
| relative RMSE tolerance | Bias tolerance | $\mathrm{EEVPD}$ |
| --- | --- | --- |
| 50% | 22% | $3.0\times 10^{-2}$ |
| 10% | 4.5% | $3.3\times 10^{-4}$ |
| 5% | 2.2% | $4.3\times 10^{-5}$ |
| 1% | 0.45% | $3.5\times 10^{-7}$ |

#### Momentum decoherence scale

A final hyperparameter of unadjusted HMC is the number of steps between momentum refreshes, and for unadjusted LMC, the amount of noise in each partial refresh.
Both can be understood in terms of a single hyperparameter, the *momentum decoherence length* $L$. For HMC, $L$ is simply the number of steps in a trajectory times the step size, since this is precisely the length along the trajectory after which momentum completely decoheres. For LMC, $L$ is given in section [3](https://arxiv.org/html/2412.08876v2#S3). In either case, we select $L$ based on the autocorrelation length of the the chain, as in (Robnik et al., 2024; Robnik
et al., 2025).

#### Demonstration of bias control

In Figure [4](https://arxiv.org/html/2412.08876v2#S6.F4) we show the posterior density for the funnel problem, obtained by the unadjusted MCLMC algorithm. We require asymptotic bias of $1\%$ and use Theorem ([4.2](https://arxiv.org/html/2412.08876v2#S4.Thmtheorem2)) to get the corresponding $\mathrm{EEVPD}=\varphi(0.01)=4\times 10^{-6}$. We use dual averaging from (Hoffman
et al., 2014) to determine the MCLMC step size, targeting this $\mathrm{EEVPD}$.
The posterior is practically indistinguishable from a very long NUTS chain, showing that the discretization error is very small, as desired.

## 7 Experiments

Our central contribution is a scheme which makes unadjusted HMC, LMC and MCLMC *black-box samplers*, in the sense of requiring no manual tuning on the part of a user. It is therefore natural to compare performance to the state of the art black-box sampler, the No-U-Turn Sampler (NUTS) (Hoffman
et al., 2014). In addition, it is of interest to compare unadjusted samplers to their adjusted counterparts.

#### Samplers

With this in mind, we report results on NUTS, unadjusted LMC (uLMC), adjusted LMC (aLMC), in particular, the version proposed in (Riou-Durand and
Vogrinc, 2023), which is termed Metropolis Adjusted Langevin Trajectories (MALT), unadjusted MCLMC (uMCLMC), and adjusted MCLMC (aMCLMC).

While the other algorithms are evaluated with their respective tuning schemes, for aLMC, we perform a grid search to demonstrate that black-box uLMC outperforms even optimal aLMC. We omit reporting of unadjusted and adjusted HMC, since the performance and implementation closely resemble that of LMC, and LMC is generally considered the preferred option (Riou-Durand and
Vogrinc, 2023).
Further details of the experiment are provided in Appendix [D.3](https://arxiv.org/html/2412.08876v2#A4.SS3), for the code, see Appendix [F](https://arxiv.org/html/2412.08876v2#A6).

For the unadjusted methods, we use the scheme from Section [6](https://arxiv.org/html/2412.08876v2#S6) and take $\mathrm{EEVPD}$ of $3\times 10^{-4}$, corresponding to the RMSE = $10\%$, which will be our notion of convergence. For MCLMC we use a slightly larger value of $5\times 10^{-4}$, as suggested by Figure [3](https://arxiv.org/html/2412.08876v2#S5.F3). In Appendix [E](https://arxiv.org/html/2412.08876v2#A5), we perform an ablation study for LMC, to examine the change in performance as desired $\mathrm{EEVPD}$ is varied. The performance does not change much in a range of reasonable $\mathrm{EEVPD}$, and the value of $3\times 10^{-4}$ is conservative, in the sense that larger values improve performance. The only exception is Stochastic Volatility, where we find that a smaller value is needed; we use $5\times 10^{-7}$ for HMC, and $2\times 10^{-8}$ for MCLMC.

#### Evaluation metric

In addition to $b^{2}_{\mathit{cov}}(\Sigma,\Sigma_{\mathrm{sampler}})$ we consider a more standard metric of convergence: following Hoffman and
Sountsov (2022a) we define the squared error of the expectation value $\mathbb{E}[f(\boldsymbol{x})]$ as

$$ $b^{2}(f)=\frac{(\mathbb{E}_{\mathrm{sampler}}[f]-\mathbb{E}[f])^{2}}{\mathrm{ Var}[f]},$ (12) $$

and consider the average second-moment error across parameters, $b^{2}_{\mathit{avg}}\equiv d^{-1}\sum_{i=1}^{d}b^{2}(x_{i}^{2})$. $b^{2}_{\mathit{avg}}$ can be interpreted as the accuracy equivalent to 100 effective samples (Hoffman and
Sountsov, 2022a). For the higher dimensional problems (Item response and Stochastic Volatility) we only show this metric, because the covariance matrix is impractically large.
In typical applications, computing the gradients $\nabla\log p(\boldsymbol{x})$ dominates the total sampling cost, so we take the number of gradient evaluations as a proxy of a wall-clock time.
As in (Hoffman and
Sountsov, 2022b), we measure the sampler’s performance as the number of gradient calls $n$ needed to achieve low error, $b^{2}_{\mathit{avg}}<0.01$ or $b^{2}_{\mathit{cov}}<0.01$.

#### Benchmarks

We use a set of benchmark problems, mostly Bayesian hierarchical inference problems adapted from the Inference Gym (Sountsov
et al., 2020).
Problems vary in dimensionality (36–2429) and are both synthetic (the first three problems in Table [2](https://arxiv.org/html/2412.08876v2#S7.T2)) and with real data (the last three problems).

#### Results

Tables [2](https://arxiv.org/html/2412.08876v2#S7.T2) and [3](https://arxiv.org/html/2412.08876v2#S7.T3) show the results. We see that the unadjusted algorithms tend to perform better than their adjusted counterparts. This is especially impressive since the hyperparameters of the adjusted samplers were found by grid search, or are otherwise near optimal Robnik
et al. (2025).
More to the point, the performance advantage of the unadjusted methods over NUTS, in particular uMCLMC, is quite striking.
We estimate the standard deviation of results by a bootstrapping procedure (see Table [4](https://arxiv.org/html/2412.08876v2#A4.T4)).

Tables [2](https://arxiv.org/html/2412.08876v2#S7.T2) and [3](https://arxiv.org/html/2412.08876v2#S7.T3) also show the results where a grid search has found optimal values of the hyperparameters ($L,\epsilon$) for unadjusted LMC. This shows that the performance of LMC with automatic tuning is close to optimal. Notable exceptions are Rosenbrock and Item response where performance is off by a factor of two, due to the conservative choice of $\mathrm{EEVPD}$ that we are using (see Appendix [E](https://arxiv.org/html/2412.08876v2#A5)).

**Table 2: Number of gradient calls needed to get $b^{2}_{\mathit{avg}}$ below 0.01. Lower is better.**
|  | aLMC | uLMC | aMCLMC | uMCLMC | uLMC (gridsearch) | NUTS |
| --- | --- | --- | --- | --- | --- | --- |
| Standard Gaussian | 803 | 563 | 436 | 246 | 568 | 2391 |
| Rosenbrock | 19,862 | 16,820 | 18,214 | 10,688 | 8410 | 27,070 |
| Brownian | 4667 | 2168 | 2876 | 1628 | 2407 | 5334 |
| German Credit | 7924 | 4730 | 6123 | 3960 | 4423 | 10,484 |
| Item Response | 5234 | 2020 | 5470 | 1612 | 930 | 6944 |
| Stochastic Volatility | 37,904 | 40,131 | 38,357 | 16,854 | 26,982 | 30,234 |

**Table 3: Number of gradient calls needed to get $b^{2}_{\mathit{cov}}$ below 0.01.**
|  | aLMC | uLMC | aMCLMC | uMCLMC | uLMC (gridsearch) | NUTS |
| --- | --- | --- | --- | --- | --- | --- |
| Standard Gaussian | 79,841 | 64,254 | 43,389 | 26,032 | 65,604 | 240,456 |
| Rosenbrock | 659,359 | 415,988 | 540866 | 348,048 | 271,825 | 852,135 |
| Brownian | 93,787 | 112,242 | 76,931 | 41,838 | 73,632 | 146,333 |
| German Credit | 462,843 | 380,569 | 462,770 | 249,674 | 306904 | 756,792 |

## 8 Conclusions

In this paper, we propose a scheme that makes unadjusted gradient-based samplers practical. We show that unadjusted versions of LMC and MCLMC considerably outperform optimally tuned versions of their adjusted counterparts and, more importantly, the state-of-the-art NUTS sampler. The advantage is particularly pronounced for high-dimensional problems (Figure [2](https://arxiv.org/html/2412.08876v2#S1.F2)).
The practical implications of our work extend to domains where unadjusted methods are already employed for sampling tasks, such as computational chemistry, enabling more informed choices of hyperparameters and minimizing bias in sampling results.
This scheme could also be used in Bayesian neural networks, even in the presence of minibatching noise (Welling and
Teh,), but this is beyond the scope of this paper.

We note that while the proposed scheme automatically selects a step size that in general performs well, it is certainly possible that there are models where bias is larger than expected. We therefore encourage the practice of validating results common in numerical analysis (Dalla Brida and
Lüscher, 2017): rerun the sampler with a smaller step size and check that results do not change significantly. Since the asymptotic bias depends strongly on the step size this provides strong evidence that the asymptotic bias is negligible.

## References

- Baydin et al. (2018)
Baydin, A. G., B. A. Pearlmutter, A. A. Radul, and J. M. Siskind (2018).
Automatic differentiation in machine learning: A survey.
Journal of Machine Learning Research 18.
Publisher: Journal of Machine Learning Research.
- Beskos et al. (2013)
Beskos, A., N. Pillai, G. Roberts, J.-M. Sanz-Serna, and A. Stuart (2013,
November).
Optimal tuning of the hybrid Monte Carlo algorithm.
Bernoulli 19(5A), 1501–1534.
Publisher: Bernoulli Society for Mathematical Statistics and
Probability.
- Betancourt (2017)
Betancourt, M. (2017).
A conceptual introduction to hamiltonian monte carlo.
arXiv preprint arXiv:1701.02434.
- Bou-Rabee and
Eberle (2023)
Bou-Rabee, N. and A. Eberle (2023, February).
Mixing time guarantees for unadjusted Hamiltonian Monte Carlo.
Bernoulli 29(1), 75–104.
Publisher: Bernoulli Society for Mathematical Statistics and
Probability.
- Bou-Rabee and
Schuh (2023)
Bou-Rabee, N. and K. Schuh (2023, January).
Convergence of unadjusted Hamiltonian Monte Carlo for
mean-field models.
Electronic Journal of Probability 28(none), 1–40.
Publisher: Institute of Mathematical Statistics and Bernoulli
Society.
- Cabezas
et al. (2024)
Cabezas, A., A. Corenflos, J. Lao, and R. Louf (2024).
Blackjax: Composable Bayesian inference in JAX.
- Campagne
et al. (2023)
Campagne, J.-E., F. Lanusse, J. Zuntz, A. Boucaud, S. Casas, M. Karamanis,
D. Kirkby, D. Lanzieri, Y. Li, and A. Peel (2023, April).
JAX-COSMO: An End-to-End Differentiable and GPU
Accelerated Cosmology Library.
The Open Journal of Astrophysics 6,
10.21105/astro.2302.05163.
arXiv:2302.05163 [astro-ph].
- Camrud
et al. (2024)
Camrud, E., A. Durmus, P. Monmarché, and G. Stoltz (2024, May).
Second order quantitative bounds for unadjusted generalized
Hamiltonian Monte Carlo.
arXiv:2306.09513 [math].
- Carpenter
et al. (2017)
Carpenter, B., A. Gelman, M. D. Hoffman, D. Lee, B. Goodrich, M. Betancourt,
M. Brubaker, J. Guo, P. Li, and A. Riddell (2017, January).
Stan: A Probabilistic Programming Language.
Journal of Statistical Software 76, 1–32.
- Chib and
Greenberg (1995)
Chib, S. and E. Greenberg (1995).
Understanding the metropolis-hastings algorithm.
The american statistician 49(4), 327–335.
- Clark and
Kennedy (2007)
Clark, M. A. and A. D. Kennedy (2007, October).
Asymptotics of fixed point distributions for inexact Monte Carlo
algorithms.
Physical Review D 76(7), 074508.
Publisher: American Physical Society.
- Dalla Brida and
Lüscher (2017)
Dalla Brida, M. and M. Lüscher (2017, May).
SMD-based numerical stochastic perturbation theory.
The European Physical Journal C 77(5), 308.
- Duane
et al. (1987a)
Duane, S., A. D. Kennedy, B. J. Pendleton, and D. Roweth (1987a, September).
Hybrid Monte Carlo.
Physics Letters B 195(2), 216–222.
- Duane
et al. (1987b)
Duane, S., A. D. Kennedy, B. J. Pendleton, and D. Roweth (1987b).
Hybrid monte carlo.
Physics letters B 195(2), 216–222.
- Durmus and
Eberle (2023)
Durmus, A. O. and A. Eberle (2023, April).
Asymptotic bias of inexact Markov Chain Monte Carlo methods
in high dimension.
arXiv:2108.00682 [cs, math, stat].
- Evans and
Holian (1985)
Evans, D. J. and B. L. Holian (1985).
The nose-hoover thermostat.
Journal of Chemical Physics 83(8), 4069–4074.
- Gattringer and
Lang (2010)
Gattringer, C. and C. B. Lang (2010).
Quantum Chromodynamics on the Lattice: An Introductory
Presentation, Volume 788 of Lecture Notes in Physics.
Berlin, Heidelberg: Springer Berlin Heidelberg.
- Gelman et al. (2013)
Gelman, A., J. B. Carlin, H. S. Stern, D. B. Dunson, A. Vehtari, and D. B.
Rubin (2013).
Bayesian data analysis third edition.
Chapman and Hall/CRC.
- Gouraud
et al. (2023)
Gouraud, N., P. L. Bris, A. Majka, and P. Monmarché (2023, June).
HMC and underdamped Langevin united in the unadjusted convex
smooth case.
arXiv:2202.00977 [math, stat].
- Griewank and
Walther (2008)
Griewank, A. and A. Walther (2008, January).
Evaluating Derivatives: Principles and Techniques of
Algorithmic Differentiation, Second Edition.
Society for Industrial and Applied Mathematics.
Edition: Second.
- Grumitt
et al. (2022)
Grumitt, R. D., B. Dai, and U. Seljak (2022).
Deterministic langevin monte carlo with normalizing flows for
bayesian inference.
arXiv preprint arXiv:2205.14240.
- Hoffman and
Gelman (2014)
Hoffman, M. D. and A. Gelman (2014).
The No-U-Turn Sampler: Adaptively Setting Path
Lengths in Hamiltonian Monte Carlo.
Journal of Machine Learning Research 15(47),
1593–1623.
- Hoffman
et al. (2014)
Hoffman, M. D., A. Gelman, et al. (2014).
The no-u-turn sampler: adaptively setting path lengths in hamiltonian
monte carlo.
J. Mach. Learn. Res. 15(1), 1593–1623.
- Hoffman and
Sountsov (2022a)
Hoffman, M. D. and P. Sountsov (2022a, May).
Tuning-Free Generalized Hamiltonian Monte Carlo.
In Proceedings of The 25th International Conference on
Artificial Intelligence and Statistics, pp.  7799–7813. PMLR.
ISSN: 2640-3498.
- Hoffman and
Sountsov (2022b)
Hoffman, M. D. and P. Sountsov (2022b).
Tuning-free generalized hamiltonian monte carlo.
In International Conference on Artificial Intelligence and
Statistics, pp.  7799–7813. PMLR.
- Horowitz (1991)
Horowitz, A. M. (1991, October).
A generalized guided Monte Carlo algorithm.
Physics Letters B 268(2), 247–252.
- Horowitz and
Lukic (2025)
Horowitz, B. and Z. Lukic (2025, February).
Differentiable Cosmological Hydrodynamics for Field-Level
Inference and High Dimensional Parameter Constraints.
arXiv:2502.02294 [astro-ph].
- Houtekamer and
Mitchell (2005)
Houtekamer, P. L. and H. L. Mitchell (2005).
Ensemble kalman filtering.
Quarterly Journal of the Royal Meteorological Society: A journal
of the atmospheric sciences, applied meteorology and physical
oceanography 131(613), 3269–3289.
- Houtekamer and
Zhang (2016)
Houtekamer, P. L. and F. Zhang (2016).
Review of the ensemble kalman filter for atmospheric data
assimilation.
Monthly Weather Review 144(12), 4489–4532.
- Kantorovich (1960)
Kantorovich, L. V. (1960, July).
Mathematical Methods of Organizing and Planning Production.
Management Science 6(4), 366–422.
Publisher: INFORMS.
- Lao and Louf (2022)
Lao, J. and R. Louf (2022).
Blackjax: Library of samplers for jax.
Astrophysics Source Code Library, ascl–2211.
- Leimkuhler and
Matthews (2015a)
Leimkuhler, B. and C. Matthews (2015a).
Molecular dynamics.
Interdisciplinary applied mathematics 36.
- Leimkuhler and
Matthews (2015b)
Leimkuhler, B. and C. Matthews (2015b).
Molecular Dynamics: With Deterministic and Stochastic
Numerical Methods, Volume 39 of Interdisciplinary Applied
Mathematics.
Cham: Springer International Publishing.
- Lüscher (2018)
Lüscher, M. (2018).
Stochastic locality and master-field simulations of very large
lattices.
EPJ Web Conf. 175, 01002.
_eprint: 1707.09758.
- Margossian (2019)
Margossian, C. C. (2019).
A review of automatic differentiation and its efficient
implementation.
WIREs Data Mining and Knowledge Discovery 9(4),
e1305.
_eprint: https://onlinelibrary.wiley.com/doi/pdf/10.1002/widm.1305.
- Margossian and
Gelman (2024)
Margossian, C. C. and A. Gelman (2024, February).
For how many iterations should we run Markov chain Monte Carlo?
arXiv:2311.02726 [stat].
- Margossian et al. (2024)
Margossian, C. C., M. D. Hoffman, P. Sountsov, L. Riou-Durand, A. Vehtari, and
A. Gelman (2024, January).
Nested Rˆ: Assessing the Convergence of Markov Chain
Monte Carlo When Running Many Short Chains.
Bayesian Analysis -1(-1), 1–28.
Publisher: International Society for Bayesian Analysis.
- Metropolis et al. (1953)
Metropolis, N., A. W. Rosenbluth, M. N. Rosenbluth, A. H. Teller, and E. Teller
(1953, June).
Equation of State Calculations by Fast Computing Machines.
The Journal of Chemical Physics 21(6), 1087–1092.
- Millea and
Seljak (2022)
Millea, M. and U. Seljak (2022).
Marginal unbiased score expansion and application to cmb lensing.
Physical Review D 105(10), 103531.
- Monmarché (2022)
Monmarché, P. (2022).
Hmc and underdamped langevin united in the unadjusted convex smooth
case.
arXiv preprint arXiv:2202.00977.
- Neal (2012)
Neal, R. M. (2012).
Bayesian learning for neural networks, Volume 118.
Springer Science & Business Media.
- Neal et al. (2011a)
Neal, R. M. et al. (2011a).
Mcmc using hamiltonian dynamics.
Handbook of markov chain monte carlo 2(11), 2.
- Neal
et al. (2011b)
Neal, R. M. et al. (2011b).
Mcmc using hamiltonian dynamics.
Handbook of markov chain monte carlo 2(11), 2.
- Nesterov (2009)
Nesterov, Y. (2009).
Primal-dual subgradient methods for convex problems.
Mathematical programming 120(1), 221–259.
- Olkin and
Pukelsheim (1982)
Olkin, I. and F. Pukelsheim (1982).
The distance between two random vectors with given dispersion
matrices.
Linear Algebra and its Applications 48, 257–263.
- Phan
et al. (2019)
Phan, D., N. Pradhan, and M. Jankowiak (2019).
Composable effects for flexible and accelerated probabilistic
programming in numpyro.
arXiv preprint arXiv:1912.11554.
- Riou-Durand and
Vogrinc (2023)
Riou-Durand, L. and J. Vogrinc (2023, December).
Metropolis Adjusted Langevin Trajectories: a robust alternative
to Hamiltonian Monte Carlo.
arXiv:2202.13230 [math, stat].
- Robnik
et al. (2025)
Robnik, J., R. Cohn-Gordon, and U. Seljak (2025, March).
Metropolis Adjusted Microcanonical Hamiltonian Monte Carlo.
arXiv:2503.01707 [stat].
- Robnik et al. (2022)
Robnik, J., G. B. De Luca, E. Silverstein, and U. Seljak (2022).
Microcanonical hamiltonian monte carlo.
arXiv preprint arXiv:2212.08549.
- Robnik et al. (2024)
Robnik, J., G. B. De Luca, E. Silverstein, and U. Seljak (2024, March).
Microcanonical Hamiltonian Monte Carlo.
The Journal of Machine Learning Research 24(1),
311:14696–311:14729.
- Robnik and
Seljak (2024)
Robnik, J. and U. Seljak (2024, July).
Fluctuation without dissipation: Microcanonical Langevin Monte
Carlo.
In Proceedings of the 6th Symposium on Advances in
Approximate Bayesian Inference, pp.  111–126. PMLR.
ISSN: 2640-3498.
- Rossky
et al. (1978)
Rossky, P. J., J. D. Doll, and H. L. Friedman (1978).
Brownian dynamics as smart monte carlo simulation.
The Journal of Chemical Physics 69(10), 4628–4633.
- Ruiz-Zapatero
et al. (2024)
Ruiz-Zapatero, J., D. Alonso, C. García-García, A. Nicola, A. Mootoovaloo,
J. M. Sullivan, M. Bonici, and P. G. Ferreira (2024, February).
LimberJack.jl: auto-differentiable methods for angular power
spectra analyses.
The Open Journal of Astrophysics 7,
10.21105/astro.2310.08306.
arXiv:2310.08306 [astro-ph].
- Simon-Onfroy et al. (2025)
Simon-Onfroy, H., F. Lanusse, and A. d. Mattia (2025, April).
Benchmarking field-level cosmological inference from galaxy redshift
surveys.
arXiv:2504.20130 [astro-ph].
- Sountsov
et al. (2020)
Sountsov, P., A. Radul, and contributors (2020).
Inference gym.
- Tuckerman et al. (2001)
Tuckerman, M. E., Y. Liu, G. Ciccotti, and G. J. Martyna (2001).
Non-hamiltonian molecular dynamics: Generalizing hamiltonian phase
space principles to non-hamiltonian systems.
The Journal of Chemical Physics 115(4), 1678–1702.
- Ver Steeg and
Galstyan (2021)
Ver Steeg, G. and A. Galstyan (2021).
Hamiltonian dynamics with non-newtonian momentum for rapid sampling.
Advances in Neural Information Processing Systems 34,
11012–11025.
- Villani (2003)
Villani, C. (2003).
Topics in optimal transportation.
Graduate studies in mathematics. Providence, RI: American
Mathematical Society.
OCLC: 51477002.
- (59)
Welling, M. and Y. W. Teh.
Bayesian Learning via Stochastic Gradient Langevin
Dynamics.
- Štrumbelj
et al. (2024)
Štrumbelj, E., A. Bouchard-Côté, J. Corander, A. Gelman, H. Rue, L. Murray,
H. Pesonen, M. Plummer, and A. Vehtari (2024, February).
Past, Present and Future of Software for Bayesian
Inference.
Statistical Science 39(1), 46–61.
Publisher: Institute of Mathematical Statistics.

## Appendix A Covariance matrix error

Let SPD($d$) be the space of symmetric positive definite matrices of size $d\times d$.

###### Lemma A.1 .

$b^{2}_{\mathit{cov}}(A,B)=\frac{1}{d}\Tr{(I-A^{-1}B)^{2}}$ is a divergence on SPD($d$), meaning that for all $A,B\in$ SPD($d$),

- 1.
$b^{2}_{\mathit{cov}}(A,B)\geq 0$
- 2.
$b^{2}_{\mathit{cov}}(A,B)=0$ if and only if $A=B$.

###### Proof.

Without loss of generality, we may assume that $A$ is diagonal with positive entries, because the trace is invariant under the change of basis and $A$ must be diagonal in some basis because it is positive-definite. $R\equiv I-A^{-1}B$ then has $R_{ii}=1-A^{-1}_{ii}B_{ii}$ on the diagonal and $R_{ij}=-A^{-1}_{ii}B_{ij}$ off the diagonal. The trace is

$$ $\Tr{R^{2}}=\sum_{i=1}^{d}(1-A^{-1}_{ii}B_{ii})^{2}+\sum_{i\neq j}A^{-1}_{ii}B_ {ij}A^{-1}_{jj}B_{ji},$ $$

where the first and the second term are the contribution from the diagonal and off-diagonal elements respectively.
Both terms are non-negative: the first term because it is a sum of squares, the second because all factors $A_{ii}$, $A_{jj}$ and $B_{ij}B_{ji}=(B_{ij})^{2}$ are non-negative.
This already proves (1).

The implication $b^{2}_{\mathit{cov}}(A,A)=0$ in (2) is trivial. To prove the other implication in (2), now suppose $b^{2}_{\mathit{cov}}(A,B)=0$. This implies that both terms in the above equation are zero, as they are both non-negative. In the second term, $A_{ii}>0$ are non-zero for any $i$, so the term can only be zero if $B_{ij}=0$ for any $i\neq j$. The first term can only be non-zero if $A_{ii}=B_{ii}$ for any $i$. We have shown that $A_{ij}=B_{ij}$ for any $i,j$, thus $A=B$.
∎

Note that the above result is not obvious from the fact that we are computing a trace of the matrix squared. There are non-zero matrices whose square is zero.

$b^{2}_{\mathit{cov}}$ has another nice property – it can be related to the effective sample size in a covariance matrix independent way:

###### Lemma A.2 .

Let $\boldsymbol{x}^{(k)}$ for $k=1,2\ldots n$ be exact i.i.d samples from $p=\mathcal{N}(0,\Sigma)$ and let
$\bar{\Sigma}=\frac{1}{n}\sum_{k=1}^{n}\boldsymbol{x}^{(k)}(\boldsymbol{x}^{(k)
})^{T}$,
be the empirical estimate for $\Sigma$. Then

$$ $\mathbb{E}_{MC}[b^{2}_{\mathit{cov}}(\Sigma,\bar{\Sigma})]=(d+1)/n,$ $$

where $\mathbb{E}_{MC}[\cdot]$ is the expectation with respect to the sample realizations.

###### Proof.

The empirical estimate is unbiased:

$$ $\mathbb{E}_{MC}[\bar{\Sigma}]=\frac{1}{n}\sum_{k=1}^{n}\mathbb{E}_{MC}[ \boldsymbol{x}^{(k)}(\boldsymbol{x}^{(k)})^{T}]=\Sigma,$ $$

and

$$ $\mathbb{E}_{MC}[\bar{\Sigma}_{ab}\bar{\Sigma}_{cd}]=\Sigma_{ab}\Sigma_{cd}+ \frac{1}{n}\big{(}\Sigma_{ac}\Sigma_{bd}+\Sigma_{ad}\Sigma_{bc}\big{)}.$ $$

The expected error of the empirical covariance matrix is then:

$$ $\displaystyle\mathbb{E}_{MC}[b^{2}_{\mathit{cov}}]=1-2\frac{1}{d}[\Sigma^{-1}] _{ij}\mathbb{E}_{MC}[\bar{\Sigma}_{ij}]+\frac{1}{d}[\Sigma^{-1}]_{ij}\mathbb{E }_{MC}[\bar{\Sigma}_{jk}\bar{\Sigma}_{li}][\Sigma^{-1}]_{kl}$ $\displaystyle=1-2+\frac{1}{d}[\Sigma^{-1}]_{ij}\bigg{(}\Sigma_{jk}\Sigma_{li}+ \frac{1}{n}\big{(}\Sigma_{jl}\Sigma_{ki}+\Sigma_{ji}\Sigma_{kl}\big{)}\bigg{)} [\Sigma^{-1}]_{kl}=\frac{d+1}{n},$ $$

where Einstein’s convention of summing over the repeated indices was used.
∎

Note that we could also take the more transparently non-negative definition $b_{F}^{2}(A,B)=\frac{1}{d}\Tr{(I-A^{-1}B)((I-A^{-1}B)^{T}}=\frac{1}{d}\norm{I-
A^{-1}B}_{F}^{2}$, where $\norm{\cdot}_{F}$ is the Frobenius norm. However, $b_{F}^{2}$ cannot be related to the effective sample size in a covariance matrix independent way. Instead,

$$ $\displaystyle\mathbb{E}_{MC}[b_{F}^{2}]$ $\displaystyle=\frac{1}{d}\mathbb{E}_{MC}[\big{(}\delta_{ij}-[\Sigma^{-1}]_{ik} \bar{\Sigma}_{kj}\big{)}\big{(}\delta_{ij}-[\Sigma^{-1}]_{il}\bar{\Sigma}_{lj} \big{)}]$ (13) $\displaystyle=1-2\frac{1}{d}[\Sigma^{-1}]_{ij}\mathbb{E}_{MC}[\bar{\Sigma}_{ij }]+\frac{1}{d}[(\Sigma^{-1})^{2}]_{ij}\mathbb{E}_{MC}[\bar{\Sigma}^{2}]_{ji}= \frac{1}{n}\big{(}1+\frac{1}{d}\Tr{\Sigma}\Tr{\Sigma^{-1}}\big{)},$ $$

so we will note use this definition here.

Finally, $b_{\mathit{cov}}$ is invariant to the linear change of basis

###### Lemma A.3 .

For an invertible matrix $A$, and a change of basis $\boldsymbol{x}^{\prime}=A\boldsymbol{x}$, the covariance matrix error does not change, that is, $b_{\mathit{cov}}^{2}(\Sigma_{p(\boldsymbol{x})},\Sigma_{q(\boldsymbol{x})})=b_
{\mathit{cov}}^{2}(\Sigma_{p(\boldsymbol{x}^{\prime})},\Sigma_{q(\boldsymbol{x
}^{\prime})})$.

###### Proof.

The covariance matrix transforms as

$$ $[\Sigma_{p(\boldsymbol{x}^{\prime})}]_{ij}=\int x^{\prime}_{i}x^{\prime}_{j}p( \boldsymbol{x}^{\prime})d\boldsymbol{x}^{\prime}=[A\,\Sigma_{p(\boldsymbol{x}) }A^{T}]_{ij},$ $$

hence the covariance matrix bias is invariant:

$$ $\displaystyle b_{\mathit{cov}}^{2}(\Sigma_{p(\boldsymbol{x})},\Sigma_{q( \boldsymbol{x})})$ $\displaystyle=\frac{1}{d}\Tr{\big{(}I-(A\Sigma_{p(\boldsymbol{x})}A^{T})^{-1}( A\Sigma_{q(\boldsymbol{x})}A^{T})\big{)}^{2}}$ $\displaystyle=\frac{1}{d}\Tr{\big{(}I-(A^{T})^{-1}\Sigma_{p(\boldsymbol{x})}^{ -1}A^{-1})(A\Sigma_{q(\boldsymbol{x})}A^{T})\big{)}^{2}}$ $\displaystyle=\frac{1}{d}\Tr{\big{(}I-(A^{T})^{-1}\Sigma_{p(\boldsymbol{x})}^{ -1}\Sigma_{q(\boldsymbol{x})}A^{T}\big{)}^{2}}$ $\displaystyle=\frac{1}{d}\Tr{(A^{T})^{-1}\big{(}I-\Sigma_{p(\boldsymbol{x})}^{ -1}\Sigma_{q(\boldsymbol{x})}\big{)}^{2}A^{T}}$ $\displaystyle=b_{\mathit{cov}}^{2}(\Sigma_{p(\boldsymbol{x}^{\prime})},\Sigma_ {q(\boldsymbol{x}^{\prime})}).$ $$

∎

###### Lemma A.1 .

###### Proof.

###### Lemma A.2 .

###### Proof.

###### Lemma A.3 .

###### Proof.

## Appendix B Proofs

### B.1 Lemma 4.1 : EEVPD

###### Proof.

We will work in the eigenbasis, where the dynamics is decoupled.
It then suffices to analyze each dimension separately. Let $x_{i}(t)$ and $u_{i}(t)$ be components of $\boldsymbol{x}(t)$ and $\boldsymbol{u}(t)$ along the dimension that we analyze and $\sigma_{i}^{2}$ the eigenvalue of the covariance matrix in that direction. The velocity Verlet integrator update ([4](https://arxiv.org/html/2412.08876v2#S3.E4)) with step size $\epsilon$
can be written compactly as
(Gouraud
et al., 2023)

$$ $\Phi_{\epsilon}(\boldsymbol{z}_{i})=\begin{bmatrix}\cos h&\alpha\sigma_{i}\sin h \\ -\alpha^{-1}\sigma_{i}^{-1}\sin h&\cos h\end{bmatrix}\boldsymbol{z}_{i}\equiv A \boldsymbol{z}_{i}.$ $$

Here, $\boldsymbol{z}_{i}(t)=(x_{i}(t),u_{i}(t))$, $\alpha=(1-y_{i}/4)^{-1/2}$, $y_{i}=\epsilon^{2}/\sigma_{i}^{2}$ and $\sin h=\sqrt{y_{i}}/\alpha$. The energy error is

$$ $\Delta_{H}^{i}\equiv\Delta_{H}(\Phi_{\epsilon}(\boldsymbol{z}_{i}),\boldsymbol {z}_{i})=\frac{1}{2}\Phi_{\epsilon}(\boldsymbol{z}_{i})^{T}D\Phi_{\epsilon}( \boldsymbol{z}_{i})-\frac{1}{2}\boldsymbol{z}_{i}^{T}D\boldsymbol{z}_{i}=\frac {1}{2}\boldsymbol{z}_{i}^{T}M\boldsymbol{z}_{i},$ $$

where $D=\mathrm{Diag}(1/\sigma_{i}^{2},1)$ and

$$ $M=A^{T}DA-D=(1-\alpha^{2})\begin{bmatrix}(\alpha\sigma_{i})^{-2}\sin^{2}h&-( \alpha\sigma_{i})^{-1}\sin h\cos h\\ -(\alpha\sigma_{i})^{-1}\sin h\cos h&-\sin^{2}h\\ \end{bmatrix}.$ $$

Denote by $\tilde{p}_{i}$ the stationary distribution for $\boldsymbol{z}_{i}$, namely
$x_{i}\sim\mathcal{N}(0,\tilde{\sigma}_{i}),\,u_{i}\sim\mathcal{N}(0,1)$, where $\tilde{\sigma}_{i}=\sigma_{i}/\alpha$ is taken from Equation ([11](https://arxiv.org/html/2412.08876v2#S4.E11)).

The contribution to the $\mathrm{EEVPD}$ from $\boldsymbol{z}_{i}$ is

$$ $\mathrm{Var}_{\tilde{p}_{i}}[\Delta_{H}^{i}]=\mathbb{E}_{\tilde{p}_{i}}[( \Delta_{H}^{i})^{2}]-\mathbb{E}_{\tilde{p}}[\Delta_{H}^{i}]^{2}.$ $$

The expectation value in the second term is

$$ $\mathbb{E}_{\tilde{p}}[\Delta_{H}^{i}]=\frac{1}{2}\big{(}M_{11}\mathbb{E}_{ \tilde{p}}[x^{2}]+M_{22}\mathbb{E}_{\tilde{p}}[u^{2}]\big{)}$ $$

and for the first

$$ $\displaystyle\mathbb{E}_{\tilde{p}}[(\Delta_{H}^{i})^{2}]$ $\displaystyle=\frac{1}{4}\mathbb{E}_{\tilde{p}}[(M_{11}x^{2}+2M_{12}xu+M_{22}u ^{2})^{2}]$ $\displaystyle=\frac{1}{4}\big{(}M_{11}^{2}\mathbb{E}_{\tilde{p}}[x^{4}]+M_{22} ^{2}\mathbb{E}_{\tilde{p}}[u^{4}]+(4M_{12}^{2}+2M_{11}M_{22})\mathbb{E}_{ \tilde{p}}[x^{2}]\mathbb{E}_{\tilde{p}}[u^{2}]\big{)}.$ $$

In both expressions we have dropped the vanishing contributions that contain $\mathbb{E}_{\tilde{p}}[xu]=0$, $\mathbb{E}_{\tilde{p}}[xu^{3}]=0$ or $\mathbb{E}_{\tilde{p}}[x^{3}u]=0$.
Combining both terms together and using $\mathbb{E}_{\tilde{p}}[x^{4}]=3\mathbb{E}_{\tilde{p}}[x^{2}]^{2}$ and $\mathbb{E}_{\tilde{p}}[u^{4}]=3\mathbb{E}_{\tilde{p}}[u^{2}]^{2}$ we get

$$ $\mathrm{Var}[\Delta_{H}^{i}]=\frac{1}{4}\big{(}2M_{11}^{2}\mathbb{E}_{\tilde{p }}[x^{2}]^{2}+2M_{22}^{2}\mathbb{E}_{\tilde{p}}[u^{2}]^{2}+4M_{12}^{2}\mathbb{ E}_{\tilde{p}}[x^{2}u^{2}]\big{)}.$ $$

Inserting
$\mathbb{E}_{\tilde{p}}[x^{2}]=\sigma^{2}\alpha^{2}$ and
$\mathbb{E}_{\tilde{p}}[u^{2}]=1$ gives

$$ $\displaystyle\mathrm{Var}[\Delta_{H}^{i}]=\frac{(1-\alpha^{2})^{2}}{2}\big{(} \sin^{4}h+\sin^{4}h+2\sin^{2}h\cos^{2}h\big{)}=(1-\alpha^{2})^{2}\sin^{2}h= \frac{y^{3}}{16(1-y/4)},$ $$

which is $E(y)$ from the statement of the theorem. $\mathrm{EEVPD}$ is thus

$$ $\mathrm{EEVPD}=\frac{1}{d}\sum_{i=1}^{d}\mathrm{Var}[\Delta_{H}^{i}]=\frac{1}{ d}\sum_{i=1}^{d}E(y_{i}).$ $$

∎

###### Proof.

### B.2 Theorem 4.2 : bias bounds

###### Proof.

Let’s start with the covariance matrix bias.
Due to Equation ([11](https://arxiv.org/html/2412.08876v2#S4.E11)), the asymptotic covariance error ([6](https://arxiv.org/html/2412.08876v2#S3.E6)) is

$$ $b^{2}_{\mathit{cov}}(\Sigma,\tilde{\Sigma})=\frac{1}{d}\sum_{i=1}^{d}\big{(}1- (1-\epsilon^{2}/4\sigma_{i}^{2})^{-1}\big{)}^{2}=\frac{1}{d}\sum_{i=1}^{d}B( \epsilon^{2}/\sigma_{i}^{2}),$ $$

where
$B(y)=\frac{y^{2}}{16(1-y/4)^{2}}$.
We thus see that $\varphi$ from the statement of the theorem is $\varphi=E\circ B^{-1}$.
Lemma [B.1](https://arxiv.org/html/2412.08876v2#A2.Thmtheorem1) shows that $\varphi(x)$ restricted to $0<x<11-4\sqrt{7}$ is a convex, monotonically increasing function. By Jensen’s inequality this implies that

$$ $\varphi(b_{\mathit{cov}}(\Sigma,\tilde{\Sigma})^{2})\leq\mathrm{EEVPD},$ $$

as long as $\varphi(b^{2}_{\mathit{cov}})<\varphi(11-4\sqrt{7})=(-134+52\sqrt{7})/9\approx
0
.397674$. The assumption of the theorem that $\mathrm{EEVPD}<0.397$ is a sufficient condition for this to hold.
Since $\varphi(x)$ is not a linear function, Jensen’s inequality becomes an equality if and only if $\sigma_{i}=\sigma_{j}$ for all $i,j$.
$\varphi$ is a monotonically increasing function so it is a bijection and its inverse is also a monotonically increasing function. The inverse can therefore be applied to both sides of the above inequality to yield the desired result.

The proof of the Wasserstein distance part of the theorem is similar. For zero-mean Gaussian distributions with covariance matrices $\Sigma_{p}$ and $\Sigma_{q}$, the Wasserstein distance reduces to (Olkin and
Pukelsheim, 1982)

$$ $\mathcal{W}_{2}(\mathcal{N}(0,\Sigma_{p}),\mathcal{N}(0,\Sigma_{q}))^{2}=\Tr{ \Sigma_{p}+\Sigma_{q}-2\big{(}\Sigma_{p}^{1/2}\Sigma_{q}\Sigma_{p}^{1/2}\big{) }^{1/2}}.$ $$

By using Equation ([11](https://arxiv.org/html/2412.08876v2#S4.E11)) for unadjusted HMC or LMC with the velocity Verlet integrator we therefore get

$$ $\mathcal{W}_{2}(p,\tilde{p})^{2}=\epsilon^{2}\sum_{i=1}^{d}W(\epsilon^{2}/ \sigma_{i}^{2}),$ $$

where

$$ $W(y)=\frac{2(1-y/8-\sqrt{1-y/4})}{y(1-y/4)}$ $$

is as in the statement of the theorem.
Lemma [B.2](https://arxiv.org/html/2412.08876v2#A2.Thmtheorem2) shows that $\varphi_{W}(x)$ is a convex, monotonically increasing function for $0<\varphi_{W}(x)<27/4=6.75$. By Jensen’s inequality this implies that

$$ $\varphi_{W}(w)\leq\mathrm{EEVPD},$ $$

as long as $\mathrm{EEVPD}<6.75$. Here $w=\mathcal{W}_{2}(p,\tilde{p})^{2}/d\epsilon^{2}$.
Since $\varphi_{W}$ is not a linear function, Jensen’s inequality becomes an equality if and only if $\sigma_{i}=\sigma_{j}$ for all $i,j$.
$\varphi_{W}$ is a monotonically increasing function so it is a bijection and its inverse is also a monotonically increasing function. The inverse can therefore be applied to both sides of the above inequality to yield

$$ $w\leq\varphi_{W}^{-1}(\mathrm{EEVPD}).$ $$

∎

###### Proof.

### B.3 Convexity

###### Lemma B.1 .

$\varphi(x)=4x^{3/2}/(1+x^{1/2})^{2}$ is monotonically increasing for $x>0$ and convex for $0<x<11-4\sqrt{7}$.

###### Proof.

$\varphi(x)$ is monotonically increasing because its derivative

$$ $\varphi^{\prime}(x)=\frac{2x^{1/2}(3+x^{1/2})}{(1+x^{1/2})^{3}}$ $$

is positive (all terms are positive).
To show that it is convex, we compute its second derivative

$$ $\varphi^{\prime\prime}(x)=\frac{3-4x^{1/2}-x}{x^{1/2}(1+x^{1/2})^{4}}$ $$

The denominator is positive for $x>0$. The numerator is a quadratic polynomial $p(y)=3-4-y^{2}$ in $y=\sqrt{x}$. Its roots are $y_{1,2}=-2\pm\sqrt{7}$. Since $p(0)=3>0$ the numerator is positive for $0<y<-2+\sqrt{7}$ corresponding to $0<x<(-2+\sqrt{7})^{2}=11-4\sqrt{7}$ so $\varphi(x)$ restricted to this interval is convex.
∎

###### Lemma B.2 .

$\varphi_{W}(x)$ from Theorem [4.2](https://arxiv.org/html/2412.08876v2#S4.Thmtheorem2) is monotonically increasing and convex for $0<\varphi_{W}(x)<27/4$.

###### Proof.

To prove that $\varphi_{W}$ is monotonically increasing we will show that its derivative is positive. We cannot solve for for $W^{-1}$ explicitly, but nontheless

$$ $\varphi_{W}^{\prime}(w)=E^{\prime}(W^{-1}(w))(W^{-1})^{\prime}(w)=\frac{E^{ \prime}(y)}{W^{\prime}(y)},$ $$

where we have denoted $y=W^{-1}(w)$. We will show that both denominator and numerator are positive for $0\leq y<4$, i.e. in the range where the velocity Verlet integrator is stable.

The numerator is

$$ $E^{\prime}(y)=\frac{3y^{2}(1-y/6)}{(1-y/4)^{2}}$ $$

which is positive for $0<y<4$.

To simplify the denominator we use the reparametrization $y=4\sin^{2}\xi$, such that $0\leq\xi<\pi/2$. In the new parametrization

$$ $W(\xi(y))=\frac{\sin^{4}(\xi/2)}{4\sin^{2}(2\xi)}.$ $$

We get

$$ $W^{\prime}(y)=\frac{1+\cos\xi-\cos^{2}\xi}{128\cos^{4}(\xi/2)\cos^{4}(\xi)} \geq\frac{\cos\xi}{128\cos^{4}(\xi/2)\cos^{4}(\xi)}>0,$ $$

where we have used that $-\cos^{2}\xi>-1$.

To prove that $\varphi_{W}$ is convex, we will show that its second derivative is positive. We have

$$ $\varphi_{W}^{\prime\prime}(w)=E^{\prime\prime}(W^{-1}(w))((W^{-1})^{\prime}(w) )^{2}+E^{\prime}(W^{-1}(w))(W^{-1})^{\prime\prime}(w)=\frac{E^{\prime\prime}(y )}{W^{\prime}(y)^{2}}-\frac{E^{\prime}(y)W^{\prime\prime}(y)}{W^{\prime}(y)^{3}}$ $$

Which after some algebra reduces to

$$ $\varphi_{W}^{\prime\prime}(y)=\frac{(2-s)(s^{2}-1)(s^{4}+4s^{3}+7s^{2}+6s-6)}{ 8(s^{2}+s-1)},$ $$

where we have first reparametrized $y(\xi)=4\sin^{2}\xi$ and then $s(\xi)=1/\cos(\xi)$.
The range $0<y<4$ corresponds to $0<\xi<\pi/2$ which in turn corresponds to $1<s<\infty$.

All the factors except for $(2-s)$ are non-negative at $s=1$ and monotonically increasing for $s>1$, so they are all positive for $s>1$. The second derivative is then positive for $1<s<2$, corresponding to $0<\xi<\pi/3$ or $0<y<3$.
$\varphi_{W}(x)$ is thus monotonically increasing and convex for $0<\varphi_{W}(x)<E(3)=27/4$.
∎

###### Lemma B.1 .

###### Proof.

###### Lemma B.2 .

###### Proof.

## Appendix C Bias-variance tradeoff

We have shown how to control the covariance matrix bias to be below some desired threshold. However, typically bias is not of direct interest, but instead we want to control the error of the expectation values, which additionally contains the variance, i.e. the second term in Equation ([2](https://arxiv.org/html/2412.08876v2#S1.E2)). The variance of a stationary chain is

$$ $\mathrm{Var}[\bar{f}]=\mathrm{Var}_{p}[f](1+2\sum_{k=1}^{n}(1-k/n)\rho_{k})/n \equiv\mathrm{Var}_{p}[f]\tau_{\mathrm{int}}/n,$ (14) $$

where the autocorrelation coefficients are $\rho_{k}=\mathbb{E}[(f(x_{i})-\mathbb{E}_{p}[f])(f(x_{i+k})-\mathbb{E}_{p}[f])
]/\mathrm{Var}_{p}[f]$.

We would therefore like to know how to optimally set the bias, given some error tolerance. Here we will provide a heuristic, based on estimating the second moment $\mathbb{E}[x^{2}]$ of a one-dimensional standard Gaussian target with velocity Verlet unadjusted HMC.
In this case, $\rho_{k}=\rho^{k}$ (Gouraud
et al., 2023), where

$$ $\rho=\cos^{2}\big{(}\frac{T}{\epsilon}\arcsin(\alpha\epsilon/\sigma)\big{)}.$ (15) $$

This makes the sum in Equation ([14](https://arxiv.org/html/2412.08876v2#A3.E14)) expressable in terms of geometric series:

$$ $\tau_{\mathrm{int}}=1+2S(\rho)-2\rho S^{\prime}(\rho)/n=\frac{1+\rho}{1-\rho} \bigg{(}1-\frac{2\rho}{n}\frac{1-\rho^{n}}{1-\rho^{2}}\bigg{)},$ (16) $$

where $S(\rho)=\sum_{k=1}^{n}\rho^{k}=\rho(1-\rho^{n})/(1-\rho)$. The variance at the lowest order in the small-$\epsilon$-expansion is then

$$ $\mathrm{Var}[x^{2}]\asymp\frac{2L}{N\epsilon}\lim_{\epsilon\xrightarrow[]{}0} \tau_{\mathrm{int}}(\epsilon)=c_{v}/\epsilon,$ (17) $$

where $c_{v}$ is constant, independent of the step size and $\asymp$ denotes asymptotic equivalence as $\epsilon\xrightarrow[]{}0$.
The bias of the second moment in this limit is

$$ $\mathrm{Bias}[x^{2}]=\frac{1}{1-\epsilon^{2}/4}-1\asymp\epsilon^{2}/4=c_{b} \epsilon^{2},$ (18) $$

where we have used Equation ([11](https://arxiv.org/html/2412.08876v2#S4.E11)) and denoted $c_{b}=1/4$.

Combining Equations ([17](https://arxiv.org/html/2412.08876v2#A3.E17)) and ([18](https://arxiv.org/html/2412.08876v2#A3.E18)) we get for the RMSE in the small step size limit:

$$ $\mathrm{RMSE}^{2}=c_{b}\epsilon^{4}+c_{v}/\epsilon,$ (19) $$

where $c_{b}$ and $c_{v}$ are defined above and are independent of the step size.
We would like to set the step size so that it minimizes the RMSE. The optimum is found at

$$ $0=\frac{d}{d\epsilon}\mathrm{RMSE}^{2}=4c_{b}\epsilon^{3}-c_{v}/\epsilon^{2},$ (20) $$

which gives the optimal step size $\epsilon_{\mathrm{opt}}=(c_{v}/4c_{b})^{1/5}$. At the optimal step size, bias squared is one fifth of the error squared:

$$ $\frac{\mathrm{Bias}^{2}}{\mathrm{RMSE}^{2}}=\frac{1}{1+\frac{c_{v}/\epsilon_{ \mathrm{opt}}}{c_{v}\epsilon_{\mathrm{opt}}^{4}}}=\frac{1}{5}.$ (21) $$

which is the prescription we use in Table [1](https://arxiv.org/html/2412.08876v2#S6.T1).

## Appendix D Experiment Details

### D.1 Benchmark Problems

The following benchmarks are used:

- •
A Standard Gaussian in $d=100$.
- •
An ill-conditioned Gaussian in $d=100$ and condition number $\kappa=1000$. The eigenvalues of the covariance matrix are equally spaced in log.
- •
A Rosenbrock function with $Q=0.1$ from Grumitt
et al. (2022). This is a banana shaped target in two dimensions, see Figure 8 in Robnik et al. (2022). We use a product of 18 independent copies, so the total dimension of the target is 36. An exception is Figure [2](https://arxiv.org/html/2412.08876v2#S1.F2) where we study the performance as a function of the number of copies.
- •
A Funnel problem in 101 dimensions: this is a hierarchical Bayesian model with a funnel shape (Grumitt
et al., 2022). The goal is to infer the hierarchical parameter $\theta$ and the latent variables $\{z_{i}\}_{i=1}^{100}$, given the noisy observations $y_{i}\sim\mathcal{N}(z_{i},1)$. The prior is Neal’s funnel (Neal
et al., 2011b):
$\theta\sim\mathcal{N}(0,3)$,
$z_{i}\sim\mathcal{N}(0,e^{\theta/2})$.
We set $\theta_{\mathrm{true}}=0$ and generate the data with the generative process described above. Given this data we then sample from the posterior for $\theta$ and $\{z_{i}\}_{i=1}^{100}$.
- •
A Brownian motion example from the Inference Gym (Sountsov
et al., 2020), where it is named BrownianMotionUnknownScalesMissingMiddleObservations.
This is a 32 dimensional hierarchical Bayesian model where Brownian motion with unknown innovation noise and measurement noise is fitted to the noisy and partially missing data.
- •
The German Credit model, also known as Sparse logistic regression (GermanCreditNumericSparseLogisticRegression) is a 51-dimensional Bayesian hierarchical model, where logistic regression is used to model the approval of the credit based on the information about the applicant.
- •
An Item Response theory model (SyntheticItemResponseTheory), which is a 501-dimensional hierarchical problem where students’ ability is inferred, given the test results.
- •
Stochastic Volatility is a 2429-dimensional hierarchical non-Gaussian random walk fit to the S&P500 returns data, adapted from NumPyro (Phan
et al., 2019).

The ground truth covariance matrix for the first two problems is known exactly. For the Rosenbrock function, we compute it by drawing exact samples from the posterior. For the other problems, we obtain the ground truth by running very long NUTS chains.

### D.2 Step size adaptation

Figure: Figure 5: The step size adaptation algorithm from Section [D.2](https://arxiv.org/html/2412.08876v2#A4.SS2), applied to the Rosenbrock target distribution in $d=36$ with MCLMC sampler. The sequential algorithm was initialized from the standard Gaussian distribution with a random initial velocity orienation. Top: the step size as a function of leapfrog integration steps. Bottom: per dimension squared energy error for each step. The algorithm quickly converges to the targeted EEVPD = 0.001, shown with a black line.
Refer to caption: x4.png

We now outline a scheme for quickly adapting the step size $\epsilon$ to achieve a desired EEVPD. We verified that it indeed yields the desired EEVPD in the experiments performed in this work, but we do not provide any convergence guarantees. More generally, one can instead use a more established dual averaging algorithm (Hoffman
et al., 2014; Nesterov, 2009).

Suppose we did a leapfrog step with size $\epsilon_{k}$ and found some energy error $\Delta_{H}^{k}$. Using only this knowledge and the scaling from Equation ([4.1](https://arxiv.org/html/2412.08876v2#S4.Thmtheorem1)) for small step sizes, $\Delta_{H}\propto\epsilon^{6}$, we could estimate the optimal step size as $\epsilon=\xi_{k}^{-1/6}$ where

$$ $\xi_{k}=\frac{(\Delta_{H}^{k})^{2}}{d}\frac{1}{\alpha\,\epsilon_{k}^{6}}$ (22) $$

and $\alpha$ is the desired EEVPD. As we do more leapfrog steps, we can improve our estimate by averaging energy errors and use the predicted optimal step size in the next step:

$$ $\epsilon_{n+1}=\bigg{(}\frac{\sum_{k=1}^{n}w(\xi_{k})\gamma^{n-k}\,\xi_{k}}{ \sum_{k=1}^{n}w(\xi_{k})\gamma^{n-k}}\bigg{)}^{-1/6}.$ (23) $$

We have introduced two types of weights:

- •
The weights $w$ parametrize our trust in the predictions from the too large and too small $\epsilon$. We take the log-normal penalty
$w(\xi)=\exp{-\frac{1}{2}(\log\xi)^{2}/\sigma_{\xi}^{2}},$
(24)
with $\sigma_{\xi}=1.5$.
- •
$\gamma$ is the forgetting factor. It is related to the effective sample size $n$ of the estimate (if $w$ were constant) by $\gamma=\frac{n-1}{n+1}$. $n$ is also the number of steps after which the weights have decayed to $e^{-2}=0.13$. In general, we don’t want $n$ to be too small, so that EEVPD is well determined and yet not too large during the burn-in such that the initially heavily biased estimates are forgotten quickly.
We find $n=50$ to work well on all benchmark tests from (Robnik et al., 2022). An example run is shown in Figure [5](https://arxiv.org/html/2412.08876v2#A4.F5).

The pseudocode for the proposed algorithm is shown in [1](https://arxiv.org/html/2412.08876v2#algorithm1).

Figure: Algorithm 1 Step size adaptation

### D.3 Details of experiments

We use an run of NUTS to find a diagonal preconditioning matrix. In practice, we would use an unadjusted sampler rather than NUTS for this purpose, but we wish the preconditioning matrix to be identical for all problems considered, since it is not the object of our study.

We determine step size by running another chain. For unadjusted methods, we use the algorithm from [D.2](https://arxiv.org/html/2412.08876v2#A4.SS2) to select a step size that ensures a desired $\mathrm{EEVPD}$. For adjusted MCLMC we use the dual averaging algorithm (Nesterov, 2009) from Hoffman and
Gelman (2014) to adapt the step size to achieve an acceptance rate of $90\%$.

We take the tuning steps as our burn-in, initializing the chain with the final state returned by the tuning procedure. For each model we run at least 128 chains, and take the median of the error across chains at each step. This reduces the error in the quantities in Tables [2](https://arxiv.org/html/2412.08876v2#S7.T2) and [3](https://arxiv.org/html/2412.08876v2#S7.T3). The errors are shown in Table [4](https://arxiv.org/html/2412.08876v2#A4.T4). They are calculated by bootstrap: for a given model, we produce a set of chains (at least 128), and calculate the bias $b_{\mathit{avg}}$ or $b_{\mathit{cov}}$ at each step of the chain. We then resample (with replacement) 100 times from this set, and compute our final metric (number of gradients to low bias) 100 times. We take the standard deviation of this list of length 100 as the estimate of the error and report it relatively to the values in Tables [2](https://arxiv.org/html/2412.08876v2#S7.T2) and [3](https://arxiv.org/html/2412.08876v2#S7.T3).

**Table 4: Relative error associated with Tables [2](https://arxiv.org/html/2412.08876v2#S7.T2) and [3](https://arxiv.org/html/2412.08876v2#S7.T3).**
|  |  | aLMC | uLMC | aMCLMC | uMCLMC | NUTS |
| --- | --- | --- | --- | --- | --- | --- |
| Model | metric |  |  |  |  |  |
| Standard Gaussian | $b_{\mathit{avg}}$ | 0.42% | 1.02% | 1.51% | 0.77% | 0.13% |
| $b_{\mathit{cov}}$ | 0.06% | 0.32% | 0.31% | 0.29% | 0.06% |  |
| Rosenbrock | $b_{\mathit{avg}}$ | 0.07% | 5.69% | 0.06% | 2.05% | 0.02% |
| $b_{\mathit{cov}}$ | 0.03% | 1.64% | 0.01% | 0.90% | 0.01% |  |
| Brownian Motion | $b_{\mathit{avg}}$ | 0.37% | 6.14% | 0.38% | 2.77% | 0.16% |
| $b_{\mathit{cov}}$ | 0.05% | 4.96% | 0.06% | 0.49% | 0.03% |  |
| German Credit | $b_{\mathit{avg}}$ | 0.06% | 7.49% | 0.10% | 2.25% | 0.04% |
| $b_{\mathit{cov}}$ | 0.05% | 11.58% | 0.08% | 0.93% | 0.05% |  |
| Item Response | $b_{\mathit{avg}}$ | 0.62% | 8.08% | 1.23% | 4.11% | 0.44% |
| Stochastic Volatility | $b_{\mathit{avg}}$ | 0.04% | 8.88% | 0.06% | 1.87% | 0.02% |

NUTS is run using the BlackJax (Cabezas
et al., 2024) implementation, with the provided window adaptation scheme.
For adjusted LMC (MALT), we perform a search over different values of trajectory length, and at each, choose $\epsilon$ to target an acceptance rate of $0.8$.

## Appendix E Ablation study

Figure: Figure 6: Performance of unadjusted LMC as a function of EEVPD (which in turn sets the step size). As can be seen, performance does not change much within a resonable range of EEVPD. The value $\mathrm{EEVPD}=3\times 10^{-4}$ that we use in section [7](https://arxiv.org/html/2412.08876v2#S7) is shown as a vertical dashed line.
Refer to caption: extracted/6448359/ablation.png

Here, we investigate how performance of unadjusted LMC, with the tuning algorithm proposed in Section [6](https://arxiv.org/html/2412.08876v2#S6) varies with the $\mathrm{EEVPD}$ value being targeted. Figure [6](https://arxiv.org/html/2412.08876v2#A5.F6) shows that our choice of $3\times 10^{-4}$ is within the safe range for problems that we consider, albeit it is somewhat conservative for the Item Response and Rosenbrock problems.

## Appendix F Reproducibility

All the samplers considered and their adaptation schemes are implemented in blackjax (Lao and Louf, 2022):

[https://blackjax-devs.github.io/blackjax/](https://blackjax-devs.github.io/blackjax/).

The code for benchmarking and numerical experiments is available here:

[https://github.com/reubenharry/sampler-benchmarks](https://github.com/reubenharry/sampler-benchmarks).

## Peer Reviews

confidence: 3
contribution: 3
initial_rating: 8
presentation: 4
questions: N/A
soundness: 4
strengths: - The proposed criterion EEVPD is intuitive and simple to compute online. Furthermore it is backed with some theory in a Gaussian setting and for the more general case the authors show experimentally that a function of EEVPD upperbounds $b_cov$, suggesting some generality to the result proved in theorem 4.2
- The paper contains broad benchmarks with high‑dimensional case and the method is winning against NUTS/adjusted baselines at fixed error. 
- Overall I found the paper to be well-written, with clear theory and experiments. The paper is definitely interesting for anyone in the ICLR  community interested in Bayesian inference and MCMC algorithms.
summary: The paper proposes a black‑box tuning scheme for unadjusted, gradient‑based samplers (HMC, underdamped LMC...). The key idea is to control asymptotic discretization bias by monitoring the Energy‑Error Variance Per Dimension (EEVPD), i.e. the variance of the one‑step Hamiltonian error normalized by dimension. This is cheap to estimate (online) quantity from the integrator’s energy changes. For Gaussian targets, the authors prove sharp, dimension‑wise bounds connecting EEVPD to (i) a covariance‑error divergence and (ii) Wasserstein‑2 distance, and they use these bounds to set a target EEVPD and adapt the step size automatically.
weaknesses: The only weakness of this paper is the scope of the theoretical results. It would have certainly been better to have more general results (at least under log-concavity).

---

confidence: 4
contribution: 1
initial_rating: 2
presentation: 2
questions: 1. In multiple instances, such as on line 175, general statements are supported by one or two references. For example:
   *"Unadjusted methods have also been analyzed theoretically, establishing a bound on the asymptotic bias and the mixing time."*
   It is unclear how the cited works were selected. Are they the earliest papers on the topic, the most recent, or those providing the sharpest results?

2. Lines 88–90 state:
   *"Noting that chains are finite in practice, one could negotiate the tradeoff differently and reduce the variance at the cost of introducing some bias. This strategy is the focus of the present work."*
   Initially, I interpreted this as proposing a method to select the step-size by balancing bias and variance. However, I now suspect this is not the intended meaning. Is my understanding correct?
soundness: 2
strengths: The paper addresses an important, timely problem relevant to the ICLR audience: selecting the step size in MCMC algorithms. This challenge has broad applications across machine learning and computational statistics.

I ike the idea of monitoring the expected change of the Hamiltonian.
summary: The paper aims to develop a practical method for selecting the discretization step-size in Hamiltonian Monte Carlo (HMC) and Underdamped Langevin Monte Carlo (ULMC) algorithms. The theoretical analysis is restricted to basic computations in a Gaussian setting. Based on these computations, an algorithm is proposed—though its derivation lacks precision—and its performance is evaluated through numerical simulations on selected examples.
weaknesses: ### 1. Theoretical Contributions
The theoretical contributions of the paper are simple and of limited interest:
The justification for minimizing EEVPD is unconvincing:
  - There is no clear evidence that EEVPD can be easily approximated in practice.
  - In non-Gaussian settings, there is no reason to assume that EEVPD controls the Wasserstein error.
  - Even if a function $\phi$ exists such that an EEVPD smaller than $\varepsilon$ implies a Wasserstein error smaller than $\phi(\varepsilon)$, the form of $\phi$ would likely depend on the target distribution, limiting the generality of the result.

---

### 2. Methodological Contributions
The main contributions of the paper are methodological, yet they lack clarity and practical guidance:
- A pseudo-code for the proposed method is missing. Section 6 provides only a high-level description, which is insufficient for understanding how the method should be implemented.

---

### 3. Unjustified Claims
There are several strong and unjustified claims; some of them are listed below:
- **Line 67:** *"In statistics, MH-adjusted versions of these dynamics are used almost exclusively."*
  This is an overgeneralization and lacks nuance.
- **Line 77:** *"Thus, unadjusted methods provably outperform adjusted methods."*
  
- **Line 79:** *"Without a principled way to do so, general practitioners cannot use unadjusted algorithms in a black-box fashion, which gravely limits their applicability."*
  This claim is inconsistent with practice: practitioners routinely use stochastic gradient descent (SGD) and its variants in deep learning, despite the lack of a universally accepted principled method for selecting the learning rate.

---

### 4. Major Weakness: Practical Estimability
A critical weakness of the paper is its reliance on the following statement from **line 253**:
*"Crucially, this is a quantity that can easily be estimated in practice."*
The quantity referred to involves an expectation with respect to the stationary distribution. For arbitrary target distributions, it is unclear when stationarity is achieved or how closely the process approaches stationarity after a finite number of steps. As a result, the claim that this quantity can be easily estimated lacks justification.

---

### 5. The fact that the initialization error is disregarded, as written on line 80, is an important limitation.

---

confidence: 3
contribution: 3
initial_rating: 6
presentation: 3
questions: --In order to better understand the strenth/weakness of this method in non-Gaussian distributions, can you provide the numerical results on some multi-modal distributions that are much different from Gaussian dsitributions (e.g., the Gaussian mixture distribution where each Gaussian distribution is far away from each other but has similar variance)?
soundness: 3
strengths: 1. The proposed method has strong theoretical fundation for Gaussian distribution. It also provides empirical results to show the effectiveness on non-Gaussian distributions.
2. The numerical results show that the proposed approach has achieved near-optimal RMSE with orders-of-magnitude fewer gradient evaluations than adjusted (MH) methods. It also outperform NUTS and manually tuned adjusted LMC/HMC.
3. Open-source codes are provided.
summary: This paper studies how to adjust the integration stepsize automatically to control the bias in HMC and Langevin Monte Carlo while without using the Metropolis-Hastings (MH) correction. The key idea behind this paper is that the energy error during numerical integration directly controls the asymptotic bias. As a result, one can adaptively selects the step size to guarantees a bias tolerance by monitoring the energy error variance per dimension (EEVPD). 

Specific contributions include:
1. Analytical results for Gaussian distribution, showing the relationship betwen EEVPD and integration stepsize (Lemma 4.1) as well as the associated covariance matrix bias.
2. Numerical impelmentation of the automatic tuning scheme;
3. Results on multiple Bayesian inference benchmarks (e.g., Gaussian, Brownian, Item Response, Stochastic Volatility) and a statistical physics model with >1 million parameters.
weaknesses: -The main weakness of this paper is the lack of theoretical justification for non-Gaussian cases. The reviewer understand that this paper still provides enough novel contribution even without theoretical justification about non-Gaussian cases. However, it will be great if the authors can provides enough explaination/clarification regarding when this framework will work or not work for non-Gaussian cases.

---

confidence: 2
contribution: 2
initial_rating: 4
presentation: 3
questions: * Can you provide a concrete example to show how $\epsilon$ is adaptively changed when navigating through different parts of the landscape (for flatter regions, a larger $\epsilon$ is desired; for sharper regions, a smaller $\epsilon$ is desired)? Can the algorithm rediscover this trend?
* Why in Figure 1, the adjusted has a bigger variance than the unadjusted, despite the fact that they are both $1/t$?
soundness: 2
strengths: * The paper is nicely written and easy to follow
* The paper is supported by theoretical analysis and numerical experiments
* The proposed method is intuitive and seems effective
summary: Although adjusted sampling methods (which use Metropolis-Hastings algorithm for correction) do not have asymptotic bias due to finite step size, they are vulnerable to the curse of dimensionality. As such, the paper's goal is to find a good balance between efficiency and accuracy for unadjusted sampling methods, by adaptively changing step size based on the running estimate of EEVPD.
weaknesses: * The automatic tuning is too simple. This is not necessarily a weakness, but I want the authors to provide more context -- why haven't people come up with something similar before? If so, what is the conceptual challenge here? If not, how is the proposed method compared to previous methods that adaptively adjust step sizes?
* Similar to the first point, I'm not totally sure about the novelty and usefulness of the proposed method. My general impression is that HMC/LMC is a well-developed field. I'm willing to increase the score if the authors can justify more the novelty and contribution to the field.


Decision: Reject