1 Introduction↩︎

Underdamped Langevin dynamics on \(\mathbb{R}^{2d}\) with constant symmetric, positive-definite mass matrix \(M \in \mathbb{R}^{d\times d}\) and friction matrix \(\Gamma \in \mathbb{R}^{d\times d}\), is defined by the following SDE system: \[\begin{align} dx &= M^{-1}p\,dt, \tag{1}\\ dp &= -\nabla U(x)\,dt - \Gamma M^{-1} p\,dt + \sqrt{2k_{\mathrm{B}}T\Gamma}\,dW_t, \tag{2} \end{align}\] where \(U:\mathbb{R}^{d} \to \mathbb{R}\) is a potential, \((W_t)_{t\geq 0}\) denotes a vector of independent Wiener processes on \(\mathbb{R}^{d}\), \(T\) is the temperature and \(k_{B}\) is Boltzmann’s constant. Under mild assumptions, the invariant measure of 1 –2 is the canonical distribution with density defined pointwise up to a normalisation constant by \(\pi(x,p) \propto \exp(-H(x,p)/k_{\mathrm{B}}T)\), where the Hamiltonian is defined via \(H(x,p) = U(x)+\tfrac12 p^T M^{-1} p\), for any constant positive-definite \(\Gamma\). Additionally, under mild assumptions, the system is ergodic with respect to the invariant measure [1]. As a consequence, it is a popular method for sampling in molecular modelling and Bayesian computation.

The practical efficiency of 1 2 depends crucially on the choice of friction: if the friction is too weak, trajectories exhibit long ballistic excursions and may explore inefficiently, whereas if friction is too strong, momentum is rapidly overdamped and the dynamics approaches the slower overdamped Langevin regime [2], [3]. The familiar scalar choice \(\Gamma=\gamma I\) damps every direction equally. Introducing a friction matrix and allowing it to vary with position, \(\Gamma: \mathbb{R}^d \rightarrow \mathbb{R}^{d\times d}\), instead lets different directions be damped independently, and in a way that may be tuned to local geometric features of the landscape under study.

Many potentials of interest involve both very stiff directions (large curvature) and flat or slowly varying directions. In case a fixed stepsize is used, it must be small enough to resolve the stiffest mode [4], sacrificing efficiency in the flat regions where much larger steps would be safe. Adaptive (variable) timestepping integrators address this by contracting the step only where necessary, but doing so in a way that preserves the invariant measure is non-trivial.

In this work we propose an algorithm that combines variable timesteps with a specific position-dependent friction tensor as follows:

1. SamAdams variable timestep [5]. The integrator continuously senses how demanding the local dynamics are through the magnitude of the force, or in tighter geometries through the speed of the motion, and smoothly contracts the integration step wherever the demand is high, relaxing it again where possible. Because the effective step never exceeds a chosen nominal value, that nominal step can be set well beyond the stability threshold of a conventional fixed-stepsize integrator. The accompanying time-rescaling is corrected for in a simple reweighting procedure, so that Monte Carlo averages from the canonical distribution can be recovered. We use two variants of the local sensor, one driven by the force magnitude (the recommended default), and one that also responds to how fast the trajectory is changing, better suited to a narrow entropic bottleneck. Both options are defined in Section 2.

2. Position-adaptive Langevin (PAL). The friction tensor is taken to be the normalized rank-one-plus-scalar form \[\Gamma(x) = \alpha I + \beta\,\hat{f}(x)\hat{f}(x)^T, \qquad \hat{f}(x) = \frac{\nabla U(x)}{|\nabla U(x)|}, \label{eq:PAL-friction}\tag{3}\] with constants \(\alpha > 0\) and \(\beta > -\alpha\). The eigenvalues of \(\Gamma\) are the spatially constant values \(\alpha\) (transverse to \(\nabla U\)) and \(\alpha+\beta\) (parallel to \(\nabla U\)); only the eigendirection rotates with \(x\). This structure admits an \(O(d)\) exact solution of the Ornstein–Uhlenbeck (OU) step and guarantees that the canonical distribution is the exact invariant measure.

The idea of using a variable friction is not entirely new even if it is rarely used in practice. The paper of Lim and Tao  [6] considered this and the specific choice of the Hessian of the potential as the friction matrix and showed theoretically that the continuous dynamics could enhance performance (even if no numerical experiments were provided in that paper).

Due to the special form (scaled identity + rank-1 for friction, SamAdams variable stepsizes), the scheme we describe is also practically free of computational overhead compared to standard integrators. The two devices are combined in the AZBOBZA palindromic integrator described herein (Section 3.3), which requires a single force evaluation per step, the same cost as other typical 2nd order Langevin integrators [7], and passes all computed quantities through each sub-step without redundant work other than a few extra inner products.

In SA-PAL, the two devices address the temporal and directional components simultaneously, and this is shown in the numerical experiments of Section 4. It is observed that the PAL contribution is most visible in landscapes with clearly curved low-energy structure, such as the Rosenbrock potential. On the thin channel, much of the gain comes from the SamAdams timestep; here the fixed-step chains barely cross the corridor at all, so SA-PAL’s main effect is to pass through it far more often, although in experiments the position adaptive friction also provides a quantifiable benefit.

On the relatively smooth Mueller-Brown potential, the picture is more intricate and reveals a further benefit of the new scheme. Compared against the fastest-mixing isotropic friction the gain in position exploration looks modest, but at that friction, it turns out that the fixed stepsize solution is substantially biased: it mixes the positions quickly but crosses between the three wells too rarely to sample their relative populations correctly. The proper comparison is against the isotropic friction that matches SA-PAL’s accuracy, and in relation to this choice SA-PAL mixes the positions about \(2.8\times\) faster while reducing the potential-energy autocorrelation several times more quickly. The benefit PAL adds on top of SA is a directional decoupling of damping rates; it damps the force direction which drives the variation of energy without over-damping the slow inter-well motion. Thus, in this example, SA-PAL achieves a more favourable balance between exploration speed and equilibrium accuracy than any of the isotropic-friction baselines considered. This improvement in both mixing and accuracy is visible whenever the landscape has sufficient complexity.

Beyond these low-dimensional examples, we also apply SA-PAL to a realistic, higher-dimensional target, namely, Bayesian linear regression under a horseshoe prior, using an example borrowed from [8], whose posterior has the sharp “funnel” geometry that is notoriously difficult for fixed-step samplers. Here the combined scheme delivers order-of-magnitude reductions in integrated autocorrelation time relative to a well-tuned fixed-step baseline while keeping the estimated posterior moments accurate.

1.0.0.1 Related work.

Adaptive-friction Langevin samplers in the literature can be grouped into several distinct families, and the rank-one-plus-scalar friction 3 adopted here does not coincide with any of them.

