.tocmtchapter mtchaptersubsection mtappendixnone

# Fast mixing of Metropolized Hamiltonian Monte Carlo: Benefits of multi-step
gradients

Yuansi Chen<EMAIL_ADDRESS>
Department of Statistics
University of California
Berkeley, CA 94720-1776, USA Raaz Dwivedi<EMAIL_ADDRESS>
Department of Electrical Engineering and Computer Sciences
University of California
Berkeley, CA 94720-1776, USA Martin J. Wainwright<EMAIL_ADDRESS>
Bin Yu<EMAIL_ADDRESS>
Department of Statistics
Department of Electrical Engineering and Computer Sciences
University of California
Berkeley, CA 94720-1776, USA

###### Abstract

Hamiltonian Monte Carlo (HMC) is a state-of-the-art Markov chain Monte Carlo
sampling algorithm for drawing samples from smooth probability densities over
continuous spaces. We study the variant most widely used in practice,
Metropolized HMC with the Störmer-Verlet or leapfrog integrator, and make two
primary contributions. First, we provide a non-asymptotic upper bound on the
mixing time of the Metropolized HMC with explicit choices of step-size and
number of leapfrog steps. This bound gives a precise quantification of the
faster convergence of Metropolized HMC relative to simpler MCMC algorithms
such as the Metropolized random walk, or Metropolized Langevin algorithm.
Second, we provide a general framework for sharpening mixing time bounds of
Markov chains initialized at a substantial distance from the target
distribution over continuous spaces. We apply this sharpening device to the
Metropolized random walk and Langevin algorithms, thereby obtaining improved
mixing time bounds from a non-warm initial distribution.

##  1 Introduction

Markov Chain Monte Carlo (MCMC) methods date back to the seminal work of
Metropolis et al. (1953), and are the method of choice for drawing samples
from high-dimensional distributions. They are widely used in practice,
including in Bayesian statistics for exploring posterior distributions
(Carpenter et al., 2017; Smith, 2014), in simulation-based methods for
reinforcement learning, and in image synthesis in computer vision, among other
areas. Since their origins in the 1950s, many MCMC algorithms have been
introduced, applied and studied; we refer the reader to the handbook by Brooks
et al. (2011) for a survey of known results and contemporary developments.

There are a variety of MCMC methods for sampling from target distributions
with smooth densities (Robert and Casella, 1999; Roberts et al., 2004; Roberts
and Stramer, 2002; Brooks et al., 2011). Among them, the method of Hamiltonian
Monte Carlo (HMC) stands out among practitioners: it is the default sampler
for sampling from complex distributions in many popular software packages,
including Stan (Carpenter et al., 2017), Mamba (Smith, 2014), and Tensorflow
(Abadi et al., 2015). We refer the reader to the papers (Neal, 2011; Hoffman
and Gelman, 2014; Durmus et al., 2017) for further examples and discussion of
the HMC method. There are a number of variants of HMC, but the most popular
choice involves a combination of the leapfrog integrator with Metropolis-
Hastings correction. Throughout this paper, we reserve the terminology HMC to
refer to this particular Metropolized algorithm. The idea of using Hamiltonian
dynamics in simulation first appeared in Alder and Wainwright (1959). Duane et
al. (1987) introduced MCMC with Hamiltonian dynamics, and referred to it as
Hybrid Monte Carlo. The algorithm was further refined by Neal (1994), and
later re-christened in the statistics community as Hamiltonian Monte Carlo. We
refer the reader to Neal (2011) for an illuminating overview of the history of
HMC and a discussion of contemporary work.

###  1.1 Past work on HMC

While HMC enjoys fast convergence in practice, a theoretical understanding of
this behavior remains incomplete. Some intuitive explanations are based on its
ability to maintain a constant asymptotic accept-reject rate with large step-
size (Creutz, 1988). Others suggest, based on intuition from the continuous-
time limit of the Hamiltonian dynamics, that HMC can suppress random walk
behavior using momentum (Neal, 2011). However, these intuitive arguments do
not provide rigorous or quantitative justification for the fast convergence of
the discrete-time HMC used in practice.

More recently, general asymptotic conditions under which HMC will or will not
be geometrically ergodic have been established in some recent papers (Durmus
et al., 2017; Livingstone et al., 2016). Other work has yielded some insight
into the mixing properties of different variants of HMC, but it has focused
mainly on _unadjusted_ versions of the algorithm. Mangoubi and Smith (2017)
and Mangoubi and Vishnoi (2018) study versions of unadjusted HMC based on
Euler discretization or leapfrog integrator (but omitting the Metropolis-
Hastings step), and provide explicit bounds on the mixing time as a function
of dimension $d$, condition number $\kappa$ and error tolerance $\epsilon>0$.
Lee et al. (2018) studied an extended version of HMC that involves applying an
ordinary differential equation (ODE) solver; they established bounds with
sublinear dimension dependence, and even polylogarithmic for certain densities
(e.g., those arising in Bayesian logistic regression). The mixing time for the
same algorithm is further refined in recent work by Chen and Vempala (2019).
In a similar spirit, Lee and Vempala (2018a) studied the Riemannian variant of
HMC (RHMC) with an ODE solver focusing on sampling uniformly from a polytope.
While their result could be extended to log-concave sampling, the practical
implementation of their ODE solver for log-concave sampling is unclear.
Moreover, it requires a regularity condition on all the derivatives of
density. One should note that such unadjusted HMC methods behave differently
from the Metropolized version most commonly used in practice. In the absence
of the Metropolis-Hastings correction, the resulting Markov chain no longer
converges to the correct target distribution, but instead exhibits a
persistent bias, even in the limit of infinite iterations. Consequently, the
analysis of such sampling methods requires controlling this bias; doing so
leads to mixing times that scale polynomially in $1/\epsilon$, in sharp
contrast with the $\log(1/\epsilon)$ that is typical for Metropolis-Hastings
corrected methods.

Most closely related to our paper is the recent work by Bou-Rabee et al.
(2018), which studies the same Metropolized HMC algorithm that we analyze in
this paper. They use coupling methods to analyze HMC for a class of
distributions that are strongly log-concave outside of a compact set. In the
strongly log-concave case, they prove a mixing-time bound that scales at least
as $d^{3/2}$ in the dimension $d$. We note that with a “warm” initialization,
this dimension dependence grows more quickly than known bounds for the MALA
algorithm (Dwivedi et al., 2019; Eberle, 2014), and so does not explain the
superiority of HMC in practice.

In practice, it is known that Metropolized HMC is fairly sensitive to the
choice of its parameters—namely the step-size $\eta$ used in the
discretization scheme, and the number of leapfrog steps $K$. At one extreme,
taking a single leapfrog step $K=1$, the algorithm reduces to the Metropolis
adjusted Langevin algorithm (MALA). More generally, if too few leapfrog steps
are taken, HMC is likely to exhibit a random walk behavior similar to that of
MALA. At the other extreme, if $K$ is too large, the leapfrog steps tend to
wander back to a neighborhood of the initial state, which leads to wasted
computation as well as slower mixing (Betancourt et al., 2014). In terms of
the step size $\eta$, when overly large step size is chosen, the
discretization diverges from the underlying continuous dynamics leading to a
drop in the Metropolis acceptance probability, thereby slowing down the mixing
of the algorithm. On the other hand, an overly small choice of $\eta$ does not
allow the algorithm to explore the state space rapidly enough. While it is
difficult to characterize the necessary and sufficient conditions on $K$ and
$\eta$ to ensure fast convergence, many works suggest the choice of these two
parameters based on the necessary conditions such as maintaining a constant
acceptance rate (Chen et al., 2001). For instance, Beskos et al. (2013) showed
that in the simplified scenario of target density with independent,
identically distributed components, the number of leapfrog steps should scale
as $d^{1/4}$ to achieve a constant acceptance rate. Besides, instead of
setting the two parameters explicitly, various automatic strategies for tuning
these two parameters have been proposed (Wang et al., 2013; Hoffman and
Gelman, 2014; Wu et al., 2018). Despite being introduced via heuristic
arguments and with additional computational cost, these methods, such as the
No-U-Turn (NUTS) sampler (Hoffman and Gelman, 2014), have shown promising
empirical evidence of its effectiveness on a wide range of simple target
distributions.

###  1.2 Past work on mixing time dependency on initialization

Many proof techniques for the convergence of continuous-state Markov chains
are inspired by the large body of work on discrete-state Markov chains; for
instance, see the surveys (Lovász et al., 1993; Aldous and Fill, 2002) and
references therein. Historically, much work has been devoted to improving the
mixing time dependency on the initial distribution. For discrete-state Markov
chains, Diaconis and Saloff-Coste (1996) were the first to show that the
logarithmic dependency of the mixing time of a Markov chain on the warmness
parameter111See equation (4) for a formal definition. of the starting
distribution can be improved to doubly logarithmic. This improvement—from
logarithmic to doubly logarithmic—allows for a good bound on the mixing time
even when a good starting distribution is not available. The innovation
underlying this improvement is the use of log-Sobolev inequalities in place of
the usual isoperimetric inequality. Later, closely related ideas such as
average conductance (Lovász and Kannan, 1999; Kannan et al., 2006), evolving
sets (Morris and Peres, 2005) and spectral profile (Goel et al., 2006) were
shown to be effective for reducing dependence on initial conditions for
discrete space chains. Thus far, only the notion of average conductance
(Lovász and Kannan, 1999; Kannan et al., 2006) has been adapted to continuous-
state Markov chains so as to sharpen the mixing time analysis of the Ball walk
(Lovász and Simonovits, 1990).

###  1.3 Our contributions

This paper makes two primary contributions. First, we provide a non-asymptotic
upper bound on the mixing time of the Metropolized HMC algorithm for smooth
densities (see Theorem 1). This theorem applies to the form of Metropolized
HMC (based on the leapfrog integrator) that is most widely used in practice.
To the best of our knowledge, Theorem 1 is the first rigorous confirmation of
the faster non-asymptotic convergence of the Metropolized HMC as compared to
MALA and other simpler Metropolized algorithms.222As noted earlier, previous
results by Bou-Rabee et al. (2018) on Metropolized HMC do not establish that
it mixes more rapidly than MALA. Other related works on HMC consider either
its unadjusted version (without accept-reject step) with different integrators
(Mangoubi and Smith, 2017; Mangoubi and Vishnoi, 2018) or the HMC based on an
ODE solver (Lee et al., 2018; Lee and Vempala, 2018a). While the dimension
dependency for these algorithms is usually better than MALA, they have
polynomial dependence on the inverse error tolerance $1/\epsilon$ while MALA’s
mixing time scales as $\log(1/\epsilon)$. Moreover, our direct analysis of the
Metropolized HMC with a leapfrog integrator provides explicit choices of the
hyper-parameters for the sampler, namely, the step-size and the number of
leapfrog updates in each step. Our theoretical choices of the hyper-parameters
could potentially provide guidelines for parameter tuning in practical HMC
implementations

Our second main contribution is formalized in Lemmas 3 and 4: we develop
results based on the conductance profile in order to prove quantitative
convergence guarantees general continuous-state space Markov chains. Doing so
involves non-trivial extensions of ideas from discrete-state Markov chains to
those in continuous state spaces. Our results not only enable us to establish
the mixing time bounds for HMC with different classes of target distributions,
but also allow simultaneous improvements on mixing time bounds of several
Markov chains (for general continuous-state space) when the starting
distribution is far from the stationary distribution. Consequentially, we
improve upon previous mixing time bounds for Metropolized Random Walk (MRW)
and MALA (Dwivedi et al., 2019), when the starting distribution is not _warm_
with respect to the target distribution (see Theorem 5).

While this high-level road map is clear, a number of technical challenges
arise en route in particular in controlling the conductance profile of HMC.
The use of multiple gradient steps in each iteration of HMC helps it mix
faster but also complicates the analysis; in particular, a key step is to
control the overlap between the transition distributions of HMC chain at two
nearby points; doing so requires a delicate argument (see Lemma 6 and Section
5.3 for further details).

Table 1 provides an informal summary of our mixing time bounds of HMC and how
they compare with known bounds for MALA when applied to log-concave target
distributions. From the table, we see that Metropolized HMC takes fewer
gradient evaluations than MALA to mix to the same accuracy for log-concave
distributions. Note that our current analysis establishes logarithmic
dependence on the target error $\epsilon$ for strongly-log-concave as well as
for a sub-class of weakly log-concave distributions.333For comparison with
previous results on unadjusted HMC or ODE based HMC refer to the discussion
after Corollary 2 and Table 7 in Appendix D.2. Moreover, in Figure 1 we
provide a sketch-diagram of how the different results in this paper interact
to provide mixing time bounds for different Markov chains.

|  Strongly log-concave |  Weakly log-concave
---|---|---
Sampling algorithm |  Assumption (B) ($\kappa\ll d$) |  Assumption (C) |  Assumption (D)
MALA (improved bound in Thm 5 in this paper) |  $\displaystyle d\kappa\log\frac{1}{\epsilon}$ Dwivedi et al. (2019) |  $\displaystyle\frac{d^{2}}{\epsilon^{\frac{3}{2}}}\log\frac{1}{\epsilon}$ Dwivedi et al. (2019) |  $\displaystyle d^{\frac{3}{2}}\log\frac{1}{\epsilon}$ Mangoubi and Vishnoi (2019)
Metropolized HMC with leapfrog integrator [this paper] |  $\displaystyle d^{\frac{11}{12}}\kappa\log\frac{1}{\epsilon}$ (Corollary 2) |  $\displaystyle\frac{d^{\frac{11}{6}}}{\epsilon}\log\frac{1}{\epsilon}$ (Corollary 18) |  $\displaystyle d^{\frac{4}{3}}\log\frac{1}{\epsilon}$ (Corollary 18)
Table 1: Comparisons of the number of gradient evaluations needed by MALA and
Metropolized HMC with leapfrog integrator from a _warm start_ to obtain an
$\epsilon$-accurate sample in TV distance from a log-concave target
distribution on d. The second column corresponds to strongly log-concave
densities with condition number $\kappa$, and the third and fourth column
correspond to weakly log-concave densities satisfying certain regularity
conditions. Theorem 1
---
general mixing-time
result for HMC
Lemma 3
---
mixing-time bound
via conductance profile
Lemma 4
---
conductance profile
via overlap bounds
Lemma 6
---
overlap bounds
for HMC
known overlap bounds
---
for MALA and MRW
Dwivedi et al. (2019)
Corollaries 2 and 18
---
HMC mixing-time for
different target distributions
Theorem 5
---
improved mixing-time
for MALA and MRW
Figure 1: High-level sketch diagram of the results and the proof techniques
developed in this paper. Lemmas 3 and 4 of this paper provide general
machinery to develop mixing time bounds for several Markov chains, provided
that we have a control over the overlap bounds of the respective transition
distributions. Establishing these overlap bounds for Metropolized HMC requires
non-trivial machinery which we develop in Lemma 6. Putting Lemmas 3, 4 and 6
together yields the general mixing time bound for HMC in Theorem 1. Moreover,
Theorem 1 easily implies a faster mixing time bound for HMC over MALA and MRW
for different target distributions as shown in Corollaries 2 and 18, and
summarized in Table 1. On the other hand, Lemmas 3 and 4 in conjunction with
overlap bounds from Dwivedi et al. (2019) immediately yield Theorem 5 that
provides improved mixing time bounds for MALA and MRW from a non-warm start.

##### Organization:

The remainder of the paper is organized as follows. Section 2 is devoted to
background on the idea of Monte Carlo approximation, Markov chains and MCMC
algorithms, and the introduction of the MRW, MALA and HMC algorithms. Section
3 contains our main results on mixing time of HMC in Section 3.2, followed by
the general framework for obtaining sharper mixing time bounds in Section 3.3
and its application to MALA and MRW in Section 3.4. In Section 4, we describe
some numerical experiments that we performed to explore the sharpness of our
theoretical predictions in some simple scenarios. In Section 5, we prove
Theorem 1 and Corollary 14, with the proofs of technical lemmas and other
results deferred to the appendices. We conclude in Section 6 with a discussion
of our results and future directions.

##### Notation:

For two real-valued sequences $\\{a_{n}\\}$ and $\\{b_{n}\\}$, we write
$a_{n}=O(b_{n})$ if there exists a universal constant $c>0$ such that
$a_{n}\leq cb_{n}$. We write $a_{n}=\widetilde{O}(b_{n})$ if $a_{n}\leq
c_{n}b_{n}$, where $c_{n}$ grows at most poly-logarithmically in $n$. We use
$[K]$ to denote the integers from the set $\left\\{1,2,\ldots,K\right\\}$. We
denote the Euclidean norm on d as $\left\|\,\cdot\,\right\|_{2}$. We use
$\mathcal{X}$ to denote the (general) state space of a Markov chain. We denote
$\mathcal{B}(\mathcal{X})$ as the Borel $\sigma$-algebra of the state space
$\mathcal{X}$. Throughout we use the notation $c$, $c_{1}$, $c_{2}$ to denote
universal constants. For a function $f:{}^{d}\rightarrow\real$ that is three
times differentiable, we represent its derivatives at $x\in{}^{d}$ by ${\nabla
f}(x)\in{}^{d}$, ${\nabla^{2}f}(x)\in{}^{d\times d}$ and
${\nabla^{3}f}(x)\in{}^{d^{3}}$. Here

$\displaystyle\left[{\nabla f}(x)\right]_{i}=\frac{\partial}{\partial
x_{i}}f(x),\quad\left[{\nabla^{2}f}(x)\right]_{i,j}=\frac{\partial^{2}}{\partial
x_{i}\partial
x_{j}}f(x),\quad\left[{\nabla^{3}f}\right]_{i,j,k}=\frac{\partial^{3}}{\partial
x_{i}\partial x_{j}\partial x_{k}}f(x).$

Moreover for a square matrix $A$, we define its $\ell_{2}$-operator norm
$\left|\\!\left|\\!\left|A\right|\\!\right|\\!\right|_{{\mbox{\tiny{op}}}}:\,=\max\limits_{\left\|v\right\|_{2}=1}\left\|Av\right\|_{2}$.

##  2 Background and problem set-up

In this section, we begin by introducing background on Markov chain Monte
Carlo in Section 2.1, followed by definitions and terminology for Markov
chains in Section 2.2. In Section 2.3, we describe several MCMC algorithms,
including the Metropolized random walk (MRW), the Metropolis-adjusted Langevin
algorithm (MALA), and the Metropolis-adjusted Hamiltonian Monte Carlo (HMC)
algorithm. Readers familiar with the literature may skip directly to the
Section 3, where we set up and state our main results.

###  2.1 Monte Carlo Markov chain methods

Consider a distribution $\Pi^{*}$ equipped with a density
$\pi^{*}:\mathcal{X}\rightarrow{}_{+}$, specified explicitly up to a
normalization constant as follows

$\displaystyle\pi^{*}(x)\propto e^{-f(x)}.$ (1)

A standard computational task is to estimate the expectation of some function
$g:\mathcal{X}\rightarrow\real$—that is, to approximate
$\Pi^{*}(g)={\mathbb{E}}_{\pi^{*}}\left[g(X)\right]=\int_{\mathcal{X}}g(x)\pi^{*}(x)dx$.
In general, analytical computation of this integral is infeasible. In high
dimensions, numerical integration is not feasible either, due to the well-
known curse of dimensionality.

