Kinetic energy choice in Hamiltonian/hybrid Monte Carlo

TL;DR

This paper investigates how different kinetic energy choices in Hamiltonian Monte Carlo influence performance, revealing that standard Gaussian momentum may not always be optimal and highlighting a trade-off between efficiency and robustness.

Contribution

It introduces two evaluative quantities for kinetic energy choices and analyzes their impact on integrator stability, convergence, and the undesirable negligible moves property.

Findings

Heavy-tailed momentum distributions can cause negligible moves.

Standard Gaussian momentum is not always optimal.

A trade-off exists between efficiency and robustness.

Abstract

We consider how different choices of kinetic energy in Hamiltonian Monte Carlo affect algorithm performance. To this end, we introduce two quantities which can be easily evaluated, the composite gradient and the implicit noise. Results are established on integrator stability and geometric convergence, and we show that choices of kinetic energy that result in heavy-tailed momentum distributions can exhibit an undesirable negligible moves property, which we define. A general efficiency-robustness trade off is outlined, and implementations which rely on approximate gradients are also discussed. Two numerical studies illustrate our theoretical findings, showing that the standard choice which results in a Gaussian momentum distribution is not always optimal in terms of either robustness or efficiency.

Figures

Tables

Table 1. Table 1 : Effective sample sizes at equilibrium and iterations until max i ⁡ | ψ i | < 2 subscript 𝑖 subscript 𝜓 𝑖 2 \max_{i}|\psi_{i}|<2 for Hamiltonian Monte Carlo on the Ginzburg–Landau model with α = 0.1 𝛼 0.1 \alpha=0.1 , λ = 0.5 𝜆 0.5 \lambda=0.5 and τ = 2 𝜏 2 \tau=2 . Step sizes were chosen based on optimizing effective sample sizes. In each simulation the number of leapfrog steps was set to 10, and results are averaged over 10 runs.

	Effective sample size			Iterations to centre
	Minimum	Mean	Maximum
Gaussian	6,251	8,748	10,000	N/A
Relativistic Power $\beta =4/3$	5,253	6,777	8,271	4.2
Relativistic	3,591	4,639	5,525	8.6
Exponential Power $\beta =4/3$	810	1,108	1,303	11.9

Equations

$$\dot{x}=\nabla K(p),\dot{p}=-\nabla U(x).$$

$$x^{\ast }=x_0+{\int }_0^T\nabla K\left\{p_0-{\int }_0^s\nabla U(x_u)du\right\}ds.$$

$$x_{\epsilon }=x_0+\epsilon \nabla K\left\{p_0-\frac{\epsilon }{2}\nabla U(x_0)\right\},$$

$$p_{\epsilon }=p_0-\frac{\epsilon }{2}\left\{\nabla U(x_0)+\nabla U(x_{\epsilon })\right\}.$$

$$\nabla K\circ \nabla U(x),$$

$$\Vert P^n(x,\cdot )-\pi (\cdot ){\Vert }_{TV}\le CV(x){\rho }^n$$

$$\underset{\Vert x\Vert \to \infty }{\lim }\frac{\Vert \nabla K\circ \nabla U(x)\Vert }{\Vert x\Vert }=\infty ,$$

$$\Vert y\Vert \ge M\Vert x\Vert ⟹\Vert \nabla U(y)\Vert \ge \varphi (M)\Vert \nabla U(x)\Vert .$$

$$△(x_0,p_0)\to \infty ⟹△H(x_0,p_0)\to \infty .$$

$$\frac{\Vert \nabla K\{\delta \nabla U(x)\}\Vert }{\Vert \nabla K\circ \nabla U(x)\Vert }\ge C_{\delta }$$

$$\Vert p\Vert \ge \Vert y\Vert ⟹\Vert \nabla K(p)\Vert \ge \Vert \nabla K(y)\Vert ,$$

$$|p_i|\ge |y_i|⟹|k^{\prime }(p_i)|\ge |k^{\prime }(y_i)|.$$

$$\Vert \nabla U(x)\Vert \le A\Vert x{\Vert }^q+B,\Vert \nabla K(p)\Vert \le C\Vert p{\Vert }^{\frac{1}{q}}+D,$$

$$\Vert x_{L\epsilon }(x_0,p_0)\Vert \le E_L\Vert x_0\Vert +F_L,\Vert p_{L\epsilon }(x_0,p_0)\Vert \le {\mathbf{G}}_L\Vert x_0{\Vert }^q+{\mathbf{H}}_L.$$

$$\underset{\Vert x\Vert \to \infty }{\lim }P\{x,B_{\delta }(x)\}=1.$$

$$x^i=\max (x^{i-1}-1+{\xi }^i,1),{\xi }^i\sim N(0,1),$$

$$H(x,p)={\alpha }^{-1}|x{|}^{\alpha }+{\beta }^{-1}|p{|}^{\beta },$$

$$\dot{x}$$

$$\dot{p}$$

$$\mathcal{P}(E)=\underset{t>0}{\inf }\left\{t:(x_t,p_t)=(x_0,p_0)\right\}$$

$$\mathcal{P}(E)=c_{\alpha ,\beta }{E}^{\eta },\eta =\frac{1-(\beta -1)(\alpha -1)}{\alpha \beta },$$

$$\nabla K(p)\approx \nabla U^{-1}\left(p\right),$$

$$U(x)=C(\alpha x^4-\beta x^2),$$

$$U(\beta )=L(\beta )+\lambda \sum _{j=1}^d|{\beta }_j{|}^q,$$

$$\nabla U(x_1,...,x_9,\nu )={\left(2x_1{e}^{-\nu },...,2x_9{e}^{-\nu },\frac{\nu }{2}+\frac{9}{2}-\frac{e^{-\nu }}{2}\sum _i{x}_i^2\right)}^T.$$

$$\nabla K\left\{p_0-\frac{\epsilon }{2}\nabla U(x_0)\right\}=\nabla K\left\{-\frac{\epsilon }{2}\nabla U(x_0)\right\}+{\nabla }^{2}K\left\{-\frac{\epsilon }{2}\nabla U(x_0)\right\}{p}_0+\cdots$$

$$\Vert {\nabla }^{2}K\left\{\frac{\epsilon }{2}\nabla U(x)-p\right\}-{\nabla }^{2}K\left\{\frac{\epsilon }{2}\nabla U(x)\right\}\Vert \to 0$$

$$\underset{\Vert p\Vert \to \infty }{\lim }\Vert {\nabla }^{2}K(p)\Vert =0.$$

$$K(p)=\sum _i{m}_i{c}_i^2{\left(1+\frac{p_i^2}{m_i^2{c}_i^2}\right)}^{\frac{1}{2}}$$

$$K(p)=\sum _i{\beta }^{-1}{\left(1+{\gamma }_i^{-1}{p}_i^2\right)}^{\beta /2},$$

$$K(p)\approx \sum _i{\beta }^{-1}(1+{\gamma }_i^{-1}{p}_i^2)⟹{\partial }_{i}K(p)\approx 2(\beta {\gamma }_i{)}^{-1}{p}_i$$

Full text

1 Introduction

Hamiltonian Monte Carlo is a Markov chain Monte Carlo method which is now both widely used in Bayesian inference, and increasingly studied and developed. The idea of the approach is to use the deterministic measure-preserving dynamics of Hamiltonian flow to promote fast exploration of a parameter space of interest. To achieve this the space is augmented with a momentum variable, and the Markov chain evolves by switching between re-sampling this momentum and solving Hamilton’s equations for a prescribed amount of time. Typically the equations cannot be solved exactly and so a time-reversible and volume-preserving numerical integrator is used, with discretisation errors controlled for using a Metropolis step. Comprehensive reviews are given in Neal (2011)* and Betancourt et al. (2017)**.*

The free variables in Hamiltonian Monte Carlo are the time for which Hamilton’s equations should be solved between momentum refreshments, the choice of numerical integrator and step-size, and the choice of distribution for the momentum. Typically the Störmer–Verlet or leapfrog numerical integrator is used, which gives a reasonable balance between energy preservation and computational cost. Guidelines for step-size choice are given in Beskos et al. (2013), which is typically tuned to reach a 65-80% Metropolis acceptance rate. Heuristics exist for the integration time (Hoffman & Gelman, 2014), and there is some justification for a stochastic choice given in Bou-Rabee et al. (2017).

