On the Geometric Ergodicity of Hamiltonian Monte Carlo

Introduction

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 . Iterative algorithms of this type are widely used in (for example) statistics and machine learning , inverse problems , and molecular dynamics . 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. ). 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, and many are well understood theoretically (e.g. ). HMC has proven an empirical success, with numerous authors noting its superior performance in a variety of settings (e.g. ) and high performance software available for its implementation . Comparatively few rigorous results, however, exist to justify this. Indeed, the absence of such analysis has been noted on more than one occasion . 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.

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 , in which the family $\{f_{\theta},\theta\in\Theta\}$ is given by

and set $r(x,y):=0$ otherwise. Then $\alpha(x,y):=1\wedge r(x,y)$. A more general definition is given in Proposition 1 of . The resulting chain $(X_{n})_{n\geq 0}$ is reversible with respect to $\pi(\cdot)$.

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 . 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 under suitable regularity conditions (see e.g. ). 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 .

The HMC method can also be considered within the above framework, as outlined in . The algorithm is designed to exploit the measure-preserving properties of Hamiltonian flow (e.g. ), which can be induced provided the state space for the chain is a symplectic manifold (e.g. ). 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

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. ). For $T=\varepsilon$ (meaning $L=1$) in fact HMC reduces to MALA. In general, however, for $L>1$ (2) 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.

Theoretical study of MCMC methods is in the main focused on two themes: convergence to equilibrium and asymptotic variance. 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 can be used to show that the distance to equilibrium decreases geometrically in $n$ (we elaborate in Section 3). 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. ). 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, 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, 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 . By contrast, we relate HMC to overdamped Langevin dynamics, as analysed in . 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. In some practical approximations are given for convergence bounds under a positive curvature assumption on the underlying chain. We discuss these further in Section 7. We comment further on connections between HMC and Langevin dynamics in the supplementary material .

Recently a HMC has been generalised to the context of sampling on spaces of infinite dimension . Due to the frequent singularity of measures in such spaces, it is often necessary to characterise distance to equilibrium here through other metrics than Total Variation. Such analysis is beyond the scope of this paper, though we note that recent work in the context of MALA in and are useful pre-cursors in this direction.

2 Notation

Overview of Main Results

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 position in the chain. In this scenario we have the following result.

If Assumptions A1 (on page 5.1), A2 (on page 5.11) and A3 (on page 5.13) hold, then a Markov chain produced by the Hamiltonian Monte Carlo method (outlined in Algorithm 1) 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. Instead of establishing $\pi$-irreducibility directly on on the multiple step HMC transition, we make a simple stochasticity assumption on the integration time parameter, which allows much of the technical difficulty to be sidestepped. Assumption A2 imposes conditions on the distribution from which expectations are desired, essentially restricting the tail behaviour to be lighter than a Laplacian but no lighter than a Gaussian distribution. 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 5.

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 (18) and (19) are satisfied

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.

For the exponential family class of models $\mathcal{E}(\beta,\alpha)$, under assumption A1 the following results hold:

(i) For $1\leq\beta\leq 2$, the Hamiltonian Monte Carlo method will produce a geometrically ergodic chain (for small enough $\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 . 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 will not 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 6, 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. ) which are commonly used in modern software packages (e.g. ). 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

then the idealised Hamiltonian Monte Carlo method introduced in Section 6 will produce a geometrically ergodic Markov chain for any choice of $\beta>0$.

The positive result in the case where $\beta>2$ is an artefact 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 7.

Preliminaries

The approach taken here to establishing geometric convergence was popularised in the monograph . A key observation first shown in that work is the following.

Recall that a set $C\in\mathcal{B}$ is called ‘small’ if there is a $t<\infty$, a measure $\nu(\cdot)$ defined on $(\mathbf{X},\mathcal{B})$ and an $\epsilon>0$ such that $\forall x\in C$ and $\forall A\in\mathcal{B}$ it holds that $P^{t}(x,A)\geq\epsilon\nu(A)$ (see e.g. ).

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

We say $P$ is of the Metropolis–Hastings type if

where $Q$ is a Markov kernel, $\alpha(x,y)$ is defined in (1) and $r(x)=1-\int\alpha(x,y)Q(x,dy)$.

The following was shown in when $P$ is of the form (4).

If $\pi(\cdot)$ and $Q(x,\cdot)$ admit Lebesgue densities $\pi(x)$ and $q(y|x)$, $\pi(x)$ is bounded away from and $\infty$ on compact sets, and there exists $\delta_{q}>0$ and $\epsilon_{q}>0$ such that, for every $x$,

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 3.3 are satisfied, then (3) holds if and only if

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