A Monte Carlo approximation to $\Pi^{*}(g)$ is based on access to a sampling
algorithm that can generate i.i.d. random variables $Z_{i}\sim\pi^{*}$ for
$i=1,\ldots,N$. Given such samples, the random variable
$\widehat{\Pi^{*}}(g):\,=\frac{1}{N}\sum_{i=1}^{N}g(Z_{i})$ is an unbiased
estimate of the quantity $\Pi^{*}(g)$, and has its variance proportional to
$1/N$. The challenge of implementing such a method is drawing the i.i.d.
samples $Z_{i}$. If $\pi^{*}$ has a complicated form and the dimension $d$ is
large, it is difficult to generate i.i.d. samples from $\pi^{*}$. For example,
rejection sampling (Gilks and Wild, 1992), which works well in low dimensions,
fails due to the curse of dimensionality.

The Markov chain Monte Carlo (MCMC) approach is to construct a Markov chain on
$\mathcal{X}$ that starts from some easy-to-simulate initial distribution
$\mu_{0}$, and converges to $\pi^{*}$ as its stationary distribution. Two
natural questions are: (i) methods for designing such Markov chains; and (ii)
how many steps will the Markov chain take to converge close enough to the
stationary distribution? Over the years, these questions have been the subject
of considerable research; for instance, see the reviews by Tierney (1994);
Smith and Roberts (1993); Roberts et al. (2004) and references therein. In
this paper, we are particularly interested in comparing three popular
Metropolis-Hastings adjusted Markov chains sampling algorithms (MRW, MALA,
HMC). Our primary goal is to tackle the second question for HMC, in particular
via establishing its concrete non-asymptotic mixing-time bound and thereby
characterizing how HMC converges faster than MRW and MALA.

###  2.2 Markov chain basics

Let us now set up some basic notation and definitions on Markov chains that we
use in the sequel. We consider _time-homogeneous_ Markov chains defined on a
measurable state space $\left(\mathcal{X},\mathcal{B}(\mathcal{X})\right)$
with a transition kernel
$\Theta:\mathcal{X}\times\mathcal{B}(\mathcal{X})\rightarrow{}_{+}$. By
definition, the transition kernel satisfies the following properties:

$\displaystyle\Theta(x,dy)\geq 0,\quad\text{for all
}x\in\mathcal{X},\quad\text{and}\quad\int_{y\in\mathcal{X}}\Theta(x,dy)dy=1\quad\text{for
all }x\in\mathcal{X}.$

The $k$-step transition kernel $\Theta^{k}$ is defined recursively as
$\Theta^{k+1}(x,dy)=\int_{z\in\mathcal{X}}\Theta^{k}(x,dz)\Theta(z,dy)dz$.

The Markov chain is _irreducible_ means that for all $x,y\in\mathcal{X}$,
there is a natural number $k>0$ such that $\Theta^{k}(x,dy)>0$. We say that a
Markov chain satisfies the _detailed balance condition_ if

$\displaystyle\pi^{*}(x)\Theta(x,dy)dx=\pi^{*}(y)\Theta(y,dx)dy\quad\mbox{for
all $x,y\in\mathcal{X}$}.$ (2)

Such a Markov chain is also called _reversible_. Finally, we say that a
probability measure $\Pi^{*}$ with density $\pi^{*}$ on $\mathcal{X}$ is
_stationary_ (or _invariant_) for a Markov chain with the transition kernel
$\Theta$ if

$\displaystyle\int_{x\in\mathcal{X}}\pi^{*}(x)\Theta(y,dx)=\pi^{*}(y)\quad\text{for
all }y\in\mathcal{X}.$

##### Transition operator:

We use $\mathcal{T}$ to denote the transition operator of the Markov chain on
the space of probability measures with state space $\mathcal{X}$. In simple
words, given a distribution $\mu_{0}$ on the current state of the Markov
chain, $\mathcal{T}(\mu_{0})$ denotes the distribution of the next state of
the chain. Mathematically, we have
$\mathcal{T}(\mu_{0})(A)=\int_{\mathcal{X}}\Theta(x,A)\mu_{0}(x)dx$ for any
$A\in\mathcal{B}(\mathcal{X})$. In an analogous fashion, $\mathcal{T}^{k}$
stands for the $k$-step transition operator. We use $\mathcal{T}_{x}$ as the
shorthand for $\mathcal{T}(\mathbf{\delta}_{x})$, the _transition distribution
at $x$_; here $\mathbf{\delta}_{x}$ denotes the Dirac delta distribution at
$x\in\mathcal{X}$. Note that by definition $\mathcal{T}_{x}=\Theta(x,\cdot)$.

##### Distances between two distributions:

In order to quantify the convergence of the Markov chain, we study the mixing
time for a class of distances denoted $\mathcal{L}_{\mathfrak{p}}$ for
$\mathfrak{p}\geq 1$. Letting $Q$ be a distribution with density $q$, its
$\mathcal{L}_{\mathfrak{p}}$-divergence with respect to the positive density
$\nu$ is defined as

$\displaystyle
d_{\mathfrak{p}}(Q,\nu)=\left(\int_{\mathcal{X}}\left|\frac{q(x)}{\nu(x)}-1\right|^{\mathfrak{p}}\nu(x)dx\right)^{\frac{1}{\mathfrak{p}}}.$
(3a) Note that for $\mathfrak{p}=2$, we get the $\chi^{2}$-divergence. For
$\mathfrak{p}=1$, the distance $d_{1}(Q,\nu)$ represents two times the total
variation distance between $Q$ and $\nu$. In order to make this distinction
clear, we use $d_{\text{\tiny{TV}}}\left(Q,\nu\right)$ to denote the total
variation distance.

##### Mixing time of a Markov chain:

Consider a Markov chain with initial distribution $\mu_{0}$, transition
operator $\mathcal{T}$ and a target distribution $\Pi^{*}$ with density
$\pi^{*}$. Its $\mathcal{L}_{\mathfrak{p}}$ mixing time with respect to
$\Pi^{*}$ is defined as follows:

$\displaystyle\tau_{\mathfrak{p}}(\epsilon;\mu_{0})=\inf\left\\{k\in\mathbb{N}\mid
d_{\mathfrak{p}}\left(\mathcal{T}^{k}(\mu_{0}),\Pi^{*}\right)\leq\epsilon\right\\}.$
(3b)

where $\epsilon>0$ is an error tolerance. Since distance
$d_{\mathfrak{p}}(Q,\Pi^{*})$ increases as $\mathfrak{p}$ increases, we have

$\displaystyle\tau_{\mathfrak{p}}(\epsilon;\mu_{0})\leq\tau_{\mathfrak{p}^{\prime}}(\epsilon;\mu_{0})\quad\text{for
any }\quad\mathfrak{p}^{\prime}\geq\mathfrak{p}\geq 1.$ (3c)

##### Warm initial distribution:

We say that a Markov chain with state space $\mathcal{X}$ and stationary
distribution $\Pi^{*}$ has a $\beta$-warm start if its initial distribution
$\mu_{0}$ satisfies

$\displaystyle\sup_{S\in\mathcal{B}(\mathcal{X})}\frac{\mu_{0}(S)}{\Pi^{*}(S)}\leq\beta,$
(4)

where $\mathcal{B}(\mathcal{X})$ denotes the Borel $\sigma$-algebra of the
state space $\mathcal{X}$. For simplicity, we say that $\mu_{0}$ is a warm
start if the warmness parameter $\beta$ is a small constant (e.g., $\beta$
does not scale with dimension $d$).

##### Lazy chain:

We say that the Markov chain is $\zeta$-lazy if, at each iteration, the chain
is forced to stay at the previous iterate with probability $\zeta$. We study
$\frac{1}{2}$-lazy chains in this paper. In practice, one is not likely to use
a lazy chain (since the lazy steps slow down the convergence rate by a
constant factor); rather, it is a convenient assumption for theoretical
analysis of the mixing rate up to constant factors.444Any lazy (time-
reversible) chain is always aperiodic and admits a unique stationary
distribution. For more details, see the survey (Vempala, 2005) and references
therein.

###  2.3 From Metropolized random walk to HMC

In this subsection, we provide a brief description of the popular algorithms
used for sampling from the space $\mathcal{X}={}^{d}$. We start with the
simpler zeroth-order Metropolized random walk (MRW), followed by the single-
step first-order Metropolis adjusted Langevin algorithm (MALA) and finally
discuss the Hamiltonian Monte Carlo (HMC) algorithm.

####  2.3.1 MRW and MALA algorithms

One of the simplest Markov chain algorithms for sampling from a density of the
form (1) defined on d is the Metropolized random walk (MRW). Given state
$x_{i}\in{}^{d}$ at iterate $i$, it generates a new proposal vector
$z_{i+1}\sim\mathcal{N}(x_{i},2\eta\mathbb{I}_{d})$, where $\eta>0$ is a step-
size parameter.555The factor $2$ in the step-size definition is a convenient
notational choice so as to facilitate comparisons with other algorithms. It
then decides to accept or reject $z_{i+1}$ using a Metropolis-Hastings
correction; see Algorithm 1 for the details. Note that the MRW algorithm uses
information about the function $f$ only via querying function values, but not
the gradients.

The Metropolis-adjusted Langevin algorithm (MALA) is a natural extension of
the MRW algorithm: in addition to the function value $f(\cdot)$, it also
assumes access to its gradient $\nabla f(\cdot)$ at any state $x\in{}^{d}$.
Given state $x_{i}$ at iterate $i$, it observes $\left(f(x_{i}),{\nabla
f}(x_{i})\right)$ and then generates a new proposal
$z_{i+1}\sim\mathcal{N}(\displaystyle x_{i}-\eta{\nabla
f}(x_{i}),2\eta\mathbb{I}_{d})$, followed by a suitable Metropolis-Hastings
correction; see Algorithm 2 for the details. The MALA algorithm has an
interesting connection to the Langevin diffusion, a stochastic process whose
evolution is characterized by the stochastic differential equation (SDE)

$\displaystyle dX_{t}=-{\nabla f}(X_{t})dt+\sqrt{2}dW_{t}.$ (5)

The MALA proposal can be understood as the Euler-Maruyama discretization of
the SDE (5).

Input: Step size $\eta>0$ and a sample $x_{0}$ from a starting distribution
$\mu_{0}$

Output: Sequence $x_{1},x_{2},\ldots$

1 for  _$i=0,1,\ldots$_ do

2 Proposal step: _Draw_ $z_{i+1}\sim\mathcal{N}(x_{i},2\eta\mathbb{I}_{d})$

3 Accept-reject step:

4 compute
$\alpha_{i+1}\leftarrow\displaystyle\min\left\\{1,\frac{\exp\left(-f(z_{i+1})\right)}{\exp\left(-f(x_{i})\right)}\right\\}$

5 With probability $\alpha_{i+1}$ _accept_ the proposal: $x_{i+1}\leftarrow
z_{i+1}$

6 With probability $1-\alpha_{i+1}$ _reject_ the proposal: $x_{i+1}\leftarrow
x_{i}$

7 end for

Algorithm 1 Metropolized Random Walk (MRW)

Input: Step size $\eta$ and a sample $x_{0}$ from a starting distribution
$\mu_{0}$

Output: Sequence $x_{1},x_{2},\ldots$

1 for  _$i=0,1,\ldots$_ do

2 Proposal step: _Draw_ $z_{i+1}\sim\mathcal{N}(\displaystyle
x_{i}-\eta{\nabla f}(x_{i}),2\eta\mathbb{I}_{d})$

3 Accept-reject step:

4 compute
$\alpha_{i+1}\leftarrow\displaystyle\min\left\\{1,\frac{\exp\left(-f(z_{i+1})-\left\|x_{i}-z_{i+1}+\eta{\nabla
f}(z_{i+1})\right\|_{2}^{2}/4\eta\right)}{\exp\left(-f(x_{i})-\left\|z_{i+1}-x_{i}+\eta{\nabla
f}(x_{i})\right\|_{2}^{2}/4\eta\right)}\right\\}$

5 With probability $\alpha_{i+1}$ _accept_ the proposal: $x_{i+1}\leftarrow
z_{i+1}$

6 With probability $1-\alpha_{i+1}$ _reject_ the proposal: $x_{i+1}\leftarrow
x_{i}$

7 end for

Algorithm 2 Metropolis adjusted Langevin algorithm (MALA)

####  2.3.2 HMC sampling

The HMC sampling algorithm from the physics literature was introduced to the
statistics literature by Neal; see his survey (Neal, 2011) for the historical
background. The method is inspired by Hamiltonian dynamics, which describe the
evolution of a state vector $q(t)\in{}^{d}$ and its momentum $p(t)\in{}^{d}$
over time $t$ based on a Hamiltonian function
$\mathcal{H}:{}^{d}\times{}^{d}\rightarrow\real$ via Hamilton’s equations:

$\displaystyle\frac{dq}{dt}(t)=\frac{\partial\mathcal{H}}{\partial
p}(p(t),q(t)),\quad\mbox{and}\quad\frac{dp}{dt}(t)=-\frac{\partial\mathcal{H}}{\partial
q}(p(t),q(t)).$ (6)

A straightforward calculation using the chain rule shows that the Hamiltonian
remains invariant under these dynamics—that is, $\mathcal{H}(p(t),q(t))=C$ for
all $t\in\real$. A typical choice of the Hamiltonian
$\mathcal{H}:{}^{d}\times{}^{d}\rightarrow\real$ is given by

$\displaystyle\mathcal{H}(p,q)=f(q)+\frac{1}{2}\left\|p\right\|_{2}^{2}.$ (7)

The ideal HMC algorithm for sampling is based on the continuous Hamiltonian
dynamics; as such, it is not implementable in practice, but instead a useful
algorithm for understanding. For a given time $T>0$ and vectors
$u,v\in{}^{d}$, let $q_{T}(u,v)$ denote the $q$-solution to Hamilton’s
equations at time $T$ and with initial conditions $(p(0),q(0))=(u,v)$. At
iteration $k$, given the current iterate $X_{k}$, the ideal HMC algorithm
generates the next iterate $X_{k+1}$ via the update rule
$X_{k+1}=q_{T}(p_{k},X_{k})$ where $p_{k}\sim N(0,\mathbb{I}_{d})$ is a
standard normal random vector, independent of $X_{k}$ and all past iterates.
It can be shown that with an appropriately chosen $T$, the ideal HMC algorithm
converges to the stationary distribution $\pi^{*}$ without a Metropolis-
Hastings adjustment (see Neal (2011); Mangoubi and Vishnoi (2018) for the
existence of such a solution and its convergence).

However, in practice, it is impossible to compute an exact solution to
Hamilton’s equations. Rather, one must approximate the solution
$q_{T}(p_{k},X_{k})$ via some discrete process. There are many ways to
discretize Hamilton’s equations other than the simple Euler discretization;
see Neal (2011) for a discussion. In particular, using the leapfrog or
Störmer-Verlet method for integrating Hamilton’s equations leads to the
Hamiltonian Monte Carlo (HMC) algorithm. It simulates the Hamiltonian dynamics
for $K$ steps via the leapfrog integrator. At each iteration, given a integer
$K\geq 1$, a previous state $q_{0}$ and fresh
$p_{0}\sim\mathcal{N}(0,\mathbb{I}_{d})$, it runs the following updates for
$k=0,1,\ldots,K-1$:

$\displaystyle p_{k+\frac{1}{2}}$ $\displaystyle=p_{k}-\frac{\eta}{2}{\nabla
f}(q_{k})$ (8a) $\displaystyle q_{k+1}$ $\displaystyle=q_{k}+\eta
p_{k+\frac{1}{2}}$ (8b) $\displaystyle p_{k+1}$
$\displaystyle=p_{k+\frac{1}{2}}-\frac{\eta}{2}{\nabla f}(q_{k+1}).$ (8c)

Since discretizing the dynamics generates discretization error at each
iteration, it is followed by a Metropolis-Hastings adjustment where the
proposal $(p_{K},q_{K})$ is accepted with probability

$\displaystyle\min\left\\{1,\frac{\exp\left(-\mathcal{H}(p_{K},q_{K})\right)}{\exp\left(-\mathcal{H}(p_{0},q_{0})\right)}\right\\}.$
(9)

See Algorithm 3 for a detailed description of the HMC algorithm with leapfrog
integrator.

Input: Step size $\eta$, number of internal leapfrog updates $K$,

and a sample $x_{0}$ from a starting distribution $\mu_{0}$

Output: Sequence $x_{1},x_{2},\ldots$

1 for  _$i=0,1,\ldots$_ do

2 Proposal step:

3 $q_{0}\leftarrow x_{i}$

4 _Draw_ $p_{0}\sim\mathcal{N}(0,\mathbb{I}_{d})$

5 for  _$k=1,\ldots,K$_ do

6 $(p_{k},q_{k})\leftarrow\text{Leapfrog}(p_{k-1},q_{k-1},\eta)$

7

8 end for

9 % $q_{K}$ is now the new proposed state

10 Accept-reject step:

11 compute
$\alpha_{i+1}\leftarrow\displaystyle\min\left\\{1,\frac{\exp\left(-\mathcal{H}(p_{K},q_{K})\right)}{\exp\left(-\mathcal{H}(p_{0},q_{0})\right)}\right\\}$

12 With probability $\alpha_{i+1}$ _accept_ the proposal: $x_{i+1}\leftarrow
q_{K}$

13 With probability $1-\alpha_{i+1}$ _reject_ the proposal: $x_{i+1}\leftarrow
x_{i}$

14 end for

15 Program  _Leapfrog(_$p$ , $q$, $\eta$_)_:

16 $\tilde{p}\leftarrow p-\frac{\eta}{2}{\nabla f}(q)$

17 $\tilde{q}\leftarrow q+\eta\tilde{p}$

18 $\tilde{p}\leftarrow\tilde{p}-\frac{\eta}{2}{\nabla f}(\tilde{q})$

19

20 return _$(\tilde{p},\tilde{q})$_

Algorithm 3 Metropolized HMC with leapfrog integrator

##### Remark:

The HMC with leapfrog integrator can also be seen as a multi-step version of a
simpler Langevin algorithm. Indeed, running the HMC algorithm with $K=1$ is
equivalent to the MALA algorithm after a re-parametrization of the step-size
$\eta$. In practice, one also uses the HMC algorithm with a modified
Hamiltonian, in which the quadratic term $\left\|p\right\|_{2}^{2}$ is
replaced by a more general quadratic form $p^{T}Mp$. Here $M$ is a symmetric
positive definite matrix to be chosen by the user; see Appendix D.1.1 for
further discussion of this choice. In the main text, we restrict our analysis
to the case $M=I$.

##  3 Main results

We now turn to the statement of our main results. We remind the readers that
HMC refers to Metropolized HMC with leapfrog integrator, unless otherwise
specified. We first collect the set of assumptions for the target
distributions in Section 3.1. Following that in Section 3.2, we state our
results for HMC: first, we derive the mixing time bounds for general target
distributions in Theorem 1 and then apply that result to obtain concrete
guarantees for HMC with strongly log-concave target distributions. We defer
the discussion of weakly log-concave target distributions and (non-log-
concave) perturbations of log-concave distributions to Appendix C. In Section
3.3, we discuss the underlying results that are used to derive sharper mixing
time bounds using conductance profile (see Lemmas 3 and 4). Besides being
central to the proof of Theorem 1, these lemmas also enable a sharpening of
the mixing time guarantees for MALA and MRW (without much work), which we
state in Section 3.4.

###  3.1 Assumptions on the target distribution

In this section, we introduce some regularity notions and state the
assumptions on the target distribution that our results in the next section
rely on.

##### Regularity conditions:

A function $f$ is called: $\displaystyle L\text{-smooth}:\quad\ \
f(y)-f(x)-{\nabla f}(x)^{\top}(y-x)$
$\displaystyle\leq\frac{L}{2}\left\|x-y\right\|_{2}^{2},$ (10a) $\displaystyle
m\text{-strongly convex}:\quad\ \ f(y)-f(x)-{\nabla f}(x)^{\top}(y-x)$
$\displaystyle\geq\frac{m}{2}\left\|x-y\right\|_{2}^{2},$ (10b) $\displaystyle
L_{\text{H}}\text{-Hessian
Lipschitz}:\quad\quad\quad\quad\quad\left|\\!\left|\\!\left|{\nabla^{2}f}(x)-{\nabla^{2}f}(y)\right|\\!\right|\\!\right|_{{\text{op}}}$
$\displaystyle\leq L_{\text{H}}\left\|x-y\right\|_{2},$ (10c) where in all
cases, the inequalities hold for all $x,y\in{}^{d}$.

A distribution $\Pi$ with support $\mathcal{X}\subset{}^{d}$ is said to
satisfy the _isoperimetric inequality_ ($\mathfrak{a}=0$) or the _log-
isoperimetric inequality_ ($\mathfrak{a}=\frac{1}{2}$) with constant
$\psi_{\mathfrak{a}}$ if given any partition $S_{1},S_{2},S_{3}$ of
$\mathcal{X}$, we have

