# Condition-number-independent convergence rate of Riemannian Hamiltonian
Monte Carlo with numerical integrators

Yunbum Kook, Yin Tat Lee, Ruoqi Shen, Santosh S. Vempala Georgia Tech,
<EMAIL_ADDRESS>of Washington and Microsoft Research,
<EMAIL_ADDRESS>of Washington<EMAIL_ADDRESS>Tech,
<EMAIL_ADDRESS>
###### Abstract

We study the convergence rate of discretized Riemannian Hamiltonian Monte
Carlo on sampling from distributions in the form of $e^{-f(x)}$ on a convex
body $\mathcal{M}\subset\mathbb{R}^{n}$. We show that for distributions in the
form of $e^{-\alpha^{\top}x}$ on a polytope with $m$ constraints, the
convergence rate of a family of commonly-used integrators is independent of
$\left\|\alpha\right\|_{2}$ and the geometry of the polytope. In particular,
the implicit midpoint method (IMM) and the generalized Leapfrog method (LM)
have a mixing time of $\widetilde{O}\left(mn^{3}\right)$ to achieve $\epsilon$
total variation distance to the target distribution. These guarantees are
based on a general bound on the convergence rate for densities of the form
$e^{-f(x)}$ in terms of parameters of the manifold and the integrator. Our
theoretical guarantee complements the empirical results of [KLSV22], which
shows that RHMC with IMM can sample ill-conditioned, non-smooth and
constrained distributions in very high dimension efficiently in practice.

###### Contents

  1. 1 Introduction
    1. 1.1 Prior work
  2. 2 RHMC with numerical integrators
    1. 2.1 Basics of RHMC
    2. 2.2 Notation and setting
    3. 2.3 Discretized RHMC
  3. 3 Our results
  4. 4 Technical overview
    1. 4.1 Mixing time via $s$-conductance: isoperimetry and one-step coupling
    2. 4.2 One-step coupling of discretized RHMC
      1. 4.2.1 Coupling of ideal RHMC
      2. 4.2.2 Coupling between ideal and discretized RHMC & Rejection probability
    3. 4.3 Analysis of numerical integrators & Parameter estimation in polytopes
  5. 5 Convergence rate of ideal RHMC
    1. 5.1 Auxiliary function and smoothness parameters
      1. 5.1.1 Working in high probability region
      2. 5.1.2 Auxiliary function $\ell$ with parameters $\ell_{0}$ and $\ell_{1}$
      3. 5.1.3 Smoothness parameters $R_{1},R_{2},R_{3}$
    2. 5.2 One-step coupling and convergence rate
  6. 6 Convergence rate of discretized RHMC
    1. 6.1 Stability via self-concordance
    2. 6.2 Coupling between ideal and discretized RHMC
      1. 6.2.1 Dynamics of ideal and discretized RHMC
      2. 6.2.2 One-step coupling
    3. 6.3 Bound on rejection probability
  7. 7 Numerical integrators
    1. 7.1 Implicit midpoint method
      1. 7.1.1 Second-order
      2. 7.1.2 Sensitivity
    2. 7.2 Generalized Leapfrog method (Störmer–Verlet)
      1. 7.2.1 Second-order
      2. 7.2.2 Sensitivity
  8. 8 Convergence rate of RHMC in polytopes
    1. 8.1 Isoperimetry of convex set
    2. 8.2 Good region $\mathbb{\mathcal{M}}_{\rho}$
    3. 8.3 Auxiliary function $\ell$ and smoothness parameters $R$
    4. 8.4 Convergence rate of RHMC with numerical integrators
      1. 8.4.1 Implicit midpoint method
      2. 8.4.2 Generalized Leapfrog method
  9. A Definitions
  10. B Lemmas

##  1 Introduction

Efficient sampling from high dimensional distributions is a fundamental
question that arises in many fields such as statistics, machine learning, and
theoretical computer science. One class of distributions that arises in many
applications is constrained distributions, where the distribution is defined
on a constrained set. Sampling from such distribution can be an efficient way
to study the geometric properties of the constrained set when direct
computation is not feasible. For instance, in systems biology, a metabolic
network is defined by a set of equalities and inequalities that represents
feasible steady state reaction rates [LNP12, TSF+13]. For large metabolic
networks, sampling from the constraint set can be an efficient way to simulate
the biochemical network and evaluate its capacity. In mathematics, computing
the volume of the Birkhoff polytope plays a key role in several areas,
including algebraic geometry, and probability. However, computing the volume
exactly using algebraic representations can take years even for a small
dimension $n=11$. On the other hand, a sampling based-algorithm can compute
the volume efficiently up to dimension half-million [KLSV22].

##### Traditional samplers

The current primary approach for sampling is Markov Chain Monte Carlo (MCMC)
method, which for many problems is the only known method with provable
efficiency guarantees. For general non-smooth distributions, the traditional
class of samplers is the zeroth-order samplers, which query the density of the
distributions to determine the algorithm’s trajectory. This class of samplers
includes Ball walk [LS93, KLS97], its affine-invariant version Dikin walk
[KN12, LLV20] and Hit-and-Run [Smi84, Lov99], which avoids an explicit step
size. However, this class of sampler is inefficient in practice because it
intrinsically needs a step size smaller than $O(1/\sqrt{n})$, where $n$ is the
dimension, to avoid stepping outside the constraint set, which leads to a
bottleneck of quadratic mixing time in dimension. Moreover, without putting
the convex body into an isotropic position, which requires expensive
computation in practice, the mixing time of Ball walk and Hit-and-Run,
$\widetilde{O}(n^{2}R^{2})$, depends on the condition number $R$ of the convex
body. The condition number of the distributions appearing in practical
applications can be large, e.g., the condition number of RECON1 [KLD+16], a
human metabolic network, can be as large as $10^{6}$. Using Hit-and-Run to
sample from metabolic networks can be over 100 times slower than the better
algorithms on this problem [CV16]. Sampling from the Birkhoff polytope can be
prohibitively expensive for any dimension higher than $n=20$ [CV16].

Another class of samplers commonly used is the first-order samplers, which
update the Markov chain based on the gradient information. The mixing time of
the continuous processes as well as the various discretization methods of this
class of samplers has been studied in a long line of recent works. The most
well-studied first-order samplers include Langevin algorithm [Dal17, DCWY18,
DMM19, VW19, CLA+21, CEL+22], its variant Underdamped Langevin algorithm
[CCBJ18, SL19], and Hamiltonian Monte Carlo (HMC) [CV22, CDWY20, LST20]. The
mixing time of this class of samplers also suffers from dependence on the
condition number of the distributions. Moreover, this class of samplers cannot
be applied to constrained distributions directly because their Markov chain
can easily step outside the constraint set. Currently, popular sampling
packages such as Stan [Sta20] and Pyro [BCJ+19] that are based on this class
of samplers are not able to handle constrained distributions, despite their
effectiveness in other settings.

##### Non-Euclidean Samplers

Given the limitations of the traditional samplers, researchers have sought to
extend these methods to non-Euclidean samplers, which leverage the local
geometry of distributions to speed up the samplers. For instance, Riemannian
Hamiltonian Monte Carlo (RHMC) extends the traditional HMC by considering the
dynamics on a Riemannian manifold that uses a non-Euclidean metric
corresponding to the distribution’s local geometry. When combined with a local
metric induced by the Hessian of a self-concordant barrier function, RHMC can
sample from ill-conditioned and non-smooth distributions efficiently. A recent
work [KLSV22] showed that RHMC can achieve a 1000-fold acceleration on the
benchmark dataset RECON3D [KLD+16], the largest published human metabolic
network, compared to previous methods. While RHMC has demonstrated superior
practical performance, the convergence rate of discretized RHMC remains open.
[LV18] bounded the convergence rate of continuous RHMC in terms of the
isoperimetry and natural smoothness parameters of the associated Riemannian
manifold. However, to implement RHMC, sophisticated integrators such as
implicit midpoint integrator (IMM) or the generalized Leapfrog integrator (LM)
are necessary to maintain measure-preservation and time reversibility. Simple
integrators, such as the naive Leapfrog method, are not suitable for RHMC as
they are no longer symplectic on general Riemannian manifolds [CBMR19]. These
sophisticated integrators provide accurate discretization and efficient
convergence in practice, but their theoretical analysis is challenging. In
particular, there is no theoretical guarantee that the convergence rate of
RHMC remains independent of the condition number after discretization, which
is the main motivation for using non-Euclidean samplers in our case.

In fact, analyzing discretized non-Euclidean samplers has been a persistent
challenge in many recent works. Another commonly studied class of non-
Euclidean samplers is the Riemannian Langevin algorithm (RLA) [GC11], which
extends the Langevin algorithm to non-Euclidean space. A closely related
process is the Mirror Langevin diffusion (MLD) [ZPFP20], which is a special
case of RLA when the metric is given by the Hessian of a Legendre-type convex
potential $\phi$. Many recent works have focused on obtaining the convergence
rate of discretized MLD or RLA, but many of them require strong assumptions or
oracles for accurate discretization. The analysis of [ZPFP20, Jia21, LTVW22]
and the empirical results in [Jia21] suggest that unless a strong regularity
assumption between the target distribution and $\phi$ is satisfied, the naive
integrators can lead to a bias term that exists even when the step size tends
to zero. This bias arises from the third-order error terms resulting from non-
Euclidean geometries and is hard to control. [AC21] circumvented this issue by
proposing an alternative discretization method that uses the exact solution to
the Brownian motion term, but it remains unclear whether the discretization is
feasible for general $\phi$. Similarly, [GV22] analyzed the convergence rate
of RLA using an oracle to sample from the natural Brownian motion on the
manifold. Given the current limitations in our understanding of the
integrators for non-Euclidean samplers, we believe it is crucial to
investigate the integrators more thoroughly and explore alternative
integrators.

##### Contribution

We provide (to our knowledge) the first convergence rate of discretized RHMC
on a class of numerical integrators. We consider a general class of
constrained distributions that can be written as

$e^{-f(x)}\text{ subject to }x\in\mathcal{M},$ (1.1)

where we assume $f$ is a convex function and
$\mathcal{M}\subset\mathbb{R}^{n}$ is a convex body with a (highly) self-
concordant barrier. We give theoretical guarantees showing that a large class
of integrators can maintain smoothness and condition number independence when
sampling from distributions in the form of $e^{-\alpha^{\top}x}$ on a polytope
with $m$ constraints. In fact, many applications can be written in this form
because any log-concave density in the form of $\eqref{eq:constrained_dist}$
can be reduced to

$e^{-t}\text{ subject to }(x,t)\in\mathcal{M}^{\prime},$ (1.2)

where $\mathcal{M}^{\prime}=\left\\{(x,t):f(x)\leq t,x\in\mathcal{M}\right\\}$
is convex in $(x,t)$. We show for distributions in the form of
$e^{-\alpha^{\top}x}$, the implicit midpoint method (IMM) and the generalized
Leapfrog method (LM) have a mixing time of $\widetilde{O}\left(mn^{3}\right)$
to achieve $\epsilon$ total variation distance to the target distribution. In
addition, we give a general convergence result on sampling from distributions
in the form of $e^{-f(x)}$ on a convex body in terms of parameters of the
manifold and the integrator, which can be useful for future works that analyze
the convergence rate on other integrators or distributions.

While numerical integration is a rich and active field [HHIL06], and the study
of the local convergence of numerical estimators is quite sophisticated, we
are not aware of global polynomial-time mixing time guarantees based on
commonly-used numerical integrators such as IMM and LM. Our convergence result
is the theoretical foundation of [KLSV22] and extends [LV18] to settings of
practical importance. Our results apply to not only IMM and LM, but also a
more general class of symplectic and time-reversible integrators that
satisfies a sensitivity condition, which advances our understanding of
integrators for RHMC and the more general non-Euclidean samplers.

Moreover, in our algorithm, we use a Metropolis filter to correct the
distribution, which is a crucial step for high-accuracy sampling. To address
the discretization issues of RLA and MLD, applying a Metropolis filter to
correct the bias is one potential solution. Nevertheless, to the best of our
knowledge, there is no general-purpose metropolized non-Euclidean Langevin
algorithm in the literature. We believe that our analysis of metropolized RHMC
can provide valuable insights into the design and analysis of future
metropolized non-Euclidean Langevin algorithms.

It is important for readers to be aware that although the convergence rate we
obtain is independent of the condition number, the convergence rate is likely
to be far from optimal due to the complicated analysis of the integrators
used. To couple the discretized and ideal RHMC in our analysis, we need a step
size much smaller than what is typically required in practice. [KLSV22]
demonstrated that RHMC with IMM can achieve sublinear mixing times in
dimension on metabolic networks and structured polytopes including hypercubes,
simplices, and Birkhoff polytopes. We believe a tighter convergence bound is
possible with more advanced analysis.

###  1.1 Prior work

The convergence rate of MCMC methods in sampling from a convex body has been a
topic of active research for decades (see [LV22] for a more detailed
discussion). The mixing time of ball walk on isotropic log-concave density is
bounded by $\widetilde{O}(n^{2})$ from a warm start [KLS97], where a convex
body can be put into a near isotropic position in $\widetilde{O}(n^{3})$
membership queries [JLLV21]. Dikin walk uses the local geometry to improve the
mixing rate to $O(mn)$ on polytopes, where $m$ is the number of constraints.
Moreover, due to its affine invariance, there is no need to put the polytope
into an isotropic position. With an LS barrier [LS14], Dikin walk can achieve
a mixing rate of $\widetilde{O}(n^{2})$ for any polytope [LLV20]. Geodesic
walk utilizes non-Euclidean geometry by taking a random walk on a manifold.
Geodesic walk with an exact exponential map and a Metropolis filter can
converge to the uniform density in $O(mn^{3/4})$ steps [LV17]. Continuous RHMC
avoids the use of a Metropolis filter due to its measure preservation and time
reversibility, which further improves the mixing time to $O(mn^{2/3})$ [LV18]
on uniform density. Our paper extends the mixing time result to discretized
RHMC with feasible integrators on more general distributions. Note that even
an extension to distribution $e^{-\alpha^{\top}x}$ needs nontrivial work to
avoid dependence on quantities such as the domain diameter.

##  2 RHMC with numerical integrators

###  2.1 Basics of RHMC

Hamiltonian Monte Carlo (HMC) is one of the most widely used MCMC methods and
is the default sampler implementation in many sampling packages ([Sta20,
SWF16, BCJ+19, KLSV22]). HMC introduces an auxiliary velocity variable $v$ in
addition to the position $x$, defines a joint density on $(x,v)$, and
determines its trajectory according to the _Hamiltonian dynamics_. The
Hamiltonian dynamics is characterized by the _Hamiltonian equations_ , the
first-order differential equations of the _Hamiltonian_ $H$ with respect to
$x$ and $v$. The Hamiltonian has a natural interpretation as the total energy
of a particle consisting of the kinetic and potential energy at position $x$
with velocity $v$.

The dynamic can be naturally generalized to the setting of Riemannian manifold
with local metric $\\{g(x)\\}_{x\in\mathbb{\mathcal{M}}}$. A natural extension
of the Hamiltonian is given by

$H(x,v)=f(x)+\frac{1}{2}v^{\top}g(x)^{-1}v+\frac{1}{2}\log\det g(x),$

with $g(x)$ viewed as a positive-definite matrix. For later use, we split $H$
into two parts $H_{1}(x,v)=f(x)+\frac{1}{2}\log\det g(x)$ and
$H_{2}(x,v)=\frac{1}{2}v^{\top}g(x)^{-1}v$. A curve
$(x(t),v(t))\in\mathbb{\mathcal{M}}\times
T_{x}\mathbb{\mathcal{M}}\subset\mathbb{R}^{n}\times\mathbb{R}^{n}$ is called
the _Hamiltonian curve_ if it is the solution to the Hamiltonian equations:

$\displaystyle\frac{dx}{dt}$
$\displaystyle=\frac{\mathcal{\partial}H}{\mathcal{\partial}v}(x,v)=g(x)^{-1}v,$
$\displaystyle\frac{dv}{dt}$
$\displaystyle=-\frac{\mathcal{\partial}H}{\mathcal{\partial}x}(x,v)=-\left(\underbrace{\nabla
f(x)+\frac{1}{2}\mathrm{Tr}\left[g(x)^{-1}Dg(x)\right]}_{\frac{\mathcal{\partial}H_{1}}{\mathcal{\partial}x}}+\underbrace{\left(-\frac{1}{2}Dg(x)\left[\frac{dx}{dt},\frac{dx}{dt}\right]\right)}_{\frac{\mathcal{\partial}H_{2}}{\mathcal{\partial}x}}\right).$
(2.1)

When clear from context, the Hamiltonian curve refers to
$x(t)\in\mathbb{\mathcal{M}}$ only. The Hamiltonian curves $(x(t),v(t))$ have
several geometric properties. For a map $F_{t}:(x,v)\mapsto(x(t),v(t))$,

  1. 1.

Hamiltonian preservation: $\frac{d}{dt}H(x(t),v(t))=0.$

  2. 2.

Symplectic: $DF_{t}(x,v)^{\top}\cdot J\cdot DF_{t}(x,v)=J$ for any $t\geq 0$
and $J=\left[\begin{array}[]{cc}0&I_{n}\\\ -I_{n}&0\end{array}\right]$.

  3. 3.