Chak et al [9] analyze the optimal choice of a constant (state-independent) friction matrix that maximizes the spectral gap of underdamped Langevin dynamics; the friction in that article does not vary with position. By contrast, the work of Lim and Tao [6] constructs a position-dependent friction tensor built from the Hessian \(\nabla^2 U\), motivated by exact preconditioning of harmonic modes.

A further line of research, distinct from those above, treats the friction as a dynamical variable governed by an auxiliary equation rather than as a chosen function of \(x\). The projective thermostatting technique of [10] applies a Nosé–Hoover-type thermostat through a projection onto a selected set of modes, this is particularly useful in multiscale settings where only a subset of modes require explicit temperature control. While related to the current work, the aim there was not to improve sampling, but to treat out-of-equilibrium physical models, and the method was in fact purely deterministic. Also close in spirit are the recent projection-based optimisation schemes discussed in [11], [12].

There has also been a substantial body of work that replaces the underdamped SDE by a preconditioned variant in which a position-dependent matrix \(B(x)\) multiplies the drift and diffusion, in close analogy with Riemannian-manifold Hamiltonian Monte Carlo [13]. A common structural difficulty of all such schemes is that preservation of the canonical measure requires the SDE drift to incorporate a divergence correction. An analogous correction is needed for the kinetic-energy preconditioner in the underdamped formulation and computing the matrix divergence is a real impediment in practice. Analytic preconditioners cost an additional \(d\) matrix-derivative entries per step.

The ensemble preconditioner of [14] attacks the same problem from a different angle: instead of computing \(\nabla\!\cdot\!B\), the authors build \(B\) from a localized empirical covariance estimated over those walkers within a distance-based cutoff of the current configuration. The effective \(B\) is then not a smooth matrix field but a discontinuous function of the joint ensemble state, depending on which walkers happen to fall inside the cutoff and on the cutoff radius itself. Divergence identities therefore do not directly apply, and the question of measure preservation is replaced by a more delicate one about the multi-walker joint dynamics, with the cutoff radius as an additional and somewhat ad hoc hyperparameter.

The friction tensor 3 avoids the divergence problem altogether, by virtue of its incorporation into an Ornstein-Uhlenbeck part of a kinetic Langevin dynamics formulation rather than via a coordinate transformation of the physical variables. The fluctuation–dissipation operator \[\mathcal{L}_{\mathrm{FD}} = -(M^{-1}p)^T \Gamma(x)\nabla_p + k_{\mathrm{B}}T\,\Gamma(x):\nabla_p^2\] preserves the canonical density \(\pi\propto e^{-H/k_{\mathrm{B}}T}\) pointwise for any smooth positive-definite \(\Gamma(x)\), with no drift correction involving \(\nabla\!\cdot\!\Gamma\) required (Section 3).

The lack of any need for a divergence correction is of course understood within the Hessian-based friction of [6]. What is specific to 3 is the specific form of \(\Gamma\) itself, involving a projection matrix based on the gradient rather than the Hessian which concentrates all position dependence in the rotation of a single eigendirection \(\hat{f}(x)\). The two eigenvalues \(\alpha\) and \(\alpha+\beta\) are spatially invariant, no derivative of \(\Gamma\) appears anywhere in the integrator, and the canonical distribution is the correct invariant measure. For rigorous treatment of the subsequent dynamics, the hypoelliptic property and the consequent ergodicity of Langevin system with general position-dependent friction is treated in [15]; that framework covers the PAL form 3 as a special case, providing existence of an invariant measure and convergence guarantees.

1.0.0.2 Organisation.

Section 2 summarizes the SamAdams adaptive-timestep device, including the relaxation ODE governing the running sensor, options for monitor functions, choice of the Sundman transformation and the time-reweighting identity by which canonical expectations are recovered from chains generated in the rescaled time. Section 3 develops the position-adaptive Langevin (PAL) dynamics, including the specific rank-one-plus-scalar form \(\Gamma(x)=\alpha I + \beta\,\hat{f}(x)\hat{f}(x)^T\), discuss its physical interpretation, and derives an exact \(O(d)\) closed-form Ornstein–Uhlenbeck update based on it that avoids matrix decomposition. Section 3.3 assembles the SamAdams scheme with the PAL friction into the AZBOBZA palindromic integrator: it defines the four elementary flows (A, Z, B, O) and explains why a particular ordering is needed. Section 4 reports the numerical experiments on three benchmark toy model potentials (each in two spatial coordinates) and discusses parameter choices. Numerical experiments are used to compute autocorrelation functions and integrated autocorrelation times at matched gradient budget. We discuss the stepsize distribution, equilibrium accuracy of certain moments and provide an ablation study that clarifies the contribution of the PAL friction on top of SamAdams. Finally in Section 5 we summarize the findings and outline future work, including bias analysis of the SamAdams device and extensions to higher-dimensional Bayesian and machine-learning targets.

2 SamAdams variable timestep↩︎

In this section, we describe the variable timestepping approach of [5] on which the SA-PAL scheme is based. Consider a system with potential \(U\) that has both stiff and soft directions. If the curvature in the stiff direction is \(\lambda_{\max}\), then stability of any Hamiltonian-based splitting integrator requires \(h < C/\sqrt{\lambda_{\max}}\) for some method-dependent constant \(C\), see for example [16]. In a region where the curvature is mild, however, the dynamics would be stable and accurate with significantly larger steps. A variable stepsize scheme can exploit this by automatically detecting local stiffness.

A local stiffness indicator, the monitor function, is a positive scalar function \(g(x,p)\) that grows in regions where the dynamics needs to be slowed down. Two specific choices have proved useful in the experiments of this paper: \[\tag{4} \begin{align} g(x) \;&=\; |\nabla U(x)|/\Omega, \qquad &&\text{(force monitor)} \tag{5}\\ g(x,p) \;&=\; \bigl(|\nabla U(x)|^2 + |M^{-1}p|^2\bigr)/\Omega^2, \qquad &&\text{(arc-length monitor)} \tag{6} \end{align}\] The force monitor 5 grows when the force is large; the scale \(\Omega>0\) is chosen near the typical magnitude of \(|\nabla U|\) over the target distribution. This is the monitor used in the original SamAdams algorithm [5], where the Langevin dynamics is coupled to an auxiliary scalar relaxation that senses stiffness through \(g\) and rescales the step accordingly.

The arc-length monitor 6 is the squared norm of the Hamiltonian part of the vector field, \((\dot{x},\dot{p})=(M^{-1}p,-\nabla U)\). To ensure the stepsize does not get too large or small, we use a function \(\psi\) and stepsizes \(\eta=\psi(g)\), where \(\psi\) is decreasing in \(g\). The effective step then contracts wherever this velocity is large, so the integrator advances in near-equal increments of arc length along the trajectory rather than in equal time. This results in a stochastic counterpart of the arc-length time-reparameterisation underlying adaptive integrators for deterministic Hamiltonian systems [17], [18]. Retaining the force term lets the step also contract in stiff regions even where the velocity is momentarily small. This control is particularly useful near entropic bottlenecks, where the force can be small yet the geometry is tight and the velocity large. In both cases \(\Omega\) provides the global sensitivity scale and a larger \(g\) signals the need for a smaller effective step. Only the force monitor was studied in [5].