$\displaystyle\Pi(S_{3})\geq\frac{1}{2\psi_{\mathfrak{a}}}\cdot
d(S_{1},S_{2})\cdot\min\left\\{\Pi(S_{1}),\Pi(S_{2})\right\\}\cdot\log^{\mathfrak{a}}\left(1+\frac{1}{\min\left\\{\Pi(S_{1}),\Pi(S_{2}\right\\}}\right).$
(10d)

where the distance between two sets $S_{1},S_{2}$ is defined as
$d(S_{1},S_{2})=\inf_{x\in S_{1},y\in
S_{2}}\left\\{\left\|x-y\right\|_{2}\right\\}$. For a distribution $\Pi$ with
density $\pi$ and a given set $\Omega$, its restriction to $\Omega$ is the
distribution $\Pi_{\Omega}$ with the density
$\pi_{\Omega}(x)=\frac{\pi(x)\mathbf{1}_{\Omega}(x)}{\Pi(\Omega)}$.

##### Assumptions on the target distribution:

We introduce two sets of assumptions for the target distribution:

  1. (A)

We say that the target distribution $\Pi^{*}$ is
$(L,L_{\text{H}},s,\psi_{\mathfrak{a}},M)$-regular if the negative log density
$f$ is $L$-smooth (10a) and has $L_{\text{H}}$-Lipschitz Hessian (10c), and
there exists a convex measurable set $\Omega$ such that the distribution
$\Pi^{*}_{\Omega}$ is $\psi_{\mathfrak{a}}$-isoperimetric (10d), and the
following conditions hold:

$\displaystyle\Pi^{*}(\Omega)\geq 1-s\quad\mbox{and}\quad\left\|{\nabla
f}(x)\right\|_{2}$ $\displaystyle\leq M,\quad\text{for all }x\in\Omega.$ (10e)

  2. (B)

We say that the target distribution $\Pi^{*}$ is
$(L,L_{\text{H}},m)$-_strongly log-concave_ if the negative log density is
$L$-smooth (10a), $m$-strongly convex (10b), and $L_{\text{H}}$-Hessian-
Lipschitz (10c). Moreover, we use $x^{\star}$ to denote the unique mode of
$\Pi^{*}$ whenever $f$ is strongly convex.

Assumption (B) has appeared in several past papers on Langevin algorithms
(Dalalyan, 2016; Dwivedi et al., 2019; Cheng and Bartlett, 2017) and the
Lipschitz-Hessian condition (10c) has been used in analyzing Langevin
algorithms with inaccurate gradients (Dalalyan and Karagulyan, 2019) as well
as the unadjusted HMC algorithm (Mangoubi and Vishnoi, 2018). It is worth
noting Assumption (A) is strictly weaker than Assumption (B), since it allows
for distributions that are not log-concave. In Appendix B (see Lemma 15), we
show how Assumption (B) implies a version of Assumption (A).

###  3.2 Mixing time bounds for HMC

We start with the mixing time bound for HMC applied to any distribution
$\Pi^{*}$ satisfying Assumption (A). Let $\text{HMC-}(K,\eta)$ denote the
$\frac{1}{2}$-lazy Metropolized HMC algorithm with $\eta$ step size and $K$
leapfrog steps in each iteration. Let $\tau_{2}^{\text{\tiny
HMC}}(\epsilon;\mu_{0})$ denote the $\mathcal{L}_{2}$-mixing time (3b) for
this chain with the starting distribution $\mu_{0}$. We use $c$ to denote a
universal constant.

###### Theorem 1

Consider an $(L,L_{\text{H}},s,\psi_{\mathfrak{a}},M)$-regular target
distribution (cf. Assumption (A)) and a $\beta$-warm initial distribution
$\mu_{0}$. Then for any fixed target accuracy $\epsilon\in(0,1)$ such that
$\epsilon^{2}\geq 3\beta s$, there exist choices of the parameters $(K,\eta)$
such that $\text{HMC-}(K,\eta)$chain with $\mu_{0}$ start satisfies

$\displaystyle\tau_{2}^{\text{\tiny
HMC}}(\epsilon;\mu_{0})\leq\begin{cases}c\cdot\max\left\\{\log\beta,\displaystyle\frac{\psi_{\mathfrak{a}}^{2}}{K^{2}\eta^{2}}\log\left(\frac{\log\beta}{\epsilon}\right)\right\\}&\text{if
$\mathfrak{a}=\frac{1}{2}\qquad$ [log-
isoperimetric~{}\eqref{eq:assumption_isoperimetric}]}\\\\[5.69054pt]
c\cdot\displaystyle\frac{\psi_{\mathfrak{a}}^{2}}{K^{2}\eta^{2}}\log\left(\frac{\beta}{\epsilon}\right)&\text{if
$\mathfrak{a}=0\qquad$
[isoperimetric~{}\eqref{eq:assumption_isoperimetric}]}.\end{cases}$

See Section 5.2 for the proof, where we also provide explicit conditions on
$\eta$ and $K$ in terms of the other parameters (cf. equation (26b)).

Theorem 1 covers mixing time bounds for distributions that satisfy
isoperimetric or log-isoperimetric inequality provided that: (a) both the
gradient and Hessian of the negative log-density are Lipschitz; and (b) there
is a convex set that contains a large mass $(1-s)$ of the distribution. The
mixing time only depends on two quantities: the log-isoperimetric (or
isoperimetric) constant of the target distribution and the effective step-size
$K^{2}\eta^{2}$. As shown in the sequel, these conditions hold for log-concave
distributions as well as certain perturbations of them. If the distribution
satisfies a log-isoperimetric inequality, then the mixing time dependency on
the initialization warmness parameter $\beta$ is relatively weak
$O(\log\log\beta)$. On the other hand, when only an isoperimetric inequality
(but not log-isoperimetric) is available, the dependency is relatively larger
$O(\log\beta)$. In our current analysis, we can establish the
$\epsilon$-mixing time bounds up-to an error $\epsilon$ such that
$\epsilon^{2}\geq 3\beta s$. If mixing time bounds up to an arbitrary accuracy
are desired, then the distribution needs to satisfy (10e) for arbitrary small
$s$. For example, as we later show in Lemma 15, arbitary small $s$ can be
imposed for strongly log-concave densities (i.e., satisfying Assumption (B)).

Let us now derive several corollaries of Theorem 1. We begin with non-
asymptotic mixing time bounds for $\text{HMC-}(K,\eta)$ chain for strongly-log
concave target distributions. Then we briefly discuss the corollaries for
weakly log-concave target and non-log-concave target distributions and defer
the precise statements to Appendix C. These results also provide a basis for
comparison of our results with prior work.

####  3.2.1 Strongly log-concave target

We now state an explicit mixing time bound of HMC for a strongly log-concave
distribution. We consider an $(L,L_{\text{H}},m)$-strongly log-concave
distribution (assumption (B)). We use $\kappa=L/m$ to denote the condition
number of the distribution. Our result makes use of the following function

$\displaystyle r(s)$
$\displaystyle:\,=1+\max\left\\{\left(\frac{\log(1/s)}{d}\right)^{1/4},\;\left(\frac{\log(1/s)}{d}\right)^{1/2}\right\\},$
(11a) for $s>0$, and involves the step-size choices
$\displaystyle\eta_{\mbox{\tiny{warm}}}=\sqrt{\frac{1}{cL\cdot
r\left(\frac{\epsilon^{2}}{3\beta}\right)d^{\frac{7}{6}}}}\;,\quad\mbox{and}\quad\eta_{\mbox{\tiny{feas}}}=\sqrt{\frac{1}{cL\cdot
r\left(\frac{\epsilon^{2}}{2\kappa^{d}}\right)}\min\left\\{\frac{1}{d\kappa^{\frac{1}{2}}},\frac{1}{d^{\frac{2}{3}}\kappa^{\frac{5}{6}}},\frac{1}{d^{\frac{1}{2}}\kappa^{\frac{3}{2}}}\right\\}}.$
(11b)

With these definitions, we have the following:

###### Corollary 2

Consider an $(L,L_{\text{H}},m)$-strongly log-concave target distribution
$\Pi^{*}$ (cf. Assumption (B)) such that $L_{\text{H}}^{2/3}=O(L)$, and any
error tolerance $\epsilon\in(0,1)$.

  1. (a)

Suppose that $\kappa=O(d^{\frac{2}{3}})$ and
$\beta=O\left(\exp\left(d^{\frac{2}{3}}\right)\right)$. Then with any
$\beta$-warm initial distribution $\mu_{0}$, hyper-parameters
$K=d^{\frac{1}{4}}$ and $\eta=\eta_{\mbox{\tiny{warm}}}$, the
$\text{HMC-}(K,\eta)$ chain satisfies

$\displaystyle\tau_{2}^{\text{\tiny HMC}}(\epsilon;\mu_{0})\leq
c\;d^{\frac{2}{3}}\;\kappa\,r\left(\frac{\epsilon^{2}}{3\beta}\right)\;\log\left(\frac{\log\beta}{\epsilon}\right).$
(12a)

  2. (b)

With the initial distribution
$\mu_{\dagger}=\mathcal{N}(x^{\star},\frac{1}{L}\mathbb{I}_{d})$, hyper-
parameters $K=\kappa^{\frac{3}{4}}$ and $\eta=\eta_{\mbox{\tiny{feas}}}$, the
$\text{HMC-}(K,\eta)$ chain satisfies

$\displaystyle\tau_{2}^{\text{\tiny HMC}}(\epsilon;\mu_{\dagger})\leq
c\;r\left(\frac{\epsilon^{2}}{3\kappa^{d}}\right)\,\max\left\\{d\log{\kappa},\max\left[d,d^{\frac{2}{3}}\kappa^{\frac{1}{3}},d^{\frac{1}{2}}\kappa\right]\log\left(\frac{d\log{\kappa}}{\epsilon}\right)\right\\}.$
(12b)

See Appendix B for the proof. In the same appendix, we also provide a more
refined mixing time of the HMC chain for a more general choice of hyper-
parameters (see Corollary 14). In fact, as shown in the proof, the assumption
$L_{\text{H}}^{2/3}=O(L)$ is not necessary in order to control mixing; rather,
we adopted it above to simplify the statement of our bounds. A more detailed
discussion on the particular choice for step size $\eta$ is provided in
Appendix D.

##### MALA vs HMC—Warm start:

Corollary 2 provides mixing time bounds for two cases. The first result (12a)
implies that given a warm start for a well-conditioned strongly log concave
distribution, i.e., with constant $\beta$ and $\kappa\ll d$, the
$\epsilon$-$\mathcal{L}_{2}$-mixing time666Note that $r(\epsilon^{2})\leq 6$
for $\epsilon\geq\frac{2}{e^{d/2}}$ and thus we can treat $r$ as a small
constant for a large range of $\epsilon$. Otherwise, if $\epsilon$ needs to be
extremely small, the results still hold with an extra
$\log^{\frac{1}{2}}\left(\frac{1}{\epsilon}\right)$ dependency. of HMC scales
$\widetilde{O}(d^{\frac{2}{3}}\log(1/\epsilon))$. It is interesting to compare
this guarantee with known bounds for the MALA algorithm. However since each
iteration of MALA uses only a single gradient evaluation, a fair comparison
would require us to track the total number of gradient evaluations required by
the $\text{HMC-}(K,\eta)$ chain to mix. For HMC to achieve accuracy
$\epsilon$, the total number of gradient evaluations is given by
$K\cdot\tau_{2}^{\text{\tiny HMC}}(\epsilon;\mu_{0})$, which in the above
setting, scales as $\widetilde{O}(d^{\frac{11}{12}}\kappa\log(1/\epsilon))$.
This rate was also summarized in Table 1. On the other hand, Theorem 1 in
Dwivedi et al. (2019) shows that the corresponding number of gradient
evaluations for MALA is $\widetilde{O}(d\kappa\log(1/\epsilon))$. As a result,
we conclude that the upper bound for HMC is $d^{\frac{1}{12}}$ better than the
known upper bound for MALA with a warm start for a well-conditioned strongly
log concave target distribution. We summarize these rates in Table 2. Note
that MRW is a zeroth order algorithm, which makes use of function evaluations
but not gradient information.

Sampling algorithm |  Mixing time |  #Gradient evaluations
---|---|---
MRW (Dwivedi et al., 2019, Theorem 2) | $d\kappa^{2}\cdot\log\frac{1}{\epsilon}$ | NA
MALA (Dwivedi et al., 2019, Theorem 1) | $d\kappa\cdot\log\frac{1}{\epsilon}$ | $d\kappa\cdot\log\frac{1}{\epsilon}$
$\text{HMC-}(K,\eta)$ [ours, Corollary 2] | $d^{\frac{2}{3}}\kappa\cdot\log\frac{1}{\epsilon}$ | $d^{\frac{11}{12}}\kappa\cdot\log\frac{1}{\epsilon}$
Table 2: Summary of the $\epsilon$-mixing time and the corresponding number of
gradient evaluations for MRW, MALA and HMC from a _warm start_ with an
$(L,L_{\text{H}},m)$-strongly-log-concave target. These statements hold under
the assumption $L_{\text{H}}^{2/3}=O(L)$, $\kappa=\frac{L}{m}\ll d$, and omit
logarithmic terms in dimension.

##### MALA vs HMC—Feasible start:

In the second result (12b), we cover the case when a warm start is not
available. In particular, we analyze the HMC chain with the feasible initial
distribution $\mu_{\dagger}=\mathcal{N}(x^{\star},\frac{1}{L}\mathbb{I}_{d})$.
Here $x^{\star}$ denotes the unique mode of the target distribution and can be
easily computed using an optimization scheme like gradient descent. It is not
hard to show (see Corollary 1 in Dwivedi et al. (2019)) that for an $L$-smooth
(10a) and $m$ strongly log-concave target distribution (10b), the distribution
$\mu_{\dagger}$ acts as a $\kappa^{d/2}$-warm start distribution. Once again,
it is of interest to determine whether HMC takes fewer gradient steps when
compared to MALA to obtain an $\epsilon$-accurate sample. We summarize the
results in Table 3, with log factors hidden, and note that HMC with
$K=\kappa^{3/4}$ is faster than MALA for as long as $\kappa$ is not too large.
From the last column, we find that when $\kappa\ll d^{\frac{1}{2}}$, HMC is
faster than MALA by a factor of $\kappa^{\frac{1}{4}}$ in terms of number of
gradient evaluations. It is worth noting that for the feasible start
$\mu_{\dagger}$, the mixing time bounds for MALA and MRW in Dwivedi et al.
Dwivedi et al. (2019) were loose by a factor $d$ when compared to the tighter
bounds in Theorem 5 derived later in this paper.

Sampling algorithm |  Mixing time |  # Gradient Evaluations
---|---|---
|  |  general $\kappa$ | $\kappa\ll d^{\frac{1}{2}}$
MRW [ours, Theorem 5] | $d\kappa^{2}$ | NA | NA
MALA [ours, Theorem 5] | $\max\left\\{d\kappa,d^{\frac{1}{2}}\kappa^{\frac{3}{2}}\right\\}$ | $\max\left\\{d\kappa,d^{\frac{1}{2}}\kappa^{\frac{3}{2}}\right\\}$ | $d\kappa$
$\text{HMC-}(K,\eta)$ [ours, Corollary 2] | $\max\big{\\{}d,d^{\frac{2}{3}}\kappa^{\frac{1}{3}},d^{\frac{1}{2}}\kappa\big{\\}}$ | $\max\big{\\{}d\kappa^{\frac{3}{4}},d^{\frac{2}{3}}\kappa^{\frac{13}{12}},d^{\frac{1}{2}}\kappa^{\frac{7}{4}}\big{\\}}$ | $d\kappa^{\frac{3}{4}}$
Table 3:  Summary of the $\epsilon$-mixing time and the corresponding number
of gradient evaluations for MRW, MALA and HMC from the _feasible start_
$\mu_{\dagger}=\mathcal{N}(x^{\star},\frac{1}{L}\mathbb{I}_{d})$ for an
$(L,L_{\text{H}},m)$-strongly-log-concave target. Here $x^{\star}$ denotes the
unique mode of the target distribution. These statements hold uner the
assumption $L_{\text{H}}=O(L^{\frac{3}{2}})$, and hide the logarithmic factors
in $\epsilon,d$ and $\kappa=L/m$.

##### Metropolized HMC vs Unadjusted HMC:

There are many recent results on the $1$-Wasserstein distance mixing of
unadjusted versions of HMC (for instance, see the papers Mangoubi and Vishnoi
(2018); Lee et al. (2018)). For completeness, we compare our results with them
in the Appendix D.2; in particular, see Table 7 for a comparative summary. We
remark that comparisons of these different results is tricky for two reasons:
(a) The $1$-Wasserstein distance and the total variation distance are not
strictly comparable, and, (b) the unadjusted HMC results always have a
polynomial dependence on the error parameter $\epsilon$ while our results for
Metropolized HMC have a superior logarithmic dependence on $\epsilon$.
Nonetheless, the second difference between these chains has a deeper
consequence, upon which we elaborate further in Appendix D.2. On one hand, the
unadjusted chains have better mixing time in terms of scaling with $d$, if we
fix $\epsilon$ or view it as independent of $d$. On the other hand, when such
chains are used to estimate certain higher-order moments, the polynomial
dependence on $\epsilon$ might become the bottleneck and Metropolis-adjusted
chains would become the method of choice.

##### Ill-conditioned target distributions:

In order to keep the statement of Corollary 2 simple, we stated the mixing
time bounds of $\text{HMC-}(K,\eta)$-chain only for a particular choice of
$(K,\eta)$. In our analysis, this choice ensures that HMC is better than MALA
only when condition number $\kappa$ is small. For Ill-conditioned
distributions, i.e., when $\kappa$ is large, finer tuning of
$\text{HMC-}(K,\eta)$-chain is required. In Appendices B and D (see Table 4
for the hyper-parameter choices), we show that HMC is strictly better than
MALA as long as $\kappa\leq d$ and as good as MALA when $\kappa\geq d$.

##### Beyond strongly-log-concave:

The proof of Corollary 2 is based on the fact that
$(L,L_{\text{H}},m)$-strongly-log-concave distribution is in fact an
$(L,L_{\text{H}},s,\psi_{1/2},M_{s})$ -regular distribution for any
$s\in(0,1)$. Here $\psi_{1/2}={1}/{\sqrt{m}}$ is fixed and the bound on the
gradient $M_{s}=r(s)\sqrt{{d}/{m}}$ depends on the choice of $s$. The result
is formally stated in Lemma 15 in Appendix B. Moreover, in Appendix C, we
discuss the case when the target distribution is weakly log concave (under a
bounded fourth moment or bounded covariance matrix assumption) or a
perturbation of log-concave distribution. See Corollary 18 for specific
details where we provide explicit expressions for the rates that appear in
third and fourth columns of Table 1.

###  3.3 Mixing time bounds via conductance profile

In this section, we discuss the general results that form the basis of the
analysis in this paper. A standard approach to controlling mixing times is via
worst-case conductance bounds. This method was introduced by Jerrum and
Sinclair (1988) for discrete space chains and then extended to the continuous
space settings by Lovász and Simonovits (1993), and has been thoroughly
studied. See the survey (Vempala, 2005) and the references therein for a
detailed discussion of conductance based methods for continuous space Markov
chains.

Somewhat more recent work on discrete state chains has introduced more refined
methods, including those based on the conductance profile (Lovász and Kannan,
1999), the spectral and conductance profile (Goel et al., 2006), as well as
the evolving set method (Morris and Peres, 2005). Here we extend one of the
conductance profile techniques from the paper by Goel et al. (2006) from
discrete state to continuous state chains, albeit with several appropriate
modifications suited for the general setting.

We first introduce some background on the conductance profile. Given a Markov
chain with transition probability
$\Theta:\mathcal{X}\times\mathcal{B}\left(\mathcal{X}\right)\rightarrow\real$,
its stationary _flow_ $\phi:\mathcal{B}(\mathcal{X})\rightarrow\real$ is
defined as

$\displaystyle\phi(S)$ $\displaystyle=\int_{x\in
S}\Theta(x,S^{c})\pi^{*}(x)dx\quad\text{for any
}S\in\mathcal{B}(\mathcal{X}).$ (13)

Given a set $\Omega\subset\mathcal{X}$, the _$\Omega$ -restricted conductance
profile_ is given by

$\displaystyle\Phi_{\Omega}(v)=\inf_{\Pi^{*}(S\cap\Omega)\in(0,v]}\frac{\phi(S)}{\Pi^{*}(S\cap\Omega)}\qquad\text{for
any }v\in\big{(}0,\;\Pi^{*}(\Omega)/2\big{]}.$ (14)

(The classical conductance constant $\Phi$ is a special case; it can be
expressed as $\Phi=\Phi_{\mathcal{X}}(\frac{1}{2})$.) Moreover, we define the
_truncated extension_ $\tilde{\Phi}_{\Omega}$ of the function $\Phi_{\Omega}$
to the positive real line as

$\displaystyle\tilde{\Phi}_{\Omega}(v)=\begin{cases}\Phi_{\Omega}(v),&v\in\left(0,\frac{\Pi^{*}(\Omega)}{2}\right]\\\
\Phi_{\Omega}(\Pi^{*}(\Omega)/2),&v\in\left[\frac{\Pi^{*}(\Omega)}{2},\infty\right).\end{cases}$
(15)

In our proofs we use the conductance profile with a suitably chosen set
$\Omega$.

##### Smooth chain assumption:

We say that the Markov chain satisfies the _smooth chain assumption_ if its
transition probability function
$\Theta:\mathcal{X}\times\mathcal{B}(\mathcal{X})\rightarrow{}_{+}$ can be
expressed in the form

$\displaystyle\Theta(x,dy)=\theta(x,y)dy+\alpha_{x}\mathbf{\delta}_{x}(dy)\quad\text{for
all }x,y\in\mathcal{X},$ (16)

where $\theta$ is the transition kernel satisfying $\theta(x,y)\geq 0$ for all
$x,y\in\mathcal{X}$. Here $\mathbf{\delta}_{x}$ denotes the Dirac-delta
function at $x$ and consequently, $\alpha_{x}$ denotes the one-step
probability of the chain to stay at its current state $x$. Note that the three
algorithms discussed in this paper (MRW, MALA and HMC) all satisfy the smooth
chain assumption (16). Throughout the paper, when dealing with a general
Markov chain, we assume that it satisfies the smooth chain assumption.

##### Mixing time via conductance profile:

We now state our Lemma 3 that provides a control on the mixing time of a
Markov chain with continuous-state space in terms of its restricted
conductance profile. We show that this control (based on conductance profile)
allows us to have a better initialization dependency than the usual
conductance based control (see Lovász and Simonovits (1990, 1993); Dwivedi et
al. (2019)). This method for sharpening the dependence is known for discrete-
state Markov chains; to the best of our knowledge, the following lemma is the
first statement and proof of an analogous sharpening for continuous state
space chains:

###### Lemma 3

Consider a reversible, irreducible, $\zeta$-lazy and smooth Markov chain (16)
with stationary distribution $\Pi^{*}$. Then for any error tolerance
$\epsilon$, and a $\beta$-warm distribution $\mu_{0}$, given a set $\Omega$
such that $\Pi^{*}(\Omega)\geq 1-\frac{\epsilon^{2}}{3\beta^{2}}$, the
$\epsilon$-$\mathcal{L}_{2}$ mixing time of the chain is bounded as

$\displaystyle\tau_{2}(\epsilon;\mu_{0})\leq\int_{4/\beta}^{8/\epsilon^{2}}\frac{8\>dv}{\zeta\cdot
v{\tilde{\Phi}_{\Omega}}^{2}(v)},$ (17)

where $\tilde{\Phi}_{\Omega}$ denotes the truncated $\Omega$-restricted
conductance profile (15).

See Appendix A.1 for the proof, which is based on an appropriate
generalization of the ideas used by Goel et al. (2006) for discrete state
chains.

The standard conductance based analysis makes use of the worst-case
conductance bound for the chain. In contrast, Lemma 3 relates the mixing time
to the conductance profile, which can be seen as point-wise conductance. We
use the $\Omega$-restricted conductance profile to state our bounds, because
often a Markov chain has poor conductance only in regions that have very small
probability under the target distribution. Such a behavior is not disastrous
as it does not really affect the mixing of the chain up to a suitable
tolerance. Given the bound (17), we can derive mixing time bound for a Markov
chain readily if we have a bound on the $\Omega$-restricted conductance
profile $\Phi_{\Omega}$ for a suitable $\Omega$. More precisely, if the
$\Omega$-restricted conductance profile $\Phi_{\Omega}$ of the Markov chain is
bounded as

$\displaystyle\Phi_{\Omega}(v)\geq\sqrt{B\log\left(\frac{1}{v}\right)}\quad\text{for
}v\in\left[\frac{4}{\beta},\frac{1}{2}\right],$

for some $\beta>0$ and $\Omega$ such that $\Pi^{*}(\Omega)\geq
1-\frac{\epsilon^{2}}{3\beta^{2}}$. Then with a $\beta$-warm start, Lemma 3
implies the following useful bound on the mixing time of the $\zeta$-lazy
Markov chain:

$\displaystyle\tau_{2}(\epsilon;\mu_{0})\leq\frac{64}{\zeta
B}\log\left(\frac{\log\beta}{2\epsilon}\right).$ (18)

We now relate our result to prior work based on conductance profile.

##### Prior work:

For discrete state chains, a result similar to Lemma 3 was already proposed by
Lovász and Kannan (Theorem 2.3 in Lovász and Kannan (1999)). Later on, Morris
and Peres (2005) and Goel et al. (2006) used the notion of evolving sets and
spectral profile respectively to sharpen the mixing time bounds based on
average conductance for discrete-state space chains. In the context of
continuous state space chains, Lovász and Kannan claimed in their original
paper (Lovász and Kannan, 1999) that a similar result should hold for general
state space chain as well, although we were unable to find any proof of such a
general result in that or any subsequent work. Nonetheless, in a later work an
average conductance based bound was used by Kannan et al. to derive faster
mixing time guarantees for uniform sampling on bounded convex sets for ball
walk (see Section 4.3 in Kannan et al. (2006)). Their proof technique is not
easily extendable to more general distributions including the general log-
concave distributions in d. Instead, our proof of Lemma 3 for general state
space chains proceeds by an appropriate generalization of the ideas based on
the spectral profile by Goel et al. (2006) (for discrete state chains).

##### Lower bound on conductance profile:

Given the bound (18), it suffices to derive a lower bound on the conductance
profile $\Phi_{\Omega}$ of the Markov chain with a suitable choice of the set
$\Omega$. We now state a lower bound for the restricted-conductance profile of
a general state space Markov chain that comes in handy for this task. We note
that a closely related logarithmic-Cheeger inequality was used for sampling
from uniform distribution of a convex body (Kannan et al., 2006) and for
sampling from log-concave distributions (Lee and Vempala, 2018b) without
explicit constants. Since we would like to derive a non-asymptotic mixing
rate, we re-derive an explicit form of their result.

Let scalars $s\in(0,1/2]$, $\omega\in(0,1)$ and $\Delta>0$ be given and let
$\mathcal{T}_{x}$ denote the one-step transition distribution of the Markov
chain at point $x$. Suppose that that chain satisfies

$\displaystyle{d_{\text{\tiny{TV}}}\left(\mathcal{T}_{x},\mathcal{T}_{y}\right)\leq
1-\omega}\qquad\mbox{ whenever $x,y\in\Omega$ and
$\left\|x-y\right\|_{2}\leq\Delta$.}$ (19)

###### Lemma 4

For a given target distribution $\Pi^{*}$, let $\Omega$ be a convex measurable
set such that the distribution $\Pi^{*}_{\Omega}$ satisfies the isoperimetry
(or log-isoperimetry) condition (10d) with $\mathfrak{a}=0$ (or
$\mathfrak{a}=\frac{1}{2}$ respectively). Then for any Markov chain satisfying
the condition (19), we have

$\displaystyle\Phi_{\Omega}(v)\geq\frac{\omega}{4}\cdot\min\left\\{1,\frac{\Delta}{16\psi_{\mathfrak{a}}}\cdot\log^{\mathfrak{a}}\left(1+\frac{1}{v}\right)\right\\},\quad\text{for
any }v\in\left[0,\frac{\Pi^{*}(\Omega)}{2}\right].$ (20)

See Appendix A.2 for the proof; the extra logarithmic term comes from the
logarithmic isoperimetric inequality ($\mathfrak{a}=\frac{1}{2}$).

##### Faster mixing time bounds:

For any target distribution satisfying a logarithmic isoperimetric inequality
(including the case of a strongly log-concave distribution), Lemma 4 is a
strict improvement of the conductance bounds derived in previous works
(Lovász, 1999; Dwivedi et al., 2019). Given this result, suppose that we can
find a convex set $\Omega$ such that $\Pi^{*}(\Omega)\approx 1$ and the
conditions of Lemma 4 are met, then with a $\beta$-warm start $\mu_{0}$, a
direct application of the bound (18) along with Lemma 4 implies the following
bound:

$\displaystyle\tau_{2}(\epsilon;\mu_{0})\leq
O\left(\frac{1}{\omega^{2}\Delta^{2}}\log\frac{\log\beta}{\epsilon}\right).$
(21)

Results known from previous work for continuous state Markov chains scale like
$\frac{\log(\beta/\epsilon)}{\omega^{2}\Delta^{2}}$; for instance, see Lemma 6
in Chen et al. (2018). In contrast, the bound (21) provides an additional
logarithmic factor improvement in the factor $\beta$. Such an improvement also
allows us to derive a sharper dependency on dimension $d$ for the mixing time
for sampling algorithms other than HMC as we now illustrate in the next
section.

###  3.4 Improved warmness dependency for MALA and MRW

As discussed earlier, the bound (21) helps derive a
$\frac{\log\log\beta}{\log\beta}$ factor improvement in the mixing time bound
from a $\beta$-warm start in comparison to earlier conductance based results.
In many settings, a suitable choice of initial distribution has a warmness
parameter that scales exponentially with dimension $d$, e.g.,
$\beta=O(e^{d})$. For such cases, this improvement implies a gain of
$O(\frac{d}{\log d})$ in mixing time bounds. As already noted the distribution
$\mu_{\dagger}=\mathcal{N}(x^{*},\frac{1}{L}\mathbb{I}_{d})$ is a feasible
starting distribution, whose warmness scales exponentially with dimension $d$.
See Section 3.2 of the paper (Dwivedi et al., 2019), where the authors show
that computing $x^{*}$ is not expensive and even approximate estimates of
$x^{*}$ and $L$ are sufficient to provide a feasible starting distribution. We
now state sharper mixing time bounds for MALA and MRW with $\mu_{\dagger}$ as
the starting distribution. In the result, we use $c_{1},c_{2}$ to denote
positive universal constants.

###### Theorem 5

Assume that the target distribution $\Pi^{*}$ satisfies the conditions (10a)
and (10b) (i.e., the negative log-density is $L$-smooth and $m$-strongly
convex). Then given the initial distribution
$\mu_{\dagger}=\mathcal{N}(x^{*},\frac{1}{L}\mathbb{I}_{d})$, the
$\frac{1}{2}$-lazy versions of MRW and MALA (Algorithms 1 and 2) with step
sizes

$\displaystyle\eta_{\text{\tiny
MRW}}=c_{2}\cdot\frac{1}{Ld\kappa},\quad\text{and}\quad\eta_{\text{\tiny
MALA}}=c_{1}\cdot\frac{1}{Ld\cdot\max\left\\{1,\sqrt{\kappa/d}\right\\}}$ (22)

respectively, satisfy the mixing time bounds

$\displaystyle\tau_{2}^{\text{\tiny MRW}}(\epsilon;\mu_{0})$
$\displaystyle=O\left(d\kappa^{2}\log{\frac{d}{\epsilon}}\right),\quad\text{and}\quad$
(23a) $\displaystyle\tau_{2}^{\text{\tiny MALA}}(\epsilon;\mu_{0})$
$\displaystyle=O\left(d\kappa\log{\frac{d}{\epsilon}}\cdot\max\left\\{1,\sqrt{\frac{\kappa}{d}}\right\\}\right).$
(23b)

The proof is omitted as it directly follows from the conductance profile based
mixing time bound in Lemma 3, Lemma 4 and the overlap bounds for MALA and MRW
provided in our prior work (Dwivedi et al., 2019). Theorem 5 states that the
mixing time bounds for MALA and MRW with the feasible distribution
$\mu_{\dagger}$ as the initial distribution scale as
$\widetilde{O}(d\kappa\log\left(1/\epsilon\right))$ and
$\widetilde{O}(d\kappa^{2}\log\left(1/\epsilon\right))$. Once again, we note
that in light of the inequality (3c) we obtain same bounds for the number of
steps taken by these algorithms to mix within $\epsilon$ total-variation
distance of the target distribution $\Pi^{*}$. Consequently, our results
improve upon the previously known mixing time bounds for MALA and MRW (Dwivedi
et al., 2019) for strongly log-concave distributions. With $\mu_{\dagger}$ as
the initial distribution, the authors had derived bounds of order
$\widetilde{O}(d^{2}\kappa\log\left(1/\epsilon\right))$ and
$\widetilde{O}(d^{2}\kappa^{2}\log\left(1/\epsilon\right))$ for MALA and MRW
respectively (cf. Corollary 1 in Dwivedi et al. (2019)). However, the
numerical experiments in that work suggested a better dependency on the
dimension for the mixing time. Indeed the mixing time bounds from Theorem 5
are smaller by a factor of $\frac{d}{\log d}$, compared to our earlier bounds
in the prior work (Dwivedi et al., 2019) for both of these chains thereby
resolving an open question. Nonetheless, it is still an open question how to
establish a lower bound on the mixing time of these sampling algorithms.

##  4 Numerical experiments

In this section, we numerically compare HMC with MALA and MRW to verify that
our suggested step-size and leapfrog steps lead to faster convergence for the
HMC algorithm. We adopt the step-size choices for MALA and MRW given in
Dwivedi et al. (2019), whereas the choices for step-size and leapfrog steps
for HMC are taken from Corollary 2 in this paper. When the Hessian-Lipschitz
constant $L_{\text{H}}$ is small, our theoretical results suggest that HMC can
be run with much larger step-size and much larger number of leapfrog steps
(see Appendix D.1.1). Since our experiments make use of multivariate Gaussian
target distribution, the Hessian-Lipschitz constant $L_{\text{H}}$ is always
zero. Consequently we also perform experiments with a more _aggressive_ choice
of parameters, i.e., larger step-size and number of leapfrog steps. We denote
this choice by HMCagg.

In this simulation, we check the dimension $d$ dependency and condition number
$\kappa$ dependency in the multivariate Gaussian case under our step-size
choices. We consider sampling from the multivariate Gaussian distribution with
density

$\displaystyle\Pi^{*}(x)\propto e^{-\frac{1}{2}x^{\top}\Sigma^{-1}x},$ (24)

for some covariance matrix $\Sigma\in{}^{d\times d}$. The log density
(disregarding constants) and its deriviatives are given by

$\displaystyle f(x)=\frac{1}{2}x^{\top}\Sigma^{-1}x,\quad{\nabla
f}(x)=\Sigma^{-1}x,\quad\text{and}\quad{\nabla^{2}f}(x)=\Sigma^{-1}.$

Consequently, the function $f$ is strongly convex with parameter
$m=1/\lambda_{\text{max}}(\Sigma)$ and smooth with parameter
$L=1/\lambda_{\text{min}}(\Sigma)$. Since
$\mathcal{L}_{\mathfrak{p}}$-divergence can not be measure with finitely many
samples, we use the error in quantiles along different directions for
convergence diagnostics. Using the exact quantile information for each
direction for Gaussians, we measure the error in the $75\%$ quantile of the
relative sample distribution and the true distribution in the _least favorable
direction_ , i.e., along the eigenvector of $\Sigma$ corresponding to the
eigenvalue $\lambda_{\text{max}}(\Sigma)$. The quantile mixing time is defined
as the smallest iteration when this relative error falls below a constant
$\delta=0.04$. We use $\mu_{0}=\mathcal{N}\left(0,L^{-1}\mathbb{I}_{d}\right)$
as the initial distribution. To make the comparison with MRW and MALA fair, we
compare the number of total function and gradient evaluations instead of
number of iterations. For HMC, the number of gradient evaluations is $K$ times
the number of outer-loop iterations.

For every case of simulation, the parameters for $\text{HMC-}(K,\eta)$ are
chosen according to the warm start case in Corollary 2 with $K=4\cdot
d^{1/4}$, and for MRW and MALA are chosen according to the paper Dwivedi et
al. (2019). As alluded to earlier, we also run the HMC chain a more aggressive
choice of parameters, and denote this chain by HMCagg. For HMCagg, both the
step-size and leapfrog steps are larger (Appendix D.1.1): $K=4\cdot
d^{1/8}\kappa^{1/4}$ where we take into account that $L_{\text{H}}$ is zero
for Gaussian distribution. We simulate $100$ independent runs of the four
chains, MRW, MALA, HMC, HMCagg, and for each chain at every iteration we
compute the quantile error across the $100$ samples from $100$ independent
runs of that chain. We compute the minimum number of total function and
gradient evaluations required for the relative quantile error to fall below
$\delta=0.04$. We repeat this computation $10$ times and report the averaged
number of total function and gradient evaluations in Figure 2. To examine the
scaling of the number of evaluations with the dimension $d$, we vary
$d\in\left\\{2,4,\ldots,128\right\\}$. For each chain, we also fit a least
squares line for the number of total function and gradient evaluations with
respect to dimension $d$ on the log-log scale, and report the slope in the
figure. Note that a slope of $\alpha$ would denote that the number of
evaluations scales with $d$ as $d^{\alpha}$.

##### (a) Dimension dependency for fixed $\kappa$:

First, we consider the case of fixed condition number. We fix $\kappa=4$ while
we vary the dimensionality $d$ of the target distribution is varied over
$\left\\{2,4,\ldots,128\right\\}$. The Hessian $\Sigma$ in the multivariate
Gaussian distribution is chosen to be diagonal and the square roots of its
eigenvalues are linearly spaced between $1.0$ to $2.0$. Figure 2(a) shows the
dependency of the number of total function and gradient evaluations as a
function of dimension $d$ for the four Markov chains on the log-log scale. The
least-squares fits of the slopes for HMC, HMCagg, MALA and MRW are $0.80(\pm
0.12)$, $0.58(\pm 0.15)$, $0.93(\pm 0.13)$ and $0.96(\pm 0.10)$, respectively,
where standard errors of the regression coefficient is reported in the
parentheses. These numbers indicate close correspondence to the theoretical
slopes (reported in Table 2 and Appendix D.1.1) are $0.92,0.63,1.0,1.0$
respectively.

|
---|---
(a) $\kappa=4$ | (b) $\kappa=d^{\frac{2}{3}}$
Figure 2: Average number of total function and gradient evaluations as a
function of dimension for four random walks on multivariate Gaussian density
(24) where the covariance has a condition number $\kappa$ that is (a) constant
$4$ and (b) scales with dimension $d$. With suggested step-size and leapfrog
steps in Corollary 2, the number of total function and gradient evaluations of
HMC has a smaller dimension dependency than that of MALA or MRW. Since the
target distributon is Gaussian and the Hessian-Lipschitz constant
$L_{\text{H}}$ is zero, larger step-size and larger number of leapfrog steps
can be chosen according to Appendix D.1.1. The plots does show that HMCagg
with larger step-size and larger number of leapfrog steps uses smaller number
of total function and gradient evaluations to achieve the same quantile
mixing.

##### (b) Dimension dependency for $\kappa=d^{2/3}$:

Next, we consider target distributions such that their condition number varies
with $d$ as $\kappa=d^{2/3}$, where $d$ is varied from $2$ to $128$. To ensure
such a scaling for $\kappa$, we choose the Hessian $\Sigma$ for the
multivariate Gaussian distribution to be diagonal and set the square roots of
its eigenvalues linearly spaced between $1.0$ to $d^{1/3}$. Figure 2(b) shows
the dependency of the number of total function and gradient evaluations as a
function of dimension $d$ for the four random walks on the log-log scale. The
least squares fits yield the slopes as $1.60(\pm 0.09)$, $1.34(\pm 0.17)$,
$1.64(\pm 0.11)$ and $2.25(\pm 0.08)$ for HMC, HMCagg, MALA and MRW,
respectively, where standard errors of the regression coefficient are reported
in the parentheses. Recall that the theoretical guarantees for HMC (Table 5),
HMCagg (Table 6), MALA and MRW (Table 2) yield that these exponent should be
close to 1.58, 1.46, 1.67 and 2.33 respectively. Once again, we observe a good
agreement of the numerical results with that of our theoretical results.

##### Remark:

We would like to caution that the aggressive parameter choices for HMCagg are
well-informed only when the Hessian-Lipschitz constant $L_{\text{H}}$ is
small—which indeed is the case for the Gaussian target distributions
considered above. When general log-concave distributions are considered, one
may use the more general choices recommended in Corollary 14. See Appendix D
for an in-depth discussion on different scenarios and the optimal parameter
choices derived from our theory.

##  5 Proofs

This section is devoted primarily to the proof of Theorem 1. In order to do
so, we begin with the mixing time bound based on the conductance profile from
Lemma 3. We then seek to apply Lemma 4 in order derive a bound on the
conductance profile itself. However, in order to do so, we need to derive
bound on the overlap between the proposal distributions of HMC at two nearby
points and show that the Metropolis-Hastings step only modifies the proposal
distribution by a relatively small amount. This control is provided by Lemma
6, stated in Section 5.1. We use it to prove Theorem 1 in Section 5.2.
Finally, Section 5.3 is devoted to the proof of Lemma 6. We provide a sketch-
diagram for how various main results of the paper interact with each other in
Figure 1.

###  5.1 Overlap bounds for HMC

In this subsection, we derive two important bounds for the Metropolized HMC
chain: (1) first, we quantify the overlap between proposal distributions of
the chain for nearby points, and, (2) second, we show that the distortion in
the proposal distribution introduced by the Metropolis-Hastings accept-reject
step can be controlled if an appropriate step-size is chosen. Putting the two
pieces together enables us to invoke Lemma 4 to prove Theorem 1.

In order to do so, we begin with some notation. Let $\mathcal{T}$ denote the
transition operator of the HMC chain with leapfrog integrator taking step-size
$\eta$ and number of leapfrog updates $K$. Let $\mathcal{P}_{x}$ denote the
proposal distribution at $x\in\mathcal{X}$ for the chain before the accept-
reject step and the lazy step. Let $\mathcal{T}^{\text{before-lazy}}_{x}$
denote the corresponding transition distribution after the proposal and the
accept-reject step, before the lazy step. By definition, we have

$\displaystyle\mathcal{T}_{x}(A)=\zeta\mathbf{\delta}_{x}(A)+(1-\zeta)\mathcal{T}^{\text{before-
lazy}}_{x}(A)\qquad\mbox{for any measurable set
$A\in\mathcal{B}(\mathcal{X})$.}$ (25)

Our proofs make use of the Euclidean ball $\mathcal{R}_{s}$ defined in
equation (29). At a high level, the HMC chain has bounded gradient inside the
ball $\mathcal{R}_{s}$ for a suitable choice of $s$, and the gradient of the
log-density gets too large outside such a ball making the chain unstable in
that region. However, since the target distribution has low mass in that
region, the chain’s visit to the region outside the ball is a rare event and
thus we can focus on the chain’s behavior inside the ball to analyze its
mixing time.

In the next lemma, we state the overlap bounds for the transition
distributions of the HMC chain. For a fixed univeral constant $c$, we require

$\displaystyle K^{2}\eta^{2}$
$\displaystyle\leq\frac{1}{4\max\left\\{d^{\frac{1}{2}}L,d^{\frac{2}{3}}L_{\text{H}}^{\frac{2}{3}}\right\\}},\qquad\mbox{and}$
(26a) $\displaystyle\eta^{2}$
$\displaystyle\leq\frac{1}{cL}\min\left\\{\frac{1}{K^{2}},\frac{1}{Kd^{\frac{1}{2}}},\frac{1}{K^{\frac{2}{3}}d^{\frac{1}{3}}\left(\frac{M^{2}}{L}\right)^{\frac{1}{3}}},\frac{1}{K\frac{M}{L^{\frac{1}{2}}}},\frac{1}{K^{\frac{2}{3}}d}\frac{L}{L_{\text{H}}^{\frac{2}{3}}},\frac{1}{K^{\frac{4}{3}}\frac{M}{L^{\frac{1}{2}}}}\left(\frac{L}{L_{\text{H}}^{\frac{2}{3}}}\right)^{\frac{1}{2}}\right\\}.$
(26b)

###### Lemma 6

Consider a $(L,L_{\text{H}},s,\psi_{\mathfrak{a}},M)$-regular target
distribution (cf. Assumption (A)) with $\Omega$ the convex measurable set
satisfying (10e). Then with the parameters $\left(K,\eta\right)$ satisfying
$K\eta\leq\frac{1}{4L}$ and condition (26a), the HMC-$\left(K,\eta\right)$
chain satisfies

$\displaystyle\sup_{\left\|q_{0}-{\tilde{q}}_{0}\right\|_{2}\leq\frac{K\eta}{4}}d_{\text{\tiny{TV}}}\left(\mathcal{P}_{q_{0}},\mathcal{P}_{{\tilde{q}}_{0}}\right)\leq\frac{1}{2}.$
(27a) If, in addition, condition (26b) holds, then we have
$\displaystyle\sup_{x\in\Omega}d_{\text{\tiny{TV}}}\left(\mathcal{P}_{x},\mathcal{T}^{\text{before-
lazy}}_{x}\right)\leq\frac{1}{8}.$ (27b)

See Appendix 5.3 for the proof.

Lemma 6 is crucial to the analysis of HMC as it enables us to apply the
conductance profile based bounds discussed in Section 3.3. It reveals two
important properties of the Metropolized HMC. First, from equation (27a), we
see that proposal distributions of HMC at two different points are close if
the two points are close. This is proved by controlling the KL-divergence of
the two proposal distributions of HMC via change of variable formula. Second,
equation (27b) shows that the accept-reject step of HMC is well behaved inside
$\Omega$ provided the gradient is bounded by $M$.

###  5.2 Proof of Theorem 1

We are now equipped to prove our main theorem. In order to prove Theorem 1, we
begin by using Lemma 4 and Lemma 6 to derive an explicit bound for on the HMC
conductance profile. Given the assumptions of Theorem 1, conditions (26a) and
(26b) hold, enabling us to invoke Lemma 6 in the proof.

Define the function $\Psi_{\Omega}:[0,1]\mapsto{}_{+}$ as

$\displaystyle\Psi_{\Omega}(v)$
$\displaystyle=\begin{cases}\displaystyle\frac{1}{32}\cdot\min\left\\{1,\frac{K\eta}{64\psi_{\mathfrak{a}}}\log^{\mathfrak{a}}\left(\frac{1}{v}\right)\right\\}&\mbox{if
$v\in\left[0,\frac{1-s}{2}\right]$.}\\\\[8.53581pt]
\displaystyle\frac{K\eta}{2048\psi_{\mathfrak{a}}},&\mbox{if
$v\in\left(\frac{1-s}{2},1\right]$.}\end{cases}$ (28)

This function acts as a lower bound on the truncated conductance profile.
Define the Euclidean ball

$\displaystyle\mathcal{R}_{s}=\mathbb{B}\left(x^{\star},r(s)\sqrt{\frac{d}{m}}\right),$
(29)

and consider a pair $(x,y)\in\mathcal{R}_{s}$ such that
$\left\|x-y\right\|_{2}\leq\frac{1}{4}K\eta$. Invoking the decomposition (25)
and applying triangle inequality for $\zeta$-lazy HMC, we have

$\displaystyle
d_{\text{\tiny{TV}}}\left(\mathcal{T}_{x},\mathcal{T}_{y}\right)$
$\displaystyle\leq\zeta+\left(1-\zeta\right)d_{\text{\tiny{TV}}}\left(\mathcal{T}^{\text{before-
lazy}}_{x},\mathcal{T}^{\text{before-lazy}}_{y}\right)$
$\displaystyle\leq\zeta+\left(1-\zeta\right)\left(d_{\text{\tiny{TV}}}\left(\mathcal{T}^{\text{before-
lazy}}_{x},\mathcal{P}_{y}\right)+d_{\text{\tiny{TV}}}\left(\mathcal{P}_{x},\mathcal{P}_{y}\right)+d_{\text{\tiny{TV}}}\left(\mathcal{P}_{x},\mathcal{T}^{\text{before-
lazy}}_{y}\right)\right)$
$\displaystyle\stackrel{{\scriptstyle(i)}}{{\leq}}\zeta+\left(1-\zeta\right)\left(\frac{1}{4}+\frac{1}{2}+\frac{1}{4}\right)$
$\displaystyle=1-\frac{1-\zeta}{4},$

where step (i) follows from the bounds (27a) and (27b) from Lemma 6. For
$\zeta=\frac{1}{2}$, substituting $\omega=\frac{1}{8}$,
$\Delta=\frac{1}{4}K\eta$ and the convex set $\Omega=\mathcal{R}_{s}$ into
Lemma 4, we obtain that

$\displaystyle\Phi_{\Omega}(v)\geq\frac{1}{32}\cdot\min\left\\{1,\frac{K\eta}{64\psi_{\mathfrak{a}}}\log^{\mathfrak{a}}\left(1+\frac{1}{v}\right)\right\\},\quad\text{for
}v\in\left[0,\frac{1-s}{2}\right].$

Here $\mathfrak{a}$ equals to $\frac{1}{2}$ or $0$, depending on the
assumption (10d). By the definition of the truncated conductance profile (15),
we have that
$\tilde{\Phi}_{\Omega}(v)\geq\frac{K\eta}{2048\psi_{\mathfrak{a}}}$ for
$v\in\left[\frac{1-s}{2},1\right]$. As a consequence, $\Psi_{\Omega}$ is
effectively a lower bound on the truncated conductance profile. Note that the
assumption (A) ensures the existence of $\Omega$ such that
$\Pi^{*}(\Omega)\geq 1-s$ for $s=\frac{\epsilon^{2}}{3\beta^{2}}$. Putting the
pieces together and applying Lemma 3 with the convex set $\Omega$ concludes
the proof of the theorem.

###  5.3 Proof of Lemma 6

In this subsection, we prove the two main claims (27a) and (27b) in Lemma 6.
Before going into the claims, we first provide several convenient properties
about the HMC proposal.

####  5.3.1 Properties of the HMC proposal

Recall the Hamiltonian Monte Carlo (HMC) with leapfrog integrator (8c). Using
an induction argument, we find that the final states in one iteration of $K$
steps of the HMC chain, denoted by $q_{K}$ and $p_{K}$ satisfy

$\displaystyle p_{K}$ $\displaystyle=p_{0}-\frac{\eta}{2}{\nabla
f}(q_{0})-\sum_{j=1}^{K-1}{\nabla f}(q_{j})-\frac{\eta}{2}{\nabla f}(q_{K}),$
(30a) $\displaystyle\quad\text{and}\quad q_{K}$ $\displaystyle=q_{0}+K\eta
p_{0}-\frac{K\eta^{2}}{2}{\nabla
f}(q_{0})-\eta^{2}\sum_{j=1}^{K-1}\left(K-j\right){\nabla f}(q_{j}).$ (30b) It
is easy to see that for $k\in\left[K\right]$, $q_{k}$ can be seen as a
function of the initial state $q_{0}$ and $p_{0}$. We denote this function as
the _forward mapping_ $F$, $\displaystyle
q_{k}=:F_{k}(p_{0},q_{0})\quad\text{and}\quad
q_{K}=:F_{K}(p_{0},q_{0})=:F\left(p_{0},q_{0}\right)$ (30c) where we
introduced the simpler notation $F:\,=F_{K}$ for the final iterate. The
forward mappings $F_{k}$ and $F$ are deterministic functions that only depends
on the gradient ${\nabla f}$, the number of leapfrog updates $K$ and the step
size $\eta$.

Denote $\mathbf{J}_{x}F$ as the Jacobian matrix of the forward mapping $F$
with respect to the first variable. By definition, it satisfies

$\displaystyle\left[\mathbf{J}_{x}F(x,q_{0})\right]_{ij}=\frac{\partial}{\partial
x_{j}}\left[F(x,q_{0})\right]_{i},\quad\text{for all}\quad
i,j\in\left[d\right].$ (30d)

Similarly, denote $\mathbf{J}_{y}F$ as the Jacobian matrix of the forward
mapping $F$ with respect to the second variable. The following lemma
characterizes the eigenvalues of the Jacobian $\mathbf{J}_{x}F$.

###### Lemma 7

Suppose the log density $f$ is $L$-smooth. For the number of leapfrog steps
and step-size satisfying $K^{2}\eta^{2}\leq\frac{1}{4L}$, we have

$\displaystyle\left|\\!\left|\\!\left|K\eta\mathbb{I}_{d}-\mathbf{J}_{x}F(x,y)\right|\\!\right|\\!\right|_{{2}}\leq\frac{1}{8}K\eta,\quad\text{for
all}\quad x,y\in\mathcal{X}\text{ and }i\in\left[d\right].$

Also all eigenvalues of $\mathbf{J}_{x}F(x,y)$ have absolute value greater or
equal to $\frac{7}{8}K\eta$.

See Appendix A.3.1 for the proof.

Since the Jacobian is invertible for $K^{2}\eta^{2}\leq\frac{1}{4L}$, we can
define the inverse function of $F$ with respect to the first variable as the
backward mapping $G$. We have

$\displaystyle F\left(G(x,y),y\right)=x,\quad\text{for all}\quad
x,y\in\mathcal{X}.$ (31)

Moreover as a direct consequence of Lemma 7, we obtain that the magnitude of
the eigenvalues of the Jacobian matrix $\mathbf{J}_{x}G(x,y)$ lies in the
interval $\left[\frac{8}{9K\eta},\frac{8}{7K\eta}\right]$. In the next lemma,
we state another set of bounds on different Jacobian matrices:

###### Lemma 8

Suppose the log density $f$ is $L$-smooth. For the number of leapfrog steps
and step-size satisfying $K^{2}\eta^{2}\leq\frac{1}{4L}$, we have

$\displaystyle\left|\\!\left|\\!\left|\mathbf{J}_{y}G(x,y)\right|\\!\right|\\!\right|_{{2}}$
$\displaystyle\leq\frac{4}{3K\eta},\quad\text{for all}\quad
x,y\in\mathcal{X},\quad\text{and}\quad$ (32a)
$\displaystyle\left|\\!\left|\\!\left|\frac{\partial F_{k}(G(x,y),y)}{\partial
y}\right|\\!\right|\\!\right|_{{2}}$ $\displaystyle\leq 3,\quad\quad\text{for
all}\quad k\in\left[K\right].$ (32b)

See Appendix A.3.2 for the proof.

Next, we would like to obtain a bound on the quantity
$\frac{\partial\log\det\mathbf{J}_{x}G(x,q_{0})}{\partial y}$. Applying the
chain rule, we find that

$\displaystyle\frac{\partial\log\det\mathbf{J}_{x}G(x,q_{0})}{\partial
y}=\begin{bmatrix}\operatorname{\text{trace}}\left([\mathbf{J}_{x}G(x,q_{0})]^{-1}\mathbf{J}_{xy_{1}}G(x,q_{0})\right)\\\
\vdots\\\
\operatorname{\text{trace}}\left([\mathbf{J}_{x}G(x,q_{0})]^{-1}\mathbf{J}_{xy_{d}}G(x,q_{0})\right)\end{bmatrix}.$
(33)

Here $\mathbf{J}_{xy}G(x,q_{0})$ is a third order tensor and we use
$\mathbf{J}_{xy_{l}}G(x,q_{0})$ to denote the matrix corresponding to the
$l$-th slice of the tensor which satisfies

$\displaystyle\left[\mathbf{J}_{xy_{l}}G(x,q_{0})\right]_{ij}=\frac{\partial\partial}{\partial
x_{j}y_{l}}\left[F(x,q_{0})\right]_{i},\ \quad\text{for all}\quad
i,j,l\in\left[d\right].$

###### Lemma 9

Suppose the log density $f$ is $L$-smooth and $L_{\text{H}}$-Hessian
Lipschitz. For the number of leapfrog steps and step-size satisfying
$K^{2}\eta^{2}\leq\frac{1}{4L}$, we have

$\displaystyle\left\|\frac{\partial\log\det\mathbf{J}_{x}G(x,q_{0})}{\partial
y}\right\|_{2}=\left\|\begin{bmatrix}\operatorname{\text{trace}}\left([\mathbf{J}_{x}G(x,q_{0})]^{-1}\mathbf{J}_{xy_{1}}G(x,q_{0})\right)\\\
\vdots\\\
\operatorname{\text{trace}}\left([\mathbf{J}_{x}G(x,q_{0})]^{-1}\mathbf{J}_{xy_{d}}G(x,q_{0})\right)\end{bmatrix}\right\|_{2}\leq
2dK^{2}\eta^{2}L_{\text{H}}.$

See Appendix A.3.3 for the proof.

As a direct consequence of the equation (30b) at $k$-th step of leapfrog
updates, we obtain the following two bounds for the difference between
successive $F_{k}$ terms that come in handy later in our proofs.

###### Lemma 10

Suppose that the log density $f$ is $L$-smooth. For the number of leapfrog
steps and step-size satisfying $K^{2}\eta^{2}\leq\frac{1}{4L}$, we have

$\displaystyle\left\|F_{k}(p_{0},q_{0})-q_{0}\right\|_{2}$ $\displaystyle\leq
2k\eta\left\|p_{0}\right\|_{2}+2k^{2}\eta^{2}\left\|{\nabla
f}(q_{0})\right\|_{2}\ \ \text{for }k\in\left[K\right],\ \text{and}$ (34a)
$\displaystyle\left\|F_{k+1}(p_{0},q_{0})-F_{k}(p_{0},q_{0})\right\|_{2}$
$\displaystyle\leq 2\eta\left\|p_{0}\right\|_{2}+2(k+1)\eta^{2}\left\|{\nabla
f}(q_{0})\right\|_{2}\ \text{for }k\in\left[K-1\right].$ (34b)

See Appendix A.3.4 for the proof.

We now turn to the proof the two claims in Lemma 6. Note that the claim (27a)
states that the proposal distributions at two close points are close; the
claim (27b) states that the proposal distribution and the transition
distribution are close.

####  5.3.2 Proof of claim (27a) in Lemma 6

In order to bound the distance between proposal distributions of nearby
points, we prove the following stronger claim: For a $L$-smooth
$L_{\text{H}}$-Hessian-Lipschitz target distribution, the proposal
distribution of the HMC algorithm with step size $\eta$ and leapfrog steps $K$
such that $K\eta\leq\frac{1}{4L}$ satisfies

$\displaystyle
d_{\text{\tiny{TV}}}\left(\mathcal{P}_{q_{0}},\mathcal{P}_{{\tilde{q}}_{0}}\right)\leq\left(\frac{2\left\|q_{0}-{\tilde{q}}_{0}\right\|_{2}^{2}}{K^{2}\eta^{2}}+3\sqrt{d}K\eta
L\left\|q_{0}-{\tilde{q}}_{0}\right\|_{2}+4dK^{2}\eta^{2}L_{\text{H}}\left\|q_{0}-{\tilde{q}}_{0}\right\|_{2}\right)^{1/2},$
(35)

for all $q_{0},{\tilde{q}}_{0}\in{}^{d}$. Then for any two points
$q_{0},{\tilde{q}}_{0}$ such that
$\left\|q_{0}-{\tilde{q}}_{0}\right\|_{2}\leq\frac{1}{4}K\eta$, under the
condition (26a), i.e.,
$K^{2}\eta^{2}\leq\frac{1}{4\max\left\\{d^{\frac{1}{2}}L,d^{\frac{2}{3}}L_{\text{H}}^{\frac{2}{3}}\right\\}}$,
we have

$\displaystyle
d_{\text{\tiny{TV}}}\left(\mathcal{P}_{q_{0}},\mathcal{P}_{{\tilde{q}}_{0}}\right)\leq\left(\frac{1}{8}+\frac{3}{64}+\frac{1}{64}\right)^{1/2}\leq\frac{1}{2},$

and the claim (27a) follows.

The proof of claim (35) involves the following steps: (1) we make use of the
update rules (30b) and change of variable formula to obtain an expression for
the density of $q_{n}$ in terms of $q_{0}$, (2) then we use Pinsker’s
inequality and derive expressions for the KL-divergence between the two
proposal distributions, and (3) finally, we upper bound the KL-divergence
between the two distributions using different properties of the forward
mapping $F$ from Appendix 5.3.1.

According to the update rule (30b), the proposals from two initial points
$q_{0}$ and ${\tilde{q}}_{0}$ satisfy respectively

$\displaystyle q_{K}=F(p_{0},q_{0}),\quad\text{ and
}\quad{\tilde{q}}_{K}=F({\tilde{p}}_{0},{\tilde{q}}_{0}),$

where $p_{0}$ and ${\tilde{p}}_{0}$ are independent random variable from
Gaussian distribution $\mathcal{N}(0,\mathbb{I}_{d})$.

Denote $\rho_{q_{0}}$ as the density function of the proposal distribution
$\mathcal{P}_{q_{0}}$. For two different initial points $q_{0}$ and
${\tilde{q}}_{0}$, the goal is to bound the total variation distance between
the two proposal distribution, which is by definition

$\displaystyle
d_{\text{\tiny{TV}}}\left(\mathcal{P}_{q_{0}},\mathcal{P}_{{\tilde{q}}_{0}}\right)=\frac{1}{2}\int_{x\in\mathcal{X}}\left|\rho_{q_{0}}(x)-\rho_{{\tilde{q}}_{0}}(x)\right|dx.$
(36)

Given $q_{0}$ fixed, the random variable $q_{K}$ can be seen as a
transformation of the Gaussian random variable $p_{0}$ through the function
$F(\cdot,q_{0})$. When $F$ is invertible, we can use the change of variable
formula to obtain an explicit expression of the density $\rho_{q_{0}}$:

$\displaystyle\rho_{q_{0}}(x)=\varphi\left(G(x,q_{0})\right)\det\left(\mathbf{J}_{x}G(x,q_{0})\right),$
(37)

where $\varphi$ is the density of the standard Gaussian distribution
$\mathcal{N}(0,\mathbb{I}_{d})$. Note that even though explicit, directly
bounding the total variation distance (36) using the complicated density
expression (37) is difficult. We first use Pinsker’s inequality (Cover and
Thomas, 1991) to give an upper bound of the total variance distance in terms
of KL-divergence

$\displaystyle
d_{\text{\tiny{TV}}}\left(\mathcal{P}_{q_{0}},\mathcal{P}_{{\tilde{q}}_{0}}\right)\leq\sqrt{2\text{KL}(\mathcal{P}_{q_{0}}\;\|\;\mathcal{P}_{{\tilde{q}}_{0}})},$
(38)

and then upper bound the KL-divergence. Plugging the density (37) into the KL-
divergence formula, we obtain that

$\displaystyle\text{KL}(\mathcal{P}_{q_{0}}\;\|\;\mathcal{P}_{{\tilde{q}}_{0}})$
$\displaystyle=\int_{{}^{d}}\rho_{q_{0}}(x)\log\left(\frac{\rho_{q_{0}}(x)}{\rho_{{\tilde{q}}_{0}}(x)}\right)dx$
$\displaystyle=\int_{{}^{d}}\rho_{q_{0}}(x)\left[\log\left(\frac{\varphi\left(G(x,q_{0})\right)}{\varphi\left(G(x,{\tilde{q}}_{0})\right)}\right)+\log\det\mathbf{J}_{x}G(x,q_{0})-\log\det\mathbf{J}_{x}G(x,{\tilde{q}}_{0})\right]dx$
$\displaystyle=\underbrace{\int_{{}^{d}}\rho_{q_{0}}(x)\left[\frac{1}{2}\left(-\left\|G(x,q_{0})\right\|_{2}^{2}+\left\|G(x,{\tilde{q}}_{0})\right\|_{2}^{2}\right)\right]dx}_{T_{1}}$
$\displaystyle\quad\quad\quad\quad+\underbrace{\int_{{}^{d}}\rho_{q_{0}}(x)\left[\log\det\mathbf{J}_{x}G(x,q_{0})-\log\det\mathbf{J}_{x}G(x,{\tilde{q}}_{0})\right]dx}_{T_{2}}$
(39)

We claim the following bounds on the terms $T_{1}$ and $T_{2}$:

$\displaystyle\left|T_{1}\right|$
$\displaystyle\leq\frac{8}{9}\frac{\left\|q_{0}-{\tilde{q}}_{0}\right\|_{2}^{2}}{K^{2}\eta^{2}}+\frac{3}{2}\sqrt{d}K\eta
L\left\|q_{0}-{\tilde{q}}_{0}\right\|_{2},\quad\text{and}$ (40a)
$\displaystyle\left|T_{2}\right|$ $\displaystyle\leq
2dK^{2}\eta^{2}L_{\text{H}}\left\|q_{0}-{\tilde{q}}_{0}\right\|_{2},$ (40b)

where the bound on $T_{2}$ follows readily from Lemma 9:

$\displaystyle\left|T_{2}\right|$
$\displaystyle=\left|\int\rho_{q_{0}}(x)\left[\log\det\mathbf{J}_{x}G(x,q_{0})-\log\det\mathbf{J}_{x}G(x,{\tilde{q}}_{0})\right]dx\right|$
$\displaystyle\leq\left\|\frac{\partial\log\det\mathbf{J}_{x}G(x,q_{0})}{\partial
y}\right\|_{2}\left\|q_{0}-{\tilde{q}}_{0}\right\|_{2}$ $\displaystyle\leq
2dK^{2}\eta^{2}L_{\text{H}}\left\|q_{0}-{\tilde{q}}_{0}\right\|_{2}.$ (41)

Putting together the inequalities (38), (5.3.2), (40a) and (40b) yields the
claim (35).

It remains to prove the bound (40a) on $T_{1}$.

##### Proof of claim (40a):

For the term $T_{1}$, we observe that

$\displaystyle\frac{1}{2}\left(\left\|G(x,{\tilde{q}}_{0})\right\|_{2}^{2}-\left\|G(x,q_{0})\right\|_{2}^{2}\right)=\frac{1}{2}\left\|G(x,q_{0})-G(x,{\tilde{q}}_{0})\right\|_{2}^{2}-\left(G(x,q_{0})-G(x,{\tilde{q}}_{0})\right)^{\top}G(x,q_{0}).$

The first term on the RHS can be bounded via the Jacobian of $G$ with respect
to the second variable. Applying the bound (32a) from Lemma 8, we find that

$\displaystyle\left\|G(x,q_{0})-G(x,{\tilde{q}}_{0})\right\|_{2}\leq\left|\\!\left|\\!\left|\mathbf{J}_{y}G(x,y)\right|\\!\right|\\!\right|_{{2}}\left\|q_{0}-{\tilde{q}}_{0})\right\|_{2}\leq\frac{4}{3K\eta}\left\|q_{0}-{\tilde{q}}_{0})\right\|_{2}.$
(42)

