Optimizing The Integrator Step Size for Hamiltonian Monte Carlo

Bounding The Cost of a Hamiltonian Monte Carlo Transition

Hamiltonian Monte Carlo transitions are generated from a Hamiltonian,

where the kinetic energy, $T\!\left(q,p\right)$, is specified by the user subject to certain constraints (Betancourt et al., 2014, Sec. 3.1.3) and the potential energy is defined by a given target distribution, $V\!\left(q\right)=-\log\varpi\!\left(q\right)$. Beginning with an initial position, $q$, each transition generates a joint state by randomly sampling the momenta,

and then producing a new state by integrating Hamilton’s equations,

When Hamilton’s equations are integrated exactly the joint state is distributed as

while the position is marginally distributed according to the desired target distribution,

In practice, however, the Hamiltonian trajectory can be integrated only approximately and the stationary distribution of the position will be biased away from $\varpi$.

This proposal is then accepted only with the probability,

where $\Delta_{\epsilon}\!\left(q,p\right)$ is the Hamiltonian error,

The number of attempts required to produce an accepted proposal follows a geometric distribution with the probability of success

With the cost of generating a proposal just the cost of simulating at trajectory,

the expected cost of generating an accepted proposal is given by averaging the expected number of rejections over the position space,

If $\tau$ is chosen independently of position, then

Following Beskos et al. (2013, Eq. 4.2), we apply Jensen’s inequality to the outer expectation to yield a lower bound on the expected cost. Jensen’s inequality, however, can also be applied on the inner expectation to give a complementary upper bound,

These bounds are particularly advantageous because they reduce to expectations of functions of the error in the Hamiltonian, $\Delta_{\epsilon}\!\left(q,p\right)$, with respect to the joint distribution, $\varpi_{H}$. These expectations admit well-behaved approximations independent of the actual form of the potential and kinetic energies and hence the particular details of the given problem.

Approximating Canonical Expectations

More formally, expectations with respect to the joint distribution, $\varpi_{H}$, in Hamiltonian Monte Carlo are canonical expectations and are readily estimated in practice with symplectic integrators. In this section we define canonical expectations and their relationship to Hamiltonian Monte Carlo, show how symplectic integrators approximate these canonical expectations and constrain the accuracy of these approximations in general, and then ultimately construct universal approximations to canonical expectations of certain functions of the Hamiltonian error.

This construction is necessarily technical and requires a strong familiarity with differential geometry and the geometric foundations of Hamiltonian Monte Carlo (Betancourt et al., 2014). We reserve the detailed proofs to Appendix A and suggest that readers interested in only the final result skip ahead to Section 3.

On a Hamiltonian system the symplectic form, $\omega$, immediately defines a canonical volume form,

known as canonical distributions. We refer to expectations of functions with respect to canonical distributions as canonical expectations.

The Hamiltonian foliates the manifold, $M$, into level sets,

and $\varpi$ naturally disintegrates into microcanonical distristributions, $\pi_{H^{-1}(E)}$, that concentrate on these submanifolds,

Combined with the symplectic form, the Hamiltonian also generates a Hamiltonian flow,

under which both the symplectic volume form and Hamiltonian, and consequently the canonical and microcanonical distributions, are invariant.

1.2 Hamiltonian Monte Carlo and Canonical Expectations

Because it preserves the canonical distribution, Hamiltonian flow can be used to construct an efficient Markov transition. The only problem is that a given probability space, $\left(Q,\mathcal{B}\!\left(Q\right),\varpi\right)$, does not have the symplectic structure necessary be a Hamiltonian system.

Hamiltonian Monte Carlo leverages Hamiltonian flow by considering not the sample space, $Q$, but rather it’s cotangent bundle, $T^{*}Q$. If $Q$ is a smooth and orientable $n$-dimensional manifold then the cotangent bundle is itself a smooth, orientable $2n$-dimensional manifold with a canonical fiber bundle structure, $\pi:T^{*}Q\rightarrow Q$, and a canonical symplectic form, $\omega$.

The target measure on $Q$, given in canonical coordinates by

is lifted onto the cotangent bundle with the choice of a disintegration,

this lift defines a Hamiltonian system, $\left(T^{*}Q,\omega,H\right)$, where the joint measure $\varpi_{H}$ is exactly the unit canonical measure with $\beta=1$. In particular, expectations with respect to $\varpi_{H}$ are all canonical expectations.