If $P$ is of Metropolis–Hastings type, then (3) 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 .

If for any $\eta>0$ there is a $\delta>0$ such that

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

Hamiltonian Monte Carlo

We give a brief introduction here. For a more detailed account see or . 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

Solving (8) results in Hamiltonian flow. To put this presentation into the framework introduced in Section 1, 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 (8) 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 (8) 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. ), meaning that for appropriate Hamiltonians the approximate solution of (8) 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 (7) 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

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. ). The full method is shown in Algorithm 1 below.

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

To use the techniques of , 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, which can be straightforwardly be derived using classical results (see e.g. ).

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

Proposition 4.2 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 (9) 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

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 6.

It is known (e.g. ) that establishing $\pi$-irreducibility is not so straightforward in the case of HMC as for Metropolis–Hastings methods based on random walks or Langevin diffusions. The canonical example where the system becomes reducible is integrating the harmonic oscillator over precisely one period (e.g. ). We show this in the supplementary material .

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

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

In fact, before the final publication of the present work, it was shown in that $\pi$-irreducibility can indeed be established without using assumption A1, but instead considering the HMC chain using a fixed number of leapfrog steps $L\geq 1$, under suitable assumptions and using appropriate techniques. We refer the interested reader to that work for details.

2 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

denote the ‘mean’ next candidate position ($x_{L\varepsilon}-L\varepsilon p_{0}$), and

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.

The HMC method produces a geometrically ergodic Markov chain if assumption A1 holds, and in addition both

where $\eta(d):=\Gamma((d+1)/2)/\Gamma(d/2)$, and

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$ and write $A(x_{0}):=R(x_{0})^{c}$. Then we can write

The last integral asymptotes to zero as $\|x_{0}\|\to\infty$ by (13). Writing $x_{L\varepsilon}(p_{0})$ to indicate that $x_{L\varepsilon}$ depends on $p_{0}$, the first integral can be written

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 (12) above then setting $\xi(x_{0}):=\sup_{\|p_{0}\|\leq\|x_{0}\|^{\delta}}(\|m_{L,\varepsilon}(x_{0},p_{0})\|-\|x_{0}\|)$ we can write

The last integral can be bounded above by the moment generating function of a Chi-distributed random variable with $d$ degrees of freedom, and so equals $e^{\log(1+s\sqrt{2}L\varepsilon\eta(d)+o(s))}\leq e^{s\sqrt{2}L\varepsilon\eta(d)+o(s)}$. Therefore by (12) the integral asymptotes to a quantity which is strictly less than one if $s>0$ is chosen to be suitably small.

Provided $\delta>1/2$, then for $\|x_{0}\|$ large enough $C\|x_{0}\|^{1-\delta}-\|x_{0}\|^{\delta}<-1$, meaning

which becomes negligibly small as $\|x_{0}\|\to\infty$, as required. ∎