This paper is concerned with the choice of momentum distribution $\nu(\cdot)$, addressing a question which has been raised previously (Barthelmé & Chopin, 2011; Stephens, 2011). We will assume throughout that $\nu(\cdot)$ possesses a density which is proportional to $\exp\{-K(p)\}$ for some $K(p)$, which we call the kinetic energy of the system. The standard choice is $K(p)=p^{T}p/2$, with the resulting $\nu(\cdot)$ a Gaussian distribution. The general requirements are simply that $K(p)$ is differentiable, symmetric about zero and that $\nu(\cdot)$ can be sampled from. Alternative options to the Gaussian have recently been suggested (Lu et al., *2017; *Zhang et al., 2016)**. Here we consider how such a choice affects the algorithm as a whole, in terms of stability and convergence, and develop guidelines for practitioners.

Our key findings are that the robustness and efficiency of the method can to some degree be controlled through the quantity $\nabla K\circ\nabla U(x)$, notation defined below, which we term the composite gradient. In particular, we propose that balancing the tails of the kinetic energy with those of the potential to make this approximately linear in $x$ should give good performance, and that no faster than linear growth is necessary for algorithm stability. When $\pi(\cdot)$ is very light-tailed, this can be done by choosing $\nu(\cdot)$ to be heavier-tailed. There are, however, serious disadvantages to choosing any heavier tails than those of a Laplace distribution, which can result in the sampler moving very slowly in certain regions of the space. We introduce a negligible moves property for Markov chains to properly characterize this behaviour. We also find that in practice considering the distribution of $\nabla K(p)$, which we term the implicit noise, is important, since this governs the behaviour of the sampler in regions where $\|\nabla U(x)\|$ is small. We suggest various choices for controlling these two quantities, and support our findings with a numerical study. We also comment on how changing the kinetic energy can affect behaviour when $\nabla U(x)$ is approximated, either through subsampling or within an exact-approximate Monte Carlo scheme.

1.1 Setup, notation and assumptions

Throughout we denote the distribution from which expectations are desired by $\pi(\cdot)$. We assume that $\pi(\cdot)$ possesses a density which is proportional to $\exp\{-U(x)\}$, for some $C^{1}(\mathbb{R}^{d})$ potential $U:\mathbb{R}^{d}\to[0,\infty)$, where $x\in\mathbb{R}^{d}$. The Hamiltonian is formed as $H(x,p)=U(x)+K(p)$, where $p\in\mathbb{R}^{d}$ and $K:\mathbb{R}^{d}\to[0,\infty)$ is a continuously differentiable kinetic energy, and the corresponding Hamiltonian dynamics are given by

[TABLE]

We restrict ourselves to separable Hamiltonians, meaning the distribution for $p$ does not depend on $x$, as in this case general purpose and well-studied explicit numerical integrators are available (Leimkuhler & Reich, 2004). We also focus on the leapfrog numerical integrator (Neal, 2011) and a choice of integration time $T=L\varepsilon$, where $\varepsilon$ is the integrator step-size and $L$ is the number of leapfrog steps. We assume that $L$ is drawn from some distribution $\Psi(\cdot)$, which is independent of the current $x$ and $p$, and denote the result of solving (1) using this integrator for $T$ units of time given the starting position and momentum $(x,p)$ as $\tilde{\varphi}_{T}(x,p)$. With these tools, the Hamiltonian Monte Carlo method is described as Algorithm 1 below.

*We work on the Borel space $(\mathbb{R}^{d},\mathcal{B})$ and product space $\mathbb{R}^{2d}$ with the product $\sigma$-field. Throughout let $\|x\|,\|x\|_{1}$ and $\|x\|_{\infty}$ be the Euclidean, $L_{1}$ and $L_{\infty}$ norms of $x\in\mathbb{R}^{d}$, $B_{r}(x)=\{y\in\mathbb{R}^{d}:\|y-x\|<r\}$ be the Euclidean ball of radius $r$ centred at $x$, and $\partial_{i}=\partial/\partial x_{i}$. For functions $f,g:\mathbb{R}^{d}\to\mathbb{R}^{d}$, let $f\circ g(x)=f(g(x))$. Whenever a distribution is referred to as having a density then this with respect to the relevant Lebesgue measure. We write $x^{i}$ for the *$i$*th point in a Markov chain produced by Hamiltonian Monte Carlo, and $x_{i\varepsilon}$ for the *$i$th step of the leapfrog integrator within a single iteration of the method.

2 General observations

The effects of changing the proposal input noise distribution are more complex in Hamiltonian Monte Carlo than in the random walk Metropolis and the Metropolis-adjusted Langevin algorithm, two methods for which this question has been studied previously (Jarner & Roberts, 2007; Stramer & Tweedie, 1999). For both of these cases the input noise is combined linearly with some deterministic function of the current point of the Markov chain. As a result, larger values of this noise will typically result in larger proposed moves, which may be advantageous if $\pi(\cdot)$ has heavy tails or possesses multiple modes. By contrast, in Hamiltonian Monte Carlo the choice of kinetic energy alters both the input noise distribution and the Hamiltonian dynamics, and so the impact of different choices is not so transparent.

If we consider how a single proposal in Hamiltonian Monte Carlo perturbs the current position $x^{i-1}$ at iteration $i$, then setting $x_{0}=x^{i-1}$, $p_{0}\sim\nu(\cdot)$ and solving (1) for $T$ units of time gives

[TABLE]

An explicit numerical method typically involves the approximation $\int_{0}^{h}\nabla U(x_{u})du\approx\nabla U(x_{0})h$ for some suitably chosen $h$. In the case of the leapfrog integrator the approximate solution over a single leapfrog step becomes

[TABLE]

with the corresponding momentum update

[TABLE]

Recalling that $K(p)$ is an even function, meaning $\nabla K(p)$ is odd, then it will turn out that in many cases the function that intuitively governs the speed of the $x$-dynamics when $\|\nabla U(x)\|$ is large is

[TABLE]

which we refer to from this point as the composite gradient. Although in practice the choice of $\varepsilon$ will influence sampler behaviour, we will demonstrate that a qualitative understanding can be developed by considering only (5). Similarly, when $\|\nabla U(x)\|$ is small, $x_{\varepsilon}$ will resemble a random walk proposal with perturbation $\varepsilon$ multiplied by $\nabla K(p_{0})$, which we refer to as the implicit noise. We argue that an appropriate choice of kinetic energy is one for which these two quantities are both suitably optimized over, and in the next sections we examine both when this can be done and what a suitable choice for each should be.

3 The composite gradient

3.1 Geometric convergence and long term stability

We begin with some general results relating the choice of kinetic energy to the composite gradient, and potential consequences for the convergence rate of the method.

Definition 1**.**

We call a distribution which has a density $f(x)\propto\exp\{-g(x)\}$ heavy-tailed if $\lim_{\|x\|\to\infty}\|\nabla g(x)\|=0$, and light-tailed if $\lim_{\|x\|\to\infty}\|\nabla g(x)\|=\infty$.

Distributions which do not fall into either category are those for which the tails are heavy in some directions and light in others, for which the sampler will behave very differently in distinct regions of the space. A detailed study of these is beyond the scope of this work, though we conjecture that the least favourable region would dictate convergence. We also do not include distributions with exponentially decaying densities, which are often treated as a special case when analysing Metropolis–Hastings methods (Jarner & Hansen, 2000).

Recall that a $\pi$-invariant Markov chain with transition kernel $P$ is called geometrically ergodic if for some positive constants $C<\infty$, $\rho<1$ and a Lyapunov function $V:\mathbb{R}^{d}\to[1,\infty)$ the bound

[TABLE]

can be constructed, where $\|\mu(\cdot)-\nu(\cdot)\|_{TV}$ denotes the total variation distance between two measures $\mu(\cdot)$ and $\nu(\cdot)$. See for example Roberts & Rosenthal (2004)* for more details.*

Proposition 1**.**

In either of the following cases Algorithm 1 will not produce a geometrically ergodic Markov chain:

(i) $\pi(\cdot)$ is heavy-tailed

(ii) $\pi(\cdot)$ is light-tailed and the following conditions hold:

The composite gradient satisfies

[TABLE] 2. 2.

There is a strictly increasing unbounded function $\phi:[0,\infty)\to[1,\infty)$ and $m_{1},m_{2}<\infty$ with $m_{1}\geq 1$ such that for every $M\geq m_{1}$ and $\|x\|\geq m_{2}$

[TABLE] 3. 3.

Setting $\triangle H(x_{0},p_{0})=H(x_{L\varepsilon},p_{L\varepsilon})-H(x_{0},p_{0})$ and $\triangle(x_{0},p_{0})=\|x_{L\varepsilon}\|-\|x_{0}\|+\|p_{L\varepsilon}\|-\|p_{0}\|$,

