2.1. RWM

The RWM algorithm (Algorithm 1) is the most basic commonly used MCMC algorithm. At each step i of the Markov chain, the RWM algorithm proposes a new value $\hat{x}_{i+1}$ obtained by adding a small random increment to the current position $x_i$ . The proposal is accepted with a probability of $\mathrm{min}\{{\pi(\hat{x}_{i+1})}/{\pi(x_{i})}, 1\}$ , in which case we set $x_{i+1} = \hat{x}_{i+1}$ . Otherwise, we say that the step is rejected and we set $x_{i+1} = x_{i}$ , so that the algorithm stays at its current position until the next time step.

Note that we describe here the simple but important special case of Gaussian increments; see, e.g., [Reference Robert and Casella49] for a more general survey.

Algorithm 1

Random walk metropolis [Reference Metropolis, Rosenbluth, Rosenbluth, Teller and Teller45]

2.2. HMC algorithms

HMC algorithms [Reference Duane, Kennedy, Pendleton and Roweth20] seek to avoid quadratic slowdowns associated with diffusive ‘random walk’ behavior. They do so by adding a notion of ‘momentum’ to the random walk, which encourages the underlying Markov chain to propose longer steps without incurring a large chance of rejection. For this reason HMC algorithms work especially well in high dimensions: concentration of the posterior measure $\pi$ causes most other MCMC algorithms to either propose steps that are very small or larger steps that are rejected with high probability (see, e.g., [Reference Beskos, Pillai, Roberts, Sanz-Serna and Stuart4, Reference Roberts and Rosenthal52] for discussion of these heuristics).

In this section we review two commonly used HMC algorithms. Both algorithms are based on Hamiltonian dynamics, which we review here. Fix a $C^{2}$ function $H\, :\, \mathbb{R}^{2d} \mapsto \mathbb{R}^{+}$ , which we call the Hamiltonian function, and starting points $p_{0},q_{0} \in \mathbb{R}^{d}$ . We then define Hamilton’s equations to be the pair of differential equations

(2.1) \begin{equation} \frac{\mathrm{d}}{\mathrm{d} t} p(t) = - \frac{\partial H(p,q)}{\partial q}, \quad \frac{\mathrm{d}}{\mathrm{d} t} q(t) = \frac{\partial H(p,q)}{\partial p},\end{equation}

with initial condition $q(0) = q_{0}$ , $p(0) = p_{0}$ . Throughout the paper, we denote by $\gamma_{p,q}(t) = q(t)$ the first component of the solution to these equations with these initial conditions, so that, e.g., $\gamma_{p,q}'(t) = p(t)$ for the ‘standard’ choice of Hamiltonian in (2.4). Throughout the paper, we consider only Hamiltonians H for which the solution $\{ \gamma_{p,q}(t)\}$ exists for all $p,q \in \mathbb{R}^{d}$ and all $t \in \mathbb{R}^{+}$ ; see Appendix C for some simple sufficient conditions under which this holds, and which are satisfied by all of our examples. Our notation elides the dependence on the associated Hamiltonian function H, which will always be clear from the context.

The solutions to Hamilton’s equations have a number of special properties. Two of the most important are conservation of energy and conservation of volume.

