---
abstract: 'We establish general conditions under which Markov chains produced by the Hamiltonian Monte Carlo method will and will not be geometrically ergodic.  We consider implementations with both position-independent and position-dependent integration times.  In the former case we find that the conditions for geometric ergodicity are essentially a non-vanishing gradient of the log-density which asymptotically points towards the centre of the space and does not grow faster than linearly.  In an idealised scenario in which the integration time is allowed to change in different regions of the space, we show that geometric ergodicity can be recovered for a much broader class of tail behaviours, leading to some guidelines for the choice of this free parameter in practice.' 
address:
- |
    School of Mathematics\
    University Walk\
    Bristol\
    BS8 1TW\
    United Kindgom\
- |
    Department of Statistics\
    University of Warwick\
    Coventry\
    CV4 7AL\
    United Kingdom\
    \
- |
    Department of Statistical Science\
    University College\
    Gower Street\
    London WC1E 6BT\
    United Kingdom\
author:
- 
- 
- 
- 
bibliography:
- 'mybibliography.bib'
title: On the Geometric Ergodicity of Hamiltonian Monte Carlo
---

,

,

Introduction {#sec:intro}
============

This paper deals with ergodic properties of Markov chains produced by the *Hamiltonian* (or *Hybrid*) Monte Carlo method (HMC), a technique for approximating high dimensional integrals through stochastic simulation \[[@duane1987hybrid]\].  Iterative algorithms of this type are widely used in (for example) statistical inference \[[@robert2011short]\], inverse problems \[[@stuart2010inverse]\], artificial intelligence \[[@andrieu2003introduction]\] and molecular dynamics \[[@lelievre2010free]\].  In many of these settings a prior distribution can be constructed for an unknown quantity, and after conditioning on some observed data, Bayes’ theorem gives a posterior — to extract relevant information from this typically high-dimensional integrals must be evaluated. 

A popular approach to such problems is to simulate a Markov chain whose limiting distribution is the posterior, and compute long-run averages \[e.g. [ @diaconis2009markov]\].  Provided the chain is *ergodic*, then a Law of Large Numbers exists for these.  Several *Markov chain Monte Carlo* (MCMC) methods of this nature have been proposed in the literature \[e.g. [ @green2015bayesian]\], and many are well understood theoretically \[e.g. [ @roberts1996geometric; @roberts1996exponential]\].  HMC has proven an empirical success, with numerous authors noting its superior performance in a variety of settings \[e.g. [ @gelman2014bayesian; @jasche2010fast; @betancourt2015hamiltonian]\] and high performance software available for its implementation \[[@carpenter2016stan]\].  Comparatively few rigorous results, however, exist to justify this.  Indeed, the absence of such analysis has been noted on more than one occasion \[[@diaconis2013some; @diaconis2014connections]\].  The major contribution of this work is to establish general scenarios under which geometric ergodicity can and cannot be established for Markov chains produced by common HMC implementations. 

Consider a Borel space $({\mathbf{X}},{\mathcal{B}})$. In the most general case of interest to us ${\mathbf{X}}= (\mathcal{M},g)$ is some geodesically complete Riemannian manifold of dimension $d$, but in this article we focus on the case ${\mathbf{X}}= {\mathbb{R}}^d$. We define a Markov chain $(X_n)_{n \geq 0}$ on $({\mathbf{X}},{\mathcal{B}})$ through an initial distribution $\delta_{x}(\cdot)$ and a family of mappings $f_\theta:{\mathbf{X}}\to {\mathbf{X}}$, indexed by an additional random variable $\theta$ defined on the Borel space $(\Theta,{\mathcal{B}}_{\theta})$ and with law $\gamma(\cdot)$ \[e.g. [ @diaconis1999iterated]\].  A transition kernel $P: {\mathbf{X}}\times {\mathcal{B}}\to [0,1]$ can then be induced through the relation $$P(x,A) = \int \mathbbm{1}_A(f_\theta(x))\gamma(d\theta),$$ for any $A \in {\mathcal{B}}$. Constructing a Markov chain for which some distribution of interest $\pi(\cdot)$ is *invariant* is not very difficult, owing to the Metropolis–Hastings algorithm \[[@metropolis1953equation; @hastings1970monte]\], in which the family $\{f_\theta, \theta \in \Theta\}$ is given by $$f_\theta(x) := \begin{cases}
g_\xi(x) & u < \alpha(x,g_\xi(x)), \\
x & \text{otherwise},
\end{cases}$$ where $\theta = \{\xi,u\}$ in this case, with $u \sim U[0,1]$, and $\{g_\xi, \xi \in \Xi\}$ is a family of ‘candidate’ maps, with $\xi \sim \mu(\cdot)$. A candidate transition kernel is induced as $Q(x,A) = \int \mathbbm{1}_A(g_\xi(x))\mu(d\xi)$ for any $A \in {\mathcal{B}}$. The ‘acceptance probability’ $\alpha: {\mathbf{X}}\times {\mathbf{X}}\to [0,1]$ is defined via the Radon–Nikodym derivative $$\label{eqn:accept}
r(x,y) := \frac{\pi(dy)Q(y,dx)}{\pi(dx)Q(x,dy)},$$ as $\alpha(x,y) = 1 \wedge r(x,y)$. The resulting chain $(X_n)_{n \geq 0}$ is reversible with respect to $\pi(\cdot)$ \[[@tierney1994markov]\]. 

Simple choices for the family $\{g_\xi, \xi \in \Xi\}$ result in Markov chains which are intuitive and convenient to analyse.  In the random walk case $g_\xi(x) = x + \xi$, with $\Xi = {\mathbf{X}}$ and $\mu(\cdot)$ a centred, symmetric distribution \[[@tierney1994markov]\].  For the Metropolis-adjusted Langevin algorithm (MALA) $g_\xi(x) = x + h\nabla\log\pi(x)/2 + \sqrt{h}\xi$, with $\mu(\cdot)$ a standard Gaussian measure on $\Xi = {\mathbf{X}}$, $h>0$ a constant, $\nabla$ the gradient operator and $\pi(x)$ the Lebesgue density of $\pi(\cdot)$. The former is in some sense a naive choice, while the latter is an Euler–Maruyama scheme for the diffusion governed by $dX_t = \nabla\log\pi(X_t)dt + \sqrt{2}dW_t$, for which $\pi(\cdot)$ is invariant.  In both cases proposals are *local* (only depending on analytic information at the current point), and $x$ is combined with $\xi$ *linearly*, with added complexity coming only through the (typically nonlinear) $\alpha$. As a result, simple bounds on $\alpha$ allow stochastic stability properties such as $\pi$-irreducibility to be deduced straightforwardly, and rates of convergence for different forms of $\pi(\cdot)$ are also well-understood in both cases \[[@mengersen1996rates; @roberts1996geometric; @roberts1996exponential]\]. 

The HMC method can also be considered within the above framework, as outlined in [@betancourt2014geometric].  The algorithm is designed to exploit the measure-preserving properties of *Hamiltonian flow* \[e.g. [ @leimkuhler2004simulating; @lee2012symplectic]\], which can be induced provided the state space for the chain is a symplectic manifold \[e.g. [ @lee2012symplectic Chapter 22]\].  The space ${\mathbf{X}}$ can be made symplectic by doubling the dimension, introducing auxiliary *momentum* variables $p$ which follow some user-specified distribution.  A Hamiltonian function can then be constructed on the resulting *phase space* which preserves a distribution for $(x,p)$, the $x$-marginal of which will be $\pi(\cdot)$. At each step of the Markov chain, a fresh value for $p$ is drawn from its conditional distribution given the current $x$ state, and then the relevant Hamiltonian flow is approximated for $T$ units of time to produce the next proposed move.  The resulting proposal map is $$\label{eqn:hmap}
g_\xi(x) = \text{Pr}_x \circ \varphi_T(x,p),$$ where $\text{Pr}_x$ denotes the projection operator onto the $x$ coordinate, $\varphi_T$ the approximate flow for $T$ units of time, and $\xi = \{T,p\}$. Typically the distribution for $p$ is chosen to be a $d$-dimensional Gaussian.  If the law of $p$ does not depend on $x$, then the Störmer–Verlet (or *leapfrog*) numerical integrator is typically used to approximate the flow, with ${\varepsilon}>0$ chosen as the integrator step-size and $L$ the number of ‘leapfrog steps’ (meaning $T = L{\varepsilon}$).  The choice of $T$ is a point of ambiguity; often it is set to be some fixed value, however heuristics have also been suggested for choosing this dynamically \[e.g. [ @hoffman-gelman:2013]\].  For $T = {\varepsilon}$ (meaning $L = 1$) in fact HMC reduces to MALA \[e.g. [ @lelievre2012langevin]\].  In general, however, for $L > 1$ (\[eqn:hmap\]) will be a non-linear function of $p$, making analysis of the method challenging, particularly in the case of a dynamic $T$.

Our main contribution is to establish conditions under which common HMC implementations produce a geometrically ergodic Markov chain.  We also establish instances where convergence will *not* be geometric, meaning the sampler may perform poorly in practice.  We first consider the case where the choice of *integration time* $T$ is chosen independently of the current position, and show that here the non-linear terms in $g_\xi(x)$ can be *bounded in probability* as the norm $\|x\| \to \infty$ under suitable assumptions, meaning that geometric convergence essentially occurs for HMC in the same scenarios as for MALA, when the tails of $\pi(x)$ are uniformly exponential or lighter, but no lighter than that of a Gaussian density.  We then consider an idealised scheme in which $T$ is chosen as a function of the current position, and show that in this case geometrically converging chains can be constructed for a much broader class of targets.  Although the latter results are in an idealised case, they do offer some practical guidelines for the choice of integration time, which can be used to examine some commonly used heuristics in the literature as well as suggest alternatives. 

Literature review
-----------------

The HMC method was first introduced in lattice field theory \[[@duane1987hybrid]\], as a *hybrid* of two differing approaches to molecular simulation introduced in [@alder1959studies] and [@metropolis1953equation] respectively.  The algorithm was introduced to Statisticians through [@neal1996bayesian]; a user-oriented review of the method is given in [@neal2011mcmc].  Several extensions have been suggested.  A *generalized* scheme in which the momentum is only partially refreshed was introduced in [@horowitz1991generalized], and see also [@ottobre2016function].  Other extensions have been proposed to allow more directed motion and reduced rejections \[[@campos2015extra; @sohl2014hamiltonian]\].  The *shadow Hamiltonian* framework was introduced in \[[@izaguirre2004shadow; @akhmatskaya2008gshmc]\], in which the properties of symplectic integrators are exploited in order to reduce rejection rates, and in some cases this can enhance sampling efficiency.  A dynamic approach to tuning the integration time parameter was introduced through the ‘No-U-Turn Sampler’ of [@hoffman-gelman:2013], which is now implemented in the Stan software \[[@carpenter2016stan]; see also [@wang2013adaptive]\].  An extension showing how to implement the sampler on a Riemannian manifold which is *globally diffeomorphic* to ${\mathbb{R}}^d$ is given in [@girolami2011riemann] \[see also [@betancourt2013general]\], and to embedded manifolds with closed form geodesics in [@byrne2013geodesic]. 

Theoretical study of MCMC methods is in the main focused on two themes: *convergence to equilibrium* and *asymptotic variance*.  When ${\mathbf{X}}$ is of finite dimension, the first is often understood through upper bounding some suitable discrepancy between the $n$th iterate of the Markov chain and its limiting distribution, as a function of $n$. When the discrepancy is taken as either the Total Variation or $V$-norm distance (for some suitable Lyapunov function $V:{\mathbf{X}}\to [1,\infty)$), then the drift and minorisation conditions popularised in \[[@meyn2012markov]\] can be exploited to envoke Harris’ theorem \[e.g. [ @hairer2011yet]\], implying that the distance to equilibrium decreases geometrically in $n$ (we elaborate in Section \[sec:prelim\]).  If such a bound holds then for reversible chains a Central Limit Theorem exists for long-run averages of $L^2(\pi)$ functionals \[e.g. [ @roberts2004general]\].  We take this approach here.  Note that such techniques rely crucially on the chain being $\psi$-irreducible for some $\sigma$-finite measure $\psi(\cdot)$.

For HMC, [@cances2007theoretical] establish that if the *potential energy* $U(x) = -\log\pi(x)$ is bounded above, continuous and has bounded derivative then the algorithm will produce a $\pi$-irreducible chain.  The result holds for both the exact flow and the leapfrog integrator variants of HMC.  Typically the boundedness assumption on $U(x)$ will only be satisfied when ${\mathbf{X}}$ is compact.  The authors also show that $\pi$-irreducibility can be established more broadly if the integration time is chosen stochastically.  More recently, [@bou2015randomized] consider a continuous-time version of HMC in which the integration step-size is randomly sampled from an Exponential distribution.  Under the assumption that Hamilton’s equations can be exactly integrated, they prove that the algorithm will produce a geometrically ergodic Markov chain whenever the tails of $\pi(x)$ decay at a Gaussian rate or faster.  The method of the authors is to relate HMC to *underdamped* Langevin dynamics, the ergodic properties of which are established in [@mattingly2002ergodicity] \[see also [@bou2010long]\].  By contrast, we relate HMC to *overdamped* Langevin dynamics, as analysed in [@roberts1996exponential].  Although at first this may seem less natural, in fact it allows us to paint a broad picture of when the algorithm as used in practice will and will not produce a geometrically ergodic Markov chain. [ @seiler2014positive] offer some practical approximations for convergence bounds under a *positive curvature* assumption on the underlying chain.  We discuss these further in Section \[sec:discussion\]. 

Asymptotic variances of long-run averages from Markov chains are often considered via analysing the expected squared jump distance $\mathbb{E}[(X_{i+1} - X_i)^2]$; at equilibrium this can then be optimised over the various parameters of the dynamics.  Careful study of this quantity can also indicate how algorithm performance depends on $d$. In the case of HMC such analysis has been performed \[[@beskos2013optimal; @betancourt2014optimizing]\], suggesting that the method scales more favourably than other approaches with dimension, and a larger optimal acceptance rate is attained. 

Recently a class of MCMC methods have been generalised to the context of sampling on spaces of infinite dimension \[e.g. [ @cotter2013mcmc]\], principally in the case where the prior is a Gaussian measure \[though see [@vollmer2015dimension]\].  The HMC method has been generalised in this way \[[@beskos2011hybrid]\].  Due to the frequent singularity of measures in such spaces, it is often necessary to characterise distance to equilibrium through more general Wasserstein metrics than the Total Variation case in these settings \[e.g. [ @hairer2014spectral]\].  Analysing this class of methods is beyond the scope of this paper, though we note that the approaches of [@eberle2014error] and [@durmus2015quantitative] in the context of MALA are a useful pre-cursor in this direction. 

Outline
-------

In the next section we give an overview of the key findings of this work.  After introducing mathematical preliminaries in Section \[sec:prelim\], in Section \[sec:hmc\] we discuss Hamiltonian Monte Carlo.  We consider the chain on ${\mathbf{X}}$-space rather than the usual phase space representation.  In Section \[sec:fixed\] we consider position-independent integration times, showing cases in which the resulting Markov chain will and will not be geometrically ergodic.  We then consider position-dependent integration times in Section \[sec:dynamic\], and present an idealised scenario in which geometrically ergodic chains can be recovered in a much broader class of settings.  After this we discuss how these results relate to practical schemes currently in use, can inform new modifications to these, and give general guidelines for the choice of integration time in Section \[sec:discussion\].  We also discuss how our results can be extended in various ways. 

Notation
--------

Let $({\mathbf{X}},{\mathcal{B}})$ denote a Borel space.  Here we restrict attention to ${\mathbf{X}}= {\mathbb{R}}^d$ (and write $\|x\|$ for the Euclidean norm of $x \in {\mathbf{X}}$).  For functions $f,g:{\mathbb{R}}_{\geq 0} \to {\mathbb{R}}_{\geq 0}$ let $f(x) \asymp g(x)$ mean that $\lim_{x \to \infty} f(x)/g(x) = c$ for some $c<\infty$. Throughout let $\pi(\cdot)$ be a finite ‘target’ measure, and $\pi(x)$ the corresponding Lebesgue density for some $x \in {\mathbf{X}}$, and let $\mathfrak{L}(\cdot)$ be a distribution defined over $\mathbb{Z}_+$ We will denote Lebesgue measure on ${\mathbb{R}}^d$ with $\mu^L(\cdot)$, the Dirac point mass at $x$ with $\delta_x(\cdot)$ and the standard Gaussian measure with $\mu^G(\cdot)$. We write $P: {\mathbf{X}}\times {\mathcal{B}}\to [0,1]$ to denote a Markov transition kernel, meaning $P(x,\cdot)$ is a probability measure for any $x \in {\mathbf{X}}$ and $P(\cdot,A)$ is measurable for any $A \in {\mathcal{B}}$. $P$ acts to the left on measures through $\mu P(\cdot) := \int \mu(dx)P(x,dy)$ and to the right on functions through $Pf(x) := \int f(y)P(x,dy)$. We let $P^n(x,\cdot) := \int P^{n-1}(x,dy)P(y,\cdot)$ and say $\pi(\cdot)$ is *invariant* for $P$ if $\pi P(\cdot) = \pi(\cdot)$.

Denote the Total Variation distance between two distributions $\mu(\cdot)$ and $\nu(\cdot)$ on $({\mathbf{X}},{\mathcal{B}})$ as $\|\mu(\cdot) - \nu(\cdot)\|_{TV} := \sup_{|f| \leq 1}|\mathbb{E}_\mu f - \mathbb{E}_\nu f|$. We say $\pi(\cdot)$ is a *limiting* distribution for $P$ if $\|P^n(x,\cdot) - \pi(\cdot)\|_{TV} \to 0$ as $n \to \infty$, for $\pi$-a.e.  $x \in {\mathbf{X}}$. Recall that an invariant distribution $\pi(\cdot)$ will be the unique limiting measure if $P$ is both $\pi$-irreducible and aperiodic \[e.g. [ @meyn2012markov; @tierney1994markov]\].  We note that the convergence results presented here could equivalently be shown under the $V$-norm distance \[see [@roberts1997geometric] for details\]. 

Overview of Main Results {#sec:main}
========================

The majority of results in this paper concern the version of HMC which is typically used in practice, in which the ‘integration time’ for a typical proposal is chosen independently of the current value of the Hamiltonian.  In this scenario we have the following result. 

If Assumptions **A1** (on page ), **A2** (on page ) and **A3** (on page ) hold, then a Markov chain produced by the Hamiltonian Monte Carlo method (outlined in Algorithm \[alg:hmc\]) will be geometrically ergodic. 

Assumption **A1** introduces a controlled degree of randomness into the integration time parameter, which ensures ergodicity of the HMC transition kernel.  Assumption **A2** imposes conditions on the distribution from which expectations are desired, essentially restricting the tail behaviour to be at least as light as that of a Laplace distribution, but no lighter than a Gaussian.  This is to ensure that when the chain is very far from the ‘centre’ of the space then typical proposals will bring it back to regions where probability mass concentrates.  Assumption **A3** relates to the Metropolis–Hastings acceptance rate, ensuring that this does not behave undesirably, in the sense that desirable proposals are often rejected.  We make these arguments precise in Section \[sec:fixed\]. 

We also present the following conditions under which Markov chains produced using HMC will not be geometrically ergodic. 

If either of the following hold then HMC will *not* produce a geometrically ergodic Markov chain:

\(i) $\lim_{\|x\| \to \infty}\frac{\|\nabla U(x)\|}{\|x\|} = \infty$ and (\[eqn:oscillations\]) and (\[eqn:oscillation2\]) are satisfied

\(ii) There is an $M<\infty$ such that $\|\nabla U(x)\| \leq M$ for all $x \in {\mathbf{X}}$, and $\mathbb{E}_\pi[e^{s\|x\|}]=\infty$ for every $s>0$.

The first of these scenarios in essence covers the case where the distribution of interest has lighter tails than those of a Gaussian distribution.  In this case explicit numerical solvers for Hamilton’s equations typically become unstable in some regions of the state space.  The second is concerned with ‘heavy tailed’ distributions, in which the resulting Hamiltonian flow can be slow, precluding a geometric rate of convergence. 

To give some intuition for these results, we apply them to the Exponential Family class of models first introduced in [@roberts1996exponential], in which $\pi(x) \in C^1({\mathbf{X}})$ and for all $\|x\|>M$ for some $M<\infty$ it holds that $$\pi(x) \propto \exp \left( -\alpha\|x\|^{\beta} \right)$$ for some $\alpha,\beta > 0$ and any $x \in {\mathbb{R}}^d$. Different choices of $\beta$ correspond to different tail behaviours, with larger values resulting in ‘lighter’ tails.  For $\beta \geq 1$ the density is log-concave, and the specific choices $\beta = 1$ and $\beta = 2$ correspond to Laplace and Gaussian distributions. 

\[crly:expfam\] For the exponential family class of models, the following results hold:

\(i) For $1 \leq \beta \leq 2$, the Hamiltonian Monte Carlo method will produce a geometrically ergodic chain (with careful choice of $L$ and ${\varepsilon}$ in the $\beta = 2$ case)

\(ii) If $\beta < 1$ or $\beta > 2$, then the Hamiltonian Monte Carlo method will *not* produce a geometrically ergodic chain

The results are analogous to those found for the Metropolis-adjusted Langevin algorithm in [@roberts1996exponential].  A key finding of this work is that when the integration time parameter is chosen in a manner which is independent of the current position, then the two methods essentially coincide in terms of presence or absence of geometric ergodicity.  In other words, taking more than a single leapfrog step in the method can never result in a chain ‘becoming’ geometrically ergodic, even though it may still improve the speed of convergence. 

We also consider an idealised version of the method in Section \[sec:dynamic\], in which the integration time is allowed to depend on the current position in a prescribed way.  This scheme was designed to mimic several more recent versions of HMC \[e.g. [ @hoffman-gelman:2013]\] which are commonly used in modern software packages \[e.g. [ @carpenter2016stan]\].  For a specific one-dimensional class of smooth exponential family models we find the following

For the one-dimensional class of distributions with densities of the form $$\pi(x) \propto \exp \left( -\beta^{-1} (1+x^2)^{\beta/2} \right),$$ then the idealised Hamiltonian Monte Carlo method introduced in Section \[sec:dynamic\] will produce a geometrically ergodic Markov chain for any choice of $\beta>0$.

The positive result in the case where $\beta>2$ is an artifact of the assumption that Hamilton’s equations can be exactly solved in the idealised scheme - this result would disappear if a typical explicit numerical solver were used instead.  However, the findings for the case $\beta<1$ suggest that there are advantages to using an position-dependent integration time in the presence of heavy tails.  We discuss this in more detail in Section \[sec:discussion\]. 

Preliminaries {#sec:prelim}
=============

The approach taken here to establishing geometric convergence was popularised in the monograph [@meyn2012markov], who built on earlier concepts from Doeblin \[e.g. [ @lindvall2002lectures]\] and Harris \[e.g. [ @harris1956existence; @baxendale2011te]\] in particular \[see also [@hairer2011yet]\].  A key observation shown in that work is that for $\pi$-irreducible and aperiodic Markov chains establishing a bound of the form $$\label{eqn:geoerg}
\|P^n(x,\cdot) - \pi(\cdot)\|_{TV} \leq M(x)\rho^n$$ for some $\rho<1$ and $\pi$-a.e.  finite $M:{\mathbf{X}}\to[0,\infty]$ is equivalent to finding a $\pi$-a.e.  finite *Lyapunov* function $V:{\mathbf{X}}\to [1,\infty]$ with ‘small’ level sets \[e.g. [ @baxendale2011te]\] such that the condition $PV(x) \leq \lambda V(x) + b\mathbbm{1}_{C_d}(x)$ holds for some $\lambda < 1$, $b<\infty$ and some set $C_d:=\{x : V(x)\leq d\}$ with $d<\infty$.

It is shown in [@roberts1997geometric] that if $P$ is reversible then (\[eqn:geoerg\]) is equivalent to the Markov chain exhibiting an $L_2(\pi)$ absolute spectral gap.  In this case then a Central Limit Theorem can be deduced for long-run averages of any $f \in L_2(\pi)$ \[e.g. [ @roberts2004general]\]. 

We are concerned here with specific forms of $P$.

We say $P$ is of the *Metropolis–Hastings* type if $$\label{eqn:mhtype}
P(x,dy) = \alpha(x,y)Q(x,dy) + r(x)\delta_x(dy),$$ where $Q$ is a Markov kernel, $\alpha(x,y)$ is defined in (\[eqn:accept\]) and $r(x) = 1- \int \alpha(x,y)Q(x,dy)$.

[@roberts1996geometric] showed the following when $P$ is of the form (\[eqn:mhtype\]). 

\[prop:mh1\] If $\pi(\cdot)$ and $Q(x,\cdot)$ admit Lebesgue densities $\pi(x)$ and $q(y|x)$, $\pi(x)$ is bounded away from $0$ and $\infty$ on compact sets, and there exists $\delta_q>0$ and $\epsilon_q>0$ such that, for every $x$, $$\|x-y\| \leq \delta_q \implies q(y|x) \geq \epsilon_q,$$ then the Metropolis–Hastings chain with candidate density $q(y|x)$ is $\pi$-irreducible and aperiodic, and all compact sets are small. 

If $P$ is of Metropolis–Hastings type and the conditions of Proposition \[prop:mh1\] are satisfied, then (\[eqn:geoerg\]) holds if and only if $$\label{eqn:drift}
\limsup_{\|x\|\to\infty}\frac{PV(x)}{V(x)} < 1,$$ for some Lyapunov function $V$.

Showing a lack of geometric ergodicity typically requires careful study of the distribution of return times to small sets.  The following result of [@roberts1996geometric], however, provides a straightforward method for doing this for Metropolis–Hastings kernels. 

\[prop:supr\] If $P$ is of Metropolis–Hastings type, then (\[eqn:geoerg\]) *fails* to hold if $\operatorname*{ess\,sup}r(x) = 1$.

Lack of geometric ergodicity can also be established in some cases using the following result of [@jarner2003necessary]. 

\[prop:tight\] If for any $\eta > 0$ there is a $\delta > 0$ such that $$\label{eqn:lackge2}
P(x,B_\delta(x)) > 1-\eta,$$ where $B_\delta(x) :=\{ y \in {\mathbf{X}}: \|x-y\| < \delta \}$, then $P$ can be geometrically ergodic *only* if $\mathbb{E}_\pi [e^{\beta \|x\|}] < \infty$ for some $\beta > 0$.

If $P$ is of Metropolis–Hastings type, it is straightforward to verify that $Q(x,B_\delta(x)) > 1- \eta$ ensures (\[eqn:lackge2\]), meaning we only need consider the candidate kernel in these cases. 

Hamiltonian Monte Carlo {#sec:hmc}
=======================

We give a brief introduction here.  For a more detailed account see [@neal2011mcmc] or [@betancourt2014geometric].  We consider probability densities of the form $\pi(x) \propto e^{-U(x)}$ for some $U:{\mathbf{X}}\to [0,\infty)$. If we view $U(x)= -\log\pi(x)$ as a ‘potential’ energy in a physical system, it is natural to consider the larger phase space and construct the *Hamiltonian* $$\label{eqn:hamiltonian}
H(x,p) = U(x) + \frac{1}{2}p^t M^{-1} p,$$ where $p$ denotes a $d$-dimensional ‘momentum’ variable, $M$ a $d \times d$ ‘mass’ matrix and $p^t M^{-1}p/2$ the ‘kinetic’ energy (other forms of kinetic energy are also possible, see e.g. [@girolami2011riemann]).  Provided $U(x)$ is differentiable, we can evolve the coordinates $(x_t,p_t)$ through time in such a way that $H(x_t,p_t) = H(x_{t+s},p_{t+s})$ for any $t,s \in \mathbb{R}$ using Hamilton’s equations $$\label{eqn:hamilton}
\frac{d p_t}{dt} = -\frac{\partial H}{\partial x}, ~~ \frac{dx_t}{dt} = \frac{\partial H}{\partial p}.$$ Solving (\[eqn:hamilton\]) results in *Hamiltonian flow*.  To put this presentation into the framework introduced in Section \[sec:intro\], we can consider constructing a measure-preserving map $f_\theta:{\mathbf{X}}\to {\mathbf{X}}$ by setting the input to be $x_0$, choosing a momentum variable $p_0$, solving (\[eqn:hamilton\]) for $T$ units of time and then projecting back down onto ${\mathbf{X}}$ to produce $x_T$. The parameters $\theta = \{p_0,T\}$ define the behaviour of a single map $f_\theta$, and how they are chosen define the behaviour of the Markov chain produced by iterating the process of randomly selecting a $\theta$ and then applying the resulting map $f_\theta$ to the current point to produce the next. 

Of course, it is often not possible to solve (\[eqn:hamilton\]) exactly, so numerical methods are needed.  Fortunately, the rich geometric structure of Hamiltonian systems allows the construction of *symplectic* integrators, which possess attractive long term numerical stability properties \[e.g. [ @leimkuhler2004simulating; @sanz1994numerical; @hairer2006geometric]\], meaning that for appropriate Hamiltonians the approximate solution of (\[eqn:hamilton\]) is such that $H(x_t,p_t) \approx H(x_0,p_0)$ for all $t<\eta$, where $\eta \gg 0$. The standard choice when the Hamiltonian is of the form (\[eqn:hamiltonian\]) is the Störmer–Verlet or *leapfrog* scheme, in which $(x_{L{\varepsilon}},p_{L{\varepsilon}})$ is generated from $(x_0,p_0)$ using $L$ steps of the recursion $$\begin{aligned}
p_{t+\frac{{\varepsilon}}{2}} &= p_t - \frac{{\varepsilon}}{2}\nabla U(x_t), \\
x_{t+{\varepsilon}} &= x_t + {\varepsilon}M^{-1}p_{t+\frac{{\varepsilon}}{2}}, \\
p_{t+{\varepsilon}} &= p_{t+\frac{{\varepsilon}}{2}} - \frac{{\varepsilon}}{2}\nabla U(x_{t+{\varepsilon}}),\end{aligned}$$ for some step-size ${\varepsilon}>0$. Although the resulting *approximate* flow map $\varphi_{L{\varepsilon}}(x_0,p_0) := (x_{L{\varepsilon}},p_{L{\varepsilon}})$ no longer preserves $\pi(\cdot)$, it can be used as a proposal mechanism within the Metropolis–Hastings framework \[e.g. [ @neal2011mcmc]\].  The full method is shown in Algorithm \[alg:hmc\] below. 

\[alg:hmc\] $x_{i-1}$, $\epsilon \geq 0$, $\mathfrak{L}(\cdot)$ Set $x_0 \leftarrow x_{i-1}$, draw $p_0 \sim N(0,M)$, $L \sim \mathfrak{L}(\cdot)$, set $T \leftarrow L\varepsilon$ Draw $u \sim U[0,1]$ Set $\delta \leftarrow H(x_0,p_0) - H\circ \varphi_T(x_0,p_0)$, Set $x_i \leftarrow \text{Pr}_x \circ \varphi_T(x_0,p_0)$ Set $x_i \leftarrow x_{i-1}$

It is known that stability conditions on ${\varepsilon}$ are required in order for the leapfrog scheme to be a proficient numerical solver \[e.g. [ @leimkuhler2004simulating]\].  We show in Section \[sec:fixed\] that such conditions, as expected, are also needed to ensure that Algorithm \[alg:hmc\] produces a geometrically ergodic chain. 

From this point forward we assume $M = I$ for ease of exposition but without loss of generality. 

The marginal chain
------------------

To use the techniques of [@meyn2012markov], it is helpful to express the HMC transition in such a way that when $\|x\|$ is large it is clear how the chain will behave.  Although it is typically presented as a map on the larger phase space, HMC can simply be thought of as a Markov chain on ${\mathbf{X}}$, and we will find this representation useful in relation to the above.  In this case the candidate map $g_\xi$ is given by the following proposition. 

\[prop:mprop\] The HMC candidate map can be written $$\label{eqn:mprop}
x_{L{\varepsilon}} = x_0 - \frac{L{\varepsilon}^2}{2}\nabla U(x_0) - {\varepsilon}^2\sum_{i=1}^{L-1} (L-i)\nabla U(x_{i{\varepsilon}}) + L{\varepsilon}p_0.$$ where $p_0 \sim N(0,I)$, $L$ is the number of leapfrog steps and ${\varepsilon}$ the integrator step-size.  With this choice, the acceptance probability will be $$\alpha(x_0,x_{L{\varepsilon}}) = 1 \wedge \frac{\pi(x_{L{\varepsilon}})}{\pi(x_0)}\exp \left( \frac{1}{2}\|p_0\|^2 - \frac{1}{2}\|p_{L{\varepsilon}}\|^2 \right),$$ where $$\label{eqn:mprop2}
p_{L{\varepsilon}} = p_0 - \frac{{\varepsilon}}{2}\nabla U(x_0) - {\varepsilon}\sum_{i=1}^{L-1} \nabla U(x_{i{\varepsilon}}) - \frac{{\varepsilon}}{2}\nabla U(x_{L{\varepsilon}}).$$

Proposition \[prop:mprop\] highlights the previously noted relationship between HMC and MALA quite explicitly, as setting $L = 1$ means the third term on the right-hand side of (\[eqn:mprop\]) disappears, leaving the MALA proposal $x_0 - {\varepsilon}^2\nabla U(x_0)/2 + {\varepsilon}p_0$. It also highlights why taking $L >1$ proposes a greater challenge, as for each $x_{i{\varepsilon}}$ with $i \geq 1$ this term will typically be a nonlinear transformation of $x_0$ and $p_0$. As $p_0$ is stochastic, then ${\varepsilon}^2\sum_{i=1}^{L-1} (L-i)\nabla U(x_{i{\varepsilon}})$ will be also, but its distribution will often be intractable. 

Results for an position-independent integration time {#sec:fixed}
====================================================

In this section we make the assumption that the distribution $\mathfrak{L}(\cdot)$ for the number of leapfrog steps $L$ does not depend on the current position.  This is relaxed in Section \[sec:dynamic\]. 

$\pi$-irreducibility
--------------------

We first present a simple (and well-known) example that shows how the proposal transition given by (\[eqn:mprop\]) can produce a method which is not $\pi$-irreducible, and hence will not be ergodic. 

\[ex:irr\] Take $\pi(x) \propto e^{-x^2/2}$, meaning $\nabla U(x) = x$, and set $L = 2$. Then the HMC proposal becomes $$\begin{aligned}
x_{2{\varepsilon}} 
&= x_0 - {\varepsilon}^2 x_0 - {\varepsilon}^2 (x_0 - {\varepsilon}^2 x_0 + {\varepsilon}p_0) + 2{\varepsilon}p_0, \\
&= x_0 - {\varepsilon}^2 x_0 - {\varepsilon}^2x_0 + {\varepsilon}^4 x_0 - {\varepsilon}^3 p_0 + 2{\varepsilon}p_0, \\
&= (1 - 2{\varepsilon}^2 + {\varepsilon}^4)x_0 + (2{\varepsilon}- {\varepsilon}^3)p_0.\end{aligned}$$ Setting ${\varepsilon}= \sqrt{2}$ means $2{\varepsilon}- {\varepsilon}^3 = 0$, so that $$x_{2{\varepsilon}} = (1 - 4 + 4)x_0 = x_0.$$ With this transition the proposal kernel is simply $Q(x,\cdot) = \delta_x(\cdot)$, so the chain is not $\pi$-irreducible unless $\pi(\cdot) = \delta_x(\cdot)$.

Although it is in some sense trivial, Example \[ex:irr\] highlights that establishing $\pi$-irreducibility is not so straightforward here as for Metropolis–Hastings methods based on random walks or overdamped Langevin diffusions.  In essence the problem is that the third and last terms on the right-hand side of (\[eqn:mprop\]) are *both* random, and if the model is of a certain type and the parameters chosen in a particular manner then these can interact in a very undesirable way.  The above case is the discrete time analogue to integrating the harmonic oscillator over precisely one period \[e.g. [ @leimkuhler2004simulating]\]. 

The observation noted here and elsewhere that HMC in the case $L=1$ corresponds to MALA, for which irreducibility is established in [@roberts1996exponential], can be exploited to alleviate these issues and establish $\pi$-irreducibility of HMC under the following assumption. 

**A1** [*The distribution $\mathfrak{L}(\cdot)$ is such that $\mathbb{P}_{\mathfrak{L}}[L=1]>0$, and for any fixed $(x_0,p_0)\in {\mathbb{R}}^{2d}$ and ${\varepsilon}>0$, and that there is an $s <\infty$ such that $\mathbb{E}_\mathfrak{L} [e^{s\|x_{L{\varepsilon}}\|}]<\infty$.*]{} \[ass:A1\]

When assumption **A1** holds then the fact that HMC produces an ergodic Markov chain can be straightforwardly invoked from existing MALA results \[[@roberts1996exponential]\].  The idea of randomising the integration time is commonly recommended for practical applications of the method \[e.g. [ @neal2011mcmc; @girolami2011riemann]\], and more theoretical motivation for doing so is given in [@betancourt2016identifying].  The finite exponential expectation condition is needed to ensure that the Lyapunov function used to prove geometric ergodicity is valid.  One simple way to ensure this in practice (under the additional assumptions imposed on the potential $U$ in the next subsection) is that $\mathbb{P}_{\mathfrak{L}}[L>l]=0$ for some fixed $l<\infty$, though weaker conditions than this are also possible. 

Assumption **A1** can be viewed as the discrete time analogue to the the exponential integration time assumption made in [@bou2015randomized], and in many respects serves a similar purpose.  Similar conditions are also exploited to prove $\pi-$irreducibility results in [@cances2007theoretical]. 

Geometric ergodicity
--------------------

We first present here some seemingly abstract conditions under which the HMC method produces a geometrically ergodic Markov chain.  We then give some natural assumptions on the potential $U(x)$ under which these hold. 

We present the results of this section conditioned on a fixed choice of the number of leapfrog steps $L$, for ease of exposition.  Note that the required drift conditions shown hold for a fixed $L$, then under **A1** they will hold when possible values for $L$ are averaged over according to $\mathfrak{L}(\cdot)$, so this does not affect the generality of the results. 

*Notation.*  We introduce some further notation for this section.  Let $I_\delta(x) := \{ y \in {\mathbf{X}}: \|y\| \leq \|x\|^\delta \}$ for some $1/2 < \delta < 1$. In the case $\delta=1$ we will simply write $I(x)$. Let $$m_{L,{\varepsilon}}(x_0,p_0) := x_0 - L{\varepsilon}^2\nabla U(x_0)/2 - {\varepsilon}^2\sum_{i=1}^{L-1} (L-i)\nabla U(x_{i{\varepsilon}})$$ denote the ‘mean’ next candidate position ($x_{L{\varepsilon}}-L{\varepsilon}p_0$), and $$\psi_{L,{\varepsilon}}(x_0,p_0):= L{\varepsilon}^2\nabla U(x_0)/2 + {\varepsilon}^2\sum_{i=1}^{L-1} (L-i)\nabla U(x_{i{\varepsilon}})$$ denote the proposal ‘drift’ (implying $m_{L,{\varepsilon}}(x_0,p_0) = x_0 - \psi_{L,{\varepsilon}}(x_0,p_0)$).  We will also sometimes write $h:={\varepsilon}^{2}/2$ in a most likely futile attempt to keep things readable. 

\[thm:geoerg\] The HMC method produces a geometrically ergodic Markov chain if assumption **A1** holds, and in addition both $$\label{eqn:mean}
\limsup_{\|x_0\|\to\infty, \|p_0\| \leq \|x_0\|^\delta} \left( \|m_{L,{\varepsilon}}(x_0,p_0)\| - \|x_0\| \right) < 0$$ and $$\label{eqn:inwards}
\lim_{\|x_0\|\to\infty} \int_{R(x_0) \cap I(x_0)} Q(x_0,dy) = 0,$$ where $R(x_0) := \{ y \in {\mathbf{X}}: \alpha(x_0,y) < 1 \}$ denotes the ‘potential rejection region’ and $I(x_0) := \{ y \in {\mathbf{X}}: \|y\| \leq \|x_0\| \}$ the ‘interior’ of $x_0$.

Take $V(x) = e^{s\|x\|}$ for some $s > 0$. It is shown in the proof of Theorem 4.1 in [@roberts1996exponential] that under (\[eqn:mean\]) and (\[eqn:inwards\]) we can write $$\frac{PV(x_0)}{V(x_0)} \leq \int e^{s(\|y\| - \|x_0\|)}Q(x_0,dy) + \int_{R(x_0) \cap I(x_0)} Q(x_0,dy).$$ The last integral asymptotes to zero as $\|x_0\| \to \infty$ by (\[eqn:inwards\]).  Writing $x_{L{\varepsilon}}(p_0)$ to indicate that $x_{L{\varepsilon}}$ depends on $p_0$, the first integral can be written $$\label{eqn:brokenupintegral}
\int_{I_\delta(x_0)} e^{s(\|x_{L{\varepsilon}}(p_0)\| - \|x_0\|)} \mu^G(dp_0) + \int_{I_\delta(x_0)^c} e^{s(\|x_{L{\varepsilon}}(p_0)\| - \|x_0\|)} \mu^G(dp_0).$$ Noting that $\|x_{L{\varepsilon}}(p_0)\| \leq \|m_{L,{\varepsilon}}(x_0,p_0)\| + L{\varepsilon}\|p_0\|$ for large enough $\|x_0\|$ and using (\[eqn:mean\]) above gives $$\int_{I_\delta(x_0)} e^{s(\|x_{L{\varepsilon}}\| - \|x_0\|)} \mu^G(dp_0) \leq e^{s(\|m_{L,{\varepsilon}}(x_0,p_0)\| - \|x_0\|)} \int_{I_\delta(x_0)^c} e^{sL{\varepsilon}\|p_0\|}\mu^G(dp_0).$$ The last integral is the absolute moment generating function for a truncated Gaussian random variable, so can be upper bounded by $2e^{(sL{\varepsilon})^2/2}\Phi(sL{\varepsilon})$, where $\Phi$ denotes the standard Gaussian cumulative distribution function.  As this is positive and finite for any $x_0$, and its size is controlled by a factor of $s^2$, then by (\[eqn:mean\]) the integral asymptotes to a quantity which is strictly less than one if $s>0$ is chosen to be suitably small. 

It remains to show that the right-hand integral in (\[eqn:brokenupintegral\]) becomes negligibly small as $\|x_0\| \to \infty$. It follows from (\[eqn:mprop\]) and (A4.4.1) that $\|x_{L{\varepsilon}}\| \in O(\max(\|x_0\|,\|p_0\|))$. This means that for some constants $C \in {\mathbb{R}}$ and for $p_0 \in I_\delta(x_0)^c$ and $\|x_0\|$ large enough we can write $$\begin{aligned}
\exp\left( s\|x_{L{\varepsilon}}\| - s\|x_0\| - \frac{1}{2}\|p_0\|^2  \right) &\leq \exp\left( c\max(\|x_0\|,\|p_0\|) - \frac{1}{2}\|p_0\|^2 \right), \\
&= \exp\left( \|p_0\|\left( C\frac{\max(\|x_0\|,\|p_0\|)}{\|p_0\|} - \|p_0\|\right) \right), \\
&\leq \exp\left( \|p_0\|\left( C\|x_0\|^{1 - \delta} - \|p_0\|\right) \right)\end{aligned}$$ Provided $\delta > 1/2$, then for $\|x_0\|$ large enough $C\|x_0\|^{1-\delta} - \|x_0\|^\delta < -1$, meaning $$\exp\left( \|p_0\|\left( C\|x_0\|^{1 - \delta} - \|p_0\|\right) \right) \leq \exp\left(-\|p_0\|\right),$$ meaning $$\int_{I_\delta(x_0)^c} e^{s(\|x_{L{\varepsilon}}\| - \|x_0\|)}\mu^G(dp_0) \leq \int_{I_\delta(x_0)^c} e^{-\|p_0\|} \mu^G(dp_0) \leq 2e^{-\|x_0\|^\delta},$$ which becomes negligibly small as $\|x_0\| \to \infty$, as required. 

Theorem \[thm:geoerg\] is a generalisation of Theorem 4.1 in [@roberts1996exponential] to the HMC case.  The nontriviality involved in this extension is essentially dealing with the randomness induced into $m_{L,{\varepsilon}}(x_0,p_0)$ from $p_0$.

### Requirements for (\[eqn:mean\]) to be satisfied. 

In the case $L=1$ (\[eqn:mean\]) corresponds to $$\label{eqn:mean1step}
\|x_{0}-h\nabla U(x_{0})\|-\|x_{0}\|<0$$ whenever $\|x_{0}\|>M$ for some $M<\infty$. The statements in this section give three simple conditions which establish this are *also* sufficient to establish that for *typical* values of $\|p_{0}\|$ then $\|m_{L,{\varepsilon}}(x_{0},p_{0})\|$ will also be strictly less than $\|x_{0}\|$ when $L\geq2$. The main result is stated below.  The crucial consequence of this is that controlling the behaviour of ‘global move’ updates produced by HMC when $L>1$ can be done through only ‘local’ knowledge, meaning analytic information at the current point $x_{0}$.

For any $L\geq1$ (\[eqn:mean\]) holds if the following conditions are met

(SC1.1) $\liminf_{\|x_{0}\|\to\infty}\|\nabla U(x_{0})\|>0$

(SC1.2) $\liminf_{\|x_{0}\|\to\infty}\frac{\langle\nabla U(x_{0}),x_{0}\rangle}{\|\nabla U(x_{0})\|\|x_{0}\|}>0$

(SC1.3) $\lim_{\|x_{0}\|\to\infty}\frac{\|\nabla U(x_{0})\|}{\|x_{0}\|}=0$.

In addition, if (SC1.3) is replaced by

(SC1.3b) $\limsup_{\|x_{0}\|\to\infty}\frac{\|\nabla U(x_{0})\|}{\|x_{0}\|}=S_{l},$

for some $S_{l}<\infty$, then there is an ${\varepsilon}_{0}\in(0,\infty)$ such that the same result holds provided ${\varepsilon}\in(0,{\varepsilon}_{0})$.

This is re-stated as Proposition \[prop:SC1.3\] and Proposition \[prop:SC1.3b\] below, which follow from the preceding Lemmas. 

The conditions of the result are intuitive.  Condition (SC1.1) ensures that the gradient is non-negligible for large values of $\|x_{0}\|$, while (SC1.3) ensures that the asymptotic growth in $\|\nabla U(x_{0})\|$ is sublinear.  Condition (SC1.3) ensures that $\nabla U(x_{0})$ asymptotically ‘points inwards’.  We begin with a straightforward observation. 

Sufficient conditions such that $$\limsup_{\|x_{0}\|\to\infty}\left(\|x_{0}-h\nabla U(x_{0})\|-\|x_{0}\|\right)<0$$ are:

(SC1.1) $\liminf_{\|x_{0}\|\to\infty}\|\nabla U(x_{0})\|>0$

(SC1.2a) $\limsup_{\|x_{0}\|\to\infty}\left(\frac{h^{2}\|\nabla U(x_{0})\|^{2}-2h\langle\nabla U(x_{0}),x_{0}\rangle}{2\|x_{0}\|}\right)<0$

(SC1.3) $\lim_{\|x_{0}\|\to\infty}\frac{\|\nabla U(x_{0})\|}{\|x_{0}\|}=0$.

First note that $\|x_{0}-h\nabla U(x_{0})\|=\sqrt{\|x_{0}\|^{2}+h^{2}\|\nabla U(x_{0})\|^{2}-2h\langle\nabla U(x_{0}),x_{0}\rangle}$. Recall the generalised Bernoulli inequality: if $y>-1$ and $r\in[0,1]$ then $\left(1+y\right)^{r}\leq1+ry$. Setting $r:=1/2$, $a(x_{0}):=\|x_{0}\|^{2}$ and $b(x_{0}):=h^{2}\|\nabla U(x_{0})\|^{2}-2h\langle\nabla U(x_{0}),x_{0}\rangle$ then we have $$\begin{aligned}
a(x_{0})^{r}\left(1+\frac{b(x_{0})}{a(x_{0})}\right)^{r} &\leq a(x_{0})^{r}\left(1+\frac{b(x_{0})}{2a(x_{0})}\right) \\ 
&=\|x_{0}\|+\frac{h^{2}\|\nabla U(x_{0})\|^{2}-2h\langle\nabla U(x_{0}),x_{0}\rangle}{2\|x_{0}\|}.\end{aligned}$$ This will be strictly less than $\|x_{0}\|$ in the limit as $\|x_{0}\|\to\infty$ provided $b(x_{0})/a(x_{0})>-1$. Noting that $$\frac{b(x_{0})}{a(x_{0})}=\frac{h^{2}\|\nabla U(x_{0})\|^{2}-2h\langle\nabla U(x_{0}),x_{0}\rangle}{\|x_{0}\|^{2}}>\frac{h^{2}\|\nabla U(x_{0})\|^{2}-2h\|\nabla U(x_{0})\|\|x_{0}\|}{\|x_{0}\|^{2}},$$ then it suffices to see that under (SC1.3) the right-hand side can be made arbitrarily close to zero by taking $\|x_{0}\|$ large enough. 

We can also recover some more intuitive sufficient conditions. 

\[crly:intuitive\]A more intuitive condition which implies (SC1.2a) conditional on (SC1.1) and (SC1.3) is

(SC1.2) $\liminf_{\|x_{0}\|\to\infty}\frac{\langle\nabla U(x_{0}),x_{0}\rangle}{\|\nabla U(x_{0})\|\|x_{0}\|}>0$

Using condition (SC1.3) gives $$\begin{aligned}
\liminf_{\|x_{0}\|\to\infty}\frac{\langle\nabla U(x_{0}),x_{0}\rangle}{\|\nabla U(x_{0})\|\|x_{0}\|} & >0\implies\limsup_{\|x_{0}\|\to\infty}\frac{h^{2}\|\nabla U(x_{0})\|}{2\|x_{0}\|}-\frac{h\langle\nabla U(x_{0}),x_{0}\rangle}{2\|\nabla U(x_{0})\|\|x_{0}\|}>0\\
 & \implies\limsup_{\|x_{0}\|\to\infty}\left(\frac{h^{2}\|\nabla U(x_{0})\|^{2}-2h\langle\nabla U(x_{0}),x_{0}\rangle}{2\|x_{0}\|}\right)<0.\end{aligned}$$

From now on we refer to (SC1.1), (SC1.2) and (SC1.3) combined as (SC1.1)-(SC1.3).  Next we show that these same broad and natural conditions are sufficient for (\[eqn:mean\]) to hold. 

\[lem:sufficient\_L&gt;1\]Under the following conditions (\[eqn:mean\]) holds:

\(i) $\liminf_{\|x_{0}\|\to\infty,\|p_{0}\|\leq\|x_{0}\|^{\delta}}\left(\frac{\|\psi_{L,{\varepsilon}}\|^{2}-2\langle\psi_{L,{\varepsilon}},x_{0}\rangle}{\|x_{0}\|^{2}}\right)>-1$

\(ii) $\limsup_{\|x_{0}\|\to\infty,\|p_{0}\|\leq\|x_{0}\|^{\delta}}\left(\frac{\|\psi_{L,{\varepsilon}}\|^{2}-2\langle\psi_{L,{\varepsilon}},x_{0}\rangle}{2\|x_{0}\|}\right)<0$.

Using the generalised Bernoulli inequality as above gives the result. 

Next we relate the conditions of Lemma \[lem:sufficient\_L&gt;1\] to *local* criteria, i.e. criteria that only depend on the current point $x_{0}$. The following lemmas give a starting point. 

\[lem:U\_lem\_1\]Provided $\|p_{0}\|\leq\|x_{0}\|^{\delta}$ and (SC1.3) holds then we have the following

\(i) For any $\eta>0$ there is an $M_{\eta<\infty}$ such that whenever $\|x_{0}\|>M_{\eta}$ it holds that $(1-\eta)\|x_{0}\|\leq\|x_{{\varepsilon}}\|\leq(1+\eta)\|x_{0}\|$

\(ii) $\|\nabla U(x_{{\varepsilon}})\|=o(\|x_{0}\|)$

\(iii) $\|p_{{\varepsilon}}\|\in o(\|x_{0}\|)$.

\(i) Noting that $\|x_{{\varepsilon}}\|=\|x_{0}-h\nabla U(x_{0})+{\varepsilon}p_{0}\|$ gives $$\begin{aligned}
(1-\delta)\|x_{0}\| &\leq\|x_{0}\|-h\|\nabla U(x_{0})\|-{\varepsilon}\|x_{0}\|^{\delta} \\
&\leq\|x_{0}-h\nabla U(x_{0})+{\varepsilon}p_{0}\| \\
&\leq\|x_{0}\|+h\|\nabla U(x_{0})\|+{\varepsilon}\|x_{0}\|^{\delta} \\
&\leq(1+\delta)\|x_{0}\|.\end{aligned}$$

\(ii) We have from (i) and (SC1.3) that for any $\gamma>0$ there is an $M_{\gamma}<\infty$ such that whenever $\|x_{0}\|>M_{\gamma}/(1-\delta$) then $\|\nabla U(x_{{\varepsilon}})\|/\|x_{{\varepsilon}}\|<\gamma$. This implies using (i) that $\|\nabla U(x_{{\varepsilon}})\|/\|x_{0}\|<\gamma(1-\delta)$, and since $\gamma(1-\delta)$ can be made arbitrarily small then the result follows. 

\(iii) $\|p_{{\varepsilon}}\|=\|p_{0}-{\varepsilon}\nabla U(x_{0})/2-{\varepsilon}\nabla U(x_{{\varepsilon}})/2\|\leq\|p_{0}\|+{\varepsilon}\|\nabla U(x_{0})\|/2+{\varepsilon}\|\nabla U(x_{{\varepsilon}})\|/2$, which is $\in o(\|x_{0}\|)$ using (i) and (ii) and the fact that $\|p_{0}\|\leq\|x_{0}\|^{\delta}$

\[lem:orders\_L\]Provided $\|p_{0}\|\leq\|x_{0}\|^{\delta}$ and (SC1.3) holds then for any $L<\infty$ and each $i\in\{0,...,L-1\}$ the following hold

\(i) For any $\eta>0$ there is an $M_{\eta<\infty}$ such that whenever $\|x_{0}\|>M_{\eta}$ it holds that $(1-\eta)\|x_{0}\|\leq\|x_{i{\varepsilon}}\|\leq(1+\eta)\|x_{0}\|$

\(ii) $\|\nabla U(x_{i{\varepsilon}})\|\in o(\|x_{0}\|)$

\(iii) $\|p_{i{\varepsilon}}\|\in o(\|x_{0}\|)$

\(iv) $\|\psi_{L,{\varepsilon}}\|\in o(\|x_{0}\|)$.

The results follow iteratively for each $i$ using the same approach as in the previous Lemma.  For the case $i=2$ then noting that $\|x_{2{\varepsilon}}\|=\|x_{{\varepsilon}}-h\nabla U(x_{{\varepsilon}})+{\varepsilon}p_{{\varepsilon}}\|$, then (i) in this case follows from Lemma \[lem:U\_lem\_1\].  It follows that $\|\nabla U(x_{2{\varepsilon}})\|\in o(\|x_{0}\|)$ and $\|p_{2{\varepsilon}}\|\in o(\|x_{0}\|)$ by an analogous argument to this Lemma.  Given this then it can be shown that (i) holds for $i=3$, and then (ii) and (iii) by the same logic, and the argument can be iterated as many times as is needed.  The last claim follows trivially from the second. 

\[prop:SC1.3\]Under (SC1.1)-(SC1.3) then for any $L<\infty$ the conditions of Lemma \[lem:sufficient\_L&gt;1\] are satisfied. 

First we show (i).  Writing $x^{*}:=\arg\max_{i\in{0,...,L-1}}\left\{ \|\nabla U(x_{i{\varepsilon}})\|\right\} $, then we have $\|\psi_{L,{\varepsilon}}\|\leq L{\varepsilon}^{2}\sum_{i=0}^{L-1}\|\nabla U(x_{i{\varepsilon}})\|\leq L^{2}{\varepsilon}^{2}\|\nabla U(x^{*})\|$, which implies $$\frac{\|\psi_{L,{\varepsilon}}\|^{2}}{\|x_{0}\|^{2}}\leq\frac{L^{4}{\varepsilon}^{4}\|\nabla U(x^{*})\|^{2}}{(1-\eta)^{2}\|x^{*}\|^{2}}$$ which can be made arbitrarily small by taking $\|x_{0}\|$ large enough using (SC1.3).  Noting that $\|\psi_{L,{\varepsilon}}\|/\|x_{0}\|\geq\langle\psi_{L,{\varepsilon}},x_{0}\rangle/\|x_{0}\|^{2}\geq-\|\psi_{L}\|/\|x_{0}\|$, then an analogous argument can be used to show that $-2\langle\psi_{L},x_{0}\rangle/\|x_{0}\|^{2}$ will also tend to zero as $\|x_{0}\|\to\infty$.

\(ii) First note from above that $\lim_{\|x_{0}\|\to\infty}\|\psi_{L,{\varepsilon}}\|^{2}/\left(L^{2}{\varepsilon}^{2}\|\nabla U(x^{*})\|\|x_{0}\|\right)=0$. By an analogous argument to that used in the proof of Corollary \[crly:intuitive\], it is clear therefore that (ii) holds if $$\limsup_{\|x_{0}\|\to\infty}\left(\frac{\|\psi_{L,{\varepsilon}}\|^{2}}{L^{2}{\varepsilon}^{2}\|\nabla U(x^{*})\|\|x_{0}\|}-\frac{\langle\psi_{L,{\varepsilon}},x_{0}\rangle}{L^{2}{\varepsilon}^{2}\|\nabla U(x^{*})\|\|x_{0}\|}\right)<0$$ which in turn is equivalent to the statement $$\liminf_{\|x_{0}\|\to\infty}\left(\frac{\langle\psi_{L,{\varepsilon}},x_{0}\rangle}{L^{2}{\varepsilon}^{2}\|\nabla U(x^{*})\|\|x_{0}\|}\right)>0.$$ The numerator can be decomposed as $$\begin{aligned}
\langle\psi_{L,{\varepsilon}},x_{0}\rangle &\geq \sum_{i=0}^{L-1}c_{i}\langle\nabla U(x_{i{\varepsilon}}),x_{0}\rangle \\
&=\sum_{i=0}^{L-1}c_{i}\langle\nabla U(x_{i{\varepsilon}}),x_{i{\varepsilon}}\rangle+\sum_{i=0}^{L-1}c_{i}\langle\nabla U(x_{i{\varepsilon}}),x_{0}-x_{i{\varepsilon}}\rangle,\end{aligned}$$ where each $c_{i}=(L-i){\varepsilon}^{2}$ for $i\geq1$ and $c_{0}=L{\varepsilon}^{2}/2$ .  The second of these terms is $o(\|x_{0}\|)$ using the Cauchy–Schwartz inequality and Lemma \[lem:orders\_L\] (which shows that both $\|\psi_{L,{\varepsilon}}\|$ and $\|x_{0}-x_{i{\varepsilon}}\|\in o(\|x_{0}\|)$), and so this term vanishes when divided by $L^{2}{\varepsilon}^{2}\|\nabla U(x^{*})\|\|x_{0}\|$. The first term divided by the same quantity will be strictly positive as each term in the sum is $\geq0$ using (SC1.2) and Lemma \[lem:orders\_L\], and at least one of them is $>0$ since it will correspond to $x^{*}$.

The condition (SC1.3) allows clarity in the proofs, but precludes the natural boundary case of distributions with Gaussian tails.  The following proposition addresses this. 

\[prop:SC1.3b\] If (SC1.1) and (SC1.2) hold, and in addition

(SC1.3b) $\limsup_{\|x_{0}\|\to\infty}\frac{\|\nabla U(x_{0})\|}{\|x_{0}\|}=S_{l}$

for some constant $S_{l}<\infty$, then there is an ${\varepsilon}_{0}\in(0,\infty)$ such that for any choice of ${\varepsilon}\in(0,{\varepsilon}_{0})$ the conditions of Lemma \[lem:sufficient\_L&gt;1\] are satisfied. 

We simply note that the term $\langle\psi_{L,{\varepsilon}},x_{0}\rangle\in O({\varepsilon}^{2})$, while $\|\psi_{L,{\varepsilon}}\|^{2}\in O({\varepsilon}^{4})$, so that the proofs of the preceding Lemmas can be straightforwardly modified when (SC1.3) is replaced by (SC1.3b) by choosing a small enough value of ${\varepsilon}$ that the inner product dominates the square norm.  We omit the details of this. 

The sensitivity to the choice of ${\varepsilon}$ in this case is well known in this scenario as a potential source of numerical instabilities, and choosing ${\varepsilon}< 1/S_l$ is recommended to alleviate such issues \[e.g. [ @leimkuhler2004simulating]\].  We conclude this subsection with the following assumption that we require for a geometrically ergodic Markov chain produced by the HMC method, which is a natural conclusion of the preceding results. 

**A2** [*The potential $U(x)$ satisfies either (SC1.1)-(SC1.3), or it satisfies (SC1.1)-(SC1.3b) and ${\varepsilon}$ is chosen to be suitably small that the conditions of Lemma \[lem:sufficient\_L&gt;1\] are satisfied.*]{}  \[ass:A2\]

It remains to consider (\[eqn:inwards\]), which reflects the role of the acceptance rate in the HMC method.  We turn to this next. 

### Discussion of (\[eqn:inwards\]). 

In [@roberts1996exponential] the authors note that (\[eqn:inwards\]) applied to the MALA transition $x_{{\varepsilon}}=x-{\varepsilon}^{2}\nabla U(x)/2+{\varepsilon}p_{0}$ can be viewed as the restriction that for $x_{{\varepsilon}}\in I(x_{0})$ $$\int_{x_{0}}^{x_{{\varepsilon}}}\nabla U(z)dz-\hat{U}_{1}\leq\frac{{\varepsilon}^2}{8}\left(\|\nabla U(x_{0})\|^{2}-\|\nabla U(x_{{\varepsilon}})\|^{2}\right),$$ where $\hat{U}_{1}:=\frac{1}{2}\left\langle x_{\varepsilon}- x_0,\nabla U(x_{\varepsilon})+\nabla U(x_0)\right\rangle $ denotes the ‘trapezium’ estimate for the line integral $\int_{x}^{y}\nabla U(z)dz$. We can extend this intuition to HMC and arrive at the following natural generalisation of the same condition. 

The acceptance rate for HMC will satisfy the ‘inwards acceptance’ property (\[eqn:inwards\]) if whenever $x_{L{\varepsilon}}\in I(x_{0})$ then in the limit as $\|x_0\|\to\infty$ it holds that $$\label{eqn:ic2}
\int_{x_{0}}^{x_{L{\varepsilon}}}\nabla U(z)dz-\hat{U}_{L}\leq\frac{1}{2L^{2}{\varepsilon}^{2}}\left(\|\psi_{L,{\varepsilon}}\|^{2}-\|\psi_{L,{\varepsilon}}^{R}\|^{2}\right),$$ where $\hat{U}_L:=\frac{1}{2L}\langle x_{L{\varepsilon}}-x_0,\nabla U(x_{0})+\nabla U(x_{L{\varepsilon}})+2\sum_{i=1}^{L-1}\nabla U(x_{i{\varepsilon}})\rangle$ denotes the quadrature rule estimate for the line integral $\int_{x_0}^{x_{L{\varepsilon}}}\nabla U(z)dz$ based on $L$ trapezia and the forward and reverse drift components are given by $$\psi_{L,{\varepsilon}}:=\frac{L{\varepsilon}^{2}}{2}\nabla U(x_{0})+{\varepsilon}^{2}\sum_{i=1}^{L-1}(L-i)\nabla U(x_{i{\varepsilon}}), ~~~~ 
\psi_{L,{\varepsilon}}^{R}:=\frac{L{\varepsilon}^{2}}{2}\nabla U(x_{L{\varepsilon}})+{\varepsilon}^{2}\sum_{i=1}^{L-1}i\nabla U(x_{i{\varepsilon}}).$$

We first note that we can write $p_{0}=\frac{1}{L{\varepsilon}}(x_{L{\varepsilon}}-x_{0}+\psi_{L,{\varepsilon}})$, and that using reversibility of the leapfrog integrator, we can also write $p_{L{\varepsilon}}=\frac{1}{L{\varepsilon}}(x_{L{\varepsilon}}-x_{0}-\psi_{L,{\varepsilon}}^{R})$. The log acceptance ratio can therefore be written $$U(x)-U(x_{L{\varepsilon}})+\frac{1}{2L^{2}{\varepsilon}^{2}}\left(\|x_{L{\varepsilon}}-x_{0}+\psi_{L,{\varepsilon}}\|^{2}-\|x_{L{\varepsilon}}-x_{0}-\psi_{L,{\varepsilon}}^{R}\|^{2}\right).$$ We require this quantity to be $\geq0$. This is equivalent to the requirement $$\int_{x_{0}}^{x_{L{\varepsilon}}}\nabla U(z)dz\leq\frac{1}{2L^{2}{\varepsilon}^{2}}\left(\|x_{L{\varepsilon}}-x_{0}+\psi_{L,{\varepsilon}}\|^{2}-\|x_{L{\varepsilon}}-x_{0}-\psi_{L,{\varepsilon}}^{R}\|^{2}\right).$$ We can re-write the right-hand side of the above expression as $$\frac{1}{2L^{2}{\varepsilon}^{2}}\left(2\langle x_{L{\varepsilon}}-x_{0},\psi_{L,{\varepsilon}}+\psi_{L,{\varepsilon}}^{R}\rangle+\|\psi_{L,{\varepsilon}}\|^{2}-\|\psi_{L,{\varepsilon}}^{R}\|^{2}\right),$$ and then note that $$\psi_{L,{\varepsilon}}+\psi_{L,{\varepsilon}}^{R}=\frac{L{\varepsilon}^{2}}{2}\left(\nabla U(x_{0})+\nabla U(x_{L{\varepsilon}})+2\sum_{i=1}^{L-1}\nabla U(x_{i{\varepsilon}})\right).$$ Substituting this into the inequality and simplifying gives the result. 

The requirement (\[eqn:ic2\]) is straightforwardly satisfied if $U(x)$ exhibits convexity in the tails, since the left-hand side of (\[eqn:ic2\]) will be negative in this case if $x_{L{\varepsilon}}\in I(x_{0})$ while the right-hand side will be positive.  Such convexity is not necessary, and there is some discussion in [@roberts1996exponential] of relaxations of (\[eqn:inwards\]) to the requirement that $\alpha(x_{0},x_{{\varepsilon}})\geq \delta$ for some $\delta >0$ if $\|x_{{\varepsilon}}\|\leq\|x_{0}\|$, which are also applicable to the HMC case and would relax the inequality (\[eqn:ic2\]) to some degree.  In essence, the key role of the ‘inwards acceptance’ property (\[eqn:inwards\]) (among the class of potentials which satisfy **A2**) is to limit the degree of oscillation in the tails of the density $e^{-U(x)}$, which can potentially mean that too many proposals $x_{L{\varepsilon}}$ for which the chosen Lyapunov function $V(x_{L{\varepsilon}})/V(x_{0})<1$ are rejected to establish a geometric bound of the form (\[eqn:geoerg\]).  Similar requirements to (\[eqn:inwards\]) are needed for many Markov chain Monte Carlo methods which rely on the Metropolis–Hastings construction \[e.g. [ @jarner2000geometric; @roberts1996geometric; @roberts1996exponential]\].  The issues are discussed in some detail in the case of the Random Walk Metropolis in [@jarner2000geometric].  It is possible that choosing the more natural (but less pliable) Lyapunov function $V(x)=e^{sU(x)}$ for some $s>0$ would remove the need for (\[eqn:inwards\]) here, owing to the ergodic nature of the proposal kernel.  We leave such explorations for future work. 

The preceding discussion leads to the following assumption that we require for geometric ergodicity here. 

**A3** [*The chain satisfies the ‘inwards acceptance’ property (\[eqn:inwards\]) which can equivalently be formulated as (\[eqn:ic2\]).*]{}  \[ass:A3\]

Assumptions **A1-A3** together are sufficient to establish a geometric bound. 

Necessary conditions for geometric ergodicity
---------------------------------------------

Next we highlight the importance of the growth assumptions we have made on the potential, by showing two general scenarios in which HMC will *not* produce geometrically ergodic Markov chains. 

### Light tails

We begin with the case where the gradient term may grow at a faster than linear rate, meaning that the resulting system of equations (\[eqn:hamilton\]) is ‘stiff’, in the sense that the derivatives can change very rapidly over small time scales, which can pose a challenge to explicit numerical integrators.  We show in Theorem \[thm:light\] that in this scenario a Markov chain produced by the HMC method can exhibit undesirable behaviour. 

\[thm:light\] If it holds that $$\label{eqn:explode}
\lim_{\|x\| \to \infty}\frac{\|\nabla U(x)\|}{\|x\|} = \infty,$$ and that there is a fixed $C <\infty$ such that whenever $\|y\|\geq 2\|x\|\geq C$ then $$\label{eqn:oscillations}
\|\nabla U(y)\| \geq 3\|\nabla U(x)\|,$$ and it also holds that $$\label{eqn:oscillation2}
\lim_{\|x\|\to\infty, \|y\| \geq 2^L\|x\|} \left( U(x) - U(y) - \frac{1}{2}\|x\|^2 \right) = -\infty,$$ then the Hamiltonian Monte Carlo method with fixed integration time $T = L{\varepsilon}$ does not produce a geometrically ergodic Markov chain for any choice $T>0$.

Lemmas \[lem:xelower\] and \[lem:xlower\] below establish that in this case $\|x_{L{\varepsilon}}\|\geq 2^L\|x_0\|$, and Lemma \[lem:arate\_0\] shows that this will result in $\alpha(x_0,x_{L{\varepsilon}})$ tending to zero as $\|x_0\|\to\infty$ provided $\|p_0\|\leq \|x_0\|^\delta$ for some $\delta<1$, allowing Proposition \[prop:supr\] to be envoked.  To conclude we simply note that $\mathbb{P}(\|p_0\|\leq \|x_0\|^\delta) \to 1$ as $\|x_0\| \to \infty$, establishing the result. 

The conditions (\[eqn:oscillations\]) and (\[eqn:oscillation2\]) limit the amount that the potential can oscillate as it approaches $\infty$, and are introduced to prevent tail oscillations in gradient from making the behaviour of the method too unpredictable to analyse sensibly.  They are very lenient and should be satisfied for the vast majority of statistical models of interest for which (\[eqn:explode\]) holds.  Below we establish several intermediate results, the first two of which relate to the values of $\|x_{L{\varepsilon}}\|$ when $\|x_0\|$ is large in this scenario. 

\[lem:xelower\] If (\[eqn:explode\]) holds then there exists an $\eta<\infty$ such that for all $\|x_0\|>\eta$ and any $\|p_0\|\leq \|x_0\|^\delta$ for some $\delta<1$, it holds that $\|x_{\varepsilon}\|>2\|x_0\|$.

Taking norms after a single leapfrog step gives $$\begin{aligned}
\|x_{\varepsilon}\| = \left\| x_0 - \frac{{\varepsilon}^2}{2}\nabla U(x_0) + {\varepsilon}p_0 \right\| &\geq \frac{{\varepsilon}^2}{2}\|\nabla U(x_0)\| - \|x_0\| - {\varepsilon}\|p_0\|.\end{aligned}$$ Dividing by $\|x_0\|$ gives $$\frac{\|x_{\varepsilon}\|}{\|x_0\|} \geq \frac{{\varepsilon}^2}{2}\frac{\|\nabla U(x_0)\|}{\|x_0\|} - 1 - {\varepsilon}\frac{\|p_0\|}{\|x_0\|}.$$ Using (\[eqn:explode\]), we can choose an $x_0$ such that the first term on the right-hand side is larger than $6/{\varepsilon}^2$, and the last term can be made negligibly small as $\|p_0\| \leq \|x_0\|^\delta$ for some $\delta < 1$, which establishes the result. 

\[lem:xlower\] If (\[eqn:explode\]) holds then there exists an $\eta<\infty$ such that for all $\|x_0\|>\eta$ and any $\|p_0\| \leq \|x_0\|^\delta$ for some $\delta < 1$, it holds that $\|x_{L{\varepsilon}}\| \geq 2^L\|x_0\|$.

We proceed iteratively.  First note that $$x_{2{\varepsilon}} = x_{\varepsilon}- {\varepsilon}^2\nabla U(x_{\varepsilon}) - \frac{{\varepsilon}^2}{2}\nabla U(x_0) + {\varepsilon}p_0.$$ Using this, we have $$\begin{aligned}
\frac{\|x_{2{\varepsilon}}\|}{\|x_{\varepsilon}\|} \geq \frac{{\varepsilon}^2}{2}\frac{\|\nabla U(x_{\varepsilon})\|}{\|x_{\varepsilon}\|} - \frac{{\varepsilon}^2}{2}\frac{\|\nabla U(x_0)\|}{\|x_{\varepsilon}\|} - 1 - {\varepsilon}\frac{\|p_0\|}{\|x_{\varepsilon}\|}.\end{aligned}$$ Showing the right-hand side is $\geq 2$ amounts to upper bounding the middle term, or equivalently lower bounding its reciprocal.  We have $$\label{eqn:inter21}
\frac{2 \|x_{\varepsilon}\|}{{\varepsilon}^2 \|\nabla U(x_0)\|} \geq \frac{ {\varepsilon}^2\|\nabla U(x_0)\|-2\|x_0\| - 2{\varepsilon}\|p_0\|}{{\varepsilon}^2\|\nabla U(x_0)\|} \geq 1 - \delta$$ for some $\delta>0$ which can be made arbitrarily small by choosing $\|x_0\|$ large enough. 

Next we have $$\frac{\|x_{3{\varepsilon}}\|}{\|x_{2{\varepsilon}}\|} \geq \frac{{\varepsilon}^2}{2}\frac{\|\nabla U(x_{2{\varepsilon}})\|}{\|x_{2{\varepsilon}}\|} - \frac{{\varepsilon}^2}{2}\frac{\|\nabla U(x_0)\|}{\|x_{2{\varepsilon}} \|} - {\varepsilon}^2\frac{\|\nabla U(x_{\varepsilon})\|}{\|x_{2{\varepsilon}}\|} - 1 - {\varepsilon}\frac{\|p_0\|}{\|x_{2{\varepsilon}}\|}.$$ Here the right-hand side will be $\geq 2$ provided the middle two terms can be bounded above.  For the first we lower bound the reciprocal, using (\[eqn:inter21\]) gives $$\frac{2}{{\varepsilon}^2}\frac{\|x_{2{\varepsilon}} \|}{\|\nabla U(x_0)\|} \geq \frac{2}{{\varepsilon}^2} \frac{ \|x_{\varepsilon}\|}{\|\nabla U(x_0)\|} \geq 1-\delta.$$ For the second we have $$\frac{\|x_{2{\varepsilon}}\|}{{\varepsilon}^2\|\nabla U(x_{\varepsilon})\|} \geq 1 - \frac{\|x_{\varepsilon}\|}{{\varepsilon}^2 \|\nabla U(x_{\varepsilon})\|} - \frac{\|\nabla U(x_0)\|}{2\|\nabla U(x_{\varepsilon})\|} - \frac{\|p_0\| }{{\varepsilon}\|\nabla U(x_{\varepsilon})\|}.$$ The second and last terms on the right hand side can be made arbitrarily small by choosing $\|x_0\|$ large enough.  Envoking (\[eqn:oscillations\]) gives $$\frac{\|\nabla U(x_0)\|}{2\|\nabla U(x_{\varepsilon})\|} \leq \frac{1}{6},$$ which therefore shows that $\|x_{3{\varepsilon}}\|\geq 2\|x_{2{\varepsilon}}\|$. An entirely analogous argument can be used to show that $\|x_{i{\varepsilon}}\|\geq 2\|x_{(i-1){\varepsilon}}\|$ for any fixed $i$, establishing the result. 

The next result shows that as a result of the fact that $\|x_{L{\varepsilon}}\|\geq 2^L\|x_0\|$ when $\|x_0\|$ is large enough, then the acceptance rate will approach $0$ in the limit as $\|x_0\| \to \infty$.

\[lem:arate\_0\] If (\[eqn:explode\]), (\[eqn:oscillations\]) and (\[eqn:oscillation2\]) hold then for any $\delta<1$ it holds that $$\lim_{\|x_0\|\to\infty,\|p_0\|\leq\|x_0\|^\delta} \alpha(x_0,x_{L{\varepsilon}}) = 0.$$

Recall that $$\begin{aligned}
\alpha(x_0,x_{L{\varepsilon}}) = 1 \wedge \exp\left( U(x_0) - U(x_{L{\varepsilon}}) + \frac{1}{2}\|p_0\|^2 - \frac{1}{2}\|p_{L{\varepsilon}}\|^2 \right).\end{aligned}$$ Note that $$\begin{aligned}
\|p_{L{\varepsilon}}\| &\geq \frac{{\varepsilon}}{2}\|\nabla U(x_{L{\varepsilon}})\| - {\varepsilon}\sum_{i=1}^{L-1} \|\nabla U(x_{i{\varepsilon}})\| - \frac{{\varepsilon}}{2}\|\nabla U(x_0)\| - \|p_0\|, \\
&\geq \frac{{\varepsilon}}{2}\|\nabla U(x_{L{\varepsilon}})\| - {\varepsilon}\sum_{i=1}^{L} \left( \frac{1}{3} \right)^i \|\nabla U(x_{L{\varepsilon}})\| - \|p_0\|, \\
&= \frac{{\varepsilon}}{2}\left( 1- 2\sum_{i=1}^L \left(\frac{1}{3}\right)^i \right) \|\nabla U(x_{L{\varepsilon}})\| - \|p_0\|,\end{aligned}$$ where (\[eqn:oscillations\]) is used for the second line.  The term inside the bracket can be bounded below by some fixed constant $\gamma_L >0$, for any fixed $L<\infty$. Squaring the result gives $$\begin{aligned}
\|p_{L{\varepsilon}}\|^2 &\geq \left( \frac{{\varepsilon}\gamma_L}{2}\|\nabla U(x_{L{\varepsilon}})\| - \|p_0\|\right)^2 \\ 
&= \frac{{\varepsilon}^2\gamma_L^2}{4}\|\nabla U(x_{L{\varepsilon}})\|^2 + \|p_0\|^2 - {\varepsilon}\delta_L\|\nabla U(x_{L{\varepsilon}})\|\|p_0\|,\end{aligned}$$ which implies that $$\|p_0\|^2 - \|p_{L{\varepsilon}}\|^2 \leq {\varepsilon}\gamma_L\|\nabla U(x_{L{\varepsilon}})\| \left( \|p_0\| - \frac{{\varepsilon}\gamma_L}{4}\|\nabla U(x_{L{\varepsilon}})\| \right).$$ Noting that $\|p_0\|\leq \|x_0\|^\delta$ and that for any $M<\infty$ we can choose an $\|x_0\|$ large enough that $$\frac{{\varepsilon}\gamma_L}{4}\|\nabla U(x_{L{\varepsilon}})\| \geq M\|x_{L{\varepsilon}}\|\geq 2^L\|x_0\|,$$ then it follows that $\|p_0\|^2 - \|p_{L{\varepsilon}}\|^2 \leq -\|x_0\|^2$. Using this, then simply envoking (\[eqn:oscillation2\]) gives the result. 

### Heavy tails

In the case where $\pi(x)$ has ‘heavier than exponential’ tails in some direction the HMC method can also exhibit slow convergence, as $\liminf_{\|x\| \to \infty}\|\nabla U(x)\| = 0$. Intuitively the problem here is that when $\|x\|$ is large then the gradient provides insufficient drift back into the ‘centre’ of the space, meaning the chain can exhibit random walk behaviour and hence convergence can be very slow.  Theorem \[thm:heavy\] makes this intuition rigorous. 

\[thm:heavy\] If $\|\nabla U(x)\| <M$ for all $x \in {\mathbf{X}}$, then a *necessary* condition for the Hamiltonian Monte Carlo method to produce a geometrically ergodic Markov chain is $$\int e^{s\|x\|}\pi(dx) < \infty$$ for some $s>0$.

From Proposition \[prop:tight\], it is sufficient to show that for any ${\varepsilon}>0$ there is a $\delta >0$ such that $Q(x,B_\delta(x)) > 1-{\varepsilon}$ for all $x \in {\mathbf{X}}$. Using equation (\[eqn:mprop\]) if $x_0$ is the current point in the chain then $$\|x_{L{\varepsilon}} - x_0\| = \left\| L{\varepsilon}p_0 - \frac{L{\varepsilon}^2}{2}\nabla U(x_0) - {\varepsilon}^2\sum_{i=1}^{L-1} (L-i)\nabla U(x_{i{\varepsilon}}) \right\|.$$ Applying the triangle inequality and then the global bound on $\|\nabla U(x)\|$ gives $$\|x_{L{\varepsilon}} - x_0\| \leq \frac{L{\varepsilon}^2}{2}M + \frac{M{\varepsilon}^2 L(L-1)}{2} + L{\varepsilon}\|p_0\|.$$ As $p_0$ follows a centred Gaussian distribution with fixed covariance then Chebyshev’s inequality gives the result. 

In fact, in this case the lack of geometric ergodicity is a property of the flow itself, rather than being a consequence of numerical instabilities as in Theorem \[thm:light\], as shown by the following result. 

\[prop:heavyexact\] Theorem \[thm:heavy\] still holds even if an exact integrator is available for Hamilton’s equations. 

Using Hamilton’s equations, we have $$x_T - x_0 = \int_0^T p_s ds = \int_0^T \left[ p_0 - \int_0^s \nabla U(x_u) du \right] ds.$$ Taking the norm and using the upper bound gives $$\begin{aligned}
\|x_T - x_0\| \leq T\|p_0\| - \int_0^T \int_0^s \|\nabla U(x_u)\| duds \leq T\|p_0\| + CT^2/2,\end{aligned}$$ and again Chebyshev’s inequality gives the result. 

The class of models for which $\|\nabla U(x)\|$ is bounded and $\int e^s\|x\|\pi(dx)<\infty$ for some $s>0$ is comparatively narrow, essentially comprising $U(x) = C\|x\| + b(x)$, where $C<\infty$ and $b:{\mathbf{X}}\to{\mathbb{R}}$ is some appropriately regular function which is bounded both above and below. 

Results for an position-dependent integration time {#sec:dynamic}
==================================================

An important free parameter in HMC is the integration time $T$, which we have previously assumed to be fixed.  The representation (\[eqn:mprop\]) does however suggest that allowing this to change can have some benefits.  If the candidate map is viewed as $$x_{L{\varepsilon}} = x_0 + \text{DRIFT}(x_0,p_0,T) + L{\varepsilon}p_0,$$ then if the ‘DRIFT’ function becomes negligible for large $\|x_0\|$ and fixed $T$, then it can be increased in magnitude by making $T$ larger.  In a trivial setting in which $\nabla U(x) = x^{-p}$ (for $x > 0$), for example, then setting $T = x^p$ would ensure that the drift no longer shrinks to zero as $\|x_0\|$ increases. 

We make this simple intuition rigorous for an idealised algorithm on the particular one-dimensional Exponential Family class of models with densities of the form $$\label{eqn:efsmooth}
\pi(x) \propto \exp \left( - \beta^{-1} (1+x^2) ^{\beta / 2} \right),$$ for some fixed $\beta>0$. Here any contour $C_{x_0,p_0} := \{ (x,p) : H(x,p) = H(x_0,p_0)\}$ consists of a single closed path, and the flow is periodic from any fixed starting point.  We additionally assume that the period length $\zeta_{x_0,p_0} >0$ is known, and that we have an exact integrator for Hamilton’s equations.  This means that we need not concern ourselves with the acceptance probability (we discuss this issue in Section \[sec:discussion\]). 

At iteration $i$ (with $x_0 = x_{i-1}$), the *dynamic* HMC implementation we consider consists of re-sampling $p_0 \sim N(0,1)$, and then setting $x_i = \text{Pr}_x \circ \varphi_{\tau}(x_0,p_0)$, where $\tau \sim U[0,\zeta_{x_0,p_0}]$. In words, we flow along the Hamiltonian for $\tau$ units of time, where $\tau$ is a uniform random variable with maximum value $\zeta_{x_0,p_0}$ (note that $\varphi_{\zeta_{x_0,p_0}}(x_0,p_0) = (x_0,p_0)$). 

Firstly, note that $\pi$-irreducibility is more straightforward to see here.  To reach any set $A \in {\mathcal{B}}$ with $\pi(A) > 0$, we first consider the single contour $C_{x_0,p_0}$, and specifically the component of this contour that is connected to $(x_0,p_0)$. Let $C_{x_0}$ be the projection of this component onto ${\mathbf{X}}$. Then any nonempty set $A' \subset C_{x_0}$ has positive probability of occuring, as the next point is chosen from a density with support all of $C_{x_0}$. As the contours are composed of single components, and cover the entire space, then for any $A$, the probability of choosing a contour for which this argument can be applied is greater than zero.  Figure \[fig:hmcperiod\] offers more intuition. 

(0,4) node \[left\] [$p$]{} – (0,0) – (5,0) node \[below\] [$x$]{}; (0,-4) – (0,0) – (-5,0); (0,0) circle \[radius=3\]; (-3,0) – (3,0); (1,2.7) – (1,0.1); (-2.1,2) – (-2.1,0.1); (-1,-2.7) – (-1,-0.1); (2.1,-2) – (2.1,-0.1); (2,-2.236) circle \[radius=0.1\]; at (2,-2.236) [$(x_0,p_0)$]{}; at (2, 2.1) [$C_{x_0,p_0}$]{}; at (-2,-0.05) [$C_{x_0}$]{};

We introduce some additional notation in this section.  We define the *microcanonical* expectation of a real-valued function $f(x_t,p_t)$, where $(x_t,p_t) = \varphi_t(x_0,p_0)$, i.e. the solution to (\[eqn:hamilton\]) for $t$ units of time initialised at $(x_0,p_0)$, as $$\left< f(x_0,p_0) \right> := \zeta_{x_0,p_0}^{-1} \int_0^{\zeta_{x_0,p_0}} f(x_s,p_s)ds.$$ This is simply the time expectation of $f$ from uniformly sampling across $C_{x_0,p_0}$.

To establish geometric ergodicity, we rely on conservation of the Hamiltonian, meaning that for any $T>0$ $$\label{eqn:conserveH0}
\int_0^T U(x_s) ds + \frac{1}{2}\int_0^T p_s^2 ds = H(x_0,p_0).$$ Averaging over initial momenta gives $$\int H(x_0,p_0) \mu^G(dp_0) = U(x_0) + 1/2.$$ We first introduce a result from the Physics literature \[e.g. [ @goldstein1965classical]\] which relates the kinetic and potential energies.  Using this we can relate the left hand side of (\[eqn:conserveH0\]) to $$PU(x_0) = \int U(y)P(x_0,dy),$$ where $P$ is the transition kernel under consideration.  Since the right-hand side of (\[eqn:conserveH0\]) relates to the current value of $U(x_0)$, then our goal will be to construct a suitable Lyapunov function from the potential energy that will help us establish the necessary drift condition. 

(Virial Theorem).  \[thm:virial\] Under Hamiltonian flow $(x_{s},p_{s}) = \varphi_s(x_0,p_0)$ we have $$\label{eqn:vt1}
\left< x_0\nabla U(x_0) \right> = \left< p_0^2 \right>.$$

Define the *virial* function $G_t = x_tp_t$. From the fundamental theorem of Calculus we have $$\left< \dot{G}_0 \right> = \frac{G_{\zeta_{x_0,p_0}} - G_0}{\zeta_{x_0,p_0}} = 0,$$ where $\dot{G}_t := dG_t/dt$. In this case $$\dot{G}_t = x_t \dot{p}_t + p_t\dot{x}_t = -x_t \nabla U(x_t) + p_t^2,$$ meaning $$\left< x_0 \nabla U(x_0) \right> = \left< p_0^2 \right>,$$ as required. 

We can now state and prove the main result of this section. 

\[thm:dynamic\] For the one-dimensional Exponential Family class of distributions with density given by (\[eqn:efsmooth\]), the dynamic Hamiltonian Monte Carlo method produces a geometrically ergodic Markov chain for any value of $\beta > 0$.

Note that by conservation of the Hamiltonian, we have $$\begin{aligned}
 \label{eqn:conserveh}
\int \left< U(x_0) + p_0^2/2 \right> \mu^G(dp_0)  &= \int H(x_0,p_0)\mu^G(dp_0) = U(x_0) + 1/2,\end{aligned}$$ Choose the Lyapunov function $V(x) = U(x) + xU'(x) + 1$. Using Theorem \[thm:virial\], we can re-write the above expression $$PV(x_0)  = U(x_0) + 3/2.$$ Note also that for any $\eta > 0$ there is an $M_\eta<\infty$ such that whenever $|x_0|>M_\eta$ $$(1+\eta)U(x_0) \geq U(x_0) + 3/2.$$ The proof will be complete if we can find a $\lambda < 1$ such that $(1+\eta)U(x_0) \leq \lambda V(x_0)$ for suitably large $|x_0|$. Now $x_0U'(x_0) \to \beta U(x_0)$ here as $|x_0| \to \infty$, meaning that there is an $M<\infty$ such that whenever $|x_0|>M$ $$x_0 U'(x_0)/2 \geq \beta U(x_0)/2 - 1.$$ Taking $|x_0|\geq \max(M_\eta,M)$ we can therefore re-write the inequality of interest $(1+\eta)U(x_0) \leq \lambda V(x_0)$ as $$(1+\eta)U(x_0) \leq \lambda (1+\beta/2) U(x),$$ which will be true if $$\lambda \geq \frac{1+\eta}{1+\beta/2}.$$ Choosing $\eta<\beta/2$ ensures $\lambda <1$ and also gives the desired inequality $PV(x_0) \leq \lambda V(x_0)$ whenever $|x_0| > \max(M_\eta,M)$, showing that the resulting Markov chain will be geometrically ergodic. 

Discussion {#sec:discussion}
==========

We have established conditions under which geometric ergodicity will and will not hold for Markov chains produced by the Hamiltonian Monte Carlo method.  Here we discuss how our results can be extended in various ways, as well as how they translate to standard implementations in widely used software \[[@carpenter2016stan]\]. 

Dynamic implementations
-----------------------

Allowing the integration time in HMC to depend on the current point in the chain without an exact integrator will typically mean that some adjustments to $\alpha(x_0,x_T)$ must be made to ensure that $\pi(\cdot)$ is still preserved.  The reason is that the approximate flow map $\varphi_T$ may no longer be reversible, as if $T_1 := T(x_0,p_0)$ and $T_2 := T(x_T,p_T)$ then $\varphi^{-1}_{T_2} \circ \varphi_{T_1}$ will typically not be the identity map if $T_1 \neq T_2$. The two possible ways of changing the integration time $T = L{\varepsilon}$ are to adjust either $L$ or ${\varepsilon}$. Increasing $L$ requires more computations per transition, while this is not necessarily true for ${\varepsilon}$. In the No-U-Turn sampler a binary tree approach is introduced to ensure preservation of detailed balance when $L$ is altered in different parts of the space \[[@hoffman-gelman:2013]\].  We are not aware of any implementations involving adjustment of ${\varepsilon}$, however it is likely that similar modifications to $\alpha(x_0,x_T)$ are possible here also.  We note that variable step-size Hamiltonian integrators are a well-explored notion in molecular dynamics \[e.g. [ @skeel1994symplectic]\].  Adjusting ${\varepsilon}$ may be a sensible option in some cases, as the leapfrog method is known to ‘almost’ preserve the modified Hamiltonian $$\tilde{H}(x,p) = H(x,p) + \left( \frac{1}{12} p^t \nabla^t\nabla U(x) p - \frac{1}{24}\nabla U(x)^t \nabla U(x) \right) {\varepsilon}^2 + O({\varepsilon}^4),$$ as shown for example in [@lopez1997explicit].  When $\pi(x)$ is not log-concave in the tails and hence the elements of $\nabla U$ and $\nabla^t \nabla U$ become negligible as $\|x\| \to \infty$, this implies that ${\varepsilon}$ can be increased for larger $\|x\|$ without compromising on numerical accuracy. 

Extension to other integrators
------------------------------

The fixed integration time results in Section \[sec:fixed\] refer specifically to the leapfrog integrator implementation of HMC (aside from Proposition \[prop:heavyexact\]).  It should be possible to use the same approach when analysing other explicit symplectic integrators, however for schemes which rely on implicit methods \[e.g. [ @girolami2011riemann]\] then composing multiple steps of the integrator as in Proposition \[prop:mprop\] cannot be done so cleanly.  Implicit methods are needed when the Hamiltonian is non-separable, and can often resolve stiffness issues such as those characterised in Theorem \[thm:light\]. 

To construct ergodicity results for the most general version of Hamiltonian Monte Carlo (i.e. without restricting attention to $\mathbb{R}^n$) we note that there are many ways to construct drift conditions in line with the purely geometric framework introduced in [@betancourt2014geometric].  We also point out that for the one-dimensional Exponential Family, choosing the Riemannian metric $G(x) = \|\nabla^2U(x)\|$ and employing the approach of [@girolami2011riemann] is mathematically equivalent to applying the transformation $x' = \text{sgn}(x)\|x\|^{\beta/2}$ with corresponding density $$\pi(x') \propto \|x'\|^{2/(\beta - 1)}\exp \left( -x'^2/\beta^2 \right).$$ This new density will have Gaussian tails for any $\beta >0$, suggesting a well-behaved sampler can be constructed.  Further discussion on the relationship between geometric Markov chain Monte Carlo methods and parameter transformations is given in [@livingstone2014information]. 

Direct comparison of HMC and MALA
---------------------------------

The fixed integration time results essentially state that HMC is geometrically ergodic whenever MALA is, which does not allow us to answer the question of which algorithm to use in a given scenario.  The results of Section \[sec:dynamic\], while idealised, give some indication that tuning the integration time dynamically when $\pi(\cdot)$ is not log-concave in the tails can be beneficial, however these do not consider the additional computational cost.  If acceptance probabilities are ignored, then the difference between the two algorithms is that $p_{\varepsilon}\sim N(0,I)$ in the case of MALA, while $p_{\varepsilon}= p_0 - {\varepsilon}\nabla U(x_0)/2 - {\varepsilon}\nabla U(x_{\varepsilon})/2$ in the case of HMC.  This translates to slightly different updates after $L$ steps of each method, but it is not clear that one strategy is more effective than another for a given class of problems and a fixed computational budget.  Establishing such an ordering of algorithms would be an interesting future research direction. 

Non-periodic orbits
-------------------

The assumption made in Section \[sec:dynamic\] that the flow is periodic is a strong one, and will in the main not hold when $d > 1$. Here flows will typically be quasi-periodic, though the Poincaré recurrence theorem can still be invoked to ensure that $x_T$ will be in some neighbourhood of $x_0$ for a finite $T$. The idea of numerical *exhaustions* introduced in [@betancourt2016identifying] is to then consider an *approximate* version of Theorem \[thm:virial\] in which $$\int_0^{T(x_0,p_0)} \dot{G}_s ds < \delta$$ for some (small) $\delta >0$. The equation can be used as a stopping criteria for $T(x_0,p_0)$, and is motivated by Theorem \[thm:dynamic\].  This approach and other similar integration time rules are possible, and these may give rise to more efficient implementations of HMC in practice. 

Theorem \[thm:dynamic\] has connections to the work of [@connor2009state] on using state-dependent drift conditions to establish subgeometric convergence rates for Markov chains.  It may be that such ideas can be extended to show that HMC with a fixed integration time produces a Markov chain that converges at a certain polynomial rate when $\pi(x)$ is not log-concave in the tails. 

Honest bounds
-------------

Geometric ergodicity is often referred to as a *qualitative* bound, as there are unknown constants which may prevent a geometrically ergodic scheme from being useful in some practical scenarios, and an explicit upper bound on the geometric rate $\rho$ is not established when using the techniques introduced in [@roberts1996geometric].  With some modifications, however, quantitative bounds can be constructed for geometrically ergodic chains \[[@jones2001honest; @rosenthal2002quantitative; @baxendale2005renewal]\].  We have refrained from doing this here, as while these bounds are ‘honest’ (in the sense of [@jones2001honest]), they are also often too conservative to be of use in practice \[[@jones2001honest; @durmus2015quantitative]\]. 

Quantitative bounds which are perhaps more informative are found using a different Wasserstein distance than the Total Variation special case in [@durmus2015quantitative] and [@eberle2014error] for a version of MALA, but these are restricted to specific forms of $\pi(\cdot)$ which consist of bounded perturbations from a Gaussian distribution.  Monte Carlo estimates for non-asymptotic quantitative bounds using the Ricci curvature approach of [@ollivier2009ricci] are applied to Hamiltonian Monte Carlo in [@seiler2014positive].  We note that the applicability of these bounds relies on the assumption of *positive curvature* in some Wasserstein distance for the underlying Markov chain.  When this distance is chosen to be Total Variation, then this is a strictly stronger condition than geometric ergodicity \[see Corollary 22 in [@ollivier2009ricci]\], so we feel that our results are a useful pre-cursor to understanding when these estimated bounds are informative in practice. 

In the case of MALA, when $\|\nabla U(x)\|$ grows at a faster than linear rate for large $\|x\|$ then it is shown in [@bou2012nonasymptotic] that useful inferences for functionals concentrated in the centre of the space can be made by setting a small enough value for ${\varepsilon}$. It is likely that the same analysis can be done with HMC, and that the result would be similar.  The disadvantage of choosing small step-sizes in MALA is that the gradient term in the proposal is often forced to be negligible in many parts of the space.  This could be counter-acted in HMC to some extent by fixing the integration time $T$ and adapting $L$ to compensate for choosing a smaller ${\varepsilon}$. We leave such explorations for future work. 

Practitioner guidelines
-----------------------

The main conclusion of our work for practitioners implementing the method in a bespoke manner is to consider the form of $\|\nabla U(x)\|$. If this term either grows very fast or becomes negligibly small when $\|x\|$ is large then it is likely that the Markov chains produced will struggle to explore the tails of $\pi(\cdot)$ effectively.  When the gradient grows at a faster than linear rate then a suitably small value for ${\varepsilon}$ must be chosen to counteract this, while when it shrinks then the integration time $T$ must be made sufficiently large.  Of course in either scenario if there is a re-parametrisation of the model that may not suffer these difficulties then this should be applied.  Users implementing the method in the Stan software \[[@carpenter2016stan]\] should note that both of these instances are captured by standard output diagnostics.  Large gradients are classed as ‘divergences’, while a failure to move far enough because of negligible gradients is recorded through the ‘maximum tree depth reached’ warning.  If this happens and $\pi(\cdot)$ is known to be proper then the user should set as large a maximum tree depth as is computationally feasible when tail exploration is of keen interest. 

Acknowledgements {#acknowledgements .unnumbered}
================

We thank Alexandros Beskos, Gareth Roberts, Krzysztof Łatuszyński and Gabriel Stoltz for useful discussions.  SL thanks Nawaf Bou–Rabee for pointing him to [@cances2007theoretical]. 

SL was supported by a PhD scholarship from Xerox Research Centre Europe and EPSRC grant EP/K014463/1 for this project.  SB was supported by EPSRC fellowship EP/K005723/1, MB is funded by EPSRC grant EP/J016934/1, and SB and MB also received a 2014 EPRSC NCSML Award for PDRA Collaboration for this project.  MG is funded by an EPSRC Established Career Research Fellowship, EP/J016934/1, and a Royal Society Wolfson Research Merit Award.