[TABLE] 4. 4.

For all $\delta<\infty$, there is a $C_{\delta}>0$ such that for all $x\in\mathbb{R}^{d}$

[TABLE] 5. 5.

Either $K(p)$ satisfies

[TABLE]

or there is a $k:\mathbb{R}\to\mathbb{R}$ such that $K(p):=\sum_{i=1}^{d}k(p_{i})$ and

[TABLE]

for all $i$.

Remark 1**.**

Condition 4 in part (ii) will be satisfied if, for example, $\nabla K$ is a homogeneous function. It is satisfied for all choices considered in this work.

For the standard quadratic choice of kinetic energy the Hamiltonian Monte Carlo method will not produce a geometrically ergodic Markov chain when $\pi(\cdot)$ has heavy-tails (Livingstone et al., 2016). Part (i) of the above states that no choice of kinetic energy which is independent of $x$ can rectify this. We explain why in subsection 5.3. Part (ii) states that provided the density is sufficiently regular that (9) holds, and that the degree of tail oscillation in $\|\nabla U(x)\|$ is restricted through (8), then the composite gradient should not grow faster than linearly in $\|x\|$ for geometric ergodicity. Within the class of kinetic energies for which this requirement is met, we motivate an optimal choice in subsection 3.3.

The composite gradient was deduced based on studying the leapfrog integrator over a single step. Next, however, we argue that controlling this quantity in fact induces stability over multiple leapfrog steps.

Proposition 2**.**

If for some $q>0$ and $A,B,C,D<\infty$

[TABLE]

then $\limsup_{\|x\|\to\infty}\|\nabla K\circ\nabla U(x)\|/\|x\|<\infty$, and furthermore if $\|p_{0}\|\leq E_{0}\|x_{0}\|^{q}+F_{0}$ for some $E_{0},F_{0}<\infty$, then writing $x_{L\varepsilon}$ and $p_{L\varepsilon}$ as functions of $x_{0}$ and $p_{0}$, it holds that

[TABLE]

for some $E_{L},F_{L},\mathbf{G}_{L},\mathbf{H}_{L}<\infty$ which do not depend on $x_{0}$.

As a contrast, if the assumption on $\nabla K(p)$ was replaced with $\|\nabla K(p)\|\leq C\|p\|+D$, then the bound would read $\|x_{L\varepsilon}(x_{0},p_{0})\|=E_{L}\|x_{0}\|^{1\vee q^{L}}+F_{L}$, which grows exponentially in $L$ when $q>1$.

3.2 Kinetic energies that induce negligible moves

Next we consider a heavy-tailed choice of $\nu(\cdot)$. In this instance the composite gradient will in fact decay as $\|\nabla U(x)\|$ grows, meaning that if $\pi(\cdot)$ is light-tailed the resulting Markov chain will move very slowly in some regions of the space. We formalize this intuition below.

Definition 2**.**

We say that a Markov chain on $\mathbb{R}^{d}$ with transition kernel $P$ possesses the negligible moves property if for every $\delta>0$

[TABLE]

Proposition 3**.**

If $\pi(\cdot)$ is light-tailed and $\nu(\cdot)$ is heavy-tailed then a Markov chain produced by the Hamiltonian Monte Carlo method possesses the negligible moves property.

In general this property seems undesirable, the implication being that if the current point $x$ has large norm then the chain is unlikely to move very far in any appreciable number of steps. It does not, however, preclude geometric ergodicity, as often such steps can be close to deterministic in a desirable direction. We give a simple example to illustrate the phenomenon.

Example 1**.**

Consider the Markov chain with transition

[TABLE]

which is geometrically ergodic as shown in section 16.1.3 of Meyn & Tweedie (1993). Then the transformed Markov chain $(y^{i})_{i\geq 1}$, where each $y^{i}=\log(x^{i})$, is also geometrically ergodic, but possesses the negligible moves property.

Although a geometrically decaying convergence bound can still be established in such cases, the slow movement of the chain when $\|x\|$ is large is accounted for in the Lyapunov function $V(x)$ present in the bound (6).

Proposition 4**.**

If (14) holds for a Hamiltonian Monte Carlo method then any Lyapunov function $V:\bm{X}\to[1,\infty)$ cannot satisfy either of the following:

(i) $\log V$ is uniformly lower semi-continuous, meaning for all $\epsilon^{\prime}>0$, there is a $\delta^{\prime}<\infty$ such that if $\|y-x\|\leq\delta^{\prime}$, uniformly in $x$, then $\log V(y)\geq\log V(x)-\epsilon^{\prime}$.

(ii) There is an $s\in(0,\infty)$ such that $\liminf_{\|x\|\to\infty}V(x)e^{-s\|x\|}=c$ for some $c\in(0,\infty)$.

The above implies that Lyapunov functions will typically exhibit faster than exponential growth in $\|x\|$, meaning the penalty in the bound (6) for starting the chain in the tails of the distribution will be very large. This reflects the fact that the chain will take a long time to move an appreciable distance and hence approach equilibrium. So although the negligible moves property does not preclude geometric convergence, we regard it as an undesirable feature.

3.3 Illustrative example

Part (ii) of Proposition 1 relates to the numerical solution of (1) using the leapfrog method. If the composite gradient grows faster than linearly in $\|x\|$ then the dynamical system exhibits stiffness, meaning gradients change very quickly and hence numerical approximations in which they are assumed to be constant over small time periods will become unstable. Two possible remedies are to use a more complex implicit numerical method or adaptive step-size control, as discussed in Okudo & Suzuki (2016)**. With the standard choice of kinetic energy $\nabla K\circ\nabla U(x)=\nabla U(x)$, meaning these instabilities occur when $\|\nabla U(x)\|$ grows faster than linearly in $\|x\|$. Since the composite gradient dictates this phenomenon in general, fast growth in $\nabla U(x)$ can be counteracted by choosing a kinetic energy which is subquadratic in $\|p\|$. Physically, choosing a slowly growing kinetic energy slows down the Hamiltonian dynamics, meaning the system (1) is no longer stiff in the $x$-coordinate and hence simple numerical schemes like the leapfrog method become stable. We can further use the same idea when $\nabla U(x)$ exhibits relatively slow growth by choosing a fast-growing kinetic energy which will speed up the flow allowing the chain to explore the state space more efficiently.

To demonstrate this phenomenon more concretely we consider the family of separable Hamiltonians $H:\mathbb{R}^{2}\to[0,\infty)$ given by

[TABLE]

for suitable choices $\alpha,\beta>1$, which is necessary to ensure continuity of the derivative. Hamilton’s equations are

[TABLE]

The flow is periodic since the contours of $H$ are closed curves in $\mathbb{R}^{2}$. These kinetic energies correspond to choosing $\nu(\cdot)$ from the exponential power family of Box & Tiao (1973), for which direct sampling schemes are provided in Mineo (2014). Zhang et al. (2016)* recommend choosing this class of momentum distributions in Hamiltonian Monte Carlo. The following proposition shows how the the period length depends on the initial value of the Hamiltonian.*

Proposition 5**.**

For the class of Hamiltonians given by (15) and dynamics given by the exact solutions to (16), the period length

[TABLE]

is related to the initial value of the Hamiltonian $E=H(x_{0},p_{0})$ through the expression

[TABLE]

where $c_{\alpha,\beta}<\infty$ does not depend on $E$. In particular, if $\alpha=1+\gamma$ and $\beta=1+\gamma^{-1}$ for some $\gamma\in(0,\infty$), then $\nabla K\circ\nabla U(x)=x$ and $\eta=0$, meaning the period length does not depend on $E$.

The above result characterizes precisely how the speed of the Hamiltonian flow depends on the growth rates of the potential and kinetic energies. If $\eta<0$ the flow will be very fast at higher energies, and so explicit numerical methods will typically become unstable when $|x|$ is large. If $\eta>0$ then the period length will increase with $|x|$, so numerical methods will be more stable, but the sampler may take some time to move from the tails back into the centre of the space.

The ideal choice in terms of finding the right balance between efficiency and numerical stability would therefore appear to be setting $\eta=0$, which for this class of models can be achieved by choosing the kinetic energy as the $L_{p}$ dual of the potential. This corresponds to setting $\alpha=1+1/(\beta-1)$, which produces an algorithm for which the flow behaves in an auto-regressive manner when $|x|$ is large, since the composite gradient will be similar in magnitude to the current point. In this case tuning the method should also be more straightforward, as a choice of integrator step-size and integration time made in the transient phase should still be reasonable at equilibrium. When $\eta<0$, for any fixed choice of step-size $\varepsilon>0$ there is an $M_{\varepsilon}<\infty$ such that when $|x|>M_{\varepsilon}$ the algorithm will become numerically unstable, and when $\eta>0$ then it would be desirable in the transient phase to choose a comparatively large value for $\varepsilon$ and take more leapfrog steps, which may result in many proposals being rejected and unnecessary computational expense when the chain reaches equilibrium.