1. 1. Conservation of energy. For $\gamma_{p,q}$ as above, we have

   (2.2) \begin{equation} \frac{\mathrm{d}}{\mathrm{d} t} H(\gamma_{p,q}'(t), \gamma_{p,q}(t)) \equiv 0.\end{equation}

2. 2. Conservation of volume. Define $\mu_{H} = {\mathrm{e}^{-H(p,q)}}/{\int_{p,q} \mathrm{e}^{-H(p,q)} \,\mathrm{d} q \,\mathrm{d} p}$ , which we call ‘Hamilton’s distribution’ (it has several names in the literature, including the more generic ‘Boltzmann–Gibbs’ distribution). If we sample $(P,Q) \sim \mu_{H}$ , then we also have

   (2.3) \begin{equation} \mathcal{L}(\gamma_{P,Q}'(t), \gamma_{P,Q}(t)) = \mu_{H}\end{equation}
   for all $t \in \mathbb{R}^{+}$ . This fact is a restatement of Liouville’s theorem in probabilistic language.

We refer the reader to [Reference Lanczos29] for proofs of these facts. The second fact motivates the standard choice of Hamiltonian

(2.4) \begin{equation} H(p,q) = -\log(\pi(q)) + \frac{1}{2} \| p \|^{2}.\end{equation}

Under this measure, if $Q \sim \pi$ and $P \sim \mathcal{N}(0,1)$ , then $\gamma_{P,Q}(t) \sim \pi$ for all $t \in \mathbb{R}$ . This choice leads to the HMC algorithm (Algorithm 2), developed in [Reference Mehlig, Heermann and Forrest44].

Algorithm 2

Isotropic-momentum HMC (idealized symplectic integrator) [Reference Mehlig, Heermann and Forrest44]

Although this is not obvious by inspection, Algorithm 2 is reversible with respect to $\mu_{H}$ .

Riemannian manifold HMC seeks to take longer steps by choosing initial momenta from a multivariate Gaussian distribution that agrees with the local geometry of the posterior density $\pi$ . This is achieved by using trajectories that evolve according to Hamiltonian dynamics on a Riemannian manifold, with metric defined by some symmetric and positive definite matrix G(q). For the special case of Algorithm 2, we have $G(q) = I_d$ . Alternatively, one may choose G(q) to be a regularization of the Hessian of $U(q) \equiv - \log(\pi(q))$ , which acts as a local pre-conditioner for U [Reference Betancourt5]. This can result in much faster mixing for distributions that are not close to isotropic; for example, it avoids the bad performance of Algorithm 2 on the nearly degenerate Gaussian that is observed in Theorem 3. The algorithm is identical to Algorithm 2, except for the choice of the Hamiltonian H.

Algorithm 3

Riemannian manifold HMC (idealized symplectic integrator) [Reference Girolami and Calderhead22, Reference Girolami, Calderhead and Chin23]

For both algorithms, we make the following assumption.

Assumption 1. Let $\lambda_{0}(q), \lambda_{1}(q)$ be the smallest and largest singular values of G(q). Assume that for any compact set $A \subset \mathbb{R}^{d}$ we have both

\begin{equation*}\sup_{q \in A} \max \left( \frac{1}{\lambda_{0}(q)}, \, \lambda_{1}(q) \right) < \infty\end{equation*}

and

\begin{equation*}\sup_{q \in A} \| \frac{ \partial H}{\partial q} \| < \infty.\end{equation*}

Assumption 1 says that the singular values of the metric G are bounded away from 0 and $\infty$ on compact sets, and also that the target distribution does not ‘blow up’ on compact sets. The former assumption is satisfied when G is continuous and positive definite, and in particular in the special case $G = Id$ used by Algorithm 2. The latter assumption is satisfied when H has continuous derivatives.

2.3. Cheeger’s inequality and the spectral gap

We recall the basic definitions used to measure the efficiency of MCMC algorithms. Let L be a reversible transition kernel with unique stationary distribution $\mu$ on $\mathbb{R}^{d}$ . We view L as an operator from $L_{2}(\pi)$ to itself:

\begin{equation*}(L f)(x) = \int_{y \in \mathbb{R}^{d}} L(x,\mathrm{d} y) f(y).\end{equation*}

The constant function is always an eigenfunction of this operator, with eigenvalue 1. We define the space $W^{\perp} = \{ f \in L_{2}(\mu)\, :\, \int_{x} f(x) \mu(\mathrm{d} x) = 0\}$ of functions that are orthogonal to the constant function, and denote by $L^{\perp}$ the restriction of the operator L to the space $W^{\perp}$ . We then define the spectral gap $\rho$ of L by the formula

(2.5) \begin{equation} \rho = \rho(L) \equiv 1 - \sup \{ \lambda \,:\, \lambda \in \mathrm{Spectrum}(L^{\perp}) \},\end{equation}

where Spectrum refers to the usual spectrum of an operator. Geometric ergodicity of HMC algorithms was proved under very general conditions in [Reference Bou-Rabee and Sanz-Serna13, Reference Livingstone, Betancourt, Byrne and Girolami34], implying existence of a non-zero spectral gap under those conditions [Reference Roberts and Rosenthal51]. We call ${1}/{\rho(L)}$ the relaxation time of L.

Cheeger’s inequality [Reference Cheeger15, Reference Lawler and Sokal30] provides bounds for the spectral gap in terms of the ability of L to move from any set to its complement in a single step. This ability is measured by the conductance $\Phi(L)$ , which is defined by the pair of equations

\begin{align*}\Phi(L) &= \inf_{S \in \mathcal{A}\, :\, 0 < \mu(S) < \frac{1}{2}} \Phi(L,S) \\\Phi(L,S) &= \frac{ \int_{x} \unicode{x1D7D9}\{x \in S\} L(x,S^{c}) \mu(\mathrm{d} x)}{\mu(S) },\end{align*}

where $\mathcal{A}= \mathcal{A}(\mathbb{R}^{d})$ denotes the usual collection of Lebesgue-measurable subsets of $\mathbb{R}^{d}$ . Note that the conductance has essentially the same behaviour as the ‘Cheeger constant’ and ‘bottleneck ratio’, which are also common in the MCMC literature.

One version of Cheeger’s inequality for Markov chains, first proved in [Reference Lawler and Sokal30], gives

(2.6) \begin{equation} \frac{\Phi(L)^2}{8} \leq \rho(L) \leq 2 \Phi(L).\end{equation}

We note that there are several small variants of Cheeger’s inequality and the conductance in the Markov chain literature, and all are essentially equivalent for our purposes. In fact, the original version of conductance from [Reference Lawler and Sokal30] is slightly different from ours: it defines a quantity k(L, S) that differs by a multiplicative factor of $1 \leq {1}/({1-\pi(S)}) \leq 2$ from our $\Phi(L,S)$ . We chose the definition $\Phi(L,S)$ because it seemed most common in the Markov chain literature (as evidenced by, e.g., the extremely popular textbook [Reference Levin, Peres and Wilmer33]), but it is straightforward to translate our results to other definitions of conductance.

3.1. Notation related to manifolds and intersections

In this section, we give some basic results about intersections of manifolds. As our target audience consists mostly of probabilists, we:

1. 1. spend extra time stating and illustrating the consequences of our definitions that will be most important, even though these illustrations are not required for the rest of the paper;

2. 2. give references for the existence of various objects under the assumptions made in this section.

Fix a set $S \subset \mathbb{R}^{d}$ whose topological boundary $\partial S$ is an embedded $(d-1)$ -dimensional smooth manifold (see, e.g., p. 17 of [Reference Guillemin and Pollack25] for the definition of an embedding; an embedded submanifold of $\mathbb{R}^{d}$ is one for which the inclusion map is an embedding). Fix also an integration time $T \in \mathbb{R}^{+}$ . In this special case, we define the (possibly infinite) number of intersections $N_{\partial S}$ as follows.

Definition 1. For fixed $T > 0$ , define the set $\mathcal{G}_{S} = \mathcal{G}_{S}(T) \subset \mathbb{R}^{2d}$ to be all pairs $(p,q) \in \mathbb{R}^{2d}$ satisfying:

1. 1. the function $\gamma_{p,q}\, :\, \mathbb{R} \mapsto \mathbb{R}^{d}$ is transverse to the manifold $\partial S$ ;

2. 2. the set $\{ t \in [0,T]\, :\, \gamma_{p,q}(t) \in \partial S\}$ is finite.

For $(p,q) \in \mathcal{G}_{S}$ , define

\begin{equation*}N_{\partial S} = N_{\partial S}(p,q) \equiv |\{ t \in [0,T]\, :\, \gamma_{p,q}(t) \in \partial S\}|.\end{equation*}

To recall a precise definition of a transverse intersection, see, e.g., p. 30 of [Reference Guillemin and Pollack25]. In this paper, we are primarily interested in the following property of transverse intersections. Let $\partial S \subset \mathbb{R}^{d}$ be a closed surface that partitions $\mathbb{R}^{d}$ into two open pieces $P_{1}'$ , $P_{2}'$ with closures $P_{1},P_{2}$ , and let $\gamma\, :\, [0,1] \mapsto \mathbb{R}^{d}$ be a curve that is transverse to $\partial S$ . Assume that $|\gamma([0,1]) \cap \partial S| < \infty$ . Then $\gamma(0) \notin \partial S, \, \gamma(1)$ are in the same closed piece of $\mathbb{R}^{d}$ if and only if $|\gamma([0,1]) \cap \partial S| < \infty$ is even.

In other words, if all intersections are transverse, we can tell which side of a curve a path will end on by simply counting the number of intersections. This is not true for non-transverse intersections. For example, consider the intersection of the circle given by $x^{2} + y^{2} = 1$ and the curve $\gamma(t)=(2t - 1, 1)$ ; they intersect at exactly one point but the curve never enters the inside of the circle. See also Figure 2.

Figure 2.

A boundary $\partial S$ is shown in black. The intersection of $\partial S$ with the top curve is transverse; the intersection with the show right-hand curve is not.

We next show that this definition of $N_{\partial S}(p,q)$ can easily be extended to all $(p,q) \in \mathbb{R}^{2d}$ , under modest assumptions about S. We need the following notation.

$\bullet$ For $q \in \partial S$ , define $\eta(q)$ to be the unit normal vector of $\partial S$ at q that points away from the interior of S, and for $p \in \mathbb{R}^{d}$ define

\begin{equation*}p_{q} = \langle G^{-1}(q) p, \eta(q) \rangle\end{equation*}

to be the component of $ G^{-1}(q) p$ in the direction orthogonal to $\partial S$ at q.

$\bullet$ For $q \in \partial S$ , define the set

\begin{equation*}\mathcal{P}^{+}_{S}(q) = \{ p \in \mathbb{R}^{d}\, :\, p_{q} > 0 \}\end{equation*}

to be the half-space of momentum vectors pointing away from S at $q \in \partial S$ .

For $x \in \partial S$ , denote by $\mathcal{T}_{x}$ the tangent space of $\partial S$ at x (see, e.g., p. 9 of [Reference Guillemin and Pollack25] for the definition of a tangent space in this context). We view this tangent space as being embedded in the same copy of $\mathbb{R}^{d}$ as S and being based at the point x (so that, for example, it will generally not include 0). Denote by $\mathbb{Proj}_{x}\, :\, \mathbb{R}^{d} \mapsto \mathcal{T}_{x}$ the usual orthogonal projection map from $\mathbb{R}^{d}$ to the tangent space. We assume.

Remark 2. We have not tried to provide the weakest-possible assumptions in Assumption 2. Since Lemma 2 implies that the conductance of most realistic HMC chains can be computed entirely in terms of sets that satisfy Assumption 2, these assumptions are weak enough for our purposes. Here, we explain why something like Assumption 2 is needed at all.

The goal is to rule out the possibility that solutions to Hamilton’s equations will pass through $\partial S$ very many times over very short intervals. To give a simple pathological example in $d=2$ dimensions, we wish to avoid sets such as

(3.1) \begin{equation} S_{\infty} = \bigcup_{n \in \mathbb{N}} [0,1] \times [(4n+1)^{-2}, (4n)^{-2}] \subset \mathbb{R}^{2}.\end{equation}

The problem with $S_{\infty}$ is that the boundary $\partial S$ contains a countably infinite union of parallel lines within the compact set $[0,1]^{2}$ , and so arbitrarily short Hamiltonian paths can cross $\partial S$ arbitrarily (or even infinitely) often.

On the other hand, ‘nice’ sets such as the usual Euclidean ball satisfy our assumption. In fact, we only need easily drawn sets, such as rectangles with rounded corners, in our proofs.

We next note that $N_{\partial S}$ can be defined $\mathbb{R}^{2d}$ -almost everywhere.

By this lemma, for all fixed $T \in \mathbb{R}^{+}$ we can define a measurable function $N_{\partial S}\, :\, \mathbb{R}^{2d} \mapsto \mathbb{N}$ that agrees with Definition 1 for almost every value of $(p,q) \in \mathbb{R}^{2d}$ . Throughout the rest of the paper, we always assume that all intersections are transverse. By this lemma, the complement of the set $\cap_{n \in \mathbb{N}} \mathcal{G}_{S}(n)$ has measure 0, and so this assumption will not influence the results of any calculations.

Next, we define a family of measures on $\mathbb{R}^{2d}$ obtained by ‘tilting’ $\mu_{H}$ according to the number $N_{\partial S}$ of crossings of $\partial S$ . Define

(3.2) \begin{equation} \Phi^+ = \frac{1}{2} T \cdot \int_{\partial S} \int_{\mathbb{R}^d}\mu_{H}(p,q) \cdot |p_q| (\mathrm{d}p\,\mathrm{d} _{\partial S}q).\end{equation}

When $0 < \Phi^{+} < \infty$ , we can make the following definition.

Definition 2. (Tilted measures.) Define $\mathbb{Q}$ to be the probability measure on $\mathbb{R}^{2d}$ with density

\begin{equation*}\mathbb{Q}(p,q) \propto \mathrm{e}^{-H(p,q)} \, N_{\partial S}(p,q).\end{equation*}

Although $\Phi^{+}$ can be infinite (and, thus, the tilted measure can fail to exist), we will see that it is always finite under the main assumptions of this paper.

3.2. Main conductance formula

Our main theorem is as follows.

Theorem 1. Let $T \in \mathbb{R}^{+}$ , let $\pi$ be any $C^{2}$ probability density (with respect to Lebesgue measure) on $\mathbb{R}^d$ , let K and H be the transition kernel and Hamiltonian of either Algorithm 2 or Algorithm 3 with these parameters, and let $\mu_{H}(p,q) \propto \mathrm{e}^{-H(p,q)}$ be the associated Hamiltonian measure. Assume that the associated solutions $\{ \gamma_{p,q}(t)\}$ to Hamilton’s equations exist for all $p, q \in \mathbb{R}^{d}$ and $t \in \mathbb{R}^{+}$ , and that the solutions are $C^{2}$ . Let $S \subset \mathbb{R}^{d}$ satisfy Assumption 2. Finally, in the case of Algorithm 3, let Assumption 1 also hold.

Then $0 < \Phi^{+} < \infty$ , and the conductance of K satisfies

(3.3) \begin{equation}\Phi(K,S) = \Phi^+ \cdot \mathrm{E}\bigg[\frac{1}{N_{\partial S}(P,Q)} \cdot \unicode{x1D7D9}\{N_{\partial S}(P,Q) \, \mathrm{odd}\}\bigg] \bigg/ \pi(S),\end{equation}

where the expectation is taken with respect to the random variables $(P,Q) \sim \mathbb{Q}$ and we use the convention that $({1}/{N_{\partial S}(P,Q)}) \cdot \unicode{x1D7D9}\{N_{\partial S}(P,Q) \, \mathrm{odd}\} = 0$ if $N_{\partial S}(P,Q)$ is almost-surely equal to 0. Furthermore, in the special case of Algorithm 2, $\Phi^{+}$ can be written as the following simpler integral:

(3.4) \begin{equation} \Phi^+ = \frac{T}{\sqrt{2 \pi}} \cdot \int_{\partial S} \pi(q) \, \mathrm{d}_{\partial S}q.\end{equation}

Recall that, in (3.2) and (3.4), the integral over $\partial S$ is taken with respect to the volume measure on the $(d-1)$ -dimensional manifold $\partial S$ , not with respect to the Lebesgue measure on the d-dimensional space $\mathbb{R}^{d}$ .

In many applications, it is more important to get a good bound on the conductance than to compute it directly. Observing that $\mathrm{E}_\mathbb{Q} [({1}/{N_{\partial S}}) \cdot \unicode{x1D7D9}\{N_{\partial S} \, \mathrm{odd}\}] \leq 1$ , we have the following useful consequence of Theorem 1.

Corollary 1. (Simple conductance bound.) Set notation as in Theorem 1. Then the conductance for Algorithm 2 is bounded by

\begin{equation*}\Phi(K,S) \leq \frac{T}{\sqrt{2 \pi}} \frac{\int_{\partial S} \pi(q) \, \mathrm{d}_{\partial S}q }{\pi(S)}.\end{equation*}

In simple examples, this bound lets us immediately check that the time it takes HMC chains to search for sub-Gaussian modes grows exponentially with both the dimension and the distance between modes. Numerical simulations for various two-mode densities that approximate Gaussian mixture models (Figure 3) suggest that this upper bound is nearly tight in many cases where T is not too large.

3.3. Removing Assumption 2

Theorem 1 only applies to sets that satisfy Assumption 2. Fortunately, it is not necessary to consider any other sets to compute the conductance of typical HMC Markov chains.

Lemma 2. (Sufficiency of compact locally good manifolds.) Let $\pi \in C^{2}$ , fix $0 < T < \infty$ , and let K be the transition kernel given by Algorithm 2 with these parameters. Let $\mathfrak{B}$ be the Lebesgue measurable subsets of $\mathbb{R}^{d}$ , and let $\mathfrak{A} \subset \mathfrak{B}$ be the subsets that also satisfy Assumption 2. Then

(3.5) \begin{equation} \inf_{S \in \mathfrak{A}} \Phi(K,S) = \inf_{S \in \mathfrak{B}} \Phi(K,S).\end{equation}

4.1. Application: HMC targeting mixture of Gaussians

For $\sigma > 0$ , define the mixture distribution

(4.1) \begin{equation} \pi_{\sigma} = \frac{1}{2} \mathcal{N}(-1,\sigma^2) + \frac{1}{2} \mathcal{N}(1,\sigma^2)\end{equation}

and denote its density by $f_{\sigma}$ . Let $K_{\sigma}$ be the transition kernel of Algorithm 2 with target distribution $\pi = \pi_{\sigma}$ and time-step $T = T_{\sigma} \equiv \sigma$ . Denote by $\lambda_{\sigma}$ the relaxation time of $K_{\sigma}$ , and denote by $\Phi_{\sigma} = \Phi(K_{\sigma}, (\!-\infty,0))$ the conductance associated with kernel $K_{\sigma}$ and set $(\!-\infty,0)$ .

For fixed $x \in (\!-\infty,0)$ , let $\{X_{t}^{(\sigma)}\}_{t \in \mathbb{N}}$ be a Markov chain with transition kernel $K_{\sigma}$ and initial point $X_{1}^{(\sigma)} = x$ . Define the hitting time

(4.2) \begin{equation} \tau_{x}^{(\sigma)} = \inf \{ t > 0\, :\, X_{t}^{(\sigma)} \notin (\!-\infty,0)\}.\end{equation}

We show the following result.

Theorem 2. (HMC for multimodal distributions.) The conductance $\Phi_{\sigma}$ satisfies

(4.3) \begin{equation} \lim_{\sigma \mapsto 0} (-2 \sigma^{2}) \log(\Phi_{\sigma}) = 1.\end{equation}

Furthermore, for all $\epsilon > 0$ and fixed $x \in (\!-\infty,0)$ , the hitting time $\tau_{x}^{(\sigma)}$ satisfies

(4.4) \begin{equation} \lim_{\sigma \mapsto 0} \mathrm{P} \left[\frac{\log(\tau_{x}^{(\sigma)})}{\log(\Phi_{\sigma})} < 1 + \epsilon \right] = 1\end{equation}

and the relaxation time satisfies

(4.5) \begin{equation} \lim_{\sigma \mapsto 0} \frac{\log(\lambda_{\sigma})}{\log(\Phi_{\sigma})} = \lim_{\sigma \mapsto 0} \frac{\log(\lambda_{\sigma})}{\log(\Phi(K_{\sigma}))} = 1.\end{equation}

We defer the proof to Appendix B.2. Note that (4.3)–(4.5) exactly match the spectral gap for an optimally tuned RWM algorithm with the same target distribution (see Theorem 3 of our companion paper [Reference Mangoubi, Pillai and Smith36]).

We also observe that this result implies the conductance $\Phi(K_{\sigma})$ of the kernel $K_{\sigma}$ is close to the conductance $\Phi_{\sigma} = \Phi(K_{\sigma}, (\!-\infty,0))$ associated with the specific test set $(\!-\infty,0)$ , at least for $\sigma$ very small. The set $(\!-\infty,0)$ is of course a natural guess for the set with the ‘worst’ conductance, though we do not know of any simple argument that would prove something like this. Deriving simple criteria to verify this fact was the motivation behind our companion paper [Reference Mangoubi, Pillai and Smith36].

4.2. High-dimensional analogue

Almost all of the work in the proof of Theorem 2 is checking that the spectral gap is not much smaller than the natural guess and upper bound $\Phi(K_{\sigma}, (\!-\infty,0))$ . If we merely wish to check that HMC is ‘slow’ for $\sigma$ small, it is straightforward to compute an upper bound on $\Phi(K_{\sigma}, (\!-\infty,0))$ that is good enough by using Theorem 1. This last fact remains true in higher-dimensional examples. For example, considering the mixture distribution

\begin{equation*}\pi_{\sigma} = \frac{1}{2} \mathcal{N}((-1,0,\ldots,0),\sigma^2 \mathrm{Id}) + \frac{1}{2}\mathcal{N}((1,0,\ldots,0),\sigma^2 \mathrm{Id})\end{equation*}

on $\mathbb{R}^{d}$ . It is natural to guess that, for the HMC algorithm with typical time step $T_{\sigma} = O(\sigma)$ , the set $\{ x \in \mathbb{R}^{d}\, :\, x[1] < 0 \}$ has (nearly) the worst conductance. Theorem 1 can be used to get a very good estimate of this conductance. In particular, the same calculation as inequality (B11) at the start of the proof of Theorem 2 immediately yields the bound

\begin{equation*}\Phi_{\sigma} \leq T_{\sigma} \left( \frac{1}{\sqrt{2 \pi} \sigma} \mathrm{e}^{-\frac{1}{2 \sigma^{2}}}\right)^{d}.\end{equation*}

4.3. Choice of integration time and computational cost

It is natural to ask: why do we study the choice $T_{\sigma} = \sigma$ , rather than some other choice of integration time? The simplest answer is that it is impossible to substantially improve the performance of the algorithm by choosing a larger integration time. To make this statement precise, we note that the computational cost of running a single step of the HMC algorithm is approximately proportional to the integration time T. Thus, the effective computational cost of a sample from a target distribution $\pi$ is roughly the ratio of the integration time to the spectral gap (this is a standard way to compare the efficiency of MCMC algorithms with widely differing costs per step; see [Reference Bornn, Pillai, Smith and Woodard10, Reference Sherlock, Thiery and Lee56] and references therein). It is this effective computational cost that is bounded below: we always have

\begin{equation*}\limsup_{\sigma \mapsto 0} 2 \sigma^{2} (\log(\lambda_{\sigma}) - \log(T_{\sigma})) \leq 1,\end{equation*}

for all $T_{\sigma} \geq \sigma$ (see the discussion in Section 5.2 for a proof of this fact). Thus, $T_{\sigma} = \sigma$ is essentially the optimal choice for $T_{\sigma}$ in this case.

Having said this, we find this simple answer slightly misleading. In practice, one never knows the optimal value of $T_{\sigma}$ . Instead, one often tunes an MCMC algorithm according to some ‘Goldilocks principle’. For HMC, one chooses $T_{\sigma}$ so that the probability of an HMC trajectory making a ‘U-turn’ (in the sense of [Reference Hoffman and Gelman28]) is not too close to 0 and not too close to 1. In this example, this means choosing $T_{\sigma} = \Theta(\sigma)$ . Note that this heuristic is very similar to the ‘Goldilocks principle’ for tuning RWM: one should choose the standard deviation of the proposal distribution so that the probability of rejecting a proposal is not too close to 0 and not too close to 1. This popular tuning choice is exactly the one that we study in our companion paper [Reference Mangoubi, Pillai and Smith36], though it is not close to optimal in that context. The difference comes from the fact that, unlike HMC, RWM can ‘jump over’ low-density regions. This means that the globally optimal choice of step size can be much larger than the ‘locally optimal’ choice that is typically made.

4.4. Application: HMC targeting degenerate multivariate Gaussian

This section is motivated by the study of the standard HMC Algorithm 2 with target distribution of the form $\mathcal{N}(0, M_{\sigma})$ , where the two-dimensional covariance matrix $M_{\sigma}$ is given by

\begin{equation*} M_{\sigma}= \left[ {\begin{array}{c@{\quad}c} 1 & 0 \\ 0 & \sigma^{2} \\ \end{array} } \right].\end{equation*}

We next observe that the target distribution $\mathcal{N}(0,M_{\sigma})$ is special: if $\{X_{t}\}_{t \in \mathbb{N}}$ is a Markov chain drawn from Algorithm 2 targeting $\mathcal{N}(0,M_{\sigma})$ , then the coordinate sequences $\{X_{t}[1]\}_{t \in \mathbb{N}}$ and $\{X_{t}[2]\}_{t \in \mathbb{N}}$ are each Markov chains as well. Furthermore, they both evolve independently, and they evolve according to Algorithm 2, with $\{X_{t}[1]\}_{t \in \mathbb{N}}$ targeting $\mathcal{N}(0,1)$ and $\{X_{t}[2]\}_{t \in \mathbb{N}}$ targeting $\mathcal{N}(0,\sigma^{2})$ . Thus, rather than analyzing the full chain $\{X_{t}\}_{t \in \mathbb{N}}$ , to calculate the spectral gap it is enough to analyze the slower-mixing marginal chain $\{X_{t}[1]\}_{t \in \mathbb{N}}$ .

Denote by $K_{\sigma}$ the transition kernel associated with target distribution $\mathcal{N}(0,1)$ and integration time $T_{\sigma}$ satisfying $T_{\sigma} = o(1)$ . Denote by $\Phi_{\sigma}$ the conductance associated with this transition kernel and the set $(\!-\infty,0)$ ; denote by $\rho_{\sigma}$ the spectral gap of $K_{\sigma}$ . We have the following result.

Theorem 3. The conductance $\Phi_{\sigma}$ satisfies

(4.6) \begin{equation} \lim_{\sigma \mapsto 0} \frac{\log(\Phi_{\sigma})}{\log(T_{\sigma})} \leq 1.\end{equation}

Furthermore, the relaxation time satisfies

(4.7) \begin{equation} \limsup_{\sigma \mapsto 0} \frac{\log(\rho_{\sigma})}{\log(T_{\sigma})} \geq \frac{1}{2}.\end{equation}

The proof is deferred to the appendix. Note that we do not give a bound for the hitting time in this example. This is not an accident: the target distribution is unimodal, and so the HMC algorithm does not exhibit metastability and the fluctuations of the hitting time remain large relative to the relaxation time as $\sigma$ goes to 0.

4.5. Numerical simulation: spectral gap of HMC

In this section we plot the spectral gap (Figure 3) of Algorithm 2 when $\pi$ is a two-mode density; the plot was generated by numerically diagonalizing an analytical solution for the transition matrix of the HMC Markov chain. The results of the calculation agree closely with the upper bounds on the spectral gap given by applying Corollary 1 with inequality (2.6).

Figure 3.

The spectral gap for the isotropic-momentum HMC algorithm with stationary distribution $\pi(q) = ({1}/{2F_{\mathcal{N}(0,1)}(a)})\max (f_{\mathcal{N}(0,1)}(q-a),f_{\mathcal{N}(0,1)}(q-a))$ , for different inter-modal distances 2a and different Hamiltonian trajectory times T. The results agree closely with the asymptotic formula in Theorem 1.

In this simulation we computed the spectral gap for Algorithm 2 with the stationary distributions of the form $\pi_{a}(q) = ({1}/{2F_{\mathcal{N}(0,1)}(a)})\max(f_{\mathcal{N}(0,1)}(q-a),f_{\mathcal{N}(0,1)}(q+a))$ , where $f_{\mathcal{N}(0,1)}$ and $F_{\mathcal{N}(0,1)}$ denote the probability distribution function (PDF) and cumulative distribution function (CDF) of the standard normal density. Note that $\pi_{a}(q)$ approximates the Gaussian mixture model $\tilde{\pi}_{a}(q) = \tfrac{1}{2}f_{\mathcal{N}(0,1)}(q-a) + \tfrac{1}{2}f_{\mathcal{N}(0,1)}(q+a)$ ; indeed $\lim_{a \mapsto \infty} \| \pi_{a} - \tilde{\pi}_{a} \|_{\mathrm{TV}} = 0$ .

As suggested by the formula for the conductance in Theorem 1, the spectral gap is bounded above by a linear function of T, and in fact increases approximately linearly with T when $a=0$ for $T\leq {\pi}/{2}$ . Note that, for fixed a small, the spectral gap $(1 - \lambda_{2})$ looks like a periodic function of T. This is due to the fact that the trajectories themselves are very close to periodic with period $\geq {\pi}/{2}$ , meaning that the expectation term in (3.3) also varies (approximately) periodically with T. The exponential decay in $a^2$ is explained by considering the set $S = (\!-\infty,0)$ and noting that $\int_{\partial S} \pi_{a}(q) \mathrm{d}_{\partial S}q = f_{\mathcal{N}(0,1)}(q-a)$ , so the corresponding term in Corollary 1 of Theorem 1 decreases exponentially in $a^2$ .

5.1. Application to numerical HMC and related methods

This paper has studied ‘ideal’ HMC, for which Hamilton’s equations (2.1) are solved exactly. In practice, HMC is always implemented with a numerical integrator such as the Verlet integrator (see discussion in [Reference Neal47]) or with a more complicated scheme that relies heavily on such an integrator (see, e.g., [Reference Hoffman and Gelman28] and recent variations such as [Reference Bou-Rabee, Carpenter, Kleppe and Liu11]). It is natural to ask which of our conclusions hold for practical HMC algorithms, and perhaps for related piecewise-deterministic Markov processes.

We first point out that a massive violation of our main qualitative conclusion is possible. For example, if we choose a very large step-size h for the Verlet integrator, it is possible to ‘jump over’ regions of very small probability (note that the step size h for the integrator is not the same as the integration length, which we have denoted T; one always has $h \leq T$ , but both in practice and in theory h is usually much smaller, as in e.g. [Reference Betancourt, Byrne and Girolami7]). This is closely analogous to the ‘jumping’ of optimal RWM mentioned in Section 4.3. We thought it was important to point out this potential violation, but we think it has little practical relevance: it would be very unusual to choose a step size that was so large compared with the local gradient.

We see, essentially, two approaches to generalizing our results for practically relevant choices of integrator and step size. The simpler option is to view numerical implementations of HMC as small perturbations of ‘ideal’ HMC. One could then (i) use continuity arguments to check that the conductance of ‘nice’ sets is similar under both algorithms (see, e.g., Lemma 6 of [Reference Chen, Dwivedi, Wainwright and Yu17] for an example of such a bound for HMC), and then (ii) check that the worst-case conductance is well-approximated by the conductance of ‘nice’ sets (as we do in our Lemma 2). This approach uses well-understood tools to approximate the conductance of an actual HMC chain, though we do not see any way to get an exact formula for the conductance this way.

A second option is to try to more closely mimic our proofs to obtain stronger results. Our arguments make crucial use of the following fact: if we sample (P, Q) from the stationary measure of HMC, then the ‘rate’ at which HMC paths cross the boundary $\partial S$ is constant throughout the full sample path (see (A3) for a precise definition of the ‘rate’). This fact is a consequence of the conservation of volume property for Hamilton’s equations. A quick calculation shows that, even for a single standard Gaussian in two dimensions, Verlet integrators do not have this property, and so directly mimicking our proof cannot possibly work. Nonetheless, we suspect that the variability is usually quite small, and so bounding this variation should allow one to obtain similar conductance bounds. We note that the Verlet integrator actually comes closer to preserving a quantity that is slightly different from the target Hamiltonian, called the shadow Hamiltonian (see the textbook [Reference Hairer, Wanner and Lubich26] for definitions). Using an appropriate shadow Hamiltonian, rather than the nominal Hamiltonian, may allow for tighter analyses.

It is natural to ask whether similar formulas exist for other piecewise-deterministic Markov processes, such as the bouncy particle sampler [Reference Bouchard-Côté, Vollmer and Doucet14]. We believe that essentially the same comments apply to these samplers as apply to numerical HMC: the ‘rate’ function is not constant and so the final conductance formula cannot be the same, but we expect perturbation arguments to yield good bounds.

5.2. Bounds on HMC improvements

One immediate consequence of Theorem 1 and Lemma 2 is that it is not possible to dramatically improve the performance of Algorithm 2 by tuning the trajectory integration time T. We give a quick discussion of this fact and some related open problems.

Fix a probability distribution with $C^{2}$ density $\pi$ and denote by $K_{T}$ the transition kernel defined by Algorithm 2 with parameters $\pi$ and $T > 0$ ; we emphasize the parameter T due to its importance in the following discussion. Define

\begin{align*}\Phi_{0}(\pi,S) &= \frac{1}{2} \frac{\int_{\partial S} \pi(q) \, \mathrm{d}_{\partial S}q }{\pi(S)}\\\Phi_{0}(\pi) &= \inf_{S \in \mathfrak{A}} \Phi_{0}(\pi,S),\end{align*}

for $S \in \mathfrak{A}$ . By Corollary 1 of Theorem 1, combined with Lemma 2, we have the upper bound

(5.1) \begin{equation} T^{-1} \Phi(K_{T}) \leq \Phi_{0}(\pi),\end{equation}

so that $\Phi(K_{T})$ is bounded by a linear function in the integration time T.

This observation is already in stark contrast to RWM as given in Algorithm 1, whose tuning parameter $\epsilon > 0$ can have very large effects on performance. As an illustration of this large change in performance, consider the family of transition kernels $\{ Q_{\sigma} \}_{0<\sigma <1}$ given by Algorithm 1 with target distribution $\pi_{\sigma}$ given in (4.1) and tuning parameter $\epsilon \equiv 10$ for all $\sigma$ . It is a straightforward exercise to check that $\log(\rho(Q_{\sigma})^{-1}) = O(\log(\sigma^{-1}))$ for $\sigma$ small. (For example, one could apply Theorem 5 of [Reference Rosenthal53] with ‘small set’ $(-5,5)$ and ‘Lyapunov function’ $V(x) = \mathrm{e}^{\|x\|}$ to obtain a bound on the rate of convergence of the walk to stationarity. One could then apply Theorem 2.1 of [Reference Roberts and Rosenthal51] to convert this bound on the convergence rate to a bound on the spectral gap.) Comparing this with the result of Theorem 3 of [Reference Mangoubi, Pillai and Smith36], we see that changing the tuning parameter $\epsilon$ of Algorithm 1 by a factor on the order of $\sigma$ will vastly decrease the size of the log of the relaxation time, from $\Theta(\sigma^{-2})$ to $O(-\log(\sigma))$ .

This dramatic performance improvement does not have obvious practical applications, since choosing a good tuning parameter for multimodal targets requires a more detailed understanding of the target than a user will typically have. Our point is just that such dramatic improvements are possible for RWM with multimodal targets, while they are not possible for HMC.

We interpret the linear bound (5.1) in terms of computational cost. Roughly speaking, the computational cost of a step of Algorithm 2 is proportional to the integration time T. Paraphrasing inequality (5.1), it is impossible to improve the cost-normalized conductance of Algorithm 2 by increasing the integration time T. This immediately implies, via Cheeger’s inequality (2.6), that the spectral gap is bounded by a quadratic in T. We point out that this behaviour is very similar to that of ‘lifted’ Markov chains (see [Reference Chen, Lovász and Pak16] for a definition). Roughly speaking, like HMC, ‘lifted’ Markov chains attempt to combat diffusive behaviour of RWM by adding an abstract notion of ‘momentum’. The central idea behind the non-improvement theorems for lifted chains in [Reference Chen, Lovász and Pak16, Reference Ramanan and Smith48] is that it is impossible to increase the conductance of a Markov chain by ‘lifting’ it. By Cheeger’s inequality, this gives an upper bound on the best-possible improvement due to lifting.

We close this discussion with an open question. We have shown that the conductance of a given set can increase at most linearly in T. Via Cheeger’s inequality (2.6), this suggests that the spectral gap should satisfy a similar quadratic inequality in T. To make this rigorous, it would be sufficient to show the following problem.

Open Problem 1. Prove the limit

(5.2) \begin{equation} \lim_{T \mapsto 0} \inf_{S \in \mathfrak{A}} T^{-1} \Phi(K_{T},S) = \inf_{S \in \mathfrak{A}} \Phi_{0}(\pi,S).\end{equation}

Note that equality (A21) already implies that, for fixed $S \in \mathfrak{A}$ ,

\begin{equation*}\lim_{T \mapsto 0} T^{-1} \Phi(K_{T},S) = \Phi_{0}(\pi,S).\end{equation*}

5.3. Multimodal targets

In Theorem 2 and Theorem 3 of the companion paper [Reference Mangoubi, Pillai and Smith36], we showed that Algorithm 2 has very similar performance to Algorithm 1 for a strongly multimodal example. However, our comparison is not direct: we prove that the two algorithms have similar spectral gaps by laboriously computing the spectral gaps of both algorithms. We suspect that this behaviour is quite general, and propose the following informal problem.

Open Problem 2. Let $\pi_{\sigma}$ be a mixture distribution of the form

\begin{equation*}\pi_{\sigma}(x) \propto \sum_{i=1}^{k} \mu_{i} f_{i}\left(\frac{x}{\sigma}\right),\end{equation*}

where $f_{1},\ldots,f_{k}$ are the densities of probability distributions and $\mu_{1},\ldots,\mu_{k} \geq 0$ are fixed weights that sum to $\sum_{i=1}^{k} \mu_{i} = 1$ .

Let $K_{\sigma}$ , $Q_{\sigma}$ be the transition kernels of Algorithms 2 and 1 respectively, with target distribution $\pi_{\sigma}$ , integration time $T_{\sigma} \propto \sigma$ , and standard deviation $\epsilon_{\sigma} \propto \sigma$ . We wish to know: what are sufficient conditions for $\{ \mu_{i}\}_{i=1}^{k}$ so that

(5.3) \begin{equation} \lim_{\sigma \mapsto 0} \frac{\log(\rho(K_{\sigma}))}{\log(\rho(Q_{\sigma}))} \geq 1?\end{equation}

Note that the constant ‘1’ in the limit (5.3) is important. We expect all algorithms to perform poorly for multimodal targets, and thus for the limit (5.3) to be strictly greater than 0. A limit of ‘1’ would suggest that HMC exhibits no asymptotic improvement over RWM; a limit of, e.g., $\frac{1}{2}$ would suggest a quadratic improvement over RWM, which is substantial.

We conjecture that inequality (5.3) holds as long as the level sets of $f_{i}$ are fairly ball-like; for example, if the sets $S_{i,C} \equiv \{x\, :\, f_{i}(x) \leq C\}$ are all convex, with the ratio of the inner and outer radii

\begin{align*}r_{\mathrm{inner}}(i,C) &\equiv \sup\{r\, :\, \exists \, x \in S_{i,C} \text{ s.t. } B_{r}(x) \subset S_{i,C} \} \\r_{\mathrm{outer}}(i,C) &\equiv \inf\{r\, :\, \exists \, x \in S_{i,C} \text{ s.t. } B_{r}(x) \supset S_{i,C} \} \end{align*}

uniformly bounded. We mention that we do not expect this to hold in complete generality. For example, if the level sets $S_{i,C}$ of $\pi_{\sigma}$ for $0 <C \ll \mathrm{e}^{\sigma^{-1}}$ look like Figure 1.

More generally, we are interested in determining which ‘momentum-based’ algorithms satisfy the limit (5.3). We note that not all similar chains satisfy this equality. In particular, lifted Markov chains can replace the ‘1’ with a ‘ $\tfrac{1}{2}$ ’ for realistic examples (see, e.g., [Reference Bierkens and Roberts9]). It would be valuable to determine if other generic algorithms, such as Algorithm 3 or [Reference Bierkens, Fearnhead and Roberts8], can also achieve this.

Finally, we suggest a more straightforward question. The level sets of the multimodal densities studied in this paper are completely disconnected below some (fairly large) threshold. Thus, the multimodality is essentially invisible to the gradient of the log-likelihood within each mode. Can HMC, and especially Riemannian manifold HMC under an appropriate implementable choice of Riemannian metric, improve on metropolitan–Hastings for multimodal examples in which the level sets are connected, and multimodality is induced by these level sets being merely very narrow? Note that this is only possible in two or more dimensions.

Open Problem 3 (Mixing time upper bounds?). Our formula for the conductance does not directly lead to mixing time upper bounds for idealized HMC. It is therefore an interesting open problem whether our conductance formula can be used as a starting point to derive mixing time upper bounds.

That being said, we note that, in subsequent works (e.g. [Reference Chen, Dwivedi, Wainwright and Yu17, Reference Chen and Gatmiry18]) appearing after the appearance of a preprint of our paper on arXiv, mixing time upper bounds were obtained for a numerical implementation of HMC via a conductance argument. However, deriving upper mixing time bounds from our conductance formula for idealized HMC would be highly non-trivial and would require overcoming a number of challenges that are beyond the scope of our paper. We describe some of the challenges below.

Recall that our formula for the conductance of a partition $\mathbb{R}^d = S \cup S^c$ of Euclidean space (Theorem 3.1) takes as input the number of times $N_{\partial S}$ that the HMC trajectory crosses the boundary $\partial S$ of the partition S. When proving a lower bound, this is not an issue as the quantity $({1}/{N_{\partial S}}) \unicode{x1D7D9}\{N_{\partial S} \textrm{ odd}\}$ which appears in our formula is trivially upper bounded by 1. This immediately implies a simple upper bound on the conductance that does not require as input the number of crossings $N_{\partial S}$ (our Corollary 1), which in turn implies a lower bound on the mixing time via Cheeger’s inequality.

Unfortunately, showing a lower bound on the conductance (which would imply an upper bound on the mixing time, via Cheeger’s inequality) requires one to show a lower bound on the (expectation of) the quantity $({1}/{N_{\partial S}}) \unicode{x1D7D9}\{N_{\partial S} \textrm{ odd}\}$ . To show an upper bound on the mixing time, one must bound the infinum of the conductance over all possible choices of partitions S. For many choices of S, the boundary $\partial S$ is highly non-smooth, which means that even a very smooth HMC trajectory can have a large number of intersections with the boundary $\partial S$ . This makes it very challenging to obtain a lower bound on the quantity $({1}/{N_{\partial S}}) \unicode{x1D7D9}\{N_{\partial S} \textrm{ odd}\}$ .