Does Hamiltonian Monte Carlo mix faster than a random walk on multimodal densities?
Hamiltonian Monte Carlo (HMC) is a very popular and generic collection of Markov chain Monte Carlo (MCMC) algorithms. One explanation for the popularity of HMC algorithms is their excellent performance as the dimension $d$ of the target becomes large: under conditions that are satisfied for many common statistical models, optimally-tuned HMC algorithms have a running time that scales like $d^{0.25}$. In stark contrast, the running time of the usual Random-Walk Metropolis (RWM) algorithm, optimally tuned, scales like $d$. This superior scaling of the HMC algorithm with dimension is attributed to the fact that it, unlike RWM, incorporates the gradient information in the proposal distribution. In this paper, we investigate a different scaling question: does HMC beat RWM for highly $\textit{multimodal}$ targets? We find that the answer is often $\textit{no}$. We compute the spectral gaps for both the algorithms for a specific class of multimodal target densities, and show that they are identical. The key reason is that, within one mode, the gradient is effectively ignorant about other modes, thus negating the advantage the HMC algorithm enjoys in unimodal targets. We also give heuristic arguments suggesting that the above observation may hold quite generally. Our main tool for answering this question is a novel simple formula for the conductance of HMC using Liouville's theorem. This result allows us to compute the spectral gap of HMC algorithms, for both the classical HMC with isotropic momentum and the recent Riemannian HMC, for multimodal targets.
\section{Introduction}
Markov chain Monte Carlo (MCMC) algorithms, and in particular the Hamiltonian Monte Carlo (HMC) algorithms, are workhorses in many scientific fields including physics \cite{Hybrid_MCMC}, statistics  and machine learning \cite{NUTS,Riemannian_HMC2, MCMC_Application_Machine_Learning,welling2011bayesian, mangoubi2018convex}, and molecular biology \cite{MCMC_Application_molecular_biology}. An important question practitioners face is how to choose an algorithm for a \textit{particular} problem? Are there any general principles that allow users to prefer one algorithm to another for \textit{generic} but \textit{simple} problems? One of the most fruitful branches of MCMC theory has focused on analyzing MCMC algorithms by studying their \textit{scaling limits.} The idea is to study how an MCMC algorithm's speed depends on some underlying parameter, often the dimension $d$ of the target distribution. Results from the theory of scaling limits suggest that HMC often has a running time of  $O(d^{0.25})$ \cite{beskos2013optimal} whereas  the Random-Walk Metropolis (RWM) algorithm has a much longer running time of $O(d)$ steps under very broad conditions \cite{roberts1997weak,mps12,bed07}. These results and other evidence, both empirical and theoretical (see \textit{e.g.,}  \cite{betancourt2017conceptual,mangoubi2017rapidp1,rabee2018HMCcoup, mangoubi2016rapid}) have lead to the widespread understanding that HMC is superior to RWM for a wide variety of problems. In this paper, we investigate the extent to which this superiority holds in another natural scaling regime: highly multimodal target distributions. There is no a priori reason to believe that algorithms superior in one regime are also superior in other regimes.

 To give a simple example, consider the mixture of Gaussians
\begin{equs} \label{eqn:gausmix}
\frac{1}{2} \mathcal{N}(-1,\sigma^2) + \frac{1}{2} \mathcal{N}(1,\sigma^2).
\end{equs} 
A natural regime to study is when $\sigma \rightarrow 0$. Let us tune both algorithms to optimize their performance within a single mode. To this end, for the RWM algorithm, results of \cite{roberts1997weak} (and also simple scale-invariance of the mixture components) imply that the proposal variance must be $O(\sigma^2)$. For the classical HMC algorithm  \cite{Hybrid_MCMC} with isotropic momentum, we set the integration length to be $O(\sigma)$; this is the first time we would expect the algorithm to take a U-turn (and is also suggested by scale-invariance of the mixture components). For the above example, with these tuning parameters, Theorem \ref{ThmHmcMultimodal} of this paper and Theorem 3 of our companion paper \cite{mangoubi2018simple} imply that the spectral gap of HMC and RWM \emph{both} decay exactly like $e^{-\frac{1}{2} \sigma^{-2}}$  as $\sigma \rightarrow 0$. Of course, we expect that the spectral gap goes to 0 quickly in both cases, just as it does in the high-dimensional scaling regime of \cite{roberts1997weak,beskos2013optimal}. The interesting fact is that they go to 0 at the \textit{same rate}, so that the asymptotic performance of the two algorithms is the same. For this example, instead of isotropic momentum, if we choose the metric to be the inverse of the Fisher information as in the Reimannian HMC of \cite{Riemannian_HMC2}, the spectral gap still decays like $e^{-\frac{1}{2} \sigma^{-2}}$  as $\sigma \rightarrow 0$.  Thus, our main conclusion is that, for this example is that HMC is \emph{not} better than RWM. 
Finally, the above results also hold for the high dimensional analog of \eqref{eqn:gausmix}; see Section \ref{sec:highdim}.

It is natural to ask if this comparison is a result of bad tuning. In fact, this is not the case. Tuning the integration length of HMC cannot improve its relative performance. The natural tuning for HMC used above is essentially optimal (in the sense of maximizing the effective sample size per unit computation; see Section \ref{sec:RemChoiceInt} for a more detailed discussion of tuning parameters and optimality). We give a further discussion of scaling results, as well as related open questions, in Section \ref{SecCompMulti}. We highlight here one fact from  Section \ref{SecCompMulti} that is particularly interesting for our simple target distribution \eqref{eqn:gausmix}: considering other tuning parameters makes HMC look even worse than RWM. More precisely, the computational cost of the HMC algorithm is \textit{always} at least on the order of $e^{-\frac{1}{2} \sigma^{-2}}$ even for much longer integration times, while well-chosen tuning parameters can make RWM much more efficient (again detailed in Section  \ref{sec:RemChoiceInt} and Section \ref{SecCompMulti}). 



Readers and practitioners familiar with the MCMC literature may be surprised at the tone of the above paragraph, since HMC is generally viewed as having much better performance than many older and simpler random walk based MCMC algorithms  \cite{betancourt2017conceptual,roberts1997weak,beskos2013optimal,mangoubi2017concave,rabee2018HMCcoup}. However, the main heuristics justifying the dominance of HMC in the high-dimensional regime do not apply in our highly-multimodal regime. In particular, HMC can use knowledge of the gradient of the log-likelihood to make larger moves than RWM in the high-dimensional regime, but the gradient cannot detect multimodality. 