4 The implicit noise

There are likely to be regions in the centre of the space in which $\|\nabla U(x_{0})\|$ is small, meaning the increment $\varepsilon\nabla K\{p_{0}-\varepsilon\nabla U(x_{0})/2\}\approx\varepsilon\nabla K(p_{0})$. This can occur when the density $\pi(x)$ contains flat regions, and we show one such case in section 6. In these regions Hamiltonian Monte Carlo resembles a random walk Metropolis with proposal $x_{L\varepsilon}\approx x_{0}+L\varepsilon\nabla K(p_{0})$. It would seem sensible to choose $K(p)$ such that the distribution of $\nabla K(p)$ would be an appropriate choice in a random walk Metropolis here.

For the standard quadratic choice $\nabla K(p)=p$, which is hence also Gaussian. Different choices of kinetic energy, which may seem sensible by analysing the composite gradient, can however result in very different distributions for the implicit noise, which can either become bi-modal or exhibit very light tails. We plot histograms of $\nabla K(p)$ for some natural choices of momentum distribution in figure 1.

The optimum choice of implicit noise will typically be problem dependent, but it would seem logical to follow the standard choices used in the random walk Metropolis. Bi-modal distributions have recently been suggested for this method (Yang & Rodríguez, 2013), but we found that something which resembles a Gaussian distribution for small values of $p$ performed favourably in experiments.

5 General guidelines and the efficiency-robustness trade off

The above results allow us to offer some general strategies and observations for choosing a kinetic energy, dependent on the objectives of the user.

If we deem the goal to be to maximize the speed of the flow when $\|x\|$ is large whilst retaining numerical stability, then an appropriate choice of kinetic energy for a given potential $U$ is one for which when $\|\nabla U(x)\|$ is large $\nabla K\{p_{0}-\varepsilon\nabla U(x_{0})/2\}\approx-c_{\varepsilon}x$ for some $c_{\varepsilon}>0$. This can be achieved when $\nabla U(x)$ is invertible through the choice

[TABLE]

which implies that $\nabla K\circ\nabla U(x)\approx\nabla U^{-1}\circ\nabla U(x)=x$. Of course often $\nabla U(x)$ will not possess an inverse, in which case the considering the leading order term in $\|\nabla U(x)\|$ is of greatest importance. Below we give some examples of how such a choice for $K(p)$ can be constructed.

Example 2**.**

If $x\sim N(0,\Sigma)$ then $\nabla U(x)=\Sigma^{-1}x$, meaning $-(\nabla U)^{-1}(-p)=\Sigma p$, implying the choice $p\sim N(0,\Sigma^{-1})$ and corresponding kinetic energy $K(p)=p^{T}\Sigma p/2$.

Remark 2**.**

This result is intuitive, as if $K(p)=p^{T}M^{-1}p/2$, then the proposal after a single leapfrog step will be $x_{\varepsilon}=x_{0}-\varepsilon^{2}M^{-1}\nabla U(x)/2+\varepsilon M^{-1/2}z$ with $z\sim N(0,I)$. This is equivalent to a Metropolis-adjusted Langevin algorithm proposal with pre-conditioning matrix $M^{-1}$, and in that setting it is known that this should be made equal to $\Sigma$ for optimal sampling (Roberts & Rosenthal, 2001; Dalalyan, 2017).

Example 3**.**

Consider a distribution with density $\pi(x)\propto e^{-U(x)}$ where

[TABLE]

for some $C,\alpha,\beta>0$. Such double well potential models are important in many areas of physics, including quantum mechanical tunnelling (Liang & Müller-Kirsten, 1992). In this case $\nabla U(x)=4C\alpha x^{3}-2C\beta x$, which is typically not invertible, but for large $|x|$ we see that $\nabla U(x)/x^{3}\approx 4C\alpha$, which motivates the choice $\nabla K(p)=p^{1/3}/(4C\alpha)$ for large $|p|$.

Example 4**.**

Park & Casella (2008)** suggest a Bayesian version of the bridge estimator for regression of Frank & Friedman (1993). In this case

[TABLE]

for some $\lambda>0$, where $L(\beta)$ is a suitable loss function. If $q>1$ and the asymptotic growth of $L(\beta)$ is at most linear in $\|\beta\|$, then the choice $K(p)=\sum_{i}|p_{i}|^{\alpha}$ with $\alpha=1+1/(q-1)$ for large $\|p\|$ would maximize the speed of convergence during the transient phase of the algorithm. Logistic regression, the Huberized loss of Park & Casella (2008) and the asymmetric pinball loss (Yu & Moyeed, 2001) are possible choices for $L(\beta)$ which fit this criterion.

Example 5**.**

Neal (2003)** proposes the hierarchical model $x_{i}|\nu\sim N(0,e^{\nu})$ for $i=1,...,9$, with $\nu\sim N(0,3^{2})$, resulting in a joint distribution that resembles a ten-dimensional funnel. Here $U(x_{1},...,x_{9},\nu)=\nu^{2}/18+9\nu/2+e^{-\nu}\sum_{i}x_{i}^{2}/2$, meaning

[TABLE]

When $\nu\ll 0$ these gradients can be exponentially large, so a choice of $\nabla K(p)$ that is at most logarithmic in $p$ would be needed, which could arise from the choice of kinetic energy $K(p_{i})\approx\{1+\log(1+|p_{i}|)\}|p_{i}|$.

A general appreciation for the consequences of different kinetic energy choices can be gained by considering a formal Taylor series expansion

[TABLE]

In the Gaussian case $\nabla^{2}K(p)=I$, meaning the refreshed momentum $p_{0}$ interacts linearly with the composite gradient term regardless of the current position in the state space. If however, the momentum distribution is chosen to have heavier than Gaussian tails, then each element of $\nabla^{2}K(p)$ will become negligibly small when $\|\nabla U(x_{0})\|$ is large, reducing the influence of $p_{0}$. A sampler with a flow of this type will move towards a high-density region of the space in an almost deterministic fashion when $\|x\|$ is large. We formalize this intuition below.

Definition 3**.**

A Hamiltonian Monte Carlo sampler is called deterministic in the tails if as $\|x\|\to\infty$

[TABLE]

a.s., where $p\sim\nu(\cdot)$.

Given the above intuition, the following result is immediate.

Proposition 6**.**

A Hamiltonian Monte Carlo sampler will be deterministic in the tails if $\pi(\cdot)$ is light-tailed and $K(p)$ is twice differentiable with

[TABLE]

Remark 3**.**

Typically if a choice of momentum distribution with heavier than Gaussian tails is made then (19) will hold.

As alluded to in subsection 3.3, there is a trade off to be made between efficiency and step-size robustness. Choosing a lighter-tailed momentum distribution will result in a faster flow, and potentially a faster mixing Markov chain as a result, but the risk of numerical instabilities is higher, and so the algorithm will be very sensitive to the choice of $\varepsilon$. Choosing a heavier-tailed momentum distribution will result in slower flows, meaning algorithm performance will be less sensitive to step-size choice, but the best-case mixing time under the optimal step-size may be slower. We give a numerical example illustrating this in the supplementary material. If the tails are too heavy then the algorithm will also exhibit the negligible moves property of section 3.2. A kinetic energy choice with Laplacian tails, for example, would result in a sampler which may be comparatively slow in the tails but is guaranteed to be numerically stable, since in this case the composite gradient is $\text{sign}\{p_{0}-\varepsilon\nabla U(x_{0})/2\}$, which is bounded from above and below regardless of how fast $\|\nabla U(x_{0})\|$ grows.

Remark 4**.**

It has recently been shown that by choosing $\nu(\cdot)$ to be a Laplace distribution one can perform Hamiltonian Monte Carlo on problems for which the standard method is not applicable, such as problems with discontinuous likelihoods or involving discrete parameters. We do not comment further on this here, but see Nishimura et al. (2017) for details.

5.1 Relativistic power kinetic energies

From the above discussion it would seem sensible to design a kinetic energy function that looks similar to the Gaussian choice when its argument is small, but also allows a robust and efficient choice of composite gradient. The relativistic kinetic energy

[TABLE]

suggested in Lu et al. (2017)**, has many such desirable features, since it behaves as a quadratic function for small $p$ and is linear when $p$ is large. We slightly generalize this choice to allow for different tail behaviours, resulting in the family of relativistic power kinetic energies