Theorem 5.3 is a generalisation of Theorem 4.1 in to the HMC case. The nontriviality involved in this extension is accounting for the randomness induced into $m_{L,\varepsilon}(x_{0},p_{0})$ from $p_{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 (12) when $L\geq 2$. 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\geq 1$ (12) holds if the following conditions are met

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

(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$.

(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 5.10 and Proposition 5.11 below, which follow from the preceding Lemmas. ∎

The conditions of the result are intuitive. Condition (SC1.2) ensures that the gradient asymptotically ‘points inwards’, while (SC1.1) and (SC1.3) ensure that $\|\nabla U(x_{0})\|$ grows but at an asymptotically sublinear rate. We begin with a straightforward observation.

(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)<-\sqrt{2}\varepsilon\eta(d)$

(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$ then $\left(1+y\right)^{r}\leq 1+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

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

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.

A more intuitive condition which implies (SC1.2a) conditional on (SC1.3) and (SC1.1) is

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

For large enough $\|x_{0}\|$ this implies

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 conditions are sufficient for (12) to hold.

Under the following conditions (12) 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 5.7 to criteria that only depend on the current point $x_{0}$. The following lemmas give a starting point.

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

(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}$ ∎

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 5.8. 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. ∎

Under (SC1.1)-(SC1.3) then for any $L<\infty$ the conditions of Lemma 5.7 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

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 5.6, it is clear therefore that (ii) holds if

where each $c_{i}=(L-i)\varepsilon^{2}$ for $i\geq 1$ and $c_{0}=L\varepsilon^{2}/2$ . The second of these terms is $o(\|\nabla U(x^{*})\|\|x_{0}\|)$ using the Cauchy–Schwartz inequality and Lemma 5.9 (which shows that $\|x_{0}-x_{i\varepsilon}\|\in o(\|x_{0}\|)$), and so this term vanishes if divided by $\|\nabla U(x^{*})\|\|x_{0}\|$. The first term divided by the same quantity will be strictly positive as each term in the sum is $\geq 0$ using (SC1.3) and Lemma 5.9, and at least one of them is $>0$ since it will correspond to $x^{*}$. Using (SC1.1) establishes that $\langle\psi_{L\varepsilon},x_{0}\rangle/\|x_{0}\|\to\infty$ as $\|x_{0}\|\to\infty$, proving the result. ∎

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.

(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 5.7 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. ). 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 5.7 are satisfied.

Condition (SC1.1) precludes densities for which $\|\nabla U(x)\|\to c$ for some $0<c<\infty$. Often geometric ergodicity will still hold in this case, as we demonstrate in Corollary 2.3, however a different argument is required to that presented above.

It remains to consider (13), which reflects the role of the acceptance rate in the HMC method. We turn to this next.

2.2 Discussion of (13).

In the authors note that (13) 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})$

where $\hat{U}_{1}:=\left\langle x_{0}-x_{\varepsilon},\nabla U(x_{\varepsilon})+\nabla U(x_{0})\right\rangle/2$ denotes the ‘trapezium’ estimate for the line integral $U(x_{0})-U(x_{\varepsilon})=\int_{x_{\varepsilon}}^{x_{0}}\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 (13) if whenever $x_{L\varepsilon}\in I(x_{0})$ then in the limit as $\|x_{0}\|\to\infty$ it holds that

where $\hat{U}_{L}:=\langle x_{0}-x_{L\varepsilon},\nabla U(x_{0})+\nabla U(x_{L\varepsilon})+2\sum_{i=1}^{L-1}\nabla U(x_{i\varepsilon})\rangle/(2L)$ 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

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

We require this quantity to be $\geq 0$. This is equivalent to the requirement

We can re-write the right-hand side of the above expression as

Substituting this into the inequality and simplifying gives the result. ∎

The requirement (16) can sometimes be established using convexity arguments. In the exponential family class of Corollary 2.3, for example, setting $x^{i}:=x_{L\varepsilon}+i(x_{0}-x_{L\varepsilon})/L$, when $1\leq\beta<4/3$ one can show as $|x_{0}|\to\infty$ that $\hat{U}_{L}\xrightarrow{a.s.}(x_{0}-x_{L\varepsilon})(\nabla U(x_{0})+\nabla U(x_{L\varepsilon})+2\sum_{i=1}^{L-1}\nabla U(x^{i}))/(2L)$, the regular trapezium rule estimate for $\int_{x_{L\varepsilon}}^{x_{0}}\nabla U(z)dz=U(x_{0})-U(x_{L\varepsilon})$. Since $\nabla U(x)$ is concave/convex for $x$ positive/negative then the trapezium rule gives an underestimate for the integral as $|x_{0}|\to\infty$, and hence the left-hand side of (16) will be positive, while it is also possible to show that the right hand side is negative in this case (using arguments given in the proof of Corollary 2.3). We omit the details of this.

There is some discussion in of relaxations of (13) 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 (16) to some degree. In essence, the key role of the ‘inwards acceptance’ property (13) (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 (3). Similar requirements to (13) are needed for many Markov chain Monte Carlo methods which rely on the Metropolis–Hastings construction (e.g. ). The issues are discussed in some detail in the case of the Random Walk Metropolis in . 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 (13) 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 (13) which can equivalently be formulated as (16).

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

Part (ii) is a direct consequence of Theorem 2.2 For part (i), we consider three cases separately.

First consider $\beta\in(1,2)$, meaning $1>\beta-1>0$. Since $\nabla U(x)=\alpha\beta\text{sgn}(x)|x|^{\beta-1}$ then A2 holds. It remains to establish A3. We let $x_{0}\to\infty$ but an analogous argument holds as $x_{0}\to-\infty$ by symmetry. Note that $x_{\varepsilon}=x_{0}-\varepsilon^{2}\alpha\beta\text{sgn}(x_{0})|x_{0}|^{\beta-1}/2+\varepsilon p_{0}$ will clearly satisfy $(1-\delta)x_{0}<x_{\varepsilon}<x_{0}$ with probability reaching one in the limit for any $\delta>0$. Similarly $(1-\delta)x_{0}<x_{2\varepsilon}=x_{\varepsilon}-\varepsilon^{2}\alpha\beta\text{sgn}(x_{0})|x_{0}|^{\beta-1}-\varepsilon^{2}\alpha\beta\text{sgn}(x_{\varepsilon})|x_{\varepsilon}|^{\beta-1}/2+\varepsilon p_{0}<x_{\varepsilon}$ in the same asymptotic regime. Iterating the argument reveals that $(1-\delta)x_{0}<x_{L\varepsilon}<...<x_{\varepsilon}<x_{0}$ with probability tending to one as $x_{0}\to\infty$. Hence a.s. the proposal will be ‘inwards’, as will each intermediate point in the trajectory. To establish geometric convergence we must show that these inwards proposals are accepted with probability tending to one as $x_{0}\to\infty$. A Taylor series expansion of the difference in Hamiltonians for large enough $x_{0}$ gives

Since the leading order term is strictly positive then the result is proved. A detailed derivation is provided in the supplementary material .

In the case $\beta=1$ then as $x_{0}\to\infty$ the proposal in fact a.s. becomes $x_{L\varepsilon}=x_{0}-L\varepsilon^{2}/2+L\varepsilon p_{0}$, which resembles that of a random walk with inwards drift. Here the acceptance rate a.s. becomes one as the leapfrog integrator becomes exact provided the zero boundary is not crossed, and hence the scheme is geometrically ergodic following Theorem 16.0.1 and the argument of Section 16.1.3 in Chapter 16 of . Again a similar argument holds as $x_{0}\to-\infty$.

In the case $\beta=2$ following Example 3.5 in , setting $\theta:=\arccos(1-\alpha\varepsilon^{2})$ the proposal becomes $x_{L\varepsilon}=\cos(\theta L)x_{0}+\sin(\theta L)p_{0}/\sqrt{2\alpha(1-\alpha\varepsilon^{2}/2)}$,which will be inwards provided $|\cos(\theta L)|<1$, which will be true for suitably small $\varepsilon$. Similarly provided $p_{0}=o(x_{0})$ then the difference in Hamiltonian values will be

The $x_{0}^{2}$ coefficient will be positive provided $(1+2\alpha-\alpha^{2}\varepsilon^{2})\sin^{2}(\theta L)>0$ which will also be true for small enough $\varepsilon$, hence as $x_{0}\to\pm\infty$ A3 holds and since A2 does also then the result is proven. ∎

3 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.

We begin with the case where the gradient term may grow at a faster than linear rate, meaning that the resulting system of equations (8) 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 5.14 that in this scenario a Markov chain produced by the HMC method can exhibit undesirable behaviour.

and that there is a fixed $C<\infty$ such that whenever $\|y\|\geq 2\|x\|\geq C$ then

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$.

The conditions (18) and (19) 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 (17) 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.

If (17) 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

Using (17), 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. ∎

If (17) 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}\|$.

Showing the right-hand side is $\geq 2$ amounts to upper bounding the middle term, or equivalently lower bounding its reciprocal. We have

for some $\delta>0$ which can be made arbitrarily small by choosing $\|x_{0}\|$ large enough.

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 (21) gives

The second and last terms on the right hand side can be made arbitrarily small by choosing $\|x_{0}\|$ large enough. Envoking (18) gives

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 in the limit as $\|x_{0}\|\to\infty$.

If (17), (18) and (19) hold then for any $\delta<1$ it holds that

where (18) 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

Noting that $\|p_{0}\|\leq\|x_{0}\|^{\delta}$ and that for any $M<\infty$ we can choose an $\|x_{0}\|$ large enough that

then it follows that $\|p_{0}\|^{2}-\|p_{L\varepsilon}\|^{2}\leq-\|x_{0}\|^{2}$. Using this, then simply envoking (19) gives the result. ∎

3.2 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 5.18 makes this intuition rigorous.

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

From Proposition 3.6, 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 (9) if $x_{0}$ is the current point in the chain then

Applying the triangle inequality and then the global bound on $\|\nabla U(x)\|$ gives

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 5.14, as shown by the following result.

Theorem 5.18 still holds even if an exact integrator is available for Hamilton’s equations.

Taking the norm and using the upper bound gives

and again Chebyshev’s inequality gives the result. ∎

Results for an position-dependent integration time

An important free parameter in HMC is the integration time $T$, which we have previously assumed to be independent of the current position. The representation (9) does however suggest that allowing this to change can have some benefits. If the candidate map is viewed as

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. 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

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 7).

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^{\prime}\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. We provide a figure in the supplementary material to offer more intuition .

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 (8) for $t$ units of time initialised at $(x_{0},p_{0})$, as

This is simply the time expectation of $f$ from uniformly sampling across $C_{x_{0},p_{0}}$.

We first introduce a result from the Physics literature (e.g. ) which relates the kinetic and potential energies.

(Virial Theorem). Under Hamiltonian flow $(x_{s},p_{s})=\varphi_{s}(x_{0},p_{0})$ we have

Define the virial function $G_{t}=x_{t}p_{t}$. From the fundamental theorem of Calculus we have

where $\dot{G}_{t}:=dG_{t}/dt$. In this case

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

For the one-dimensional Exponential Family class of distributions with density given by (23), 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

Choose the Lyapunov function $V(x)=U(x)+xU^{\prime}(x)+1$. Using Theorem 6.1, we can re-write the above expression

Note also that for any $\eta>0$ there is an $M_{\eta}<\infty$ such that whenever $|x_{0}|>M_{\eta}$

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_{0}U^{\prime}(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$

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

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

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 .

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 . 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. Adjusting $\varepsilon$ may be a sensible option in some cases, as the leapfrog method is known to ‘almost’ preserve the modified Hamiltonian

as shown for example in . 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.

2 Extension to other integrators

The fixed integration time results in Section 5 refer specifically to the leapfrog integrator implementation of HMC (aside from Proposition 5.19). 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. ) then composing multiple steps of the integrator as in Proposition 4.2 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 5.14.

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 .

3 Honest bounds

Geometric ergodicity is often called a qualitative bound, as an explicit upper bound on the geometric rate $\rho$ is not established when using the techniques of . With some modifications, however, quantitative bounds can be constructed (e.g. ). We have refrained from doing this here, as these bounds are also often too conservative to be of use in practice .

Monte Carlo estimates for non-asymptotic quantitative bounds using the Ricci curvature approach of are applied to Hamiltonian Monte Carlo in . 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 ], 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 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, but we leave such explorations for future work.

