{'fraction_non_alphanumeric': 0.06448, 'fraction_numerical': 0.02384, 'mean_word_length': 3.7462573226296376, 'pattern_counts': {'":': 0, '<': 8, '<?xml version=': 0, '>': 27, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 8, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Coupling and Convergence for Hamiltonian Monte Carlo

Based on a new coupling approach, we prove that the transition step of the Hamiltonian Monte Carlo algorithm is contractive w.r.t. a carefully designed Kantorovich (L1 Wasserstein) distance. The lower bound for the contraction rate is explicit. Global convexity of the potential is not required, and thus multimodal target distributions are included. Explicit quantitative bounds for the number of steps required to approximate the stationary distribution up to a given error are a direct consequence of contractivity. These bounds show that HMC can overcome diffusive behaviour if the duration of the Hamiltonian dynamics is adjusted appropriately.


Introduction
Markov Chain Monte Carlo (MCMC) is a family of methods to approximately sample from an arbitrary probability distribution. In conjunction with Bayesian methods, MCMC has revolutionized statistics and enabled applications of statistical inference to machine learning, pattern recognition, and artificial intelligence [2,15,31,3,7]. Much of the classical research activity related to MCMC considered techniques based on random walks. Regrettably the meandering behavior of these random walks leads to sampling methods that are slow [16,8]. Therefore, a more recent focus of MCMC research activity is to develop faster methods by overcoming this random walk or diffusive behavior. Hamiltonian Monte Carlo (HMC), in principle, provides one way to do this [9,24,28,29,5]. The basic idea is to "give the walker momentum". With this momentum, the walker intermixes periods of fast running with slower running to efficiently explore features of a probability distribution. However, beyond this intuition, the mathematical properties of HMC are not well understood.
In this paper, we use a coupling technique [21,Ch. 14] to analyze the HMC algorithm. Let µ denote a probability measure on R d with non-normalized density e −U (x) . In HMC, the function U (x) is viewed as a potential energy. The algorithm simulates a Markov chain with state space R d × R d and invariant probability measureμ with non-normalized density e −H (x,v) where H(x, v) = 1 2 |v| 2 + U (x) is the Hamiltonian function. We are only considering the first component which is a Markov chain on R d with invariant probability measure µ. In its simplest form, a transition step of HMC inputs an initial position x ∈ R d and a duration parameter T > 0, and outputs a final position by taking the following steps: Step 1. Draw an initial velocity ξ ∼ N (0, I d ).
Step 2. Run the Hamiltonian dynamics associated to the Hamiltonian function H for a duration T with initial position x and initial velocity ξ.
Step 3. Output the final position of this Hamiltonian dynamics.
We call the algorithm with this transition step exact HMC. In practice, the Hamiltonian dynamics is approximated by a geometric numerical integrator, and an accept/reject step is added to remove the bias due to time discretization error. In order to distinguish this from the idealized version above, we call the resulting algorithm numerical HMC.
Below, we introduce a new coupling between the transition steps of two copies of exact HMC, or numerical HMC, respectively. The approach we use is based on the framework introduced in [13], and the specific coupling of the velocities is strongly inspired by a recently developed coupling for second-order Langevin dynamics [14]. Essentially, the underlying idea is to couple two copies of HMC at different positions x and y by coupling their velocities ξ and η such that the event x + ξT = y + ηT happens with maximal probability, and to apply a reflection coupling to the velocities otherwise. In particular, the coupling is designed such that in the free case where U = 0, the positions after the transition step coincide with maximal probability. We leverage this contractive property of the coupling to obtain an explicit contraction rate for HMC in a specially designed Kantorovich (L 1 Wasserstein) metric.
To be more specific, we state a simpified version of one of our main results, which will later be reformulated rigorously as Corollary 2.6 -a corollary of Theorem 2.4. Let π(x, dy) denote the one-step transition kernel of exact HMC and let W 1 denote the standard L 1 -Wasserstein distance. Assuming sufficient regularity on the potential energy function U (see Assumption 2.1) including that ∇U is globally Lipschitz with Lipschitz constant L, and U is strongly convex outside a Euclidean ball of diameter R with strict convexity constant K, we prove that if then for all initial distributions ν and η, and for all n ≥ 0, More precisely, we prove that the transition kernel π is even contractive with contraction rate c w.r.t. the L 1 Wasserstein distance W ρ based on an explicit metric ρ that is equivalent to the Euclidean distance. This statement can be used to quantify the speed of convergence of HMC to equilibrium, and it also directly implies completely explicit bias and variance bounds, as well as concentration inequalities for ergodic averages, see e.g. [19]. A remarkable feature of the contraction rate c is that under our hypothesis on LT 2 , it only depends on K and R/T . Hence, if we choose T proportional to R, and assume that K and LR 2 are fixed (which excludes the possibility of high energy barriers), then the rate does not deteriorate as R increases. Noting that the Hamiltonian dynamics is run for time T during each transition step, we can conclude that a given approximation accuracy can be obtained after running the dynamics for a total time of kinetic order O(R), where R basically is the diameter of a ball where the target distribution concentrates in. On the other hand, a Random Walk based method would require a time of diffusive order O(R 2 ). Hence if T is chosen adequately then HMC can indeed overcome diffusive behaviour.
In Theorem 2.10, we extend our results to numerical HMC with a velocity Verlet integrator. The corresponding result is more involved than for exact HMC, but the bound c for the contraction rate is the same provided the time discretization step size h is chosen sufficiently small depending on the other parameters.
Several recent works have studied ergodic properties of HMC methods. In [4], geometric ergodicity has been proven for a variant of exact HMC (called randomized HMC) where the lengths of the durations of the Hamiltonian dynamics at the different transitions of the Markov chain are independent and identically distributed exponential random variables with mean T . The proof relies on Harris' theorem, which requires a (local) version of Doeblin's condition: a minorization condition for the transition probabilities at a finite time and in a compact set. Unfortunately, given the complicated form of these transition probabilities, the minorization condition involves non-explicit constants, and in particular, the dependence of the convergence rate on parameters in HMC is unclear. We remark that randomized HMC is related to Anderson's dynamics, which describes a molecular system interacting with a heat bath [1,22,11]. Convergence of Anderson's dynamics on an n-torus was proven in Ref. [11] by showing that Doeblin's condition holds. Very recently, geometric ergodicity for Hamiltonian Monte Carlo, but without explicit rates, has been shown by Durmus, Moulines and Saksman [10], cf. also [25] for a related work. Closely related to our results is a recent preprint by Mangoubi and Smith [27] that extends significantly ideas from [30]. In [27], coupling techniques are applied in order to analyse the properties of HMC in high dimension under the assumption of strict convexity of U ; see also [6] for related work on second-order Langevin dynamics. The major difference to our approach is that these works rely on synchronous couplings of the initial velocities in HMC, i.e., they set η = ξ. This simplifies the analysis considerably, but as a consequence, the couplings are contractive only if the stationary distribution is strongly log-concave. Another difference is that in [27] and [6], the coupling is only applied to the exact dynamics, whereas the numerical discretization is controlled by a perturbative approach. In contrast, the coupling introduced below is contractive both for exact and numerical HMC. Its superiority to synchronous couplings is supported both by theoretical results and by numerical simulations. In connection with [17], the coupling may also be useful to parallelize HMC.
Let us finally remark that the Hamiltonian flow is what, in principle, enables HMC to make large moves in state space that reduce correlations in the resulting Markov chain. One might hope that, by increasing the duration T further, the final position moves even further away from the initial position, thus reducing correlation. However, simple examples show that this outcome is far from assured. For example, for a standard normal target distribution, the corresponding Hamiltonian flow is a planar rotation with period 2π. It is easy to see that, if the initial position is taken from the target distribution, as T increases from 0 to π/2, the correlation between the initial and final positions decreases and for T = π/2, the initial and final positions are independent. However increasing T beyond π/2 will cause an increase in the correlation and for T = π, the chain is not even ergodic. For general distributions, it is likely that a small T will lead to a highly correlated chain, while choosing T too large may cause the Hamiltonian trajectory to make a U-turn and fold back on itself, thus increasing correlation [18]. Generally speaking the performance of HMC may be very sensitive to changes in T as first noted by Mackenzie in [26]. This sensitivity is reflected in our conditions on the duration parameter T .

Setting
We fix a potential energy function U : R d → R, and we denote by the corresponding Hamiltonian for unit mass. For our main results we impose the following regularity condition on U : Assumption 2.1. U is a function in C 4 (R d ) satisfying the following conditions: (A1) U has a local minimum at 0, and U (0) = 0.
(A2) U has bounded second, third and fourth derivatives. We set (A3) U is strongly convex outside a Euclidean ball, i.e., there exist constants Notice that (A3) implies that U has a local minimum. Hence if (A3) holds then (A1) can always be satisfied by centering the coordinate system appropriately and subtracting a constant from U . Conditions (A1) and (A2) imply that It is possible to replace (A3) by a Lyapunov type drift condition but this requires a slightly different approach that will be considered in a forthcoming work on randomized HMC. For some of the results stated below only a part of the assumptions is required.

Hamiltonian Monte Carlo
Hamiltonian Monte Carlo (HMC) is an MCMC method for approximate sampling from probability measures of the form We consider HMC as a Markov chain on R d (not on the phase space R d ×R d ). The transition step from x is given by x → X (x) with Here the duration T : Ω → R + is in general a random variable with a given distribution ν (e.g. ν = δ s or ν = Exp(λ −1 )), ξ ∼ N (0, I d ) and U ∼ Unif(0, 1) are independent random variables, and the acceptance event for a proposed transition is We will only consider the case where T ∈ (0, ∞) is a given deterministic constant. Moreover, is either the exact Hamiltonian flow or a numerical approximation of the Hamiltonian flow. The exact Hamiltonian flow is the solution of the ODE The corresponding Markov chain with transition step determined by (6), (7) and (8) is called exact HMC. Notice that for exact HMC, H(q T (x, ξ), p T (x, ξ)) = H(x, ξ). Hence all proposed transitions are accepted, and the transition step is simply given by In practice, the Hamiltonian flow has to be approximated by a numerical integrator. Here, we focus on the velocity Verlet integrator with discretization step size h > 0. In this case, φ t = (q t , p t ) is the solution of the equation The corresponding Markov chain with transition step determined by (6), (7) and (10) is called numerical HMC. Whenever h > 0 is fixed, we briefly write t and t instead of t h and t h , respectively. Since the velocity Verlet integrator does not preserve the Hamiltonian exactly, the rejection event A(x) C is not empty in general for numerical HMC. However, the rejection probability goes to 0 as h ↓ 0.
The HMC algorithm induces a time-homogeneous Markov chain on R d with transition kernel Here 1 − P [A(x)] is the rejection probability for a proposed transition from x. The probability measure µ defined by (5) is invariant for π, cf. e.g. [5,28].

Coupling
We now introduce a coupling for the transition steps of two copies of the HMC chain starting at different initial conditions x and y. The coupling is defined in a different way depending on whether x and y are far apart or sufficiently close.

Synchronous coupling for |x − y| ≥ 2R
The easiest way to couple the transition probabilities π(x, ·) and π(y, ·) for two states x, y ∈ R d is to use the same random variables ξ and U in both cases for the momentum refreshment and to decide whether a proposed move is accepted. The corresponding coupling transition is given by ( with A(x), A(y) defined as in (7) above. We will apply synchronous coupling for |x − y| ≥ 2R. Here we can exploit the strict convexity condition (A3) to ensure contractivity for the coupling transition.

A contractive coupling for |x − y| < 2R
For |x − y| < 2R we use a different coupling that enables us to derive a weak form of contractivity even in the absence of convexity. Let γ > 0 be a positive constant. The precise value of the parameter γ will be chosen in an appropriate way below. The coupling transition step is now given by with the event A(x) defined as in (7) above, and Here the same random variable U as in (7) is used to decide whether the proposed move to q T (y, η) is accepted. Moreover, we set where z = x − y, e = z/ |z|, ϕ 0,1 denotes the density of the standard normal distribution, and U ∼ Unif(0, 1) is independent of T , ξ and U.
This coupling is partially motivated by a coupling for second order Langevin diffusions introduced in [14]. It is defined in such a way that ξ − η = −γz holds with the maximal possible probability, and a reflection coupling is applied otherwise. As illustrated in Figure 1, the reason for this choice is that the difference process q t (x, ξ) − q t (y, η) is contracting in a time interval [0, t 0 ] if the difference ξ − η of the initial velocities is negatively proportional to the difference of the initial positions.
In order to verify that (X (x, y), Y (x, y)) is indeed a coupling of the transition probabilities π(x, ·) and π(y, ·), we remark that the distribution of η is N (0, I d ) since, by definition of η in (15) and a change of variables, for any measurable set B. Here a ∧ b denotes the minimum of real numbers a and b, and we have used that ϕ 0,I d (y − 2(e · y)e) = ϕ 0,I d (y) = ϕ 0,I d (−y). As a byproduct of this calculation, note also that where d TV is the total variation distance. Hence, by the coupling characterization of the total variation distance, ξ − η = −γz does indeed hold with maximal possible probability. z x y ξ η = ξ + γz x + T ξ = y + T η

Numerical Illustration of couplings
Before stating our theoretical results, we test the coupling defined by (13) numerically on the following two examples: • A multimodal distribution that is a mixture of twenty two-dimensional Gaussian distributions with covariance matrix given by the 2 × 2 identity matrix and with mean vector given by 20 independent samples from the uniform distribution over the rectangle [0, 10]×[0, 10]. The energy barriers between modes are not large. This example is adapted from [23,20]. • A banana-shaped distribution whose associated potential energy U : R 2 → R is given by the Rosenbrock function U (x, y) = (1 − x) 2 + 10(y − x 2 ) 2 . This function is highly non-convex and unimodal with a global minimum at the point (1, 1) where U (1, 1) = 0. This minimum lies in a long, narrow, banana shaped valley.
For simplicity, we apply the coupling globally and choose the step size h to integrate the Hamiltonian dynamics small enough to ensure that essentially all proposed moves are accepted. Realizations of the coupling process with T = 1 and γ = 1 are shown in Figure 2. We chose these parameters only for visualization purposes. The different components of the coupling are shown as different color dots. The insets of the figures show the distance between the components of the coupling as a function of the number of steps. Figure 3 shows the average time after which the distance between the components of the coupling is for the first time within 10 −9 . To produce this figure, we generated 10 5 samples of the coupled process for one hundred different values of the duration parameter T . We chose the coupling parameter γ equal to either T −1 , or equal to zero which corresponds to a synchronous coupling. The former choice is motivated by Figure 1.

Contractivity
We now state our main contraction bounds for the coupling introduced above. For given x, y ∈ R d let denote the coupling distance before and after the transition step. For exact HMC we set h := 0, whereas for numerical HMC, h > 0 is the discretization step size.

Contractivitiy by strict convexity
The assumed strict convexity of U outside of a euclidean ball directly implies contractivity of a transition step for exact HMC for initial values x and y that are sufficiently far apart.
Theorem 2.1 (Contractivity for exact HMC, strongly convex case). Suppose that Assumption 2.1 is satisfied, and let h = 0. Then for any x, y ∈ R d and T ∈ R + such that |x − y| ≥ 2R and LT 2 ≤ K/L, The proof is a direct consequence of Lemma 3.4 below, see Section 5. A similar result is proven in [27].
Notice that contractivity is only guaranteed for LT 2 smaller than the conditioning number K/L. Sometimes, contraction bounds for longer durations can be obtained. However, as discussed in the introduction, due to possible periodicity of the Hamiltonian flow, in general these do hot hold for arbitrary T .
with σ max ≥ σ min > 0. In this case, Assumption 2.1 is satisfied with R = 0, L = σ −2 min and K = σ −2 max . Theorem 2.1 gives a global contraction for synchronous coupling with rate KT 2 /2 provided that T 2 ≤ σ 4 min /σ 2 max . In particular, a necessary condition is that T is no greater than σ min , which avoids periodicities in the Hamiltonian dynamics [26].
Next we consider numerical HMC. We fix an upper bound h 1 > 0 for the discretization step size h. We assume that Under similar conditions as in Theorem 2.1 we obtain contractivity on average for coupled HMC transition steps: Theorem 2.2 (Contractivity for numerical HMC, strongly convex case). Suppose that Assumption 2.1 is satisfied, and fix T, R 2 , h 1 ∈ (0, ∞) such that (17) The estimated average is plotted as a function of the duration T of the Hamiltonian dynamics for γ = 0 (black) and γ = T −1 (gray). The latter choice is motivated by Figure 1. From (a), note that the minimum of the function is smaller and occurs at a smaller value of T when γ = T −1 . This difference is more pronounced in (b) because the underlying potential is highly nonconvex. The kinks in the graphs are due to statistical error.

holds. Then there exists
Moreover, for fixed K, L, M and N , h 0 can be chosen such that h −1 0 is of order . The proof is given in Section 5. Key ingredients are the bound for contractivity of the proposal in Lemma 3.4 and a bound for the probability that the proposal move gets accepted for one of the components of the coupling and rejected for the other component, cf. Theorem 3.8.

Contractivity without convexity
Even if we do not assume convexity, we can still obtain contractivity for x, y at a bounded distance if we replace r(x, y) = |x − y| by a modified metric. To this end we consider a distance function of the form where with parameters a, R 1 ∈ (0, ∞) to be specified below.
We again fix an upper bound h 1 > 0 for the discretization step size h. We now replace (17) by the more stringent assumption where Λ := 16LR 2 . Under this condition we obtain contractivity on average w.r.t. the metric ρ if the parameters γ, a and R 1 defining the coupling and the metric are adjusted appropriately. Explicitly, we set Theorem 2.3 (Contractivity for numerical HMC, general case). Suppose that Assumption 2.1 is satisfied, and fix T, R 2 , h 1 ∈ (0, ∞) such that (21) holds. Let γ, a and R 1 be given by (22), (23) and (24), respectively. Then there exists Moreover, for fixed K, L, M and N , h can be chosen such that h −1 is of order . For exact HMC, the corresponding contraction bound is valid for all x, y ∈ R d . Theorem 2.4 (Contractivity for exact HMC, general case). Suppose that Assumption 2.1 is satisfied, and fix T ∈ (0, ∞) such that LT 2 ≤ min K L , 1 4 , 1 16 Λ . Let γ, a and R 1 be given by (22), (23) and (24), respectively. Then (25) holds for any x, y ∈ R d with the contraction rate c given by (26).
The proofs of Theorems 2.3 and 2.4 are given in Section 5.
Remark 2.5 (Dimension dependence). The lower bound c for the contraction rate in Theorems 2.3 and 2.4 does not depend explicitly on the dimension. However, in applications, the parameter R may depend on the dimension which would cause an implicit (and possibly exponential) dimension dependence of c. This can not be avoided in the general setup considered here. In concrete models (e.g. perturbations of product measures), modifications of the approach might avoid a possible dimension dependence, see e.g. [12,13,32] for related results. For numerical HMC, the dimension also affects the upper bound h for the discretization step size which is relevant for the computational complexity. Under restrictive assumptions on the potential U (strong convexity and product structure), more precise results on the dimension dependence are proven in [27].
The following simple example demonstrates that our results provide applicable bounds in multimodal situations if T is adjusted appropriately: . Consider a mixture of two Gaussians with means ±2σ and variances σ 2 where σ > 0. The corresponding potential is In this case, for all |x| > σ, which allows us to choose K = 2/(3σ 2 ) and R = 2σ. Hence, the condition on LT 2 in Theorem 2.4 reduces to LT 2 ≤ 1/3072, and the rate in (26) reduces to If we choose T proportional to σ, then this rate is constant.

Quantitative bounds for distance to the invariant measure
For exact HMC, Theorem 2.4 establishes global contractivity of the transition kernel π(x, dy) w.r.t. the Kantorovich (L 1 Wasserstein) distance on probability measures ν, η on R d . Here the infimum is over all couplings γ of ν and η. Since the metric ρ is comparable to the Euclidean Distance on R d , contractivity w.r.t. W ρ immediately implies a quantitative bound on the standard L 1 -Wasserstein distance between the law of the HMC chain after n steps and the invariant probability measure µ.
Corollary 2.6. Suppose that Assumption 2.1 is satisfied, and let T ∈ (0, ∞) such that Then for any n ∈ N and for any probability measures ν, η on R d , where c is given by (26), and M = exp 5 2 (1 + R/T ) . In particular, for a given constant ∈ (0, ∞), the standard L 1 Wasserstein distance ∆(n) = W 1 (νπ n , µ) w.r.t. µ after n steps of the chain with initial distribution ν satisfies ∆(n) ≤ provided n ≥ 1 c The corollary is a rather direct consequence of Theorem 2.4. The proof is given in Section 6.
Remark 2.7 (Kinetic bounds). One remarkable feature of the result in Corollary 2.6 is that for a given initial error ∆(0), the number of steps required to stay below a certain error bound can be chosen universally provided T is chosen proportional to R. Notice, however, that by Condition (27), it is only possible to choose T proportional to R with a fixed proportionality constant if LR 2 is bounded by a fixed constant ! Remark 2.8 (Quantitative bounds for ergodic averages). MCMC methods are often applied to approximate expectation values w.r.t. the target distribution by ergodic averages of the Markov chain. Our results (e.g. 28) directly imply completely explicit bounds for bias and variances, as well as explicit concentration inequalities for these ergodic averages in the case of HMC. Indeed, the general results by Joulin and Ollivier [19] show that such bounds follow directly from an L 1 Wasserstein contraction w.r.t. an arbitrary metric ρ, which is precisely the statement shown above.
We now return to numerical HMC. Here, our main result in Theorem 2.3 only establishes contractivity w.r.t. W ρ on a ball of given radius R 2 . In order to derive bounds for the distance to the invariant measure of the law after n steps, we additionally have to control exit probabilities from the ball. This is achieved by a Lyapunov bound that we first state in a general form. Suppose that π(x, dy) is the transition kernel of a Markov chain on a complete separable metric space (S, ρ), and let W ρ denote the corresponding Kantorovich distance on probability measures on S.
(C3) Local contractivity: There are a measurable map (X , Y ) : S × S × Ω → S × S defined on a probability space (Ω, A, P ) and a constant c ∈ (0, ∞) such that for any x, y ∈ S, (X (x, y, ·), Y (x, y, ·)) is a realization of a coupling of π(x, ·) and π(y, ·) satisfying The proof of the following theorem is given in Section 6.
Theorem 2.9. Suppose that Assumption 2.2 is satisfied. Then for any n ∈ N and for any probability measures ν, η on (S, B(S)), where Theorem 2.9 can be applied to bound the distance to the invariant measure µ after n steps of numerical HMC. Suppose that π is the corresponding transition kernel for a given discretization step size h > 0, and let ρ denote the metric on S = R d defined by (19), (20), (23) and (24). We will then apply Theorem 2.9 with Lyapunov functions of the form ϕ(x) = 2T d 1/2 + |x| and ψ(x) = exp(U (x) 2/3 ). As a consequence of Theorem 2.3 and Theorem 2.9, we can prove that a similar number of steps as for exact HMC is also sufficient for an approximation of the invariant measure by numerical HMC, provided h is chosen sufficiently small. Theorem 2.10. Suppose that Assumption 2.1 is satisfied. Let T, h 1 ∈ (0, ∞) such that (21) holds, let ν be a probability measure on R d , and let ∆(n) = W 1 (νπ n , µ) denote the standard L 1 Wasserstein distance to the invariant probability measure after n steps with initial distribution ν. Let ∈ (0, ∞) and n ∈ N such that n ≥ 1 c where c is given by (26). Then there exists h > 0 depending only on K, L, M , N , R, T , d, ν and n, such that for any h ∈ (0, min(h , h 1 )) with T /h ∈ Z, Furthermore, h can be chosen such that for fixed values of K, L, M , N , h −1 is of order where A(ν) = log exp(U 2/3 ) dν.

A priori estimates
In this section we state several bounds for the Hamiltonian flow, for the coupling, and for acceptance-rejection probabilities that will be crucial in the proof of our main result. The proofs of all the results stated in this section are included in Section 4.

Bounds for the Hamiltonian flow and for velocity Verlet
In the following, we consider t ∈ [0, ∞) and h ∈ [0, 1] such that t/h ∈ Z if h > 0. We assume throughout that Assumption 2.1 is satisfied, and Recall that φ t = (q t , p t ) denotes the Hamiltonian flow for h = 0, and the flow of the velocity Verlet integrator for h > 0. The proofs of the following statements are provided in Section 4.
In particular, Remark 3.3. The lemma shows that on sufficiently short time intervals, the first variation of velocity Verlet can be controlled by that of the corresponding motion with constant velocity. In particular, contractivity for small times holds if u − v = −γ(x − y) for some γ > 0.
We will show next that in the region of strict convexity, the bounds in Lemma 3.2 can be improved if the initial velocities coincide. For such initial conditions, (40) and (41) imply For |x − y| ≥ 2R, the bound in (42) can be improved considerably: Lemma 3.4. There exists a finite constant C ∈ (0, ∞), depending only on L and M , such that the bound holds for any t, h as above such that and for any x, y, v ∈ R d such that Remark 3.5. The lemma does not provide a bound if |v| is very large. However, in this case we still have the upper bound that follows from (42). Hence if |v| is large with small probability, then we still get a contraction on average. For the exact Hamiltonian dynamics, there is no corresponding restriction on |x| and |v|. Here, the lemma immediately yields a contraction result for synchronous coupling.
In the case of the exact Hamiltonian flow, i.e. for h = 0, we have We are now going to quantify the error in (48) in the case where the exact flow is replaced by the flow of the velocity Verlet integrator. This is crucial to quantify the acceptance-rejection probabilities.
Lemma 3.6. There exist finite constants C 1 , C 2 ∈ (0, ∞) that depend only on L, M and N such that the bounds hold for any x, v, z, w ∈ R d and t, h as above satisfying (35).

Bounds for acceptance-rejection probabilities
We now provide some crucial bounds for probabilities and expectations that involve acceptance-rejection events and the coupling. Recall that the coupling that we consider for |x − y| < 2R ensures that ξ − η = −γz with the maximal possible probability, where z = x − y. The following lemma enables us to control probabilities and expectations when ξ − η = −γz.
Lemma 3.7. For any p ≥ 1 there exist finite constants C p and C p such that for any choice of γ, The a priori bounds (49) and (50) for velocity Verlet can be used to obtain a rather precise control for the rejection probabilities in HMC, and for the probability that in a coupling for HMC, the proposal is accepted for one component and rejected for the other. The resulting bounds are crucial to prove contractivity on average for the coupling.

Proofs of a priori bounds
If h > 0 then we define t = t h and t = t h by (11). For h = 0 we set t = t = t. In both cases, (q t , p t ) solves (10). ) and v s = p s (x, v). By (10), we have for any s ∈ [0, t] that By (4) and since t ∈ hZ, By (35), we obtain: We now derive bounds for v s . By (10) and (4), Since t ∈ hZ, we obtain by (62)  Proof of Lemma 3.2. The proof can be carried out in a similar way to the proof of Lemma 3.1, where instead of (4), we directly apply the Lipschitz bound |∇U (x) − ∇U (y)| ≤ L|x − y| for x, y ∈ R d .
Proof of Lemma 3.4. Notice that we are in the case where the initial velocities coincide. We fix x, y, v ∈ R d such that (46) holds true and set In particular, z 0 = x − y and w 0 = 0. Let By Lemma 3.2, we have w t ≤ Ltz t ≤ 2Lt for any t ∈ hZ + s.t. L(t 2 + ht) ≤ 1.
The following computations are valid for t ∈ R + such that |z s | ≥ R for s ∈ [0, t].
Recall that by (10), Let a(t) := |z t | 2 and b(t) := 2z t · w t . Our goal is to derive an upper bound for a(t). To this end we note that a and b satisfy the following differential equations: where We see thatḃ (t) = −2Ka(t) + β(t) with a function β satisfying The initial value probleṁ has a unique solution that is given by We now bound the terms δ, and β. Note first that the assumptions imply that Kt 2 ≤ Lt 2 ≤ 1 ≤ π 2 /2. Hence t ≤ π/ √ 2K, and thus sin( √ 2K(t − r)) ≥ 0 for any r ∈ [0, t]. Moreover, by (3), Hence the first term in the second line of equation (68) is negative. Lett : In particular, for t ∈ hZ + , where z ,2 t := (z t ) 2 . Moreover, δ 2 (t) is given by and hence by (43), for t ∈ hZ + , t 0 cos( In order to control δ 3 in an efficient way note that Therefore, we obtain for t ∈ hZ + , by (43), Next, we derive bounds for β(t). We first observe that by (67) and (43), and hence The terms 1 , 2 and 3 can be controlled similarly to δ 1 , δ 2 and δ 3 . Analogously to (72), we obtain and thus by (43), for t ∈ hZ + , 3 (t) ≤ Lhw t z t ≤ 2L 2 htz ,2 t . Thus in total, we obtain by (74), By (68), (69), (70), (71) and (73), we obtain By Lemma 3.1, there is a finite constant C that depends only on K, L and M such that for any t ∈ hZ + satisfying (45), we have Noting that cos( √ 2Kt) ≤ 1 − Kt 2 + K 2 t 4 /6, and K 2 t 4 ≤ L 2 t 4 ≤ Kt 2 by (45), we obtain by (75) and (46), for any t ∈ h · Z + s.t. (45) holds. This inequality then implies for t as before. Indeed, suppose first that h > 0. Then (79) follows directly from (78) if |z t | ≤ |z 0 | holds for any t > 0 satisfying (45). Now suppose for a contradiction that there exists t > 0 s.t. (45) holds and |z t | > |z 0 |. Since z t is linear on each partition interval, we may assume that t ∈ hZ + . Let t 0 denote the smallest s ∈ hZ + for which |z s | > |z 0 |. Then z t0 = |z t0 |, and hence by (78), |z t0 | ≤ |z 0 | in contradiction to the definition of t 0 . Thus (79) holds for all t as above. For h = 0, we can argue similarly by the intermediate value theorem.
Summarizing, we have shown that (79) holds for t ∈ hZ + satisfying (45) provided |z s | ≥ R for all s ∈ [0, t] and h satisfies (46). To conclude the proof suppose that |z 0 | ≥ 2R. We claim that then |z t | ≥ R holds for all t satisfying (45). Indeed let t 1 := inf{t : |z t | < R}. Then |z s | ≥ R on [0, t 1 ]. Suppose for a contradiction that L(t 2 1 + t 1 h) < K/L ≤ 1. Then by (79), |z s | ≤ |z 0 | for s ∈ [0, t 1 ]. Hence by (65) and (64), and thus |z t1 | > 1 2 |z 0 | ≥ R in contradiction to the definition of t 1 . Proof of Lemma 3.6. Fix x, v, z, w ∈ R d . We set where where In particular, for any t ∈ hZ + , t 0 I a s ds = 0, and Similarly, we obtain where In particular, for any t ∈ hZ + , t 0 III a s ds = 0, and By combining the bounds, we obtain for t ∈ hZ + : This implies the first claim (49), since for t ∈ hZ + satisfying (35), both x t and v t are bounded by a constant multiple of max(|x|, |v|).
Next, we consider the derivative flow where the derivatives are taken w.r.t. the initial condition. We have with initial condition (x 0 , v 0 ) = (z, w). In particular, for We now first derive a priori bounds for x t and v t . By (90) and (91), Similarly, by (91) and (92), Now we can derive bounds for H t . We have Similarly as above, we bound the terms I t , II t , III t and IV t individually. By (80), Similarly to the decomposition of I t above, we have For t ∈ hZ + , t 0 VII a s ds = 0. Moreover, similarly to (85) and (86), we have Furthermore, by (90), we have For the second term we have II t = IX t + X t + XI t where For t ∈ hZ + , we obtain by (90) and (91), Furthermore, III t = (III a t ) + (III b t ) + (III c t ) . By (87) and the chain rule, t 0 (III a t ) ds = 0 for t ∈ hZ + . Moreover, by (88) and the chain rule, Finally, a similar computation as for (III b t ) shows that Collecting all the bounds derived above, we eventually obtain for t ∈ hZ, where Q 1 and Q 2 are explicit quadratic forms. We can conclude that Proof of Lemma 3.7. Let p = 0 or p ≥ 1. Then by definition of η, For p = 0, we directly obtain (51), and for p ≥ 1, 2π, 2pm p−1 ) with m p denoting the p th moment of the standard normal distribution. Finally, since ξ ∼ N (0, I d ) and the event {ξ − η = −γz} is measurable w.r.t. σ(e · ξ), we obtain Proof of Theorem 3.8. Recall that Therefore, and since U is independent of ξ, we obtain by Lemma 3.6, SinceÂ(y) is defined similarly to A(x) with x, ξ replaced by y, η and η ∼ ξ, we obtain corresponding bounds for P [Â(y) C |η] and P [Â(y) C ].
Next, we derive the corresponding bounds for the probabilities that the proposed move is accepted for one of the two components and rejected for the other. By independence of U from ξ and η and by Lemma 3.6, we have where x u = ux + (1 − u)y, ξ u = uξ + (1 − u)η, z = x − y and W = ξ − η. This proves (60), and (58) can be shown similarly with η replaced by ξ.

Proofs of main results
Proof of Theorem 2.1. We fix x, y ∈ R d such that |x−y| ≥ 2R and max(|x|, |y|) ≤ R 2 . Since synchronous coupling is applied for |x−y| ≥ 2R, we have η = ξ. Hence by (40) and Lemma 3.4 with h = 0, we obtain Proof of Theorem 2.2. Without loss of generality, we may assume that R 2 is chosen sufficiently large such that We fix x, y ∈ R d such that |x − y| ≥ 2R and max(|x|, |y|) ≤ R 2 . Since synchronous coupling is applied for |x − y| ≥ 2R, we have η = ξ, and hence Moreover, on A(x) ∩ A(y) C we have Y = y, and thus Similarly, on A(x) C ∩ A(y), Therefore, we obtain In order to control the first term, we choose a constant C ∈ (0, ∞) as in Lemma 3.4, and we assume h ≤ min(h 1 , h 2 ) where h 2 := K C(1+2R2) . Then by Lemma 3.4, Therefore, by (47), and since K ≤ L, Therefore, by Theorem 3.8, for |x| , |y| ≤ R 2 , We choose h 3 > 0 such that for h ≤ h 3 , the expression on the r.h.s. is smaller than 1/5. Because of (17), this can be achieved with h −2 3 of order O((R 3 2 +d 3/2 )). For h ≤ h 3 , we obtain Therefore, and by (99), In order to control II, we note that by Lemma 3.1, Hence in this case we obtain A corresponding bound with x and y interchanged holds for III. Hence by the bound (58) for the conditional AR probability given ξ, We choose a strictly positive constant h 4 such that for h ≤ h 4 , the right hand side is smaller than 1 24 KT 2 r. By (17), this can be achieved with h −1 Proof of Theorem 2.3. The parameters γ, a and R 1 have been chosen in (22), (23) and (24) such that the following conditions are satisfied: aT ≥ 1, Indeed, (103) and (105) hold by (22), (104) holds by (22) and (21), (106) holds by (23), (107) holds by (24) and (22), and (108) holds by (24) and (23). The bounds (103)-(108) will be essential in the following arguments. We have chosen γ and a as large resp. small as possible such that (103), (105) and (106) hold. Then (104) implies the additional constraints on T in (21), and R 1 is chosen such that (107) and (108) are satisfied.
To prove contractivity, we fix x, y ∈ R d such that max(|x|, |y|) ≤ R 2 . Since x and y are fixed, we briefly write r and R instead of r(x, y) and R (x, y). We consider separately the cases where |x − y| ≥ 2R and |x − y| < 2R.
(i) Contractivity for |x − y| ≥ 2R. For |x − y| ≥ 2R, we can apply the result of Theorem 2.2. Indeed, choose h 0 as in Theorem 2.2. Then for h ≤ h 0 , by concavity of f and by (18), where the lower bound c 1 for the contraction rate is given by Recall that f is concave with f (0) = 0, and f is linear for r ≥ R 1 . Hence the function r → rf (r)/f (r) attains its minimum at R 1 , where f (R 1 ) = e −aR1 and f (R 1 ) = R1 0 e −as ds ≤ min(R 1 , a −1 ). Therefore, by (23) and (24), (ii) Contractivity for |x − y| < 2R. For |x − y| < 2R, we apply the coupling defined by (13) and (14). where Only the first term is responsible for contractivity. The other terms are perturbations that have to be controlled. We will now derive upper bounds for each of the four terms. We remark at first that on A(x) ∩Â(y), by Lemma 3.2.
Here we have used that by (15), Moreover, we have used that by (107), R 1 ≥ 5 4 (1 + γT )r. By concavity of f and by (108), we obtain IV. By a similar argument as in the proof of Theorem 2.2, we obtain Here we have used (61) in the last step. By concavity of f we obtain for h ≤ h 6 where h 6 is a positive constant depending only on L, M , N , R, R 2 and d that by (21) can be chosen such that h −2 6 is of order O((1 + R)(R 4 2 + d 2 )).
Combining the bounds for the terms I, II, III and IV in (113), (114), (115) and (116), we obtain for h ≤ min(h 5 , h 6 ), where the contraction rate c 2 satisfies  (110) and (117), we see that for h ≤ min(h 1 , h ) and for any x, y ∈ R d with max(|x|, |y|) ≤ R 2 , where c := min(c 1 , c 2 ). Moreover, by (111) and (118), Proof of Theorem 2.4. The contraction bound for exact HMC can be derived similarly to the proof of Theorem 2.3. In this case, instead of Theorem 2.2, we apply Theorem 2.1 in Step (i). Furthermore, the rejection events A(x) C and A(y) C , x, y ∈ R d , are empty for exact HMC. Therefore, the corresponding terms do not have to be taken into account in Step (ii). Consequently, the resulting bound (25) is valid for all x, y ∈ R d with the same rate c as above.
6. Proofs of results in Section 2.6 All bounds in Section 2.6 are based on the following observation: Lemma 6.1 (Locally contractive couplings and supermartingales). Let π(x, dy) be a Markov transition kernel on a complete separable metric space (S, ρ). Suppose that there exist a constant c ∈ (0, ∞), a measurable subset A ⊆ S, a probability space (Ω, A, P ), and a measurable map from S × S × Ω to S × S such that for any x, y ∈ S, (X (x, y), Y (x, y)) is a realization of a coupling of π(x, ·) and π(y, ·), and Then, for any probability measure γ on S×S, there is a Markov chain (X n , Y n ) n≥0 defined on a probability space ( Ω, A, P ) such that (X 0 , Y 0 ) ∼ γ, both marginal processes (X n ) n≥0 and (Y n ) n≥0 are Markov chains on S with transition kernel π, and such that the process is a non-negative supermartingale w.r.t. the filtration generated by (X n , Y n ) n≥0 .
Now let (X n , Y n ) n≥0 be a time-homogeneous Markov chain on a probability space ( Ω, A, P ) with initial distribution (X 0 , Y 0 ) ∼ γ and transition kernel k, and let F n = σ((X i , Y i ) : 0 ≤ i ≤ n). Then for any n ≥ 0, holds P -almost surely on {(X n , Y n ) ∈ A × A}. Therefore, the process (M n ) defined by (120) is a non-negative (F n )-supermartingale.
The error bound for exact HMC in Corollary 2.6 is a direct consequence of Theorem 2.4 and Lemma 6.1: Proof of Corollary 2.6. For exact HMC, by Theorem 2.4, the local contractivity condition (119) in Lemma 6.1 is satisfied for S = A = R d , ρ and c given by (19) and (26), and the coupling (X (x, y), Y (x, y)) introduced above. Now let ν and η be probability measures on R d , and let γ be an arbitrary coupling of ν and η. Then by Lemma 6.1, there is a Markov chain (X n , Y n ) n≥0 on a probability space ( Ω, A, P ) such that (X 0 , Y 0 ) ∼ γ, both (X n ) and (Y n ) are Markov chains with transition kernel π and initial laws ν and η, respectively, and M n = e cn ρ(X n , Y n ) is a non-negative supermartingale. Hence for any n ∈ N, Taking the infimum over all couplings γ ∈ Π(ν, η), we see that (28) holds. Furthermore, by (19) and (20), Therefore, (28) implies W 1 (νπ n , ηπ n ) ≤ e aR1 e −cn W 1 (ν, η).
Choosing η = µ, we have ηπ n = µ for all n. Hence The second part of the assertion now follows, because by (23) and (24), aR 1 = 5 2 (1 + R/T ). Now suppose again that π(x, dy) is an arbitrary Markov transition kernel on a complete separable metric space (S, ρ). For proving Theorem 2.9, we combine Lemma 6.1 with a Lyapunov bound for exit probabilities: Proof of Theorem 2.9. Let ν and η be probability measures on S, and let γ be a coupling of ν and η. By (C3) in Assumption 2.2, the conditions in Lemma 6.1 are satisfied with A = {ψ > C}. Hence on some probability space ( Ω, A, P ), there is a coupling (X n , Y n ) n≥0 of Markov chains with transition kernel π and joint initial law (X 0 , Y 0 ) ∼ γ such that M n = e c(n∧T ) ρ(X n∧T , Y n∧T ) is a non-negative supermartingale stopped at T = min {n ≥ 0 : ψ(X n ) > C or ψ(Y n ) > C} .
Here we have used that by (C2), for any k ≤ n; and a corresponding inequality holds for E[ϕ(Y n ); T = k].
From now on, we consider numerical HMC. Let π denote the transition kernel for a given step size h > 0, and let ρ be the metric defined by (19), (20), (23) and (24). To be able to apply Theorem 2.9, we first identify appropriate Lyapunov functions.

Condition (135) holds if and only if
log C ≥ log(u + v log C) + log w + n log(βλ).
In particular, since log(u + v log C) ≤ log + (2u) + log + (2v log C) ≤ log + u + log + v + 2 + log log C, there is a universal finite constant C 0 such that (135) is satisfied if C ≥ C 0 and log C ≥ log + u + log + v + log w + n log(βλ).