The SamAdams algorithm uses a relaxation approach, in particular, it introduces a sensor variable \(Z\) obeying \[\frac{dZ}{dt} \;=\; -\theta\bigl(Z - g(x,p)\bigr), \label{eq:sensor-ode}\tag{7}\] so that \(Z\) relaxes towards the current value of the monitor at rate \(\theta > 0\). A large \(\theta\) makes \(Z\) track \(g\) closely; a small \(\theta\) gives a time-averaged memory and smoother adaptation. The averaging also helps to improve the robustness of the method in the presence of gradient noise.

For fixed \(x,p\), 7 can be solved over a half-step of length \(h/2\) resulting in the SamAdams update: \[Z \;\leftarrow\; e^{-\theta h/2}\,Z \;+\; (1-e^{-\theta h/2})\,g(x,p), \label{eq:sensor-update}\tag{8}\] which is a moving average. The effective scale factor, \(\eta\), is then defined through a Sundman kernel \(\psi\) (see [19]) such that \(\eta \;=\; \psi(Z)\), and the nominal stepsize \(h\) is replaced by the effective stepsize \(h_\mathrm{eff} = h\,\eta\). We use two kernels, the first of which is defined with the force monitor; the bounded SamAdams kernel [5], \[\psi^{(2)}(Z) \;=\; m_{\mathrm{lo}}\,\frac{Z^{r} + m_{\mathrm{hi}}/m_{\mathrm{lo}}}{Z^{r} + 1} \;\in\; [m_{\mathrm{lo}},\, m_{\mathrm{hi}}], \qquad 0 < m_{\mathrm{lo}} \le m_{\mathrm{hi}} \le 1,\;\; r > 0, \label{eq:psi-bounded}\tag{9}\] with bounds \(\psi^{(2)}(0)=m_{\mathrm{hi}}\) and \(\psi^{(2)}(Z{\to}\infty)=m_{\mathrm{lo}}\), and a transition centered on \(Z=1\) (i.e.on \(|\nabla U|=\Omega\) at equilibrium) of sharpness controlled by \(r\). The parameter \(m_{\mathrm{hi}}\) caps the maximum step in quiescent regions (\(m_{\mathrm{hi}} h\) should respect the stability threshold), while \(m_{\mathrm{lo}}\) guarantees the chain never freezes, i.e. \(h\eta \ge m_{\mathrm{lo}} h>0\) always.

For the “arc-length” monitor we use instead the simple unbounded form \[\psi(Z) \;=\; (1+Z)^{-1/2} \;\in\;(0,1], \label{eq:psi-classic}\tag{10}\] which decays smoothly as \(Z\) grows and imposes no upper bound, so here the nominal \(h\) itself plays the role of the maximum step. In either case \(\eta\le 1\), so \(h\) need only satisfy a stability condition “on average” and may be chosen far larger than the fixed-step threshold.

3 Position-adaptive Langevin dynamics↩︎

The underdamped Langevin dynamics with position-dependent friction and constant symmetric positive-definite mass matrix \(M\in\mathrm{Sym}^{d\times d}_{+}\) are given by \[\begin{align} dx &= M^{-1}p\,dt, \\ dp &= -\nabla U(x)\,dt-\Gamma(x)M^{-1}p\,dt +\sqrt{2k_{\mathrm{B}}T\,\Gamma(x)}\,dW_t, \label{eq:pal-sde} \end{align}\tag{11}\] a minor modification of 1 2 for any smooth, uniformly positive-definite friction tensor \(\Gamma:\mathbb{R}^d\to\mathrm{Sym}^{d\times d}_{+}\). Under the usual regularity and confinement assumptions, these dynamics preserve the canonical distribution with density \[\pi(x,p) = Z^{-1}\exp\left(-\frac{H(x,p)}{k_{\mathrm{B}}T}\right), \qquad H(x,p)=U(x)+\frac{1}{2} p^T M^{-1}p,\] for a normalisation constant \(Z > 0\).

To verify invariance, consider the fluctuation–dissipation part of the forward Kolmogorov operator. Acting on a density \(\rho\), it is \[\mathcal{L}_{\mathrm{FD}}^*\rho = \nabla_p\cdot \left[ \Gamma(x) \left( M^{-1}p\,\rho+k_{\mathrm{B}}T\nabla_p\rho \right) \right].\] Since \(\nabla_p\pi =-(k_{\mathrm{B}}TM)^{-1}p\,\pi,\) it follows immediately that \(\mathcal{L}_{\mathrm{FD}}^*\pi=0\) pointwise. In particular, no derivative or divergence correction involving \(\Gamma(x)\) is required and the diffusion acts only in the momentum variables, while \(x\) is fixed under the fluctuation–dissipation dynamics. The Hamiltonian part also preserves \(\pi\), since its flow is divergence-free and conserves \(H\). Hence the full forward operator satisfies \(\mathcal{L}^*\pi=0\).

Under suitable growth conditions on \(U\), together with uniform non-degeneracy and sufficient regularity of \(\Gamma\), the resulting process is hypoelliptic and ergodic with respect to \(\pi\); see [15]. Although the choice \(M=I\) is often made for simplicity, a non-trivial constant mass matrix may be used to precondition anisotropic target distributions.

3.1 The normalized rank-one friction↩︎

We adopt and make use of the following choice of position-dependent friction \[\Gamma(x) = \alpha I+\beta\,\hat{f}(x)\hat{f}(x)^T, \qquad \hat{f}(x)=\frac{\nabla U(x)}{|\nabla U(x)|}, \label{eq:Gamma-norm}\tag{12}\] with \(\alpha>0\) and \(\alpha+\beta>0\)³. Since \(\hat{f}(x)\hat{f}(x)^T\) is the orthogonal projector onto the force direction, \(\Gamma(x)\) has eigenvalue \(\alpha+\beta\) along \(\hat{f}(x)\) and eigenvalue \(\alpha\), with multiplicity \(d-1\), in the orthogonal directions. Thus, the anisotropy ratio \((\alpha+\beta)/\alpha\) is spatially constant, while the principal dissipation direction rotates with the force.

For \(\beta>0\), momentum parallel to the force is damped more strongly than momentum tangent to the local level sets of \(U\). This promotes rapid thermalisation in the energy-changing direction while retaining weaker damping along equipotential contours.

3.2 Exact \(O(d)\) Ornstein–Uhlenbeck update↩︎