For the second part, we claim that there exists a deterministic function $C$
of $q_{0}$ and ${\tilde{q}}_{0}$ and independent of $x$, such that

$\displaystyle\left\|G(x,q_{0})-G(x,{\tilde{q}}_{0})-C(q_{0},{\tilde{q}}_{0})\right\|_{2}\leq\frac{3}{2}K\eta
L\left\|q_{0}-{\tilde{q}}_{0}\right\|_{2}.$ (43)

Assuming the claim (43) as given at the moment, we can further decompose the
second part of $T_{1}$ into two parts:

$\displaystyle\left(G(x,q_{0})-G(x,{\tilde{q}}_{0})\right)^{\top}G(x,q_{0})=\left(G(x,q_{0})-G(x,{\tilde{q}}_{0})-C(q_{0},{\tilde{q}}_{0})\right)^{\top}G(x,q_{0})+C(q_{0},{\tilde{q}}_{0})^{\top}G(x,q_{0})$
(44)

Applying change of variables along with equation (37), we find that

$\displaystyle\int\rho_{q_{0}}(x)G(x,q_{0})dx=\int\varphi(x)xdx=0.$

Furthermore, we also have

$\displaystyle\int_{x\in\mathcal{X}}\rho_{q_{0}}(x)\left\|G(x,q_{0})\right\|_{2}dx$
$\displaystyle=\int_{x\in\mathcal{X}}\varphi(x)\left\|x\right\|_{2}dx$
$\displaystyle\stackrel{{\scriptstyle(i)}}{{\leq}}\left[\left(\int_{x\in\mathcal{X}}\varphi(x)\left\|x\right\|_{2}^{2}dx\right){\left(\int_{x\in\mathcal{X}}\varphi(x)dx\right)}\right]^{1/2}=\sqrt{d},$