1.3 Computing Canonical Expectations

can be computed by taking expectations with respect to the microcanonical distributions on the level sets,

In particular the microcanonical expectations are readily computed using the Hamiltonian flow. Birkhoff’s ergodic theorem (Petersen, 1989) states that given certain ergodicity conditions the expectation of any function with respect to the microcanonical distribution is equal to its expectation along the Hamiltonian flow,

2 Approximating Canonical Expectations with Symplectic Integrators

The only problem with using Hamiltonian flow to compute expectations is that the Hamiltonian flow itself requires the solution to a system of $2n$ first-order ordinary differential equations. For all but the simplest systems, analytical solution are unfeasible and we must instead resort to simulating the flow numerically.

Fortunately, there exist a family of numerical integrators that leverage the underlying symplectic geometry to conserve many of the properties of the exact flow (Hairer, Lubich and Wanner, 2006; Leimkuhler and Reich, 2004). These symplectic integrators exactly preserve the symplectic volume form with only small variations in the Hamiltonian along the simulated flow.

In fact, symplectic integrators simulate some flow exactly, just not the flow corresponding to $H$. Using backwards error analysis one can show that a $k$-th order symmetric symplectic integrator exactly simulates the flow for some modified Hamiltonian, given by an even, asymptotic expansion with respect to the integrator step size, $\epsilon$,

Because it is exponentially small in the step size, the asymptotic error is typically neglected and the leading-order behavior of is given by

As in the exact case, the modified Hamiltonian foliates the manifold and we can define level sets,

and a corresponding transverse vector field,

Provided that the asymptotic error is indeed negligible and the symplectic integrator is topologically stable (McLachlan, Perlmutter and Quispel, 2004), the modified level sets will have the same topology as the exact level sets. In particular, when the exact foliation defines a well-behaved disintegration into microcanonical distributions we can define a corresponding modified density of states,

Using the flow from a numerical integrator to compute averages yields expectations with respect to these modified measures,

The ultimate utility of a symplectic integrator and its modified Hamiltonian system is in the accuracy of its expectations relative to the true canonical expectations. Fortunately, the geometric structure of symplectic integrators ensures that the approximation error of both microcanonical and canonical expectations computed with a symplectic integrator is well-controlled.

Note that, when $H^{-1}\!\left(E\right)$ and $\widetilde{H}^{-1}\!\left(E\right)$ intersect at some initial point, this reduces to the calculation in Arizumi and Bond (2012). Moreover, to leading-order we can replace to the expectations over $\widetilde{H}^{-1}\!\left(E\right)$ on the RHS with expectations over $H^{-1}\!\left(E\right)$ to give

This is convenient for numerical experiments when the canonical expectations can be computed analytically.

Given the decomposition of the canonical distributions over level sets, the result for the accuracy of microcanonical expectations immediately carries over to a canonical expectations.

3 Approximating Canonical Expectations of the Hamiltonian Error

The expectations necessary for bounding the cost of a basic Hamiltonian Monte Carlo transition are not just any canonical expectations but canonical expectations of functions of the Hamiltonian error,

is the Metropolis proposal. By constraining how the moments and then the cumulants of the Hamiltonian error scale with the symplectic integrator step size we can construct universal approximations to these particular expectations.

Let $\left(M,\omega,H\right)$ be a Hamiltonian system and consider a $k$-th order symmetric symplectic integrator with the corresponding modified Hamiltonian $\widetilde{H}=H+\epsilon^{k}\,G+\mathcal{O}(\epsilon^{k+2})$ and Metropolis proposal $\Phi^{\widetilde{H}}_{\epsilon,\tau}$. If the integrator is topologically stable and the asymptotic error is negligible, then the moments of the Hamiltonian error with respect to the unit canonical distribution scale as

The mean of the Hamiltonian error is particularly interesting because it can be computed analytically (Appendix A.3). In the case of a Gaussian target distribution, a Euclidean kinetic energy, and a second-order leapfrog integrator we have the Hamiltonian,

the sub-leading contribution to the modified Hamiltonian,

in agreement with numerical experiments (Figure 1).

3.2 Cumulants of the Hamiltonian Error

Let $\left(M,\omega,H\right)$ be a Hamiltonian system and consider a $k$-th order symmetric symplectic integrator with the corresponding modified Hamiltonian $\widetilde{H}=H+\epsilon^{k}\,G+\mathcal{O}(\epsilon^{k+2})$ and Metropolis proposal $\Phi^{\widetilde{H}}_{\epsilon,\tau}$. If the integrator is topologically stable and the asymptotic error is negligible, then the cumulants of the Hamiltonian error with respect to the unit canonical distribution scale as