For \(\Gamma\) of the form 12 with \(\hat{u}\) a fixed unit vector (to denote \(\hat{f}(x)\) at position \(x \in \mathbb{R}^{d}\)) and isotropic mass \(M = m I\), the OU process given by the solution to \[dp = -\Gamma (p/m)\,dt + \sqrt{2k_{\mathrm{B}}T\Gamma}dW_t,\] can be solved exactly (in the weak sense) over time \(\Delta t > 0\). We then use the following convenient decomposition of the momentum \(p = p_\parallel \hat{u} + p_\perp\), with \(p_\parallel = \hat{u}\cdot p\), \(p_\perp = p - p_\parallel\hat{u}\) whose solutions satisfy \[p_\parallel(\Delta t) = e^{-(\alpha+\beta)\Delta t/m}\,p_\parallel(0) + \sqrt{mk_{\mathrm{B}}T\bigl(1 - e^{-2(\alpha+\beta)\Delta t/m}\bigr)}\,\xi_\parallel, \label{eq:OU-par}\tag{13}\] \[p_\perp(\Delta t) = e^{-\alpha\Delta t/m}\,p_\perp(0) + \sqrt{mk_{\mathrm{B}}T\bigl(1 - e^{-2\alpha\Delta t/m}\bigr)}\,\xi_\perp, \label{eq:OU-perp}\tag{14}\] where \(\xi_\parallel \sim \mathcal{N}(0,1)\) and \(\xi_\perp \sim \mathcal{N}(0,I_{d-1})\) are independent. Due to the fact that \(\alpha/m\) and \((\alpha+\beta)/m\) are global constants, the four scalar factors in 13 –14 can be precomputed once for a given parameter selection \((\alpha,\beta,m,\Delta t)\); only the direction \(\hat{u}\) changes at each step.

For a general mass matrix \(M\), the OU update can still be performed exactly in the basis of \(M^{-1}\Gamma\); see Appendix 6 for a general construction.

3.3 The AZBOBZA integrator↩︎

We combine the Langevin system 1 –2 with the PAL friction 12 and the SamAdams sensor 8 into four elementary flows, each of which can be executed exactly. The A step which advances the position using the momentum, the B step which advances the momentum using the force, the O-step as an exact Ornstein-Uhlenbeck solve as described in Section 3.2, all of which are present in classical ABO-based splitting integrators [7]. Then additionally, we denote the Z-step which is used for adaptive stepsizing, which given \(\nu = e^{-\theta h/2}\) and the current monitor value \(g(x,p)\) (either 5 or 6 ) is defined via 8 , where \(\psi\) is either the unbounded kernel 10 or the bounded kernel 9 as appropriate. The Z step uses the force \(F=-\nabla U(x)\) already computed and generates no additional gradient evaluation.

We then use the splitting method ordering AZBOBZA, where the position \(x\) is advanced only by the two outer \(A\) half-steps; the intervening \(Z\), \(B\), \(O\), \(B\), \(Z\) operators all leave \(x\) fixed. A single force \(F=-\nabla U(x)\), evaluated once after the first \(A\) half-step, therefore serves both \(B\) momentum kicks, the \(O\) update (through \(\hat{F}=F/|F|\)), and both \(Z\) sensor updates. A comparison of integrated autocorrelation times per step is thus a comparison per force evaluation. As in [5] one can show that under appropriate assumptions the AZBOBZA integrator has weak error of order two by using [20].

The complete single-step update for AZBOBZA, starting from state \((x,p,Z,\eta)\) with the pre-computed constant \(\nu=e^{-\theta h/2}\), is given as pseudocode as follows.

4 Numerical experiments↩︎

Before turning to a problem in higher dimensions, we test the method on three two-dimensional model potentials, chosen to exercise distinct pathologies: curved metastable basins (Rosenbrock), strongly anisotropic geometry with an entropic bottleneck (thin channel), and multiple metastable states with rare inter-well transitions at elevated temperature (Mueller–Brown).

For each model we obtain reference marginal densities by direct numerical quadrature of the Gibbs distribution with density \(\pi(x_1,x_2)\propto\exp(-U(x_1,x_2)/k_{\mathrm{B}}T)\) rather than by simulation.

4.1 Model Problems↩︎

4.1.1 Rosenbrock (\(k_{\mathrm{B}}T= 1\)).↩︎

For \(x_{1}\), \(x_{2} \in \mathbb{R}\), we define our Rosenbrock potential function pointwise in \(\mathbb{R}^{2}\) via \[U(x_1, x_2) = (1-x_1)^2 + 100(x_2 - x_1^2)^2.\] The distribution is a narrow ridge following the Rosenbrock valley. Fixed friction and stepsize integrators such as BAOAB and OBABO typically require a moderately large friction to suppress oscillations along the ridge.

4.1.2 Thin entropic channel (\(k_{\mathrm{B}}T= 1\)).↩︎

For \(x_{1}\), \(x_{2} \in \mathbb{R}\), we define our thin entropic channel potential function pointwise in \(\mathbb{R}^{2}\) via \[U(x_1, x_2) = \frac{100 x_2^2}{1 + 10 x_1^4} + 0.001(x_1^2 - 9)^2.\] This is constructed such that the channel narrows dramatically near \(x_1 = 0\), where the transverse curvature is \(200\). Fixed stepsize Langevin integrators, like BAOAB or OBABO are slow to traverse the narrow neck, often failing to pass through at all when the stepsize is large.

4.1.3 Mueller–Brown (\(k_{\mathrm{B}}T= 15\)).↩︎

The Mueller–Brown potential function from [21] is defined for \(x_{1}\), \(x_{2} \in \mathbb{R}\) by the sum of four shifted anisotropic Gaussians \[U(x_1,x_2) \;=\; \sum_{k=1}^{4} A_k \exp\!\Bigl(a_k(x_1-x_{1,k})^2 + b_k(x_1-x_{1,k})(x_2-x_{2,k}) + c_k(x_2-x_{2,k})^2\Bigr),\] with the commonly used parameters \[\begin{array}{r|rrrr} k & 1 & 2 & 3 & 4 \\ \hline A_k & -200 & -100 & -170 & 15 \\ a_k & -1 & -1 & -6.5 & 0.7 \\ b_k & 0 & 0 & 11 & 0.6 \\ c_k & -10 & -10 & -6.5 & 0.7 \\ x_{1,k} & 1 & 0 & -0.5 & -1.0 \\ x_{2,k} & 0 & 0.5 & 1.5 & 1.0 \end{array}\] This surface has three local minima connected by two saddle points; at \(k_{\mathrm{B}}T= 15\) inter-well transitions occur at a physically meaningful rate.

4.2 Parameter choices↩︎

There are a collection of parameters both in the fixed stepsize, fixed friction integrators and the adaptive methods. We detail our choices of parameters below.