[TABLE]

where $\beta\geq 1$ and each $\gamma_{i}>0$. When each $p_{i}^{2}\ll\gamma_{i}$ then

[TABLE]

and when each $p_{i}^{2}\gg\gamma_{i}$

[TABLE]

Drawing samples from this family is a requirement for implementation. In our experiments we used adaptive rejection sampling  (Gilks & Wild, 1992), which we found to be reasonably efficient. However, since the distribution is fairly regular in shape it is likely that specialized methods can be designed to be even more efficient. In addition, other strategies such a dependent sampling using Metropolis–Hastings could be employed, as in Akhmatskaya & Reich (2008)**.

5.2 Approximate gradients and doing away with Metropolis–Hastings

One example in which the robustness-efficiency trade off outlined above can be informative is when estimates are used in place of $\nabla U(x)$, which may be either intractable or very expensive to compute. Examples are given in Chen et al. (2014); Strathmann et al. (2015); Lindsten & Doucet (2016)**.

If we assume that at each iteration the approximate gradient $\tilde{\nabla}U(x)$ satisfies

[TABLE]

where $\eta_{x}$ follows a distribution that may depend on $x$, but such that $E(\eta_{x})=0$ for all $x$, then a similar Taylor expansion gives

[TABLE]

Therefore, if the momentum distribution is chosen to have heavy enough tails that for large $\|x\|$ the last term becomes negligibly small, then the effects of such an approximation are mitigated and the the resulting approximate composite gradient will closely resemble $\nabla K\left\{-\frac{\varepsilon}{2}\nabla U(x_{0})\right\}$.

In some of the above mentioned approximate implementations, the Metropolis step is also omitted to reduce computational costs. When the gradient term $\|\nabla U(x)\|$ grows at a faster than linear rate, then often the resulting Markov chains become transient when this is done, as shown for the unadjusted Langevin algorithm in Roberts & Tweedie (1996a)**. In such cases it is difficult to give long time guarantees on the level of bias induced from a finite simulation. Ensuring that the composite gradient grows no faster than linearly should mean that such transience is averted, meaning a stronger grasp of the degree of approximation can be established. We leave a detailed exploration of this for future work.

5.3 Heavy-tailed models

Based on the composite gradient intuition, it may seem desirable when $\pi(\cdot)$ is heavy-tailed to choose a kinetic energy for which $\|\nabla K(p)\|\to\infty$ as $\|p\|\to 0$. In this case when $\|x\|$ is large the potential gradient will be very small, so choosing the kinetic energy in this way one could still make $\nabla K\circ\nabla U(x)$ linear in $\|x\|$. However, the actual leapfrog perturbation is

[TABLE]

When $\pi(\cdot)$ is light-tailed the $p_{0}$ term becomes safe to ignore when $\|x\|$ is large, but if it is heavy-tailed then this is no longer the case as for large $\|x\|$ typical proposals will be of the form $x^{*}\approx x+L\varepsilon\nabla K(p)$, where $p\sim\nu(\cdot)$, which resembles a random walk. It is because of this that, as shown in Part (i) of Proposition 1, no choice of kinetic energy which is independent of $x$ can produce a geometrically ergodic chain.

It will likely be the case that choosing a heavier-tailed momentum distribution will be of benefit in this scenario in terms of rates of convergence, since this is true of the random walk Metropolis, as shown in Jarner & Roberts (2007)**. However, since it is the implicit noise $\nabla K(p)$ which drives the sampler, then care must also be taken with its form. In the case where $\nu(\cdot)$ is a Cauchy distribution, for example, then $\nabla K(p)=2p/(1+p^{2})$, which is lighter-tailed than $p$ itself, so here the benefits of choosing heavier tails are not as strong as those for the random walk algorithm.

6 Examples

6.1 Quantile regression

We consider a Bayesian quantile regression as introduced in Yu & Moyeed (2001). The goal is to estimate the $\tau$th quantile of a response $y\in\mathbb{R}$ conditioned on a collection of covariates $x\in\mathbb{R}^{m}$, written $\mu(x)=F^{-1}(\tau|x)$. Given $n$ data points $\{(x_{i},y_{i})\}_{i=1}^{n}$, we take the natural choice $\mu(x_{i},\beta)=\sum_{j=1}^{m}x_{ij}\beta_{j}$, and follow the approach of Fasiolo et al. (2017) by estimating a posterior distribution for $\beta$ as $\pi(\beta|y,x)\propto\exp\left\{-\sum_{i=1}^{n}L(\beta,x_{i},y_{i})\right\}\pi_{0}(\beta)$, using the general Bayesian updating framework of Bissiri et al. (2016), with $L(\beta,x_{i},y_{i})$ a smoothed version of the pinball loss introduced in Fasiolo et al. (2017)**, given by

[TABLE]

Here $g(\xi,\sigma,\tau)=\log\left\{\xi\sigma\text{Beta}\left[\xi(1-\tau),\xi\tau\right]\right\}$. As $\xi\to 0$ then the non-smooth pinball loss is recovered. This is linear in $\beta$, meaning that if $L_{q}$ priors are chosen as in bridge regression (Park & Casella, 2008) then the corresponding potential is

[TABLE]

for $1<q\leq 2$. The dominant term in $U(\beta)$ is the contribution from the prior, meaning that each $|\partial_{j}U(\beta)|=O(|\beta_{j}|^{q-1})$ for $j=1,...,m$.

We performed two studies on 20 simulated data points with $m=2$, fixing $\sigma=1$, $\lambda=1$ and $\xi=0.01$. In the first we set $q=2$ for each $\beta_{j}$, and in the second $q=1.5$. Four different momentum distributions were tested: an exponential power family with shape parameter 3, Gaussian, Laplacian, and a t distribution with 4 degrees of freedom. Although the set up is simple it still enables a demonstration of both the negligible moves property and how the composite gradient dictates sampler performance. The number of leapfrog steps was set to 1 and the step size was tuned based on achieving a 65-75% acceptance rate at equilibrium. The samplers were then initialised far from the region of high probability and convergence speed was assessed, and also initialised at equilibrium to assess mixing in favourable areas. Results are shown in Figure 2. In the Gaussian prior study the exponential power choice is not shown since the resulting sampler did not move, as explained by Proposition 1. As expected, the Gaussian choice, which results in a linearly growing composite gradient, reaches equilibrium quicker than the others. The Laplacian choice converges in a straight-line fashion, as the size of the proposed jump is always $\varepsilon\surd 2$ regardless of the current position. The Student’s t choice exhibits the negligible moves property as outlined in Proposition 3, with convergence visibly slow during the first 5000 iterations. When the heavier-tailed prior is chosen, the exponential power choice is now numerically stable and produces faster convergence than the Gaussian. Since the potential growth is $O(|\beta|^{1.5})$ then this option results in a linear composite gradient, so the performance is to be expected. The speed of the remaining choices is dictated by how the composite gradient grows, with the Student’s t distribution again performing the worst. In this study the difference is not so pronounced, as in probabilistic terms the sampler is initialised in a region that is not so far into the tails. The slower convergence overall compared to the Gaussian study reflects the fact that $\pi(\cdot)$ is heavier-tailed. In both studies all samplers mixed similarly well regardless of the kinetic energy when initialised at equilibrium, with slightly worse performance observed for the Student’s t distribution. As discussed in Section 4, the optimal choice of implicit noise is likely to be problem specific, with the Gaussian shape sensible when no other information is available.

6.2 The Ginzburg–Landau model

As a second example we take the model of phase transitions in condensed-matter physics proposed by Ginzburg and Landau (Goldenfeld, 1992, Chapter 5). We consider a three-dimensional $10^{3}$ lattice, where each site represents a random variable $\psi_{ijk}\in\mathbb{R}$. The corresponding potential for the 1000-dimensional parameter $\psi$ is

[TABLE]

where $\alpha,\lambda,\tau>0$ and $\tilde{\nabla}\psi_{ijk}=(\psi_{i_{+}jk}-\psi_{ijk},\psi_{ij_{+}k}-\psi_{ijk},\psi_{ijk_{+}}-\psi_{ijk})$, where $i_{+}=(i+1)\mod(10)$, and $j_{+},k_{+}$ are defined analogously. The modular structure imposes periodic boundary conditions on the lattice. When the parameter $\tau<1$, the conditional distribution for each $\psi_{ijk}$ looks Gaussian in the centre of the space but with lighter tails. The model exhibits a phase transition at $\tau=1$, so that when $\tau>1$ each conditional distribution is bi-modal. The parameter $\alpha$ controls the strength of correlation between neighbouring lattice sites, with larger values making the sampling problem more challenging. When $\tau<1$, each $\psi_{ijk}$ is likely to be close to zero, while in the bi-modal phase they are more likely to be non-negligible in magnitude. When $\alpha$ is large, neighbouring parameters are likely to have the same sign. The bi-modal phase therefore represents the system in its ordered state, whereas the system is disordered for $\tau<1$.