where step (i) follows from Cauchy-Schwarz’s inequality. Combining the
inequalities (42), (43) and (44) together, we obtain the following bound on
term $T_{1}$:

$\displaystyle\left|T_{1}\right|$
$\displaystyle=\left|\int\rho_{q_{0}}(x)\left[-\frac{1}{2}\left\|G(x,q_{0})\right\|_{2}^{2}+\frac{1}{2}\left\|G(x,{\tilde{q}}_{0})\right\|_{2}^{2}\right]dx\right|$
$\displaystyle\leq\frac{1}{2}\left|\int\rho_{q_{0}}(x)\left\|G(x,q_{0})-G(x,{\tilde{q}}_{0})\right\|_{2}^{2}dx\right|$
$\displaystyle\qquad\qquad\qquad+\left|\int\rho_{q_{0}}(x)\left\|G(x,q_{0})-G(x,{\tilde{q}}_{0})-C(q_{0},{\tilde{q}}_{0})\right\|_{2}\left\|G(x,q_{0})\right\|_{2}dx\right|$
$\displaystyle\leq\frac{8}{9}\frac{\left\|q_{0}-{\tilde{q}}_{0}\right\|_{2}^{2}}{K^{2}\eta^{2}}+\frac{3}{2}\sqrt{d}K\eta\left\|q_{0}-{\tilde{q}}_{0}\right\|_{2},$
(45)

which yields the claimed bound on $T_{1}$.

We now prove our earlier claim (43).

##### Proof of claim (43):

For any pair of states $q_{0}$ and ${\tilde{q}}_{0}$, invoking the definition
(31) of the map $G(x,\cdot)$, we obtain the following implicit equations:

$\displaystyle x$ $\displaystyle=q_{0}+K\eta
G(x,q_{0})-K\frac{\eta^{2}}{2}{\nabla
f}(q_{0})-\eta^{2}\sum_{j=1}^{K-1}(K-j){\nabla
f}(F_{j}(G(x,q_{0}),q_{0})),\quad\text{and}$ $\displaystyle x$
$\displaystyle={\tilde{q}}_{0}+K\eta
G(x,{\tilde{q}}_{0})-K\frac{\eta^{2}}{2}{\nabla
f}({\tilde{q}}_{0})-\eta^{2}\sum_{j=1}^{K-1}(K-j){\nabla
f}(F_{j}(G(x,{\tilde{q}}_{0}),{\tilde{q}}_{0})).$

Taking the difference between the two equations above, we obtain

$\displaystyle G(x,q_{0})-G(x,{\tilde{q}}_{0})$
$\displaystyle-\frac{q_{0}-{\tilde{q}}_{0}}{K\eta}-\frac{\eta}{2}\left({\nabla
f}(q_{0})-{\nabla f}({\tilde{q}}_{0})\right)$
$\displaystyle=\frac{\eta^{2}}{K\eta}\sum_{k=1}^{K-1}(K-j)\left({\nabla
f}(F_{k}(G(x,q_{0}),q_{0}))-{\nabla
f}(F_{k}(G(x,{\tilde{q}}_{0}),{\tilde{q}}_{0}))\right).$

Applying $L$-smoothness of $f$ along with the bound (32b) from Lemma 8, we
find that

$\displaystyle\left\|{\nabla f}(F_{k}(G(x,q_{0}),q_{0}))-{\nabla
f}(F_{k}(G(x,{\tilde{q}}_{0}),{\tilde{q}}_{0}))\right\|_{2}$
$\displaystyle\leq L\left|\\!\left|\\!\left|\frac{\partial
F_{k}(G(x,y),y)}{\partial
y}\right|\\!\right|\\!\right|_{{2}}\left\|q_{0}-{\tilde{q}}_{0}\right\|_{2}$
$\displaystyle\leq 3L\left\|q_{0}-{\tilde{q}}_{0}\right\|_{2}.$

Putting the pieces together, we find that

$\displaystyle\left\|G(x,q_{0})-G(x,{\tilde{q}}_{0})-\frac{q_{0}-{\tilde{q}}_{0}}{K\eta}-\frac{1}{2}\left({\nabla
f}(q_{0})-{\nabla f}({\tilde{q}}_{0})\right)\right\|_{2}\leq\frac{3K\eta
L}{2}\left\|q_{0}-{\tilde{q}}_{0}\right\|_{2},$

which yields the claim (43).

####  5.3.3 Proof of claim (27b) in Lemma 6

We now bound the distance between the one-step proposal distribution
$\mathcal{P}_{x}$ at point $x$ and the one-step transition distribution
$\mathcal{T}^{\text{before-lazy}}_{x}$ at $x$ obtained after performing the
accept-reject step (and no lazy step). Using equation (30a), we define the
forward mapping $E$ for the variable $p_{K}$ as follows

$\displaystyle p_{K}=E(p_{0},q_{0}):\,=p_{0}-\frac{\eta}{2}{\nabla
f}(q_{0})-\eta\sum_{j=1}^{K-1}{\nabla f}(q_{j})-\frac{\eta}{2}{\nabla
f}(q_{K}).$

Consequently, the probability of staying at $x$ is given by

$\displaystyle\mathcal{T}^{\text{before-
lazy}}_{x}(\left\\{x\right\\})=1-\int_{\mathcal{X}}\min\left\\{1,\frac{\exp(-\mathcal{H}(E(z,x),F(z,x)))}{\exp(-\mathcal{H}(z,x))}\right\\}\varphi_{x}(z)dz,$

where the Hamiltonian
$\mathcal{H}(q,p)=f(q)+\frac{1}{2}\left\|p\right\|_{2}^{2}$ was defined in
equation (7). As a result, the TV-distance between the proposal and transition
distribution is given by

$\displaystyle
d_{\text{\tiny{TV}}}\left(\mathcal{P}_{x},\mathcal{T}^{\text{before-
lazy}}_{x}\right)$
$\displaystyle=1-\int_{\mathcal{X}}\min\left\\{1,\frac{\exp(-\mathcal{H}(E(z,x),F(z,x)))}{\exp(-\mathcal{H}(z,x))}\right\\}\varphi_{x}(z)dz$
$\displaystyle=1-{\mathbb{E}}_{z\sim\mathcal{N}(0,\mathbb{I}_{d})}\left[\min\left\\{1,\frac{\exp(-\mathcal{H}(E(z,x),F(z,x)))}{\exp(-\mathcal{H}(z,x))}\right\\}\right].$
(46)

An application of Markov’s inequality yields that

$\displaystyle{\mathbb{E}}_{z\sim\mathcal{N}(0,\mathbb{I}_{d})}\left[\min\left\\{1,\frac{\exp(-\mathcal{H}(E(z,x),F(z,x)))}{\exp(-\mathcal{H}(z,x))}\right\\}\right]$
$\displaystyle\qquad\geq\alpha{\mathbb{P}}_{z\sim\mathcal{N}(0,\mathbb{I}_{d})}\left[\frac{\exp(-\mathcal{H}(E(z,x),F(z,x)))}{\exp(-\mathcal{H}(z,x))}\geq\alpha\right],$
(47)