Let $\left(M,\omega,H\right)$ be a Hamiltonian system and consider a $k$-th order symmetric symplectic integrator with the corresponding modified Hamiltonian $\widetilde{H}=H+\epsilon^{k}\,G+\mathcal{O}(\epsilon^{k+2})$ and Metropolis proposal $\Phi^{\widetilde{H}}_{\epsilon,\tau}$. If the integrator is topologically stable and the asymptotic error is negligible, then the cumulant generating function of the Hamiltonian error vanishes

Let $\left(M,\omega,H\right)$ be a Hamiltonian system and consider a $k$-th order symmetric symplectic integrator with the corresponding modified Hamiltonian $\widetilde{H}=H+\epsilon^{k}\,G+\mathcal{O}(\epsilon^{2k})$ and Metropolis proposal $\Phi^{\widetilde{H}}_{\epsilon,\tau}$. If the integrator is topologically stable and the asymptotic error is negligible, then to leading-order in $\epsilon$ the first two cumulants of the unit canonical distribution satisfy

or taking $t=1$ and appealing to Lemma 5,

But from Lemma 4 we know that to leading-order only the first two cumulants contribute,

Explicitly introducing the scaling from Lemma 4 gives

At this point we can note that if the target distribution composes into $d$ independently and identically distributions components then the cumulants scale as

If we scale the step size as $\epsilon^{2k}=\epsilon_{0}^{2k}/d$ then the cumulants scaling becomes

Consequently, in the infinite limit $d\rightarrow\infty$ all of the cumulants beyond second-order vanish, $\varpi_{\Delta_{\epsilon}}$ converges to a $\mathcal{N}\!\left(-\frac{1}{2}\alpha\epsilon^{2k},\alpha\epsilon^{2k}\right)$ in distribution, and the desired expectations simply to Gaussian integrals. Extending this argument to independently but not necessarily identically distributed distributions corresponds to the results in Beskos et al. (2013) generalized to any symplectic integrator.

Fortunately, even outside of the limit of infinite independently distributed distributions the expectations are remarkably well-behaved.

3.3 Expectations of the Hamiltonian Error

Together, these Lemmas imply that canonical expectations of any smooth function of the Hamiltonian error, as well as the Metropolis acceptance probability with its single cusp, are well-approximated by straightforward Gaussian integrals.

Let $\left(M,\omega,H\right)$ be a Hamiltonian system and consider a $k$-th order symmetric symplectic integrator with the corresponding modified Hamiltonian $\widetilde{H}=H+\epsilon^{k}\,G+\mathcal{O}(\epsilon^{k+2})$ and Metropolis proposal $\Phi^{\widetilde{H}}_{\epsilon,\tau}$. If the integrator is topologically stable and the asymptotic error is negligible, then the expectation of any smooth function of the Hamiltonian error is given by

Let $\left(M,\omega,H\right)$ be a Hamiltonian system and consider a $k$-th order symmetric symplectic integrator with the corresponding modified Hamiltonian $\widetilde{H}=H+\epsilon^{k}\,G+\mathcal{O}(\epsilon^{k+2})$ and Metropolis proposal $\Phi^{\widetilde{H}}_{\epsilon,\tau}$. If the integrator is topologically stable and the asymptotic error is negligible, then the expectation of the Metropolis acceptance probability is given by

Approximating Bounds and the Step Size Optimization Criterion

The approximation expectations in Theorem 7 and 8 immediately admit universal, approximation bounds on the cost of a basic Hamiltonian Monte Carlo transition, and minimizing these bounds provides a correspondingly universal strategy for tuning the integrator step size.

Approximations for the both the lower and upper bounds in (1) are given immediately by carrying out the Gaussian integrals analytically (Roberts, Gelman and Gilks, 1997),

Following previous work (Roberts, Gelman and Gilks, 1997; Beskos et al., 2013) we now consider the cost as a function of not the step size but rather the average acceptance probability,

and subsequently substituting into the bounds gives

Provided that the approximations hold, we can determine an optimal average acceptance probability, and hence a criterion for tuning the integrator step size, by minimizing these bounds. The dependence on the particular problem is isolated to the $\alpha^{k/2}$ scaling common to both bounds and hence does not effect the resulting optimimum; consequently the optimal average acceptance probability is the same for all choices of the potential and kinetic energies and hence defines a universal tuning strategy.