The AZBOBZA nominal stepsize \(h_{\mathrm{sa}}\) is chosen significantly larger than the stability-limited BAOAB/OBABO stepsizes; the SamAdams sensor automatically reduces the effective step in stiff regions. PAL parameters \((\alpha,\beta)\) are selected via a short two-stage grid sweep (described in the supplementary material); the Mueller–Brown values reflect the choice \(\gamma_{\mathrm{opt}} \approx 10\) for Mueller–Brown with \(k_{\mathrm{B}}T= 15\). The Rosenbrock force-monitor row is retained for the monitor comparison in Table 1, whereas the arc-length row is the Rosenbrock configuration used in the main Rosenbrock figures and in the SA-only versus SA-PAL ablation. When a figure uses a fixed-step friction or stepsize sweep, the exact fixed-step values are given in the corresponding caption or table; the entries above are the default model settings.

The SamAdams parameters (sensor relaxation rate \(\theta\), kernel form, bounds \(m_{\mathrm{lo}},m_{\mathrm{hi}}\), sharpness \(r\), scale \(\Omega\)) are summarized below. The force monitor 5 with the bounded kernel 9 is used by default; the channel and the Rosenbrock arc-length configuration use the monitor 6 with the unbounded kernel 10 . On the channel this is needed because the rate-limiting feature is the entropic narrowing at \(x_1=0\) (where \(|\nabla U|\approx0\)), which the force-only monitor cannot detect. The arc-length monitor is tested in the Mueller-Brown example see its effect on accuracy as well as convergence rate.

4.3 Rosenbrock function↩︎

The Rosenbrock potential is a narrow, curved valley flanked by stiff walls. It exposes a dilemma intrinsic to fixed-step Langevin integrators: the largest stable stepsize depends on the friction as well as the conservative force. At low friction the dynamics are near-ballistic, so the chain gains energy and probes the stiff walls, forcing a small stable \(h\); heavier friction damps these excursions and permits a larger stepsize, but reduces the diffusion along the valley. Figure 1 (a) shows this stability boundary \(h_{\mathrm{crit}}(\gamma)\) rising with \(\gamma\): \(\gamma=0.5\) is unstable already at \(h=0.02\) (it requires \(h=0.01\)), whereas \(\gamma=5\) allows \(h=0.02\). BAOAB is therefore squeezed, meaning that low friction values improve mixing, but limit the stepsize stability threshold. Figure 1 (b,c) and Table 1 demonstrate the consequence: BAOAB’s best achievable per-gradient IACT for the valley coordinate is about \(480\), whereas the retuned arc-length SA-PAL run gives \(\tau_{x_1}\approx291\). The moderately damped stable choice \(\gamma=5,h=0.02\) is nearly six times slower on \(x_1\) despite taking larger fixed steps. For the stiff transverse coordinate the arc-length SA-PAL run improves on the best low-friction fixed-step baseline and is much faster than the stable \(\gamma=5,h=0.02\) row.

Rosenbrock also illustrates the difference between the two Sundman monitors of Section 2. With the force monitor \(\psi^{(2)}\), which contracts the step only where \(|\nabla U|\) is large, SA-PAL accelerates the valley coordinate \(x_1\) and the energy but leaves the stiff transverse mode \(x_2\) essentially at the BAOAB rate: the force is modest during the fast transverse oscillations, so a force-only monitor does not register them. The arc-length monitor 6 also responds to large velocities, and it is the configuration used in the Rosenbrock figures and ablation study. The monitor comparison in Table 1 should therefore be read as a sensitivity check: both choices improve the slow valley coordinate substantially, while the transverse mode remains close to the best low-friction fixed-step baseline.

In light of Figure 1 (d,e), the \(\gamma=5,h=0.02\) BAOAB run is stable but visibly more damped than the low-friction BAOAB rows, and its valley-coordinate IACT remains high (\(\tau_{x_1}\approx 1673\)). Running BAOAB at \(\gamma=5\) with \(h=0.01\) makes the physical step comparable to the adaptive run but slows the valley coordinate further.

SA-PAL keeps the isotropic friction low and uses the adaptive step to handle stiff excursions, so the velocities retain useful memory without destabilizing the trajectory. The manuscript uses \(h_{\mathrm{sa}}=0.040\), for which the ordinary mean effective step \(N^{-1}\sum_n h_{\mathrm{sa}}\eta_n\) is about \(0.0106\), only modestly larger than BAOAB’s \(h=0.01\) low-friction runs and still below the stable \(\gamma=5\) fixed step. Thus the improvement is not merely a hidden larger-step comparison: it comes from combining underdamped motion with local step contraction in stiff regions.

4.4 The thin channel↩︎

The thin channel isolates the regime where the SamAdams adaptive timestep is essential. The two basins at \(x_1=\pm3\) communicate only through a narrow entropic neck at \(x_1=0\), where the transverse confinement stiffens sharply. A fixed-step integrator must pick \(h\) small enough to stay stable in the neck, yet that same \(h\) is wastefully small throughout the wide basins, so each gradient evaluation advances the dynamics very little. Figure 2 shows that, at an equal number of gradient evaluations, BAOAB explores poorly at every friction. In the overdamped regime (\(\gamma=10\)), BAOAB is trapped in one basin; in the underdamped regime, it diffuses slowly, whereas SA-PAL takes large steps in the basins and automatically shrinks \(\eta\) in the neck, allowing much more frequent passes through the channel. Table 2 quantifies this: at equal cost SA-PAL makes about \(16\times\) as many neck crossings as the best-tuned BAOAB and is the only scheme whose \(\tau_{x_1}\) actually converges. Fixed-step BAOAB crosses the neck so rarely (a few hundred times in \(2\times10^7\) steps, and even less often at the smaller \(h=0.001\)) that its \(\tau_{x_1}\) is only a loose lower bound.

4.5 Mueller–Brown↩︎