for any $\alpha\in(0,1]$. Thus, to bound the distance
$d_{\text{\tiny{TV}}}\left(\mathcal{P}_{x},\mathcal{T}^{\text{before-
lazy}}_{x}\right)$, it suffices to derive a high probability lower bound on
the ratio ${\exp(-\mathcal{H}(E(z,x),F(z,x)))}/{\exp(-\mathcal{H}(z,x))}$ when
$z\sim\mathcal{N}(0,\mathbb{I}_{d})$.

We now derive a lower bound on the following quantity:

$\displaystyle\exp\left(-f(F(p_{0},q_{0}))+f(q_{0})-\frac{1}{2}\left\|E(p_{0},q_{0})\right\|_{2}^{2}+\frac{1}{2}\left\|p_{0}\right\|_{2}^{2}\right),\quad\mbox{when
$p_{0}\sim\mathcal{N}(0,\mathbb{I}_{d})$.}$

We derive the bounds on the two terms $-f(F(p_{0},q_{0}))+f(q_{0})$ and
$\left\|E(p_{0},q_{0})\right\|_{2}^{2}$ separately.

Observe that

$\displaystyle
f(F(p_{0},q_{0}))-f(q_{0})=\sum_{j=0}^{K-1}\left[f(F_{j+1}(p_{0},q_{0}))-f(F_{j}(p_{0},q_{0}))\right].$

The intuition is that it is better to apply Taylor expansion on closer points.
Applying the third order Taylor expansion and using the smoothness assumptions
(10a) and (10c) for the function $f$, we obtain

$\displaystyle f(x)-f(y)\leq\frac{(x-y)^{\top}}{2}\left({\nabla f}(x)+{\nabla
f}(y)\right)+L_{\text{H}}\left\|x-y\right\|_{2}^{3}.$

For the indices $j\in\left\\{0,\ldots,K-1\right\\}$, using $F_{j}$ as the
shorthand for $F_{j}(p_{0},q_{0})$, we find that

$\displaystyle f(F_{j+1})-f(F_{j})$
$\displaystyle\leq\frac{(F_{j+1}-F_{j})^{\top}}{2}\left({\nabla
f}(F_{j+1})+{\nabla
f}(F_{j})\right)+L_{\text{H}}\left\|F_{j+1}-F_{j}\right\|_{2}^{3}$
$\displaystyle=\frac{1}{2}\eta p_{0}^{\top}\left({\nabla f}(F_{j+1})+{\nabla
f}(F_{j})\right)$
$\displaystyle\quad-\frac{\eta^{2}}{2}\bigg{[}{\frac{1}{2}{\nabla
f}(p_{0})+\sum_{k=1}^{j}{\nabla f}(F_{k})}\bigg{]}^{\top}\left({\nabla
f}(F_{j+1})+{\nabla
f}(F_{j})\right)+L_{\text{H}}\left\|F_{j+1}-F_{j}\right\|_{2}^{3},$ (48)

where the last equality follows by definition (30c) of the operator $F_{j}$.

Now to bound the term $E(p_{0},q_{0})$, we observe that

$\displaystyle\frac{\left\|E(p_{0},q_{0})\right\|_{2}^{2}}{2}$
$\displaystyle=\frac{\left\|p_{0}-\frac{\eta}{2}{\nabla
f}(q_{0})-\eta\sum_{j=1}^{K-1}{\nabla f}(F_{j})-\frac{\eta}{2}{\nabla
f}(F_{K})\right\|_{2}^{2}}{2}$
$\displaystyle=\frac{\left\|p_{0}\right\|_{2}^{2}}{2}-\eta
p_{0}^{\top}\bigg{(}{\frac{1}{2}{\nabla f}(q_{0})+\sum_{j=1}^{K-1}{\nabla
f}(F_{j})+\frac{1}{2}{\nabla f}(F_{K})}\bigg{)}$
$\displaystyle\qquad\qquad\qquad+\frac{\eta^{2}}{2}\big{\|}{\frac{1}{2}{\nabla
f}(q_{0})+\sum_{j=1}^{K-1}{\nabla f}(F_{j})+\frac{1}{2}{\nabla
f}(F_{K})}\big{\|}_{2}^{2}.$ (49)

Putting the equations (5.3.3) and (49) together leads to cancellation of many
gradient terms and we obtain

$\displaystyle-f(F(p_{0},q_{0}))+f(q_{0})-\frac{1}{2}\left\|E(p_{0},q_{0})\right\|_{2}^{2}+\frac{1}{2}\left\|p_{0}\right\|_{2}^{2}$
$\displaystyle\quad\geq\frac{\eta^{2}}{8}\left({\nabla f}(q_{0})-{\nabla
f}(F_{K})\right)^{\top}\left({\nabla f}(q_{0})+{\nabla
f}(F_{K})\right)-L_{\text{H}}\sum_{j=0}^{K-1}\left\|F_{j+1}-F_{j}\right\|_{2}^{3}$
$\displaystyle\quad\geq-\frac{\eta^{2}L}{4}\left\|q_{0}-F(p_{0},q_{0})\right\|_{2}\left\|{\nabla
f}(q_{0})\right\|_{2}-\frac{\eta^{2}L^{2}}{2}\left\|q_{0}-F(p_{0},q_{0})\right\|_{2}^{2}-L_{\text{H}}\sum_{j=0}^{K-1}\left\|F_{j+1}-F_{j}\right\|_{2}^{3}$
(50)

The last inequality uses the smoothness condition (10a) for the function $f$.
Plugging the bounds (34a) and (34b) in equation (5.3.3), we obtain a lower
bound that only depends on $\left\|p_{0}\right\|_{2}$ and $\left\|{\nabla
f}(q_{0})\right\|_{2}$:

$\displaystyle\text{RHS of
}\eqref{eq:accept_rate_lower_bound}\geq-2K^{2}\eta^{4}L^{2}\left\|p_{0}\right\|_{2}^{2}-2K\eta^{3}L\left\|p_{0}\right\|_{2}\left\|{\nabla
f}(q_{0})\right\|_{2}-2K^{2}\eta^{4}L\left\|{\nabla f}(q_{0})\right\|_{2}^{2}$
$\displaystyle-
L_{\text{H}}\left(32K\eta^{3}\left\|p_{0}\right\|_{2}^{3}+8K^{4}\eta^{6}\left\|{\nabla
f}(q_{0})\right\|_{2}^{3}\right).$ (51)

According to assumption (A), we have bounded gradient in the convex set
$\Omega$. For any $x\in\Omega$, we have $\left\|{\nabla f}(x)\right\|_{2}\leq
M$. Standard Chi-squared tail bounds imply that

$\displaystyle{\mathbb{P}}\left[\left\|p_{0}\right\|_{2}^{2}\leq
d\alpha_{1}\right]\geq 1-\frac{1}{16},\quad\mbox{for
$\alpha_{1}=1+2\sqrt{\log(16)}+2\log(16)$}.$ (52)

Plugging the gradient bound and the bound (52) into equation (5.3.3), we
conclude that there exists an absolute constant $c\leq 2000$ such that for
$\eta^{2}$ satisfying equation (26b), namely

$\displaystyle\eta^{2}$
$\displaystyle\leq\frac{1}{cL}\min\left\\{\frac{1}{K^{2}},\frac{1}{Kd^{\frac{1}{2}}},\frac{1}{K^{\frac{2}{3}}d^{\frac{1}{3}}\left(\frac{M^{2}}{L}\right)^{\frac{1}{3}}},\frac{1}{K\frac{M}{L^{\frac{1}{2}}}},\frac{1}{K^{\frac{2}{3}}d}\frac{L}{L_{\text{H}}^{\frac{2}{3}}},\frac{1}{K^{\frac{4}{3}}\frac{M}{L^{\frac{1}{2}}}}\left(\frac{L}{L_{\text{H}}^{\frac{2}{3}}}\right)^{\frac{1}{2}}\right\\},$

we have

$\displaystyle{\mathbb{P}}\left[-f(F(p_{0},q_{0}))+f(q_{0})-\frac{1}{2}\left\|E(p_{0},q_{0})\right\|_{2}^{2}+\frac{1}{2}\left\|p_{0}\right\|_{2}^{2}\geq-1/16\right]\geq
1-\frac{1}{16}.$

Plugging this bound in the inequality (47) yields that

$\displaystyle{\mathbb{E}}_{z\sim\mathcal{N}(0,\mathbb{I}_{d})}\left[\min\left\\{1,\frac{\exp(-\mathcal{H}(E(z,x),F(z,x)))}{\exp(-\mathcal{H}(z,x))}\right\\}\right]\geq
1-\frac{1}{8},$

which when plugged in equation (5.3.3) implies that
$d_{\text{\tiny{TV}}}\left(\mathcal{P}_{x},\mathcal{T}^{\text{before-
lazy}}_{x}\right)\leq 1/8$ for any $x\in\mathcal{R}_{s}$, as claimed. The
proof is now complete.

##  6 Discussion

In this paper, we derived non-asymptotic bounds on mixing time of Metropolized
Hamiltonian Monte Carlo for log-concave distributions. By choosing appropriate
step-size and number of leapfrog steps, we obtain mixing-time bounds for HMC
that are smaller than the best known mixing-time bounds for MALA. This
improvement can be seen as the benefit of using multi-step gradients in HMC.
An interesting open problem is to determine whether our HMC mixing-time bounds
are tight for log-concave sampling under the assumptions made in the paper.

Even though, we focused on the problem of sampling only from strongly and
weakly log-concave distribution, our Theorem 1 can be applied to general
distributions including nearly log-concave distributions as mentioned in
Appendix C.2. It would be interesting to determine the explicit expressions
for mixing-time of HMC for more general target distributions. The other main
contribution of our paper is to improve the warmness dependency in mixing
rates of Metropolized algorithms that are proved previously such as MRW and
MALA (Dwivedi et al., 2019). Our techniques are inspired by those used to
improve warmness dependency in the literature of discrete-state Markov chains.
It is an interesting future direction to determine if this warmness dependency
can be further improved to prove a convergence sub-linear in $d$ for HMC with
generic initialization even for small condition number $\kappa$.

### Acknowledgements

We would like to thank Wenlong Mou for his insights and helpful discussions
and Pierre Monmarché for pointing out an error in the previous revision. We
would also like to thank the anonymous reviewers for their helpful comments.
This work was supported by Office of Naval Research grant DOD
ONR-N00014-18-1-2640, and National Science Foundation Grant NSF-DMS-1612948 to
MJW and by ARO W911NF1710005, NSF-DMS-1613002, NSF-IIS-174134 and the Center
for Science of Information (CSoI), US NSF Science and Technology Center, under
grant agreement CCF-0939370 to BY.

1Appendix .tocmtappendix mtchapternone mtappendixsubsection
mtappendixsubsubsection

###### Contents

  1. 1 Introduction
    1. 1.1 Past work on HMC
    2. 1.2 Past work on mixing time dependency on initialization
    3. 1.3 Our contributions
  2. 2 Background and problem set-up
    1. 2.1 Monte Carlo Markov chain methods
    2. 2.2 Markov chain basics
    3. 2.3 From Metropolized random walk to HMC
      1. 2.3.1 MRW and MALA algorithms
      2. 2.3.2 HMC sampling
  3. 3 Main results
    1. 3.1 Assumptions on the target distribution
    2. 3.2 Mixing time bounds for HMC
      1. 3.2.1 Strongly log-concave target
    3. 3.3 Mixing time bounds via conductance profile
    4. 3.4 Improved warmness dependency for MALA and MRW
  4. 4 Numerical experiments
  5. 5 Proofs
    1. 5.1 Overlap bounds for HMC
    2. 5.2 Proof of Theorem 1
    3. 5.3 Proof of Lemma 6
      1. 5.3.1 Properties of the HMC proposal
      2. 5.3.2 Proof of claim (27a) in Lemma 6
      3. 5.3.3 Proof of claim (27b) in Lemma 6
  6. 6 Discussion
  7. A Proof of Lemmas 3, 4 and 6
    1. A.1 Proof of Lemma 3
      1. A.1.1 Proof of Lemma 11
      2. A.1.2 Proof of Lemma 12
    2. A.2 Proof of Lemma 4
    3. A.3 Proofs related to Lemma 6
      1. A.3.1 Proof of Lemma 7
      2. A.3.2 Proof of Lemma 8
      3. A.3.3 Proof of Lemma 9
      4. A.3.4 Proof of Lemma 10
  8. B Proof of Corollary 2
    1. B.1 Proof of Lemma 15
    2. B.2 Proof of Lemma 16
    3. B.3 Proof of Lemma 17
  9. C Beyond strongly log-concave target distributions
    1. C.1 Weakly log-concave target
    2. C.2 Non-log-concave target
  10. D Optimal choice for HMC hyper-parameters
    1. D.1 Optimal choices for Corollary 14
      1. D.1.1 Faster mixing time bounds
    2. D.2 Comparison with guarantees for unadjusted versions of HMC

##  A Proof of Lemmas 3, 4 and 6

In this appendix, we collect the proofs of Lemmas 3, and 4, as previously
stated in Section 3.3, that are used in proving Theorem 1. Moreover, we
provide the proof of auxiliary results related to HMC proposal that were used
in the proof of Lemma 6.

###  A.1 Proof of Lemma 3

In order to prove Lemma 3, we begin by adapting the spectral profile technique
(Goel et al., 2006) to the continuous state setting, and next we relate
conductance profile with the spectral profile.

First, we briefly recall the notation from Section 2.2. Let
$\Theta:\mathcal{X}\times\mathcal{B}(\mathcal{X})\rightarrow{}_{+}$ denote the
transition probability function for the Markov chain and let $\mathcal{T}$ be
the corresponding transition operator, which maps a probability measure to
another according to the transition probability $\Theta$. Note that for a
Markov chain satisfying the smooth chain assumption (16), if the distribution
$\mu$ admits a density then the distribution $\mathcal{T}(\mu)$ would also
admits a density. We use $\mathcal{T}_{x}$ as the shorthand for
$\mathcal{T}(\mathbf{\delta}_{x})$, the transition distribution of the Markov
chain at $x$.

Let $L_{2}(\pi^{*})$ be the space of square integrable functions under
function $\pi^{*}$. The _Dirichlet form_ $\mathcal{E}:L_{2}(\pi^{*})\times
L_{2}(\pi^{*})\rightarrow\real$ associated with the transition probability
$\Theta$ is given by

$\displaystyle\mathcal{E}(g,h)=\frac{1}{2}\int_{(x,y)\in\mathcal{X}^{2}}\left(g(x)-h(y)\right)^{2}\Theta(x,dy)\pi^{*}(x)dx.$
(53)

The expectation ${\mathbb{E}}_{\pi^{*}}:L_{2}(\pi^{*})\rightarrow\real$ and
the variance $\operatorname{Var}_{\pi^{*}}:L_{2}(\pi^{*})\rightarrow\real$
with respect to the density $\pi^{*}$ are given by

$\displaystyle{\mathbb{E}}_{\pi^{*}}(g)=\int_{x\in\mathcal{X}}g(x)\pi^{*}(x)dx\quad\text{and}\quad\operatorname{Var}_{\pi^{*}}(g)=\int_{x\in\mathcal{X}}\left(g(x)-{\mathbb{E}}_{\pi^{*}}(g)\right)^{2}\pi^{*}(x)dx.$
(54a)

Furthermore, for a pair of measurable sets $(S,\Omega)\subset\mathcal{X}^{2}$,
the $\Omega$-_restricted spectral gap for the set $S$_ is defined as

$\displaystyle\lambda_{\Omega}(S):\,=\inf_{g\in
c_{0,\Omega}^{+}(S)}\frac{\mathcal{E}(g,g)}{\operatorname{Var}_{\pi^{*}}(g\cdot\mathbf{1}_{\Omega})},$
(55)

where

$\displaystyle c_{0,\Omega}^{+}(S):\,=\left\\{g\in
L_{2}(\pi^{*})\mid\text{supp}(g)\subset S,g\geq
0,g\cdot\mathbf{1}_{\Omega}\neq\text{constant},{\mathbb{E}}[g^{2}\cdot\mathbf{1}_{\Omega^{c}}]\leq{\mathbb{E}}[g^{2}\cdot\mathbf{1}_{\Omega}]\right\\}.$
(56)

Finally, the $\Omega$-_restricted spectral profile_ $\Lambda_{\Omega}$ is
defined as

$\displaystyle\Lambda_{\Omega}(v):\,=\inf_{\Pi^{*}(S\cap\Omega)\in[0,v]}\lambda_{\Omega}(S),\quad\text{
for all }v\in\big{[}0,\infty).$ (57)

Note that we restrict the spectral profile to the set $\Omega$. Taking
$\Omega$ to be $\mathcal{X}$, our definition agrees with the standard
definition definitions of the restricted spectral gap and spectral profile in
the paper (Goel et al., 2006) for finite state space Markov chains to
continuous state space Markov chains.

We are now ready to state a mixing time bound using spectral profile.

###### Lemma 11

Consider a reversible irreducible $\zeta$-lazy Markov chain with stationary
distribution $\Pi^{*}$ satifying the smooth chain assumption (16). Given a
$\beta$-warm start $\mu_{0}$, an error tolerance $\epsilon\in(0,1)$ and a set
$\Omega\subset\mathcal{X}$ with $\Pi^{*}(\Omega)\geq
1-\frac{\epsilon^{2}}{3\beta^{2}}$, the $L_{2}$-mixing time is bounded as

$\displaystyle\tau_{2}(\epsilon;\mu_{0})\leq\left\lceil\int_{4/\beta}^{8/\epsilon^{2}}\frac{dv}{\zeta\cdot
v\Lambda_{\Omega}(v)}\right\rceil,$ (58)

where $\Lambda_{\Omega}$ denotes the $\Omega$-restricted spectral profile (57)
of the chain.

See Appendix A.1.1 for the proof.

In the next lemma, we state the relationship between the $\Omega$-restricted
spectral profile (57) of the Markov chain to its $\Omega$-restricted
conductance profile (14).

###### Lemma 12

For a Markov chain with state space $\mathcal{X}$ and stationary distribution
$\Pi^{*}$, given any measurable set $\Omega\subset\mathcal{X}$, its
$\Omega$-restricted spectral profile (57) and $\Omega$-restricted conductance
profile (14) are related as

$\displaystyle\Lambda_{\Omega}(v)\geq\begin{cases}\displaystyle\frac{\Phi_{\Omega}^{2}(v)}{4}\quad&\text{for
all }v\in\left[0,\frac{\Pi^{*}(\Omega)}{2}\right]\\\
\displaystyle\frac{\Phi_{\Omega}^{2}(\Pi^{*}(\Omega)/2)}{8}\quad&\text{for all
}v\in\big{(}\frac{\Pi^{*}(\Omega)}{2},\infty).\end{cases}$ (59)

See Appendix A.1.2 for the proof.

Lemma 3 now follows from Lemmas 11 and 12 as well as the definition (15) of
$\tilde{\Phi}_{\Omega}$.

####  A.1.1 Proof of Lemma 11

We need the following lemma, proved in for the case of finite state Markov
chains in Goel et al. (2006), which lower bounds the Dirichlet form in terms
of the spectral profile.

###### Lemma 13

For any measurable set $\Omega\subset\mathcal{X}$, any function
$g:\mathcal{X}\rightarrow{}_{+}$ such that $g\in L_{2}(\pi^{*})$,
$g\cdot\mathbf{1}_{\Omega}$ is not constant and
${\mathbb{E}}[g^{2}\cdot\mathbf{1}_{\Omega}]\geq
2{\mathbb{E}}[g^{2}\cdot\mathbf{1}_{\Omega^{c}}]$, we have