For example, with a second-order symplectic integrator the lower bound is minimized at $a\!\left(\epsilon\right)=0.651$ while the upper bound is minimized at $a\!\left(\epsilon\right)=0.801$. Because the bounds are relatively flat between these two optima any target acceptance probability between $0.6\lesssim a\!\left(\epsilon\right)\lesssim 0.9$ essentially yields equivalent results (Figure 2).

Limitations of the Step Size Optimization Criterion

When applying this optimization criterion we have to be careful to account for both the fundamental limitations in its construction and the possibility that the underlying assumptions may fail.

For example, although the cost function is applicable to both a constant integration time and an integration time chosen uniformly over some static distribution it is not applicable to an integration time that varies with the initial position, as would be necessary for a dynamically optimized integration time (Betancourt, 2013). Technically this precludes implementations of Hamiltonian Monte Carlo like the No-U-Turn sampler (Hoffman and Gelman, 2014), although in practice it has performed well as the default tuning mechanism for Stan (Stan Development Team, 2014a).

Similarly, the optimization criterion is only as good as the approximate bounds from which is it constructed. One source of error in these bounds are the high-order contributions beyond the Gaussian integral, although empirically these appear to be small for simple models (Figure 3). Because more complex models typically require smaller step sizes to achieve the same average acceptance probability, the higher-order contributions should continue to be negligible.

A more serious concern with the approximations is not in the neglected higher-order contributions but rather the assumption of topological stability and vanishing asymptotic error. In complex models the step sizes necessary for these conditions to hold can be much smaller than the step size motivated by the optimization strategy. Fortunately, when these conditions do not hold the integrator becomes unstable, manifesting almost immediately in numerical divergences that pull the state towards infinity and are readily incorporated into user-facing diagnostics. Consequently our initial optimization strategy can be made robust by monitoring these diagnostics and increasing the target average acceptance probability until no divergences occur (Figure 4). This more robust strategy has proven especially effective for hierarchical models that are particularly sensitive to this pathology (Betancourt and Girolami, 2015).

Conclusions and Future Work

By appealing to the underlying geometry of Hamiltonian Monte Carlo we have constructed a robust, universal scheme for the optimal tuning of the integrator step size, for for simple implementations of Hamiltonian Monte Carlo that use an approximate Hamiltonian flow to construct only a single Metropolis proposal.

The constructing of formal optimization criteria for more elaborate implementations of the algorithm, including windowed samplers (Neal, 1994), such as the No-U-Turn Sampler (Hoffman and Gelman, 2014), that subsample from each approximate trajectory, can also be placed into this geometric framework and utilize the theorems proved in this paper. Similarly, the understanding of canonical expectations we have built in this paper is applicable to Rao-Blackwellization schemes that keep all points along each approximate trajectory, using weights to correct for the error in the symplectic integrator.

Acknowledgements

We warmly thank Chris Wendl for illuminating discussions on various aspects of symplectic geometry used in the proofs and Elena Akhmatskaya for helpful comments. Michael Betancourt is supported under EPSRC grant EP/J016934/1, Simon Byrne is a EPSRC Postdoctoral Research Fellow under grant EP/K005723/1, and Mark Girolami is an EPSRC Established Career Research Fellow under grant EP/J016934/1.

A Proofs

Here we present proofs for the Lemmas and Theorems appearing in Sections 2.2 and 2.3, as well as a rigorous calculation of the average Hamiltonian error for the example in Section 2.3.

Provided that the integrator is stable and at the topologies of the level sets are the same, we can compare the two expectations by perturbing the true level set into the modified one by dragging it along $G\,\vec{v}$ (Figure 5). In particular, dragging the integrand gives

Before we can pull these terms onto $\widetilde{H}^{-1}\!\left(E\right)$ we have to relate them to the proper volume form, $\vec{u}\,\lrcorner\,\Omega$. Because

we have, to leading-order in the step size,

Pulling back onto the modified level set gives

Consequently the microcanonical expectation becomes

Let $\left(M,\omega,H\right)$ be a Hamiltonian system and consider a $k$-th order symmetric symplectic integrator with the corresponding modified Hamiltonian $\widetilde{H}=H+\epsilon^{k}\,G+\mathcal{O}(\epsilon^{k+2})$. If the integrator is topologically stable and the asymptotic error is negligible then