The inference problem is to estimate expectations with respect to the distribution with density proportional to $\exp\{-U(\psi)\}$. We generate samples using Hamiltonian Monte Carlo with four different choices of kinetic energy: the standard quadratic choice, the relativistic choice of Lu et al. (2017), the exponential power choice of Zhang et al. (2016) and the relativistic power choice introduced in section 5.1. In the latter two cases we set the shape parameter $\beta=4/3$, which results in a composite gradient that is linear when $\|\psi\|$ is large.

We perform two studies. In the first we initialize samplers in the centre of the space and compute the effective sample sizes from 10,000 iterations. In the second we initialize each $\psi_{ijk}\sim U[-10,10]$, and compute the number of iterations until $\max|\psi_{ijk}|\leq 2$. Results are shown in Table 1. The kinetic energies with Gaussian implicit noise produce larger effective sample sizes in this example. The reason is that there is a sharp drop in density when the quartic term begins to dominate the potential. As the exponential power choice results in a bi-modal implicit noise, it is more likely to propose larger jumps which fall into this low density region, and are hence rejected. To prevent this from happening too often a smaller step-size must be chosen than for the other methods. Among the other choices the Gaussian has the highest effective sample size, though this cannot necessarily be relied upon as the sampler is not geometrically ergodic. The slower speed for larger values of $\|\psi\|$ enforced in the relativistic and to a lesser degree relativistic power cases appears to slightly reduce efficiency here, which is sacrificed in favour of robustness. We simply set each tuning parameter $\gamma_{i}=1$ in our experiments. This could be tuned in order to improve performance, but we preferred to use default values to limit the extra tuning required for this choice. Sensitivity checks with $\gamma_{i}=2$ and $3$ did not yield appreciably different results. In terms of iterations to the centre, the Gaussian choice performs poorly as expected. The remaining samplers perform similarly. The relativistic power choice outperforms the exponential power choice because a smaller step-size was needed in the latter case for adequate mixing at equilibrium. The relativistic choice, though slower, is not critically so, because $10$ is not very large in absolute terms, meaning the differences between a drift which is $O(1)$ and $O(\|\psi\|)$, while visible, is not substantial here.

7 Discussion

We have described how changing the kinetic energy in Hamiltonian Monte Carlo affects performance. In practice, several other strategies could be considered, such as mixtures of kinetic energies, or employing a delayed rejection approach as in Mira et al. (2001)**. In addition, other stochastic processes that utilise Hamiltonian dynamics, such as underdamped Langevin diffusions, can also be considered (Stoltz & Trstanova, 2016).

In this work kinetic energies of the form $K(p)=\sum_{i}K_{i}(p_{i})$ were used. The reason for this is that $\partial_{i}K(p)=\partial_{i}K_{i}(p_{i})$, meaning that the when $\pi(x)$ is of the product form the Hamiltonian flows become independent. If a different choice were made, for example $K(p)=\beta^{-1}(1+\gamma^{-1}p^{T}p)^{\beta/2}$, then $\partial_{i}K(p)=(1+\gamma^{-1}p^{T}p)^{\beta/2-1}p_{i}$, so this does not always happen, here only when $\beta=2$. It is unclear which of these approaches is preferable from a theoretical perspective, but evaluating $p^{T}p$ is certainly a computational burden, and previous empirical evidence favours our approach (Lu et al., 2017).

Outside the realm of separable Hamiltonians, in Girolami & Calderhead (2011)* the choice $K(x,p)=\{\log|G(x)|+p^{T}G^{-1}(x)p\}/2$ is advised. Since an implicit integrator is required, the composite gradient intuition no longer directly applies. Here instead*

[TABLE]

where $G(x)$ is some Riemannian metric and $2f(x,p)=\log|G(x)|+p^{T}G^{-1}(x)p$. Key drivers of the dynamics are the terms $G^{-1}(x_{0})\nabla U(x_{0})$ and $G^{-1}(x_{\varepsilon})\nabla U(x_{0})$, which are often called the natural gradient (Amari & Nagaoka, 2007). Given recent results (Taylor, 2015), this approach may be advantageous when $\pi(\cdot)$ has heavy-tails. The log-determinant term can also be beneficial (Betancourt & Girolami, 2015).

Another line of further study is to assess how kinetic energy choice affects dynamic implementations such as the No-U-Turn Sampler (Hoffman & Gelman, 2014). Results in simple cases suggest that such implementations are favourable when $\pi(\cdot)$ is heavy-tailed (Livingstone et al., 2016).

Acknowledgement

We thank the reviewers for helping to improve the paper. SL thanks Michael Betancourt for sharing a calculation related to Proposition 5, Matteo Fasiolo for advice on quantile regression, and Gabriel Stoltz, Zofia Trstanova and Paul Fearnhead for useful discussions. The authors acknowledge support from the Engineering and Physical Sciences Research Council for support through grants EP/K014463/1, i-like (SL and GOR) and EP/D002060/1, CRiSM (GOR).

Appendix A Proofs of propositions

Proof of proposition 1.

(i) In Lemma B1 we show that $\int e^{\gamma\|x\|-U(x)}dx=\infty$ for any $\gamma>0$, and in Lemma B2 that for any $\eta>0$ there is an $r<\infty$ such that $P\{x,B_{r}(x)\}>1-\eta$. Theorem 2.2 of Jarner & Tweedie (2003) establishes that if these two conditions hold then the resulting Markov chain cannot be geometrically ergodic.

(ii) Lemmas B4, B5 and B6 show that when (7)-(12) hold, with probability one $\lim_{\|x_{0}\|\to\infty}\triangle(x_{0},p_{0})=\infty$, where $\triangle(x_{0,}p_{0})=(\|x_{L\varepsilon}\|+\|p_{L\varepsilon}\|)-(\|x_{0}\|+\|p_{0}\|)$. Under (9) this implies that with probability one $\lim_{\|x_{0}\|\to\infty}\triangle H(x_{0},p_{0})=\infty$, where $\triangle H(x_{0},p_{0})=H(x_{L\varepsilon},p_{L\varepsilon})-H(x_{0},p_{0})$. This in turn implies that with probability one $\lim_{\|x_{0}\|\to\infty}\alpha(x_{0},x_{L\varepsilon})=0$, which, using Proposition 5.1 of Roberts & Tweedie (1996b), establishes the result. ∎

Proof of proposition 2.

For the first part, note that the assumptions imply

[TABLE]

which implies $\limsup_{\|x\|\to\infty}\|\nabla K\circ\nabla U(x)\|/\|x\|<\infty$ as required. We prove the second part by induction. Precisely, we show that assuming $\|p_{i\varepsilon}\|\leq E_{i}\|x_{i\varepsilon}\|^{q}+F_{i}$ for some $E_{i},F_{i}<\infty$ implies $\|p_{(i+1)\varepsilon}\|\leq E_{i+1}\|x_{i\varepsilon}\|^{q}+F_{i+1}$ and $\|x_{(i+1)\varepsilon}\|\leq G_{i}\|x_{i\varepsilon}\|+H_{i}$ for $E_{i+1},F_{i+1},G_{i},H_{i}<\infty$. These in turn imply the result. First note that

[TABLE]

Using $\|\nabla U(x_{i\varepsilon})\|\leq A\|x_{i\varepsilon}\|^{q}+B$ gives

[TABLE]

Given this we can choose $G_{i}=\varepsilon C(E_{i}+\varepsilon A/2+\varepsilon B/2)^{1/q}+1$ and $H_{i}=\varepsilon C(E_{i}+\varepsilon A/2+\varepsilon B/2)^{1/q}+\varepsilon D$ to see that

[TABLE]

Iterating gives

[TABLE]

where $\mathbf{G}_{L}=G_{L-1}G_{L-2}...G_{0}$ and $\mathbf{H}_{L}=H_{L-1}+G_{L-1}H_{L-2}+G_{L-1}G_{L-2}H_{L-3}+...+G_{L-1}...G_{1}H_{0}$. Next recall that

[TABLE]

Combining with the assumption that $\|p_{i\varepsilon}\|\leq E_{i}\|x_{i\varepsilon}\|^{q}+F_{i}$, gives