$\displaystyle\frac{\mathcal{E}(g,g)}{\operatorname{Var}_{\pi^{*}}(g\cdot\mathbf{1}_{\Omega})}\geq\frac{1}{2}\Lambda_{\Omega}\left(\frac{4\left({\mathbb{E}}_{\pi^{*}}(g\cdot\mathbf{1}_{\Omega})\right)^{2}}{\operatorname{Var}_{\pi^{*}}(g\cdot\mathbf{1}_{\Omega})}\right).$
(60)

The proof of Lemma 13 is a straightforward extension of Lemma 2.1 from Goel et
al. (2006), which deals with finite state spaces, to the continuous state
Markov chain. See the end of Section A.1.1 for the proof.

We are now equipped to prove Lemma 11.

##### Proof of Lemma 11:

We begin by introducing some notations. Recall that for any Markov chain
satisfying the smooth chain assumption (16), given an initial distribution
$\mu_{0}$ that admits a density, the distribution of the chain at any step $n$
also admits a density. As a result, we can define the ratio of the density of
the Markov chain at the $n$-th iteration
$h_{\mu_{0},n}:\mathcal{X}\rightarrow\real$ with respect to the target density
$\pi^{*}$ via the following recursion

$\displaystyle
h_{\mu_{0},0}(x)=\frac{\mu_{0}(x)}{\pi^{*}(x)}\quad\text{and}\quad
h_{\mu_{0},n+1}(x)=\frac{\mathcal{T}\left(\pi^{*}\cdot
h_{\mu_{0},n}\right)(x)}{\pi^{*}(x)},$

where we have used the notation $\mathcal{T}(\mu)(x)$ to denote the density of
the distribution $\mathcal{T}(\mu)$ at $x$. Note that

$\displaystyle{\mathbb{E}}_{\pi^{*}}(h_{\mu_{0},n})=1\quad\text{and}\quad{\mathbb{E}}_{\pi^{*}}(h_{\mu_{0},n}\cdot\mathbf{1}_{\Omega})\leq
1\quad\text{for all }n\geq 0,$ (61)

where $\Omega\subset\mathcal{X}$ is a measurable set.

We also define the quantity
$J(n):\,=\operatorname{Var}_{\pi^{*}}(h_{\mu_{0},n})$ (we prove the existence
of this variance below in Step (1)) and also
$\tilde{J}(n):\,=\operatorname{Var}_{\pi^{*}}(h_{\mu_{0},n}\cdot\mathbf{1}_{\Omega})$.
Note that the $L_{2}$-distance between the distribution of the chain at step
$n$ and the target distribution is given by

$\displaystyle
d_{2,\pi^{*}}(\mathcal{T}^{n}(\mu_{0}),\Pi^{*})=\left(\int_{x\in{}^{d}}\left(h_{\mu_{0},n}(x)-1\right)^{2}\pi^{*}(x)dx\right)^{1/2}=\operatorname{Var}_{\pi^{*}}(h_{\mu_{0},n}).$

Consequently, to prove the $\epsilon$-$L_{2}$ mixing time bound (58), it
suffices to show that for any measurable set $\Omega\subset\mathcal{X}$, with
$\Pi^{*}(\Omega)\geq 1-\frac{\epsilon^{2}}{3\beta^{2}}$, we have

$\displaystyle J(n)\leq\epsilon^{2}\quad\text{ for
}n\geq\left\lceil\int_{4/\beta}^{8/\epsilon^{2}}\frac{dv}{\zeta\cdot
v\Lambda_{\Omega}(v)}\right\rceil$ (62)

We now establish the claim (62) via a three-step argument: (1) we prove the
existence of the variance $J(n)$ for all $n\in\mathbb{N}$ and relate $J(n)$
with $\tilde{J}(n)$. (2) then we derive a recurrence relation for the
difference $J(n+1)-J(n)$ in terms of Dirichlet forms that shows the $J$ is a
decreasing function, and (3) finally, using an extension of the variance $J$
from natural indices to real numbers, we derive an explicit upper bound on the
number of steps taken by the chain until $J$ lies below the required
threshold.

##### Step (1):

Using the reversibility (2) of the chain, we find that

$\displaystyle
h_{\mu_{0},n+1}(x)dx=\frac{\int_{y\in\mathcal{X}}\Theta(y,dx)h_{\mu_{0},n}(y)\pi^{*}(y)dy}{\pi^{*}(x)}$
$\displaystyle=\frac{\int_{y\in\mathcal{X}}\Theta(x,dy)h_{\mu_{0},n}(y)\pi^{*}(x)dx}{\pi^{*}(x)}$
$\displaystyle=\int_{y\in\mathcal{X}}\Theta(x,dy)h_{\mu_{0},n}(y)dx$ (63)

Applying an induction argument along with the relationship (63) and the
initial condition $h_{\mu_{0},0}(x)\leq\beta$, we obtain that

$\displaystyle h_{\mu_{0},n}(x)\leq\beta,\quad\text{for all }n\geq 0.$ (64)

As a result, the variances of the functions $h_{\mu_{0},0}$ and
$h_{\mu_{0},n}\cdot\mathbf{1}_{\Omega}$ under the target density $\pi^{*}$ are
well-defined and

$\displaystyle J(n)=\int_{\mathcal{X}}h^{2}_{\mu_{0},n}(x)\pi^{*}(x)dx-1.$
(65)

Then we show that $J$ can be upper bounded via $\tilde{J}$ as follows

$\displaystyle J(n)-\tilde{J}(n)$
$\displaystyle=\operatorname{Var}_{\pi^{*}}(h_{\mu_{0},n})-\operatorname{Var}_{\pi^{*}}(h_{\mu_{0},n}\cdot\mathbf{1}_{\Omega})$
$\displaystyle=\int_{x\in\mathcal{X}\setminus\Omega}h^{2}_{\mu_{0},n}(x)\pi^{*}(x)dx-\left(\int_{x\in\mathcal{X}}h_{\mu_{0},n}(x)\pi^{*}(x)dx\right)^{2}$
$\displaystyle\quad\quad\quad+\left(\int_{x\in\Omega}h_{\mu_{0},n}(x)\pi^{*}(x)dx\right)^{2}$
$\displaystyle\leq\beta^{2}\left(1-\Pi^{*}(\Omega)\right)\leq\frac{\epsilon^{2}}{3}=:B,$
(66)

where the last inequality follows from the fact that $\Omega$ satisfies
$\Pi^{*}(\Omega)\geq 1-\epsilon^{2}/(3\beta^{2})$. Similarly, we have

$\displaystyle{\mathbb{E}}_{\pi^{*}}(h_{\mu_{0},n}^{2})-{\mathbb{E}}_{\pi^{*}}(h_{\mu_{0},n}^{2}\cdot\mathbf{1}_{\Omega})\leq
B$ (67)

##### Step (2):

First, note the following bound on $J(0)$:

$\displaystyle
J(0)=\int_{x\in\mathcal{X}}\frac{\mu_{0}(x)^{2}}{\pi^{*}(x)}dx-1\leq\beta\int_{x\in\mathcal{X}}\mu_{0}(x)dx-1\leq\beta-1.$
(68)

Define the two step transition kernel $\Theta\circ\Theta$ as

$\displaystyle\Theta\circ\Theta(y,dz)=\int_{x\in\mathcal{X}}\Theta(y,dx)\Theta(x,dz).$

We have

$\displaystyle J(n+1):\,=\operatorname{Var}_{\pi^{*}}(h_{\mu_{0},n+1})$
$\displaystyle=\int_{x\in\mathcal{X}}h^{2}_{\mu_{0},n+1}(x)\pi^{*}(x)dx-1$
$\displaystyle\stackrel{{\scriptstyle(i)}}{{=}}\int_{x\in\mathcal{X}}\int_{y\in\mathcal{X}}\Theta(y,dx)h_{\mu_{0},n}(y)\pi^{*}(y)dy\int_{z\in\mathcal{X}}\Theta(x,dz)h_{\mu_{0},n}(z)-1$
$\displaystyle=\int_{y,z\in\mathcal{X}^{2}}\Theta\circ\Theta(y,dz)h_{\mu_{0},n}(y)h_{\mu_{0},n}(z)\pi^{*}(y)dy-1,$

where step (i) follows from the relation (63). Using the above expression for
$J(n+1)$ and the expression from equation (65) for $J(n)$, we find that

$\displaystyle J(n+1)-J(n)$
$\displaystyle=\int_{\mathcal{X}^{2}}\Theta\circ\Theta(y,dz)h_{\mu_{0},n}(y)h_{\mu_{0},n}(z)\pi^{*}(y)dy-\int_{\mathcal{X}}h^{2}_{\mu_{0},n}(x)\pi^{*}(x)dx,$
$\displaystyle\stackrel{{\scriptstyle(a)}}{{=}}-\mathcal{E}_{\Theta\circ\Theta}(h_{\mu_{0},n},h_{\mu_{0},n}),$
(69)

where $\mathcal{E}_{\Theta\circ\Theta}$ is the Dirichlet form (53) with
transition probability $\Theta$ being replaced by $\Theta\circ\Theta$. We come
back to the proof of equality (a) at the end of this paragraph. Assuming it as
given at the moment, we proceed further. Since the Markov chain is
$\zeta$-lazy, we can relate the two Dirichlet forms
$\mathcal{E}_{\Theta\circ\Theta}$ and $\mathcal{E}_{\Theta}$ as follows: For
any $y,z\in\mathcal{X}$ such that $y\neq z$, we have

$\displaystyle\Theta\circ\Theta(y,dz)=\int_{x\in\mathcal{X}}\Theta(y,dx)\Theta(x,dz)$
$\displaystyle\geq\Theta(y,dy)\Theta(y,dz)+\Theta(y,dz)\Theta(z,dz)$
$\displaystyle\geq 2\zeta\Theta(y,dz).$ (70)

If
${\mathbb{E}}_{\pi^{*}}(h_{\mu_{0},n}^{2}\cdot\mathbf{1}_{\Omega})<2{\mathbb{E}}_{\pi^{*}}(h_{\mu_{0},n}^{2}\cdot\mathbf{1}_{\Omega^{c}})$,
then according to Equation (67), we have

$\displaystyle
J(n)\leq{\mathbb{E}}(h_{\mu_{0},n}^{2})\leq{\mathbb{E}}(h_{\mu_{0},n}^{2}\cdot\mathbf{1}_{\Omega})+B\leq
3B\leq\epsilon^{2},$ (71)

and we are done. If $h_{\mu_{0},n}^{2}\cdot\mathbf{1}_{\Omega}$ is constant,
then
$\operatorname{Var}_{\pi^{*}}(h_{\mu_{0},n}^{2}\cdot\mathbf{1}_{\Omega})=0$
and

$\displaystyle J(n)\leq\operatorname{Var}_{\pi^{*}}(h_{\mu_{0},n}^{2})\leq
0+B\leq\epsilon^{2}.$

Otherwise, we meet the assumptions of Lemma 13 and we have

$\displaystyle
J(n+1)-J(n)=-\mathcal{E}_{\Theta\circ\Theta}(h_{\mu_{0},n},h_{\mu_{0},n})$
$\displaystyle\stackrel{{\scriptstyle(i)}}{{\leq}}-2\zeta\mathcal{E}_{\Theta}(h_{\mu_{0},n},h_{\mu_{0},n})$
$\displaystyle\stackrel{{\scriptstyle(ii)}}{{\leq}}-\zeta\operatorname{Var}_{\pi^{*}}(h_{\mu_{0},n}\cdot\mathbf{1}_{\Omega})\Lambda_{\Omega}\left(\frac{4\left[{\mathbb{E}}_{\pi^{*}}(h_{\mu_{0},n}\cdot\mathbf{1}_{\Omega})\right]^{2}}{\operatorname{Var}_{\pi^{*}}(h_{\mu_{0},n}\cdot\mathbf{1}_{\Omega})}\right)$
$\displaystyle\stackrel{{\scriptstyle(iii)}}{{\leq}}-\zeta\cdot\left(J(n)-B\right)\Lambda_{\Omega}\left(\frac{4}{J(n)-B}\right).$
(72)

where step (i) follows from inequality (A.1.1), step (ii) follows from Lemma
13, and finally step (iii) follows from inequality (A.1.1) which implies that
$\operatorname{Var}_{\pi^{*}}(h_{\mu_{0},n}\cdot\mathbf{1}_{\Omega})\geq
J(n)-B$, and the fact that the spectral profile $\Lambda_{\Omega}$ is a non-
increasing function.

##### Proof of equality (a) in equation (69):

Since the distribution $\Pi^{*}$ is stationary with respect to the kernel
$\Theta$, it is also stationary with respect to the two step kernel
$\Theta\circ\Theta$. We now prove a more general claim: For any transition
kernel $K$ which has stationary distribution $\Pi^{*}$ and any measurable
function $h$, the Dirichlet form $\mathcal{E}_{K}$, defined by replacing
$\Theta$ with $K$ in equation (53), we have

$\displaystyle\mathcal{E}_{K}(h,h)=\int_{\mathcal{X}}h^{2}(x)\pi^{*}(x)dx-\int_{\mathcal{X}}\int_{\mathcal{X}}h(x)h(y)K(x,dy)\pi^{*}(x)dx.$
(73)

Note that invoking this claim with $K=\Theta\circ\Theta$ and $h=h_{\mu_{0},n}$
implies step (a) in equation (69). We now establish the claim (73). Expanding
the square in the definition (53), we obtain that

$\displaystyle\mathcal{E}_{K}(h,h)$
$\displaystyle=\frac{1}{2}\int_{\mathcal{X}}\int_{\mathcal{X}}h^{2}(x)K(x,dy)\pi^{*}(x)dx+\frac{1}{2}\int_{\mathcal{X}}\int_{\mathcal{X}}h^{2}(y)K(x,dy)\pi^{*}(x)dx$
$\displaystyle\quad\quad\quad\quad-\int_{\mathcal{X}}\int_{\mathcal{X}}h(x)h(y)K(x,dy)\pi^{*}(x)dx$
$\displaystyle\stackrel{{\scriptstyle(i)}}{{=}}\frac{1}{2}\int_{\mathcal{X}}h^{2}(x)\pi^{*}(x)dx+\frac{1}{2}\int_{\mathcal{X}}h^{2}(x)\pi^{*}(x)dx-\int_{\mathcal{X}}\int_{\mathcal{X}}h(x)h(y)K(x,dy)\pi^{*}(x)dx,$

where equality (i) follows from the following facts: For the first term, we
use the fact that $\int_{\mathcal{X}}K(x,dy)=1$ since $K$ is a transition
kernel, and, for the second term we use the fact that
$\int_{\mathcal{X}}K(x,dy)\pi^{*}(x)dx=\pi^{*}(y)dy$, since $\Pi^{*}$ is the
stationary distribution for the kernel $K$. The claim now follows.

##### Step (3):

Consider the domain extension of the function $J$ from $\mathbb{N}$ to the set
of non-negative real numbers + by piecewise linear interpolation. We abuse
notation and denote this extension also by $J$. The extended function $J$ is
continuous and is differentiable on the set ${}_{+}\backslash\mathbb{N}$. Let
$n^{*}\in{}_{+}\cup\left\\{\infty\right\\}$ denote the first index such that
$J(n^{*})\leq 3B$. Since $\Lambda_{\Omega}$ is non-increasing and $J$ is non-
increasing, we have

$\displaystyle
J^{\prime}(t)\leq-\zeta\cdot\left(J(t)-B\right)\Lambda_{\Omega}\left(\frac{4}{J(t)-B}\right)\quad\text{for
all }t\in{}_{+}\backslash\mathbb{N}\text{ such that }t\leq n^{*}.$ (74)

Moving the $J$ terms on one side and integrating for $t\leq n^{*}$, we obtain

$\displaystyle\int_{J(0)}^{J(t)}\frac{dJ}{\left(J-B\right)\cdot\Lambda_{\Omega}\left(\frac{4}{J-B}\right)}\leq-\zeta
t.$

Using the change of variable $v=4/\left(J-B\right)$, we obtain

$\displaystyle\zeta
t\leq\int_{4/\left(J(0)-B\right)}^{4/\left(J(t)-B\right)}\frac{dv}{v\Lambda_{\Omega}(v)}$
(75)

Furthermore, equation (75) implies that for
$T\geq\frac{1}{\zeta}\int_{4/\beta}^{8/\epsilon^{2}}\frac{dv}{v\Lambda_{\Omega}(v)}$,
we have

$\displaystyle\int_{4/\beta}^{8/\epsilon^{2}}\frac{dv}{v\Lambda_{\Omega}(v)}\leq\int_{4/\left(J(0)-B\right)}^{4/\left(J(T)-B\right)}\frac{dv}{v\Lambda_{\Omega}(v)}.$

The bound (68) and the fact that $B=\epsilon^{2}/3$ imply that
$4/(J(0)-B)>4/\beta$. Using this observation, the fact that
$0\leq\Lambda_{\Omega}(v)<\infty$ for $v\geq 4/\beta$ and combining with the
case in Equation (71), we conclude that $J$ satisfies

$\displaystyle J(T)\leq 3B=\epsilon^{2}\text{ or
}\frac{4}{J(T)-B}\geq\frac{8}{\epsilon^{2}}\quad\text{for
}T\geq\frac{1}{\zeta}\int_{4/\beta}^{8/\epsilon^{2}}\frac{dv}{v\Lambda(v)},$

which implies the claimed bound (62).

Finally, we turn to the proof of Lemma 13.

##### Proof of Lemma 13:

Fix a function $g:\mathcal{X}\rightarrow{}_{+}$ such that $g\in
L_{2}(\pi^{*})$ and $g\cdot\mathbf{1}_{\Omega}$ is not constant and
${\mathbb{E}}[g^{2}\cdot\mathbf{1}_{\Omega}]\geq
2{\mathbb{E}}[g^{2}\cdot\mathbf{1}_{\Omega^{c}}]$. Note that for any constant
$\gamma>0$, we have

$\displaystyle\mathcal{E}(g,g)$
$\displaystyle=\frac{1}{2}\int_{(x,y)\in\mathcal{X}^{2}}\left(g(x)-g(y)\right)^{2}\Theta(x,dy)\Pi^{*}(x)dx$
$\displaystyle=\frac{1}{2}\int_{(x,y)\in\mathcal{X}^{2}}\left((g(x)-\gamma)-(g(y)-\gamma)\right)^{2}\Theta(x,dy)\Pi^{*}(x)dx$
$\displaystyle=\mathcal{E}\left((g-\gamma),(g-\gamma)\right).$

Let
$\gamma=\operatorname{Var}_{\pi^{*}}(g\cdot\mathbf{1}_{\Omega})/4{\mathbb{E}}_{\pi^{*}}\left(g\cdot\mathbf{1}_{\Omega}\right)$.
We have

$\displaystyle\mathcal{E}(g,g)=\mathcal{E}\left((g-\gamma),(g-\gamma)\right)$
$\displaystyle\stackrel{{\scriptstyle(i)}}{{\geq}}\mathcal{E}\left((g-\gamma)_{+},(g-\gamma)_{+}\right)$
$\displaystyle\stackrel{{\scriptstyle(ii)}}{{\geq}}\operatorname{Var}_{\pi^{*}}\left((g-\gamma)_{+}\cdot\mathbf{1}_{\Omega}\right)\inf_{f\in
c_{0,\Omega}^{+}(\left\\{g>\gamma\right\\})}\frac{\mathcal{E}(f,f)}{\operatorname{Var}_{\pi^{*}}\left(f\cdot\mathbf{1}_{\Omega}\right)}$