The Mueller–Brown potential (with \(kT=15\)) is the classic three-well benchmark, with curved saddle paths connecting the basins. Here the fixed-step baseline has a genuine optimal friction: BAOAB mixes best at intermediate \(\gamma\approx10\) (Table 3), since low friction is slow to cross the saddles (under-damped, ballistic recrossing) and high friction over-damps the diffusion. Figure 3 (top) shows all methods explore the three wells at this temperature, while the per-variable autocorrelations (bottom) quantify the mixing: force-monitor SA-PAL (\(\alpha=0.5,\beta=10\)) decorrelates \(x_1\), \(x_2\) and \(U\) faster than BAOAB at any friction. Relative to the best-tuned baseline (\(\gamma=10\)), the per-gradient IACT improves by about \(2.1\times\) for \(x_1\), \(2.0\times\) for \(x_2\), and \(1.7\times\) for the potential energy \(U\). The gain for \(U\) is consistent with the structure of the dynamics: since \(\dot{U} = \nabla U\cdot\dot{x}\), fluctuations of the potential energy are carried by the velocity component along the force direction \(\hat{f}=\nabla U/\lVert\nabla U\rVert\), and SA-PAL’s rank-one friction \(\beta\,\hat{f}\hat{f}^{\top}\) damps precisely that component (here with \(\beta=10\)), driving the energy-carrying modes toward critical damping. The configurational coordinates \(x_1,x_2\) instead depend on the slow inter-well transitions. The comparison is at matched accuracy: both methods exhibit a comparable discretisation bias at their operating stepsizes (about \(2\%\) on the second moments, verified to be independent of chain length rather than a sampling artefact), so SA-PAL is not trading accuracy for speed. We also find the anisotropic friction is important for accuracy as well as mixing: the SamAdams-only limit \(\beta=0\) samples the well populations poorly (bias \(>10\%\)), which the rank-one friction \(\beta>0\) removes. For comparison we also include the arc-length SamAdams monitor used in the Rosenbrock and channel examples. On Mueller–Brown it is accurate but takes smaller effective steps, and consequently mixes more slowly than the force monitor. In some sense it represents a useful compromise, offering better accuracy than force-based SA-PAL while still giving better mixing than the accurate fixed-stepsize methods with low friction.

4.6 Effective-stepsize distributions↩︎

A key diagnostic for the SamAdams device is the distribution of the effective stepsize \(h_{\mathrm{sa}}\eta\) taken along converged SA-PAL chains, where \(\eta = \psi(Z)\) is the adaptive scale factor. Figure 4 shows these distributions for the three tuned models, with both Mueller–Brown monitor choices overlaid in the right panel. The arc-length kernel 10 (Rosenbrock, channel, and the Mueller–Brown comparison row) produces a continuous spread, while the bounded force-monitor kernel 9 (the main Mueller–Brown row) confines \(\eta\) to \([m_{\mathrm{lo}},m_{\mathrm{hi}}]\). The contrast with the fixed BAOAB step (dashed line) is informative: on the channel the mean effective step is roughly \(24\times h_{\mathrm{BAOAB}}\), exploiting the wide basins, whereas on Rosenbrock it is actually smaller than \(h_{\mathrm{BAOAB}}\). SA-PAL’s advantage there appears to be due to improved in-basin mixing, not a larger step.

4.7 The contribution of PAL on top of SamAdams↩︎

The results so far combine the two devices. To isolate the contribution of PAL, we repeat the runs with \(\beta=0\), which reduces the friction to the isotropic form \(\gamma I\) and leaves only the SamAdams timestep active. All other settings including the nominal step, the sensor parameters, and the reweighting are unchanged. We call this variant SA-only and run it at three friction values per model. Table 4 reports the results.

There is an important distinction between the fixed-step and adaptive-step comparisons. The fixed-step BAOAB baselines are primarily stability-limited in these examples: once the step is below the stability threshold, their equilibrium bias remains small until the step is pushed close to instability. By contrast, SamAdams and SA-PAL are often accuracy-limited before they are stability-limited. If the nominal adaptive step is made sufficiently small, both SA-only and SA-PAL produce low bias; near the largest useful adaptive steps, however, SA-PAL is noticeably more accurate than SA-only. Thus the fair comparison is not at the same aggressive nominal step when one row is biased, but at comparable sampling accuracy. Under this matched-accuracy convention, SA-PAL converges much more rapidly than the corresponding SA-only operating point.

Most of the speedup over fixed-step methods is already present in SA-only: the adaptive timestep alone accounts for the bulk of the gain. On Rosenbrock the best SA-only run is slightly faster than the best BAOAB, and on the channel it reduces the within-well autocorrelation times by an order of magnitude. PAL adds a further improvement on top when the comparison is restricted to rows with acceptable bias. On Rosenbrock this is clear-cut in mixing: SA-PAL has the lowest IACTs in Table 4, while the very high friction SA-only row has the smallest bias but mixes slowly. The curved valley gives the force-aligned friction a well-defined direction to act along, which is the case PAL is designed for.

Mueller–Brown is more instructive, because here an isotropic friction must cope with a trade-off: small constant friction mixes the positions quickly but samples the basins inaccurately, because the chain crosses between the wells too rarely. A large friction corrects the accuracy but slows the position mixing. Table 4 shows that at low and intermediate friction SA-only has visibly high bias (these rows are colored red and are not used as baselines), and bringing the bias down requires a friction at which the position autocorrelation times are several times larger.

PAL with either monitor addresses this trade-off. It damps strongly along the force direction, which carries the potential energy, while leaving the slow inter-well motion lightly damped. SA-PAL therefore combines the position mixing of the low-friction case with a bias as low as the high-friction case. Against the one isotropic friction that matches its accuracy, SA-PAL mixes the slow coordinates \(2.7\times\) faster. The arc-length monitor has better bias than the force-based SA-PAL on Mueller–Brown (as good as fixed stepsize), but its effective steps are smaller and its IACTs lie between the force-monitor SA-PAL row and the high-friction SA-only row.

An intermediate friction does not help: it is faster but still too biased to use. For this reason the low-friction SA-only curves in Figure 5 are drawn dashed in pink/red, since their faster decay does not correspond to correct sampling; the only accurate (solid) SA-only option has much slower mixing than SA-PAL. We set \(\beta\) to the largest value that keeps the bias below the fixed-step reference; a smaller \(\beta\) gives slightly faster position mixing at the cost of more bias.

4.8 Equilibrium accuracy at fixed budget↩︎

Having established the mixing efficiency on each model and the separate contribution of each device, we finally report the accuracy of ensemble-averaged observables. All methods are run for the same number of gradient evaluations, so the comparison is at matched computational budget.

Figure 6 summarizes these estimates. On Rosenbrock SA-PAL is the most accurate scheme and also the least variable: its boxes lie closest to the target lines and are the narrowest of any method. In the thin channel there is little gain in the invariant measure accuracy for SA-PAL over SA-only. In the Mueller-Brown model the accuracy of SA-PAL roughly matches that of BAOAB at sufficiently high friction (of course achieving this with much less computation due to the lower IACT). The fixed-step integrators behave as one expects when sampling the invariant measure. BAOAB is generally accurate given enough computation, and is typically most accurate at the higher frictions, that is, where it mixes slowest. Therefore its accuracy and its mixing speed pull in opposite directions. OBABO is the least accurate overall: its bias is largest at low friction, and conspicuously so on Mueller–Brown, where the low-friction chains miss the target by a wide margin. SA-only generally improves on both fixed-step methods, most clearly on the thin channel, where the adaptive step samples the slow coordinate much better at the same budget, but does not help much in the Mueller–Brown. The added Mueller–Brown arc-length box is close to the target lines, confirming that it has improved invariant-measure bias compared to SA-PAL but at the cost of some mixing efficiency.

4.9 Bayesian linear regression with horseshoe prior↩︎