4 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 should note that both of these instances are captured by standard output diagnostics. Numerical trajectories that become unstable due to 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

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

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, a Royal Society Wolfson Research Merit Award, and EPSRC grants EP/P020720/1, EP/J016934/3, EP/K034154/1.

References

Appendix A Examples of π𝜋\pi-irreducibility

The below example shows how the Hamiltonian Monte Carlo proposal transition can produce a method which is not $\pi$-irreducible, and hence will not be ergodic.

Take $\pi(x)\propto e^{-x^{2}/2}$, meaning $\nabla U(x)=x$, and set $L=2$. Then the HMC proposal becomes

Setting $\varepsilon=\sqrt{2}$ means $2\varepsilon-\varepsilon^{3}=0$, so that

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)$.

The following diagram give intuition for the $\pi$-irreducibility argument of the idealised Hamiltonian Monte Carlo method.

Appendix B Connections between HMC and Langevin dynamics

Hamiltonian Monte Carlo is based on interspersing Hamiltonian dynamics given by the equations

with intermittent re-sampling of $p\sim N(0,I)$ to inject stochasticity into the system. After a small time period the $x$-coordinate will be

Taking $\sqrt{\delta t}\ll 1$ then for suitably regular $\nabla U(x)$ one can approximate this with the expression

Hence, if the dynamics are only performed for a short period before momentum re-sampling, the dynamics will be close to those of an overdamped Langevin diffusion described by the stochastic differential equation