[TABLE]

Setting $E_{i+1}=\left[E_{i}+\varepsilon A\left\{1+(G_{i}+H_{i})^{q}/2\right\}\right]$ and $F_{i+1}=\varepsilon\left\{A(G_{i}+H_{i})^{q}+2B\right\}/2+F_{i}$ then gives $\|p_{(i+1)\varepsilon}\|\leq E_{i+1}\|x_{i\varepsilon}\|^{q}+F_{i+1}$. Iterating then gives $\|p_{L\varepsilon}\|\leq E_{L}\|x_{0}\|^{q}+F_{L}$. Recalling that $\|p_{0}\|\leq E_{0}\|x_{0}\|^{q}+F_{0}$ by assumption completes the proof. ∎

Proof of proposition 3.

Consider the event $B=\{4\|p_{0}\|\leq\varepsilon\|\nabla U(x_{0})\|\}$, and note that $\lim_{\|x\|\to\infty}\text{pr}(B)=1$. We use the facts that $\|x_{L\varepsilon}-x_{0}\|\leq\sum_{i=1}^{L-1}\|x_{(i+1)\varepsilon}-x_{i\varepsilon}\|$, and that for any $i\in\{0,...,L-1\}$

[TABLE]

Taking $i=0$ gives

[TABLE]

Since $4\|p_{0}-\varepsilon\nabla U(x_{0})/2\|\geq\varepsilon\|\nabla U(x_{0})\|$ under $B$, it follows from the fact that $\pi(\cdot)$ is light-tailed and $\nu(\cdot)$ heavy-tailed that for every $\delta>0$ there is an $M<\infty$ such that whenever $\|x_{0}\|>M$ then $\|\nabla K(p_{\varepsilon/2})\|<\delta/\varepsilon$. Thus $\|x_{\varepsilon}-x_{0}\|$ can be made arbitrarily small by choosing an $x_{0}$ with large enough norm.

Recall that $\nabla U(x)$ is continuous by assumption. It follows from the preceding argument that for any $\gamma_{1}>0$ we can choose an $x_{0}$ with large enough norm that $\|\nabla U(x_{\varepsilon})-\nabla U(x_{0})\|<\gamma_{1}$ under $B$.

To complete the proof we show that if $\sum_{j=1}^{i}\|x_{j\varepsilon}-x_{(j-1)\varepsilon}\|<\delta/2$ then $\|x_{(i+1)\varepsilon}-x_{i\varepsilon}\|\leq\delta/2$ under $B$. Combining this with the previous paragraphs establishes that for any $\delta>0$ then there is an $x_{0}$ with large enough norm that $\|x_{L\varepsilon}-x_{0}\|<\delta$ if event $B$ holds, establishing the result.

From equation (21) the key factor in controlling $\|x_{(i+1)\varepsilon}-x_{i\varepsilon}\|$ is $\|p_{(i+1/2)\varepsilon}\|$, which can be lower bounded using

[TABLE]

If for any $\delta>0$ we can choose an $x_{0}$ with large enough norm that $\sum_{j=1}^{i}\|x_{j\varepsilon}-x_{(j-1)\varepsilon}\|<\delta/2$ then $\sum_{j=1}^{i}\|\nabla U(x_{j\varepsilon})-\nabla U(x_{0})\|$ can be made arbitrarily small through the same continuity argument made above. Thus, under $B$ it holds that $\|p_{(i+1/2)\varepsilon}\|\geq i\varepsilon\|\nabla U(x_{0})\|$, from which it follows that $\|x_{(i+1)\varepsilon}-x_{i\varepsilon}\|$ can be made arbitrarily small by choosing $\|x_{0}\|$ large enough. ∎

Proof of proposition 4.

It is shown in chapter 16 of Meyn & Tweedie (1993) that a geometric convergence bound is equivalent to the drift condition $\int V(y)P(x,dy)\leq\lambda V(x)$ whenever $x$ is outside some small set $C$, where $\lambda<1$. Lemma B7 establishes that if (14) holds then any small set must be bounded. Hence if a geometric bound holds here then

[TABLE]

For any $\delta>0$ we can write

[TABLE]

where $\epsilon=P\{x,B_{\delta}^{c}(x)\}$. If (i) holds then we can choose a $\delta<\delta^{\prime}$, so that

[TABLE]

Noting that both $\epsilon$ and $\epsilon^{\prime}$ can be made arbitrarily small as $\|x\|\to\infty$, this expression tends to $1$ in the same limit, proving the result. If (ii) holds, note that $\liminf_{\|x\|\to\infty}V(x)e^{-s\|x\|}=c$ implies that $\forall\epsilon^{\prime}>0$ there is an $M<\infty$ such that $V(x)e^{-s\|x\|}\geq c-\epsilon^{\prime}$ whenever $\|x\|\geq M$. This means that when $\|x\|>M$

[TABLE]

Condition (ii) also implies that for all $\epsilon^{\prime}>0$, there is a sequence $\{x_{i}\}_{i\geq 1}$ for which $\|x_{i}\|\to\infty$ as $i\to\infty$ such that whenever $i\geq N$ for some $N<\infty$ then $\|x_{i}\|>M$ and the condition $V(x_{i})e^{-s\|x_{i}\|}\leq c+\epsilon^{\prime}$ holds. Combining gives that for all $i\geq N$

[TABLE]

Since $\epsilon,\epsilon^{\prime}$ and $\delta$ can all be made arbitrarily small and $V(x_{i})\to\infty$ as $\|x_{i}\|\to\infty$, then this proves the result. ∎

Proof of proposition 5.

Assume $H(x_{0},p_{0})=E$ and $x_{0}=0$, $p_{0}=\left(\beta E\right)^{\frac{1}{\beta}}$. Take $4T$ to be the period length, and note that by the symmetry of the Hamiltonian in question this implies that $p_{T}=0$ and $x_{T}=\left(\alpha E\right)^{\frac{1}{\alpha}}$. Then

[TABLE]

Setting $b=(1-\beta)/\beta$, $c_{\beta}=\beta^{b}$ and noting that $p_{t}^{1-\beta}=c_{\beta}(E-\alpha^{-1}x_{t}^{\alpha})^{b}$ for $t\in[0,T]$, then the expression can be written

[TABLE]

Applying the change of variables $y_{t}=(\alpha E)^{-1/\alpha}x_{t}$ and setting $c_{\alpha}=\alpha^{1/\alpha}$ gives

[TABLE]

where . Now, we have that $\mathcal{P}(E)=f(E^{\eta})$, for some function $f$, where

[TABLE]

Setting $\alpha=1+\gamma$ and $\beta=1+\gamma^{-1}$ for some $\gamma>0$ gives

[TABLE]

as required. ∎

Proof of proposition 6.

Set $\gamma(x)=\min\left[\frac{\varepsilon}{4}\|\nabla U(x)\|,\left\|\nabla^{2}K\left\{\frac{\varepsilon}{4}\nabla U(x)\right\}\right\|^{-1/2}\right]$, and note that $\lim_{\|x\|\to\infty}\text{pr}_{\nu}\left\{\|p\|\leq\gamma(x)\right\}=1$. For $\|p\|\leq\gamma(x)$, as a direct consequence of the mean value inequality (Dieudonné, 1961)

[TABLE]

where $M(x)=\sup_{\{4\|p\|\geq\varepsilon\|\nabla U(x)\|\}}\|\nabla^{2}K(p)\|$. As the right-hand side tends to [math] as $\|x\|\to\infty$, then the result follows. ∎

Appendix B Technical results

Lemma 1**.**

If $\pi(\cdot)$ is heavy-tailed then for every $\gamma>0$

[TABLE]

Proof.

Choose $\delta<\gamma$. Let $B$ be a Euclidean ball centred at the origin such that $\|\nabla U(x)\|\leq\delta$ whenever $x\not\in B$. By continuity of $U(x)$, there is an $M<\infty$ such that $U(x)\leq M$ for all $x\in\partial B$. Then for all $x\not\in B$ the integrand is bounded below by $e^{(\gamma-\delta)\|x\|-M}$, which diverges uniformly and hence is not integrable. ∎

Lemma 2**.**

If $\pi(x)$ is heavy-tailed then for any $\eta>0$ there is an $r<\infty$ such that

[TABLE]

Proof.

We need to show that $Q\{x,B_{r}(x)\}>1-\eta$, for any $x$. After one leapfrog step we have

[TABLE]