Our last example is a Bayesian linear regression with a horseshoe prior [8]. This is a higher-dimensional target, and a harder one. The horseshoe is a shrinkage prior, meaning it pulls small coefficients towards zero while leaving the large ones essentially alone, and the price is a posterior that is sharply curved near the origin. We use \(N=100\) observations and \(K=25\) predictors, with true coefficients \(\beta_{1:4} = [3, -2, 1.5, 1]\) and \(\beta_{5:25} = 0\), noise scale \(\sigma=0.5\), and a design matrix \(X\) with i.i.d.Gaussian entries \(X_{ij} \sim \mathcal{N}(0, 1/N)\). Writing the state as \(w = [\beta, \eta, \xi]^T \in \mathbb{R}^{2K+1}\), the potential is \[U(\beta, \eta, \xi) = \frac{\|y - X\beta\|^2}{2\sigma^2} + \sum_{j=1}^K \frac{\beta_j^2}{2 e^{2(\xi + \eta_j)}} + (K-1)\xi + \sum_{j=1}^K \log(1 + e^{2\eta_j}) + \log(1 + e^{2\xi}),\] where \(X \in \mathbb{R}^{N \times K}\) is the design matrix, \(y \in \mathbb{R}^N\) is the response, and \(\sigma\) is the noise scale. The difficulty lives in the \((\beta_j, \xi+\eta_j)\) plane, where the posterior has the familiar “funnel” shape. When a coefficient sits near zero its scale \(\xi+\eta_j\) can shrink as well, leaving a narrow neck of high curvature; further out the funnel opens up and large steps are safe. A single fixed step cannot serve both ends: a step small enough for the neck is wasteful in the mouth, and a step sized for the mouth is unstable in the neck. We give SamAdams a step-size floor of \(m=10^{-4}\) so the step can drop far enough to stay stable in the neck.

We compared BAOAB, SA-only (\(\beta=0\)), and SA-PAL, using 20 independent seeds for each configuration, looking at three distinct fidelity regimes. For SA-PAL we used \(h_{\mathrm{sa}}=0.1\) in the high-fidelity regime, \(h_{\mathrm{sa}}=0.2\) in the moderate regime, and \(h_{\mathrm{sa}}=0.28\) in the high-performance regime, reduced slightly from \(0.3\) to keep the bias below \(2.5\%\) tolerance. For SA-only we chose \(h_{\mathrm{sa}}\) by the same bias rule: within each regime we swept the nominal SA-only step downward and kept the largest value whose estimated bias in \(\langle|\beta|^2\rangle\) was below \(2.5\%\). This gives \(h_{\mathrm{sa}}=0.1\), \(0.16\), and \(0.20\) for the high-fidelity, moderate, and high-performance SA-only rows, respectively. BAOAB was run at \(h \in \{10^{-3}, 5\cdot 10^{-3}, 10^{-2}\}\); on this funnel its largest reliably stable step is \(5\cdot 10^{-3}\), while \(h=10^{-2}\) already sits at the edge of stability and is seed-sensitive (some runs diverge). That fixed-step BAOAB cannot be pushed to such a step is precisely the regime the adaptive scheme is built for.

Figure 7 collects the results. The top panel reports the effective IACT of the aggregate coefficient observable \(|\beta|^2\), the same quantity used for the accuracy check in the bottom panel. For reweighted SamAdams chains this effective IACT is computed from the long-run variance of the delta-method ratio-estimator influence sequence and normalized by the weighted target variance of \(|\beta|^2\), as described in Appendix 7. With this convention, SA-PAL is the fastest aggregate mixer in all three regimes among rows satisfying the \(2.5\%\) bias tolerance. The gains over BAOAB are large because the fixed-step chain must remain stable in the narrow funnel neck. After bias matching the SA-only rows, SA-PAL reduces \(\tau(|\beta|^2)\) by factors of about \(1.88\), \(2.70\), and \(2.04\) in the high-fidelity, moderate, and high-performance regimes, respectively. The bottom panel shows why the bias convention matters: increasing the SA-only step beyond the selected values gives faster but biased estimates of \(\langle|\beta|^2\rangle\), whereas the reported SA-PAL and SA-only rows remain close to the reference value \(\langle|\beta|^2\rangle \approx 20.83\).

The numerical IACT and accuracy values corresponding to Figure 7 are reported in Table 5. Here the mixing observable is \(|\beta|^2\), matching the posterior moment whose accuracy is reported in the same table. Red rows have bias above \(2.5\%\) and are not used as accuracy-matched baselines. With this convention, SA-PAL is the best usable aggregate mixer in each fidelity regime. The SA-only steps in the last two regimes are smaller than the nominal aggressive values precisely because they were retuned to satisfy the same bias cutoff; without this retuning SA-only appears faster but no longer gives a fair accuracy-matched comparison. Retuning the comparison for individual signal coordinates \(\beta_1,\beta_2\) largely removes the SA-PAL advantage over SA-only; the claim here is therefore deliberately about the aggregate second-moment observable, and it remains sensitive to the tuning criterion.

5 Conclusions↩︎

We have introduced the AZBOBZA integrator, combining SamAdams adaptive timesteps with position-adaptive Langevin (PAL) dynamics. The scheme retains the \(O(d)\) cost per step of BAOAB/OBABO integrators while adapting simultaneously to local stiffness (via \(\eta\)) and local force direction (via the PAL friction tensor). Numerical experiments on various distinct benchmark potentials demonstrate substantially reduced integrated autocorrelation times and accurate reproduction of the canonical distribution.

The free parameters \((\alpha, \beta, \theta, h_{\mathrm{sa}})\) have interpretable roles: \(\alpha\) and \(\alpha+\beta\) set the transverse and parallel friction; \(\theta^{-1}\) sets the memory of the stiffness sensor; and \(h_{\mathrm{sa}}\) is a nominal stepsize that is automatically reduced in stiff regions. Nevertheless, the numerical experiments also show that the tuning problem is not entirely trivial. Fixed-step BAOAB is usually limited by stability, whereas SamAdams-based schemes can become limited by equilibrium accuracy at large nominal steps even while the trajectories remain stable. Consequently the meaningful comparison between SA-only and SA-PAL is at matched bias: all convergent methods have low bias when the step is small, but at larger useful adaptive steps SA-PAL is more accurate than SA-only, and when SA-only is retuned to match that accuracy SA-PAL converges substantially faster. These parameters are coupled through the discrete stability threshold, the SamAdams sensor response, and the observable used to measure mixing. In the horseshoe example, for instance, SA-PAL gives a clear advantage for the aggregate second-moment observable, while retuning around individual signal coordinates can largely remove the advantage over SA-only. Developing a principled, problem-independent tuning procedure for the combined SA-PAL scheme therefore remains an important open issue.