If instead $T$ is typically large in between momentum refreshments, then one can instead consider underdamped Langevin dynamics, described by the system

for some $\gamma>0$. This system more obviously relates to HMC, since it consists of a Hamiltonian part combined with some stochasticity, which is also injected into the momentum variable $V_{t}$. The stochasticity is however in this case continuously injected in the form of an Ornstein–Uhlenbeck (OU) process. As such, as the dynamics evolve the conservative Hamiltonian flow is constantly perturbed by small random adjustments to the momentum. This is closely connected to the behaviour of generalized HMC, in which integration times are typically smaller and the momentum is only partially refreshed using $p_{new}\sim N(\xi p_{old},(1-\xi^{2})I)$. Indeed, if $\xi:=e^{-\gamma t}$ then this step is an exact solution to the OU part of the system.

Appendix C Proof of inwards convergence for the Exponential family model class

We provide a verbose proof in the case $\beta\in(1,2)$, to elaborate on the short version provided in the main text. In what follows $U(x):=\alpha|x|^{\beta}$ for some $\alpha>0$, and $U^{(k)}(x):=d^{k}U(x)/dx^{k}.$ We precede the main result with a technical lemma.

For every $i\in\{1,...,L\}$, if $p_{0}=o(x_{0}^{\beta-1})$ then