\subsection{Conductance Results}
One of the main tools we use in our calculations is \textit{Cheeger's inequality}. This well known tool provides a coarse estimate of the efficiency of an MCMC algorithm in terms of a geometric quantity called the \textit{conductance} (see\cite{cheeger1970lower, lawler1988bounds} and a survey of variants \cite{montenegro2006mathematical}). It is useful primarily because it can be easier to estimate than more precise measures of efficiency. 

Unfortunately, the conductance can still be difficult to compute. One of the main contributions of this paper is a simple exact formula for the conductance of both HMC and RHMC (see Theorem \ref{thm:1}). This formula relates the conductance of a set $S$ to the integral of a function, which we interpret as its ``crossing rate," on the boundary $\partial S$ of $S$.  For the `standard' HMC algorithm that uses isotropic momentum, Theorem \ref{thm:1} (see Corollary \ref{cor:1}) yields the following bound for conductance of a set $S$:
\begin{equs} 
\mathrm{Conductance}(S) \leq \frac{1}{2}T { \int_{\partial S} \pi(q) \mathrm{d}q \over \pi(S)},
\end{equs} 
where $T$ is the integration length and the integral is over the boundary $\partial S$ of the set $S$ and $\pi$ is the target density. In particular, the above formula
yields that the conductance can increase \emph{at most} linearly with the integration time for the HMC algorithm. This is also true for RHMC. In Section \ref{sec:RemChoiceInt}, we use this observation to show that changing the integration length cannot vastly improve the performance of the HMC algorithm in multimodal regimes.

 Markov chains do not typically have conductance formulas that can easily be expressed in terms of surface integrals. As seen from above, our new formula is much easier to use.  Having a simple formula for HMC is particularly useful because HMC algorithms are some of the most widely-used \cite{Riemannian_HMC2, cheung2009bayesian, mehlig1992hybrid} modern MCMC algorithms, but their theoretical properties are not yet well-understood (though again see \cite{holmes2014curvature, livingstone2016geometric, mangoubi2017concave, rabee2018HMCcoup}). Understanding the conductance of HMC algorithms can help us understand when these algorithms perform poorly, and often suggest what has gone wrong. 

As with \textit{e.g.,} \cite{holmes2014curvature}, we begin our quantitative study of mixing times by studying ``idealized" versions of HMC. That is, we ignore the error introduced by solving Hamilton's equations with a numerical integrator such as the leapfrog integrator and instead assume that we can solve the Hamilton's ODEs exactly. We expect similar qualitative conclusions to hold for appropriate practical implementations of HMC, and there is substantial work on relating ``ideal" Monte Carlo schemes to their numerical implementations (see \textit{e.g.,}  \cite{durmus2016sampling2} \cite{livingstone2016geometric}, \cite{mangoubi2017rapidp1}, \cite{mangoubi2017rapidp2}, \cite{lee2017convergence} \cite{mangoubi2018dimensionally}).



\subsection*{When can HMC beat RWM for multimodal targets?}
\begin{figure}
\begin{tikzpicture}[>=latex]
\coordinate (A) at (0,0) ;
\coordinate (B) at (0,-5) ;
\coordinate (C) at (4,0) ;
\coordinate (D) at (4,-5) ;

\filldraw[fill=gray!20]
 (A) to[out=-100,in=100] 
  (A)
 
 (C) to[out=-80,in=80] 
   coordinate[pos=0.05] (auxru)
   coordinate[pos=0.95] (auxrl)
 (D) to[out=60,in=-60] 
 (C);
\node[circle,draw = black, fill = gray!20,inner sep=0.5cm] 
  at (1.8,-2.5) (circle) {};    
  \end{tikzpicture}
  \caption{We see a small deep mode next to a `banana' shaped mode. HMC can often explore long, skinny modes much more quickly than RWM. Thus, it is possible to tune the length of the banana shaped mode, in relation to the distance between the centers of the banana and the circle, so that the time for RWM to mix on the long mode is much larger than the time to escape the small mode, while the time for HMC to mix on the long mode is much smaller than the time to escape the small mode. In this case, the HMC algorithm can exhibit metastability while the RWM algorithm does not; we would expect HMC to mix more quickly than RWM in this situation.} \label{FigBanana}
  \end{figure}

We envision that there may be some classes of highly multimodal target densities where the HMC algorithm is superior to random walk based algorithms. One such example is shown in Figure \ref{FigBanana}. The basic idea is that HMC can often explore long, skinny modes much more quickly than RWM. If we consider an example with several modes, at least one of which is very long and skinny, the running time for  RWM may be determined by the time it takes to traverse the long and skinny mode, whereas the running time for HMC will be determined by the time it takes to travel between the modes. This `energy-entropy' competition could lead to different performances of the two algorithms, and certainly merits further study. In particular, the HMC algorithm in this example could exhibit `metastability' while the RWM algorithm does not. In our companion paper \cite{mangoubi2018simple}, we derive simple conditions for checking the metastability of commonly used MCMC algorithms.


It is also natural to ask if any of the popular HMC variants can substantially improve the performance of HMC on multimodal targets. We believe that this is an important research question, but there are many HMC variants and a careful survey of their performance is beyond the scope of the present paper; instead we give a quick summary in relation to our results. Our main results do not apply as stated to the popular NUTS algorithm \cite{NUTS}, though we expect very similar bounds to be true. Our main conductance bounds, Theorem \ref{thm:1}  and Corollary \ref{cor:1}, do apply as stated to the Riemannian HMC (RHMC) algorithm of \cite{Riemannian_HMC2}. For the toy example \eqref{eqn:gausmix}, these bounds can be used to show that the RHMC algorithm with its ``usual" tuning (the inverse of the Fisher information as suggested in \cite{Riemannian_HMC2}) does not offer substantial improvements on the RWM algorithm with the ``usual" tuning discussed above. On the other hand, it does appear possible to use extremely unusual tuning parameters for RHMC to substantially improve performance, just as this is possible for RWM (see Section \ref{SecHmcImp} for detailed discussion of this tuning for RWM). As with RWM, there does not seem any obvious way to ``guess" these good tuning parameters without unrealistically good knowledge of the locations of the modes, and so it is not clear if the existence of good tuning parameters has any practical importance. 



\subsection{Guide to Paper}

In Section \ref{SecAlgDefs}, we review basic notation and the definitions of the MCMC algorithms that we study in this paper: an ``ideal" version of the standard and Riemannian HMC algorithms, as well as a simple Metropolis-Hastings algorithm for comparison. We then prove bounds on the conductance for HMC in Section \ref{sec:HMC}. Finally, we illustrate the uses of our results by analyzing the performance of HMC for two simple but important examples in Section \ref{SecAppl}: a mixture of Gaussians and a highly degenerate Gaussian. Finally, we discuss the consequences of our work and give related conjectures in Section \ref{SecCompMulti}. The longer proofs are deferred to the appendices.

\subsection{Basic Notation}
We review the standard ``big-O" notation, which is used throughout the paper. For two nonnegative functions or sequences $f,g$, we write $f = O(g)$ as shorthand for the statement: there exist constants $0 < C_{1},C_{2} < \infty$ so that for all $x > C_{1}$, we have $f(x) \leq C_{2} \, g(x)$. We write $f = \Omega(g)$ for $g = O(f)$, and we write $f = \Theta(g)$ if both $f= O(g)$ and $g=O(f)$. Relatedly, we write $f = o(g)$ as shorthand for the statement: $\lim_{x \rightarrow \infty} \frac{f(x)}{g(x)} = 0$.  We write $f = \tilde{O}(g)$ if there exist constants $0 < C_{1},C_{2}, C_{3} < \infty$ so that for all $x > C_{1}$, we have $f(x) \leq C_{2} \, g(x) \log(x)^{C_{3}}$, and write $f = \tilde{\Omega}(g)$ for $g = \tilde{O}(f)$. Finally, we say that a function $f$ is ``bounded by a polynomial" if there exists $0 < c < \infty$ so that $f(x) = O(x^{c})$. 
\section{Algorithms and Notation} \label{SecAlgDefs}


In this section we review two important Hamiltonian Monte Carlo algorithms, as well as the commonly-used Random Walk Metropolis algorithm.  We also review some important definitions for MCMC algorithms, including careful definitions of the conductance, spectral gap and Cheeger inequality. 

\subsection{Basic Notation}
Throughout the remainder of the paper, we denote by $\pi$ the smooth density function of a probability distribution on $\mathbb{R}^{d}$. We denote by $\mathcal{L}(X)$ the distribution of a random variable $X$. Similarly, if $\mu$ is a probability measure, we write ``$X \sim \mu$" for ``$X$ has distribution $\mu$."  Throughout, we will generically let $Q \sim \pi$ and $P \sim \mathcal{N}(0, \mathrm{Id})$ be independent random variables, where $\mathrm{Id}$ is the $d$-dimensional identity matrix.


\subsection{Random Walk Metropolis}
The Random Walk Metropolis (RWM) algorithm (Algorithm \ref{alg:RWM}) is the most basic commonly-used MCMC algorithm.  At each step $i$ of the Markov chain, the RWM algorithm proposes to take the next step $\hat{x}_{i+1}$ in a random direction and distance from the current position $x_i$.  The step is accepted with a probability of $\mathrm{min}\{\frac{\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. 


\begin{algorithm}[H]
\caption{Random Walk Metropolis \cite{metropolis1953equation}}\label{alg:RWM}
\flushleft
\textbf{Input:} Tuning parameter $\epsilon > 0$, starting point $x_0$, target density $\pi: \mathbb{R}^{d} \rightarrow [0, \infty)$. \\
\textbf{Output:} Markov chain $x_1, x_2, \ldots$, with stationary distribution $\pi$. \\
\begin{algorithmic}[1]
\FOR{$i = 1,2,\ldots$,}
\STATE Sample $p_{i} \sim \mathcal{N}(0,\epsilon^{2})^{d}$ and $U \sim \mathrm{Unif}([0,1])$. 
\STATE Set $\hat{x}_{i+1} = x_{i} + p_{i}$.
\IF{$U < \mathrm{min}\{\frac{\pi(\hat{x}_{i+1})}{\pi(x_i)}, 1\},$}
\STATE Set $x_{i+1} = \hat{x}_{i+1}$.
\ELSE  
\STATE Set $x_{i+1} = x_{i}$.
\ENDIF
\ENDFOR
\end{algorithmic}
\end{algorithm}



\subsection{Hamiltonian Monte Carlo algorithms}


Hamiltonian Monte Carlo (HMC) algorithms \cite{Hybrid_MCMC} 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 \textit{e.g.,} \cite{roberts2001optimal,beskos2013optimal} 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 smooth function $H \, : \, \mathbb{R}^{2d} \mapsto \mathbb{R}^{+}$, called the \textit{Hamiltonian function}, and starting points $p_{0},q_{0} \in \mathbb{R}^{d}$. We then define \textit{Hamilton's equations} to be the pair of differential equations
\begin{equs} \label{EqHamEq}
\frac{d}{dt} p(t) = - \frac{\partial H(p,q)}{\partial q}, \quad \frac{d}{dt} q(t) = \frac{\partial H(p,q)}{\partial p},
\end{equs}
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 \textit{e.g.}, $\gamma_{p,q}'(t) = p(t)$ for the ``standard" choice of Hamiltonian in Equation \eqref{EqStandardHam}.  The associated Hamiltonian function $H$ 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 \textit{conservation of energy} and \textit{conservation of volume}: 



\begin{enumerate}
\item \textbf{Conservation of Energy:} For $\gamma_{p,q}$ as above, we have 
\begin{equs} \label{EqConsEn}
\frac{d}{dt} H(\gamma_{p,q}'(t), \gamma_{p,q}(t)) \equiv 0.
\end{equs}
\item \textbf{Invariant distribution:} Assume that $\int_{p,q} e^{-H(p,q)} dq dp = 1$, and write $\mu_{H} = e^{-H(p,q)}$. If we sample $(P,Q) \sim \mu_{H}$, then we also have 
\begin{equs}  \label{EqConsVol}
\mathcal{L}(\gamma_{P,Q}'(t), \gamma_{P,Q}(t)) = \mu_{H}
\end{equs} for all $t \in \mathbb{R}^{+}$. This fact is a restatement of Liouville's theorem in probabilistic language.
\end{enumerate}

We refer the reader to \cite{lanczos1949variational} for proofs of these facts. The second fact motivates the standard choice of Hamiltonian
\begin{equs} \label{EqStandardHam}
H(p,q) = -\log(\pi(q)) + \frac{1}{2} \| p \|^{2}.
\end{equs} 
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 isotropic-momentum HMC algorithm (Algorithm \ref{alg:isotropic_HMC}), developed in \cite{mehlig1992hybrid}:


\begin{algorithm}[H]
\caption{Isotropic-Momentum HMC (idealized symplectic integrator) \cite{mehlig1992hybrid}} \label{alg:isotropic_HMC}
\flushleft
\textbf{Input:} Starting point $q_0 \in \mathbb{R}^{d}$, integration time $T \in \mathbb{R}^{+}$, smooth target density $\pi: \mathbb{R}^n \rightarrow [0, \infty)$. \\
\textbf{Output:} Markov chain $q_1, q_2, \ldots$, with stationary distribution $\pi$. \\
\begin{algorithmic}[1]
\STATE Define $H(p,q) := -\mathrm{log}(\pi(q)) + \frac{1}{2} \| p \|^{2}$.
\FOR{$i = 1,2,\ldots$,}
\STATE Sample $p_i \sim \mathcal{N}(0,1)^d$. 
\STATE Set $q_{i+1} = \gamma_{p_{i}, q_{i}}(T)$.
\ENDFOR
\end{algorithmic}
\end{algorithm}


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 Riemmanian manifold, with metric defined by some positive definite matrix $G(q)$.  For the special case of isotropic-momentum HMC, we have $G(q) = I_d$.  Alternatively, one may choose $G(q)$ to be a regularization of the Hessian of $U(q)$, which acts as a local pre-conditioner for $U$ \cite{betancourt2013general}. This can result in much faster mixing for distributions that are not close to isotropic; for example, it avoids the bad performance of Algorithm \ref{alg:isotropic_HMC} on the nearly-degenerate Gaussian that is observed in Theorem \ref{ThmHmcDegenerate}. The algorithm is identical to Algorithm \ref{alg:isotropic_HMC}, except for the choice of the Hamiltonian $H$: 



\begin{algorithm}[H]
\caption{Riemannian Manifold HMC (idealized symplectic integrator) \cite{Riemannian_HMC1, Riemannian_HMC2}} \label{alg:Riemannian_algorithm}
\flushleft
\textbf{Input:} Starting point $q_0 \in \mathbb{R}^{d}$, Integration time $T \in \mathbb{R}^{+}$, smooth target density $\pi: \mathbb{R}^d \rightarrow [0, \infty)$. \\
\textbf{Input:} Positive definite matrix $G(q)$ defining a Riemannian metric at every $q\in \mathbb{R}^d$. \\
\textbf{Output:} Markov chain $q_1, q_2, \ldots$, with stationary distribution $\pi$. \\
\begin{algorithmic}[1]
\STATE Define $H(p,q) := -\mathrm{log}(\pi(q)) +   \, \frac{1}{2}\log\left((2\pi)^d\det(G(q))\right) + \frac{1}{2} \langle p^\top G^{-1}(q), p \rangle$
\FOR{$i = 1,2,\ldots$,}
\STATE Sample $p_i \sim \mathcal{N}(0, G^{-1}(q_i))$. 
\STATE Set $q_{i+1} = \gamma_{p_{i}, q_{i}}(T)$.
\ENDFOR
\end{algorithmic}
\end{algorithm}

We will usually make the following additional assumption on the Riemannian metric matrix $G(q)$ in this case:

\begin{assumption} \label{AssumptionBoundFIM}
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  
\begin{equs} 
\sup_{q \in A} \max \left( \frac{1}{\lambda_{0}(q)}, \, \lambda_{1}(q) \right) < \infty.
\end{equs} 
\end{assumption}

Assumption \ref{AssumptionBoundFIM} says that the singular values of the 
metric $G$ are bounded away from $0$ and $\infty$ on compact sets. This is satisfied when $G$ is continuous and positive definite. More effort is required
when $G$ is allowed to be positive semi-definite; we do not study this case.

\subsection{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{equs}
(L  f)(x) = \int_{y \in \mathbb{R}^{d}} L(x,dy) f(y).
\end{equs}
The constant function is always an eigenfuncton of this operator, with eigenvalue 1. We define the space $W^{\perp} = \{ f \in L_{2}(\mu) \, : \, \int_{x} f(x) \mu(dx) = 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 \textit{spectral gap} $\rho$ of $L$ by the formula 
\begin{equs}
\rho = \rho(L) \equiv 1 - \sup \{ |\lambda| \, : \, \lambda \in \mathrm{Spectrum}(L^{\perp}) \},
\end{equs}
where $\mathrm{Spectrum}$ refers to the usual spectrum of an operator. If $L^{\perp}$ has a largest eigenvalue $|\lambda_2|$, then $\rho = 1-| \lambda_2 |$.  Geometric ergodicity of HMC algorithms was proved under very general conditions in \cite{livingstone2016geometric, Nawaf}, implying existence of a non-zero spectral gap under those conditions \cite{roberts1997geometric}.


Cheeger's inequality \cite{cheeger1970lower,lawler1988bounds} 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{equs} 
\Phi(L) &= \inf_{S  \in \mathcal{A} \, : \, 0 < \mu(S) < \frac{1}{2}} \Phi(L,S) \\
\Phi(L,S) &= \frac{ \int_{x} \mathbbm{1}\{x \in S\} L(x,S^{c}) \mu(dx)}{\mu(S) },
\end{equs} 
where $\mathcal{A}= \mathcal{A}(\mathbb{R}^{d})$ denote the usual collection of Lebesgue-measurable subsets of $\mathbb{R}^{d}$.

Cheeger's inequality for Markov chains, first proved in \cite{lawler1988bounds}, gives:
\begin{equation} \label{IneqCheegPoin}
\frac{\Phi(L)^2}{2} \leq \rho(L) \leq 2 \Phi(L).
\end{equation}




\section{Conductance of Hamiltonian Monte Carlo} \label{sec:HMC}


We derive equations for the conductance of HMC. Our main result, Theorem \ref{thm:1}, is based on the following two heuristics about HMC started at stationarity:
\begin{enumerate}
\item The conductance of a Markov chain measures the probability that a Markov chain will jump from some set $S$ to its complement in a single step. In the special case of HMC, this is roughly equivalent to the probability that the path $\gamma_{P,Q}$ will cross the boundary $\partial S$ an odd number of times.
\item By Liouville's theorem and reversibility, the \textit{rate} $\frac{d}{dt} | \gamma_{P,Q}([0,t]) \cap (\partial S)|$ at which the path $\gamma_{P,Q}$ intersects $\partial S$ is constant.
\end{enumerate}


Theorem \ref{thm:1} describes the conductance of a set purely in terms of this constant ``crossing rate" and the probability that the total number of crossings will be odd. Unfortunately, defining ``crossings" in a precise way requires a substantial amount of additional notation, which we give in Section \ref{SubsecNotationManifolds}. Readers primarily interested in upper bounds on the conductance may skip to Corollary \ref{cor:1}, which gives an easier-to-state upper bound on the conductance that is often close to sharp.

Although we introduce some restrictions on the set $S$ before giving our formula for the conductance of a \textit{particular set} $S$, we later prove that is possible to compute the conductance of \textit{the entire Markov chain} exactly using only the sets that satisfy these assumptions; see Lemma \ref{LemmSuffLoc} for details. Thus, our formula can in fact be used to compute the conductance exactly.

\subsection{Notation Related to Manifolds and Intersections} \label{SubsecNotationManifolds}

Fix a set $S \subset \mathbb{R}^{d}$ whose topological boundary $\partial S$  is an embedded $(d-1)$-dimensional manifold. Fix also an integration time $T \in \mathbb{R}^{+}$. In this special case, we define the \textit{number of intersections} $N_{\partial S}$ as follows:

\begin{defn} \label{DefNumInt}
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:
\begin{enumerate}
\item The function $\gamma_{p,q} \, : \, \mathbb{R} \mapsto \mathbb{R}^{d}$ 
is \textit{transverse} to the manifold $\partial S$.
\item The set $\{ t \in [0,T] \, : \, \gamma_{p,q}(t) \in \partial S\}$ is finite.
\end{enumerate}

For $(p,q) \in \mathcal{G}_{S}$, define
\begin{equs}
N_{\partial S} = N_{\partial S}(p,q) \equiv |\{ t \in [0,T] \, : \, \gamma_{p,q}(t) \in \partial S\}|.
\end{equs}
\end{defn} 


To recall a precise definition of a \textit{transverse} intersection, see \textit{e.g.,} \cite{guillemin2010differential}. 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 \ref{FigTransverse}.

\begin{figure}[H]
\includegraphics[scale=0.4]{TransverseNot}
\caption{A boundary $\partial S$ is shown in black. The intersection of $\partial S$ with the blue curve is transverse; the intersection with the red curve is not.\label{FigTransverse}}
\end{figure}




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 simple notation:

\begin{itemize}
\item 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{equs}
p_{q} = \langle G^{-1}(q) p, \eta(q) \rangle
\end{equs}

to be the component of $p$ in the direction orthogonal to $\partial S$ at q.
\item For $q \in \partial S$, define the set 

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

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

For $x \in \partial S$, denote by $\mathcal{T}_{x}$ the tangent space of $\partial S$ at $x$. 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 $\mathbbm{Proj}_{x} \, : \, \mathbb{R}^{d} \mapsto \mathcal{T}_{x}$ the projection map from $\mathbb{R}^{d}$ to the tangent space. We assume:



\begin{assumption} [Locally Well-Behaved Manifold] \label{DefLocallyWellBehaved}
We say that $S$ is \textit{locally well-behaved} if, for every compact set $A \subset \mathbb{R}^{d}$, there exist constants $0 < \epsilon, C < \infty$ so that:
\begin{enumerate}
\item For every ball $A' \subset A$ of diameter less than $\epsilon$, the set $A' \cap \partial S$ has a single component, and
\item For every $x, y \in (A \cap \partial S)$,
\begin{equs}
\| y - \mathbbm{Proj}_{x}(y) \| \leq C \|x - \mathbbm{Proj}_{x}(y) \|^{2}
\end{equs}
and 
\begin{equs}
\| \eta(x) - \eta(y) \| \leq C \|x-y \|.
\end{equs}
\end{enumerate}
\end{assumption}



\begin{remark}
We have not tried to provide the weakest-possible assumptions in Assumption \ref{DefLocallyWellBehaved}. Since Lemma \ref{LemmSuffLoc} implies that the conductance of most realistic HMC chains can be computed entirely in terms of sets that satisfy Assumption \ref{DefLocallyWellBehaved}, these assumptions are weak enough for our purposes. Here, we explain why something like Assumption \ref{DefLocallyWellBehaved} 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, we wish to avoid sets such as
\begin{equs}  \label{BadSetSDef}
S_{\infty} = \cup_{n \in \mathbb{N}} [0,1] \times [(4n)^{-2}, (4n+1)^{-2}].
\end{equs}

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. 
\end{remark}



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

\begin{lemma} \label{LemCountingPossible}
Set notation as above. Then the set $\mathcal{G}_{S}^{c}(T)$ has measure 0 for any $T \in \mathbb{R}^{+}$.
\end{lemma}

\begin{proof}
The proof is deferred to Appendix \ref{SecCountingPossible}.
\end{proof}


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 \ref{DefNumInt} for almost every value of $(p,q) \in \mathbb{R}^{2d}$. Throughout the rest of the paper, we will 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 two families of tilted measures on $\mathbb{R}^{2d}$:

\begin{defn} [Tilted Measures]

Let $\tilde{\mathbb{P}}$ be the probability measure on $\mathbb{R}^{2d}$ with density
\begin{equs}
\tilde{\mathbb{P}}(p,q) \propto e^{-H(p,q)} \, \mathbbm{1} \{ N_{\partial S}(p,q) \geq 1 \}.
\end{equs}

Define $\mathbb{Q}$ to be the probability measure on $\mathbb{R}^{2d}$ with density
\begin{equs} 
\mathbb{Q}(p,q) = \tilde{\mathbb{P}}(p,q) \cdot N_{\partial S}(p,q).
\end{equs}
\end{defn}







For a set $A \subset \mathbb{R}^{d}$ and a constant $c > 0$, define the $c$-thickening of $A$ to be 
\begin{equs} 
A_{c} = \{x \in \mathbb{R}^{d} \, : \, \inf_{a \in A} \, \|x-a\| \leq c\}.
\end{equs} 


\subsection{Main Conductance Formula}

Our main theorem is: 


\begin{thm} \label{thm:1} 

Let $T \in \mathbb{R}^{+}$, let $\pi(q)$ be any smooth probability density on $\mathbb{R}^d$, let $K$ and $H$ be the transition kernel and Hamiltonian of either Algorithm \ref{alg:isotropic_HMC} or Algorithm \ref{alg:Riemannian_algorithm} with these parameters, and let $\mu_{H}(p,q) \propto e^{-H(p,q)}$ be the associated Hamiltonian measure. Let $S \subset \mathbb{R}^{d}$ be any subset whose boundary $\partial S$ is a smooth manifold satisfying Assumption \ref{DefLocallyWellBehaved}. Finally, in the case of  Algorithm \ref{alg:Riemannian_algorithm}, let Assumption \ref{AssumptionBoundFIM} also hold.

Then the conductance of $K$ satisfies
\begin{equs}\label{eq:HMC_traditional_Cheeger}
\Phi(K,S) = \Phi^+ \cdot \mathbb{E}\bigg[\frac{1}{N_{\partial S}(P,Q)} \cdot \mathbbm{1}\{N_{\partial S}(P,Q) \, \mathrm{odd}\}\bigg] \bigg/ (\pi(S)),
\end{equs}

where the expectation is taken with respect to the random variables $(P,Q) \sim \mathbb{Q}$ and the total positive flux $\Phi^+$ is given by
\begin{equs}  \label{EqPosFluxConc1}
\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}q.
\end{equs}

In the case of Algorithm \ref{alg:isotropic_HMC}, this formula for $\Phi^+$ reduces to the simpler expression

\begin{equs}  \label{EqPosFluxConc2}
\Phi^+ = \frac{1}{2}T \cdot \int_{\partial S} \pi(q) \mathrm{d}q.
\end{equs} 

\end{thm}

\begin{proof}
The proof of this result is given in Appendix \ref{SecPfThm1}.
\end{proof}



Note that, in Equations \eqref{EqPosFluxConc1} and \eqref{EqPosFluxConc2}, 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 $\mathbb{E}_\mathbb{Q}\bigg[\frac{1}{N_{\partial S}} \cdot \mathbbm{1}\{N_{\partial S} \, \mathrm{odd}\}\bigg] \leq 1$, we have the following useful consequence of Theorem \ref{thm:1}:


\begin{cor}[Simple Conductance Bound] \label{cor:1} 

Set notation as in Theorem \ref{thm:1}. Then the conductance for the Isotropic-Momentum HMC algorithm is bounded by 

\begin{equs} 
\Phi(K,S) \leq \frac{1}{2}T\frac{\int_{\partial S} \pi(q) \mathrm{d}q \cdot }{\pi(S)}. 
\end{equs} 
\end{cor}

An immediate consequence of this bound is that the time it takes energy-conserving Hamiltonian Markov 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 \ref{figure:simulation}) suggest that this upper bound is nearly tight in many cases where $T$ is not too large. 


\subsection{Removing Assumption \ref{DefLocallyWellBehaved}}

Theorem \ref{thm:1} only applies to sets that satisfy Assumption \ref{DefLocallyWellBehaved}. Fortunately, it is not necessary to consider any other sets to compute the conductance of typical HMC Markov chains:



\begin{lemma} [Sufficiency of Compact Locally-Good Manifolds] \label{LemmSuffLoc}
Let $\pi$ be a distribution with twice-differentiable density, let $T > 0$, and let $K_{T}$ be the transition kernel given by Algorithm \ref{alg:isotropic_HMC} 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 \ref{DefLocallyWellBehaved}. Then
\begin{equs} \label{IneqSuffLoc1}
\inf_{S \in \mathfrak{A}} \Phi(K_{T},S) = \inf_{S \in \mathfrak{B}} \Phi(K_{T},S).
\end{equs} 
\end{lemma}

\begin{proof}
The proof is deferred to Appendix \ref{SecLemmSuffLoc}.
\end{proof}




\section{Examples and Applications} \label{SecAppl}

We illustrate our conductance bounds with two examples:


\begin{enumerate}
\item In Section \ref{SubsecAppMix}, we use Theorem \ref{thm:1} to compute the Cheeger constant of an HMC algorithm targetting a mixture of Gaussians. We also compute the spectral gap of this Markov chain, showing that it is roughly equal to the Cheeger constant. In particular, the \textit{upper} bound of Inequality \eqref{IneqCheegPoin} is close to sharp in this example.


\item In Section \ref{SubsecAppDeg}, we use Theorem \ref{thm:1} to compute the Cheeger constant of an HMC algorithm targetting a multivariate Gaussian with nearly-singular covariance matrix. We also compute the spectral gap of this Markov chain, showing that it is roughly equal to the \textit{square} of the Cheeger constant.  In particular, the \textit{lower} bound of Inequality \eqref{IneqCheegPoin} is close to sharp.
\item In Section \ref{SubsecAppNum}, we numerically compute the spectral gap and Cheeger constant of an HMC algorithm to illustrate Theorem \ref{thm:1}.
\end{enumerate}


\subsection{Application: HMC Targetting Mixture of Gaussians} \label{SubsecAppMix}



For $\sigma > 0$, define the mixture distribution 
\begin{equs}  \label{EqDefMixDistBasic}
\pi_{\sigma} = \frac{1}{2} \mathcal{N}(-1,\sigma^2) + \frac{1}{2} \mathcal{N}(1,\sigma^2)
\end{equs}  
and denote its density by $f_{\sigma}$. Let $K_{\sigma}$ be the transition kernel of Algorithm \ref{alg:isotropic_HMC} 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 Cheeger constant 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 
\begin{equs} \label{EqDefTauSigmaX}
\tau_{x}^{(\sigma)} = \inf \{ t > 0 \, : \, X_{t}^{(\sigma)} \notin (-\infty,0)\}.
\end{equs}


We show:

\begin{thm} [HMC for Multimodal Distributions] \label{ThmHmcMultimodal}
The Cheeger constant $\Phi_{\sigma}$ satisfies 
\begin{equs}  \label{EqMultiCheegAsym}
\lim_{\sigma \rightarrow 0} (-2 \sigma^{2}) \, \log(\Phi_{\sigma}) = 1.
\end{equs}
Furthermore, for all $\epsilon > 0$ and fixed $x \in (-\infty,0)$, the hitting time $\tau_{x}^{(\sigma)}$ satisfies 
\begin{equs} \label{EqMultiCheegHitting}
\lim_{\sigma \rightarrow 0} \mathbb{P}[\frac{\log(\tau_{x}^{(\sigma)})}{\log(\Phi_{\sigma})} < 1 + \epsilon] = 1
\end{equs} 
and the relaxation time satisfies 
\begin{equs}   \label{IneqRelMulti}
\lim_{\sigma \rightarrow 0} \frac{\log(\lambda_{\sigma})}{\log(\Phi_{\sigma})} = \lim_{\sigma \rightarrow 0} \frac{\log(\lambda_{\sigma})}{\log(\Phi(K_{\sigma}))} = 1.
\end{equs} 
\end{thm}

We defer the proof to Appendix \ref{AppSubsecPfHmcMulti}. Note that this result exactly matches the spectral gap for optimally-tuned random-walk Metropolis algorithm with the same target distribution (see Theorem 3 of the our companion paper \cite{mangoubi2018simple}).


We also observe that this result implies the Cheeger constant $\Phi(K_{\sigma})$ of $K_{\sigma}$ is close to the bottleneck ratio $\Phi_{\sigma} = \Phi(K_{\sigma}, (-\infty,0))$ associated with the 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. In fact, deriving simple criteria to verify this fact was the motivation behind our companion paper \cite{mangoubi2018simple}.


\subsection{High Dimensional Analog} \label{sec:highdim} Almost all of the work in the proof of Theorem \ref{ThmHmcMultimodal} is checking that, in fact, the spectral gap is not much \textit{smaller} than the natural guess and upper bound $\Phi(K_{\sigma}, (-\infty,0))$; computing a sharp upper bound $\Phi(K_{\sigma}, (-\infty,0))$ itself is straightforward using Theorem \ref{thm:1}. This remains true in higher-dimensional examples. For example, considering the mixture distribution
\begin{equs} 
\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{equs}  
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 \ref{thm:1} can be used to get a very good estimate of this conductance, using the same calculation as the start of Theorem \ref{ThmHmcMultimodal}.

\subsection{Choice of Integration Time and Computational Cost}\label{sec:RemChoiceInt}
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  \cite{bornn2017use,sherlock2017pseudo} and references therein). It is this effective computational cost that is bounded below: we always have
\begin{equs} 
\limsup_{\sigma \rightarrow 0} 2 \sigma^{2} (\log(\lambda_{\sigma}) - \log(T_{\sigma})) \leq 1,
\end{equs} for all $T_{\sigma} \geq \sigma$ (see discussion in Section \ref{SecHmcImp} 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."\footnote{The first use of Goldilocks analogy in this context is by Jeff Rosenthal.} For HMC, one chooses $T_{\sigma}$  so that the probability of an HMC trajectory making a ``U-turn" (in the sense of \cite{NUTS}) is not too close to 0 and not too close to 1. In this example, this means choosig $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 \cite{mangoubi2018simple}, even though it is \textit{not} close to optimal in that context.



\subsection{Application: HMC Targetting Degenerate Multivariate Gaussian} \label{SubsecAppDeg}

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

\begin{equs} 
   M_{\sigma}=
  \left[ {\begin{array}{cc}
   1 & 0 \\
   0 & \sigma^{2} \\
  \end{array} } \right].
\end{equs}


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 \ref{alg:isotropic_HMC} targetting $\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 \ref{alg:isotropic_HMC}, with $\{X_{t}[1]\}_{t \in \mathbb{N}}$ targetting $\mathcal{N}(0,1)$ and $\{X_{t}[2]\}_{t \in \mathbb{N}}$ targetting $\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 Cheeger constant associated with this transition kernel and the set $(-\infty,0)$; denote by $\rho_{\sigma}$ the spectral gap of $K_{\sigma}$. We have:



\begin{thm} \label{ThmHmcDegenerate}
The Cheeger constant $\Phi_{\sigma}$ satisfies 
\begin{equs}  \label{EqMultiCheegDeg}
\lim_{\sigma \rightarrow 0} \frac{\log(\Phi_{\sigma})}{\log(T_{\sigma})} \leq 1.
\end{equs} 
Furthermore, the relaxation time satisfies 
\begin{equs}  \label{IneqRelDeg}
\limsup_{\sigma \rightarrow 0} \frac{\log(\rho_{\sigma})}{\log(T_{\sigma})} \geq \frac{1}{2}.
\end{equs}
\end{thm}

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. 

\begin{remark}
This example is studied in Example 2.1 of \cite{rabee2018HMCcoup}. We observe that our results are sharper in a few ways (we allow larger integration times; we obtain the correct order of the dependence of the spectral gap on $\sigma$), but are also the results of exact computations rather than general theorems.

We mention that our own previous work \cite{mangoubi2017concave} and that of \cite{holmes2014curvature} gives similar (non-sharp) estimates to \cite{rabee2018HMCcoup} in this example. We think that finding general bounds that give sharp answers in this prototypical case to be an interesting problem.
\end{remark}

\subsection{Numerical Simulation: Spectral Gap of Hamiltonian Monte Carlo} \label{SubsecAppNum}



In this section we plot the spectral gap (Figure \ref{figure:simulation}) of an Isotropic-Momentum HMC algorithm sampling 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 \ref{cor:1} with Inequality \eqref{IneqCheegPoin}.

In this simulation we computed the spectral gap for the Isotropic-Momentum HMC algorithm with the stationary distributions of the form $\pi_{a}(q) = \frac{1}{2F_{\mathcal{N}(0,1)}(a)}\mathrm{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 CDF and PDF of the standard normal density.  Note that $\pi_{a}(q)$ approximates the Gaussian mixture model $\tilde{\pi}_{a}(q) =  \frac{1}{2}f_{\mathcal{N}(0,1)}(q-a) + \frac{1}{2}f_{\mathcal{N}(0,1)}(q+a)$; indeed $\lim_{a \rightarrow \infty} \| \pi_{a} - \tilde{\pi}_{a} \|_{\mathrm{TV}} = 0$.

As suggested by the formula for the conductance in Theorem \ref{thm: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 \frac{\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\frac{\pi}{2}$, meaning that the expectation term in Equation \eqref{eq:HMC_traditional_Cheeger} 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}q = f_{\mathcal{N}(0,1)}(q-a)$, so the corresponding term in Corollary \ref{cor:1} of Theorem \ref{thm:1}  decreases exponentially in $a^2$.

\begin{figure}[H]
\includegraphics[scale=0.3]{Spectral_gap}
\caption{The spectral gap for the Isotropic-Momentum HMC algorithm with stationary distribution $\pi(q) = \frac{1}{2F_{\mathcal{N}(0,1)}(a)}\mathrm{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 bound in Theorem \ref{thm:1}.\label{figure:simulation}}
\end{figure}

\section{Discussion and Open Problems} \label{SecCompMulti} 

We give some consequences of Theorem \ref{thm:1}, and discuss open questions related to the performance of HMC.


\subsection{Bounds on HMC Improvements} \label{SecHmcImp}

One immediate consequence of Theorem \ref{thm:1} and Lemma \ref{LemmSuffLoc} is that it is not possible to dramatically improve the performance of Algorithm \ref{alg:isotropic_HMC} 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 twice-differentiable density $\pi$ and denote by $K_{T}$ the transition kernel defined by Algorithm \ref{alg:isotropic_HMC} with parameters $\pi$ and $T > 0$. Define 
\begin{equs} 
\Phi_{0}(\pi,S) &= \frac{1}{2} \frac{\int_{\partial S} \pi(q) \mathrm{d}q  }{\pi(S)}\\
\Phi_{0}(\pi) &= \inf_{S \in \mathfrak{A}}  \Phi_{0}(\pi,S), \\
\end{equs} 
for $S \in \mathfrak{A}$. By Corollary \ref{cor:1} of Theorem \ref{thm:1}, combined with Lemma \ref{LemmSuffLoc}, we have the upper bound
\begin{equs} \label{IneqLinIncCond}
T^{-1} \Phi(K_{T}) \leq  \Phi_{0}(\pi),
\end{equs} 
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 \ref{alg:RWM}, whose tuning parameter $\epsilon > 0$ can have very large impacts 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 \ref{alg:RWM} with target distribution $\pi_{\sigma}$ given in Equation \eqref{EqDefMixDistBasic} and tuning parameter $\epsilon \equiv 10$ for all $\sigma$. It is a straightforward exercise \footnote{ For example, one could apply Theorem 5 of \cite{rosenthal1995minorization} with ``small set" $(-5,5)$ and ``Lyapunov function" $V(x) = e^{\|x\|}$ to obtain a bound on the rate of convergence of the walk to stationarity. One could then apply Theorem 2.1 of \cite{roberts1997geometric} to convert this bound on the convergence rate to a bound on the spectral gap.} to check that $\log(\rho(Q_{\sigma})^{-1}) = O(\log(\sigma^{-1}))$ for $\sigma$ small. Comparing this with the result of Theorem 3 of \cite{mangoubi2018simple}, we see that changing the tuning parameter $\epsilon$ of Algorithm \ref{alg:RWM} 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 mode 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 \eqref{IneqLinIncCond} in terms of computational cost. Roughly speaking, the computational cost of a step of Algorithm \ref{alg:isotropic_HMC} is  proportional to the integration time $T$. Paraphrasing Inequality \eqref{IneqLinIncCond}, it is impossible to improve the cost-normalized conductance of Algorithm \ref{alg:isotropic_HMC} by increasing the integration time $T$. This immediately implies, via Cheeger's inequality \eqref{IneqCheegPoin}, 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 \cite{chen1999lifting} 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 \cite{chen1999lifting,ramanan2017bounds} 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 \textit{of a given set} can increase at most \textit{linearly} in $T$. Via Cheeger's inequality \eqref{IneqCheegPoin}, this suggests that the spectral gap should satisfy a similar \textit{quadratic} inequality in $T$. To make this rigorous, it would be sufficient to show:

\textbf{Open Problem 1:} Prove the limit
\begin{equs} \label{EqLimOpenQuestion}
\lim_{T \rightarrow 0} \inf_{S \in \mathfrak{A}} T^{-1} \Phi(K_{T},S) = \inf_{S \in \mathfrak{A}} \Phi_{0}(\pi,S). 
\end{equs} 


Note that Equality \eqref{IneqE4Small} already implies that, for fixed $S \in \mathfrak{A}$, 
\begin{equs} 
\lim_{T \rightarrow 0}  T^{-1} \Phi(K_{T},S) =  \Phi_{0}(\pi,S). 
\end{equs} 


\subsection{Multimodal Targets}



In Theorem \ref{ThmHmcMultimodal} and Theorem 3 of the companion paper \cite{mangoubi2018simple}, we showed that Algorithm \ref{alg:isotropic_HMC} has very similar performance to Algorithm \ref{alg:RWM} 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:

\textbf{Open Question 2:} Let $\pi_{\sigma}$ be a mixture distribution of the form 
\begin{equs} 
\pi_{\sigma}(x) \propto \sum_{i=1}^{k} \mu_{i} f_{i}(\frac{x}{\sigma}),
\end{equs} 
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$, $a_{1},\ldots,a_{k}$.

Let $K_{\sigma}$, $Q_{\sigma}$ be the transition kernels of Algorithms \ref{alg:isotropic_HMC} and \ref{alg:RWM} 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 
\begin{equs} \label{LimConj2Strong}
\lim_{\sigma \rightarrow 0} \frac{\log(\rho(K_{\sigma}))}{\log(\rho(Q_{\sigma}))} \geq 1?
\end{equs} 

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

We conjecture that Inequality \eqref{LimConj2Strong} 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{equs} 
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{equs} 
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 e^{\sigma^{-1}}$ look like Figure \ref{FigBanana},






More generally, we are interested in determining which ``momentum-based" algorithms satisfy the limit \eqref{LimConj2Strong}. We note that not all similar chains satisfy this equality. In particular, lifted Markov chains can replace the ``1" with a ``$\frac{1}{2}$" for realistic examples (see \textit{e.g.} \cite{bierkens2017piecewise}). It would be valuable to determine if other generic algorithms, such as Algorithm \ref{alg:Riemannian_algorithm} or \cite{bierkens2016zig}, can also achieve this.

Finally, we suggest a more straightforward question. The level sets of the multimodal densities studied in this paper are \textit{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 MH for multimodal examples in which the level sets are \textit{connected}, and multimodality is induced by these level sets being merely \textit{very narrow}? Note that this is only possible in two or more dimensions.

\section*{Acknowledgement}
We would like to thank Neil Shephard for asking us the title question.
NSP thanks Gareth Roberts and Andrew Stuart for introducing him to the optimal scaling of MCMC algorithms, and Jeff Rosenthal for references. Part of this work was done when OM was a Ph.D. student at MIT working on this project under the supervision of NSP at Harvard University, and later when he was a postdoctoral researcher with AS at the University of Ottawa. 
We thank these institutions for their hospitality.


\bibliographystyle{alpha}

2018-09-05T02:32:46.000
1808
1808.03230
en
https://arxiv.org/abs/1808.03230
Hamiltonian Monte Carlo (HMC) is a very popular and generic collection of Markov chain Monte Carlo (MCMC) algorithms. One explanation for the popularity of HMC algorithms is their excellent performance as the dimension $d$ of the target becomes large: under conditions that are satisfied for many common statistical models, optimally-tuned HMC algorithms have a running time that scales like $d^{0.25}$. In stark contrast, the running time of the usual Random-Walk Metropolis (RWM) algorithm, optimally tuned, scales like $d$. This superior scaling of the HMC algorithm with dimension is attributed to the fact that it, unlike RWM, incorporates the gradient information in the proposal distribution. In this paper, we investigate a different scaling question: does HMC beat RWM for highly $\textit{multimodal}$ targets? We find that the answer is often $\textit{no}$. We compute the spectral gaps for both the algorithms for a specific class of multimodal target densities, and show that they are identical. The key reason is that, within one mode, the gradient is effectively ignorant about other modes, thus negating the advantage the HMC algorithm enjoys in unimodal targets. We also give heuristic arguments suggesting that the above observation may hold quite generally. Our main tool for answering this question is a novel simple formula for the conductance of HMC using Liouville's theorem. This result allows us to compute the spectral gap of HMC algorithms, for both the classical HMC with isotropic momentum and the recent Riemannian HMC, for multimodal targets.
Probability (math.PR); Machine Learning (cs.LG); Computation (stat.CO); Methodology (stat.ME); Machine Learning (stat.ML)
Does Hamiltonian Monte Carlo mix faster than a random walk on multimodal densities?
Qwen/Qwen-72B
1. YES
2. YES