Future directions include: a rigorous analysis of the bias introduced by the SamAdams sensor; a systematic theory and practical protocol for tuning the coupled parameters; extension to higher-dimensional targets from Bayesian machine learning; and the combination of AZBOBZA with splitting methods for stochastic gradient estimators.

6 Langevin dynamics with a general mass matrix↩︎

Consider the underdamped Langevin system with a constant, symmetric positive-definite mass matrix \(M\) and Hamiltonian \(H(x,p) = U(x) + \tfrac{1}{2} p^T M^{-1} p\). The Ornstein–Uhlenbeck (OU) process for the momentum \(p\) is given by the SDE \[dp = -\Gamma M^{-1} p \, dt + \sqrt{2k_{\mathrm{B}}T\Gamma} \, dW_t,\] which preserves the marginal canonical density \(\pi_p(p) \propto \exp(-\tfrac{1}{2} p^T M^{-1} p / k_{\mathrm{B}}T)\) for any smooth positive-definite friction tensor \(\Gamma\). To maintain cheap computational cost when \(M \neq I\), we define the friction tensor in a way that is compatible with the kinetic energy structure: \[\Gamma(x) = \alpha M + \beta \hat{f}(x) \hat{f}(x)^T, \qquad \hat{f}(x) = \frac{\nabla U(x)}{|\nabla U(x)|}. \label{eq:Gamma-mass}\tag{15}\] The exact integration of the OU process can be performed efficiently by transforming to the decorrelated variable \(q = M^{-1/2} p\), in which the kinetic energy is simply \(\tfrac{1}{2} |q|^2\). Precomputing \(M^{1/2}\) and \(M^{-1/2}\) (which is done only once), the SDE for \(q\) becomes \[dq = -\mathcal{A}(x) q \, dt + \sqrt{2k_{\mathrm{B}}T\mathcal{A}(x)} \, dW_t, \qquad \mathcal{A}(x) = M^{-1/2} \Gamma(x) M^{-1/2}.\] Substituting 15 into the expression for \(\mathcal{A}\) yields: \[\mathcal{A}(x) = \alpha I + \beta \hat{v}(x) \hat{v}(x)^T, \qquad \hat{v}(x) = M^{-1/2} \hat{f}(x).\] This is a rank-one update to the identity matrix. Let \(w(x) = \hat{v}(x)/|\hat{v}(x)|\) be the normalized direction and \(\lambda(x) = \beta |\hat{v}(x)|^2\). Then \(\mathcal{A} = \alpha I + \lambda w w^T\), which has eigenvalues \(\alpha\) (with multiplicity \(d-1\)) and \(\alpha + \lambda\) (with multiplicity 1). The process for \(q\) is thus an isotropic OU process with a single directional enhancement, which can be solved exactly using 13 –14 with parameters \((\alpha, \lambda)\) and projection vector \(w\). The momentum is then recovered as \(p = M^{1/2} q\).

7 Conventions for the results tables↩︎

Unless stated otherwise, Tables 1, 2, 3, 4, and 5 follow the conventions below.

7.0.0.1 Mixing.

We report the integrated autocorrelation time (IACT) \(\tau\) per gradient evaluation, estimated via an automatically truncated self-consistent automatic windowing rule [22] and averaged over independent seeds (\(8\) seeds, or \(24\) for Rosenbrock). For fixed-step chains, the ACF is the usual normalized ACF of the centered observable. For SamAdams-reweighted chains, the target estimator is the self-normalized ratio \[\hat{\mu}_\phi=\frac{\sum_n w_n\phi_n}{\sum_n w_n}, \qquad w_n=\eta_n,\] so the reported ACF is computed from the ratio-estimator influence sequence \[Y_n = w_n(\phi_n-\hat{\mu}_\phi).\] This is the delta-method linearisation of the ratio estimator, due to the fact that the random denominator contributes through the centering term. The normalized ACF curves therefore show the correlation structure of \(Y_n\). For the reported scalar IACT, however, we first estimate the long-run variance \[\operatorname{LRV}(Y)=\operatorname{Var}(Y_n)+2\sum_{\ell\ge1} \operatorname{Cov}(Y_n,Y_{n+\ell})\] using the same automatically truncated window, and then report the estimator effective value \[\tau_{\mathrm{eff}}(\phi) = \frac{\operatorname{LRV}(Y)}{\bar w^2\,\widehat{\operatorname{Var}}_\pi(\phi)}, \qquad \widehat{\operatorname{Var}}_\pi(\phi) = \frac{\sum_n w_n(\phi_n-\hat{\mu}_\phi)^2}{\sum_n w_n}.\] This reduces to the ordinary IACT when \(w_n\equiv1\), while for reweighted chains it retains the denominator fluctuations through \(Y_n\) and normalizes by the target variance of the original observable rather than by \(\operatorname{Var}(Y_n)\). This convention is standard in Markov-chain error analysis for reweighted ratio estimators [23]; the long-run variance is estimated from the autocovariance sum as in standard spectral/batch-means output analysis [24]. Every integrator uses one gradient per step, so IACT per gradient is a cost-normalized measure of mixing; lower is better, and bold marks the best value in each column.

7.0.0.2 Accuracy.

The bias column is the largest relative error, over the moments \(\langle x_1\rangle, \langle x_2\rangle, \langle x_1^2\rangle\) and \(\langle U\rangle\), between the sampled estimate and the exact value obtained by quadrature; first-moment errors are measured in units of the coordinate standard deviation. On the Rosenbrock and channel marginals the quartic moment \(\langle x_2^2\rangle=\langle x_1^4\rangle\) is excluded, because at these chain lengths its estimate is dominated by sampling variance rather than discretisation error. In the ablation and horseshoe tables, red rows have bias above \(2.5\%\) and are excluded from accuracy-matched baseline comparisons. We use this same \(2.5\%\) cutoff throughout the paper when deciding whether a row is an acceptable operating point for matched-accuracy comparisons. The precise value is necessarily somewhat arbitrary, as any finite-sample bias cutoff would be. We chose it as a pragmatic compromise: small enough to rule out visibly biased aggressive parameter choices, but large enough that the required ensembles finish in a reasonable amount of computation time.

7.0.0.3 Stepsizes.

BAOAB and OBABO run at their stability-limited fixed step; SA-PAL runs at a larger nominal timestep whose effective value \(h_{\mathrm{sa}}\langle\eta\rangle\) is automatically reduced in stiff regions. All per-model parameters are listed in Section 4.2.

Acknowledgements↩︎

The work described in this article was initiated while both authors were guests of the Erwin Schrödinger Institute (ESI) at the University of Vienna in May of 2026; the work was completed during the workshop “Scalable MCMC Sampling” held at the Forschungsinstitut für Mathematik (FIM) at ETH Zürich in June of 2026. We are greatly indebted to these two institutes for making this work possible. We also thank Christoph Dellago, Jianfeng Lu, and Gilles Vilmart for helpful discussions.

References↩︎