Measure-preservation: $\det(DF_{t}(x,v)=1$ for any $t\geq 0$. Note that
measure-preservation immediately follows from symplecticity.

  4. 4.

Time-reversible: $F_{t}(x(t),-v(t))=(x,-v)$.

Just as the Hamiltonian dynamics can be extended to the Riemannian setting,
[GC11] extended HMC to a Riemannian version called Riemannian Hamiltonian
Monte Carlo (RHMC); see Algorithm 1 for its one-step description. In fact, HMC
can be recovered from RHMC using the Euclidean metric (i.e., $g(x)=I$).

Our goal is to sample from a probability density proportional to $e^{-f(x)}$
supported on a convex body. To this end, we use RHMC with the Hamiltonian
$H:\mathbb{\mathcal{M}}\times\mathbb{R}^{n}\subset\mathbb{R}^{n}\times\mathbb{R}^{n}\to\mathbb{R}$,
viewing the convex body as a Riemannian manifold $\mathbb{\mathcal{M}}$ with a
local metric $g$.

Input: Initial point $x$, step size $h$

// Step 1: Sample an initial velocity $v$

Sample $v\sim\mathcal{N}(0,g(x))$.

// Step 2: Solve the Hamiltonian equations

Solve the Hamiltonian equations (2.1) to obtain $(x(t),v(t))$.

// Step 3: Metropolis-filter (skipped for ideal RHMC)

Accept $x(h)$ with probability
$\min\left(1,\frac{e^{-H(x(h),v(h))}}{e^{-H(x,v)}}\right)$. Otherwise, stay at
$x$.

Algorithm 1 Riemannian Hamiltonian Monte Carlo

###  2.2 Notation and setting

We use $\left(\mathbb{\mathcal{M}},g\right)$ to denote a connected and compact
Riemannian manifold with a boundary and a metric $g$ on which a target
distribution is supported. For a function
$f:\mathbb{\mathcal{M}}\subset\mathbb{R}^{n}\to\mathbb{R}$, we denote a target
distribution by $\pi(x)$ whose density is proportional to $e^{-f(x)}$ (i.e.,
$\frac{d\pi}{dx}\sim e^{-f(x)}$). We use $T_{x}\mathbb{\mathcal{M}}$ to denote
the tangent space of $\mathbb{\mathcal{M}}$ at $x\in\mathbb{\mathcal{M}}$. We
denote by $\pi_{x}$ the projection map onto $x$-space (i.e.,
$\pi_{x}(x,v)\stackrel{{\scriptstyle\mathrm{{\scriptscriptstyle
def}}}}{{=}}x$) and by $i_{x}$ the inclusion map (i.e.,
$i_{x}(v)\stackrel{{\scriptstyle\mathrm{{\scriptscriptstyle
def}}}}{{=}}(x,v)$). We reserve $h$ for the step size of RHMC.

With both manifold $\mathbb{\mathcal{M}}$ and tangent space
$T_{x}\mathbb{\mathcal{M}}$ endowed with the Euclidean metric, we define a map
$F_{t}:\mathbb{\mathcal{M}}\times
T_{x}\mathbb{\mathcal{M}}\to\mathbb{\mathcal{M}}\times\bigcup_{z\in\mathbb{\mathcal{M}}}T_{z}\mathbb{\mathcal{M}}$
by $F_{t}(x,v)\stackrel{{\scriptstyle\mathrm{{\scriptscriptstyle
def}}}}{{=}}(x(t),v(t))$, where $(x(t),v(t))$ is the solution to the
Hamiltonian equations at time $t$ with an initial condition $(x,v)$. In
particular, we define
$T_{x,h}:T_{x}\mathbb{\mathcal{M}}\to\mathbb{\mathcal{M}}$ by
$T_{x,h}(v)\stackrel{{\scriptstyle\mathrm{{\scriptscriptstyle
def}}}}{{=}}(\pi_{x}\circ F_{h}\circ i_{x})(v)=x(h)$. When both
$\mathbb{\mathcal{M}}$ and $T_{x}\mathbb{\mathcal{M}}$ are endowed with the
local metric $g$, we instead use
$\mathrm{Ham}_{x,t}:T_{x}\mathbb{\mathcal{M}}\to\mathbb{\mathcal{M}}$ defined
by $\mathrm{Ham}_{x,t}(v)\stackrel{{\scriptstyle\mathrm{{\scriptscriptstyle
def}}}}{{=}}x(t)$.

When a numerical integrator with step size $h$ outputs
$(\bar{x}_{h},\bar{v}_{h})$ by solving the Hamiltonian equations with an
initial condition $(x,v)$, we denote
$\overline{F}_{h}(x,v)\stackrel{{\scriptstyle\mathrm{{\scriptscriptstyle
def}}}}{{=}}(\bar{x}_{h},\bar{v}_{h})$ for a function
$\overline{F}_{h}:\mathbb{\mathcal{M}}\times
T_{x}\mathbb{\mathcal{M}}\to\mathbb{\mathcal{M}}\times\bigcup_{z\in\mathbb{\mathcal{M}}}T_{z}\mathbb{\mathcal{M}}$,
where the domain and range are endowed with the Euclidean metric. We define
$\overline{T}_{x,h}:T_{x}\mathbb{\mathcal{M}}\to\mathbb{\mathcal{M}}$ (endowed
with the Euclidean metric) by
$\overline{T}_{x,h}(v)=(\pi_{x}\circ\overline{F}_{h}\circ
i_{x})(v)=\bar{x}_{h}$. We drop $h$ from $T_{x,h},\overline{F}_{h}$ and
$\overline{T}_{x,h}$ if the step size is clear from context.

We assume that the domain $\mathbb{\mathcal{M}}\subset\mathbb{R}^{n}$ with a
boundary is convex and has a (highly) self-concordant barrier
$\phi:\mathbb{\mathcal{M}}\subset\mathbb{R}^{n}\to\mathbb{R}$ (Definition 54),
and that the metric $g$ is induced by the Hessian of the barrier (i.e.,
$g(x)=\nabla^{2}\phi(x)$). We denote the local norm of a vector $v$ by
$\left\|v\right\|_{x}$ or $\left\|v\right\|_{g(x)}$, and the Riemannian
distance by $d_{\phi}$ (Definition 56). We use $a\lesssim b$ to indicate that
$a\leq cb$ for some universal constant $c>0$.

###  2.3 Discretized RHMC

We use _ideal_ RHMC to denote the algorithm when the Hamiltonian equations in
Step 2 is accurately solved without any error. However, we cannot expect such
an accurate ODE solver to always exist in reality, so numerical integrators
with solutions that approximate the accurate ODE solutions are necessary. We
use _discretized_ RHMC to denote the algorithm when Step 2 of RHMC is solved
by a numerical integrator and a Metropolis-filter is used to correct the
distribution.

We now define a condition of numerical integrators that plays an important
role in our convergence-rate analysis.

######  Definition 1.

For a numerical integrator $\overline{F}$ and
$(x,v)\in\mathbb{\mathcal{M}}\times
T_{x}\mathbb{\mathcal{M}}\subset\mathbb{R}^{n}\times\mathbb{R}^{n}$, we call
$\overline{F}$ _sensitive_ at $(x,v)$ if there exists step size $h_{0}(x,v)$
such that the numerical integrator with step size $h$ less than $h_{0}$
satisfies

$\frac{\left|D\overline{T}_{x,h}(v)\right|}{\left|DT_{x,h}(v^{\prime})\right|}\geq
0.998,$

where $v^{\prime}$ satisfies $T_{x,h}(v^{\prime})=\overline{T}_{x,h}(v)$ and
the Jacobian $DT$ is taken with respect to the velocity variable. In other
words, the solution of the numerical integrator changes almost as fast as the
ideal solution does. Unless specified otherwise, a sensitive integrator is
additionally assumed to be measure-preserving (i.e.,
$\det(D\overline{F}_{h}(x,v))=1$) and time-reversible (i.e.,
$\overline{F}_{h}(\overline{x}_{h},-\overline{v}_{h})=(x,-v)$).

As a time-reversible numerical integrator is even-order, second-orderness
automatically follows. That is, for sufficiently small step size $h>0$,
$d_{g}(\bar{x}_{h},x_{h})\leq C_{x}(x,v)h^{2}$ and
$\left\|\bar{v}_{h}-v_{h}\right\|_{g(x)^{-1}}\leq C_{v}(x,v)h^{2}$ for some
functions of $x$ and $v$, $C_{x}$ and $C_{v}$. In other words, the errors of
the numerical integrator $\overline{F}_{h}$ with respect to the exact ODE
solver $F_{h}$ grow at most quadratically in the step size $h$.

This family of numerical integrators turns out to cover many commonly used
integrators in practice. For example, the implicit midpoint integrator (IMM)
(Algorithm 2) and the generalized Leapfrog integrator (LM) (Algorithm 3)
satisfy symplecticity, time-reversibility, and sensitivity (as shown in
Section 7). Measure-preservation gives the simple formula of the acceptance
probability in Step 3 of Algorithm 1. Measure-preservation together with time-
reversibility plays an important role in showing that the discretized RHMC
converges to its stationary distribution with density proportional to
$e^{-f(x)}$ (see Theorem 8 in [KLSV22]).

##  3 Our results

We analyze the mixing time of RHMC discretized by numerical integrators
commonly used in practice, with the Hamiltonian set to be
$H(x,v)=f(x)+\frac{1}{2}v^{\top}g(x)^{-1}v+\frac{1}{2}\log\det g(x).$ Previous
analysis of RHMC was based on high accuracy numerical integrators, which are
not always achievable in practice [LV18], and the complexity bounds were
derived for uniform density on a polytope. We extend the setting to sampling
exponential densities with practically feasible integrators. In the next
theorem, we denote
$\mathbb{\mathcal{M}}_{\rho}:=\left\\{x\in\mathbb{\mathcal{M}}:\left\|\alpha\right\|_{g(x)^{-1}}^{2}\leq
10n^{2}\log^{2}\frac{1}{\rho}\right\\}$ for $\rho>0$.

######  Theorem 2.

Let $\pi$ be a target distribution on a polytope with $m$ constraints in
$\mathbb{R}^{n}$ such that $\frac{d\pi}{dx}\sim e^{-\alpha^{\top}x}$ for
$\alpha\in\mathbb{R}^{n}$. Let $\mathbb{\mathcal{M}}$ be the Hessian manifold
of the polytope induced by the logarithmic barrier of the polytope. Let
$\Lambda=\sup_{S\subset\mathbb{\mathcal{M}}}\frac{\pi_{0}(S)}{\pi(S)}$ be the
warmness of the initial distribution $\pi_{0}$. Let $\pi_{T}$ be the
distribution obtained after $T$ steps of RHMC discretized by a sensitive
integrator on $\mathbb{\mathcal{M}}$. For any $\varepsilon>0$, if for
$x\in\mathbb{\mathcal{M}}_{\frac{\varepsilon}{2\Lambda}}$ and
$v\in\mathbb{R}^{n}$ randomly drawn from $\mathcal{N}(0,g(x))$, we have that
with probability at least $0.99,$ step size $h\leq h_{0}(x,v)$,

$h\leq\frac{10^{-20}}{n^{7/12}\log^{1/2}\frac{\Lambda}{\varepsilon}},\
hC_{x}(x,v)\leq\frac{10^{-20}}{\sqrt{n}},\
h^{2}C_{x}(x,v)\leq\frac{10^{-10}}{n\log\frac{\Lambda}{\varepsilon}}\text{ and
}h^{2}C_{v}(x,v)\leq\frac{10^{-10}}{\sqrt{n\log\frac{\Lambda}{\varepsilon}}},$

then $d_{\text{TV}}(\pi_{T},\pi)\leq\varepsilon$ for
$T=O\left(mh^{-2}\log\frac{\Lambda}{\varepsilon}\right)$.

By setting $C_{x}=C_{v}=0$, we can obtain the following corollary for the
mixing time of the ideal RHMC in this setting.

######  Corollary 3.

Let $\pi$ be a target distribution on a polytope with $m$ constraints in
$\mathbb{R}^{n}$ such that $\frac{d\pi}{dx}\sim e^{-\alpha^{\top}x}$ for
$\alpha\in\mathbb{R}^{n}$. Let $\mathbb{\mathcal{M}}$ be the Hessian manifold
of the polytope induced by the logarithmic barrier of the polytope. Let
$\Lambda=\sup_{S\subset\mathbb{\mathcal{M}}}\frac{\pi_{0}(S)}{\pi(S)}$ be the
warmness of the initial distribution $\pi_{0}$. Let $\pi_{T}$ be the
distribution obtained after $T$ iterations of the ideal RHMC on
$\mathbb{\mathcal{M}}$. For any $\varepsilon>0$ and step size
$h=O\left(\frac{1}{n^{7/12}\log^{1/2}\frac{\Lambda}{\varepsilon}}\right)$,
there exists $T=O\left(mn^{7/6}\log^{2}\frac{\Lambda}{\varepsilon}\right)$
such that $d_{\text{TV}}(\pi_{T},\pi)\leq\varepsilon$.

After we compute the parameters $C_{x}$ and $C_{v}$ of IMM and LM (see Section
7 and 8), and identify the sufficient conditions on the step size for their
sensitivity, the following mixing times of RHMC discretized by IMM or LM
immediately follow.

######  Corollary 4.

Let $\pi$ be a target distribution on a polytope with $m$ constraints in
$\mathbb{R}^{n}$ such that $\frac{d\pi}{dx}\sim e^{-\alpha^{\top}x}$ for
$\alpha\in\mathbb{R}^{n}$. Let $\mathbb{\mathcal{M}}$ be the Hessian manifold
of the polytope induced by the logarithmic barrier of the polytope. Let
$\Lambda=\sup_{S\subset\mathbb{\mathcal{M}}}\frac{\pi_{0}(S)}{\pi(S)}$ be the
warmness of the initial distribution $\pi_{0}$. Let $\pi_{T}$ be the
distribution obtained after $T$ iterations of RHMC discretized by IMM on
$\mathbb{\mathcal{M}}$. For any $\varepsilon>0$ and step size
$h=O\left(\frac{1}{n^{3/2}\log\frac{\Lambda}{\varepsilon}}\right)$, there
exists $T=O\left(mn^{3}\log^{3}\frac{\Lambda}{\varepsilon}\right)$ such that
$d_{\text{TV}}(\pi_{T},\pi)\leq\varepsilon$.

######  Corollary 5.

Let $\pi$ be a target distribution on a polytope with $m$ constraints in
$\mathbb{R}^{n}$ such that $\frac{d\pi}{dx}\sim e^{-\alpha^{\top}x}$ for
$\alpha\in\mathbb{R}^{n}$. Let $\mathbb{\mathcal{M}}$ be the Hessian manifold
of the polytope induced by the logarithmic barrier of the polytope. Let
$\Lambda=\sup_{S\subset\mathbb{\mathcal{M}}}\frac{\pi_{0}(S)}{\pi(S)}$ be the
warmness of the initial distribution $\pi_{0}$. Let $\pi_{T}$ be the
distribution obtained after $T$ iterations of RHMC discretized by LM on
$\mathbb{\mathcal{M}}$. For any $\varepsilon>0$ and step size
$h=O\left(\frac{1}{n^{3/2}\log\frac{\Lambda}{\varepsilon}}\right)$, there
exists $T=O\left(mn^{3}\log^{3}\frac{\Lambda}{\varepsilon}\right)$ such that
$d_{\text{TV}}(\pi_{T},\pi)\leq\varepsilon$.

In fact, Theorem 2 comes from a general result on the mixing time of RHMC for
density $e^{-f}$ on a convex body $\mathbb{\mathcal{M}}\subset\mathbb{R}^{n}$.
We provide its informal version here and defer its full statement (Theorem 24)
to Section 6.

######  Theorem.

(Informal) Let $\pi$ be a target distribution on a convex set
$\mathbb{\mathcal{M}}\subset\mathbb{R}^{n}$ and
$\Lambda=\sup_{S\subset\mathbb{\mathcal{M}}}\frac{\pi_{0}(S)}{\pi(S)}$ be the
warmness of the initial distribution $\pi_{0}$. Let $\mathbb{\mathcal{M}}$ be
the Hessian manifold with its metric induced by the Hessian of a highly self-
concordant barrier and $\pi_{T}$ the distribution obtained after $T$ steps of
RHMC discretized by a sensitive integrator on $\mathbb{\mathcal{M}}$. For any
$\varepsilon>0$, let
$\mathbb{\mathcal{M}}_{\frac{\varepsilon}{2\Lambda}}\subset\mathbb{\mathcal{M}}$
be a convex subset of measure at least $1-\frac{\varepsilon}{2\Lambda}$. There
is a step size bound $h_{0}$, defined in terms of smoothness parameters of the
manifold and the integrator, so that for any step size $h\leq h_{0}$, there
exists
$T=O\left(\left(h\psi_{\mathbb{\mathcal{M}}_{\frac{\varepsilon}{2\Lambda}}}\right)^{-2}\log\frac{\Lambda}{\varepsilon}\right)$
where $\psi_{\mathbb{\mathcal{M}}_{\frac{\varepsilon}{2\Lambda}}}$ is the
isoperimetry of $\mathbb{\mathcal{M}}_{\frac{\varepsilon}{2\Lambda}}$, such
that $d_{\text{TV}}(\pi_{T},\pi)\leq\varepsilon$.

##  4 Technical overview

RHMC with
---
numerical integrators
for $e^{-f}$
(Theorem 24)
Coupling of ideal
---
(Lemma 14)
Coupling between
---
ideal and discretized
(Lemma 29)
Rejection probability
---
(Lemma 32)
IMM & LM
---
(Section 7)
Parameters in
---
polytope setting
(Section 8)
RHMC with
---
numerical integrators
for $e^{-\alpha^{\top}x}$ in polytope
(Theorem 2)
RHMC with
---
IMM & LM
for $e^{-\alpha^{\top}x}$ in polytope
(Corollary 4, 5)
Figure 4.1: Proof outline

In this section, we provide a summary of the key proof ingredients that gives
the convergence rate of RHMC with numerical integrators to samples from
density $e^{-f(x)}$ on a convex body; see Figure 4.1 for the roadmap. In
Section 4.1, we review a general technique using $s$-conductance for bounding
the mixing time of a Markov chain. In Section 4.2, we summarize a refined
analysis of the ideal RHMC (Section 5) and the technique to couple the ideal
and discretized RHMC (Section 6). Finally in Section 4.3, we describe the
high-level ideas of our analysis of the numerical integrators (Section 7), IMM
and LM, and how to get the results for sampling from $e^{-\alpha^{\top}x}$ on
the Hessian manifolds of polytopes (Section 8).

###  4.1 Mixing time via $s$-conductance: isoperimetry and one-step coupling

Consider a Markov chain with a state space $\mathbb{\mathcal{M}}$, a
transition distribution $\mathbb{\mathcal{T}}_{x}$ and stationary distribution
$\pi$. We consider a _lazy_ Markov chain to avoid a uniqueness issue of the
stationary distribution. At each step, the lazy version of the Markov chain
does nothing with probability $\frac{1}{2}$ (i.e., stays at where it is and
does not move). Note that this change for the purpose of proof worsens the
mixing time only by a factor of $2$.

We use a standard conductance-based argument in [Vem05] to bound the mixing
time, which consists of two main ingredients – _the isoperimetry_ and _the
total variation (TV) distance coupling of one-step distributions_(Definition
57) staring from two close points.

######  Definition 6 ($s$-conductance).

Consider a Markov chain with a state space $\mathbb{\mathcal{M}}$, a
transition distribution $\mathbb{\mathcal{T}}_{x}$ and stationary distribution
$\pi$. For any $s\in[0,1/2)$, the $s$-_conductance_ of the Markov chain is

$\Phi_{s}\stackrel{{\scriptstyle\mathrm{{\scriptscriptstyle
def}}}}{{=}}\inf_{\pi(S)\in(s,1-s)}\frac{\int_{S}\mathbb{\mathcal{T}}_{x}(S^{c})\pi(x)dx}{\min(\pi(S)-s,\pi(S^{c})-s)}.$

As shown by [LS93], a lower bound on the $s$-conductance of a Markov chain
leads to an upper bound on the mixing time of the Markov chain.

######  Lemma 7 ([LS93]).

Let $\pi_{t}$ be the distribution of the points obtained after $t$ steps of a
lazy reversible Markov chain with the stationary distribution $\pi$. Let
$\Lambda=\sup_{S\subset\mathbb{\mathcal{M}}}\frac{\pi_{0}(S)}{\pi(S)}$ be the
warmness of an initial distribution $\pi_{0}$. For
$H_{s}=\sup\left\\{\left|\pi_{0}(A)-\pi(A)\right|:A\subset\mathbb{\mathcal{M}},\,\pi(A)\leq
s\right\\}$ with $0<s\leq\frac{1}{2}$, it follows that

$d_{\text{TV}}(\pi_{t},\pi)\leq
H_{s}+\frac{H_{s}}{s}\left(1-\frac{\Phi_{s}^{2}}{2}\right)^{t}.$

We now define the isoperimetry of a subset of $\mathbb{\mathcal{M}}$.

######  Definition 8 (Isoperimetry).

Let $(\mathbb{\mathcal{M}},g)$ be a Riemannian manifold and
$\mathbb{\mathcal{M}}^{\prime}$ a measurable subset of $\mathbb{\mathcal{M}}$
with $\pi(\mathbb{\mathcal{M}}^{\prime})>\frac{1}{2}$. The _isoperimetry_
$\psi$ of the subset with stationary distribution $\pi$ is defined by

$\psi_{\mathbb{\mathcal{M}}^{\prime}}=\inf_{S\subset\mathbb{\mathcal{M}}^{\prime}}\frac{\lim_{\delta\to
0^{+}}\frac{1}{\delta}\int_{\\{x\in\mathbb{\mathcal{M}}^{\prime}\,:\,0<d_{g}(S,x)\leq\delta\\}}\pi(x)dx}{\min(\pi(S),\pi(\mathbb{\mathcal{M}}^{\prime}\backslash
S))}.$

The following illustrates how one-step coupling with the isoperimetry leads to
a lower bound on the $s$-conductance. It can be proved similarly as Lemma 13
in [LV18].

######  Proposition 9.

For a Riemannian manifold $(\mathbb{\mathcal{M}},g)$, let $\pi$ be the
stationary distribution of a reversible Markov chain on $\mathbb{\mathcal{M}}$
with a transition distribution $\mathbb{\mathcal{T}}_{x}$. Let
$\mathbb{\mathcal{M}}^{\prime}\subset\mathbb{\mathcal{M}}$ be a subset with
$\pi(\mathbb{\mathcal{M}}^{\prime})\geq 1-\rho$ for some $\rho<\frac{1}{2}$.
We assume the following one-step coupling: if
$d_{g}(x,x^{\prime})\leq\Delta\leq 1$ for
$x,x^{\prime}\in\mathbb{\mathcal{M}}^{\prime}$, then
$d_{\text{TV}}(\mathbb{\mathcal{T}}_{x},\mathbb{\mathcal{T}}_{x^{\prime}})\leq
0.9$. Then for any $\rho\leq s<\frac{1}{2}$, the $s$-conductance is bounded
below by

$\Phi_{s}\geq\Omega\left(\psi_{\mathbb{\mathcal{M}}^{\prime}}\Delta\right).$

$\delta_{x}$$\delta_{y}$$d_{\text{TV}}(\mathcal{P}_{x},\mathcal{P}_{y})$:
Between ideal RHMCs (Lemma
14)$d_{\text{TV}}(\mathcal{P}_{x},\overline{\mathcal{P}}_{x})$: Between ideal
and discretized (Lemma
29)$d_{\text{TV}}(\mathcal{P}_{y},\overline{\mathcal{P}}_{y})$ Figure 4.2: An
illustration of our approach to one-step coupling. The thick line indicates
the ideal RHMC, and the dashed line indicates the discretized RHMC.

###  4.2 One-step coupling of discretized RHMC

In light of Proposition 9, we can focus on coupling the one-step distributions
of the discretized RHMC starting from two close-by points. Let
$\mathcal{P}_{x}$ and $\overline{\mathcal{P}}_{x}$ be the one-step
distributions on $\mathbb{\mathcal{M}}$ of the ideal and the discretized RHMC
starting from $x$, respectively. We use $\overline{\mathcal{P}}_{x}^{\prime}$
to denote the discretized RHMC without the Metropolis filter. As illustrated
in Figure 4.2, for two close points $x$ and $y$, the triangle inequality leads
to

$\displaystyle
d_{\text{TV}}(\overline{\mathcal{P}}_{x},\overline{\mathcal{P}}_{y})$
$\displaystyle\leq
d_{\text{TV}}(\overline{\mathcal{P}}_{x},\mathcal{P}_{x})+d_{\text{TV}}(\mathcal{P}_{x},\mathcal{P}_{y})+d_{\text{TV}}(\mathcal{P}_{y},\overline{\mathcal{P}}_{y})$
$\displaystyle\leq\left(d_{\text{TV}}(\overline{\mathcal{P}}_{x}^{\prime},\mathcal{P}_{x})+d_{\text{TV}}(\mathcal{P}_{x},\mathcal{P}_{y})+d_{\text{TV}}(\mathcal{P}_{y},\overline{\mathcal{P}}_{y}^{\prime})\right)+\left(d_{\text{TV}}(\overline{\mathcal{P}}_{x}^{\prime},\overline{\mathcal{P}}_{x})+d_{\text{TV}}(\overline{\mathcal{P}}_{y}^{\prime},\overline{\mathcal{P}}_{y})\right).$

Hence, it suffices to bound
$d_{\text{TV}}(\mathcal{P}_{x},\mathcal{P}_{y}),d_{\text{TV}}(\overline{\mathcal{P}}_{x}^{\prime},\mathcal{P}_{x})$
and
$d_{\text{TV}}(\overline{\mathcal{P}}_{x}^{\prime},\overline{\mathcal{P}}_{x})$,
respectively. We bound in Section 5 the first term
$d_{\text{TV}}(\mathcal{P}_{x},\mathcal{P}_{y})$, the TV distance of one-step
distributions of the ideal RHMC. For the remaining terms, when numerical
integrators do not preserve the Hamiltonian, a Metropolis filter is necessary
to ensure that the discretized RHMC converges to a target distribution. Due to
the filter, we need to handle a point-mass distribution at $x$. We address
this by first bounding the second term
$d_{\text{TV}}(\overline{\mathcal{P}}_{x}^{\prime},\mathcal{P}_{x})$, the TV
distance between the ideal and discretized RHMC _without_ the Metropolis
filter in Section 6.2. We then separately bound the rejection probability
$d_{\text{TV}}(\overline{\mathcal{P}}_{x}^{\prime},\overline{\mathcal{P}}_{x})$
in Section 6.3.

####  4.2.1 Coupling of ideal RHMC

We summarize how to bound $d_{\text{TV}}(\mathcal{P}_{x}{},\mathcal{P}_{y})$
here (see Section 5 for the full version).

######  Lemma.

(Informal, Lemma 14) For most of $x$ and $y$, and step size $h$ small enough,
if $d_{\phi}(x,y)\leq\frac{1}{100}$, then
$d_{\text{TV}}(\mathcal{P}_{x}{},\mathcal{P}_{y})\leq
O\left(\frac{1}{h}\right)d_{\phi}(x,y)+\frac{1}{25}$.

##### Previous approach

[LV18] provided a general framework for computing the mixing rate of RHMC on a
manifold embedded in $\mathbb{R}^{n}$, in terms of the isoperimetry and
smoothness parameters depending on the manifold and step size. One of the
major proof ingredients is one-step coupling: for two close points $x$ and
$y$, the one-step distributions at $x$ and $y$ have large overlap.

They use the notion of a ‘regular’ Hamiltonian curve, which enables them to
handle this task in low level, where the regularity can be understood as
_average_ behavior of Hamiltonian curves with high probability and is
quantified by some auxiliary functions. As the starting point of a regular
Hamiltonian curve changes from $x$ to $y$ along a length-minimizing geodesic
$c(s)$ joining $x$ and $y$, they find a one-to-one correspondence between
regular Hamiltonian curves started at $x$ and $y$, and bound
$\left|\frac{d}{ds}d_{\text{TV}}(\mathcal{P}_{x},\mathcal{P}_{c(s)})\right|$
over $s$. They achieve this by quantifying the rate of changes of the
probability density (see (5.1)).

It is daunting to directly work with the exact density function, so they make
use of the following techniques: (1) Show that the determinant of Jacobian is
close to $h^{n}$ up to small step size by applying a matrix-ODE theory to the
second-order ODE of the Hamiltonian equation (see Lemma 59). It allows them to
work with an approximate but simpler density with the Jacobian replaced by
$h^{n}$ (see (5.3)). (2) Establish the one-to-one correspondence along
variations of Hamiltonian curves by the implicit function theorem; for a given
endpoint $z$, as the starting point of a Hamiltonian curve moves along $c(s)$,
there exists a unique initial velocity $v_{c(s)}$ at each point on $c(s)$ that
brings $c(s)$ to the endpoint $z$ in step size $h$ (i.e.,
$\mathrm{Ham}_{c(s),h}(v_{c(s)})=z$). At the same time, by using the matrix-
ODE theory again they show that the regularity of Hamiltonian curves does not
blow up along $c(s)$ and quantify how much the proper initial velocity
changes.

##### Refined analysis

[LV18] bounded
$\left|\frac{d}{ds}d_{\text{TV}}(\mathcal{P}_{x},\mathcal{P}_{c(s)})\right|$
in terms of smoothness parameters, sumpremum bounds on some quantities defined
over the regular Hamiltonian curves starting from _any_ point in
$\mathbb{\mathcal{M}}$. However, considering all starting points leads to a
weaker coupling in the end. In fact, this makes sampling from an exponential
density have dependence on then condition number, since one of the smoothness
parameters requires the supremum bound on $\left\|\alpha\right\|_{g(x)^{-1}}$
over $x\in\mathbb{\mathcal{M}}$, which can be as large as
$\left\|\alpha\right\|_{2}$ times the diameter of the convex body.

To achieve a condition-number independent mixing time, we work in a convex
subset $\mathbb{\mathcal{M}}_{\rho}$ (call a good region) instead of
$\mathbb{\mathcal{M}}$, which requires refinement of the framework by
generalizing the smoothness parameters (Section 5.1) and theorems in their
paper accordingly. It allows us to obtain a stronger coupling by only
considering Hamiltonian curves starting from $\mathbb{\mathcal{M}}_{\rho}$.
This region is the region $\mathbb{\mathcal{M}}^{\prime}$ in Proposition 9.

This simple change, however, yields technical difficulties in following how
[LV18] proceeds with the original parameters. Recall in the one-step coupling,
they consider a Hamiltonian variation along a geodesic joining two points, but
the geodesic might step out of the good region. To address this issue, we use
the straight line between the points instead of the geodesic, as the straight
line is contained in $\mathbb{\mathcal{M}}_{\rho}$ due to the convexity. We
elaborate on how the technical details of the previous approach can be
modified accordingly under the redefined parameters and new variation curve in
order to get valid one-step coupling on this smaller region in Section 5.2.

####  4.2.2 Coupling between ideal and discretized RHMC & Rejection
probability

We provide a summary of Section 6, where we prove the following lemma and
Theorem 24.

######  Lemma.

(Informal, Lemma 29 and Lemma 32) For most of $(x,v)$, if step size $h$ is
small enough and falls under a sensitivity regime at $(x,v)$ of a numerical
integrator, then
$d_{\text{TV}}(\overline{\mathcal{P}}_{x}^{\prime},\mathcal{P}_{x})<\frac{1}{10}$
and
$d_{\text{TV}}(\overline{\mathcal{P}}_{x}^{\prime},\overline{\mathcal{P}}_{x})<\frac{1}{10^{3}}$.

For the former (bound on
$d_{\text{TV}}(\overline{\mathcal{P}}_{x}^{\prime},\mathcal{P}_{x})$), we show
that the densities of the ideal and discretized RHMC are similar by relating
two velocities $v$ and $v^{*}$, where $\overline{T}_{x}(v)=T_{x}(v^{*})$. It
can be reduced to establishing a constant lower bound on
$\frac{p_{x}^{\ast}(v^{*})}{p_{x}^{\ast}(v)}\frac{\left|D\overline{T}_{x}(v)\right|}{\left|DT_{x}(v^{*})\right|}$
for the probability density $p_{x}^{\ast}$ of Gaussian $\mathcal{N}(0,g(x))$.

We first define numerical integrators’ analogues of the smoothness parameters.
Then, we elaborate the idea above in Section 6.2, where we study the dynamics
of the ideal and discretized RHMC. In particular, we show the existence of
$v^{*}$ for a given $v$ by the Banach fixed-point theorem and a one-to-one
correspondence between them, together with an upper bound on
$\left\|v-v^{*}\right\|_{g^{-1}}$. This upper bound allows us to bound the
ratio of $p_{x}^{\ast}(v^{*})/p_{x}^{\ast}(v)$. We note that these results
heavily rely on the stability of local norm (see Section 6.1), which follows
from that the local metric is given by the Hessian of self-concordant
barriers. The ratio of the Jacobian follows from the sensitivity of the
integrators.

For the latter (bound on
$d_{\text{TV}}(\overline{\mathcal{P}}_{x}^{\prime},\overline{\mathcal{P}}_{x})$),
we observe that the acceptance probability comes down to bounding the
difference of the Hamiltonian at ideal and numerical solutions. To bound this,
we heavily use the stability of local norm as well as the quantitative
relationships between the ideal and discretized RHMC established above.
Putting these pieces together, we can obtain the mixing-time bound of the
discretized RHMC in Theorem 24.

###  4.3 Analysis of numerical integrators & Parameter estimation in polytopes

To apply the framework established so far, we analyze in Section 7 two
practical numerical integrators, IMM (Section 7.1) and LM (Section 7.2), by
estimating second-order parameters $C_{x}$ and $C_{v}$ and quantifying their
sensitivity regimes, and then apply these estimations to sampling from
distribution $e^{-\alpha^{\top}x}$ on a polytope in Section 8.

######  Lemma.

(Informal, adapted to polytope setting) For most of $(x,v)$, both IMM and LM
have $C_{x}(x,v)=\widetilde{O}(n)$ and $C_{v}(x,v)=\widetilde{O}(n^{3/2})$ for
step size $h=\widetilde{O}(1/\sqrt{n})$. The sensitivity region is
$h=\widetilde{O}(1/n)$ for IMM and $h=\widetilde{O}(1/n^{3/2})$ for LM.

To analyze one-step process of each numerical integrator, we need to keep
track all the quantities explicitly to obtain the condition-number
independence. However, the implicit nature of both integrators lead to
_coupled_ equations for $x$ and $v$, making the analysis complicated. To
address this, we handle these coupled equations parallelly by moving back and
forth between local norms at different points, where we use the stability of
local norm due to self-concordance. We remark that our approach to analyze
each integrator depends on the specific implementation of the integrator, so
each integrator requires slightly different techniques in this task.

For the sensitivity, we apply implicit differentiation to the one-step
equation of each integrator, obtaining a matrix equation in the form of
$(I-hE)D\overline{F}_{h}=hC$ for matrices $E,C\in\mathbb{R}^{2n\times 2n}$. We
use matrix-perturbation theory to quantify a sufficient condition on the step
size $h$ that ensures the invertibility of $(I-hE)$, obtaining an equation of
the form $D\overline{F}_{h}=\sum_{i=0}^{\infty}(hC^{\prime})^{i}$ for a matrix
$C^{\prime}\in\mathbb{R}^{2n\times 2n}$. By extracting the upper-right
$n\times n$ block matrix, it follows that $D\overline{T}_{h}=I+E^{\prime}$ for
$E^{\prime}=\sum_{i=1}^{\infty}(hC^{*})^{i}$ with a matrix
$C^{*}\in\mathbb{R}^{n\times n}$. Using the self-concordance of the local
metric, we get upper bounds on matrix quantities of $C^{*}$ including the
trace, two-norm and Frobenius norm. With $E^{\prime}$ viewed as perturbation,
we apply matrix-perturbation theory again to estimate a lower bound on
$\left|D\overline{T}_{h}\right|$.

Lastly in Section 8, we show that
$\mathbb{\mathcal{M}}_{\rho}=\left\\{x\in\mathbb{\mathcal{M}}:\left\|\alpha\right\|^{2}_{g(x)^{-1}}\leq
10n^{2}\log^{2}\frac{1}{\rho}\right\\}$ is convex by checking the second-order
condition and that $\mathbb{\mathcal{M}}_{\rho}$ has large measure by using a
functional inequality. Then we compute all parameters discussed so far –
isoperimetry, smoothness parameters of the manifold and numerical integrator –
for the polytope setting, putting them together to obtain the results in
Section 3.

##  5 Convergence rate of ideal RHMC

[LV18] provided a general framework for computing the mixing rate of RHMC on a
manifold embedded in $\mathbb{R}^{n}$. They represent the mixing rate in terms
of the isoperimetry and smoothness parameters depending on the manifold and
step size. In particular, they explicitly compute those parameters and
isoperimetry for the uniform distribution on a polytope with $m$ constraints,
concluding that the mixing rate of RHMC on the Hessian manifold induced by the
logarithmic barrier of the polytope is $O\left(mn^{2/3}\right)$. Notably, this
mixing rate is independent of the condition number of the polytope.
Independence of the condition number is desirable in practice, since real-
world instances are highly skewed and thus make it challenging for sampling
algorithms to sample efficiently.

Going beyond uniform sampling, we would like to obtain the condition-number-
independence of RHMC for more densities. However, even an extension to an
exponential density needs care to avoid dependence on a condition number (such
as the diameter of the domain).

In this section, we refine this framework by working on a subset
$\mathbb{\mathcal{M}}_{\rho}$ instead of $\mathbb{\mathcal{M}}$ and extending
the smoothness parameters and theorems developed in their paper accordingly.
It enables us to couple the one-step distributions of the ideal RHMC starting
at two close points by bounding the TV distance in terms of the smoothness
parameters.

###  5.1 Auxiliary function and smoothness parameters

We redefine those smoothness parameters in [LV18] that depend on a subset
$\mathbb{\mathcal{M}}_{\rho}$ of manifold (internally parameterized by
$\rho>0$) and step size $h$, pointing out how ours differ from the original
ones. We then develop the theory for one-step coupling based on the new
parameters.

####  5.1.1 Working in high probability region

When defining smoothness parameters, [LV18] pays attention to “well-behaved”
Hamiltonian curves $\gamma$ starting at _any_ point in $\mathbb{\mathcal{M}}$,
where the well-behavedness may be viewed as the average behavior of
Hamiltonian curves with _high probability_ and is quantified by some auxiliary
function. Then the smoothness parameters are estimated by bounding some
quantities along the curves. To do so, they should give supremum bounds on
those parameters over all points in $\mathbb{\mathcal{M}}$, which lead to a
weaker mixing rate in the end.

For a refined analysis, we apply a high-probability idea once again to
starting points of curves this time. In other words, we consider well-behaved
Hamiltonian curves starting only from a good region that has high probability.
Then we couple the one-step distributions at two close-by points only in this
region. This region will serve as $\mathbb{\mathcal{M}}^{\prime}$ in
Proposition 9.

This simple change, however, turns out to yield technical difficulties in
following how [LV18] proceeds with the original parameters. In bounding the
overlap of the one-step distributions, they deal with Hamiltonian curves and
Hamiltonian variations, starting points of which are on a geodesic between two
points, but the geodesic might step out of the good region. Hence, it leads to
us considering a different path joining two points instead of the geodesic. We
choose the straight line between two points instead and carefully check if the
original approach to one-step coupling still goes through. In addition to
this, we have to redefine each of the smoothness parameters and modify most of
the statements proven in [LV18] accordingly, as we work in the region smaller
than the entire domain. We now formalize this approach.

######  Definition 10.

Let $\pi$ be a target distribution on $\mathbb{\mathcal{M}}$ such that
$\frac{d\pi}{dx}\sim e^{-f(x)}$. Given $\rho>0$, we call a measurable subset
$\mathbb{\mathcal{M}}_{\rho}$ of $\mathbb{\mathcal{M}}$ a _good region_ if it
is convex and has measure $\pi(\mathbb{\mathcal{M}}_{\rho})\geq 1-\rho$.

##### Good region for exponential density.

As mentioned earlier, the necessity of a refined analysis naturally arises in
attempts to obtain a condition-number-independent mixing rate of RHMC for
$f(x)=e^{-\alpha^{\top}x}$. One of parameters in [LV18] depends on
$\sup_{x\in\mathbb{\mathcal{M}}}\left\|\nabla
f(x)\right\|_{g(x)^{-1}}^{2}=\sup_{x\in\mathbb{\mathcal{M}}}\left\|\alpha\right\|_{g(x)^{-1}}^{2}$,
but this supremum bound can be worsened by scaling up $\alpha$, and even for
fixed $\alpha$ it can be as large as the diameter of $\mathbb{\mathcal{M}}$.
To address this issue, for given $\rho>0$ we work on a smaller convex region
that has probability at least $1-\rho$, in which the quantity only depends on
the dimension $n$ and $\rho$, and set it to be a good region. To be precise,
we will take
$\mathbb{\mathcal{M}}_{\rho}=\left\\{x\in\mathbb{\mathcal{M}}:\left\|\nabla
f(x)\right\|_{g(x)^{-1}}^{2}\leq 10n^{2}\log^{2}\frac{1}{\rho}\right\\}$ for
the exponential densities.

####  5.1.2 Auxiliary function $\ell$ with parameters $\ell_{0}$ and
$\ell_{1}$

Initial velocities of Hamiltonian trajectories drawn from
$\mathcal{N}(0,g(x)^{-1})$ can be large even though it rarely happens, as seen
in the standard concentration inequality for Gaussian distributions. Since
those worst-case trajectories lead to a weaker coupling, [LV18] focuses on
“well-behaved” Hamiltonian trajectories rather than all trajectories. They
formalize this idea by defining an auxiliary function $\ell$, which measures
how regular a Hamiltonian trajectory is, along with two parameters $\ell_{0}$
and $\ell_{1}$.

######  Definition 11.

An auxiliary function $\ell$ with parameters $\ell_{0}$ and $\ell_{1}$ is a
function that assigns a non-negative real value to any Hamiltonian curve with
step size $h$, such that

  * •

For any $x\in\mathbb{\mathcal{M}}_{\rho}$, we have

$\mathbf{P}_{\gamma}\left(\ell(\gamma)>\frac{1}{2}\ell_{0}\right)<\frac{1}{100}\min\left(1,\frac{\ell_{0}}{\ell_{1}h}\right),$

where $\gamma$ is a Hamiltonian trajectory starting at $x$ with an initial
random velocity drawn from $\mathcal{N}(0,g(x)^{-1})$.

  * •

For any variations $\gamma_{s}$ starting from $\mathbb{\mathcal{M}}_{\rho}$
with $\ell(\gamma_{s})\leq\ell_{0}$, we have

$\left|\frac{d}{ds}\ell(\gamma_{s})\right|\leq\ell_{1}\cdot\left(\left\|\frac{d}{ds}\gamma_{s}(0)\right\|_{\gamma_{s}(0)}+\delta\left\|D_{s}\gamma_{s}^{\prime}(0)\right\|_{\gamma_{s}(0)}\right),$

where the variations $\gamma_{s}$ satisfy the Hamiltonian equations, and
$D_{s}$ denotes the covariant derivative of the velocity field
$\gamma_{s}^{\prime}(0)$ along a curve of starting points of the variations.

In the original definitions, the parameter $\ell_{0}$ is defined over the
Hamiltonian curves starting from any $x\in\mathbb{\mathcal{M}}$, and the
parameter $\ell_{1}$ is defined over the variations starting from any
$x\in\mathbb{\mathcal{M}}$ with $\ell(\gamma_{s})\leq\ell_{0}$.

Intuitions behind these parameters can be understood in the following way. The
auxiliary function $\ell$ measures how regular Hamiltonian trajectories are,
and $\ell_{0}$ serves as a threshold that allows us to consider only
Hamiltonian curves with regularity below the threshold, while it is large
enough to capture most trajectories.

To see the role of $\ell_{1}$, we run through a high-level idea for one-step
coupling. For a given endpoint $z$, we consider the set of regular Hamiltonian
curves $\gamma_{x}$ stating at $x$ with
$\ell(\gamma_{x})\leq\frac{1}{2}\ell_{0}$, which takes into account most
trajectories due to the definition of $\ell_{0}$. Along the straight line
joining $x$ and $y$, we smoothly vary the starting point of the Hamiltonian
curve to obtain a Hamiltonian curve $\gamma_{y}$ starting at $y$ with the same
endpoint $z$ and then find a correspondence between $\gamma_{x}$ and
$\gamma_{y}$. In doing so, it is desirable to maintain the regularity of
Hamiltonian curves. In other words, the auxiliary function should not change
rapidly so that $\ell(\gamma_{y})$ is still bounded by $\ell_{0}$. We enforce
this situation via the parameter $\ell_{1}$ that bounds the rate of change of
the auxiliary function, $\frac{d}{ds}\ell(\gamma_{s})$, along the straight
line.

####  5.1.3 Smoothness parameters $R_{1},R_{2},R_{3}$

In relating the regular Hamiltonian curves $\gamma_{x}$ and $\gamma_{y}$, some
quantities naturally arise from the proof. We begin with the definition of
Riemannian curvature tensor and then define three important parameters that
govern those quantities.

######  Definition 12.

The Riemannian curvature tensor is a map $R:V(\mathbb{\mathcal{M}})\times
V(\mathbb{\mathcal{M}})\times V(\mathbb{\mathcal{M}})\to
V(\mathbb{\mathcal{M}})$ for $V(\mathbb{\mathcal{M}})$, the collection of
vector fields on $\mathbb{\mathcal{M}}$, defined by

$R(u,v)w=\nabla_{u}\nabla_{v}w-\nabla_{v}\nabla_{u}w-\nabla_{[u,v]}w\quad\text{for
}u,v,w\in V(\mathbb{\mathcal{M}}),$

where $\nabla$ is the Levi-Civita connection on $\mathbb{\mathcal{M}}$, and
$[u,v]\stackrel{{\scriptstyle\mathrm{{\scriptscriptstyle
def}}}}{{=}}\nabla_{u}v-\nabla_{v}u$ is the Lie bracket of the vector fields
$u$ and $v$.

######  Definition 13.

Given an auxiliary function $\ell$ with parameters $\ell_{0}$ and $\ell_{1}$
and the operator $\Phi(\gamma,t):V(\mathbb{\mathcal{M}})\to
V(\mathbb{\mathcal{M}})$ defined by
$\Phi(\gamma,t)u\stackrel{{\scriptstyle\mathrm{{\scriptscriptstyle
def}}}}{{=}}D_{u}\mu(\gamma(t))-R(u,\gamma^{\prime}(t))\gamma^{\prime}(t)$,

  * •

$R_{1}$ is a parameter such that for any $t\in[0,h]$ and any Hamiltonian
curves $\gamma$ starting from $\mathbb{\mathcal{M}}_{\rho}$ with step size $h$
and $\ell(\gamma)\leq\ell_{0}$

$\left\|\Phi(\gamma,t)\right\|_{F,\gamma(t)}\stackrel{{\scriptstyle\mathrm{{\scriptscriptstyle
def}}}}{{=}}\sqrt{\mathbb{E}_{v,w\sim\mathcal{N}(0,g(\gamma(t))^{-1})}\left\langle
v,\Phi(\gamma,t)w\right\rangle_{g(\gamma(t))}^{2}}\leq R_{1}.$

  * •

$R_{2}$ is a parameter such that for any $t\in[0,h]$, any Hamiltonian curves
$\gamma$ starting from $\mathbb{\mathcal{M}}_{\rho}$ with step size $h$ and
$\ell(\gamma)\leq\ell_{0}$, any curve $c(s)$ starting from $\gamma(t)$ and any
vector field $v(s)$ along the curve $c(s)$ with $v(0)=\gamma^{\prime}(t)$,

$\left|\frac{d}{ds}\mathrm{Tr}\Phi(v(s))\bigg{|}_{s=0}\right|\leq
R_{2}\cdot\left(\left\|\frac{dc}{ds}\bigg{|}_{s=0}\right\|_{\gamma(t)}+h\left\|D_{s}v(s)|_{s=0}\right\|_{\gamma(t)}\right),$

where $v(s)$ in $\Phi$ indicates a Hamiltonian curve at time $t$ starting with
an initial condition $(c(s),v(s))$.

  * •

$R_{3}$ is a parameter such that for any Hamiltonian curves $\gamma$ starting
from $\mathbb{\mathcal{M}}_{\rho}$ with step size $h$ and
$\ell(\gamma)\leq\ell_{0}$, if $\zeta(t)\in T_{\gamma(t)}\mathbb{\mathcal{M}}$
is the parallel transport of the vector $\gamma^{\prime}(0)$ along $\gamma$,
then

$\sup_{t\in[0,h]}\left\|\Phi(\gamma,t)\zeta(t)\right\|_{\gamma(t)}\leq R_{3}.$

###  5.2 One-step coupling and convergence rate

In this section, we bound the TV distance of two one-step distributions of the
ideal RHMC starting at two close points in terms of the redefined parameters
and step size. The following result is a slight tweak of Theorem 29 in [LV18].

######  Lemma 14.

For $x,y\in\mathbb{\mathcal{M}}_{\rho}$ and step size
$h\leq\min\left(\frac{1}{10^{5}R_{1}^{1/2}},\left(\frac{\ell_{0}}{10^{3}R_{1}^{2}\ell_{1}}\right)^{1/5}\right)$,
if
$d_{\phi}(x,y)\leq\frac{1}{100}\min\left(1,\frac{\ell_{0}}{\ell_{1}}\right)$,
then

$d_{\text{TV}}(\mathcal{P}_{x},\mathcal{P}_{y})\leq
O\left(\frac{1}{h}+h^{2}R_{2}+hR_{3}\right)d_{\phi}(x,y)+\frac{1}{25}.$

This provides the convergence rate of the ideal RHMC for a general density
$e^{-f}$, which is a slight generalization of Theorem 30 in [LV18].

######  Proposition 15.

Let $\pi_{T}$ be the distribution obtained after $T$ steps of a lazy ideal
RHMC with the stationary distribution $\pi$ satisfying $\frac{d\pi}{dx}\sim
e^{-f(x)}$. Let
$\Lambda=\sup_{S\subset\mathbb{\mathcal{M}}}\frac{\pi_{0}(S)}{\pi(S)}$ be the
warmness of an initial distribution $\pi_{0}$. For any $\varepsilon>0$, let
$\rho=\frac{\varepsilon}{2\Lambda}$ and $\mathbb{\mathcal{M}}_{\rho}$ a good
region. If step size $h$ satisfies

$h^{2}\leq\frac{1}{10^{10}R_{1}},\
h^{5}\leq\frac{\ell_{0}}{10^{3}R_{1}^{2}\ell_{1}},\
h^{3}R_{2}+h^{2}R_{3}\leq\frac{1}{10^{10}}\text{ and
}h\leq\frac{1}{10^{10}}\min\left(1,\frac{\ell_{0}}{\ell_{1}}\right),$

where the parameters are defined in Definition 11 and 13, then for the
isoperimetry $\psi_{\mathbb{\mathcal{M}}_{\rho}}$ of
$\mathbb{\mathcal{M}}_{\rho}$ there exists
$T=O\left(\left(h\psi_{\mathbb{\mathcal{M}}_{\rho}}\right)^{-2}\log\frac{1}{\rho}\right)$
such that $d_{\text{TV}}(\pi_{T},\pi)\leq\varepsilon$.

Toward this result, we walk through how each lemma and theorem should change
so that they can be put together well, along with the modified smoothness
parameters and auxiliary function. We start with the formula of the
probability density of the one-step distribution at $x$.

######  Lemma 16 ([LV18], Lemma 10).

The probability density of one-step distribution of RHMC at
$x\in\mathbb{\mathcal{M}}\subset\mathbb{R}^{n}$ is

$p_{x}(z)=\sum_{v_{x}:\mathrm{Ham}_{x,h}(v_{x})=z}\underbrace{\left|D\mathrm{Ham}_{x,h}(v_{x})\right|^{-1}\sqrt{\frac{\left|g(z)\right|}{(2\pi)^{n}}}\exp\left(-\frac{1}{2}\left\|v_{x}\right\|_{x}^{2}\right)}_{\stackrel{{\scriptstyle\mathrm{{\scriptscriptstyle
def}}}}{{=}}p_{x}^{0}(v_{x})}.$ (5.1)

Note that the velocity $v_{x}$ is normalized by $g(x)^{-1}$, since the domain
of $\mathrm{Ham}_{x,h}$ is endowed with the local metric $g(x)$. In the
Euclidean coordinate, the density can be rewritten as

$\displaystyle p_{x}(z)$
$\displaystyle=\sum_{v_{x}^{\prime}:T_{x,h}(v_{x}^{\prime})=z}\left|DT_{x,h}(v_{x}^{\prime})\right|^{-1}\frac{1}{\sqrt{(2\pi)^{n}\left|g(x)\right|}}\exp\left(-\frac{1}{2}\left\|v_{x}^{\prime}\right\|_{g(x)^{-1}}^{2}\right)$
$\displaystyle=\sum_{v_{x}^{\prime}:T_{x,h}(v_{x}^{\prime})=z}\left|DT_{x,h}(v_{x}^{\prime})\right|^{-1}p_{x}^{*}(v_{x}^{\prime}),$
(5.2)

where $p_{x}^{*}$ is the probability density of the Gaussian distribution
$\mathcal{N}(0,g(x))$. This relation follows from $v_{x}^{\prime}=g(x)v_{x}$
and
$\left|DT_{x,h}(v_{x}^{\prime})\right|=\frac{\left|D\mathrm{Ham}_{x,h}(v_{x})\right|}{\sqrt{\left|g(x)\right|\left|g(x^{\prime})\right|}}$
(see the proof of Proposition 30) for
$x^{\prime}=\mathrm{Ham}_{x,h}(v_{x})=T_{x,h}(v_{x}^{\prime})$.

We can derive (5.2) in the following way as well. Intuitively, the probability
of moving from $x$ to $z$ through one step of RHMC is the summation of the
probability of choosing a proper initial velocity that brings $x$ to $z$,
which is the probability density function of $\mathcal{N}(0,g(x))$ divided by
$\left|DT_{x,h}(v_{x}^{\prime})\right|$. This Jacobian term comes from the
change of variables used when moving back from the position space $z$ to the
velocity space $v_{x}^{\prime}$.

##### High-level idea.

We are ready to run through a high-level idea of the one-step coupling proof
in [LV18]. For two close-by points, it is plausible that the probability
densities at $x$ and $y$ are similar, and one should relate those two
densities to quantify how close they are, which in turn results in a bound on
the overlap of two one-step distributions. It is a Hamiltonian curve that
enables them to handle this task in low level. Then they find a one-to-one
correspondence between the set of regular Hamiltonian curves from $x$ and $y$.

As one varies a starting point of a regular Hamiltonian curve from $x$ to $y$
along a curve $c(s)$ joining $x$ and $y$, one should quantify how fast each
term in (5.1) changes. To this end, they first prove that the determinant of
Jacobian is close to $h^{n}$ and that $\mathrm{Ham}_{x,h}$ is locally
injective, which makes it possible to work with an approximate but simpler
density

$\tilde{p}_{x}(z)\stackrel{{\scriptstyle\mathrm{{\scriptscriptstyle
def}}}}{{=}}\sum_{v_{x}:\mathrm{Ham}_{x,h}(v_{x})=z}\sqrt{\frac{\left|g(z)\right|}{(2\pi
h^{2})^{n}}}\exp\left(-\frac{1}{2}\left\|v_{x}\right\|_{x}^{2}\right).$ (5.3)

Next, they prove the following on variations of Hamiltonian curves; for a
given endpoint $z$, as the starting point of a Hamiltonian curve moves along
$c(s)$, there exists a unique initial velocity $v_{c(s)}$ at each point on
$c(s)$ that brings $c(s)$ to the fixed endpoint $z$ in step size $h$ (i.e.,
$\mathrm{Ham}_{c(s),h}(v_{c(s)})=z$). At the same time, they prove that the
regularity of Hamiltonian curves, $\ell(\gamma_{c(s)})$, does not change too
rapidly along $c(s)$, quantifying how much the proper initial velocity changes
as well. These are enough to achieve one-step coupling in terms of only
$R_{1}$. For further improvement, a more accurate estimate (5.4) of the
determinant of Jacobian is used instead, leading to an improved bound via
$R_{2}$ and $R_{3}$.

Following this approach, we elaborate on how each of the proof ingredients can
be formalized under the redefined parameters. The first ingredient about the
local injectivity of $\mathrm{Ham}_{x,h}$ and an approximation of its Jacobian
follows from Lemma 22 in [LV18] by restricting starting points of Hamiltonian
curves to $\mathbb{\mathcal{M}}_{\rho}$ in the statement.

######  Lemma 17.

Let $\gamma(t)=\mathrm{Ham}_{x,t}(v_{x})$ be a Hamiltonian curve starting at
$x\in\mathbb{\mathcal{M}}_{\rho}$ with $\ell(\gamma)\leq\ell_{0}$ and step
size $h$ satisfying $h^{2}\leq 1/R_{1}$. Then $D\mathrm{Ham}_{x,h}$ is
invertible and

$\left|\log\left|\frac{1}{h}D\mathrm{Ham}_{x,h}(v_{x})\right|-\int_{0}^{h}\frac{t(h-t)}{h}\mathrm{Tr}\Phi(t)dt\right|\leq\frac{(h^{2}R_{1})^{2}}{10}.$

As a corollary, we obtain the following estimate on the Jacobian of the
Hamiltonian map.

######  Corollary 18.

Let $(x(t),v(t))$ be the Hamiltonian curve starting with
$(x,v)\in\mathbb{\mathcal{M}}_{\rho}\times T_{x}\mathbb{\mathcal{M}}$, where
$T_{x}\mathbb{\mathcal{M}}\subset\mathbb{R}^{n}$ is endowed with the local
metric $g(x)$. For step size $h$ with
$h^{2}\leq\frac{1}{10^{5}\sqrt{n}R_{1}}$, and $v\in T_{x}\mathbb{\mathcal{M}}$
with $\ell(\mathrm{Ham}_{x,t}(v))\leq\ell_{0}$, we have

$h^{n}e^{-\frac{1}{600}}\leq\left|D\mathrm{Ham}_{x,h}(v)\right|\leq
h^{n}e^{\frac{1}{600}}.$

Namely, $\left|\frac{1}{h^{n}}\left|D\mathrm{Ham}_{x,h}(v)\right|-1\right|\leq
0.002$

The next one is the local uniqueness and existence of Hamiltonian variations,
obtained by adjusting Lemma 23 in [LV18].

######  Lemma 19.

Let $\gamma(t)=\mathrm{Ham}_{x,t}(v_{x})$ be a Hamiltonian curve starting at
$x\in\mathbb{\mathcal{M}}_{\rho}$ with $\ell(\gamma)\leq\ell_{0}$ and step
size $h$ satisfying $h^{2}\leq 1/R_{1}$. Let $x=\gamma(0)$ and $z=\gamma(h)$
be its endpoints.

  * •

For a neighborhood $U$ of $x\in\mathbb{\mathcal{M}}_{\rho}$ and neighborhood
$V$ of $v_{x}$, there exists a unique smooth invertible vector field $v:U\to
V$ such that $v(x)=v_{x}$ and $z=\mathrm{Ham}_{x,h}(v(y))$ for any $y\in U$.

  * •

For $\eta\in T_{x}\mathbb{\mathcal{M}}$, we have that
$\left\|\nabla_{\eta}v(x)\right\|_{x}\leq\frac{5}{2h}\left\|\eta\right\|_{x}$
and
$\left\|\frac{1}{h}\eta+\nabla_{\eta}v(x)\right\|_{x}\leq\frac{3}{2}R_{1}h\left\|\eta\right\|_{x}$.

  * •

Let $\gamma_{s}(t)=\mathrm{Ham}_{c(s),h}(v(c(s))$ be a variation of $\gamma$
along a path $c(s)$ in $U$ with $c(0)=x$ and $c^{\prime}(0)=\eta$. For
$t\in[0,h]$, we have
$\left\|\frac{\mathcal{\partial}\gamma_{s}(t)}{\mathcal{\partial}s}\big{|}_{s=0}\right\|_{\gamma(t)}\leq
5\left\|\eta\right\|_{x}$ and
$\left\|D_{s}\gamma_{s}^{\prime}(t)\big{|}_{s=0}\right\|_{\gamma(t)}\leq\frac{10}{h}\left\|\eta\right\|_{x}$.

The first item reveals the local uniqueness and existence of proper initial
velocities at any starting point around $x$. The second item bounds how fast
the initial velocity at $x$ change in a given direction $\eta$. The last item
extends the second result to each point on $\gamma(t)$.

The corresponding result in [LV18] is made for all regular Hamiltonian curves
from $\mathbb{\mathcal{M}}$. Since its proof relies on Lemma 17 to apply the
implicit function theorem to $f(y,w)=\mathrm{Ham}_{y,h}(w)$ and also on the
definition of $R_{1}$ for the second result, the statement should be
restricted to Hamiltonian curves starting from $\mathbb{\mathcal{M}}_{\rho}$.

We can now prove that regular Hamiltonian curves starting at $x$ with an
endpoint $z$ can be smoothly varied along the straight line between $x$ and
$y$, with the regularity of variations almost preserved.

######  Lemma 20.

Let $\gamma(t)=\mathrm{Ham}_{x,t}(v_{x})$ be a Hamiltonian curve starting at
$x\in\mathbb{\mathcal{M}}_{\rho}$ with $\ell(\gamma)\leq\frac{1}{2}\ell_{0}$
and step size $h$ satisfying $h^{2}\leq 1/R_{1}$. Let $x=\gamma(0)$ and
$z=\gamma(h)$ be its endpoints. For $y\in\mathbb{\mathcal{M}}_{\rho}$ and
$\beta=\frac{y-x}{\left\|y-x\right\|_{x}}$, let $c(s)=s\beta+x$ be a straight
line joining $x$ and $y$ with $c(0)=x$ and $c(s^{\prime})=y$. Let
$s^{\prime}\leq\frac{1}{100}\min\left(1,\frac{\ell_{0}}{\ell_{1}}\right)$.

  * •

There exists a unique velocity field $v$ along $c$ such that
$z=\mathrm{Ham}_{c(s),h}(v(c(s)))$. Furthermore, this vector field is also
uniquely determined by $c(s)$ and $v(c(s))$ on $c(s)$.

  * •

$\ell(\mathrm{Ham}_{c(s),h}(v(c(s)))\leq\ell_{0}$ for all $s$.

Compared to the original result, we use a straight line instead of a unit-
speed geodesic between $x$ and $y$, as the geodesic might escape the good
region $\mathbb{\mathcal{M}}_{\rho}$. The first item implies that there exists
an initial velocity $v(c(s))$ at each point $c(s)$ such that we can reach the
fixed endpoint $z$ via the Hamiltonian trajectory with the initial condition
$(c(s),v(c(s))$. The second item indicates that the regularity of such
Hamiltonian trajectories is preserved up to constant along the straight line.

###### Proof.

The first result can be proven similarly as in [LV18]. For the second, we
denote by $\gamma_{s}$ the Hamiltonian trajectory starting at $c(s)$ with the
proper initial velocity $v(c(s))$. We note that
$s^{\prime}=\left\|y-x\right\|_{x}$ and $\frac{dc(s)}{ds}=\beta$. By self-
concordance of $g$, we have that
$\left\|\beta\right\|_{c(s)}\leq\left(1+\left\|x-c(s)\right\|_{x}\right)\left\|\beta\right\|_{x}\leq
1+s^{\prime}$. Thus by Lemma 19,

$\left\|D_{s}v(s)\right\|_{c(s)}\leq\frac{5}{2h}\left\|\beta\right\|_{c(s)}\leq\frac{5}{2h}\left(1+\left\|x-c(s)\right\|_{x}\right),$

and

$\displaystyle\ell(\gamma_{y})$
$\displaystyle\leq\ell(\gamma_{x})+\int_{0}^{s^{\prime}}\left|\frac{d}{ds}\ell(\gamma_{s})\right|ds\leq\frac{1}{2}\ell_{0}+\ell_{1}\int_{0}^{s^{\prime}}\left(\left\|\beta\right\|_{c(s)}+\frac{5}{2}\left\|\beta\right\|_{c(s)}\right)ds$
$\displaystyle\leq\frac{1}{2}\ell_{0}+\ell_{1}s^{\prime}\cdot\frac{7}{2}(1+s^{\prime})\leq\frac{1}{2}\ell_{0}+\ell_{1}s^{\prime}\cdot
4$ $\displaystyle\leq\ell_{0},$

where we used that $1+s^{\prime}\leq 1.01$ and
$s^{\prime}\leq\frac{1}{100}\frac{\ell_{0}}{\ell_{1}}$. ∎

The next two lemmas provide bounds on some quantities via $R_{2}$ and $R_{3}$,
which are modifications of Lemma 34 and 32 in [LV18].

######  Lemma 21.

Let $\gamma_{s}$ be a family of Hamiltonian curves joining $c(s)$ and $z$
defined in Lemma 20 with $\ell(\gamma_{s})\leq\ell_{0}$ and step size $h$
satisfying $h^{2}\leq 1/R_{1}$. Then,

$\left|\int_{0}^{h}\frac{t(h-t)}{h}\frac{d}{ds}\mathrm{Tr}\Phi\left(\gamma_{s}^{\prime}(t)\right)dt\right|\leq
O\left(h^{2}R_{2}\right).$

Recall that $\gamma_{s}$ given in Lemma 20 has a starting point in
$\mathbb{\mathcal{M}}_{\rho}$ with $\ell(\gamma_{s})\leq\ell_{0}$. Since its
original proof uses the definition of $R_{2}$ and Lemma 19, and they are
applicable to regular Hamiltonian curves starting from
$\mathbb{\mathcal{M}}_{\rho}$ with $\ell(\gamma_{s})\leq\ell_{0}$, the
original proof of this lemma still works with our new definitions of the
parameters.

######  Lemma 22.

Let $\gamma(t)=\mathrm{Ham}_{x,t}(v_{x})$ be a Hamiltonian curve starting at
$x\in\mathbb{\mathcal{M}}_{\rho}$ with $\ell(\gamma)\leq\ell_{0}$ and step
size $h$ satisfying $h^{2}\leq 1/R_{1}$. Then,

$\frac{h}{2}\left|\nabla_{\eta}\left\|v(x)\right\|_{x}^{2}\right|\leq\left|\left\langle
v_{x},\eta\right\rangle_{x}\right|+3h^{2}R_{3}\left\|\eta\right\|_{x}.$

We can follow its original proof by using Lemma 19 and the definition of
$R_{3}$, as the regular Hamiltonian curve considered starts at
$x\in\mathbb{\mathcal{M}}_{\rho}$. We are now ready to prove Lemma 14.

###### Proof of Lemma 14.

Let $c(s)$ be the straight line joining $x$ and $y$, contained in
$\mathbb{\mathcal{M}}_{\rho}$ due to the convexity of
$\mathbb{\mathcal{M}}_{\rho}$. We denote
$\tilde{\ell}\stackrel{{\scriptstyle\mathrm{{\scriptscriptstyle
def}}}}{{=}}\min\left(1,\frac{\ell_{0}}{\ell_{1}h}\right)$. For
$x\in\mathbb{\mathcal{M}}_{\rho}$, let $V_{x}$ be the set of velocities
$v_{x}$ such that $\ell(\mathrm{Ham}_{x,h}(v_{x}))\leq\frac{1}{2}\ell_{0}$.
Note that $\mathcal{P}_{x}^{*}(V_{x}^{c})\leq\frac{1}{100}\tilde{\ell}$ by the
definition of $\ell_{0}$, where $\mathcal{P}_{x}^{*}$ is the one-step
distribution over velocities (not position) at $x$. Since $c(s)$ is contained
in $\mathbb{\mathcal{M}}_{\rho}$ and $\gamma(t)=\mathrm{Ham}_{x,t}(v_{x})$ has
regularity at most $\frac{1}{2}\ell_{0}$, Lemma 20 guarantees the existence of
a family of Hamiltonian variations $\gamma_{s}(t)$ joining $c(s)$ and
$\gamma(h)$ with $\ell(\gamma_{s})\leq\ell_{0}$ for all
$s\in[0,\left\|y-x\right\|_{x}]$.

We define an approximate probability density $\tilde{p}_{c(s)}$ of $p_{c(s)}$,
where $p_{c(s)}$ is the probability density of $\mathcal{P}_{c(s)}$. Driven by
Lemma 17, for $z\in\mathbb{\mathcal{M}}$ we define

$\tilde{p}_{c(s)}(z)\stackrel{{\scriptstyle\mathrm{{\scriptscriptstyle
def}}}}{{=}}\sum_{v:\,\mathrm{Ham}_{c(s),h}(v)=z}\underbrace{\sqrt{\frac{\left|g\left(\mathrm{Ham}_{c(s),h}(v)\right)\right|}{(2\pi
h^{2})^{n}}}\exp\left(-\int_{0}^{h}\frac{t(h-t)}{h}\mathrm{Tr}\Phi(\gamma_{s},t)dt\right)\cdot\exp\left(-\frac{1}{2}\left\|v\right\|_{c(s)}^{2}\right)}_{\stackrel{{\scriptstyle\mathrm{{\scriptscriptstyle
def}}}}{{=}}\tilde{p}_{c(s)}^{0}(v)},$ (5.4)

which is obtained by using
$\exp\left(-\int_{0}^{h}\frac{t(h-t)}{h}\mathrm{Tr}\Phi(t)dt\right)$ in place
of $\left|D\mathrm{Ham}_{c(s),h}(v)\right|^{-1}$ in $p_{x}(z)$ (see (5.1)).

We now relate $\tilde{p}_{c(s)}$ to $p_{c(s)}$. Note that the ratio of the
summand of $\tilde{p}_{c(s)}(z)$ and $p_{c(s)}(z)$ is equal to
$\tilde{p}_{c(s)}^{0}(v)/p_{c(s)}^{0}(v)=\frac{\frac{1}{h^{n}}\exp\left(-\int_{0}^{h}\frac{t(h-t)}{h}\mathrm{Tr}\Phi(\gamma_{s},t)dt\right)}{\left|D\mathrm{Ham}_{c(s),h}(v)\right|^{-1}}$
(see (5.1)). Due to $\ell(\gamma_{s})\leq\ell_{0}$, we can apply Lemma 17 to
$\gamma_{s}(t)$, obtaining

$\exp\left(-\frac{(h^{2}R_{1})^{2}}{10}\right)\leq\frac{\frac{1}{h^{n}}\exp\left(-\int_{0}^{h}\frac{t(h-t)}{h}\mathrm{Tr}\Phi(\gamma_{s},t)dt\right)}{\left|D\mathrm{Ham}_{c(s),h}(v)\right|^{-1}}\leq\exp\left(\frac{(h^{2}R_{1})^{2}}{10}\right).$

Using the conditions on the step size $h$, we can show that for
$C\stackrel{{\scriptstyle\mathrm{{\scriptscriptstyle
def}}}}{{=}}1+\frac{1}{10^{3}}\tilde{\ell}$

$\displaystyle\exp\left(\frac{(h^{2}R_{1})^{2}}{10}\right)$ $\displaystyle\leq
1+2\frac{(h^{2}R_{1})^{2}}{10}\leq
1+2\min\left(\frac{1}{10^{10}},\frac{\ell_{0}}{10^{3}\ell_{1}}\right)\leq C.$

Thus, the ratio is bounded below by $C^{-1}$ and above by $C$, and it implies
that

$C^{-1}\cdot p_{c(s)}^{0}(v)\leq\tilde{p}_{c(s)}^{0}(v)\leq C\cdot
p_{c(s)}^{0}(v).$ (5.5)

By Lemma 20, for each $v_{x}\in V_{x}$ with $\mathrm{Ham}_{x,h}(v_{x})=z$
there is a one-to-one correspondence between $v_{x}$ and $v_{y}$, where
$v_{y}$ satisfies $\mathrm{Ham}_{y,h}(v_{y})=z$. For this $v_{y}$, (5.5) leads
to

$\displaystyle p_{x}^{0}(v_{x})-p_{y}^{0}(v_{y})$ $\displaystyle\leq
C\cdot\tilde{p}_{x}^{0}(v_{x})-C^{-1}\cdot\tilde{p}_{y}^{0}(v_{y})$
$\displaystyle=\left(C^{2}-1\right)C^{-1}\tilde{p}_{x}^{0}(v_{x})+C^{-1}\left(\tilde{p}_{x}^{0}(v_{x})-\tilde{p}_{y}^{0}(v_{y})\right)$
$\displaystyle\leq\left(C^{2}-1\right)p_{x}^{0}(v_{x})+C^{-1}\left(\tilde{p}_{x}^{0}(v_{x})-\tilde{p}_{y}^{0}(v_{y})\right).$
(5.6)

In a similar way, we can show that

$\left(C^{-2}-1\right)p_{x}^{0}(v_{x})+C\left(\tilde{p}_{x}^{0}(v_{x})-\tilde{p}_{y}^{0}(v_{y})\right)\leq
p_{x}^{0}(v_{x})-p_{y}^{0}(v_{y}).$ (5.7)

Using this,

$\displaystyle p_{x}(z)-p_{y}(z)$
$\displaystyle=\sum_{v_{x}:\mathrm{Ham}_{x,h}(v_{x})=z}p_{x}^{0}(v_{x})-\sum_{v_{y}:\mathrm{Ham}_{y,h}(v_{y})=z}p_{y}^{0}(v_{y})$
$\displaystyle\leq\sum_{v_{x}\notin
V_{x}:\mathrm{Ham}_{x,h}(v_{x})=z}p_{x}^{0}(v_{x})+\sum_{v_{x}\in
V_{x}:\mathrm{Ham}_{x,h}(v_{x})=z}\left(p_{x}^{0}(v_{x})-p_{y}^{0}(v_{y})\right),$
(5.8)

where in the inequality we only left $v_{y}$ such that
$\mathrm{Ham}_{y,h}(v_{y})=z$ and that $v_{y}$ is the counterpart of $v_{x}\in
V_{x}$ given by the one-to-one correspondence.

We now bound the TV distance between $\mathcal{P}_{x}$ and $\mathcal{P}_{y}$
as follow:

$\displaystyle d_{\text{TV}}(\mathcal{P}_{x},\mathcal{P}_{y})$
$\displaystyle=\frac{1}{2}\int\left|p_{x}(z)-p_{y}(z)\right|dz$
$\displaystyle\underset{\eqref{eq:ineq1}}{\leq}\int_{z}\sum_{v_{x}\notin
V_{x}:\mathrm{Ham}_{x,h}(v_{x})=z}p_{x}^{0}(v_{x})dz+\int_{z}\sum_{v_{x}\in
V_{x}:\mathrm{Ham}_{x,h}(v_{x})=z}\left|p_{x}^{0}(v_{x})-p_{y}^{0}(v_{y})\right|dz$
$\displaystyle\underset{\eqref{eq:ineq2},\,\eqref{eq:ineq3}}{\leq}\mathcal{P}_{x}^{*}(V_{x}^{c})+\left(C^{2}-1\right)\int_{z}\sum_{v_{x}\in
V_{x}:\mathrm{Ham}_{x,h}(v_{x})=z}p_{x}^{0}(v_{x})dz$
$\displaystyle\qquad\qquad+2\int_{z}\sum_{v_{x}\in
V_{x}:\mathrm{Ham}_{x,h}(v_{x})=z}\left|\tilde{p}_{x}^{0}(v_{x})-\tilde{p}_{y}^{0}(v_{y})\right|dz$
$\displaystyle\leq\frac{\tilde{\ell}}{100}+\frac{\tilde{\ell}}{100}\int_{V_{x}}p_{x}^{*}(v)dv+2\int_{z}\sum_{v_{x}\in
V_{x}:\mathrm{Ham}_{x,h}(v_{x})=z}\int_{s}\left|\frac{d}{ds}\tilde{p}_{c(s)}^{0}(v_{c(s)})\right|dsdz$
$\displaystyle\leq\frac{\tilde{\ell}}{50}+2\int_{s}\underbrace{\int_{z}\sum_{v_{x}\in
V_{x}:\mathrm{Ham}_{x,h}(v_{x})=z}\left|\frac{d}{ds}\tilde{p}_{c(s)}^{0}(v_{c(s)})\right|dz}_{\stackrel{{\scriptstyle\mathrm{{\scriptscriptstyle
def}}}}{{=}}F_{s}}ds,$ (5.9)

where we used that $\int_{V_{x}}p_{x}^{*}(v)dv\leq 1$ in the last inequality,
and $v_{c(s)}$ is the initial velocity at $c(s)$ corresponding to $v_{x}\in
V_{x}$ (via the one-to-one correspondence).

Let us bound $F_{s}$ in terms of the parameters. From direct computation

$\frac{d}{ds}\tilde{p}_{c(s)}^{0}(v_{c(s)})=\left(-\int_{0}^{h}\frac{t(h-t)}{h}\frac{d}{ds}\mathrm{Tr}\Phi(\gamma_{s}^{\prime}(t))dt-\frac{1}{2}\frac{d}{ds}\left\|v_{c(s)}\right\|_{c(s)}^{2}\right)\tilde{p}_{c(s)}^{0}(v_{c(s)}).$

Due to $\tilde{p}_{c(s)}^{0}(v_{c(s)})\leq 2p_{c(s)}^{0}(v_{c(s)})$, we have

$\left|\frac{d}{ds}\tilde{p}_{c(s)}^{0}(v_{c(s)})\right|\leq
2\left(\left|\int_{0}^{h}\frac{t(h-t)}{h}\frac{d}{ds}\mathrm{Tr}\Phi(\gamma_{s}^{\prime}(t))dt\right|+\frac{1}{2}\left|\frac{d}{ds}\left\|v_{c(s)}\right\|_{c(s)}^{2}\right|\right)p_{c(s)}^{0}(v_{c(s)}).$

As $\ell(\gamma_{s})\leq\ell_{0}$ due to Lemma 20, it follow from Lemma 21
that

$\displaystyle F_{s}$ $\displaystyle\leq 4\int_{z}\sum_{v_{x}\in
V_{x}:\mathrm{Ham}_{x,h}(v_{x})=z}\left(\left|\int_{0}^{h}\frac{t(h-t)}{h}\frac{d}{ds}\mathrm{Tr}\Phi(\gamma_{s}^{\prime}(t))dt\right|+\frac{1}{2}\left|\frac{d}{ds}\left\|v_{c(s)}\right\|_{c(s)}^{2}\right|\right)p_{c(s)}^{0}(v_{c(s)})dz$
$\displaystyle\leq
O\left(h^{2}R_{2}\right)\int_{V_{x}}p_{v}^{*}(v)dv+2\int_{z}\sum_{v_{x}\in
V_{x}:\mathrm{Ham}_{x,h}(v_{x})=z}\left|\frac{d}{ds}\left\|v_{c(s)}\right\|_{c(s)}^{2}\right|p_{c(s)}^{0}(v_{c(s)})dz$
$\displaystyle\leq
O\left(h^{2}R_{2}\right)+2\underbrace{\int_{\left\\{v\,:\,\ell\left(\mathrm{Ham}_{c(s),h}(v)\right)\leq\ell_{0}\right\\}}\left|\frac{d}{ds}\left\|v\right\|_{c(s)}^{2}\right|p_{c(s)}^{*}(v)dv}_{\stackrel{{\scriptstyle\mathrm{{\scriptscriptstyle
def}}}}{{=}}S},$

where we used that $\int_{V_{x}}p_{v}^{*}(v)dv\leq 1$ again for the first term
and that $\ell(\gamma_{s})\leq\ell_{0}$ as well as the change of variable with
$z=\mathrm{Ham}_{c(s),h}(v(c(s)))$ for the second term.

We now bound $S$ in terms of $R_{3}$. As
$\ell(\gamma_{s})=\ell(\mathrm{Ham}_{c(s),h}(v(c(s))))\leq\ell_{0}$, we use
Lemma 22 to show that

$\displaystyle S$
$\displaystyle=\mathbb{E}_{\ell(\gamma_{s})\leq\ell_{0}}\left|\frac{d}{ds}\left\|v\right\|_{c(s)}^{2}\right|$
$\displaystyle\leq\frac{2}{h}\mathbb{E}_{\ell(\gamma_{s})\leq\ell_{0}}\left|\left\langle
v,\frac{d}{ds}c(s)\right\rangle_{c(s)}\right|+6hR_{3}\mathbb{E}_{\ell(\gamma_{s})\leq\ell_{0}}\left\|\frac{d}{ds}c(s)\right\|_{c(s)}.$

We recall from the proof of Lemma 20 that
$\left\|\frac{d}{ds}c(s)\right\|_{c(s)}\leq 1.01$. In addition to this, as $v$
is a Gaussian vector with respect to the local metric, $\left|\left\langle
v,\frac{d}{ds}c(s)\right\rangle_{c(s)}\right|=O(1)$ with high probability,
which easily follows from the standard concentration inequality for the
Gaussian distributions. Therefore,

$S\leq O\left(\frac{1}{h}\right)+7hR_{3}.$

Substituting this back to the inequality for $F_{s}$, we have

$F_{s}\leq O\left(h^{2}R_{2}+\frac{1}{h}+hR_{3}\right).$

Putting this to (5.9), it follows that

$d_{\text{TV}}(\mathcal{P}_{x},\mathcal{P}_{y})\leq
O\left(h^{2}R_{2}+\frac{1}{h}+hR_{3}\right)\left\|x-y\right\|_{x}+\frac{\tilde{\ell}}{50}.$

Due to $\left\|x-y\right\|_{g(x)}\leq 2d_{\phi}(x,y)$ by Lemma 64 and
$\tilde{\ell}\leq\frac{1}{100}$, it follows that

$d_{\text{TV}}(\mathcal{P}_{x},\mathcal{P}_{y})\leq
O\left(h^{2}R_{2}+\frac{1}{h}+hR_{3}\right)d_{\phi}(x,y)+\frac{1}{5000}.$

∎

Using this one-step coupling, we can prove Proposition 15 on the mixing time
of the ideal RHMC for a general density $e^{-f}$.

###### Proof.

Due to the assumptions on the step size $h$, Lemma 14 implies that if
$d_{\phi}(x,y)\lesssim h$, then
$d_{\text{TV}}(\mathcal{P}_{x},\mathcal{P}_{y})\leq\frac{1}{1000}$. By
Proposition 9 with $\rho=s=\frac{\epsilon}{2\Lambda}$, we can obtain the
following lower bound on the $s$-conductance:

$\Phi_{s}=\Omega\left(h\psi_{\mathbb{\mathcal{M}}_{\rho}}\right).$

By Lemma 7, we have

$d_{\text{TV}}(\pi_{t},\pi)\leq
s\Lambda+\Lambda\left(1-\frac{\Phi_{s}^{2}}{2}\right)^{t}.$

Therefore, it suffices to choose
$T=O\left(\left(h\psi_{\mathbb{\mathcal{M}}_{\rho}}\right)^{-2}\log\frac{1}{\rho}\right)$
to ensure $d_{\text{TV}}(\pi_{T},\pi)\leq\varepsilon$. ∎

##  6 Convergence rate of discretized RHMC

We bound the remaining two terms,
$d_{\text{TV}}(\overline{\mathcal{P}}_{x}^{\prime},\mathcal{P}_{x})$ in
Section 6.2 and
$d_{\text{TV}}(\overline{\mathcal{P}}_{x},\overline{\mathcal{P}}_{x}^{\prime})$
in Section 6.3, obtaining a result on the one-step coupling of RHMC
discretized by a numerical integrator with parameters $C_{x}$ and $C_{v}$. To
analyze the convergence rate of the discretized RHMC, we define additional
parameters.

######  Definition 23.

Given an auxiliary function $\ell$, a good region
$\mathbb{\mathcal{M}}_{\rho}$ and step size $h$, we define new parameters
$M_{1},M_{1}^{*},M_{2},M_{2}^{*}$ and
$\bar{\ell}_{0},\bar{\ell}_{1},\bar{R}_{1}$.

  * •

$M_{1}$ is a parameter such that for any $t\in[0,h]$ and any Hamiltonian curve
$\gamma$ starting at $x\in\mathbb{\mathcal{M}}_{\rho}$ with step size $h$ and
$\ell(\gamma)\leq\ell_{0}$

$n\leq M_{1}\quad\text{and}\quad\left\|\nabla
f(\gamma(t))\right\|_{g(x)^{-1}}^{2}\leq M_{1}.$

  * •

$M_{2}$ is a parameter such that for any $t\in[0,h]$ and any two Hamiltonian
curves $\gamma_{1},\gamma_{2}$ starting at $x\in\mathbb{\mathcal{M}}_{\rho}$
with step size $h$ and $\ell(\gamma_{i})\leq\ell_{0}$ for $i=1,2$

$\frac{\left\|\nabla f(\gamma_{1}(t))-\nabla
f(\gamma_{2}(t))\right\|_{g(x)^{-1}}}{\left\|\gamma_{1}(t)-\gamma_{2}(t)\right\|_{x}}\leq
M_{2}.$

  * •

Let $\gamma$ be any Hamiltonian curve $\gamma$ starting from
$(x,v)\in\mathbb{\mathcal{M}}_{\rho}\times T_{x}\mathbb{\mathcal{M}}$ with
step size $h$ and $\ell(\gamma)\leq\ell_{0}$. Let $\bar{x}_{j}$’s be
intermediate points produced by a numerical integrator with step size $h$ and
an initial condition $(x,v)$. We define $M_{1}^{*}$ to be the smallest number
such that for any $t\in[0,h]$

$\frac{\left|f(\gamma(t))-f(\bar{x}_{j})\right|}{\left\|\gamma(t)-\bar{x}_{j}\right\|_{x}}\leq\sqrt{M_{1}^{*}}\
\text{for all }j.$

We define $M_{2}^{*}$ to be the smallest number such that for any $t\in[0,h]$

$\frac{\left\|\nabla f(\gamma(t))-\nabla
f(\bar{x}_{j})\right\|_{g(x)^{-1}}}{\left\|\gamma(t)-\bar{x}_{j}\right\|_{x}}\leq
M_{2}^{*}\ \text{for all }j.$

  * •

Let $\overline{\mathbb{\mathcal{M}}_{\rho}}$ be a convex subset of
$\mathbb{\mathcal{M}}$ that contains $\bar{x}_{h}$ and $\gamma(h)$. We call an
auxiliary function $\bar{\ell}$ symmetric if
$\bar{\ell}(\mathrm{Ham}_{x,h}(v))=\bar{\ell}(\mathrm{Ham}_{x^{\prime},h}(-v^{\prime}))$
for $F_{h}(x,v)=(x^{\prime},v^{\prime})$. For a symmetric auxiliary function
$\bar{\ell}$, the parameters $\bar{\ell}_{0},\bar{\ell}_{1}$ and $\bar{R}_{1}$
are defined as in Definition 11 and 13 with
$\overline{\mathbb{\mathcal{M}}_{\rho}}$ in place of
$\mathbb{\mathcal{M}}_{\rho}$.

Note that such $\overline{\mathbb{\mathcal{M}}_{\rho}}$ always exists, as
$\mathbb{\mathcal{M}}$ is convex. We are now ready to formalize the informal
statement on the convergence rate of RHMC with a sensitive integrator for a
density $e^{-f}$ on the Hessian manifold induced by the highly self-concordant
barrier of $\mathbb{\mathcal{M}}$.

######  Theorem 24.

Let $\pi$ be a target distribution on a convex set
$\mathbb{\mathcal{M}}\subset\mathbb{R}^{n}$ and
$\Lambda=\sup_{S\subset\mathbb{\mathcal{M}}}\frac{\pi_{0}(S)}{\pi(S)}$ be the
warmness of the initial distribution $\pi_{0}$. Let $\mathbb{\mathcal{M}}$ be
the Hessian manifold with its metric induced by the Hessian of a strongly
self-concordant barrier and $\pi_{T}$ the distribution obtained after $T$
steps of RHMC discretized by a numerical integrator on $\mathbb{\mathcal{M}}$.
For any $\varepsilon>0$, let $\rho=\frac{\varepsilon}{2\Lambda}$ and
$\mathbb{\mathcal{M}}_{\rho}$ any good region. If step size $h$ guarantees the
sensitivity of the integrator and

$\displaystyle h^{2}$
$\displaystyle\leq\frac{10^{-10}}{\max(R_{1},\bar{R}_{1})},\
h^{5}\leq\frac{\ell_{0}}{10^{3}R_{1}^{2}\ell_{1}},\ h^{3}R_{2}+h^{2}R_{3}\leq
1,\ h\leq\frac{1}{10^{10}}\min\left(1,\frac{\ell_{0}}{\ell_{1}}\right),\
h^{2}\leq\frac{10^{-10}}{n+\sqrt{M_{1}}+M_{2}},$ $\displaystyle h$
$\displaystyle C_{x}(x,v)\leq\frac{10^{-10}}{\sqrt{n}},\ h^{2}C_{x}(x,v)\leq
10^{-10}\min\left(1,\frac{\bar{\ell}_{0}}{\bar{\ell}_{1}},\frac{1}{n+\sqrt{M_{1}}+\sqrt{M_{1}^{*}}}\right),\
h^{2}C_{v}(x,v)\leq\frac{10^{-10}}{\sqrt{n+\sqrt{M_{1}}}}$

for $x\in\mathbb{\mathcal{M}}_{\rho}$ and $v\in
V_{\text{good}}^{x}=\left\\{v\in\mathbb{R}^{n}:\left\|v\right\|_{g^{-1}}\leq
128\sqrt{n},\
\bar{\ell}(\mathrm{Ham}_{x,t}(g(x)^{-1}v))\leq\frac{1}{2}\bar{\ell}_{0}\right\\}$
(see (6.1)), where the parameters are defined in Definition 11, 13 and 23,
then for the isoperimetry $\psi_{\mathbb{\mathcal{M}}_{\rho}}$ of
$\mathbb{\mathcal{M}}_{\rho}$ there exists
$T=O\left(\left(h\psi_{\mathbb{\mathcal{M}}_{\rho}}\right)^{-2}\log\frac{1}{\rho}\right)$
such that $d_{\text{TV}}(\pi_{T},\pi)\leq\varepsilon$.

###  6.1 Stability via self-concordance

We summarize computational lemmas used in coupling one-step distributions and
bounding rejection probability. Going forward, the self-concordance of $g$ is
repetitively used to relate local metrics $g$ at two close points (see Lemma
25). We recall that $(1-\left\|x-y\right\|_{g(x)})^{2}g(x)\preceq
g(y)\preceq\frac{1}{(1-\left\|x-y\right\|_{g(x)})^{2}}g(x)$ for the local
metric $g$ induced by the Hessian of a self-concordant barrier when
$\left\|x-y\right\|_{g(x)}<c<1$. It implies that the local norm of a vector
with respect to $g(x)$ is within a constant factor of the local norm with
respect to $g(y)$ (and vice versa). Namely, for a vector $v$ we have
$\left\|v\right\|_{g(x)}\leq O(1)\cdot\left\|v\right\|_{g(y)}$ and
$\left\|v\right\|_{g(y)}\leq O(1)\cdot\left\|v\right\|_{g(x)}$. It enables us
to move back and forth between the local metric $g(x)$ and $g(y)$ whenever $x$
and $y$ are sufficiently close in the local metric $g(x)$ or $g(y)$.

######  Lemma 25.

Let $g(x)=\nabla^{2}\phi(x)$ for some highly self-concordant barrier $\phi$.

  * •

$(1-\|y-x\|_{g(x)})^{2}g(x)\preceq
g(y)\preceq\frac{1}{(1-\|y-x\|_{g(x)})^{2}}g(x).$

  * •

$\|Dg(x)[v,v]\|_{g(x)^{-1}}\leq 2\|v\|_{g(x)}^{2}.$

  * •

$\|Dg(x)[v,v]-Dg(y)[v,v]\|_{g(x)^{-1}}\leq\frac{6}{(1-\|y-x\|_{g(x)})^{3}}\|v\|_{g(x)}^{2}\|y-x\|_{g(x)}.$

  * •

$\left\|Dg(x)[v,v]-Dg(x)[w,w]\right\|_{g(x)^{-1}}\leq
2\left\|v-w\right\|_{g(x)}\left\|v+w\right\|_{g(x)}.$

###### Proof.

The first fact follows from Theorem 4.1.6 in [Nes03]. The second fact follows
from Lemma 4.1.2 in [Nes03]. To be precise,

$\left\|Dg(x)[v,v]\right\|_{g(x)^{-1}}=\max_{\left\|v\right\|_{g(x)}=1}Dg(x)[v,v,u]\leq
2\left\|v\right\|_{g(x)}^{2}.$

The third fact is from the following calculation:

$\displaystyle\|Dg(y)[v,v]-Dg(x)[v,v]\|_{g(x)^{-1}}$ $\displaystyle\leq$
$\displaystyle\int_{0}^{1}\|D^{2}g(x+t(y-x))[v,v,y-x]\|_{g(x)^{-1}}dt$
$\displaystyle\leq$
$\displaystyle\int_{0}^{1}\frac{1}{1-t\|y-x\|_{g(x)}}\|D^{2}g(x+t(y-x))[v,v,y-x]\|_{g(x+t(y-x))^{-1}}dt$
$\displaystyle\leq$
$\displaystyle\int_{0}^{1}\frac{6}{1-t\|y-x\|_{g(x)}}\|v\|_{g(x+t(y-x))}^{2}\|y-x\|_{g(x+t(y-x))}dt$
$\displaystyle\leq$
$\displaystyle\int_{0}^{1}\frac{6}{(1-t\|y-x\|_{g(x)})^{4}}dt\cdot\|v\|_{g(x)}^{2}\|y-x\|_{g(x)}$
$\displaystyle\leq$
$\displaystyle\frac{6}{(1-\|y-x\|_{g(x)})^{3}}\|v\|_{g(x)}^{2}\|y-x\|_{g(x)}.$

where the third and fifth line above follow from the first fact, and the
fourth line follows from Proposition 9.1.1 in [NN94].

The fourth fact is from the following calculation:

$\displaystyle\left\|Dg(x)[v,v]-Dg(x)[w,w]\right\|_{g(x)^{-1}}$
$\displaystyle=$
$\displaystyle\max_{\left\|u\right\|_{g(x)}=1}Dg(x)[v,v,u]-Dg(x)[w,w,u]$
$\displaystyle=$
$\displaystyle\max_{\left\|u\right\|_{g(x)}=1}Dg(x)[v-w,v,u]+Dg(x)[w,v-w,u]$
$\displaystyle=$
$\displaystyle\max_{\left\|u\right\|_{g(x)}=1}Dg(x)[v-w,v,u]+Dg(x)[v-w,w,u]$
$\displaystyle=$
$\displaystyle\max_{\left\|u\right\|_{g(x)}=1}Dg(x)[v-w,v+w,u]$
$\displaystyle\leq$
$\displaystyle\max_{\left\|u\right\|_{g(x)}=1}2\left\|v-w\right\|_{g(x)}\left\|v+w\right\|_{g(x)}\left\|u\right\|_{g(x)}$
$\displaystyle\leq$ $\displaystyle
2\left\|v-w\right\|_{g(x)}\left\|v+w\right\|_{g(x)}.$

∎

######  Lemma 26.

For $x,x^{\prime}\in\mathbb{\mathcal{M}}$, let $g=g(x)$ and
$g^{\prime}=g(x^{\prime})$. Let
$\delta_{x}:=\left\|x-x^{\prime}\right\|_{g}<0.99$ and
$\delta_{v}:=\left\|v-v^{\prime}\right\|_{g^{-1}}$.

  1. 1.

$(1-O(\delta_{x}))g\preceq g^{\prime}\preceq(1+O(\delta_{x}))g.$

  2. 2.

$(1-O(\delta_{x}))g^{\prime}\preceq g\preceq(1+O(\delta_{x}))g^{\prime}.$

  3. 3.

$(1-O(\delta_{x}))g^{-1}\preceq g^{\prime-1}\preceq(1+O(\delta_{x}))g^{-1}.$

  4. 4.

$(1-O(\delta_{x}))g^{\prime-1}\preceq
g^{-1}\preceq(1+O(\delta_{x}))g^{\prime-1}.$

  5. 5.

$-O(\delta_{x})I\preceq I-g^{\frac{1}{2}}g^{\prime-1}g^{\frac{1}{2}}\preceq
O(\delta_{x})I.$

  6. 6.

$-O(\delta_{x})I\preceq
I-g^{\prime\frac{1}{2}}g^{-1}g^{\prime\frac{1}{2}}\preceq O(\delta_{x})I.$

  7. 7.

$\left\|g^{\prime-\frac{1}{2}}g^{\frac{1}{2}}\right\|_{2}\leq
1+O(\delta_{x})\quad\&\quad\left\|g^{\prime\frac{1}{2}}g^{-\frac{1}{2}}\right\|_{2}\leq
1+O(\delta_{x}).$

  8. 8.

$\left\|g^{\frac{1}{2}}g^{\prime-\frac{1}{2}}\right\|_{2}\leq
1+O(\delta_{x})\quad\&\quad\left\|g^{-\frac{1}{2}}g^{\prime\frac{1}{2}}\right\|_{2}\leq
1+O(\delta_{x}).$

  9. 9.

$\left\|(g^{-1}-g^{\prime-1})p\right\|_{g}\lesssim\delta_{x}\left\|p\right\|_{g^{-1}}.$

  10. 10.

$\left\|g^{-1}p-g^{\prime-1}q\right\|_{g}\leq\left\|p-q\right\|_{g^{-1}}+O(\delta_{x})\left\|q\right\|_{g^{-1}}.$

  11. 11.

$\left\|\frac{\mathcal{\partial}H}{\mathcal{\partial}v}(x,v)-\frac{\mathcal{\partial}H}{\mathcal{\partial}v}(x^{\prime},v^{\prime})\right\|_{g}\leq\delta_{v}+O(\delta_{x})\left\|v^{\prime}\right\|_{g^{-1}}.$

  12. 12.

$\left\|\frac{\mathcal{\partial}H}{\mathcal{\partial}x}(x,v)-\frac{\mathcal{\partial}H}{\mathcal{\partial}x}(x^{\prime},v^{\prime})\right\|_{g^{-1}}\lesssim(\delta_{v}+\delta_{x}\left\|v\right\|_{g^{-1}})(\left\|v\right\|_{g^{-1}}+\left\|v^{\prime}\right\|_{g^{-1}})+n\delta_{x}+\left\|\nabla
f(x)-\nabla f(x^{\prime})\right\|_{g^{-1}}.$

###### Proof.

The first four lemmas follow from Lemma 25-1. For 5 (and 6, 7, 8 similarly),
using 3

$(1-O(\delta_{x}))I\preceq
g^{\frac{1}{2}}g^{\prime-1}g^{\frac{1}{2}}\preceq(1+O(\delta_{x}))I.$

Thus $-O(\delta_{x})I\preceq
I-g^{\frac{1}{2}}g^{\prime-1}g^{\frac{1}{2}}\preceq O(\delta_{x})I$. Also by
the definition of two-norm, it follows that

$\left\|g^{\prime-\frac{1}{2}}g^{\frac{1}{2}}\right\|_{2}\leq
1+O(\delta_{x}).$

Fact 9 follows from the following computation:

$\displaystyle\left\|(g^{-1}-g^{\prime-1})p\right\|_{g}$
$\displaystyle=\left\|(I-g^{\frac{1}{2}}g^{\prime-1}g^{\frac{1}{2}})g^{-\frac{1}{2}}p\right\|_{2}\leq
O(\delta_{x})\left\|p\right\|_{g^{-1}}.\quad(\text{Fact 5})$

Fact 10 follows from the following computation:

$\displaystyle\left\|g^{-1}p-g^{\prime-1}q\right\|_{g}$
$\displaystyle\leq\left\|g^{-1}(p-q)+(g^{-1}-g^{\prime-1})q\right\|_{g}\leq\left\|p-q\right\|_{g^{-1}}+\underbrace{O(\delta_{x})\left\|q\right\|_{g^{-1}}}_{\text{Fact
9}}.$

Fact 11 follows from the following computation and Fact 10:

$\displaystyle\left\|\frac{\mathcal{\partial}H}{\mathcal{\partial}v}(x,v)-\frac{\mathcal{\partial}H}{\mathcal{\partial}v}(x^{\prime},v^{\prime})\right\|_{g}$
$\displaystyle=\left\|g^{-1}v-g^{\prime-1}v^{\prime}\right\|_{g}\leq\delta_{v}+O(\delta_{x})\left\|v^{\prime}\right\|_{g^{-1}}.$

For Fact 12, we note that

$\displaystyle\frac{\mathcal{\partial}H}{\mathcal{\partial}x}(x,v)-\frac{\mathcal{\partial}H}{\mathcal{\partial}x}(x^{\prime},v^{\prime})$
$\displaystyle=\left(\nabla f(x)-\nabla f(x^{\prime})\right)$ $\displaystyle\
\
-\frac{1}{2}\left(Dg\left[g^{-1}v,g^{-1}v\right]-Dg^{\prime}\left[g^{\prime-1}v,g^{\prime-1}v\right]\right)+\frac{1}{2}\left(\mathrm{Tr}(g^{-1}Dg)-\mathrm{Tr}(g^{\prime-1}Dg^{\prime})\right)$
$\displaystyle=-\frac{1}{2}\left(\underbrace{Dg\left[g^{-1}v,g^{-1}v\right]-Dg^{\prime}\left[g^{-1}v,g^{-1}v\right]}_{F}+\underbrace{Dg^{\prime}\left[g^{-1}v,g^{-1}v\right]-Dg^{\prime}\left[g^{\prime-1}v,g^{\prime-1}v\right]}_{S}\right)$
$\displaystyle\leavevmode\nobreak\ \leavevmode\nobreak\
+\frac{1}{2}\left(\underbrace{\mathrm{Tr}(g^{-1}Dg-g^{\prime-1}Dg)}_{T}+\underbrace{\mathrm{Tr}(g^{\prime-1}Dg-g^{\prime-1}Dg^{\prime})}_{R}\right)+\left(\nabla
f(x)-\nabla f(x^{\prime})\right).$

For $F$, by the third fact in Lemma 25

$\displaystyle\left\|F\right\|_{g^{-1}}$
$\displaystyle\lesssim\frac{1}{(1-\delta_{x})^{3}}\left\|g^{-1}v\right\|_{g}^{2}\left\|x-x^{\prime}\right\|_{g}=\frac{1}{(1-\delta_{x})^{3}}\left\|v\right\|_{g^{-1}}^{2}\delta_{x}\lesssim\delta_{x}\left\|v\right\|_{g^{-1}}^{2}.$

For $S$, by the fourth fact in Lemma 25

$\displaystyle\left\|S\right\|_{g^{-1}}$
$\displaystyle\lesssim\left\|S\right\|_{g^{\prime-1}}\lesssim\left\|g^{-1}v-g^{\prime-1}v^{\prime}\right\|_{g^{\prime}}\left\|g^{-1}v+g^{\prime-1}v^{\prime}\right\|_{g^{\prime}}$
$\displaystyle\lesssim\left(\left\|v-v^{\prime}\right\|_{g^{\prime-1}}+O(\delta_{x})\left\|v\right\|_{g^{\prime-1}}\right)\left(\left\|v\right\|_{g^{\prime-1}}+\left\|v^{\prime}\right\|_{g^{\prime-1}}\right)$
$\displaystyle\lesssim\left(\delta_{v}+\delta_{x}\left\|v\right\|_{g^{-1}}\right)\left(\left\|v\right\|_{g^{\prime-1}}+\left\|v^{\prime}\right\|_{g^{\prime-1}}\right).$

For $T$, using the stochastic estimator of trace

$\displaystyle\left\|\mathrm{Tr}\left(g^{-1}Dg-g^{\prime-1}Dg\right)\right\|_{g^{-1}}$
$\displaystyle=\max_{\left\|u\right\|_{g}=1}\mathrm{Tr}\left((g^{-1}-g^{\prime-1})Dg[u]\right)$
$\displaystyle=\max_{\left\|u\right\|_{g}=1}\mathrm{Tr}\left(g^{\frac{1}{2}}(g^{-1}-g^{\prime-1})Dg[u]g^{-\frac{1}{2}}\right)$
$\displaystyle=\max_{\left\|u\right\|_{g}=1}\mathbb{E}_{z\sim\mathcal{N}(0,I)}\left[z^{\top}g^{\frac{1}{2}}(g^{-1}-g^{\prime-1})Dg[u]g^{-\frac{1}{2}}z\right]$
$\displaystyle=\max_{\left\|u\right\|_{g}=1}\mathbb{E}Dg\left[u,g^{-\frac{1}{2}}z,(g^{-1}-g^{\prime-1})g^{\frac{1}{2}}z\right]$
$\displaystyle\leq
2\mathbb{E}\max_{\left\|u\right\|_{g}=1}\left\|u\right\|_{g}\left\|g^{-\frac{1}{2}}z\right\|_{g}\underbrace{\left\|\left(g^{-1}-g^{\prime-1}\right)g^{\frac{1}{2}}z\right\|_{g}}_{\text{Fact
9}}$ $\displaystyle\leq
O(\delta_{x})\mathbb{E}\left\|z\right\|_{2}\left\|g^{\frac{1}{2}}z\right\|_{g^{-1}}=O(\delta_{x})\mathbb{E}\left\|z\right\|_{2}^{2}$
$\displaystyle=O(n\delta_{x}).$

For $R$, in a similar way that we bounded $\left\|T\right\|_{g^{-1}}$

$\displaystyle\left\|\mathrm{Tr}\left(g^{\prime-1}Dg-g^{\prime-1}Dg^{\prime}\right)\right\|_{g^{-1}}$
$\displaystyle=\left\|\mathbb{E}_{z\sim\mathcal{N}(0,I)}\left(Dg[g^{\prime-\frac{1}{2}}z,g^{\prime-\frac{1}{2}}z]-Dg^{\prime}[g^{\prime-\frac{1}{2}}z,g^{\prime-\frac{1}{2}}z]\right)\right\|_{g^{-1}}$
$\displaystyle\leq\mathbb{E}\underbrace{\left\|Dg[g^{\prime-\frac{1}{2}}z,g^{\prime-\frac{1}{2}}z]-Dg^{\prime}[g^{\prime-\frac{1}{2}}z,g^{\prime-\frac{1}{2}}z]\right\|_{g^{-1}}}_{\text{Use
Lemma\leavevmode\nobreak\ \ref{lem:sc_facts}}}$ $\displaystyle\leq
O(\delta_{x})\mathbb{E}\left\|g^{\prime-\frac{1}{2}}z\right\|_{g^{\prime}}^{2}\leq
O(\delta_{x})\mathbb{E}\left\|z\right\|_{2}^{2}$
$\displaystyle=O(n\delta_{x}).$

By adding up these bounds, we obtain

$\left\|\frac{\mathcal{\partial}H}{\mathcal{\partial}x}(x,v)-\frac{\mathcal{\partial}H}{\mathcal{\partial}x}(x^{\prime},v^{\prime})\right\|_{g^{-1}}\lesssim\left(\delta_{v}+\delta_{x}\left\|v\right\|_{g^{-1}}\right)\left(\left\|v\right\|_{g^{-1}}+\left\|v^{\prime}\right\|_{g^{-1}}\right)+n\delta_{x}+\left\|\nabla
f(x)-\nabla f(x^{\prime})\right\|_{g^{-1}}.$

∎

We now bound the partial derivatives of $H$ with respect to $x$ and $v$. For
$H_{1}$ and $H_{2}$ given by

$H_{1}(x,v)=f(x)+\frac{1}{2}\log\det g(x)\ \text{ and }\
H_{2}(x,v)=\frac{1}{2}v^{\top}g(x)^{-1}v,$

we recall from (2.1) that

$\displaystyle\frac{\mathcal{\partial}H_{1}}{\partial x}(x,v)$
$\displaystyle=\nabla f(x)+\frac{1}{2}\mathrm{Tr}(g^{-1}Dg),$
$\displaystyle\frac{\partial H_{2}}{\partial x}(x,v)$
$\displaystyle=-\frac{1}{2}Dg[g^{-1}v,g^{-1}v]\ \text{ and }\ \frac{\partial
H_{2}}{\partial v}(x,v)=g^{-1}v.$

######  Lemma 27.

For $x\in\mathbb{\mathcal{M}}$ and $g:=g(x)$, the following inequalities hold.

$\displaystyle\left\|\frac{\partial H_{1}(x,v)}{\partial x}\right\|_{g^{-1}}$
$\displaystyle\leq\left\|\nabla f(x)\right\|_{g^{-1}}+n,$
$\displaystyle\left\|\frac{\partial H_{2}(x,v)}{\partial
v}\right\|_{g}\leq\left\|v\right\|_{g^{-1}}\quad$
$\displaystyle\&\quad\left\|\frac{\partial H_{2}(x,v)}{\partial
x}\right\|_{g^{-1}}\leq\left\|v\right\|_{g^{-1}}^{2}$

###### Proof.

For $\frac{\mathcal{\partial}H_{1}(x,v)}{\partial x}$,

$\displaystyle\left\|\frac{\mathcal{\partial}H_{1}(x,v)}{\partial
x}\right\|_{g^{-1}}$ $\displaystyle\leq\left\|\nabla
f(x)+\frac{1}{2}\mathrm{Tr}(g^{-1}Dg)\right\|_{g^{-1}}\leq\left\|\nabla
f(x)\right\|_{g^{-1}}+\left\|\frac{1}{2}\mathrm{Tr}(g^{-1}Dg)\right\|_{g^{-1}}.$

Note that

$\displaystyle\left\|\frac{1}{2}\mathrm{Tr}(g^{-1}Dg)\right\|_{g^{-1}}$
$\displaystyle=\frac{1}{2}\max_{\left\|u\right\|_{g}=1}\mathrm{Tr}(g^{-1}Dg[u])$

By self-concordance, for any $h\in\mathbb{R}^{n}$ we have $h^{\top}Dg[u]h\leq
2\left\|h\right\|_{g}^{2}$ and thus $Dg[u]\preceq 2g$, resulting in
$g^{-\frac{1}{2}}Dg[u]g^{-\frac{1}{2}}\preceq 2I$. Then

$\displaystyle\mathrm{Tr}(g^{-1}Dg[u])$ $\displaystyle\leq 2\mathrm{Tr}(I)\leq
2n.$

For $\frac{\mathcal{\partial}H_{2}(x,v)}{\partial v}$,

$\displaystyle\left\|\frac{\mathcal{\partial}H_{2}(x,v)}{\partial
v}\right\|_{g}$
$\displaystyle=\left\|g^{-1}v\right\|_{g}=\left\|v\right\|_{g^{-1}}.$

For $\frac{\partial{H}_{2}(x,v)}{\partial x}$,

$\displaystyle\left\|\frac{\mathcal{\partial}H_{2}(x,v)}{\partial
x}\right\|_{g^{-1}}$
$\displaystyle\leq\left\|\frac{1}{2}Dg\left[g^{-1}v,g^{-1}v\right]\right\|_{g^{-1}}\leq\left\|g^{-1}v\right\|_{g}^{2}=\left\|v\right\|_{g^{-1}}^{2},$

where the second step follows from Lemma 25. ∎

###  6.2 Coupling between ideal and discretized RHMC

We bound
$d_{\text{TV}}(\overline{\mathcal{P}}_{x}^{\prime},\,\mathcal{P}_{x})$, the TV
distance between the one-step distributions of the ideal RHMC and the
discretized RHMC without the rejection step. We use
$\overline{\mathcal{P}}_{x}$ to indicate $\overline{\mathcal{P}}_{x}^{\prime}$
for simplicity in this section only. We denote by $p_{x}$ and
$\overline{p}_{x}$ the probability density functions of $\mathcal{P}_{x}$ and
$\overline{\mathcal{P}}_{x}$ respectively. We let $g=g(x)$ and
$g_{t}=g(x_{t})$.

Let us elaborate on our approach. We work with the Euclidean metric this time,
as we find it easier to handle numerical integrators with the Euclidean
representation. As mentioned in (5.2), the one-step distributions
$\mathcal{P}_{x}$ and $\overline{\mathcal{P}}_{x}$ of the ideal and
discretized RHMC on $\mathbb{\mathcal{M}}$ are the pushforwards by $T_{x}$ and
$\overline{T}_{x}$ of the Gaussian distribution of initial velocities on the
tangent space $T_{x}\mathbb{\mathcal{M}}$. Thus, it follows by the change of
variables that for $z=T_{x}(v^{*})=\overline{T}_{x}(v)$ these two probability
densities on the different spaces (one on $\mathbb{\mathcal{M}}$ and another
on $T_{x}\mathbb{\mathcal{M}}$) are related as follows. For $p_{x}^{\ast}$ the
probability density function of $\mathcal{N}(0,g(x))$,

$p_{x}(z)=\sum_{v^{*}:T_{x}(v^{*})=z}\frac{p_{x}^{\ast}(v^{*})}{\left|DT_{x}(v^{*})\right|}\quad\text{and}\quad\overline{p}_{x}(z)=\sum_{v:\overline{T}_{x}(v)=z}\frac{p_{x}^{\ast}(v)}{\left|D\overline{T}_{x}(v)\right|}.$

We aim to couple these $v^{*}$ and $v$ on $T_{x}\mathbb{\mathcal{M}}$. In this
coupling, we can exclude ‘bad’ velocities, as long as such velocities have
small measure. To see this, let $V_{\text{bad}}^{x}$ be a set of bad initial
velocities of measure $\varepsilon<1$ and $V_{\text{good}}^{x}$ be the rest.
Assuming a one-to-one correspondence between $v$ and $v^{*}$ for $v\in
V_{\text{good}}^{x}$, we have that for $x\in\mathbb{\mathcal{M}}_{\rho}$

$\displaystyle d_{\text{TV}}(\overline{\mathcal{P}}_{x},\mathcal{P}_{x})$
$\displaystyle=\sup_{A\subset\mathbb{\mathcal{M}}}\int_{A}(\overline{p}_{x}(z)-p_{x}(z))dz$
$\displaystyle\leq\sup_{A\subset\mathbb{\mathcal{M}}}\int_{A}\left(\sum_{v:\overline{T}_{x}(v)=z}\frac{p_{x}^{\ast}(v)}{\left|D\overline{T}_{x}(v)\right|}-\sum_{v^{*}:T_{x}(v^{*})=z}\frac{p_{x}^{\ast}(v^{*})}{\left|DT_{x}(v^{*})\right|}\right)dz$
$\displaystyle\leq\int_{V_{\text{bad}}^{x}}p_{x}^{*}(v)dv+\sup_{A\subset\mathbb{\mathcal{M}}}\int_{A}\sum_{v\in
V_{\text{good}}^{x}:\overline{T}_{x}(v)=z}\left(\frac{p_{x}^{\ast}(v)}{\left|D\overline{T}_{x}(v)\right|}-\frac{p_{x}^{\ast}(v^{*})}{\left|DT_{x}(v^{*})\right|}\right)dz$
$\displaystyle=P_{x}^{\ast}(V_{\text{bad}}^{x})+\sup_{A\subset\mathbb{\mathcal{M}}}\int_{A}\sum_{v\in
V_{\text{good}}^{x}:\overline{T}_{x}(v)=z}\frac{p_{x}^{\ast}(v)}{\left|D\overline{T}_{x}(v)\right|}\left(1-\frac{p_{x}^{\ast}(v^{*})}{p_{x}^{\ast}(v)}\frac{\left|D\overline{T}_{x}(v)\right|}{\left|DT_{x}(v^{*})\right|}\right)dz,$

where $P_{x}^{\ast}$ is $\mathcal{N}(0,g(x))$, and in the third line we used
the one-to-one correspondence between $v$ and $v^{*}$ to pair them in the
summation. If we show that on $v\in V_{\text{good}}^{x}$ the term of
$1-\frac{p_{x}^{\ast}(v^{*})}{p_{x}^{\ast}(v)}\frac{\left|D\overline{T}_{x}(v)\right|}{\left|DT_{x}(v^{*})\right|}$
is bounded by a small constant (say, $\eta$), then

$\displaystyle d_{\text{TV}}(\overline{\mathcal{P}}_{x},\mathcal{P}_{x})$
$\displaystyle\leq
P_{x}^{\ast}(V_{\text{bad}}^{x})+\eta\sup_{A\subset\mathbb{\mathcal{M}}}\int_{A}\sum_{v\in
V_{\text{good}}^{x}:\overline{T}_{x}(v)=z}\frac{p_{x}^{\ast}(v)}{\left|D\overline{T}_{x}(v)\right|}dz$
$\displaystyle\leq
P_{x}^{\ast}(V_{\text{bad}}^{x})+\eta\int_{V_{\text{good}}^{x}}p_{x}^{\ast}(v)dv$
$\displaystyle=P_{x}^{\ast}(V_{\text{bad}}^{x})+\eta
P_{x}^{\ast}(V_{\text{good}}^{x})$
$\displaystyle\leq\eta+(1-\eta)\varepsilon.$

By taking $\varepsilon$ sufficiently small, we can bound
$d_{\text{TV}}(\overline{\mathcal{P}}_{x},\mathcal{P}_{x})$ smaller then
$1/10$.

For each $x\in\mathbb{\mathcal{M}}$, our bad set $V_{\text{bad}}^{x}$ of
velocities is the union of the following sets:

$\displaystyle V_{1}$
$\displaystyle=\left\\{v\in\mathbb{R}^{n}:\left\|v\right\|_{g^{-1}}>128\sqrt{n}\right\\},$
$\displaystyle V_{2}$
$\displaystyle=\left\\{v\in\mathbb{R}^{n}:\bar{\ell}(\mathrm{Ham}_{x,t}(g(x)^{-1}v))>\frac{1}{2}\bar{\ell}_{0}\right\\},$

and thus

$V_{\text{good}}^{x}=\left\\{v\in\mathbb{R}^{n}:\left\|v\right\|_{g^{-1}}\leq
128\sqrt{n},\
\bar{\ell}(\mathrm{Ham}_{x,t}(g(x)^{-1}v))\leq\frac{1}{2}\bar{\ell}_{0}\right\\}.$
(6.1)

We remark that a velocity $v\in\mathbb{R}^{n}$ should be normalized by
$g(x)^{-1}$ before feeding into $\mathrm{Ham}_{x,t}$, since the domain
$T_{x}\mathbb{\mathcal{M}}$ of $\mathrm{Ham}_{x,t}$ is endowed with the local
metric. Since the standard concentration inequality for the Gaussian
distributions implies that $P_{x}^{*}(V_{1})<\frac{1}{100}$, and the
definition $\bar{\ell}_{0}$ implies that $P_{x}^{\ast}(V_{2})<\frac{1}{100}$,
it follows that $P_{x}^{\ast}(V_{\text{bad}}^{x})<0.02$.

####  6.2.1 Dynamics of ideal and discretized RHMC

We study the dynamics of the ideal and discretized RHMC.

######  Proposition 28.

For $x\in\mathbb{\mathcal{M}}_{\rho}$ and $v,v^{\prime}\in
V_{\text{good}}^{x}$ , let $g:=g(x)$ and $h$ step size satisfying

$\underbrace{h^{2}\leq\frac{10^{-10}}{n+\sqrt{M_{1}}+M_{2}}}_{\small{1}⃝},\underbrace{h^{2}\leq\frac{10^{-10}}{\bar{R}_{1}}}_{\small{2}⃝},\underbrace{C_{x}(x,v)h^{2}\leq\frac{1}{10^{10}}\min\left(1,\frac{\bar{\ell}_{0}}{\bar{\ell}_{1}}\right)}_{\small{3}⃝},\underbrace{C_{x}(x,v)h\leq\frac{\sqrt{n}}{10^{10}}}_{\small{4}⃝}.$

For $t\in[0,h]$, we let $(x_{t},v_{t})$ and $(x_{t}^{\prime},v_{t}^{\prime})$
be the Hamiltonian curves of the ideal RHMC at time $t$ with initial
conditions $(x,v)$ and $(x,v^{\prime})$, respectively. Let
$(\overline{x},\overline{v})$ be the point obtained from RHMC with a sensitive
second-order numerical integrator with the step size $h$ and initial condition
$(x,v)$. Let $\phi\stackrel{{\scriptstyle\mathrm{{\scriptscriptstyle
def}}}}{{=}}\sup_{t\in[0,h]}\left\|x_{t}-x_{t}^{\prime}\right\|_{g(x_{t})}$,
$\psi\stackrel{{\scriptstyle\mathrm{{\scriptscriptstyle
def}}}}{{=}}\sup_{t\in[0,h]}\left\|v_{t}-v_{t}^{\prime}\right\|_{g(x_{t})^{-1}}$
and $\Gamma_{t}(v)\stackrel{{\scriptstyle\mathrm{{\scriptscriptstyle
def}}}}{{=}}g(x_{t})^{-1}(v_{t}-v)$.

  1. 1.

$\left\|x-x_{t}\right\|_{g}=O\left(t\sqrt{n}+t^{2}(n+\sqrt{M_{1}})\right)<\frac{1}{4}$
and $\left\|v-v_{t}\right\|_{g^{-1}}=O\left(t(n+\sqrt{M_{1}})\right)$.

  2. 2.

$(1-o(1))\left\|v_{h}-v_{h}^{\prime}\right\|_{g^{-1}}\leq\psi\leq(1+o(1))\left\|v-v^{\prime}\right\|_{g^{-1}}$.

  3. 3.

$(1-o(1))\left\|T_{x}(v)-T_{x}(v^{\prime})\right\|_{g}\leq\phi\leq(1+o(1))h\psi\leq(1+o(1))h\left\|v-v^{\prime}\right\|_{g}$.

  4. 4.

$\left\|\Gamma_{t}(v)-\Gamma_{t}(v^{\prime})\right\|_{g}\leq
L\left\|v-v^{\prime}\right\|_{g^{-1}}$ for some $L<1/10$.

  5. 5.

For $z=\overline{T}_{x}(v)$, there exists $v^{*}\in\mathbb{R}^{n}$ with
$z=T_{x}(v^{*})$ such that
$\bar{\ell}(\mathrm{Ham}_{x,t}(g(x)^{-1}v^{*}))\leq\bar{\ell}_{0}$ and
$\left\|v-v^{*}\right\|_{g^{-1}}=O\left(\frac{1}{h}\right)\left\|T_{x}(v)-\overline{T}_{x}(v)\right\|_{g}$.
Moreover, there is a one-to-one correspondence between $v$ and $v^{*}$.

Proof of 1. For $0\leq t\leq h$, let us define
$\phi(t):=\left\|x-x_{t}\right\|_{g}$ and
$\psi(t):=\left\|v-v_{t}\right\|_{g^{-1}}$. Note that

$\displaystyle 2\left\|x-x_{t}\right\|_{g}\frac{d\phi(t)}{dt}$
$\displaystyle=\frac{d\phi^{2}(t)}{dt}=\frac{d}{dt}(x_{t}-x)^{\top}g(x_{t}-x)$
$\displaystyle=2\left(\frac{\partial H}{\partial
v}(x_{t},v_{t})\right)^{\top}g(x_{t}-x).$

Hence,

$\left|\left\|x-x_{t}\right\|_{g}\phi^{\prime}(t)\right|=\left|\left(\frac{\partial
H}{\partial
v}(x_{t},v_{t})\right)^{\top}g(x_{t}-x)\right|\leq\left\|\frac{\partial
H}{\partial v}(x_{t},v_{t})\right\|_{g}\left\|x_{t}-x\right\|_{g},$

and $|\phi^{\prime}(t)|\leq\left\|\frac{\partial H}{\partial
v}(x_{t},v_{t})\right\|_{g}$. When $\left\|x-x_{t}\right\|_{g}<\frac{1}{4}$
for $0\leq t\leq h$, since the local norms at $x$ and $x_{t}$ are within a
small constant factor as follows, we have

$\displaystyle\left\|\frac{\partial H}{\partial v}(x_{t},v_{t})\right\|_{g}$
$\displaystyle\leq 2\left\|\frac{\partial H}{\partial
v}(x_{t},v_{t})\right\|_{g_{t}}=2\left\|v_{t}\right\|_{g_{t}^{-1}}\quad(\text{Lemma
\ref{lem:partH_dist}})$ $\displaystyle\leq 4\left\|v_{t}\right\|_{g^{-1}}\leq
4(\left\|v_{t}-v\right\|_{g^{-1}}+\left\|v\right\|_{g^{-1}}),$

and thus

$\phi^{\prime}(t)\leq 10^{3}\sqrt{n}+4\psi(t)\quad\text{if
}\left\|x-x_{t}\right\|_{g}<\frac{1}{4}.$ (6.2)

Similarly, we can obtain

$\displaystyle 2\left\|v-v_{t}\right\|_{g^{-1}}\frac{d\psi}{dt}$
$\displaystyle=\frac{d\psi^{2}(t)}{dt}=\frac{d}{dt}(v_{t}-v)^{\top}g^{-1}(v_{t}-v)$
$\displaystyle=-2\left(\frac{\partial H}{\partial
x}(x_{t},v_{t})\right)^{\top}g^{-1}(v_{t}-v),$

and thus $|\psi^{\prime}(t)|\leq\left\|\frac{\partial H}{\partial
x}(x_{t},v_{t})\right\|_{g^{-1}}$. If $\left\|x-x_{t}\right\|_{g}<\frac{1}{4}$
for $0\leq t\leq h$, then by Lemma 27

$\displaystyle\left\|\frac{\partial H}{\partial
x}(x_{t},v_{t})\right\|_{g^{-1}}$ $\displaystyle\leq
2(\left\|v_{t}\right\|_{g_{t}^{-1}}^{2}+\sqrt{M_{1}}+n)\quad(\because
x\in\mathbb{\mathcal{M}}_{\rho})$ $\displaystyle\leq
2\left(4\left\|v_{t}\right\|_{g^{-1}}^{2}+\sqrt{M_{1}}+n\right)$
$\displaystyle\leq
2\left(8(\left\|v\right\|_{g^{-1}}^{2}+\left\|v_{t}-v\right\|_{g^{-1}}^{2})+\sqrt{M_{1}}+n\right)\quad(\because(a+b)^{2}\leq
2(a^{2}+b^{2}))$ $\displaystyle\leq 10^{6}n+16\psi^{2}(t)+2\sqrt{M_{1}},$

and thus

$\psi^{\prime}(t)\leq 10^{6}n+16\psi^{2}(t)+2\sqrt{M_{1}}\quad\text{if
}\left\|x-x_{t}\right\|_{g}<\frac{1}{4}.$ (6.3)

Now let us solve the coupled inequalities (6.2) and (6.3). When $\psi(t)\leq
1000t(n+\sqrt{M_{1}})$, (6.3) becomes

$\psi^{\prime}(t)\leq 10^{6}n+2\sqrt{M_{1}}+16\cdot
10^{6}t^{2}(n+\sqrt{M_{1}})^{2},$

and this inequality holds up until $h$ satisfying
$\int_{0}^{h}\left(10^{6}n+2\sqrt{M_{1}}+16\cdot
10^{6}t^{2}(n+\sqrt{M_{1}})^{2}\right)dt\leq 1000h(n+\sqrt{M_{1}})$. We can
check that for any $h\leq\frac{10^{-10}}{\sqrt{n+\sqrt{M_{1}}}}$ (i.e.,
condition $\small{1}⃝$) this integral inequality is satisfied. Recall that we
have to ensure that $\phi(t)<\frac{1}{4}$ for $t\leq h$. By substituting
$\psi(t)\leq 1000t(n+\sqrt{M_{1}})$ into (6.2), we have

$\phi^{\prime}(t)\leq 10^{3}\sqrt{n}+4000t(n+\sqrt{M_{1}}),$

as long as $\phi(t)=\left\|x_{t}-x\right\|_{g}<\frac{1}{4}$. It is
straightforward to see that for $t\leq h$ one has

$\displaystyle\phi(t)$ $\displaystyle\leq
10^{3}\sqrt{n}t+2000t^{2}(n+\sqrt{M_{1}})\leq
10^{-7}\sqrt{\frac{n}{n+\sqrt{M_{1}}}}+\frac{2000}{10^{10}}$
$\displaystyle<\frac{1}{10^{6}}.$

Proof of 2. From the first item,
$\left\|x_{t}-x_{t}^{\prime}\right\|_{g}\leq\left\|x-x_{t}\right\|_{g}+\left\|x-x_{t}^{\prime}\right\|_{g}<10^{-5}$.
Due to Lemma 26, we can switch the local norms among $g,g_{t}=g(x_{t})$ and
$g_{t}^{\prime}=g(x_{t}^{\prime})$ by losing a multiplicative constant like
$1+10^{-4}$.

For $\delta_{x}=\left\|x_{t}-x_{t}^{\prime}\right\|_{g_{t}}$ and
$\delta_{v}=\left\|v_{t}-v_{t}^{\prime}\right\|_{g_{t}^{-1}}$,

$\displaystyle\left\|v_{t}-v_{t}^{\prime}\right\|_{g_{t}}$
$\displaystyle=\left\|v-v^{\prime}+\int_{0}^{t}\left(\frac{\mathcal{\partial}H}{\mathcal{\partial}x}(x_{s},v_{s})-\frac{\mathcal{\partial}H}{\mathcal{\partial}x}(x_{s}^{\prime},v_{s}^{\prime})\right)ds\right\|_{g_{t}^{-1}}$
$\displaystyle\leq\left\|v-v^{\prime}\right\|_{g_{t}^{-1}}+O(h)\sup_{t\in[0,h]}\left\|\frac{\mathcal{\partial}H}{\mathcal{\partial}x}(x_{t},v_{t})-\frac{\mathcal{\partial}H}{\mathcal{\partial}x}(x_{t}^{\prime},v_{t}^{\prime})\right\|_{g_{t}^{-1}}$
$\displaystyle\leq\left\|v-v^{\prime}\right\|_{g_{t}^{-1}}+O(h)\sup_{t\in[0,h]}\left((\delta_{v}+\delta_{x}\left\|v_{t}\right\|_{g^{-1}})(\left\|v_{t}\right\|_{g^{-1}}+\left\|v_{t}^{\prime}\right\|_{g^{-1}})+(n+M_{2})\delta_{x}\right),$

where the last step follows from Lemma 26-12. By the first item and
$\small{1}⃝$, we have
$\left\|v_{t}\right\|_{g^{-1}},\left\|v_{t}^{\prime}\right\|_{g^{-1}}\leq
7\sqrt{n+\sqrt{M_{1}}}$ and thus

$\displaystyle\left\|v_{t}-v_{t}^{\prime}\right\|_{g_{t}^{-1}}$
$\displaystyle\leq\left\|v-v^{\prime}\right\|_{g_{t}^{-1}}+O(h)\left(\psi\sqrt{n+\sqrt{M_{1}}}+\phi\left(n+\sqrt{M_{1}}+M_{2}\right)\right)$
$\displaystyle\leq(1+o(1))\left\|v-v^{\prime}\right\|_{g^{-1}}+O(h)\psi\left(\sqrt{n+\sqrt{M_{1}}}+h\left(n+\sqrt{M_{1}}+M_{2}\right)\right),$

where we used $\phi\leq(1+o(1))O(h)\psi$ that we prove in the next item.
Taking the supremum over $t\in[0,h]$, we obtain

$\left(1-O(h)\left(\sqrt{n+\sqrt{M_{1}}}+h\left(n+\sqrt{M_{1}}+M_{2}\right)\right)\right)\psi\leq(1+o(1))\left\|v-v^{\prime}\right\|_{g^{-1}}.$

Taking a sufficiently small constant in $h$ and using $\small{1}⃝$, it follows
that $\psi\leq(1+o(1))\left\|v-v^{\prime}\right\|_{g^{-1}}$.

Proof of 3. By Lemma 26-11,

$\displaystyle\left\|x_{t}-x_{t}^{\prime}\right\|_{g_{t}}$
$\displaystyle=\left\|x+\int_{0}^{t}\frac{\mathcal{\partial}H}{\mathcal{\partial}v}(x_{s},v_{s})ds-\left(x+\int_{0}^{t}\frac{\mathcal{\partial}H}{\mathcal{\partial}v}(x_{s}^{\prime},v_{s}^{\prime})ds\right)\right\|_{g_{t}}$
$\displaystyle\leq
O(h)\sup_{t\in[0,h]}\left\|\frac{\mathcal{\partial}H}{\mathcal{\partial}v}(x_{t},v_{t})-\frac{\mathcal{\partial}H}{\mathcal{\partial}v}(x_{t}^{\prime},v_{t}^{\prime})\right\|_{g_{t}}$
$\displaystyle\leq O(h)\left(\psi+\phi\sqrt{n+\sqrt{M_{1}}}\right).$

By using $\small{1}⃝$ and taking the supremum over $t\in[0,h]$ and a
sufficiently small constant in $h$, we obtain the inequality of
$\phi\leq(1+o(1))O(h)\psi$ as we promised, and the second item implies that

$\phi\leq(1+o(1))O(h)\left\|v-v^{\prime}\right\|_{g^{-1}}.$

Proof of 4. By Lemma 26-12,

$\displaystyle\left\|\Gamma_{t}(v)-\Gamma_{t}(v^{\prime})\right\|_{g_{t}}=$
$\displaystyle\left\|g_{t}^{-1}(v_{t}-v)-g_{t}^{\prime-1}(v_{t}^{\prime}-v^{\prime})\right\|_{g_{t}}$
$\displaystyle\leq$
$\displaystyle\left\|v_{t}-v-v_{t}^{\prime}+v^{\prime}\right\|_{g_{t}^{-1}}+O(1)\underbrace{\left\|x_{t}-x_{t}^{\prime}\right\|_{g_{t}}}_{\leq\phi}\underbrace{\left\|v_{t}^{\prime}-v^{\prime}\right\|_{g_{t}^{-1}}}_{\leq
O(h(n+\sqrt{M_{1}}))}$ $\displaystyle\leq$
$\displaystyle\left\|\int_{0}^{t}\left(\frac{\mathcal{\partial}H}{\mathcal{\partial}x}(x_{s},v_{s})-\frac{\mathcal{\partial}H}{\mathcal{\partial}x}(x_{s}^{\prime},v_{s}^{\prime})\right)ds\right\|_{g_{t}^{-1}}+O\left(h\left(n+\sqrt{M_{1}}\right)\right)\phi$
$\displaystyle\leq$ $\displaystyle
O(h)\sup_{t\in[0,h]}\left\|\frac{\mathcal{\partial}H}{\mathcal{\partial}x}(x_{t},v_{t})-\frac{\mathcal{\partial}H}{\mathcal{\partial}x}(x_{t}^{\prime},v_{t}^{\prime})\right\|_{g_{t}^{-1}}+O\left(h\left(n+\sqrt{M_{1}}\right)\right)\phi.$

We can bound the first term by
$O(h^{2})\left(\sqrt{n+\sqrt{M_{1}}}+h\left(n+\sqrt{M_{1}}+M_{2}\right)\right)\psi$
by following the proof of the second item. Using the second and third items
with the condition $\small{1}⃝$, and taking a sufficiently small constant in
$h$, for some $L<1/10$

$\displaystyle\left\|\Gamma_{t}(v)-\Gamma_{t}(v^{\prime})\right\|_{g}$
$\displaystyle\lesssim
h^{2}\sqrt{n+\sqrt{M_{1}}}+h^{3}\left(n+\sqrt{M_{1}}+M_{2}\right)+h^{2}\left(n+\sqrt{M_{1}}\right)\left\|v-v^{\prime}\right\|_{g^{-1}}$
$\displaystyle\leq L\left\|v-v^{\prime}\right\|_{g^{-1}}.$

Proof of 5. Let $z=\overline{T}_{x}(v)$ for $v\in V_{\text{good}}^{x}$. We
show that the map defined on $u\in
V_{\text{dom}}=\\{v^{\prime}\in\mathbb{R}^{n}:\left\|v-v^{\prime}\right\|_{g^{-1}}\leq
4\sqrt{n}\\}$ by

$\Upsilon(u)=u-\frac{1}{h}gT_{x}(u)+\frac{1}{h}gz,$

is Lipschitz in $u$ with respect to the local norm $g^{-1}$, and then apply
the Banach fixed-point theorem to obtain the unique fixed-point $v^{*}$. Note
that it satisfies $g(T_{x}(v^{*})-\overline{T}_{x}(v))=0$ and thus
$T_{x}(v^{*})=\overline{T}_{x}(v)$.

For Lipschitzness, let $(x_{t},u_{t})$ and $(x_{t}^{\prime},u_{t}^{\prime})$
be the Hamiltonian curves of the ideal RHMC starting from $(x,u)$ and
$(x,u^{\prime})$ for $u,u^{\prime}\in V_{\text{dom}}$, respectively. Observe
that

$\displaystyle\left\|\Upsilon(u)-\Upsilon(u^{\prime})\right\|_{g^{-1}}=\left\|u-u^{\prime}-\frac{1}{h}g\left(T_{x}(u)-T_{x}(u^{\prime})\right)\right\|_{g^{-1}}$
(6.4)
$\displaystyle=\left\|u-u^{\prime}-\frac{1}{h}\int_{0}^{h}g\left(g_{t}^{-1}u_{t}-g_{t}^{-1}u_{t}^{\prime}\right)dt\right\|_{g^{-1}}$
$\displaystyle=\left\|\left(I-\frac{1}{h}g\int_{0}^{h}g_{t}^{-1}dt\right)u-\left(I-\frac{1}{h}g\int_{0}^{h}g_{t}^{\prime-1}dt\right)u^{\prime}-\frac{1}{h}\int_{0}^{h}g\left(\Gamma_{t}(u)-\Gamma_{t}(u^{\prime})\right)dt\right\|_{g^{-1}}$
$\displaystyle=\left\|\underbrace{\frac{1}{h}\left(\int_{0}^{h}(I-gg_{t}^{-1})dt\right)}_{I_{u}}u-\underbrace{\frac{1}{h}\left(\int_{0}^{h}(I-gg_{t}^{\prime-1})dt\right)}_{I_{u^{\prime}}}u^{\prime}-\frac{1}{h}\int_{0}^{h}g\left(\Gamma_{t}(u)-\Gamma_{t}(u^{\prime})\right)dt\right\|_{g^{-1}}$
$\displaystyle\leq\left\|I_{u}(u-u^{\prime})+(I_{u}-I_{u^{\prime}})u^{\prime}\right\|_{g^{-1}}+\frac{1}{h}\left\|\int_{0}^{h}g\left(\Gamma_{t}(u)-\Gamma_{t}(u^{\prime})\right)dt\right\|_{g^{-1}}$
$\displaystyle\leq\left\|I_{u}(u-u^{\prime})\right\|_{g^{-1}}+\left\|(I_{u}-I_{u^{\prime}})u^{\prime}\right\|_{g^{-1}}+\sup_{t\in[0,h]}\left\|\Gamma_{t}(u)-\Gamma_{t}(u^{\prime})\right\|_{g}$
$\displaystyle\leq\underbrace{\left\|I_{u}(u-u^{\prime})\right\|_{g^{-1}}}_{F}+\underbrace{\left\|(I_{u}-I_{u^{\prime}})u^{\prime}\right\|_{g^{-1}}}_{S}+L\left\|u-u^{\prime}\right\|_{g^{-1}},$
(6.5)

where the last inequality follows from the fourth item.

For $F$, let $p=u-u^{\prime}$ and observe that

$\displaystyle\left\|I_{u}p\right\|_{g^{-1}}$
$\displaystyle\leq\frac{1}{h}\int_{0}^{h}\left\|(I-gg_{t}^{-1})p\right\|_{g^{-1}}dt\leq\sup_{t\in[0,h]}\left\|(I-gg_{t}^{-1})p\right\|_{g^{-1}}$
$\displaystyle\leq
O(1)\sup_{t\in[0,h]}\left\|I-g^{\frac{1}{2}}g_{t}^{-1}g^{\frac{1}{2}}\right\|_{2}\left\|p\right\|_{g^{-1}}\leq
O(1)\sup_{t\in[0,h]}\left\|x-x_{t}\right\|_{g}\left\|p\right\|_{g^{-1}}$
$\displaystyle\lesssim\left(h\sqrt{n}+h^{2}\left(n+\sqrt{M_{1}}\right)\right)\left\|u-u^{\prime}\right\|_{g^{-1}}.\quad(\text{First
item})$

For $S$, we can bound it as follows:

$\displaystyle\left\|(I_{u}-I_{u^{\prime}})u^{\prime}\right\|_{g^{-1}}$
$\displaystyle\leq\frac{1}{h}\int_{0}^{h}\left\|g(g_{t}^{\prime-1}-g_{t}^{-1})u^{\prime}\right\|_{g^{-1}}dt\leq\sup_{t\in[0,h]}\left\|(g_{t}^{\prime-1}-g_{t}^{-1})u^{\prime}\right\|_{g}$
$\displaystyle\lesssim\sup_{t\in[0,h]}\left\|(g_{t}^{\prime-1}-g_{t}^{-1})u^{\prime}\right\|_{g_{t}^{\prime}}\lesssim\sup_{t\in[0,h]}\left\|x_{t}-x_{t}^{\prime}\right\|_{g_{t}}\left\|u^{\prime}\right\|_{g_{t}^{\prime}}$
$\displaystyle\lesssim\phi\sqrt{n+\sqrt{M_{1}}}\lesssim
h\sqrt{n+\sqrt{M_{1}}}\left\|u-u^{\prime}\right\|_{g^{-1}}.\quad(\text{Third
item})$

Substituting the bounds on $F$ and $S$ into (6.5) with a small constant in $h$
taken, we can conclude that

$\left\|\Upsilon(u)-\Upsilon(u^{\prime})\right\|_{g^{-1}}<\frac{1}{3}\left\|u-u^{\prime}\right\|_{g^{-1}}.$

Next, we show that the image of $\Upsilon$ is included in $V_{\text{dom}}$.
For $u\in V_{\text{dom}}$,

$\displaystyle\left\|\Upsilon(u)-v\right\|_{g^{-1}}$
$\displaystyle=\left\|u-v-\frac{1}{h}g(T_{x}(u)-\overline{T}_{x}(v))\right\|_{g^{-1}}$
$\displaystyle=\left\|u-v-\frac{1}{h}g(T_{x}(u)-T_{x}(v))+\frac{1}{h}g(T_{x}(v)-\overline{T}_{x}(v))\right\|_{g^{-1}}$
$\displaystyle\leq\left\|u-v-\frac{1}{h}g(T_{x}(u)-T_{x}(v))\right\|_{g^{-1}}+\frac{1}{h}\left\|T_{x}(v)-\overline{T}_{x}(v)\right\|_{g}.$

Repeating the proof for the first item111For $u\in V_{\text{dom}}$, we might
have $\left\|u\right\|_{g^{-1}}\geq 128\sqrt{n}$ though, it is still bounded
above by $132\sqrt{n}$. The proofs of the second to fifth items can be exactly
reproduced for
$V_{\text{relaxed}}\stackrel{{\scriptstyle\mathrm{{\scriptscriptstyle
def}}}}{{=}}\\{v^{\prime}\in\mathbb{R}^{n}:\left\|v^{\prime}\right\|_{g^{-1}}\leq
132\sqrt{n}\\}$, leading to a similar conclusion like
$\left\|\Upsilon(u)-\Upsilon(u^{\prime})\right\|_{g^{-1}}<(\frac{1}{3}+\epsilon)\left\|u-u^{\prime}\right\|_{g^{-1}}$
for a small constant $\epsilon>0$. (see (6.4)), we can bound the first term by
$\frac{1}{2}\left\|u-v\right\|_{g^{-1}}$ and thus by $2\sqrt{n}$, due to $u\in
V_{\text{dom}}$. By Lemma 64 and $\small{4}⃝$, we can bound the second term by

$\displaystyle\frac{1}{h}\left\|T_{x}(v)-\overline{T}_{x}(v)\right\|_{g}$
$\displaystyle\leq\frac{2}{h}d_{g}(T_{x}(v),\overline{T}_{x}(v))\leq\frac{2}{h}C_{x}(x,v)h^{2}\leq
2\sqrt{n}.$

Putting them together, we obtain $\left\|\Upsilon(u)-v\right\|_{g^{-1}}\leq
4\sqrt{n}$.

By the Banach fixed-point theorem, there is a unique fixed point $v^{*}$ of
$\Upsilon$ such that $T_{x}(v^{*})=\overline{T}_{x}(v)$ and

$\left\|\Upsilon(v)-v^{*}\right\|_{g^{-1}}<\frac{1}{3}\left\|v-v^{*}\right\|_{g^{-1}}.$

Moreover,
$\left\|\Upsilon(v)-v^{*}\right\|_{g^{-1}}=\left\|v-v^{*}-\frac{1}{h}g(T_{x}(v)-z)\right\|_{g^{-1}}\geq\left\|v-v^{*}\right\|_{g^{-1}}-\frac{1}{h}\left\|T_{x}(v)-\overline{T}_{x}(v)\right\|_{g}$.
Relating these two inequalities, we obtain

$\left\|v-v^{*}\right\|_{g^{-1}}\lesssim\frac{1}{h}\left\|T_{x}(v)-\overline{T}_{x}(v)\right\|_{g^{-1}}.$

We now show a one-to-one correspondence between $v\in V_{\text{good}}^{x}$ and
$v^{*}$. Let $F_{h}(x,v)=(x_{2},v_{2})$ and $F_{h}(x,v^{*})=(z,v^{\prime})$.
By the reversibility of the Hamiltonian trajectories, we have a one-to-one
correspondence between $v^{*}$ and $v^{\prime}$ in a sense that
$F_{h}(x,v^{*})=(z,v^{\prime})$ and $F_{h}(z,-v^{\prime})=(x,-v^{*})$.
Similarly, we also have a one-to-one correspondence between $v$ and $v_{2}$.
Thus, it suffices to show a one-to-one correspondence between $v_{2}\in
T_{x_{2}}\mathbb{\mathcal{M}}$ and $v^{\prime}\in T_{z}\mathbb{\mathcal{M}}$.

Consider the straight line between $z$ and $x_{2}$. We have that
$\left\|z-x_{2}\right\|_{g}=\left\|\overline{T}_{x}(v)-T_{x}(v)\right\|_{g}\leq
2C_{x}(x,v)h^{2}\leq
10^{-9}\min\left(1,\frac{\bar{\ell}_{0}}{\bar{\ell}_{1}}\right)$ by
$\small{3}⃝$ and
$\bar{\ell}(\mathrm{Ham}_{x_{2},h}(-g(x_{2})^{-1}v_{2}))=\bar{\ell}(\mathrm{Ham}_{x,h}(g(x)^{-1}v))\leq\bar{\ell}_{0}/2$
by the symmetry of $\bar{\ell}$. Due to $\small{2}⃝$, we can apply Lemma 20 to
$x_{2}$ with an initial velocity $-g(x_{2})^{-1}v_{2}$. Thus, a one-to-one
correspondence between $v^{\prime}$ and $v_{2}$ follows, and we also have that
$\bar{\ell}_{0}\geq\bar{\ell}(\mathrm{Ham}_{z,h}(-g(z)^{-1}v^{\prime}))=\bar{\ell}(\mathrm{Ham}_{x,h}(g(x)^{-1}v^{*})$.

####  6.2.2 One-step coupling

As elaborated in Section 6.2, it suffices to prove that for $v\in
V_{\text{good}}^{x}$ the term of
$1-\frac{p_{x}^{\ast}(v^{*})}{p_{x}^{\ast}(v)}\frac{\left|D\overline{T}_{x}(v)\right|}{\left|DT_{x}(v^{*})\right|}$
is bounded by a constant smaller than $1$.

######  Lemma 29.

For $x\in\mathbb{\mathcal{M}}_{\rho}$ and $v\in V_{\text{good}}^{x}$, let step
size $h$ guarantee the sensitivity of a numerical integrator at $(x,v)$, and
satisfy

$\underbrace{h^{2}\leq\frac{10^{-10}}{n+\sqrt{M_{1}}+M_{2}}}_{\small{1}⃝},\underbrace{h^{2}\leq\frac{10^{-10}}{\bar{R}_{1}}}_{\small{2}⃝},\underbrace{C_{x}(x,v)h^{2}\leq\frac{1}{10^{10}}\min\left(1,\frac{\bar{\ell}_{0}}{\bar{\ell}_{1}}\right)}_{\small{3}⃝},\underbrace{C_{x}(x,v)h\leq\frac{1}{10^{10}\sqrt{n}}}_{\small{4}⃝}.$

Then
$d_{\text{TV}}(\overline{\mathcal{P}}_{x},\mathcal{P}_{x})\leq\frac{1}{10}$.

###### Proof.

For given $v\in V_{\text{good}}^{x}$, Proposition 28-5
($\small{1}⃝\sim\small{4}⃝$ required) and the order of the numerical
integrator ensure that there exists $v^{*}\in T_{x}\mathbb{\mathcal{M}}$ such
that $T_{x}(v^{*})=\overline{T}_{x}(v)$ and
$\left\|v-v^{*}\right\|_{g^{-1}}\lesssim
C_{x}(x,v)h\leq\frac{1}{10^{10}\sqrt{n}}$ by $\small{4}⃝$. As $p_{x}^{*}$ is
the probability density function of $\mathcal{N}(0,g(x))$,

$\left|\log\left(\frac{p_{x}^{\ast}(v^{*})}{p_{x}^{\ast}(v)}\right)\right|=\left|\left\|v^{*}\right\|_{g^{-1}}^{2}-\left\|v\right\|_{g^{-1}}^{2}\right|\leq\left\|v^{*}-v\right\|_{g^{-1}}\left(\left\|v\right\|_{g^{-1}}+\left\|v^{*}\right\|_{g^{-1}}\right)\leq\frac{1}{10^{5}},$

and thus the ratio of $\frac{p_{x}^{\ast}(v^{*})}{p_{x}^{\ast}(v)}$ is bounded
below by $0.999$. Also, the sensitivity of the numerical integrator yields
$\frac{\left|D\overline{T}_{x}(v)\right|}{\left|DT_{x}(v^{*})\right|}\geq
0.998$. Hence, for any $v\in V_{\text{good}}^{x}$

$1-\frac{p_{x}^{\ast}(v^{*})}{p_{x}^{\ast}(v)}\frac{\left|D\overline{T}_{x}(v)\right|}{\left|DT_{x}(v^{*})\right|}\leq
0.003,$

and the claim follows. ∎

We finish this section by providing a sufficient condition on the step size
for the sensitivity of a numerical integrator, which we find useful later when
checking the sensitivity of IMM and LM.

######  Proposition 30.

Let $x\in\mathbb{\mathcal{M}}_{\rho}$ and $v\in V_{\text{good}}^{x}$. Let step
size $h$ satisfy $h^{2}\leq\frac{1}{10^{5}\sqrt{n}\bar{R}_{1}}$ in addition to
the step-size conditions in Proposition 28. A numerical integrator
$\overline{T}_{x,h}$ is sensitive at $(x,v)$ if
$\left|D\overline{T}_{x}(v)\right|\geq\frac{(1-10^{-6})h^{n}}{\sqrt{\left|g(x^{\prime})\right|\left|g(x)\right|}}$
for $x^{\prime}=\overline{T}_{x}(v)$.

###### Proof.

By Proposition 28-5, there exists $v^{*}$ such that
$T_{x}(v^{*})=\overline{T}_{x}(v)$ and
$\bar{\ell}(\mathrm{Ham}_{x,t}(g(x)^{-1}v^{*}))\leq\bar{\ell}_{0}$. Let us
estimate $\left|DT_{x}(v^{*})\right|$. Recall that $\mathrm{Ham}_{x,h}$ is the
Hamiltonian map from $T_{x}\mathbb{\mathcal{M}}$ to $\mathbb{\mathcal{M}}$,
where both spaces are endowed with the local metric $g$. Even though $T_{x}$
has the same domain and range, these spaces are endowed with the Euclidean
metric. Therefore, we can relate $T_{x}$ to $\mathrm{Ham}_{x,h}$ by

$T_{x}(v^{*})=\left(\mathrm{id}_{\mathbb{\mathcal{M}}\to\mathbb{R}^{n}}\circ\mathrm{Ham}_{x,h}\circ\mathrm{id}_{\mathbb{R}^{n}\to
T_{x}\mathbb{\mathcal{M}}}\right)(g(x)^{-1}v^{*}),$

where $\mathrm{id}_{\mathbb{\mathcal{M}}\to\mathbb{R}^{n}}$ is the embedding
with transition of metric from $g(x)$ to the Euclidean, and
$\mathrm{id}_{\mathbb{R}^{n}\to T_{x}\mathbb{\mathcal{M}}}$ is the embedding
with transition of metric from the Euclidean to $g(x)$. Note that we have to
normalize $v^{*}$ by $g(x)^{-1}$ before $\mathrm{Ham}_{x,h}$ takes it as
input. Using this formula and the chain rule,

$\displaystyle\left|DT_{x}(v^{*})\right|$
$\displaystyle=\left|D\mathrm{id}_{\mathbb{\mathcal{M}}\to\mathbb{R}^{n}}(x^{\prime})\right|\left|D\mathrm{Ham}_{x,h}(g(x)^{-1}v^{*})\right|\left|D\mathrm{id}_{\mathbb{R}^{n}\to
T_{x}\mathbb{\mathcal{M}}}(g(x)^{-1}v^{*})\right|$
$\displaystyle\leq\left|g(x^{\prime})\right|^{-\frac{1}{2}}\cdot
h^{n}\left(1+\frac{2}{1000}\right)\cdot\left|g(x)\right|^{-1}\cdot\left|g(x)\right|^{\frac{1}{2}}$
$\displaystyle=\frac{h^{n}}{\sqrt{\left|g(x^{\prime})\right|\left|g(x)\right|}}\left(1+\frac{2}{1000}\right),$

where we used Corollary 18 for $\bar{\ell}$ in the second line. Hence, if
$\left|D\overline{T}_{x}(v)\right|\geq\frac{(1-10^{-6})h^{n}}{\sqrt{\left|g(x^{\prime})\right|\left|g(x)\right|}}$,
then

$\frac{\left|D\overline{T}_{x}(v)\right|}{\left|DT_{x}(v^{*})\right|}\geq\frac{1-10^{-6}}{1+0.002}\geq
0.998.$

∎

###  6.3 Bound on rejection probability

Lastly, we bound the rejection probability
$d_{\text{TV}}(\overline{\mathcal{P}}_{x}^{\prime},\overline{\mathcal{P}}_{x})$.

######  Lemma 31.

For $x,x^{\prime}\in\mathbb{\mathcal{M}}$, let $g=g(x)$ and
$g^{\prime}=g(x^{\prime})$. If for some $0<\delta_{x}<1$ and $0<\delta_{v}$ we
have $\left\|x-x^{\prime}\right\|_{g}\leq\delta_{x}$ and