Taking $f=1$ in the derivation Theorem 1 we have

Substituting these into the definition of the canonical expectation gives

A.2 Approximating Canonical Expectations of the Hamiltonian Error

In order to compute moments of the Hamiltonian error we first have to be able to manipulate the Metropolis proposal operator, $\Phi^{\widetilde{H}}_{\epsilon,\tau}$. Fortunately, the operator is a diffeomorphism, $\Phi^{\widetilde{H}}_{\epsilon,\tau}:M\rightarrow M$ which admits a variety of convenient manipulations.

Because $\Phi^{\widetilde{H}}_{\epsilon,\tau}$ is a diffeomorphism we can push-forward integrals; in particular the numerator becomes

where we have used the fact that $\Phi^{\widetilde{H}}_{\epsilon,\tau}$ is idempotent. Substituting the definition of the modified Hamiltonian into the exponent and expanding then gives

But the modified Hamiltonian is invariant to the modified flow, $\widetilde{H}\circ\Phi^{\widetilde{H}}_{\epsilon,\tau}=\widetilde{H}$ and the integrals on the RHS become

Finally, we apply Corollary 10 to replace the modified normalization on the LHS with the true normalization and the necessary correction,

Now can compute the scaling of the moments of the Hamiltonian error, $\Delta_{\epsilon}$.

Let $\left(M,\omega,H\right)$ be a Hamiltonian system and consider a $k$-th order symmetric symmetric symplectic integrator with the corresponding modified Hamiltonian $\widetilde{H}=H+\epsilon^{k}\,G+\mathcal{O}(\epsilon^{k+2})$ and Metropolis proposal $\Phi^{\widetilde{H}}_{\epsilon,\tau}$. If the integrator is topologically stable and the asymptotic error is negligible, then the moments of the Hamiltonian error with respect to the unit canonical distribution scale as

First we explicitly introduce the modified Hamiltonian into error,

Because $\widetilde{H}$ is invariant to the flow of the integrator we can drag the second $\widetilde{H}$ along the flow to give

For odd $n$ there are an even number of terms in the expansion which pair up as

Even powers do not benefit from a similar cancelation and we’re left with the nominal scaling which gives

Given the moments the scaling of the cumulants immediately follows.

Let $\left(M,\omega,H\right)$ be a Hamiltonian system and consider a $k$-th order symmetric symplectic integrator with the corresponding modified Hamiltonian $\widetilde{H}=H+\epsilon^{k}\,G+\mathcal{O}(\epsilon^{k+2})$ and Metropolis proposal $\Phi^{\widetilde{H}}_{\epsilon,\tau}$. If the integrator is topologically stable and the asymptotic error is negligible, then the cumulants of the Hamiltonian error with respect to the unit canonical distribution scale as

Cumulants can be constructed from the moments via the recursion relation

and we use this relationship to proceed inductively.

Provided that the scaling holds up to $\kappa_{n-1}$, then if $n$ is even each term in the sum is a product of terms with like parity whereas if $n$ is odd then each term is the sum is a product of terms with odd parity. Consequently each term scales with $\epsilon^{km},\,m\geq n$ and to leading order $\kappa_{n}\propto\epsilon^{kn}$. The base case is confirmed immediately as $\kappa_{1}=\mu_{1}$.

Moreover, the symplectic structure provides a global constraint on the cumulants.

Let $\left(M,\omega,H\right)$ be a Hamiltonian system and consider a $k$-th order symmetric symplectic integrator with the corresponding modified Hamiltonian $\widetilde{H}=H+\epsilon^{k}\,G+\mathcal{O}(\epsilon^{k+2})$ and Metropolis proposal $\Phi^{\widetilde{H}}_{\epsilon,\tau}$. If the integrator is topologically stable and the asymptotic error is negligible, then the cumulant generating function of the Hamiltonian error vanishes

for the canonical distribution with $\beta=1$,

Pushing back the numerator against $\Phi^{\widetilde{H}}_{\epsilon,\tau}$ finally gives

In order to compute expectations of functions of the Hamiltonian error in the general case we appeal to the Gram-Charlier expansion (Barndorff-Nielsen and Cox, 1989).