Write $\|x\|_{\infty}$ for the supremum norm, and note that by equivalence of norms in finite dimensions we can write $\|x\|_{\infty}\leq C\|x\|$ for all $x$, for some $C<\infty$. We have that $\nabla U(x)\in C_{0}(\mathbb{R}^{d})$, which implies $\|\nabla U(x)\|<M/C$ for some $M<\infty$ which does not depend on $x$, so that $\|\nabla U(x)\|_{\infty}<M$. The class of distributions for $\{p_{0}-\varepsilon\nabla U(x_{0})/2\}$ is therefore tight. Now recall that if $f$ is a locally bounded function, and $\mathcal{F}$ a tight family of probability measures, then the resulting family of probability measures induced by pushing forward each element of $\mathcal{F}$ through $f$ is also tight. So since $\nabla K$ is continuous and hence locally bounded, the result follows. ∎

Lemma 3**.**

If (7) and (10) hold then

[TABLE]

Proof.

First we re-write the expression

[TABLE]

Now, (10) implies that the first term will be bounded below by a finite positive constant, while (7) ensures that the second will have an infinite limit, proving the result. ∎

Lemma 4**.**

If $\pi(\cdot)$ is light-tailed, (7) and either of (11) or (12) hold and $\|p_{0}\|\leq\frac{\varepsilon}{4}\min\{\|\nabla U(x_{0})\|,\|\nabla U(x_{0})\|_{\infty}\}$ then there is a $\gamma_{M}<\infty$ such that, provided $\|x_{0}\|\geq\gamma_{M}$, it holds that $\|x_{\varepsilon}\|\geq M\|x_{0}\|$, for any $M<\infty$.

Proof.

Note

[TABLE]

It is therefore sufficient to show that for any $M<\infty$ we can choose an $\|x_{0}\|$ large enough that

[TABLE]

Under (11), note that

[TABLE]

which implies

[TABLE]

By (24), therefore, if $\|x_{0}\|$ is chosen to be large enough then this can be made $\geq(M+1)\|x_{0}\|/\varepsilon$, for any finite $M$, proving the result.

Under (12), recall that there exists global constants $C,c>0$ such that $C\|\nabla U(x)\|\geq\|\nabla U(x)\|_{\infty}\geq c\|\nabla U(x)\|$ for all $x\in\mathbb{R}^{d}$. It suffices in this setting therefore to show that we can choose an $\|x_{0}\|$ large enough that

[TABLE]

We have

[TABLE]

Write $i^{*}$ and $j^{*}$ to denote the indices for which $\|p_{0}-\nabla U(x_{0})\|_{\infty}=|p_{0}(i^{*})-\partial_{i^{*}}U(x_{0})|$ and $\|\nabla U(x_{0})\|_{\infty}=|\partial_{j^{*}}U(x_{0})|$. We have:

[TABLE]

Now, $|p_{0}(i^{*})-(\varepsilon/2)\partial_{i^{*}}U(x_{0})|\geq(\varepsilon/4)|\partial_{j^{*}}U(x_{0})|$ implies that $|k^{\prime}\{p_{0}(i^{*})-(\varepsilon/2)\partial_{i^{*}}U(x_{0})\}|\geq|k^{\prime}\{(\varepsilon/4)\partial_{j^{*}}U(x_{0})\}|=\|\nabla K\{(\varepsilon/4)\nabla U(x_{0})\}\|_{\infty}$. Using the global bounds then we see that for any $M<\infty$ we can choose an $\|x_{0}\|$ large enough that

[TABLE]

establishing the result. ∎

Lemma 5**.**

*If $\pi(\cdot)$ is light-tailed and (7)-(10) and one of (11) and (12) hold, and provided that for any fixed $i\geq 0$

(i) $\|x_{0}\|\geq\gamma_{M}$ for some $\gamma_{M}<\infty$,

(ii) $\|p_{0}\|\leq(\varepsilon/4)\min\{\|\nabla U(x_{0})\|,\|\nabla U(x_{0})\|_{\infty}\}$,

(iii) $M$ is large enough that $\phi(M)\geq 7/3$ with $\phi$ as in (8),

(iv) $\|x_{j\varepsilon}\|\geq M\|x_{(j-1)\varepsilon}\|$ for all $j\leq i$,

it holds for any finite $M<\infty$ that $\|x_{(i+1)\varepsilon}\|\geq M\|x_{i\varepsilon}\|$.*

Proof.

We first show

[TABLE]

To show (25), first note by iterating (4) and noting $p_{(i+1/2)\varepsilon}=p_{i\varepsilon}-\varepsilon\nabla U(x_{i\varepsilon})/2$ that

[TABLE]

Using the stated assumption that $\|p_{0}\|\leq(\varepsilon/4)\|\nabla U(x_{0})\|$ then gives

[TABLE]

Now, (8) implies that if $\|x_{(j-1)\varepsilon}\|\leq M^{-1}\|x_{j\varepsilon}\|$, $\|\nabla U(x_{(j-1)\varepsilon})\|\leq\phi(M)^{-1}\|\nabla U(x_{j\varepsilon})\|$. Substituting into the above expression gives

[TABLE]

To finish the argument we need to therefore show that $\{1-\sum_{j=0}^{i-1}\phi(M)^{j-i}\}\geq 1/4$. First note that

[TABLE]

Using simple geometric series identities gives

[TABLE]

Hence $\{1-\sum_{j=0}^{i-1}\phi(M)^{j-i}\}\geq 1/4$ if

[TABLE]

which will hold if $\phi(M)\geq 7/3$ for any $i\geq 1$. The stated assumption that $\phi(M)\geq 7/3$ therefore establishes (25).

If (11) holds, then (25) can be combined with (24) directly to show that for any $M<\infty$ and $\varepsilon\in(0,\infty)$, if $\|x_{i\varepsilon}\|\geq\gamma_{M}$ for some suitably chosen $\gamma_{M}<\infty$ it will hold that

[TABLE]

If (12) holds instead of (11), then the same result can be established using a similar argument to that given in the last paragraph of the proof of Lemma B4. Hence, provided that $\|x_{0}\|\geq\gamma_{M}$, for all $i\geq 0$ it holds that

[TABLE]

which proves the result.

∎

Lemma 6**.**

Under the same conditions as Lemma B5 and provided $M$ is such that $\phi(M)\geq 5$, then $\|p_{L\varepsilon}\|\geq\|p_{0}\|$ for any $L\geq 1$.

Proof.

We have

[TABLE]

by recalling that as in Lemma B5 $\|p_{0}\|\leq(\varepsilon/4)\|\nabla U(x_{0})\|$. Using (8) and the stated assumptions we have for any $i\leq L$ that

[TABLE]

which implies

[TABLE]

The stated assumptions $\|p_{0}\|\leq(\varepsilon/4)\|\nabla U(x_{0})\|$ and $\phi(M)^{L}\|\nabla U(x_{0})\|\leq\|\nabla U(x_{L\varepsilon})\|$ lead to the bound

[TABLE]

which, when combined with the above inequality, give

[TABLE]

As shown in the proof of Lemma B5 $\sum_{i=0}^{L-1}\phi(M)^{i-L}=\{1-\phi(M)^{-L}\}/\{\phi(M)-1\}$. Hence the result is proven if

[TABLE]

which will indeed be true under the stated assumption that $\phi(M)>5$. ∎

Lemma 7**.**

If (14) holds for a Hamiltonian Monte Carlo method then any small set must be bounded.

Proof.

Since $\nabla K(p)$ and $\nabla U(x)$ are continuous and $\nabla K\circ\nabla U(x)$ is vanishing at infinity, then $\nabla K\circ\nabla U(x)$ is also bounded, implying that the collection of Hamiltonian Monte Carlo increments $\{P(x,\cdot)-x\}$ is uniformly tight, using a similar argument to that employed in the proof of Lemma B2. Hence, any small set must be bounded, following Lemma 2.2 of Jarner & Hansen (2000). ∎

Appendix C Numerical example of the efficiency-robustness trade-off

We consider the distribution $\pi(x)\propto\exp\left(-\beta^{-1}|x|^{\beta}\right)$ with $\beta=1.5$, and compare the quadratic and relativistic kinetic energy choices. We test 80 evenly spaced choices of step-size $\varepsilon$ spanning from 0.1 to 5, and in each case begin the sampler at equilibrium and compute the empirical expected squared jump distance from a chain of length 200,000, with the number of leapfrog steps randomly selected uniformly between 1 and 5 at each iteration. The results are shown in Figure 3.

As can be seen, the quadratic choice leads to a higher optimal value, but when step-sizes are chosen to be too large the jump distance drops quickly. The relativistic choice, by comparison, exhibits a larger degree of robustness to bigger step-sizes.