Let $\left(M,\omega,H\right)$ be a Hamiltonian system and consider a $k$-th order symmetric symplectic integrator with the corresponding modified Hamiltonian $\widetilde{H}=H+\epsilon^{k}\,G+\mathcal{O}(\epsilon^{k+2})$ and Metropolis proposal $\Phi^{\widetilde{H}}_{\epsilon,\tau}$. If the integrator is topologically stable and the asymptotic error is negligible, then the expectation of any smooth function of the Hamiltonian error is given by

The Gram-Charlier series defines an expansion of a target density $\varpi\!\left(x\right)$ in terms of derivatives of a reference Gaussian density,

where $\kappa_{n}$ are the cumulants of $\pi$ and $\gamma_{n}$ are the cumulants of the Gaussian. By matching the Gaussian to the first two moments of $\varpi$, the first two terms in the sum vanish leaving

Given this expansion, expectations with respect to $x$ can be written as

Given the internal summation, expanding the exponential is no easy task. The action of each term, however, is relatively easy to deduce: because there is no $y$ dependence in the exponential, any term in the expansion will reduce to

with $C_{n}$ some product of coefficients, $\left(-1\right)^{m}\,\kappa_{m}\kappa_{2}^{-m/2}/m!$, whose order sums to $n$. Integrating this term by parts yields

Upon repeated integration by parts this eventually reduces to

Consequently $C_{n}$ will only introduce addition factors of $\kappa_{2}$,

If $f$ is smooth then the derivatives $f^{\left(n\right)}\!\left(\sqrt{\kappa_{2}}y+\kappa_{1}\right)$ will not introduce any addition factors of $\kappa_{1}$ nor $\kappa_{2}$ as we have already incorporated the contributions from the Jacobian. Consequently the first contribution from the expansion beyond unity will be given by

and to leading-order the expectation becomes

or substituting the results of Corollary 6,

At leading-order in $\epsilon$, the expectation of any smooth function of $\Delta_{\epsilon}$ becomes a straightforward Gaussian integral, equivalent to the infinite independently distributed limit.

Unfortunately the expectation in which we are mainly interested, the Metropolis acceptance probability, is not smooth because of a cusp at $\Delta_{\epsilon}=0$. The cusp introduces non-trivial boundary terms that then induce additional leading-order contributions to the expectations beyond those found in the smooth case.

Let $\left(M,\omega,H\right)$ be a Hamiltonian system and consider a $k$-th order symmetric symplectic integrator with the corresponding modified Hamiltonian $\widetilde{H}=H+\epsilon^{k}\,G+\mathcal{O}(\epsilon^{k+2})$ and Metropolis proposal $\Phi^{\widetilde{H}}_{\epsilon,\tau}$. If the integrator is topologically stable and the asymptotic error is negligible, then the expectation of the Metropolis acceptance probability is given by

In order to understand the effect of the cusp it is easiest to proceed as with Theorem 7 up until each term is integrating by parts,

On the first integration by parts the boundary term vanishes as before,

Continued applications, however, introduce nontrivial boundary contributions,

The second term is exactly the result we would have if the acceptance probability were smooth, with the contributions from the cusp isolated to the first term. In particular, the Hermite polynomials introduce terms proportional to $\sqrt{\kappa_{2}}$ and $\kappa_{2}$ so that at best we have

Although the expectation of the Metropolis acceptance probability is not equivalent to the infinite independently distributed limit, the deviations are isolated into two terms, $D_{1}$ and $D_{2}$, which admit further study. Indeed, empirically these terms appear to be small indicating that there may be a means of constraining their values in general and improving this result.

A.3 Approximating the Average Hamiltonian Error

Taking $n=1$ in Lemma 3 gives the full leading-order result for the average Hamiltonian error,

where $D$ depends on the second-order contribution to the modified Hamiltonian.

In the Gaussian case the integrals can be computed analytically and provide us with a means of validating Lemma 3 against numerical experiments. Here we consider second-order leapfrog integrators, encompassing both the explicit Stromer-Verlet integrator and the implicit midpoint integrator common in Hamiltonian Monte Carlo implementations. Here

and all higher-order contributions to the modified Hamiltonian vanish so that $D=0$.

The Gaussian target distribution induces the Hamiltonian

with the sub-leading contribution to the modified Hamiltonian given by

Computing the canonical expectation requires a transverse vector field, $\vec{v}$, and given the underlying Euclidean geometry with unit metric $\delta^{ij}$ an immediate choice is

Finally we use the fact that the action of the exact flow is simply a rotation in phase space,

References