Convergence Rate of Riemannian Hamiltonian Monte Carlo and
Faster Polytope Volume Computation

arXiv:1710.06261v1 [] 17 Oct 2017

Yin Tat Lee∗, Santosh S. Vempala†
October 18, 2017

Abstract
We give the first rigorous proof of the convergence of Riemannian Hamiltonian Monte Carlo,
a general (and practical) method for sampling Gibbs distributions. Our analysis shows that
the rate of convergence is bounded in terms of natural smoothness parameters of an associated
Riemannian manifold. We then apply the method with the manifold defined by the log barrier
function to the problems of (1) uniformly sampling a polytope and (2) computing its volume, the
latter by extending Gaussian cooling to the manifold setting. In both cases, the total number
2
of steps needed is O∗ (mn 3 ), improving the state of the art. A key ingredient of our analysis is
a proof of an analog of the KLS conjecture for Gibbs distributions over manifolds.

Contents
.
.
.
.
.

3
3
5
6
7
7

2 Basics of Hamiltonian Monte Carlo
2.1 Hamiltonian Monte Carlo on Riemannian manifolds . . . . . . . . . . . . . . . . . . .

8
10

3 Convergence of Riemannian Hamiltonian Monte Carlo
3.1 Basics of geometric Markov chains . . . . . . . . . . . . .
3.2 Overlap of one-step distributions . . . . . . . . . . . . . .
3.2.1 Proof Outline . . . . . . . . . . . . . . . . . . . . .
3.2.2 Variation of Hamiltonian curve . . . . . . . . . . .
3.2.3 Local Uniqueness of Hamiltonian Curves . . . . . .
3.2.4 Smoothness of one-step distributions . . . . . . . .
3.3 Convergence bound . . . . . . . . . . . . . . . . . . . . . .

13
13
15
15
16
19
21
23

1 Introduction
1.1 Results . . . .
1.2 Approach and
1.3 Practicality .
1.4 Notation . .
1.5 Organization

. . . . . . . .
contributions
. . . . . . . .
. . . . . . . .
. . . . . . . .

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

4 Improved analysis of the convergence
24
4.1 Improved one-to-one correspondence for Hamiltonian curve . . . . . . . . . . . . . . 25
4.2 Improved smoothness of one-step distributions . . . . . . . . . . . . . . . . . . . . . . 26
∗
†

University of Washington and Microsoft Research, yintat@uw.edu
Georgia Tech, vempala@gatech.edu

1

5 Gibbs sampling on manifolds
29
5.1 Isoperimetry for Hessian manifolds . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29
5.2 Sampling with the log barrier . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
6 Polytope volume computation: Gaussian cooling on manifolds
6.1 Algorithm: cooling schedule . . . . . . . . . . . . . . . . . . . . .
6.2 Correctness of the algorithm . . . . . . . . . . . . . . . . . . . . .
6.2.1 Initial and terminal conditions . . . . . . . . . . . . . . .
6.2.2 Variance of Yx . . . . . . . . . . . . . . . . . . . . . . . .
6.2.3 Main lemma . . . . . . . . . . . . . . . . . . . . . . . . . .
6.3 Volume computation with the log barrier . . . . . . . . . . . . . .

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

32
32
32
32
35
37
38

7 Logarithmic barrier
7.1 Riemannian geometry on ML (G2 ) . . . .
7.2 Hamiltonian walk on ML . . . . . . . . .
7.3 Randomness of the Hamiltonian flow (ℓ0 )
7.4 Parameters R1 , R2 and R3 . . . . . . . . .
7.5 Stability of L2 + L4 + L∞ norm (ℓ1 ) . . .
7.6 Mixing Time . . . . . . . . . . . . . . . .

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

39
39
40
42
46
53
55

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

A Matrix ODE

59

B Concentration

61

C Calculus

61

D Basic definitions of Riemannian geometry
62
D.1 Curvature . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64
D.2 Hessian manifolds . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65

2

1

Introduction

Hamiltonian dynamics provide an elegant alternative to Newtonian mechanics. The Hamiltonian
H, which captures jointly the potential and kinetic energy of a particle, is a function of its position
and velocity. First-order differential equations describe the change in both.
∂H(x, v)
dx
=
,
dt
∂v
dv
∂H(x, v)
=−
.
dt
∂x
As we review in Section 2, these equations preserve the Hamiltonian H.
Riemannian Hamiltonian Monte Carlo (or RHMC) [26, 25][6, 7] is a Markov Chain Monte Carlo
method for sampling from a desired distribution. The target distribution is encoded in the definition
of the Hamiltonian. Each step of the method consists of the following: At a current point x,
1. Pick a random velocity y according to a local distribution defined by x (in the simplest setting,
this is the standard Gaussian distribution for every x).
2. Move along the Hamiltonian curve defined by Hamiltonian dynamics at (x, y) for time (distance) δ.
For a suitable choice of H, the marginal distribution of the current point x approaches the desired
target distribution. Conceptually, the main advantage of RHMC is that it does not require a
Metropolis filter (as in the Metropolis-Hastings method) and its step sizes are therefore not severely
limited even in high dimension.
Over the past two decades, RHMC has become very popular in Statistics and Machine Learning,
being applied to Bayesian learning, to evaluate expectations and likelihood of large models by sampling from the appropriate Gibbs distribution, etc. It has been reported to significantly outperform
other known methods [2, 25] and much effort has been made to make each step efficient by the use
of numerical methods for solving ODEs.
In spite of all these developments and the remarkable empirical popularity of RHMC, analyzing
its rate of convergence and thus rigorously explaining its success has remained an open question.

1.1

Results

In this paper, we analyze the mixing rate of Hamiltonian Monte Carlo for a general function f as a
Gibbs sampler, i.e., to generate samples from the density proportional to e−f (x) . The corresponding
Hamiltonian is H(x, v) = f (x) + 21 log((2π)n det g(x)) + 12 v T g(x)−1 v for some metric g. We show
that for x in a compact manifold, the conductance of the Markov chain is bounded in terms of
a few parameters of the metric g and the function f . The parameters and resulting bounds are
given in Corollary 28 and Theorem 30. Roughly speaking, the guarantee says that Hamiltonian
Monte Carlo mixes in polynomial time for smooth Hamiltonians. We note that these guarantees
use only the smoothness and Cheeger constant (expansion) of the function, without any convexity
type assumptions. Thus, they might provide insight in nonconvex settings where (R)HMC is often
applied.
We then focus on logconcave densities in Rn , i.e., f (x) is a convex function. This class of
functions appears naturally in many contexts and is known to be sampleable in polynomial-time
given access to a function value oracle. For logconcave densities, the current fastest sampling
algorithms use n4 function calls, even for uniform sampling [20, 23], and n2.5 oracle calls given a
3

warm start after appropriate rounding (linear transformation) [14]. In the prototypical setting of
uniform sampling from a polytope Ax ≥ b, with m inequalities, the general complexity is no better,
with each function evaluation taking O(mn) arithmetic operations, for an overall complexity of
n4 · mn = mn5 in the worst case and n2.5 · mn after rounding from a warm start. The work of
Kannan and Narayanan [11] gives an algorithm of complexity mn2 · mnω−1 from an arbitrary start
and mn · mnω−1 from a warm start (here ω is the matrix multiplication exponent), which is better
than the general case when the number of facets m is not too large. This was recently improved
to mn0.75 · mnω−1 from a warm start [13]; the subquadratic complexity for the number of steps is
significant since all known general oracle methods cannot below a quadratic number of steps. The
leading algorithms and their guarantees are summarized in Table 1.
Year
1997 [10]
2003 [21]
2009 [11]
2016 [13]
2016 [15]
This paper

Algorithm
Ball walk#
Hit-and-run#
Dikin walk
Geodesic walk
Ball walk#
RHMC

Steps
n3
n3
mn
3
mn 4
n2.5
2
mn 3

Cost per step
mn
mn
mnω−1
mnω−1
mn
mnω−1

Table 1: The complexity of uniform polytope sampling from a warm start, where each step of
e
every algorithm uses O(n)
bit of randomness. The entries marked # are for general convex bodies
presented by oracles, while the rest are for polytopes.
In this paper, using RHMC, we improve the complexity of sampling polytopes. In fact we do
this for a general family of Gibbs distributions, of the form e−αφ(x) where φ(x) is a convex function
over a polytope. When φ(x) is the standard logarithmic barrier function and g(x) is its Hessian,
1

1

we get a sampling method that mixes in only n 6 m 2 +

1

2

2

n5 m5

1
1
α 5 +m− 5

+

n3
α+m−1

steps from a warm start!

When α = 1/m, the resulting distribution is very close to uniform over the polytope.
Theorem 1. Let φ be the logarithmic barrier for a polytope M with m constraints and n variables.
Hamiltonian Monte Carlo applied to the function f = exp(−αφ(x)) and the metric given by ∇2 φ
with appropriate step size mixes in
!
1
2
1
3
3n3
1
1
n
m
e
+ m2 n6
O
+
α + m−1 α 31 + m− 13
steps where each step is the solution of a Hamiltonian ODE.

In recent independent work, Mangoubi and Smith [24] analyze Euclidean HMC in the oracle setting, i.e., assuming an oracle for evaluating φ. Their analysis formally gives a dimension-independent
convergence rate based on certain regularity assumptions such as strong convexity and smoothness
of the Hamiltonian H. Unfortunately, these assumptions do not hold for the polytope sampling
problem.
An important application of sampling is integration. The complexity of integration for general
logconcave functions is also n4 oracle calls. For polytopes, the most natural questions is computing
its volume. For this problem, the current best complexity is n4 · mn, where the factor of O(mn) is
the complexity of checking membership in a polytope. Thus, even for explicitly specified polytopes,
the complexity of estimating the volume from previous work is asymptotically the same as that
4

for a general convex body given by a membership oracle. Here we obtain a volume algorithm with
2
complexity mn 3 · mnω−1 , improving substantially on previous algorithms. The volume algorithm
is based using Hamiltonian Monte Carlo for sampling from a sequence of Gibbs distributions over
polytopes. We remark that in the case when m = O(n)1 , the final complexity is o(n4 ) arithmetic
operations, improving by more than a quadratic factor in the dimension over the previous best
complexity of Õ(n6 ) operations for arbitrary polytopes. These results and prior developments are
given in Table 2.
Year
1989 [5]
1989-93 [17, 4, 1, 18, 19]
1997 [10]
2003 [22]
2015 [3]
This paper

Algorithm
DFK
many improvements
DFK, Speedy walk, isotropy
Annealing, hit-and-run
Gaussian Cooling*
RHMC + Gaussian Cooling

Steps
n23
n7
n5
n4
n3
2
mn 3

Cost per step
mn
mn
mn
mn
mn
mnω−1

e
Table 2: The complexity of volume estimation, each step uses O(n)
bit of randomness, all except
*
the last for general convex bodies (the result marked is for well-rounded convex bodies). The
current paper applies to general polytopes, and is the first improvement utilizing their structure.
Theorem 2. For any polytope P = {x : Ax ≥ b} with m constraints and n variables, and any
ε > 0, the Hamiltonian
volume algorithm estimates the volume of P to within 1 ± ε multiplicative


2
−2
e
3
factor using O mn ε
steps where each step consists of solving a first-order ODE and takes time


e mnω−1 LO(1) logO(1) 1 and L is the bit complexity2 of the polytope.
O
ε
A key ingredient in the analysis of RHMC is a new isoperimetric inequality for Gibbs distributions over manifolds. This inequality can be seen as an evidence of a manifold version of the
KLS hyperplane conjecture. For the family of Gibbs distributions induced by convex functions with
convex Hessians, the expansion is within a constant factor of that of a hyperplane cut. This result
might be of independent interest.

1.2

Approach and contributions

Traditional methods to sample from distributions in Rn are based on random walks that take straight
line steps (grid walk, ball walk, hit-and-run). While this leads to polynomial-time convergence for
logconcave distributions, the length of each step has to be small due to boundary effects, and a
Metropolis filter (rejection sampling) has to be applied to ensure the limiting
 distribution is the
desired one. These walks cannot afford a step of length greater than δ = O √1n for a distribution
in isotropic position, and take a quadratic number of steps even for the hypercube. The Dikin walk
for polytopes [11], which explicitly takes into account the boundary of polytope at each step, has a
varying step size, but still runs into similar issues and the bound on its convergence rate is O(mn)
for a polytope with m facets.
1

We suspect that the LS barrier [12] might be used to get a faster algorithm even in the regime even if m is
sub-exponential. However, our proof requires a delicate estimate of the fourth derivative of the barrier functions.
Therefore, such a result either requires a new proof or a unpleasantly long version of the current proof.
2
L = log(m + dmax + kbk∞ ) where dmax is the largest absolute value of the determinant of a square sub-matrix of
A.

5

In a recent paper [13], we introduced the geodesic walk. Rather than using straight lines in
Euclidean space, each step of the walk is along a geodesic (locally shortest path) of a Riemannian
metric. More precisely, each step first makes a deterministic move depending on the current point
(drift), then moves along a geodesic in a random initial direction and finally uses a Metropolis filter.
Each step can be computed by solving a first-order ODE. Due to the combination of drift and
geodesic, the local 1-step distributions are smoother than that of the Dikin walk and larger steps
can be taken while keeping a bounded rejection probability for the filter. For sampling polytopes,
3
the manifold/metric defined by the standard log barrier gives a convergence rate of mn 4 , going
below the quadratic (or higher) bound of all previous sampling methods.
One major difficulty with geodesic walk is ensuring the stationary distribution is uniform. For
high dimensional problems, this necessitates taking a sufficiently small step size and then rejecting
some samples according to the desired transition probabilities according to Metropolis filter. Unfortunately, computing these transition probabilities can be very expensive. For the geodesic walk,
it entails solving an n × n size matrix ODE.
Hamiltonian Monte Carlo bears some similarity to the geodesic walk — each step is a random
(non-linear) curve. But the Hamiltonian-preserving nature of the process obviates the most expensive ingredient, Metropolis filter. Due to this, the step size can be made longer, and as a result we
2
obtain a faster sampling algorithm for polytopes that mixes in mn 3 steps (the per-step complexity
remains essentially the same, needing the solution of an ODE).
To get a faster algorithm for volume computation, we extends the analysis to a general family of
Gibbs distributions, including f (x) = e−αφ(x) where φ(x) is the standard log-barrier and α > 0. We
show that the smoothness we need for the sampling corresponding to a variant of self-concordance
defined in Definition 46. Furthermore, we establish an isoperimetric inequality for this class of
functions. This can be viewed as an extension of the KLS hyperplane conjecture from Euclidean to
Riemannian metrics (the analogous case in Euclidean space to what we prove here is the isoperimetry
of the Gaussian density function multiplied by any logconcave function, a case for which the KLS
conjecture holds). The mixing rate for this family of functions is sublinear for α = Ω(1).
Finally, we study the Gaussian Cooling schedule of [3]. We show that in the manifold setting,
2
the Gaussian distribution e−kxk /2 can be replaced by e−αφ(x) . Moreover, the speed of Gaussian
Cooling depends on the “thin-shell” constant of the manifold and classical self-concordance of φ.
Combining all of these ideas, we obtain a faster algorithm for polytope volume computation.
The resulting complexity of polytope volume computation is the same as that of sampling uniformly
2
from a warm start: mn 3 steps. To illustrate the improvement, for polytopes with m = O(n) facets,
5
the new bound is n 3 while the previous best bound was n4 .

1.3

Practicality

From the experiments, the ball walk/hit-and-run seem to mix in n2 steps, the geodesic walk seems
to mix in sublinear number of steps (due to the Metropolis filter bottleneck) and RHMC seems to
mix in only polylogarithmic number of steps. One advantage of RHMC compared to the geodesic
walk is that it does not require the expensive Metropolis filter that involves solving n × n matrix
ODEs. In the future, we plan to do an empirical comparison study of different sampling algorithms.
We are hopeful that using RHMC we might finally be able to sample from polytopes in millions of
dimensions after more than three decades of research on this topic!

6

1.4

Notation

Throughout the paper, we use lowercase letter for vectors and vector fields and uppercase letter for
d
matrices and tensors. We use ek to denote coordinate vectors. We use dt
for the usual derivative,
df (c(t))
∂
is the derivative of some function f along a curve c parametrized by t, we use ∂v
for
e.g.
dt
k
th
the usual partial derivative. We use D f (x)[v1 , v2 , · · · , vk ] for the k directional derivative of f at
x along v1 , v2 , · · · , vk . We use ∇ for the usual gradient and the connection (manifold derivative,
defined in Section D which takes into account the local metric), Dv for the directional derivative
of a vector with respect to the vector (or vector field) v (again, defined in Section D), and Dt if
the curve v(t) is clear from the context. We use g for the local metric. Given a point x ∈ M , g
is a matrix with entries gij . Its inverse has entries g ij . Also, n is the dimension, m the number of
inequalities. We use dT V for the total variation (or L1 ) distance between two distributions.

1.5

Organization

In Section 2, we define the Riemannian Hamiltonian Monte Carlo and study its basic properties
such as time-reversibility. In Section 3, we give the the first convergence rate analysis of RHMC.
However, the convergence rate is weak for the sampling applications (it is polynomial, but not better
than previous methods). In Section 4, we introduce more parameters and use them to get a tighter
analysis of RHMC. In Section 5, we study the isoperimetric constant of f (x) = e−αφ(x) under the
metric ∇2 φ(x). In Section 6, we study the generalized Gaussian Cooling schedule and its relation
to the thin-shell constant. Finally, in Section 7, we compute the parameters we need for the log
barrier function.

7

2

Basics of Hamiltonian Monte Carlo

In this section, we define the Hamiltonian Monte Carlo method for sampling from a general distribution e−H(x,y) . Hamiltonian Monte Carlo uses curves instead of straight lines and this makes
the walk time-reversible even if the target distribution is not uniform, with no need for a rejection
sampling step. In contrast, classical approaches such as the ball walk require an explicit rejection
step to converge to a desired stationary distribution.
Definition 3. Given a continuous, twice-differentiable function H : M×Rn ⊂ Rn ×Rn → R (called
the Hamiltonian, which often corresponds to the total energy of a system) where M is the x domain
of H, we say (x(t), y(t)) follows a Hamiltonian curve if it satisfies the Hamiltonian equations
∂H(x, y)
dx
=
,
dt
∂y
dy
∂H(x, y)
=−
.
dt
∂x

(2.1)

def

We define the map Tδ (x, y) = (x(δ), y(δ)) where the (x(t), y(t)) follows the Hamiltonian curve with
the initial condition (x(0), y(0)) = (x, y).
Hamiltonian Monte Carlo is the result of a sequence of randomly generated Hamiltonian curves.
Algorithm 1: Hamiltonian Monte Carlo
Input: some initial point x(1) ∈ M.
for i = 1, 2, · · · , T do
R
1
(k)
Sample y (k+ 2 ) according to e−H(x ,y) /π(x(k) ) where π(x) = Rn e−H(x,y) dy.
1
With probability 12 , set (x(k+1) , y (k+1) ) = Tδ (x(k) , y (k+ 2 ) ).
1
Otherwise, (x(k+1) , y (k+1) ) = T−δ (x(k) , y (k+ 2 ) ).
end
Output: (x(T +1) , y (T +1) ).

Lemma 4 (Energy Conservation). For any Hamiltonian curve (x(t), y(t)), we have that
d
H(x(t), y(t)) = 0.
dt
Proof. Note that
∂H dx ∂H dy
∂H ∂H
∂H ∂H
d
H(x(t), y(t)) =
+
=
−
= 0.
dt
∂x dt
∂y dt
∂x ∂y
∂y ∂x

Lemma 5 (Measure Preservation). For any t ≥ 0, we have that
det (DTt (x, y)) = 1
where DTt (x, y) is the Jacobian of the map Tt at the point (x, y).

8

Proof. Let (x(t, s), y(t, s)) be a family of Hamiltonian curves given by Tt (x + sdx , y + sdy ). We write
u(t) =

∂
∂
x(t, s)|s=0 , v(t) =
y(t, s)|s=0 .
∂s
∂s

By differentiating the Hamiltonian equations (2.1) w.r.t. s, we have that
du
∂ 2 H(x, y)
∂ 2 H(x, y)
=
u+
v,
dt
∂y∂x
∂y∂y
dv
∂ 2 H(x, y)
∂ 2 H(x, y)
=−
u−
v,
dt
∂x∂x
∂x∂y
(u(0), v(0)) = (dx , dy ).
This can be captured by the following matrix ODE
∂ 2 H(x(t),y(t))
∂y∂x
2
− ∂ H(x(t),y(t))
∂x∂x

dΦ
=
dt

∂ 2 H(x(t),y(t))
∂y∂y
2
− ∂ H(x(t),y(t))
∂x∂y

!

Φ(t)

dx
dy



Φ(0) = I
using the equation
DTt (x, y)



dx
dy



=



u(t)
v(t)



= Φ(t)



.

Therefore, DTt (x, y) = Φ(t). Next, we observe that


d
−1 d
log det Φ(t) = Tr Φ(t)
Φ(t) = Tr
dt
dt

∂ 2 H(x(t),y(t))
∂y∂x
2
− ∂ H(x(t),y(t))
∂x∂x

∂ 2 H(x(t),y(t))
∂y∂y
2
− ∂ H(x(t),y(t))
∂x∂y

!

= 0.

Hence,
det Φ(t) = det Φ(0) = 1.

Using the previous two lemmas, we now show that Hamiltonian Monte Carlo indeed converges
to the desired distribution.
Lemma 6 (Time reversibility). Let px (x′ ) denote the probability density of one step of the Hamiltonian Monte Carlo starting at x. We have that
π(x)px (x′ ) = π(x′ )px′ (x)
R
for almost everywhere in x and x′ where π(x) = Rn e−H(x,y) dy.

Proof. Fix x and x′ . Let Fδx (y) be the x component of Tδ (x, y). Let V+ = {y : Fδx (y) = x′ } and
x (x) = x′ )}. Then,
V− = {y : F−δ
π(x)px (x′ ) =

1
2

Z

y∈V+

1
e−H(x,y)
 +
x
2
det DFδ (y)

Z

y∈V−

e−H(x,y)
 .
x (y)
det DF−δ

We note that this formula assumed that DFδx is invertible. Sard’s theorem showed that Fδx (N ) is
def
measure 0 where N = {y : DFsx (y) is not invertible}. Therefore, the formula is correct except for
a measure zero subset.
9

By reversing time for the Hamiltonian curve, we have that for the same V± ,
Z
Z
′ ′
′ ′
1
e−H(x ,y )
e−H(x ,y )
1


+
π(x′ )px′ (x) =
x′ (y ′ )
2 y∈V+ det DF−δ
2 y∈V− det DFδx′ (y ′ )

(2.2)

where y ′ denotes the y component of Tδ (x, y) and T−δ (x, y) in the first
 and second
 sum respectively.
A B
We compare the first terms in both equations. Let DTδ (x, y) =
. Since Tδ ◦ T−δ = I
C D
and Tδ (x, y) = (x′ , y ′ ), the inverse function theorem shows that DT−δ (x′ , y ′ ) is the inverse map of
DTδ (x, y). Hence, we have that


−1 
A B
· · · −A−1 B(D − CA−1 B)−1
′ ′
DT−δ (x , y ) =
.
=
···
···
C D
′

x (y ′ ) = −A−1 B(D − CA−1 B)−1 . Hence, we have that
Therefore, we have that Fδx (y) = B and F−δ



−1
x′
=
det DF−δ
(y ′ ) = det A−1 det B det D − CA−1 B

|det B|

 .
A B
det
C D




A B
= 1 (Lemma 5), we have that
C D


x′
det DF−δ
(y ′ ) = |det (DFδx (y))| .

Using that det (DTt (x, y)) = det

Hence, we have that
1
2

Z

y∈V+

1
e−H(x,y)
 =
x
2
det DFδ (y)
=

1
2

−H(x′ ,y ′ )

Z
Z

y∈V+

e−H(x,y)

x′ (y ′ )
det DF−δ
′

y∈V+

′

e−H(x ,y )

x′ (y ′ )
det DF−δ

where we used that e−H(x,y) = e
(Lemma 4) at the end.
For the second term in (2.2), by the same calculation, we have that
Z
Z
′ ′
1
e−H(x,y)
e−H(x ,y )
1

 =
x (y)
2 y∈V− det DF−δ
2 y∈V+ det DFδx′ (y ′ )
Combining both terms we have the result.

The main challenge in analyzing Hamiltonian Monte Carlo is to bound its mixing time.

2.1

Hamiltonian Monte Carlo on Riemannian manifolds

Suppose we want to sample from the distribution e−f (x) . We define the following energy function
H:
1
1
def
(2.3)
H(x, v) = f (x) + log((2π)n det g(x)) + v T g(x)−1 v.
2
2
One can view x as the location and v as the velocity. The following lemma shows that the first
variable x(t) in the Hamiltonian curve satisfies a second-order differential equation. When we view
the domain M as a manifold, this equation is simply Dt dx
dt = µ(x), namely, x acts like a particle
under the force field µ. (For relevant background on manifolds, we refer the reader to Appendix D).
10

Lemma 7. In Euclidean coordinates, The Hamiltonian equation for (2.3) can be rewritten as
Dt

dx
=µ(x),
dt

dx
(0) ∼N (0, g(x)−1 )
dt


where µ(x) = −g(x)−1 ∇f (x) − 12 g(x)−1 Tr g(x)−1 Dg(x) and Dt is the Levi-Civita connection on
the manifold M with metric g.
Proof. From the definition of the Hamiltonian curve, we have that
dx
= g(x)−1 v
dt
 1 dx T
dv
1 
dx
= −∇f (x) − Tr g(x)−1 Dg(x) +
Dg(x) .
dt
2
2 dt
dt

Putting the two equations together, we have that

dx
dv
d2 x
= − g(x)−1 Dg(x)[ ]g(x)−1 v + g(x)−1
2
dt
dt
dt

 1
dx dx
1
dx T
dx
= − g(x)−1 Dg(x)[ ]
− g(x)−1 ∇f (x) − g(x)−1 Tr g(x)−1 Dg(x) + g(x)−1
Dg(x) .
dt dt
2
2
dt
dt

Hence,

T


dx dx 1
1
d2 x
dx
−1
−1 dx
+
g(x)
Dg(x)[
]
−
g(x)
= − g(x)−1 ∇f (x) − g(x)−1 Tr g(x)−1 Dg(x) .
Dg(x)
2
dt
dt dt
2
dt
dt
2
(2.4)

Using the formula of Christoffel symbols
Dt

dx
d2 x X dxi dxj k
= 2 +
Γ ek
dt
dt
dt dt ij

where

Γkij =

1 X kl
g (∂j gli + ∂i glj − ∂l gij ),
2
l

ijk

we have that
Dt

X dxi dxj
d2 x 1
dx
= 2 + g(x)−1
(∂j gli + ∂i glj − ∂l gij )el
dt
dt
2
dt dt
ijl

=

T
d2 x
dx dx 1
dx
−1
−1 dx
+
g(x)
Dg(x)[
]
−
g(x)
Dg(x) .
dt2
dt dt
2
dt
dt

Putting this into (2.4) gives
Dt



dx
1
= − g(x)−1 ∇f − g(x)−1 Tr g(x)−1 Dg(x) .
dt
2

Motivated by this, we define the Hamiltonian map as the first component of the Hamiltonian
dynamics operator T defined earlier. For the reader familiar with Riemannian geometry, this is
similar to the exponential map (for background, see Appendix D).
11

Definition 8. Let Hamx,δ (vx ) = γ(δ) where γ(t) be the solution of the Hamiltonian equation
′
Dt dγ
dt = µ with initial conditions γ(0) = x and γ (0) = vx . We also denote Hamx,1 (vx ) by Hamx (vx ).
We now give two examples of Hamiltonian Monte Carlo.
Example 9. When g(x) = I, the Hamiltonian curve acts like stochastic gradient descent for the
function f with each random perturbation drawn from a standard Gaussian.
Dt

dx
= −∇f (x).
dt

When g(x) = ∇2 f (x), the Hamiltonian curve acts like a stochastic Newton curve for the function
f + ψ:
−1
dx
Dt
= − ∇2 f (x)
∇(f (x) + ψ(x))
dt

where the volumetric function ψ(x) = log det ∇2 f (x).

Next we derive a formula for the transition probability in Euclidean coordinates.
Lemma 10. For any x ∈ M ⊂ Rn and s > 0, the probability density of the 1-step distribution from
x is given by
s


X
1
det (g(y))
2
−1
exp − kvx kx
(2.5)
px (y) =
|det(DHamx,δ (vx ))|
(2π)n
2
vx :Hamx,δ (vx )=y

where DHamx,δ (vx ) is the Jacobian of the Hamiltonian map Hamx,δ .
Proof. We prove the formula by separately considering each vx ∈ Tx M s.t. Hamx,δ (vx ) = y, then
summing up. In the tangent space Tx M, the point vx follows a Gaussian step. Therefore, the
probability density of vx in Tx M is as follows:


1
1
2
exp
−
kv
k
pTxx M (vx ) =
x x .
2
(2π)n/2
Let y = Hamx,δ (vx ) and F : Tx M → Rn be defined by F (v) = idM→Rn ◦ Hamx,δ (v). Here Rn is
the same set as M but endowed with the Euclidean metric. Hence, we have
DF (vx ) = DidM→Rn (y)DHamx,δ (vx ).
The result follows from px (y) = |det(DF (vx ))|−1 pxTx M (vx ) and
det DF (vx ) = det (DidM→Rn (y)) det (DHamx,δ (vx ))
= det(g(y))−1/2 det (DHamx,δ (vx )) .

12

3

Convergence of Riemannian Hamiltonian Monte Carlo

Hamiltonian Monte Carlo is a Markov chain on a manifold whose stationary stationary distribution
has density q(x) proportional to exp(−f (x)). We will bound the conductance of this Markov chain
and thereby its mixing time to converge to the stationary distribution. Bounding conductance
involves showing (a) the induced metric on the state space satisfies a strong isoperimetric inequality
and (b) two points that are close in metric distance are also close in probabilistic distance, i.e., the
one-step distributions from two nearby points have large overlap. In this section and the next, we
present general conductance bounds using parameters determined by the associated manifold. In
Section 7, we bound these parameters for the manifold corresponding to the logarithmic barrier in
a polytope.

3.1

Basics of geometric Markov chains

For completeness, we will discuss some standard techniques in geometric random walks in this subsection. For a Markov chain with state space M, stationary distribution q and next step distribution
pu (·) for any u ∈ M, the conductance of the Markov chain is
R
def
S pu (M \ S)dq(u)
φ = inf
.
S⊂M min {q(S), q(M \ S)}
The conductance of an ergodic Markov chain allows us to bound its mixing time, i.e., the rate of
convergence to its stationary distribution, e.g., via the following theorem of Lovász and Simonovits.
Theorem 11 ([19]). Let qt be the distribution of the current point after t steps of a Markov chain
with stationary distribution q and conductance at least φ, starting from initial distribution q0 . For
any ε > 0,
s

t
dq0 (x)
φ2
1
.
Ex∼q0
dT V (qt , q) ≤ ε +
1−
ε
dq(x)
2
Definition 12. The isoperimetry of a metric space M with target distribution q is
R
d(S,x)≤δ q(x)dx − q(S)
def
ψ = inf min
δ>0 S⊂M δ min {q(S), q(M \ S)}
where d is the shortest path distance in M.
The proof of the following theorem follows the standard outline for geometric random walks (see
e.g., [29]).
Lemma 13. Given a metric space M and a time-reversible Markov chain p on M with stationary
distribution q. Fix any r > 0. Suppose that for any x, y ∈ M with d(x, z) < r, we have that
dT V (px , py ) ≤ 0.9. Then, the conductance of the Markov chain is Ω(rψ).
Proof. Let S be any measurable subset of M. Then our goal is to bound the conductance of the
Markov chain
R
S px (M \ S) dq(x)
= Ω (rψ) .
min {q(S), q(M \ S)}
R
R
Since the Markov chain is time-reversible (For any two subsets A, B, A px (B) dq(x) = B px (A) dq(x)),
we can write the numerator of the left hand side above as
!
Z
Z
1
px (S) dq(x) .
px (M \ S) dq(x) +
2
M\S
S
13

Define
S1 = {x ∈ S : px (M \ S) < 0.05}

S2 = {x ∈ M \ S : px (S) < 0.05}

S3 = M \ S1 \ S2 .

Without
loss of generality, we can assume that q(S1 ) ≥ (1/2)q(S) and q(S2 ) ≥ (1/2)q(M \ S) (if
R
not, S px (M \ S) dq(x) = Ω(1) and hence the conductance is Ω(1).)
Next, we note that for any two points x ∈ S1 and y ∈ S2 , dT V (px , py ) > 0.9. Therefore, by the
assumption, we have that d(x, y) ≥ r. Therefore, by the definition of ψr , we have that
Z
q(x)dx − q(S1 )
q(S3 ) ≥
d(S1 ,x)≤r

≥ rψ min {q(S1 ), q(M \ S1 )}

≥ rψ min {q(S1 ), q(S2 )} .
Going back to the conductance,
1
2

Z

S

px (M \ S) dq(x) +

Z

!

px (S) dq(x)
M\S

1
≥
2

Z

(0.05)dq(x)

S3

= Ω (rψ) min{q(S1 ), q(S2 )}
= Ω (rψ) min{q(S), q(M \ S)}.
Therefore, the conductance of the Markov chain is Ω(rψ).
Combining Theorem 13 and Lemma 13 gives the following result for bounding mixing time of
general geometric random walk.
Lemma 14. Given a metric space M and a time-reversible Markov chain p on M with stationary
distribution q. Suppose that there exist r > 0 and ψ > 0 such that
1. For any x, y ∈ M with d(x, z) < r, we have that dT V (px , py ) ≤ 0.9.
2. For any S ⊂ M, we have that
Z

0<d(S,x)≤r

q(x)dx ≥ rψ min {q(S), q(M \ S)} .

Let qt be the distribution of the current point after t steps of a Markov chain with stationary distribution q starting from initial distribution q0 . For any ε > 0,
s
t
1
dq0 (x)
Ex∼q0
1 − Ω(r 2 ψ 2 ) .
dT V (qt , q) ≤ ε +
ε
dq(x)

14

3.2

Overlap of one-step distributions

The mixing of the walk depends on smoothness parameters of the manifold and the functions f, g
used to define the Hamiltonian. Since each step of our walk involves a Gaussian vector, many
smoothness parameters depend on choices of the random vector. Formally, let γ be the Hamiltonian curve used in a step of Hamiltonian Monte Carlo. In the analysis, we need a large fraction
of Hamiltonian curves from any point on the manifold to be well-behaved. A Hamiltonian curve
can be problematic when its velocity or length is too large and this happens with non-zero probability. Rather than using supremum bounds for our smoothness parameters, it suffices to use large
probability bounds, where the probability is over the random choice of Hamiltonian curve at any
point x ∈ Ω. To capture the notion that “most Hamiltonian curves are well-behaved”, we use an
auxiliary function ℓ(γ) ≥ 0 which assigns a real number to each Hamiltonian curve γ and measures
how “good” the curve is. The smoothness parameters assume that this function ℓ is bounded and
Lipshitz. One possible choice of such ℓ is ℓ(γ) = kγ ′ (0)kγ(0) which measures the initial velocity, but
this will give us a weaker bound. Instead, we use the following which jointly bounds the change in
position (first term) and change in velocity (second term).
Definition 15. An auxiliary function ℓ is a non-negative real-valued function on the set of Hamiltonian curves, i.e., maps γ : [0, δ] → M , with bounded parameters ℓ0 , ℓ1 such that
1. For any variation γs of a Hamiltonian curve (see Definition 18) with ℓ(γs ) ≤ ℓ0 , we have
!
d
d
′
ℓ(γs ) ≤ ℓ1
γs (0)
+ δ Ds γs (0) γs (0) .
ds
ds
γs (0)


1
min 1, ℓℓ10δ where γ ∼ x indicates a ran2. For any x ∈ M , Pγ∼x (ℓ(γ) ≤ 21 ℓ0 ) ≥ 1 − 100
dom Hamiltonian curve starting at x, chosen by picking a random Gaussian initial velocity
according to the local metric at x.
3.2.1

Proof Outline

To bound the conductance of HMC, we need to show that one-step distributions from nearby points
have large overlap for reasonably large step size δ. To this end, recall that the probability density
of going from x to y is given by the following formula
s


X
1
det (g(y))
2
−1
exp − kvx kx .
px (y) =
|det (DHamx,δ (vx ))|
(2π)n
2
vx :Hamx,δ (vx )=y

In Section 3.2.2, we introduce the concept of variations of Hamiltonian curves and use it to bound
|det (DHamx,δ (vx ))|−1 . We can show that px (y) is in fact close to
s


X
1
det (g(y))
1
2
·
exp − kvx kx .
(3.1)
pex (y) =
δn
(2π)n
2
vx :Hamx,δ (vx )=y

To compare px (y) with pz (y), we need to relate vx and vz that map x and z to y respectively. In
Section 3.2.3, we shows that if x and z are close enough, for every vx , there is a unique vz such that
vx is close to vz and that Hamz,δ (vz ) = Hamx,δ (vx ). Combining these facts, we obtain our main
theorem for this section, stated in Subsection 3.2.4.
15

In the analysis, we use three important operators from the tangent space to itself. The motivation for defining these operators comes directly from Lemma 19, which studies the variation in
Hamiltonian curves as the solution of a Jacobi equation. In words, the operator R(.) below allows
us to write the change in the Hamiltonian curve as an ODE.
Definition 16. Given a Hamiltonian curve γ, let R(γ, t), M (γ, t) and Φ(γ, t) be the operators from
T M to T M defined by
R(t)u = R(u, γ ′ (t))γ ′ (t),
M (t)u = Du µ(γ(t)),
Φ(t)u = M (t)u − R(t)u.
When γ is explicit from the context, we simply write them as R(t), M (t) and Φ(t).
The key parameter R1 we use in this section is a bound on the Frobenius norm of Φ formally
defined as follows.
Definition 17. Given a manifold M with metric g and an auxiliary function ℓ with parameters
ℓ0 , ℓ1 , we define the smoothness parameter R1 depending only on M and the step size δ such that
kΦ(γ, t)kF,γ(t) ≤ R1
for any γ such that ℓ(γ) ≤ ℓ0 and any 0 ≤ t ≤ δ where the Frobenius norm kAkF,γ(t) is defined by
kAk2F,γ(t) = Eα,β∼N (0,g(x)−1 ) (αT Aβ)2 .

The above definitions are related to but different from our previous paper analyzing the geodesic
walk [13].
3.2.2

Variation of Hamiltonian curve

To bound the determinant of the Jacobian of Hamx , we study variations of Hamiltonian curves.
Definition 18. We call γs (t) a Hamiltonian variation if γs (·) satisfies the Hamiltonian equation for
s
every s. We call ∂γ
∂s a Jacobi field.
The following lemma shows that a Jacobi field satisfies the following Jacobi equation.
Lemma 19. Given a path c(s), let γs (t) = Hamc(s) (t(v + sw)) be a Hamiltonian variation. The
def

Jacobi field ψ(t) =

∂
∂s γs (t)|s=0

satisfies the following Jacobi equation
Dt2 ψ(t) = Φ(t)ψ(t)

(3.2)

Let Γt parallel transport from Tγ(t,0) M to Tγ(0,0) M and ψ(t) = Γt ψ(t). Then, ψ(t) satisfies the
following ODE on the tangent space Tγ(0,0) M:
′′

ψ (t) = Γt Φ(t)Γ−1
t ψ(t)
′

ψ (0) = w,
ψ(0) = Ds c(0).

16

∀t ≥ 0,

(3.3)

Proof. Taking derivative Ds on both sides of Dt ∂γ
∂t = µ(γ), and using Fact 67, we get
∂γ
∂t
∂γ ∂γ ∂γ
∂γ
+ R( ,
)
= Dt Ds
∂t
∂s ∂t ∂t
∂γ
∂γ ∂γ ∂γ
= Dt2
+ R( ,
) .
∂s
∂s ∂t ∂t

Ds µ(γ) = Ds Dt

In short, we have Dt2 ψ(t) = Φ(t)ψ(t). This shows (3.2).
Equation (3.3) follows from the fact that
Dt v(t) = Γt


d
Γ−1
t v(t)
dt

′

for any vector field on γ0 (t) (see Definition 11 in the appendix) applied to v(t) = ψ (t).
We now proceed to estimate the determinant of the Jacobian of Hamx . For this we will use the
following elementary lemmas describing the solution of the following second-order matrix ODE:
d2
Ψ(t) = Φ(t)Ψ(t),
dt2
d
Ψ(0) = B,
dt
Ψ(0) = A.

(3.4)

Lemma 20. Consider the matrix ODE (3.4). Let λ = max0≤t≤ℓ kΦ(t)k2 . For any t ≥ 0, we have
that
√
√
kBk
kΨ(t)k2 ≤ kAk2 cosh( λt) + √ 2 sinh( λt).
λ
Lemma 21. Consider the matrix ODE (3.4). Let λ = max0≤t≤ℓ kΦ(t)kF . For any 0 ≤ t ≤
have that


t3
2
kΨ(t) − A − BtkF ≤ λ t kAk2 + kBk2 .
5

√1 ,
λ

we

In particular, this shows that
Ψ(t) = A + Bt +

with kE(s)kF ≤ λ s2 kAk2 +

s3
5


kBk2 .

Z

t
0

(t − s)Φ(s)(A + Bs + E(s))ds

The proofs of these lemmas are in Appendix A. We continue with the main proof here.

Lemma 22. Let γ(t) = Hamx (tvx ) be a Hamiltonian curve and step size δ satisfy 0 < δ2 ≤ R11
where R1 = max0≤t≤h kΦ(t)kF,γ(t) . Then DHamx,δ is invertible with kDHamx,δ − δIkF,γ(δ) ≤ 5δ .
Also, we have,
log det



2
 Z δ
δ 2 R1
t(δ − t)
1
DHamx,δ (vx ) −
TrΦ(t)dt ≤
.
δ
δ
10
0
17

(3.5)

Proof. We want to compute DHamx,δ (vx )[w] for some w ∈ Tx M . By definition, we have that
DHamx,δ (vx )[w] =

∂
γ(t, s)|t=δ,s=0
∂s

(3.6)

where γ(t, s) = Hamx (t(vx + sw)). Define ψ(t) as in Lemma 19 with c(s) = x. So, Ds c(0) = 0, i.e.,
ψ(0) = 0. Then, by the lemma,
′′

ψ (t) = Γt Φ(t)Γ−1
t ψ(t),

′

ψ (0) = w,

ψ(0) = 0.

Now, we define Ψ be the solution of the matrix ODE
Ψ′′ (t) = Γt Φ(t)Γ−1
t Ψ(t) ∀t ≥ 0,
Ψ′ (0) = I,
Ψ(0) = 0.
By the definition of Ψ, we see that ψ(t) = Ψ(t)w. Therefore, we have that
∂
−1
γ(t, s)|s=0 = Γ−1
t ψ(t) = Γt Ψ(t)w.
∂s
Combining it with (3.6), we have that DHamx,δ (vx ) = Γ−1
t Ψ(δ). Since Γt is an orthonormal matrix,
we have that
log det (DHamx,δ (vx )) = log det Ψ(δ).
(3.7)
Note that Γt Φ(t)Γ−1
t
that

F,γ(t)

= kΦ(t)kF,γ(t) ≤ R1 for all 0 ≤ t ≤ δ. Using this, Lemma 21 shows
1
Ψ(δ) − I
δ

F,x

≤ R1



δ2
kIk2
5



≤

1
5

(3.8)

Hence, Ψ(δ) is invertible, and so is DHamx .
By Lemma 64, we have that
1
log det( Ψ(δ)) − Tr
δ



1
Ψ(δ) − I
δ



≤



1 2
δ R1
5

2

.

(3.9)

Now we need to estimate Tr(Ψ(δ) − δI). Lemma 21 shows that
Z δ
Z δ
(δ − t)Φ(t)E(t)dt
t(δ − t)Φ(t)dt +
Ψ(δ) = δI +
0

0

3

with kE(t)kF,γ(t) ≤ t 5R1 for all 0 ≤ t ≤ δ. Since

Z δ
Z δ
δ5
(δ − t) kΦ(t)kF,x kE(t)kF,x dt ≤ R12 ,
(δ − t)Φ(t)E(t)dt ≤
Tr
20
0
0
we have that
Tr



1
Ψ(δ) − I −
δ

Z

δ
0

t(δ − t)
Φ(t)dt
δ



≤

δ4 2
R .
20 1

Combining (3.9) and (3.10), we have
2

2
Z δ
δ 2 R1
t(δ − t)
1 2
δ4 2
1
TrΦ(t)dt ≤
δ R1 + R1 ≤
.
log det( Ψ(δ)) −
δ
δ
5
20
10
0
Applying (3.7), we have the result.
18

(3.10)

3.2.3

Local Uniqueness of Hamiltonian Curves

Next, we study the local uniqueness of Hamiltonian curves. We know that for every pair x, y, there
can be multiple Hamiltonian curves connecting x and y. Due to this, the probability density px
at y in M is the sum over all possible Hamiltonian curves connecting x and y. The next lemma
establishes a 1-1 map between Hamiltonian curves connecting x to y as we vary x.
Lemma 23. Let γ(t) = Hamx (tvx ) be a Hamiltonian curve and let the step size δ satisfy 0 < δ2 ≤
1
R1 , where R1 = max0≤t≤δ kΦ(t)kF,γ(t) . Let the end points be x = γ(0) and y = γ(δ). Then there is
an unique smooth invertible function v : U ⊂ M → V ⊂ T M such that
y = Hamz,δ (v(z))
for any z ∈ U where U is a neighborhood of x and V is a neighborhood of vx = v(x). Furthermore,
5
we have that k∇η v(x)kx ≤ 2δ
kηkx and
1
η + ∇η v(x)
δ

x

≤

3
R1 δ kηkx .
2

Let γs (t) = Hamc(s) (t · v(c(s))) where c(s) is any path with c(0) = x and c′ (0) = η. Then, for
all 0 ≤ t ≤ δ, we have that
∂
∂s

γs (t)
s=0

γ(t)

≤ 5 kηkx

and

Ds γs′ (t)

s=0 γ(t)

≤

10
kηkx .
δ

Proof. Consider the smooth function f (z, w) = Hamz,δ (w). From Lemma 22, the Jacobian of w
at (x, vx ) in the w variables, i.e., DHamx,δ (vx ), is invertible. Hence, the implicit function theorem
shows that there is a open neighborhood U of x and a unique function v on U such that f (z, v(z)) =
f (x, vx ), i.e. Hamz,δ (v(z)) = Hamx,δ (vx ) = y.
To bound ∇η v(x), we let γs (t) = Hamc(s) (t · v(c(s))) and c(s) be any path with c(0) = x and
′
c (0) = η. Let Γt be the parallel transport from Tγ(t) M to Tγ(0) M. Define
ψ(t) = Γt

∂
∂s

γs (t).
s=0

Lemma 19 shows that ψ(t) satisfies the following ODE
′′

ψ (t) = Γt Φ(t)Γ−1
t ψ(t) ∀t ≥ 0,
′

ψ (0) = ∇η v(x),
ψ(0) = η.

Moreover, we know that ψ(δ) = 0 because γs (δ) = Hamc(s) (δ · v(c(s))) = Hamc(s),δ (v(c(s))) = y for
small enough s.
To bound k∇η v(x)kx , we note that Lemma 21 shows that


t
2
(3.11)
ψ(t) − η − t · ∇η v(x) x ≤ R1 t kηkx + k∇η v(x)kx .
5
Since ψ(δ) = 0 and δ2 ≤

1
R1 ,

we have that
kη + δ · ∇η v(x)kx ≤ kηkx +
19

δ
k∇η v(x)kx
5

which implies that
δ k∇η v(x)kx − kηkx ≤ kηkx +
Therefore, k∇η v(x)kx ≤

5
2δ

δ
k∇η v(x)kx .
5

kηkx . More precisely, from (3.11), we have that
1
η + ∇η v(x)
δ

x

≤

3
R1 δ kηkx .
2

Putting this into (3.11), for t ≤ δ we get
ψ(t)

x



5
1
2
≤ kηkx + kηkx + R1 t kηkx + kηkx ≤ 5 kηkx .
2
2

Now, apply the conclusion of Lemma 21 after taking a derivative, we have that
Z t
′
Φ(s)(η + s · ∇η v(x) + E(s))ds
ψ (t) = ∇η v(x) +
0



where kE(s)kx ≤ R1 s2 kηkx +
that t ≤ δ, we have
′

ψ (t)

x

s3
5

k∇η v(x)kx



≤

3
2

kηkx . Hence, bounding each term and noting



5
10
5 3
≤
kηkx + δR1 1 + +
kηkx .
kηkx ≤
2δ
2 2
δ

When we vary x, the Hamiltonian curve γ from x to y varies and we need to bound ℓ(γ) over
the variation.
Lemma 24. Given a Hamiltonian curve γ(t) = Hamx (t · vx ) with step size δ satisfying δ2 ≤ R11 , let
c(s) be any geodesic starting at γ(0). Let x = c(0) = γ(0) and y = γ(δ). Suppose that the auxiliary
ℓ0
dc
function ℓ satisfies ds
≤ 7ℓ
and ℓ(γ) ≤ 21 ℓ0 . Then, there is a unique vector field v on c such
c(0)
1
that
y = Hamc(s),δ (v(s)).
Moreover, this vector field is uniquely determined by the geodesic c(s) and any v(s) on this vector
field. Also, we have that ℓ(Hamc(s),δ (v(s))) ≤ ℓ0 for all s ≤ 1.
Proof. Let smax be the supremum of s such that v(s) can be defined continuously such that y =
Hamc(s),δ (v(s)) and ℓ(γs ) ≤ ℓ0 where γs (t) = Hamc(s) (t · v(s)). Lemma 23 shows that there is a
neighborhood N at x and a vector field u on N such that for any z ∈ N , we have that
y = Hamz,δ (u(z)).
Also, this lemma shows that u(s) is smooth and hence the parameter ℓ1 shows that ℓ(γs ) is Lipschitz
in s. Therefore, ℓ(γs ) ≤ ℓ0 in a small neighborhood of 0. Hence smax > 0.
Now, we show smax > 1 by contradiction. By the definition of smax , we have that ℓ(γs ) ≤ ℓ0 for
dc
5
5
any 0 ≤ s < smax . Hence, we can apply Lemma 23 to show that kDs v(s)kγ(s) ≤ 2δ
ds γ(s) = 2δ L
where L is the length of c up to s = 1 (since the speed is constant on any geodesic and the curve
is defined over [0, 1]). Therefore, the function v is Lipschitz and hence v(smax ) is well-defined and
ℓ(γsmax ) ≤ ℓ0 by continuity. Hence, we can apply Lemma 23 at ℓ(smax ) and extend the domain of
ℓ(s) beyond smax .
20

To bound ℓ(γs ) beyond smax , we note that kDs γs′ kγ(s) = kDs v(s)kγ(s) ≤

d
ds c γ(0)

= L. Hence,

d
ds ℓ(γs )

5
2δ L

and

≤ (L + 25 L)ℓ1 by the definition of ℓ1 . Therefore, if L ≤

that ℓ(γs ) ≤ ℓ(γ) + 21 ℓ0 ≤ ℓ0 for all s ≤ 1.01 wherever v(s) is
the assumption that smax is the supremum. Hence, smax > 1.

d
ds c γ(s)
ℓ0
7ℓ1 ,

=

we have

defined. Therefore, this contradicts

The uniqueness follows from Lemma 23.
3.2.4

Smoothness of one-step distributions

Lemma 25. For δ2 ≤
px , pz from x, z satisfy

1
√
100 nR1

and δ3 ≤

√ℓ0
,
100 nR1 ℓ1

the one-step Hamiltonian walk distributions

 
1
1
dTV (px , pz ) = O
d(x, z) + .
δ
25

ℓ0
. Let c(s) be a unit speed geodesic connecting x and
Proof. We first consider the case d(x, y) < 7ℓ
1
ℓ0
z of length L < 7ℓ1 .


def
Let ℓe = min 1, ℓℓ10δ . By the definition of ℓ0 , with probability at least 1 − ℓe in paths γ start at

x, we have that ℓ(γ) ≤ 21 ℓ0 . Let Vx be the set of vx such that ℓ(Hamx (t · vx )) ≤ 21 ℓ0 . Since the
ℓ0
and ℓ(γ) ≤ 12 ℓ0 , for those vx , Lemma 24 shows there is a family
distance from x to z is less than 7ℓ
1
of Hamiltonian curves γs that connect c(s) to y, and ℓ(γs ) ≤ ℓ0 for each of them.
For any vx ∈ V , we have that ℓ(γs ) ≤ ℓ0 . When ℓ(γs ) ≤ ℓ0 , by the definition of R1 , we indeed
have that kΦ(t)k ≤ R1 and hence Lemma 22 shows that
2


Z δ
δ 2 R1
t(δ − t)
1
DHamx,δ (vx ) ≤
TrΦ(t)dt +
log det
δ
δ
10
0
2
2
2
δ R1
1 e
δ √
nR1 +
≤
ℓ
≤
6
10
600

where we used our assumption on δ. We use p(vx ) to denote the probability density of choosing vx
and
s


det (g(Hamx (δ · vx )))
1
def
2
exp − kvx kx .
pe(vx ) =
(2πδ2 )n
2
Hence, we have that

C −1 · p(vx ) ≤ pe(vx ) ≤ C · p(vx )

(3.12)

1 e
ℓ. As we noted, for every vx , there is a corresponding vz such that Hamz (δ·vz ) = y
where C = 1+ 600
and ℓ(Hamz (t · vz )) ≤ ℓ0 . Therefore, we have that
def

(1 − C 2 )p(vx ) + C(e
p(vx ) − pe(vz )) ≤ p(vx ) − p(vz ) ≤(1 − C −2 )p(vx ) + C −1 (e
p(vx ) − pe(vz )).


R
1
Since Vx p(vx )dvx ≥ 1 − 100
min 1, ℓℓ10δ , we have that
Z
ℓe
|p(vx ) − p(vz )| dvx
+
100
Vx
Z
Z
ℓe
≤
|e
p(vx ) − pe(vz )| dvx
p(vx )dvx ) + 2
(1 +
100
Vx
Vx
Z Z
d
ℓe
+2
pe(vc(s) ) dsdvx
≤
50
Vx s ds

dT V (px , pz ) ≤

21

(3.13)

Note that

d
pe(v ) =
ds c(s)



−


1 d
kv(s)k2c(s) pe(vc(s) ).
2 ds

Using (3.12), we have that pe(vc(s) ) ≤ 2 · p(vc(s) ) and hence
Z

Vx

d
pe(v ) dy ≤
ds c(s)

Z

d
kv(s)k2c(s) p(vc(s) )dy
Vx ds
d
kv(s)k2c(s)
≤Eℓ(γs )≤ℓ0
ds

(3.14)

Using that k∂s c(s)k = 1, Lemma 23 shows that
1
∂s c(s) + Ds v(s)
δ

x

3
≤ R1 δ.
2

Therefore, we have that
d
kv(s)k2c(s) = 2 hv(s), Ds v(s)ic(s)
ds
2
hv(s), ∂s c(s)ic(s) + 3R1 δ kv(s)kc(s) .
≤
δ
Since v(s) is a random Gaussian vector from the local metric, we have that hv(s), ∂s c(s)ic(s) = O(1)
√
and kv(s)kc(s) = O( n) with high probability. Putting it into (3.14), we have that
Eℓ(γs )≤ℓ0

d
kv(s)k2c(s) = O
ds




√
1
+ δR1 n .
δ

Putting this into 3.14, we have that


Z
√
d
1
pe(vc(s) ) dy = O
+ δR1 n .
δ
Vx ds
Putting this into (3.13), we get

dTV (px , pz ) = O




√
ℓe
1
+ δR1 n L +
δ
50

ℓ0
for any L < 7ℓ
. By taking a minimal length geodesic, and summing over segment of length
1
any x and z, we have


√
1
1
+ δR1 n d(x, z) +
dTV (px , pz ) = O
δ
25
 
1
1
=O
d(x, z) + .
δ
25

22

ℓ0
8ℓ1 ,

for

3.3

Convergence bound

Combining Lemma 25 and Lemma 14, we have the following result.
Theorem 26. Given a manifold M. Let ℓ0 , ℓ1 , R1 be the parameters of the Hamiltonian Monte
Carlo defined in Definition 15 and 17. Let qt be the distribution of the current point after t steps
Hamiltonian Monte Carlo with step size δ satisfying
δ2 ≤

1
√
100 nR1

δ3 ≤

and

ℓ
√0
,
100 nR1 ℓ1

starting from initial distribution q0 . Let q be the distribution proportional to e−f . For any ε > 0,
we have that
s

t
dq0 (x)
(δψ)2
1
Ex∼q0
1−
dT V (qt , q) ≤ ε +
ε
dq(x)
2
where ψ is the conductance of the manifold defined in Definition 12.
For ℓ(γ) = kγ ′ (0)kγ(0) , we can bound ℓ0 and ℓ1 as follows:

√
Lemma 27. For the auxiliary function ℓ(γ) = kγ ′ (0)kγ(0) , we have that ℓ0 = 10 n and ℓ1 = O( 1δ ).
Furthermore, we have that
γ ′ (t) γ(t) ≤ γ ′ (0) γ(0) + R0 t
where R0 = supx∈M kµ(x)kx .

√
Proof. For ℓ0 , we note that γ ′ (0) ∼ N (0, g(x)−1 ) and hence kγ ′ (0)kγ(0) ≤ 5 n with probability
√

e−O( n) .
For ℓ1 , we note that

d
d
γ ′ (0)
ℓ(γs ) =
ds
ds s
=

γs (0)

≤

1
2

hDs γs′ (0), γs′ (0)iγs (0)

d
ds

kγs′ (0)k2γs (0)

kγs′ (0)kγs (0)

kγs′ (0)kγs (0)

≤ Ds γs′ (0)

γs (0)

.

Hence, we have that ℓ1 = 1δ .
Next, we note that
d
γ ′ (t)
dt

2
γ(t)

= 2 Dt γ ′ (t), γ ′ (t)

γ(t)

≤ 2 kµ(γ(t))kγ(t) γ ′ (t)

γ(t)

.

Therefore, we get the last result.
Combining Lemma 27 and Theorem 26, we have the following result. This result might be
more convenient to establish an upper bound on the rate of convergence, as it depends on only two
worst-case smoothness parameters. In the next section, we will see a more refined bound that uses
the randomness of Hamiltonian curves via additional parameters.

23

Corollary 28. Given a manifold M. Let qt be the distribution of the current point after t steps
Hamiltonian Monte Carlo with step size δ and q be the distribution proportional to e−f . Let R0 and
R1 be parameters such that
1. kµ(x)kx ≤ R0 for any x ∈ M where µ is defined in Lemma 7.
2. Eα,β∼N (0,g(x)−1 ) hDα µ(x), βi2x ≤ R12 for any x ∈ M.
3. Eα,β∼N (0,g(x)−1 ) hR(α, v)v, βi 2x ≤ R12 for any x ∈ M and any kvkx ≤
Riemann curvature tensor of M.
Suppose that δ ≤

√

n
R0

and δ2 ≤

1
√
,
100 nR1

√

n where R is the

then for any ε > 0, we have that

dT V (qt , q) ≤ ε +

s

1
dq0 (x)
Ex∼q0
ε
dq(x)



1−

(δψ)2
2

t

where ψ is the conductance of the manifold defined in Definition 12. In short, the mixing time of
Hamiltonian Monte Carlo is



2
√
R
−2
0
e ψ
.
nR1 +
O
n
Proof. The statement is basically restating the√ definition of R1 and R0 used in Lemma 27 and
Theorem 26. The only difference is that if δ ≤ Rn0 , then we know that
γ ′ (t)

γ(t)

≤ γ ′ (0)

γ(0)

√
+ R0 t = O( n)

for all 0 ≤ t ≤ δ. Therefore, we can relax the constraints ℓ(γ) ≤ ℓ0 in the definition of R1 to simply
√
kvkx ≤ n. It allows us to use R1 without mentioning the auxiliary function ℓ.

4

Improved analysis of the convergence

Corollary 28 gives a polynomial mixing time for the log barrier function. There are two bottlenecks
to improving the bound. First, the auxiliary function ℓ it used does not capture the fact each curve
γ in Hamiltonian Monte Carlo follows a random initial direction. Second, Lemma 23 also does not
take full advantage of the random initial direction. In this section, we focus on improving Lemma
23.
Our main theorem for convergence can be stated as follows in terms of ψ and additional parameters ℓ0 , ℓ1 , R1 , R2 , R3 (see Definitions 12, 15, 16, 17, 33 and 31). It uses the following key
lemma.
Theorem 29. For δ2 ≤
x, z satisfy

1
R1

and δ5 ≤

ℓ0
,
R21 ℓ1

the one-step Hamiltonian walk distributions px , pz from


1
1
dTV (px , pz ) = O δ R2 + + δR3 d(x, z) + .
δ
25


2

Remark. The constant term at the end is an artifact that comes from bounding the probability of
some bad events of the Hamiltonian walk, and can be made to arbitrary small.

24

We now prove this theorem. In a later section, we specialize to sampling distributions over
polytopes using the logarithmic barrier, by defining a suitable auxiliary function and bounding all
the parameters. The key ingredient is Theorem 29 about the overlap of one-step distributions, which
we prove in the next section.
Here is the consequence of Theorem 29 and Theorem 13.
Theorem 30. Given a manifold M. Let ℓ0 , ℓ1 , R1 , R2 , R3 be the parameters of the Hamiltonian
Monte Carlo defined in Definition 15, 17, 33 and 31. Let qt be the distribution of the current point
after t steps Hamiltonian Monte Carlo starting from initial distribution q0 . Let q be the distribution
proportional to e−f . Suppose that the step size δ satisfies
δ2 ≤

1
ℓ0
, δ5 ≤ 2 and δ3 R2 + δ2 R3 ≤ 1.
R1
R1 ℓ 1

For any ε > 0, we have that
s

dT V (qt , q) ≤ ε +

dq0 (x)
1
Ex∼q0
ε
dq(x)



(δψ)2
1−
2

t

where ψ is the conductance of the manifold.

4.1

Improved one-to-one correspondence for Hamiltonian curve

In the previous section, we only used R1 to analyze how much a Hamiltonian curves change as one
end point varies. Here we derive a more refined analysis of Lemma 23 using an additional parameter
R3 .
Definition 31. For a manifold M and auxiliary function ℓ, R3 is a constant such that for any
Hamiltonian curve γ(t) of step size δ with ℓ(γ0 ) ≤ ℓ0 , if ζ(t) is the parallel transport of the vector
γ ′ (0) along γ(t), then we have
sup kΦ(t)ζ(t)kγ(t) ≤ R3 .
0≤t≤δ

Lemma 32. Under the same assumptions as Lemma 23, we have that
δ
∇η kv(x)k2x ≤ |hvx , ηix | + 3δ2 R3 kηkx .
2
Proof. Let χ = ∇η v(x). Define ψ(t) as in the proof of Lemma 23. Using Lemma 21 we get that
0 = ψ(δ) = η + δχ +

Z

0

i.e.,
−δχ = η +
with

Z

δ
0

kE(s)kx ≤ R1

δ

(δ − s)Γs Φ(s)Γ−1
s (η + sχ + e(s))ds

(δ − s)Γs Φ(s)Γ−1
s (η + sχ + E(s))ds


δ3
kχkx
δ kηkx +
5
2

for all 0 ≤ s ≤ δ where we used k∇η v(x)kx ≤

5
2δ



≤ 2R1 δ2 kηkx

kηkx at the end (by Lemma 23).
25

Therefore,
δ |hv(x), χix | ≤ |hv(x), ηix | +

δ2
sup hv(x), Γs Φ(s)Γ−1
s (η + sχ + E(s))ix .
2 0≤s≤δ

Noting that
kη + sχ + E(s)kx ≤ kηkx + δ kχkx + kE(s)kx ≤ (1 +

5
+ 2R1 δ2 ) kηkx ≤ 6 kηkx ,
2

δ |hv(x), χix | ≤ |hv(x), ηix | + 3δ2 sup

x

we have
0≤s≤δ

= |hv(x), ηix | + 3δ2 sup

0≤s≤δ

γ ′ (0)T Γs Φ(s)Γ−1
s
Φ(s)Γs γ ′ (0)

γ(s)

kηkx

kηkx

≤ |hv(x), ηix | + 3δ2 R3 kηkx .
Finally, we note that

4.2

δ
∇η kv(x)k2x = δ hv(x), ∇η v(x)ix .
2

Improved smoothness of one-step distributions

The proof of Theorem 29 is pretty similar to Lemma 25. First, we show that px (y) is in fact close
to
s
 Z δ

X
t(δ − t)
1
det (g(y))
2
exp −
TrΦ(t)dt − kvx kx .
(4.1)
pex (y) =
(2πδ2 )n
δ
2
0
vx :Hamx,δ (vx )=y

Note that this is a more refined estimate than (3.1). We use Lemma 32 to bound the change of
kvx k2x . For the change of TrΦ(t), we defer the calculation until the end of this section.

Theorem 29. For δ2 ≤
x, z satisfy

1
R1

and δ5 ≤

ℓ0
,
R21 ℓ1

the one-step Hamiltonian walk distributions px , pz from


1
1
dTV (px , pz ) = O δ R2 + + δR3 d(x, z) + .
δ
25


Proof. We first consider the case d(x, y) <
δ2 ≤ R11 and δ5 ≤ Rℓ20ℓ , we have that

2

ℓ0
7ℓ1 .

By a similar argument as Lemma 25, using that

1 1

Z Z
d
ℓe
+2
pe(vc(s) ) dsdvx
dT V (px , pz ) ≤
50
Vx s ds



def
where ℓe = min 1, ℓℓ10δ . By direct calculation, we have
d
pe(v ) =
ds c(s)

 Z
−

δ

0


t(δ − t) d
1 d
2
′
TrΦ(γs (t))dt −
kv(s)kc(s) pe(vc(s) ).
δ
ds
2 ds

By similar argument as (3.12), we have that pe(vc(s) ) ≤ 2 · p(vc(s) ) and hence

Z δ
t(δ − t) d
1 d
d
2
′
pe(v ) ≤ 2
TrΦ(γs (t))dt +
kv(s)kc(s) p(vc(s) ).
ds c(s)
δ
ds
2 ds
0
26

(4.2)

Since ℓ(γs ) ≤ ℓ0 , we can use Lemma 34 to get
Z

δ

0

Hence,


t(δ − t) d
TrΦ(γs′ (t))dt ≤ O δ2 R2 .
δ
ds

Z

Z
Z
d
2
pe(vc(s) ) dvx ≤O(δ R2 )
p(vc(s) )dvx +
Vx ds
Vx
Vx
R
For the first term, we note that Vx p(vc(s) )dvx ≤ 1.
For the second term, we have that
Z
d
kv(s)k2c(s) p(vc(s) )dvx ≤Eℓ(γs )≤ℓ0
Vx ds

d
kv(s)k2c(s) p(vc(s) )dvx .
ds

(4.3)

d
kv(s)k2c(s) .
ds

(4.4)

By Lemma 32, we have that
δ d
kv(s)k2c(s) ≤ |hv(s), ∂s cix | + 3δ2 R3 k∂s ckx .
2 ds
Since k∂s ckx = 1 and |hv(s), ∂s cix | = O(1) with high probability since v(s) is a Gaussian vector
from the local metric. We have
δ d
kv(s)k2c(s) ≤ O(1) + 3δ2 R3 .
2 ds
Putting it into (4.4) and (4.3), we have that


Z
d
1
2
pe(vc(s) ) dvx = O δ R2 + + δR3 .
δ
Vx ds
Putting this into 4.2, we get




1
ℓe
dTV (px , pz ) = O δ R2 + + δR3 L +
δ
50
2

ℓ0
. By taking a minimal length geodesic, and summing over segment of length
for any L < 7ℓ
1
any x and z, we have


1
1
2
dTV (px , pz ) = O δ R2 + + δR3 d(x, z) + .
δ
25

ℓ0
8ℓ1 ,

for

Definition 33. Given a Hamiltonian curve γ(t) with ℓ(γ0 ) ≤ ℓ0 . Let R2 be a constant depending
on the manifold M and the step size δ such that for any 0 ≤ t ≤ δ, any curve c(s) starting from
γ(t) and any vector field v(s) on c(s) with v(0) = γ ′ (t), we have that
!
dc
d
+ δ k Ds v|s=0 kγ(t) R2 .
TrΦ(v(s))|s=0 ≤
ds
ds s=0 γ(t)

27

Lemma 34. For δ2 ≤

1
R1

and ℓ(γs ) ≤ ℓ0 , we have
Z

δ
0


t(δ − t) d
TrΦ(γs′ (t))dt ≤ O δ2 R2
δ
ds

where γs is a family of Hamiltonian curve that connect c(s) to y defined in Lemma 24.
Proof. By Definition 33, we have that
d
TrΦ(γs′ (t)) ≤
ds

∂
∂s

γs (t)
s=0

+δ
γ(t)

Ds γs′ (t) s=0 γ(t)

!

· R2 .

(4.5)

d
d
d
γs (0) = ds
Hamc(s) (0) = ds
c(s) is a unit vector and hence
By definition of γs , we have that ds
d
ds γs (0) γ(0) = 1. Since Hamc(s),δ (v(s)) = y, Lemma 23 shows that

∂
∂s

γs (t)
s=0

γ(t)

≤ 5 and

Ds γs′ (t)

Therefore, we have that
d
TrΦ(γs′ (t)) ≤ 15R2 .
ds

28

s=0 γ(t)

≤

10
.
δ

5
5.1

Gibbs sampling on manifolds
Isoperimetry for Hessian manifolds

Here we derive a general isoperimetry bound, assuming that the manifold is defined by the Hessian
of a convex function, and that the directional fourth directive is non-negative, a property satisfied,
e.g., by the standard logarithmic barrier.
Lemma 35. Let φ : [a, b] → R be a convex function such that φ′′ is also convex. For any x ∈ [a, b],
we have that

Z b
Z x
p
−φ(t)
−φ(t)
−φ(x)
′′
e
dt .
e
dt,
e
≥ 0.372 φ (x) min
a

x

Let f (x) be the logconcave density proportional to
VarX∼f (x) X ≤

e−φ(x) .

Then, we have that

O(1)
.
minx∈[a,b] φ′′ (x)

p
Proof. Case 1)p|φ′ (x)| ≥ a φ′′ (x) (we will pick a at the end). Without loss of generality, we have
that φ′ (x) ≥ a φ′′ (x). Since φ is convex, we have that
φ(t) ≥ φ(x) + φ′ (x)(t − x).

Therefore, we have that
Z

Case 2) |φ′ (x)| ≤ a
we have that

p

b
x

e−φ(t) dt ≤

Z

∞

′

e−φ(x)−φ (x)(t−x) dt

x

e−φ(x)
e−φ(x)
= ′
≤ p
.
φ (x)
a φ′′ (x)

φ′′ (x). Without loss of generality, we have that φ′′′ (x) ≥ 0. Since φ′′ is convex,
φ′′ (t) ≥ φ′′ (x) + φ′′′ (x)(t − x) ≥ φ′′ (x)

for t ≥ x. Hence, we have that
′

φ(t) = φ(x) + φ (x)(t − x) +

Z

t
x

(t − s)φ′′ (t)ds

1
≥ φ(x) + φ (x)(t − x) + φ′′ (x)(t − x)2
2
′

for all t ≥ x. Therefore, we have that
Z
Z b
−φ(t)
e
dt ≤
x

∞

′

1

′′

2

e−φ(x)−φ (x)(t−x)− 2 φ (x)(t−x) ds
−∞
Z ∞
1 ′′
2
′
= e−φ(x)
e−φ (x)s− 2 φ (x)s ds
−∞
s
2
2π −φ(x)+ 21 φφ′′′(x)
(x)
=
e
φ′′ (x)
s
2π −φ(x)+ a2
2 .
e
≤
′′
φ (x)
29

Hence, we have

Z

−φ(x)
√
a2 e
.
2πe 2 p
φ′′ (x)
x

2 
− a2
e√
Combining both cases, the isoperimetric ratio is min a, 2π . Setting a to be the solution of
p
a2
ae 2 = √12π , this minimum is achieved at W (1/2π) > 0.372 where W is the inverse Lambert
b

e−φ(t) dt ≤

1
. This proves the first result.
function, i.e., W (x)eW (x) = 2π
The variance of f follows from the fact that f is logconcave.

This generalizes to higher dimension with no dependence on the dimension using localization,
which we review next. Define an exponential needle E = (a, b, γ) as a segment [a, b] ⊆ Rn and γ ∈ R
corresponding to the weight function eγt applied the segment [a, b]. The integral of an n-dimensional
function h : Rn → R over this one dimensional needle is
Z |b−a|
Z
b−a
.
h(a + tu)eγt dt
where u =
h=
|b − a|
0
E
Theorem 36 (Theorem 2.7 in [9]). Let f1 , f2 , f3 , f3 be four nonnegative continuous functions defined
on Rn , and c1 , c2 > 0. Then, the following are equivalent:
1. For every logconcave function F defined on Rn with compact support,
 c 1 Z
 c 2 Z
 c 1 Z
Z
F (x)f3 (x) dt
≤
F (x)f2 (x) dx
F (x)f1 (x) dx

Rn

Rn

Rn

Rn

F (x)f4 (x) dx

 c2

2. For every exponential needle E,
c 2
c 1 Z
c 2 Z
c 1 Z
Z
f4
f3
≤
f2
f1
E

E

E

E

Lemma 37. Let φ : K ⊂ Rn → R be a convex function defined over a convex body K such
that D 4 φ(x)[h, h, h, h] ≥ 0 for all x ∈ K and h ∈ Rn . Given any partition S1 , S2 , S3 of K with
d = minx∈S1 ,y∈S2 d(x, y), i.e., the minimum distance between S1 and S2 in the Riemannian metric
induced by φ. For any α > 0, we have that
R
−αφ(x) dx
√
S3 e
o = Ω( α · d)
nR
R
min S1 e−αφ(x) dx, S2 e−αφ(x) dx

Proof. By rescaling φ, we can assume α = 1. We write the desired inequality as follows, for a
constant C, with χS being the indicator of set S:
Z
Z
Z
Z
−φ(x)
−φ(x)
−φ(x)
e−φ(x) χS3 (x) dx.
e
dx ·
e
χS2 (x) dx ≤
e
χS1 (x) dx ·
Cd
Rn

Rn

Rn

Rn

Using the localization lemma for exponential needles (Theorem 36), with fi (x) being Cde−φ(x) χS1 (x),
e−φ(x) χS2 (x), e−φ(x) and e−φ(x) χS3 (x) with respectively, it suffices to prove the following onedimensional inequality for functions φ defined on an interval and shifted by a linear term:
Z 1
Z 1
e−φ((1−t)a+tb) e−ct χS2 ((1 − t)a + tb) dt
e−φ((1−t)a+tb) e−ct χS1 ((1 − t)a + tb) dt ·
Cd
0
0
Z 1
Z 1
e−φ((1−t)a+tb) e−ct χS3 ((1 − t)a + tb) dt.
e−φ((1−t)a+tb) e−ct dt ·
≤
0

0

30

Each Ti = {t : (1 − t)a + tb ∈ Si } is a union of intervals. By a standard argument (see [19]),
it suffices to consider the case when each Si is a single interval and add up over all intervals in
S3 . Thus it suffices to prove the statement in one dimension for all convex φ with convex φ′′ . In
one-dimension, we have
Z yp
d(x, y) =
φ′′ (t) dt.
x

Taking T3 =
is

[a′ , b′ ]

⊂ [a, b], the inequality we need to prove is that for any convex φ with convex φ′′
R a′ −φ(t) R b −φ(t)
−φ(t) dt
e
dt b′ e
a′ e
≥ Ω(1) a
Rb
R b′ p
−φ(t) dt
φ′′ (t) dt
a e
a′
R b′ −φ(t)
e
dt
e−φ(x)
≥ minx∈[a′ ,b′ ] √
and
that R ba′′ √ ′′
φ′′ (x)
φ
(t)
dt
a′
R b′

which is implied by noting

5.2

dt

applying Lemma 35.

Sampling with the log barrier

For any polytope M = {Ax > b}, the logarithmic barrier function φ(x) is defined as
φ(x) = −

m
X
i=1

log(aTi x − bi ).

Theorem 1. Let φ be the logarithmic barrier for a polytope M with m constraints and n variables.
Hamiltonian Monte Carlo applied to the function f = exp(−αφ(x)) and the metric given by ∇2 φ
with appropriate step size mixes in
!
1
2
1
3
3n3
1
1
n
m
e
+ m2 n6
O
+
α + m−1 α 31 + m− 13
steps where each step is the solution of a Hamiltonian ODE.

√
Proof. Lemma 37 shows that the isoperimetric coefficient ψ is Ω( α). Also, we know that isoperi1
metric coefficient ψ is at worst Ω(m− 2 ) [16]. Lemma 61 shows that the condition of Theorem 30 is
satisfied, thus implying the bound


√
e max( α, m− 12 )−2 δ−2
O



√
e max( α, m− 12 )−2 max n 32 , α 23 m 13 n 31 , αm 12 n 61
=O
!
2
1
1
1
1
2
3 + α 3 m 3 n 3 + αm 2 n 6
n
e
=O
α + m−1
!
1
2
1
3
3n3
1
1
n
m
e
+ m2 n6
=O
+
α + m−1 α 13 + m− 13

To implement the walk, we solve this ODE using the collocation method as described in [13].
The key lemma used in that paper is that the slack of the geodesic does not change by more than
1
). Similarly in Lemma 53, we proved
a constant multiplicative factor up to the step size δ = O( n1/4
the the Hamiltonian flow does not change by more than a constant multiplicative factor up to
−1/3
n−1/4
√ ). Since the step size we use is Θ( n √ ), we can apply the
the step size δ = O( 11/4 ) = O( 1+
α
1+ α
M1

31

collocation method as described in [13] and obtain an algorithm to compute the Hamiltonian flow
ω−1 logO(1) 1 ) time with additive error η distance in the local norm. Due to the exponential
e
in O(mn
η
convergence, it suffices for the sampling purpose with only a polylogarithmic overhead in the total
running time.

6

Polytope volume computation: Gaussian cooling on manifolds

The volume algorithm is essentially the Gaussian cooling algorithm introduced in [3]. Here we apply
it to a sequence of Gibbs distributions rather than a sequence of Gaussians. More precisely, for a
convex body K and a convex barrier function φ : K → R, we define
(

exp −σ −2 φ(x)
if x ∈ K
2
f (σ , x) =
0
otherwise
and
2

F (σ ) =
x∗

Z

f (σ 2 , x) dx
Rn

where
is the minimizer of φ (the center of K). Let µi be the probability distribution proportional
to f (σi2 , x) where σi is the tempature of the Gibbs distribution to be fixed. The algorithm estimates
each ratio in the following telescoping product:
e−σ

−2 φ(x∗ )

vol(K) ≈ F (σk ) = F (σ0 )

k
2 )
Y
F (σi+1
i=1

F (σi2 )

for some large enough σk .
2 , x)/f (σ 2 , x). Then,
Let x be a random sample point from µi and let Yx = f (σi+1
i
Ex∼µi (Yx ) =

6.1

Algorithm: cooling schedule

6.2

Correctness of the algorithm

2 )
F (σi+1
.
F (σi2 )

In this subsection, we prove the correctness of the algorithm. We separate the proof into two parts.
In the first part, we estimate how small σ0 we start with should be and how large σk we end with
should be. In the second part, we estimate the variance of the estimator Yx .
6.2.1

Initial and terminal conditions

First, we need the following lemmas about self-concordance functions and logconcave functions:
Lemma 38 ([27, Thm 2.1.1]). Let φ be a self-concordant function and x∗ be its minimizer. For any
x such that φ(x) ≤ φ(x∗ ) + r 2 with r ≤ 12 , we have that
φ(x) = φ(x∗ ) +

1 ± Θ(r)
(x − x∗ )T ∇2 φ(x∗ )(x − x∗ ).
2

32

Algorithm 2: Volume(M , ε)
Let σ02 = Θ(ε2 n−3 log−3 (n/ε)),
2
σi+1

 

σ 2 1 + √1
if ϑ ≤ nσi2
i
n


=
.
σ 2 1 + min( √σi , 1 )
otherwise
i
2
ϑ

and

ki =

( √
Θ( ε2n log nε )
Θ((

√

ϑ
σi

if ϑ ≤ nσi2

+ 1)ǫ−2 log nε ) otherwise

.

Set i = 0. Compute x∗ = arg min φ(x).
Assume that φ(x∗ ) = 0 by shifting the barrier function.
Sample k0 points {X1 , . . . , Xk0 } from the Gaussian distribution centered at the minimizer x∗
−1
.
of φ with covariance σ02 ∇2 φ(x∗ )
nϑ
ϑ
2
while σi ≤ Θ(1) ε log ε do
Sample ki points {X1 , . . . , Xk } using Hamiltonian Monte Carlo with target density
f (σi2 , X) and the previous {X1 , . . . , Xk } as warm start.
Compute the ratio
ki
2 ,X )
f (σi+1
1 X
j
Wi+1 =
·
2, X ) .
ki
f
(σ
j
i
j=1
Increment i.
end
n

1

Output: (2πσ02 ) 2 det(∇2 φ(x∗ ))− 2 W1 . . . Wi as the volume estimate.

33

Lemma 39 ([27, Prop 2.3.2]). For any ϑ-self-concordance barrier function φ on convex K, for any
interior point x and y in K, we have that
φ(x) ≤ φ(y) + ϑ ln

1
1 − πy (x)

where πy (x) = inf{t ≥ 0|y + t−1 (x − y) ∈ K}.
Lemma 40 ([23, Lem 5.16]). For any logconcave distribution f on Rn and any β ≥ 2, we have
Px∼f (f (x) ≤ e−βn max f (y)) ≤ e−O(βn) .
y

3

3

Lemma 41 (Large and small σ 2 ). Let φ be a ϑ-self concordant barrier function for K. If σn 2 log 2
1
1
1 and σn 2 log 2 σ1 ≤ 1, then we have that

1
σ

≤



3
1
3 1
n
−2
∗
ln F (σ) − ln e−σ φ(x ) (2πσ 2 ) 2 det(∇2 φ(x∗ ))− 2 ≤ O(σn 2 log 2 ).
σ

If σ 2 ≥ ϑ, then



−2
∗
ln F (σ) − ln e−σ φ(x ) vol(K) ≤ O(σ −2 ϑ ln(σ 2 n/ϑ)).

Proof. We begin with the first inequality. Let S be the set of x such that φ(x) ≤ φ(x∗ ) + βσ 2 n for
some β ≥ 2 to be determined. Therefore,
f (σ 2 , x) ≥ e−βn f (σ 2 , x∗ )
for all x ∈ S. Since f (σ 2 , x) is a logconcave function, Lemma 40 shows that Px∼f (S) ≤ e−O(βn) .
Therefore,
Z
Z
Z
f (σ 2 , x) dx ≤ f (σ 2 , x) dx ≤ (1 + e−O(βn) ) f (σ 2 , x) dx.
S

S

In short, we have

−O(βn)

F (σ) = (1 ± e

)

Z

f (σ 2 , x) dx.
S

By the construction of S, Lemma 38 and the fact that φ is self-concordant, if βσ 2 n ≤ 14 , we have
that
√
1 ± Θ(σ βn)
φ(x) = φ(x∗ ) +
(x − x∗ )T ∇2 φ(x∗ )(x − x∗ )
2
for all x ∈ S. Hence,
Z
√
σ −2
−2
∗
∗ T 2
∗
∗
F (σ) = (1 ± e−O(βn) )e−σ φ(x ) e−(1±O(σ βn)) 2 (x−x ) ∇ φ(x )(x−x ) dx.
S

Now we note that
Z
Z
σ −2
σ −2
∗ T 2
∗
∗
∗ T 2
∗
∗
e−Θ( 2 (x−x ) ∇ φ(x )(x−x )) dx = e−O(βn) e− 2 (x−x ) ∇ φ(x )(x−x ) dx
Sc

because

S

σ −2
(x − x∗ )T ∇2 φ(x∗ )(x − x∗ ) = Ω(βn)
2

34

outside S. Therefore,
F (σ) = (1 ± e−O(βn) )e−σ
−O(βn)

= (1 ± e

−2 φ(x∗ )

Z
3
2

e−(1±O(σ

√

Rn
−σ−2 φ(x∗ )

± O(σ(βn) ))e
3

= (1 ± e−O(βn) ± O(σ(βn) 2 ))e−σ

−2 φ(x∗ )

−2

βn)) σ 2 (x−x∗ )T ∇2 φ(x∗ )(x−x∗ )

Z

e−

dx

σ −2
(x−x∗ )T ∇2 φ(x∗ )(x−x∗ )
2

dx

Rn
n

1

(2πσ 2 ) 2 det(∇2 φ(x∗ ))− 2

1

3

where we used that σ(βn) 2 = O(1) in first sentence and σ(βn) 2 = O(1) in the second sentence.
Setting β = Θ(log σ1 ), we get the first result.
For the second inequality, for any 0 ≤ t < 1 and any x ∈ x∗ + t(K − x∗ ), we have that πx∗ (x) ≤ t
(πx∗ is defined in Lemma 39). Therefore, Lemma 39 shows that
φ(x) ≤ φ(x∗ ) + ϑ ln

1
.
1−t

Note that Px∼µ (x ∈ x∗ + t(K − x∗ )) = tn where µ is the uniform distribution in K. Therefore, for
any 0 < β < 1, we have that
Px∼µ (φ(x) ≤ φ(x∗ ) + ϑ ln

1
}) ≥ (1 − β)n .
β

Hence,
vol(K) · e−σ

−2 φ(x∗ )

≥ F (σ 2 )
n



≥ vol(K) · (1 − β) exp −σ

−2


1
(φ(x ) + ϑ ln ) .
β
∗

Setting β = σ −2 n−1 ϑ, we get the second result.
6.2.2

Variance of Yx

Our goal is to estimate Ex∼µi (Yx ) within a target relative error. The algorithm estimates the
quantity Ex∼µi (Yx ) by taking random sample points x1 , . . . , xk and computing the empirical estimate
for Ex∼µi (Yx ) from the corresponding Yx1 , . . . , Yxk . The variance of Yxi divided by its expectation
squared will give a bound on how many independent samples xi are needed to estimate Ex∼µi (Yx )
within the target accuracy. We have


R
σ2 σi2
2φ(x)
φ(x)
−
dx
)
exp
F ( 2σ2i+1
2
2
2
K
σ
σ
i
i+1
i −σi+1


Ex∼µi (Yx2 ) =
=
R
φ(x)
F (σi2 )
dx
K exp − σ2
i

and

σ2

σ2

i
)
F (σi2 )F ( 2σ2i+1
2
Ex∼µi (Yx2 )
i −σi+1
=
2 )2
Ex∼µi (Yx )2
F (σi+1

2
If we let σ 2 = σi+1
and σi2 = σ 2 /(1 + r), then we can further simplify as
 2   2 
σ
σ
2
F
1+r F 1−r
Ex∼µi (Yx )
=
.
Ex∼µi (Yx )2
F (σ 2 )2

35

(6.1)

Lemma 42. For any 1 > r ≥ 0, we have that
  2   2 
σ
σ
Z r Z 1+t
F 1−r
F 1+r
1


= 4
Varx∼µs φ(x)dsdt
ln
σ 0 1−t
F (σ 2 )2
2

where µs be the probability distribution proportional to f ( σs , x).
2

Proof. Fix σ 2 . Let g(t) = ln F ( σt ). Then, we have that
  2   2 
  2   2 
σ
σ
σ
σ
Z r
F 1−r
F 1+r
d  F 1+t F 1−t 


=
dt
ln
ln
F (σ 2 )2
F (σ 2 )2
0 dt
Z r
d
d
=
g(1 + t) − g(1 − t)dt
dt
0 dt
Z r Z 1+t 2
d
=
g(s)dsdt.
2
ds
0
1−t
For

d2
g(s),
ds2

(6.2)

we have that
Z

 s

exp − 2 φ(x) dx
σ
K
R

1 d K φ(x) · exp − σs2 φ(x) dx
R

=− 2 ·
s
σ ds
K exp − σ2 φ(x) dx

 2 !
R
 2 R
s
2 (x) · exp − s φ(x) dx
φ(x)
dx
φ(x)
·
exp
−
φ
1
2
2
K R
σ 
σ
− KR
=
 2
s
2
s
σ
K exp − σ2 φ(x) dx
K exp − σ2 φ(x) dx


1
1
= 4 Ex∼µs φ2 (x) − (Ex∼µs φ(x))2 = 4 Varx∼µs φ(x).
σ
σ

d2
d2
g(s)
=
ln
ds2
ds2

Putting it into (6.2), we have the result.

Now, we bound Varx∼µs φ(x). This can be viewed as a manifold version of the thin shell or
variance hypothesis estimate.
Lemma 43 (Thin shell estimates). Let φ be a ϑ-self concordant barrier function for K. Then, we
have that
 2 
σ
Varx∼µs φ(x) = O
ϑ .
s
def

Proof. Let Kt = {x ∈ K such that φ(x) ≤ t} and m be the number such that µs (Km ) = 21 . Let
Km,r = {x such that d(x, y) ≤ r and y ∈ Km }. By repeatedly applying Lemma 37, we have that
µs (Km,r ) = 1 − e−Ω(

√

s
r)
σ

.

√
By our assumption on φ, for any x and y, we have that |φ(x) − φ(y)| ≤ ϑd(x, y). Therefore,√for
√
s
any x ∈ Km,r , we have that φ(x) ≤ m + ϑr. Therefore, with probability at least 1 − e−Ω( √σ r)
√
√
s
φ(x) ≥ m − ϑr. Hence, with 1 − e−Ω( σ r)
in µs , it follows that φ(x) ≤ m + ϑr. Similarly,
√
probability in µs , we have that |φ(x) − m| ≤ ϑr. The bound on the variance follows.
36

Now we are ready to prove the key lemma.
Lemma 44. Let φ be a ϑ-self concordant barrier function for K. For any

Ex∼µi (Yx2 )
= O r 2 min
2
Ex∼µi (Yx )



1
2


ϑ
,n .
σi2

> r ≥ 0, we have that

Proof. Using Lemma 43 and Lemma 42, we have that
  2   2 
σ
σ
 2 
F 1−r
F 1+r
r ϑ


=O
.
ln
2
σ2
F (σ 2 )

(6.3)

This bounds is useful when σ 2 is large.
case σ 2 is small, we recall that for any logconcave function f , the function a →
RFor the
n
a
a
f (x) dx is logconcave (Lemma 3.2 in [8]). In particular, this shows that an F a1 is logconcave
1 1−r
in a. Therefore, with a = 1+r
σ2 , σ2 , σ2 ,we have
1
F (σ 2 )2 ≥
σ 4n



1+r
σ2

n

F

Rearranging the term, we have that
 2   2 
σ
σ
F 1−r
F 1+r
F (σ 2 )2

Therefore, we have that



ln 

F



σ2
1+r

F



F



σ2
1+r

≤




σ2
1−r

(σ 2 )2

1−r
σ2

F



n

.

n

1
(1 + r)(1 − r)

σ2
1−r



.



Combining (6.1), (6.3) and (6.4), we have the result.
6.2.3




 = O nr 2 .

(6.4)

Main lemma

Lemma 45. Given any ϑ-self-concordance barrier φ on a convex set K and 0 < ε < 12 , the algorithm
Volume(M , ε) outputs the volume of K to within a 1 ± ε multiplicative factor.
−2

∗

n

1

Proof. By our choice of ε, Lemma 41 shows that e−σ0 φ(x ) (2πσ02 ) 2 det(∇2 φ(x∗ ))− 2 is an 1± 4ε mul−2
∗
tiplicative approximation of F (σ0 ) and that e−σk φ(x ) vol(K) is a 1± 4ε multiplicative approximation
of F (σk ). Note that we shifted the function φ such that φ(x∗ ) = 0. Therefore,
k
2
n
1 Y F (σi+1 )
ε
vol(K) = (1 ± )(2πσ02 ) 2 det(∇2 φ(x∗ ))− 2
2) .
2
F
(σ
i
i=1
2 , X)/f (σ 2 , X) is upper
In Lemma 44, we showed that the variance of the estimator Y = f (σi+1
i
q
√
2
bounded by O(1)(EY ) . Note that the algorithm takes O( n) iterations to double σi if ϑn ≤ σi
√
and O( ϑσi−1 ) iterations otherwise. By a simple analysis of variance, to have relative error ε, it
e i ) samples in each phase.
suffices to have O(k

37

6.3

Volume computation with the log barrier

In this section, we prove the Theorem 2, restated below for convenience.
Theorem 2. For any polytope P = {x : Ax ≥ b} with m constraints and n variables, and any
ε > 0, the Hamiltonian
volume algorithm estimates the volume of P to within 1 ± ε multiplicative


2
−2
e mn 3 ε
factor using O
steps where each step consists of solving a first-order ODE and takes time


e mnω−1 LO(1) logO(1) 1 and L is the bit complexity3 of the polytope.
O
ε
1
m
e
n− 3 ). Since the number
n , the mixing time of HMC is O(m ·√
√
e 2n ), the total number of
σ 2 is O( n) and since we samples O(
ǫ

Proof. In the first part, when σ 2 ≤

of sampling phases to double such
steps of HMC is
√ 


√ 
n
− 31
e
e
e
O mn
=O
×O n ×O
2
ǫ
In the second part, when σ 2 ≥

m
n,

2

m · n3
ǫ2

!

.

√

1

2

e −2n 3 −1 +
the mixing time of HMC is O(
σ +m

1

m3 n3

σ

−2
3

1

1

3
+m−√

1

+ m 2 n 6 ).

m
−2
e
Since the number of sampling phases to double σ 2 is O(1+ σm ) and since we sample O((
σ +1)ε )
in each phase, the total number of steps of HMC is
!
!
 √


√ 
1
2
2
1
3
3 n3
3
1
1
m
m
·
n
n
m
m
−2
e (
e
e
+ m2 n6 × O 1 +
+
+ 1)ε
.
×O
=O
O
σ −2 + m−1 σ − 32 + m− 13
σ
σ
ǫ2

Combining both parts, the total number of steps of HMC is
!
2
3
m
·
n
e
.
O
ǫ2

3

A.

L = log(m + dmax + kbk∞ ) where dmax is the largest absolute value of the determinant of a square sub-matrix of

38

7

Logarithmic barrier

For any polytope M = {Ax > b}, the logarithmic barrier function φ(x) is defined as
φ(x) = −

m
X
i=1

log(aTi x − bi ).

We denote the manifold induced by the logarithmic barrier on M by ML . The goal of this section
is to analyze Hamiltonian Monte Carlo on ML . In Section 7.1, we give explicit formulas for various
Riemannian geometry concepts on ML . In Section 7.2, we describe the HMC specialized to ML .
In Sections 7.3 to 7.5, we bound the parameters required by Theorem 30, resulting in Theorem 1.
The following parameters that are associated with barrier functions will be convenient.
Definition 46. For a convex function f , let M1 , M2 and M3 be the smallest numbers such that
−1
1. M1 ≥ maxx∈M (∇f (x))T ATx Ax
∇f (x) and M1 ≥ n.
2. ∇2 f  M2 · ATx Ax .

3. Tr((ATx Ax )−1 ∇3 f (x)[v]) ≤ M3 kvkx for all v.
For the case f = φ are the standard logarithmic barrier, these parameters are n, 1,
tively.

7.1

√

n respec-

Riemannian geometry on ML (G2 )

We use the following definitions throughout this section.
Definition 47. For any matrix A ∈ Rm×n and vectors b ∈ Rm and x ∈ Rn , define
1. sx = Ax − b, Sx = Diag(sx ), Ax = Sx−1 A.
2. sx,v = Ax v, Sx,v = Diag(Ax v).


(2)
= (Px )2ij .
3. Px = Ax (ATx Ax )−1 ATx , σx = diag(Px ), Σx = Diag(P ), Px
ij

4. Gradient of φ: φi = −

P

ℓ


eTℓ Ax ei .

−1



P T
ej .
eℓ Ax ei eTℓ Ax ej , gij = eTi ATx Ax
5. Hessian of φ and its inverse: gij = φij = ATx Ax ij =



P
6. Third derivatives of φ: φijk = −2 ℓ eTℓ Ax ei eTℓ Ax ej eTℓ Ax ek .
7. For brevity (overloading notation), we define sγ ′ = sγ,γ ′ , sγ ′′ = sγ,γ ′′ , Sγ ′ = Sγ,γ ′ and
Sγ ′′ = Sγ,γ ′′ for a curve γ(t).

In this section, we will frequently use the following identities derived from elementary calculus
(using only the chain/product rules and the formula for derivative of the inverse of a matrix). For
reference, we include proofs in Appendix C.
Fact 48. For any matrix A and any curve γ(t), we have
dAγ
= −Sγ ′ Aγ ,
dt
dPγ
= −Sγ ′ Pγ − Pγ Sγ ′ + 2Pγ Sγ ′ Pγ ,
dt
dSγ ′
= Diag(−Sγ ′ Aγ γ ′ + Aγ γ ′′ ) = −Sγ2′ + Sγ ′′ ,
dt
39

We also use these matrix inequalities: Tr(AB) = Tr(BA), Tr(P AP ) ≤ Tr(A) for any psd
matrix A; Tr(ABAT ) ≤ Tr(AZAT ) for any B  Z; the Cauchy-Schwartz, namely, Tr(AB) ≤
1
1
Tr(AAT ) 2 Tr(BB T ) 2 . We note Px2 = Px because Px is a projection matrix.
Since the manifold ML is naturally embedded in Rn , we can identify Tx ML with Euclidean
coordinates. We have that
hu, vix = uT ∇2 φ(x)v = uT ATx Ax v.
We will use the following two lemmas proved in [13].
Lemma 49. Let w(t) be a vector field defined on a curve z(t) in ML . Then,
∇z ′ w =

−1 T
−1 T
dw
dw
Az Sz ′ sz,w =
Az Sz,w sz ′ .
− ATz Az
− ATz Az
dt
dt

In particular, the equation for parallel transport on a curve γ(t) is given by
−1 T
d
Aγ Sγ ′ Aγ v.
v(t) = ATγ Aγ
dt

(7.1)

Lemma 50. Given u, v, w, x ∈ Tx ML , the Riemann Curvature Tensor at x is given by
R(u, v)w =

ATx Ax
def

−1

ATx (Sx,v Px Sx,w − Diag(Px sx,v sx,w )) Ax u

and the Ricci curvature Ric(u) = TrR(u, u) is given by
Ric(u) = sTx,u Px(2) sx,u − σxT Px s2x,u
where R(u, u) is the operator defined above.

7.2

Hamiltonian walk on ML

We often work in Euclidean coordinates. In this case, the Hamiltonian walk is given by the formula
in the next lemma. To implement the walk, we solve this ODE using the collocation method as
described in [13], after first reducing it to a first-order ODE. The resulting complexity is Õ(mnω−1 )
per step.
Lemma 51. The Hamiltonian curve at a point x in Euclidean coordinates is given by the following
equations
γ ′′ (t) =

ATγ Aγ

′

γ (0) = w,

−1

ATγ s2γ ′ + µ(γ(t))

∀t ≥ 0

γ(0) = x.
where µ(x) = ATx Ax

−1

ATx σx − (ATx Ax )−1 ∇f (x) and w ∼ N (0, (ATγ Aγ )−1 ).

Proof. Recall from Lemma 7 that the Hamiltonian walk is given by
Dt

dγ
=µ(γ(t)),
dt

dγ
(0) ∼N (0, g(x)−1 )
dt

40



where µ(x) = −g(x)−1 ∇f (x)− 12 g(x)−1 Tr g(x)−1 Dg(x) . By Lemma 49, applied with w(t) = γ ′ (t),
z(t) = γ(t), we have
−1 T 2
dγ
Dt
Aγ sγ ′ .
= γ ′′ (t) − ATγ Aγ
dt
For the formula of µ, we note that
 ∂


1 
1X
(ATx Ax )−1 ij
Tr g(x)−1 Dg(x) k =
ATx Ax ji
2
2
∂xk
ij
X
X




(ATx Ax )−1 ij
eTℓ Ax ei eTℓ Ax ej eTℓ Ax ek
by Defn.47(6) = −
ij

=−

=

X

ℓ

ATx (ATx Ax )−1 Ax

ℓ
−ATx σx .



ℓℓ

eTℓ Ax ek

Therefore,
µ(x) = −(ATx Ax )−1 ∇f (x) + ATx Ax

−1



ATx σx .

Many parameters for Hamiltonian walk depends on the operator Φ(t). Here, we give a formula
of Φ(t) in Euclidean coordinates.
Lemma 52. Given a curve γ(t), in Euclidean coordinates, we have that
Φ(t) = M (t) − R(t)
where
−1


ATγ Sγ ′ Pγ Sγ ′ Aγ − ATγ Diag(Pγ s2γ ′ )Aγ


−1  T 
Ax Sx,µ − 3Σx + 2Px(2) Ax − ∇2 f (x) .
M (t) = ATγ Aγ
R(t) = ATγ Aγ

Proof. Lemma 50 with v = w = γ ′ , gives the formula for R(t).
For M (t), Lemma 49 with w = µ(x), z ′ = u shows that

Du µ(x) = ∇u µ(x) − (ATx Ax )−1 ATx Sx,µ Ax u.
For the first term ∇u µ(x), we note that

−1

−1 T
ATx Sx,u Ax ATx Ax
Ax σx

−1
ATx Sx,u σx
− 3 ATx Ax
−1 T
+ 2 ATx Ax
Ax diag(Px Sx,u Px )
−1

−1
∇f (x)
ATx Sx,u Ax ATx Ax
− 2 ATx Ax

∇u µ(x) =2 ATx Ax

− (ATx Ax )−1 ∇2 f (x)u.

41

Therefore, we have that
Du µ(x)

−1 T (2)
−1 T
Ax Px Ax u
Ax Σx Ax u + 2 ATx Ax
ATx Diag(Px σx )Ax u − 3 ATx Ax




−1
−1 T
∇f (x) Ax u − (ATx Ax )−1 ∇2 f (x)u − (ATx Ax )−1 ATx Sx,µ Ax u
Ax Diag Ax ATx Ax
− 2 ATx Ax




−1
−1  T 
∇f (x) − Sx,µ − 3Σx + 2Px(2) Ax − ∇2 f (x) u
Ax 2Diag(Px σx ) − 2Diag Aγ ATx Ax
= ATx Ax


−1  T 
= ATx Ax
Ax 2Sx,µ − Sx,µ − 3Σx + 2Px(2) Ax − ∇2 f (x) u.


−1  T 
= ATx Ax
Ax Sx,µ − 3Σx + 2Px(2) Ax − ∇2 f (x) u.
=2 ATx Ax

−1

where we used the facts that

µ = ATx Ax

−1

ATx σx − (ATx Ax )−1 ∇f (x)

Sx,µ = Ax µ = Diag(Px σx − Ax (ATx Ax )−1 ∇f (x)).
Remark. Note that R(t) and M (t) is symmetric in h·, ·iγ , but not in h·, ·i2 . That is why the formula
does not look symmetric.

7.3

Randomness of the Hamiltonian flow (ℓ0 )

Many parameters of a Hessian manifold relate to how fast a Hamiltonian curve approaches the
boundary of the polytope. Since the initial velocity of the
curve is drawn from a

 Hamiltonian
1
Gaussian distribution, one can imagine that sγ ′ (0) ∞ = O √m sγ ′ (0) 2 (each coordinate of sγ ′
measures the relative rate at which the curve is approaching the corresponding facet). So the walk
initial approaches/leaves every facet of the polytope at roughly the same slow pace. If this holds
for the entire walk, it would allow us to get very tight bounds on various parameters. Although we
are not able to prove that sγ ′ (t) ∞ is stable throughout 0 ≤ t ≤ δ, we will show that sγ ′ (t) 4 is
stable and thereby obtain a good bound on sγ ′ (t) ∞ .
Throughout this section, we only use the randomness of the walk to prove that both sγ ′ (t) 4 and
sγ ′ (t) ∞ are small with high probability. Looking ahead, we will show that sγ ′ (t)
√
√
and sγ ′ (t) ∞ = O( log n + M1 δ) (Lemma 55), we define
sγ ′ (t)

def

ℓ(γ) = max

0≤t≤δ

n1/2

+

2
1/4
M1

+

sγ ′ (t)
1/4
M1

4

sγ ′ (t)

√∞
+√
+
log n + M1 δ

sγ ′ (0)
n1/2

2

+

sγ ′ (0)
n1/4

4

1/4

4

= O(M1 )

sγ ′ (0) ∞
+ √
log n

!

to capture this randomness involves in generating the geodesic walk. This allows us to perturb the
geodesic (Lemma 24) without worrying about the dependence on randomness.
We first prove the the walk is stable in the L4 norm and hence ℓ(γ) can be simply approximated
by sγ ′ (0) 4 and sγ ′ (0) ∞ .
Lemma 53. Let γ be a Hamiltonian flow in ML starting at x. Let v4 =
1
0≤t≤
1/4 , we have that
12(v4 +M1

)

42

sγ ′ (0)

4

. Then, for

1.

sγ ′ (t)

1/4

4

≤ 2v4 + M1 .

2. kγ ′′ (t)k2γ(t) ≤ 128v44 + 30M1 .
Proof. Let u(t) = sγ ′ (t)

4

. Then, we have (using Holder’s inequality in the first step),

d
du
≤
Aγ γ ′
dt
dt
≤ Aγ γ

′′

4

4
2

= Aγ γ ′′ − Aγ γ ′

+ u (t).

2

4

(7.2)

Under the Euclidean coordinates, by Lemma 51 the Hamiltonian flow is given by
γ ′′ (t) = ATγ Aγ
with µ(x) = ATx Ax

−1

γ ′′

−1

ATγ s2γ ′ + µ(γ(t))

ATx σx − (ATx Ax )−1 ∇f (x). Hence, we have that

2
γ

−1 T  T −1 T 2
Aγ sγ ′
Aγ Aγ Aγ Aγ
Aγ ATγ Aγ
−1 T  T −1 T
+ 3σγT Aγ ATγ Aγ
Aγ Aγ Aγ Aγ
Aγ σγ


−1
−1
+ 3(∇f (x))T ATγ Aγ
ATγ Aγ ATγ Aγ
∇f (x)
X
X

−1
≤3
(s4γ ′ )i + 3
(σγ2 )i + 3(∇f (x))T ATγ Aγ
∇f (x)
≤3 s2γ ′

T

i
4

i

≤3u (t) + 3(n + M1 ) ≤ 3u4 (t) + 6M1
Therefore, we have
Aγ γ ′′

4

≤ Aγ γ ′′

2

Plugging it into (7.2), we have that

1/4

(7.3)

p
≤ 2u2 (t) + 3 M1 .

p
du
≤ 3u2 (t) + 3 M1 .
dt

Note that when u ≤ 2v4 + M1 , we have that
p
du
1/4
≤ 12v42 + 9 M1 ≤ 12(v4 + M1 )2 .
dt

Since u(0) = v4 , for 0 ≤ t ≤

1
1/4 ,
12(v4 +M1 )

1/4

we have that u(t) ≤ 2v4 + M1

and this gives the first

inequality.
Using (7.3), we get the second inequality.

We can now prove that ℓ(γ) is small with high probability.
Lemma 54. Assume that δ ≤

1
1/4 ,
36M1

ℓ1 = Ω(n1/4 δ) and n is large enough, we have that



ℓ0
1
min 1,
Pγ∼x (ℓ(γ) ≥ 128) ≤
.
100
ℓ1 δ
Therefore, we have ℓ0 ≤ 256.
43

Proof. From the definition of the Hamiltonian curve (Lemma 7), we have that
Aγ γ ′ (0) = Bz
−1/2
and z ∼ N (0, I).
where B = Aγ ATγ Aγ
First, we estimate kAγ γ ′ (t)k4 . Lemma 65 shows that




Pz∼N (0,I) kBzk44 ≤  3

Note that

P

eTi B

i

4
2

=

P

2
i (σγ )i

Pγ ′ (0)
def

X

eTi B

i

4
2

!1/4

4 

2

s
+ kBk2→4 s  ≤ 1 − exp(− ).
2

≤ n and kBk2→4 ≤ kBk2→2 = 1. Putting s =



′

Aγ γ (0)

4
4

n1/4
2 ,

we have that

√
n
).
≤ 11n ≤ 1 − exp(−
8


Therefore, we have that v4 = kAγ γ ′ (0)k4 ≤ 2n1/4 with probability at least 1 − exp(−
apply Lemma 53 to get that
1/4
1/4
sγ ′ (t) 4 ≤ 2v4 + M1 ≤ 5M1

√
n
8 ).

Now, we

1
1/4 .
12(v4 +M1 )

for all 0 ≤ t ≤

Next, we estimate kAγ γ ′ (t)k∞ . Since eTi Aγ γ ′ (0) = eTi Bx ∼ N (0, σi ), we have
 2
√ 
t
T
′
.
Pγ ′ (0) ei Aγ γ (0) ≥ σi t ≤ 2 exp −
2

Hence, we have that

Pγ ′ (0)


′

Aγ γ (0)

∞




X
p
2 log n
≥ 2 log n ≤ 2
exp −
σi
i



P
P
n
log n
Since i exp − 2 log
is
concave
in
σ,
the
maximum
of
exp
−
on the feasible set {0 ≤
i
σi
σi
P
σ ≤ 1, σi = n} occurs on its vertices. Hence, we have that


p
2
Pγ ′ (0) Aγ γ ′ (0) ∞ ≥ 2 log n ≤ 2n exp (−2 log n) = .
n
√
√
Lemma 53 shows that kAγ γ ′′ k∞ ≤ kγ ′′ kγ(t) ≤ 46 n + 6 M1 . Hence, for any 0 ≤ t ≤ δ, we have
that

sγ ′ (t)





∞

Z

t

Aγ(t) γ ′′ (r) ∞ dr
≤ Aγ(t) γ (0) ∞ +
0


Z δ
sγ(t),i
′
≤ max
Aγ(r) γ ′′ (r)
Aγ(0) γ (0) ∞ +
i,0≤s≤t sγ(s),i
0



p
sγ(t),i  p
√
≤ max
2 log n + (46 n + 6 M1 )δ
i,0≤s≤t sγ(s),i



p
sγ(t),i  p
≤ max
2 log n + 52 M1 δ .
i,0≤s≤t sγ(s),i
′

44

∞

dr



(7.4)

sγ(t),i
sγ(s),i

Let z(t) = maxi,0≤s≤t
sγ(t),i = sγ(s),i exp

Z

. Note that


sγ ′ (α),i dα

t

r

since

(Aγ(t)−b)i = (Aγ(r)−b)i exp

Z

t
r


ai · γ ′ (α)
dα
(Aγ(α) − b)i

Hence, we have that
z ′ (t) ≤ z(t) sγ ′ (t)

∞

Solving this, since z(0) = 1, we get
z(t) ≤
Since t ≤ δ ≤

1
1/4 ,
36M1

1
 .
√
√
1 − 2 log n + 52 M1 δ t

we have that z(t) ≤ 1.05. Putting this into (7.4), we have that
sγ ′ (t)

Finally, we estimate sγ ′ (t)

2

X

Pz∼N (0,I) kBzk22 ≤ 

P

i

eTi B

2
2

=

P

∞

p
p
≤ 3 log n + 55 M1 δ.

= kAγ γ ′ (t)k2 . Lemma 65 shows that





Note that


 p
p
≤ z 2 (t) 2 log n + 52 M1 δ .

eTi B

i

2
2

!1/2

2 

2

r
+ kBk2→2 r   ≤ 1 − exp(− ).
2
1/2

≤ n and kBk2→2 ≤ 1. Putting s = n 3 , we have that


n
2
Pγ ′ (0) Aγ γ ′ (0) 2 ≤ 2n ≤ 1 − exp(− ).
18

i (σγ )i

Therefore, kAγ γ ′ (0)k22 ≤ 2n with high probability. Next, we note that
d
Aγ γ ′
dt

2
2

= Aγ γ ′ , Aγ γ ′′ − s2γ ′
E
D
−1 T 2
Aγ sγ ′ + Aγ µ(γ(t)) − Aγ s2γ ′
= Aγ γ ′ , Aγ ATγ Aγ
E
D
−1 T 2
Aγ sγ ′ − s2γ ′ + Aγ γ ′ , Aγ µ(γ(t))
= Aγ γ ′ , Aγ ATγ Aγ
= Aγ γ ′ , Aγ µ(γ(t)) .

Using that µ(x) = ATx Ax
d
Aγ γ ′
dt

2
2

=

−1

ATx σx − (ATx Ax )−1 ∇f (x), we have that

X
X
(γ ′ )i (∇f )i
(sγ ′ )i (σγ )i −
i

i

sX
sX
qX
q
2
2
≤
(sγ ′ )i
(σγ )i +
(sγ ′ )2i (∇f )T (ATγ Aγ )−1 (∇f )
i

≤ 2 Aγ γ ′

2

p

i

M1 .
√
≤ M1 . Since δ ≤

d
1
′
Therefore, we have that dt
kAγ γ ′ k2
1/4 , we have that kAγ γ (t)k2 ≤
36M1
1/4
1/4
√
√
M1
M1
1/4
kAγ γ ′ (0)k2 + 36
≤ 2n + 36
. Therefore, we have that kAγ γ ′ (t)k2 ≤ 2 n + M1 with
n
probability at least 1 − exp(− 18 ).

45

Combining our estimates on sγ ′ (t)
we have that
sγ ′ (t)

2
1/4
1/2
n + M1

max

P

0≤t≤δ

+

sγ ′ (t)

4

1/4

M1

2

, sγ ′ (t)

+√

√
n
2
n
≤ exp(− ) + exp(−
)+ .
18
8
n

sγ ′ (t)
log n +

4

and sγ ′ (t)

√∞

M1 δ

+

∞

sγ ′ (0)
n1/2

and using the assumption on δ,

2

+

sγ ′ (0)
n1/4

4

sγ ′ (0) ∞
+ √
log n

!

!

≥ 128

In Lemma 59, we indeed have that ℓ1 = Ω(n1/4 δ) and hence ℓℓ10δ = Ω(n−1/4 δ−2 ) = Ω(1). Therefore,


1
min 1, ℓℓ10δ when n is large enough.
the probability is less than 100
Here, we collect some simple consequences of small ℓ(γ) that we will use later.

Lemma 55. Given a Hamiltonian flow γ on ML with ℓ(γ) ≤ ℓ0 ≤ 256. For any 0 ≤ t ≤ δ,
1/4

1. kAγ γ ′ (t)k2 ≤ 256(n1/2 + M1 ), kAγ γ ′ (0)k2 ≤ 256n1/2 .
1/4

2. kAγ γ ′ (t)k4 ≤ 256M1 , kAγ γ ′ (0)k4 ≤ 256n1/4 .

√
√
√
3. kAγ γ ′ (t)k∞ ≤ 256 log n + M1 δ , kAγ γ ′ (0)k∞ ≤ 256 log n.
4. kγ ′′ (t)k2γ ≤ 1013 M1 .

Proof. The first three inequalities simply follow from the definition of ℓ(γ). Since kAγ γ ′ (0)k4 ≤
1/4

256M1 , Lemma 53 shows the last inequality.

7.4

Parameters R1 , R2 and R3

Lemma 56. For a Hamiltonian curve γ on ML with ℓ(γ) ≤ ℓ0 , we have that
sup kΦ(t)kF,γ(t) ≤ R1

0≤t≤ℓ

√
√
with R1 = O( M1 + M2 n).
Proof. Note that Φ(t) = M (t) − R(t) where

ATγ Sγ ′ Pγ Sγ ′ Aγ − ATγ Diag(Pγ s2γ ′ )Aγ ,

−1  T
M (t) = ATγ Aγ
Aγ (Sγ,µ − 3Σγ + 2Pγ(2) )Aγ − ∇2 f .
R(t) = ATγ Aγ

−1

We bound the Frobenius norm of Φ(t) separately.
For kR(t)kF,γ , we note that
kR(t)k2F,γ
≤2 (ATγ Aγ )−1/2 ATγ Sγ ′ Pγ Sγ ′ Aγ (ATγ Aγ )−1/2
=2TrPγ S Pγ S Pγ S Pγ S +
γ′

≤4 sγ ′

γ′

γ′

γ′

2
F

+ 2 (ATγ Aγ )−1/2 ATγ Diag(Pγ s2γ ′ )Aγ (ATγ Aγ )−1/2

2TrPγ Diag(Pγ s2γ ′ )Pγ Diag(Pγ s2γ ′ )

4
.
4

46

2
F

For kM (t)kF,γ , we note that
kM (t)k2F,γ
≤2 (ATγ Aγ )−1/2 ATγ (Sγ,µ − 3Σγ + 2Pγ(2) )Aγ (ATγ Aγ )−1/2

2
F

+ 2 (ATγ Aγ )−1/2 ∇2 f (ATγ Aγ )−1/2

2
F

≤6TrPγ Sγ,µ Pγ Sγ,µ + 54TrPγ Σγ Pγ Σγ + 24TrPγ Pγ(2) Pγ Pγ(2) + 2M22 n

2
≤6TrSγ,µ
+ 54TrΣ2γ + 24Tr(Pγ(2) )2 + 2M22 n.
(2)

where we used that diag(Pγ )

2

≤ kσγ k2 ≤ n and

TrPγ M Pγ M = TrPγ M Pγ M Pγ = TrPγ M M Pγ
= TrM Pγ Pγ M ≤ TrM 2 .
Note that
2
2

ksγ,µ k22 = kµ(γ)k2γ = Pγ σγ − Aγ (ATγ Aγ )−1 ∇f (γ)
and
Aγ ATγ Aγ
Therefore, we have that

−1

∇f (γ(t))

2
2

= ∇f (γ(t))T ATγ Aγ

≤ 2n + 2M1 ≤ 4M1

−1

∇f (γ(t)) = M1 .

kM (t)k2F,γ ≤ 102M1 + 2M22 n.

The claim follows from Lemma 55.
Lemma 57. Let γ be a Hamiltonian curve on ML with ℓ(γ) ≤ ℓ0 . Assume that δ2 ≤ √1n . For any
0 ≤ t ≤ δ, any curve c(r) starting from γ(t) and any vector field v(r) on c(r) with v(0) = γ ′ (t), we
have that
!
dc
d
+ δ k Dr v|r=0 kγ(t) .
TrΦ(v(r))|r=0 ≤ R2
dr
dr r=0 γ(t)
where

1/4
p
√
M
R2 = O( nM1 + nM1 δ2 + 1 +
δ

√

n log n √
+ nM2 + M3 ).
δ

Proof. We first bound TrR(t). By Lemma 50, we know that
(2)

T
Ric(v(r)) = sTc(r),v(r) Pc(r) sc(r),v(r) − σc(r)
Pc(r) s2c(r),v(r)

= Tr(Sc(r),v(r) Pc(r) Sc(r),v(r) Pc(r) ) − Tr(Diag(Pc(r) s2c(r),v(r) )Pc(r) ).

For simplicity, we suppress the parameter r and hence, we have
Ric(v) = Tr(Sc,v Pc Sc,v Pc ) − Tr(Diag(Pc s2c,v )Pc ).

47

d
d
c = c′ and dr
v = v ′ (in Euclidean coordinates). Since
We write dr
d
and dr
Sc,v = −Sc′ Sc,v + Sc,v′ , we have that

d
dr Pc

= −Sc′ Pc − Pc Sc′ + 2Pc Sc′ Pc

d
Ric(v)
dr
= −2Tr(Sc,v Sc′ Pc Sc,v Pc ) − 2Tr(Sc,v Pc Sc′ Sc,v Pc ) + 4Tr(Sc,v Pc Sc′ Pc Sc,v Pc )
−2Tr(Sc′ Sc,v Pc Sc,v Pc ) + 2Tr(Sc,v′ Pc Sc,v Pc )

+Tr(Diag(Pc s2c,v )Sc′ Pc ) + Tr(Diag(Pc s2c,v )Pc Sc′ ) − 2Tr(Diag(Pc s2c,v )Pc Sc′ Pc )

+Tr(Diag(Pc Sc′ s2c,v )Pc ) + Tr(Diag(Sc′ Pc s2c,v )Pc ) − 2Tr(Diag(Pc Sc′ Pc s2c,v )Pc )

+2Tr(Diag(Pc Sc,v Sc′ sc,v )Pc ) − 2Tr(Diag(Pc Sc,v sc,v′ )Pc )

= −6Tr(Sc,v Sc′ Pc Sc,v Pc ) + 4Tr(Sc,v Pc Sc′ Pc Sc,v Pc ) + 2Tr(Sc,v′ Pc Sc,v Pc )
+3Tr(Diag(Pc s2c,v )Sc′ Pc ) − 2Tr(Diag(Pc s2c,v )Pc Sc′ Pc )

+3Tr(Diag(Pc Sc′ s2c,v )Pc ) − 2Tr(Diag(Pc Sc′ Pc s2c,v )Pc )

−2Tr(Diag(Pc Sc,v sc,v′ )Pc ).

d
Let dr
Ric(v) = (1) + (2) where (1) is the sum of all terms not involving v ′ and (2) is the sum of
other terms.
For the first term (1), we have that

|(1)| ≤ 6 |Tr(Sc,v Sc′ Pc Sc,v Pc )| + 4 |Tr(Sc,v Pc Sc′ Pc Sc,v Pc )|

+3 Tr(Diag(Pc s2c,v )Sc′ Pc ) + 2 Tr(Diag(Pc s2c,v )Pc Sc′ Pc )

+3 Tr(Diag(Pc Sc′ s2c,v )Pc ) + 2 Tr(Diag(Pc Sc′ Pc s2c,v )Pc )
sX
sX
2
(sc,v )i
(sc,v )2i + 4 ksc′ k∞ |Tr(Pc Sc,v Pc Sc,v Pc )|
≤ 6 ksc′ k∞
i

i

sX
sX
sX
4
4
(sc,v )i kSc′ k2 + 2
(sc,v )i
(Pc Sc′ Pc )2ii
+3
i

i

i

sX
sX
sX
sX
(sc,v )4i
(Pc )2ii + 2 ksc′ k∞
(sc,v )4i
(Pc )2ii
+3 ksc′ k∞
i

≤

≤

i

i

10 ksc′ k∞ ksc,v k22 + 3 ksc,v k24 ksc′ k2
√
20 ksc′ k2 ksc,v k24 n.

+

√
7 ksc′ k∞ ksc,v k24 n

i

1/2

Since sc,v = sγ ′ at r = 0, we have that ksc,v k24 = O(M1 ) and hence
p

|(1)| = O
nM1 ksc′ k2 .
For the second term (2), we have that

|(2)| ≤ 2 Tr(Sc,v′ Pc Sc,v Pc ) + 2 Tr(Diag(Pc Sc,v sc,v′ )Pc )
sX
√
(sc,v′ sc,v )2i
≤ 2 sc,v′ 2 ksc,v k2 + 2 n
i




p
p
1/4
n log n + nM1 δ sc,v′
≤ O n1/2 + M1
sc,v′ 2 + O


p
p
1/4
= O M1 + n log n + nM1 δ sc,v′ 2
48

2

where we used ksc,v k∞ = sγ ′
1/4
M1 )

∞

=O

√

log n +

√


M1 δ and at r = 0, we have ksc,v k2 = Aγ(0) γ ′ (0)

sγ ′ 2 = O(n1/2 +
in the second-to-last line.
Note that at r = 0, by Lemma 49, we have
Dr v =

Therefore,

−1 T
dv
− ATc Ac
Ac Sc′ sc,v .
dr

sc,v′ = Ac v ′ = Ac (Dr v) − Ac ATc Ac
and hence
sc,v′

2

≤ kDr vk + Ac ATc Ac
≤ kDr vk + sγ ′

−1

−1

ATc Sc′ sc,v

ATc Sc′ sc,v

2

ksc′ k2
∞

Therefore,


p


p
p
p
1/4
|(2)| = O M1 + n log n + nM1 δ kDr vk +
log n + M1 δ ksc′ k2 .

Therefore, we have

d
Ric(v(r))|s=0
dr
p



p
p
1/4
=O
nM1 ksc′ k2 + O M1 + n log n + nM1 δ kDr vk
 p


p
p
p
1/4
log n + M1 δ ksc′ k2
+ O M1 + n log n + nM1 δ
p

p
p
√
√
1/4
1/4
=O
nM1 + M1
log n + M1
M1 δ + n log n + nM1 δ2 ksc′ k2
!
√
1/4
M1
n log n p
+
+ nM1 δ kDr vk
+O
δ
δ
!
√
1/4
p
√
n log n
M1
2
=O
(ksc′ k2 + δ kDr vk) .
+
nM1 + nM1 δ +
δ
δ

√
√
3/4
where we used M1 δ = O( nM1 δ2 + nM1 ) at the last line.
Next, we bound TrM (t). Lemma 52 shows that

TrM (r) =Tr((ATc Ac )−1 ATc (Sc,µ − 3Σc + 2Pc(2) )Ac − (ATc Ac )−1 ∇2 f )
=Tr(Pc (Sc,µ − 3Σc + 2Pc(2) )) − Tr((ATc Ac )−1 ∇2 f )

=σcT Pc σc − σcT Ac (ATc Ac )−1 ∇f − 3Tr(Σ2c ) + 2Tr(Pc(3) ) − Tr((ATc Ac )−1 ∇2 f ).
where in the last step we used
Sc,µ = Ac µ = Diag(Pc σc − Ac (ATc Ac )−1 ∇f ).

49

2

=

Since

d
dr Pc

= −Sc′ Pc − Pc Sc′ + 2Pc Sc′ Pc and

d
dr Ac

= −Sc′ Ac , we have that

d
TrM (r)
dr
= − σcT Sc′ Pc σc − σcT Pc Sc′ σc + 2σcT Pc Sc′ Pc σc
− 4σcT Pc Sc′ σc + 4σcT Pc diag(Pc Sc′ Pc )

+ 3σcT Sc′ Ac (ATc Ac )−1 ∇f − 2diag(Pc Sc′ Pc )T Ac (ATc Ac )−1 ∇f

− 2σcT Ac (ATc Ac )−1 ATc Sc′ Ac (ATc Ac )−1 ∇f − σcT Ac (ATc Ac )−1 ∇2 f · c′

+ 12Tr(Sc′ Σ2c ) − 12Tr(Σc diag(Pc Sc′ Pc ))

− 6Tr(Pc(2) Sc′ Pc + Pc(2) Pc Sc′ ) + 12Tr(Pc(2) Pc Sc′ Pc )

− 2Tr((ATc Ac )−1 ATc Sc′ Ac (ATc Ac )−1 ∇2 f ) − Tr((ATc Ac )−1 ∇3 f [c′ ]).
Simplifying it, we have
d
TrM (r)
dr
= − 6σcT Sc′ Pc σc + 2σcT Pc Sc′ Pc σc + 4σcT Pc Pc(2) sc′

+ 3σcT Sc′ Ac (ATc Ac )−1 ∇f − 2sTc′ Pc(2) Ac (ATc Ac )−1 ∇f

− 2σcT Ac (ATc Ac )−1 ATc Sc′ Ac (ATc Ac )−1 ∇f − σcT Ac (ATc Ac )−1 ∇2 f · c′

+ 12Tr(Sc′ Σ2c ) − 12σcT Pc(2) sc′

− 6Tr(Pc(2) Sc′ Pc + Pc(2) Pc Sc′ ) + 12Tr(Pc(2) Pc Sc′ Pc )

− 2Tr((ATc Ac )−1 ATc Sc′ Ac (ATc Ac )−1 ∇2 f ) − Tr((ATc Ac )−1 ∇3 f [c′ ]).
d
TrM (r) = (3) + (4) where (3) is the sum of all terms not involving f and (4) is the sum of
Let dr
other terms with f .
For the first term (3), we have that

|(3)| ≤6 σcT Sc′ Pc σc + 2 σcT Pc Sc′ Pc σc + 4 σcT Pc Pc(2) sc′ + 12 Tr(Sc′ Σ2c ) + 12 σcT Pc(2) sc′
+ 6 Tr(Pc(2) Sc′ Pc + Pc(2) Pc Sc′ ) + 12 Tr(Pc(2) Pc Sc′ Pc )
q
√
√
√
√
≤6 σcT Sc2′ σc n + 2 ksc′ k∞ n + 4 n ksc′ k2 + 12 ksc′ k2 n + 12 kSc′ k2 n
+ 6 diag(Pc Pc(2) )

2

ksc′ k2 + 6 diag(Pc(2) Pc )

≤36n ksc′ k2 .

50

2

ksc′ k2 + 12 diag(Pc Pc(2) Pc )

2

ksc′ k2

For the second term (4), we have that
(4) ≤3 σcT Sc′ Ac (ATc Ac )−1 ∇f + 2 sTc′ Pc(2) Ac (ATc Ac )−1 ∇f
+ 2 σcT Ac (ATc Ac )−1 ATc Sc′ Ac (ATc Ac )−1 ∇f + σcT Ac (ATc Ac )−1 ∇2 f · c′

+ 2 Tr((ATc Ac )−1 ATc Sc′ Ac (ATc Ac )−1 ∇2 f ) + Tr((ATc Ac )−1 ∇3 f [c′ ])
q
q
q
q
(2)
(2)
T
T
T
T
−1
≤3 sc′ Σc sc′ ∇f (Ac Ac ) ∇f + sc′ Pc Pc Pc sc′ ∇f T (ATc Ac )−1 ∇f
q
q
q
q
T
T
T
−1
T
+ 2 σc Pc Sc′ Pc Sc′ Pc σc ∇f (Ac Ac ) ∇f + σc Pc σc c′ ∇2 f (ATc Ac )−1 ∇2 f · c′
+ 2 diag(Ac (ATc Ac )−1 ∇2 f (ATc Ac )−1 ATc ) 2 ksc′ k2 + Tr((ATc Ac )−1 ∇3 f [c′ ])
p
p
√
≤4 ksc′ k2 M1 + 2 ksc′ k∞ nM1 + 3 nM2 ksc′ k2 + M3 ksc′ k2
 p

√
≤ 6 M1 + 3 nM2 + M3 ksc′ k2 .

Therefore,



p
√
d
TrM (s) ≤ O n + M1 + nM2 + M3 ksc′ k2 .
dr

Lemma 58. Let γ be a Hamiltonian curve on ML with ℓ(γ) ≤ ℓ0 . Assume that δ ≤
ζ(t) be the parallel transport of the vector γ ′ (0) on γ(t). Then,

1
1/4 .
36M1

sup kΦ(t)ζ(t)kγ(t) ≤ R3

0≤t≤δ

where



1
2

R3 = O M1

p

3
4

1
4

log n + M1 n δ + M2 n

1
2



.

Proof. By Lemma 52, we have that

−1  T
Aγ (Sγ,µ − 3Σγ + 2Pγ(2) − Sγ ′ Pγ Sγ ′ + Diag(Pγ s2γ ′ ))Aγ − ∇2 f
Φ(t) = ATγ Aγ
= (1) + (2)

where (2) is the last term − ATγ Aγ
For the first term, we have that

−1

∇2 f .

k(1)ζkγ = Pγ (Sγ,µ − 3Σγ + 2Pγ(2) − Sγ ′ Pγ Sγ ′ + Diag(Pγ s2γ ′ ))sγ,ζ
≤ (Sγ,µ − 3Σγ +

2Pγ(2)

− Sγ ′ Pγ Sγ ′ +

Diag(Pγ s2γ ′ ))sγ,ζ

Pγ(2) sγ,ζ

2

2

+ Sγ ′ Pγ Sγ ′ sγ,ζ 2 + Diag(Pγ s2γ ′ )sγ,ζ
≤ kSγ,µ sγ,ζ k2 + 3 kΣγ sγ,ζ k2 + 2
2


≤ ksγ,ζ k∞ ksγ,µ k2 + kσγ k2 + Pγ s2γ ′ 2 + 2 Pγ(2) sγ,ζ + sγ ′ ∞ Sγ ′ sγ,ζ 2
2
 p

√
2
(2)
≤ ksγ,ζ k∞ 2 M1 + n + sγ ′ 4 + 2 Pγ sγ,ζ + sγ ′ ∞ sγ ′ 4 ksγ,ζ k4
2

51

2

Let

where we used that ksγ,µ k22 = kµ(x)k2x ≤ 4M1 , kσγ k2 ≤ n. Now, we note that
Pγ(2) sγ,ζ

2
2

=

X
i


2
X
 (Pγ )2ij (sγ,ζ )j 
j

≤ ksγ,ζ k2∞
= ksγ,ζ k2∞

X
i

X



X
 (Pγ )2ij 
j

σγ2

i



i

= n ksγ,ζ k2∞ .

√
√
1/4
′
log
n
+
M1 δ), we have that
=
O(M
)
and
s
=
O(
γ
1
4
∞
 p

2
k(1)ζkγ ≤ ksγ,ζ k∞ 5 M1 + sγ ′ 4 + sγ ′ ∞ sγ ′ 4 ksγ,ζ k4


p
1/2
1/4
3/4
≤ O M1
ksγ,ζ k∞ + O( log nM1 + M1 δ) ksγ,ζ k4 .

Using also that sγ ′

For the second term, we have that

k(1)ζkγ = Pγ ∇2 f ζ

2
2

≤ ζ T ∇2 f · ζ ≤ M2 ksγ,ζ k2 .

Combining both terms, we have that


p
1/4
3/4
1/2
kΦ(t)ζkγ = O M1
ksγ,ζ k∞ + O( log nM1 + M1 δ) ksγ,ζ k4 + M2 ksγ,ζ k2 .

(7.5)

Now, we bound ksγ,ζ k2 , ksγ,ζ k4 and ksγ,ζ k∞ . (7.1) shows that
−1 T
d
ζ(t) = ATγ Aγ
Aγ Sγ ′ Aγ ζ.
dt

Let wp (t) = kAγ ζ(t)kp . Then, we have that

d
d
wp (t) ≤
Aγ ζ(t)
dt
dt

p

d
ζ(t)
dt

≤ Sγ ′ Aγ ζ(t)

p

+ Aγ

≤ Sγ ′ Aγ ζ(t)

p

+ Pγ Sγ ′ Aγ ζ

≤ sγ ′
d
w2 (t) ≤ 2 sγ ′
For p = 2, we have that dt
1
that t ≤ δ ≤
1/4 , we have

∞

w (t) + Sγ ′ Aγ ζ
∞ p

w2 (t) ≤ e512(
= O(n

where we used that ζ(0) = γ ′ (0) and t ≤
O(n1/2 ).

√

√
log n+ M1 δ)t

1/2

p

.
2

w2 (t). Using that sγ ′

36M1

p

∞

≤ 256

√

log n +

√


M1 δ and

w2 (0)

)

1
1/4
12(v4 +M1 )

52

at the end. Therefore, we have that ksγ,ζ k2 =

For p = 4, we note that
d
w4 (t) ≤ sγ ′ ∞ w4 (t) + sγ ′ 4 ksγ,ζ k4
dt
1/4
≤ 2 sγ ′ 4 w4 (t) = O(M1 w4 (t)).
1
1/4 , we have again
36M1
O(n1/4 ), we have that

that w4 (t) = O(w4 (0)). Since w4 (0) = kAγ ζ(0)k4 =

Since t ≤ δ ≤
kAγ γ ′ (0)k4 =

w4 (t) = O(n1/4 ).

For p = ∞, we note that
d
w∞ (t) ≤ sγ ′ ∞ w∞ (t) + sγ ′ 4 ksγ,ζ k4
dt
1/4
1/4
≤ O(M1 w∞ (t)) + O(M1 n1/4 ).
Again using t ≤ δ ≤

1
1/4 ,
36M1

√
1/4
we have that w∞ (t) ≤ O( log n + M1 n1/4 δ).

Combining our bounds on w2 , w4 , w∞ to (7.5), we get

 1


3
3
1
p
p
1
1
1
1
2
4
kΦ(t)ζkγ = O M1 log n + M1 n 4 δ + O( log nM14 n 4 + M14 n 4 δ) + O M2 n 2
 1

3
p
1
1
2
4
4
2
= O M1 log n + M1 n δ + M2 n
.

7.5

Stability of L2 + L4 + L∞ norm (ℓ1 )

Lemma 59. Given a family of Hamiltonian curves γr (t) on ML with ℓ(γ0 ) ≤ ℓ0 , for δ2 ≤
1 √
√
, we have that
M +M n
1

2



d
1
1/4
ℓ(γr ) ≤ O M1 δ + √
dr
δ log n

1/4
Hence, ℓ1 = O M1 δ +
Proof. For brevity, all

d
dr

√1
δ log n

d
γr (0)
dr


.

are evaluated at r = 0. Since
d
sγ ′
dr r

2

≤ sγr′

∞

Aγr

d
γr
dr

d
′
dr sγr

+ δ Dr γr′ (0)
γr (0)

γr (0)

!

(7.6)

.

d ′
= −Sγr , d γr sγr′ + Aγr dr
γr , we have
dr

+ Aγr
2

d ′
γ
dr r

.
2

For the last term, we note that
Dr γr′ =
Hence, we have
Aγr

d ′
γ
dr r

2

≤ Dr γr′

γr

−1 T
d ′
γr − ATγr Aγr
Aγr Sγr , d γr sγr′ .
dr
dr
+ Sγr , d γr sγ ′
dr

53

2

≤ Dr γr′

γr

+ sγr′

∞

sγr , d γr
dr

2

.

Therefore, we have that
d
sγ ′
dr r

2

d
+ Dr γ ′
γr
dr
γr
p
 d
p
=O
log n + M1 δ
γr
dr

≤ 2 sγr′

∞

γ

+ Dr γr′ .

(7.7)

γr

Since γr is a family of Hamiltonian curves, Lemma 19 shows that
′′

ψ (t) = Γt Φ(t)Γ−1
t ψ(t)
where ψ(t) is the parallel transport of
d
ψ(t) = γr (0) + Dr γr′ (0)t +
dr

d
dr γr (t)

Z

0

t

from γr (t) to γr (0). By Lemma 49, we have that

(t − r)Γr Φ(r)Γ−1
r (

d
γr (0) + Dr γr′ (0)r + E(r))dr
dr

def

d
with kE(r)kF ≤ O(1)∆ and ∆ = dr
γr (0) γr (0) + δ kDr γr′ (0)kγr (0) where we used that kΦ(t)kF,γ =
√
√
1 √
O( M1 + M2 n) (Lemma 56) and that s2 ≤ δ2 ≤ √M +M
.
1
2 n
Therefore, we have that
Z t
d
γr (t)
= ψ(t) γr (0) ≤ ∆ + O(∆) (t − s) Γr Φ(s)Γ−1
dr
r
γr (0)
dr
0
γr (t)

≤ O(∆)

√
√
where we used again kΦ(s)kF,γ = O( M1 + M2 n) and s2 ≤ δ2 ≤
′

Similarly, we have that kDr γr′ (t)kγr (t) = ψ (t)
√1 ,
n

Putting these into (7.7) and using h ≤
d
sγ ′
dr r

We write
ℓ(γr ) = max

0≤t≤δ

sγr′ (t)
n1/2 +

2
1/4
M1

+

sγr′ (t)

1 √
.
M1 +M2 n

≤ O( ∆
δ ).

we have


∆
M1 δ ∆ +
δ
p

1
M1 δ +
=O
∆.
δ

=O
2

γr (0)

√

4

1/4

M1

p

log n +

sγr′ (t)

+√
log n +

p

√∞

M1 δ

+

sγr′ (0)

(7.8)

2

n1/2

+

According to same calculation as (7.7), we can improve the estimate on
d
sγ ′ (0)
dr r

2

≤

p

log n

d
γr (0)
dr

∆
≤ .
δ

sγr′ (0)
n1/4
d
′
dr sγr 2

4

sγ ′ (0) ∞
+ √r
log n

!

.

for t = 0 and get

+ Dr γr′ (0)
γr

(7.9)

54

Using (7.8) and (7.9), we have that
d
d
′
′
dr sγr (t) 2
dr sγr (0) 2
√
max
+
1/4
0≤t≤δ
log n
M1
!
√
M1 δ + 1δ
1
∆
+ √
1/4
δ log n
M1

d
ℓ(γr ) = O
dr
=O
=O

7.6



1/4
M1 δ

1
+ √
δ log n



!

∆.

Mixing Time

Lemma 60. If f (x) = α · φ(x) (logarithmic barrier), then, we have that M1 = n + α2 m, M2 = α
√
and M3 = 2α · n.
−1
Proof. For M1 , we note that (∇f (x))T ATx Ax
∇f (x) = α2 1T Ax (ATx Ax )−1 ATx 1 ≤ α2 m. Hence,
M1 = n + α2 m.
For M2 , it directly follows from the definition.
For M3 , we note that
Tr((ATx Ax )−1 ∇3 f (x)[v]) = −2αTr((ATx Ax )−1 ATx Sx,v Ax )
X
σx,i (sx,v )i .
= −2α
i

Hence, we have
Tr((ATx Ax )−1 ∇3 f (x)[v])

sX
sX
√
2
σx,i
(sx,v )2i ≤ 2α n kvkx .
≤ 2α
i

i

Using Theorem 29, we have the following
Lemma 61. There is a universal constant c > 0 such that if the step size
 1

1
1
1
1
1
1
δ ≤ c · min n− 3 , α− 3 m− 6 n− 6 , α− 2 m− 4 n− 12 ,
then, all the δ conditions for Theorem 30 are satisfied.

Proof. In the previous section, we proved that if δ ≤

1
1/4
36M1

and n is large enough,

1. ℓ0 = 256 (Lemma 54)


1/4
1
(Lemma 59)
2. ℓ1 = O M1 δ + δ√log
n
√
√
3. R1 = O( M1 + M2 n) (Lemma 56)
√
√
4. R2 = O( nM1 + nM1 δ2 +

1/4

M1
δ

+

√

n log n
δ

+

55

√

nM2 + M3 ) (Lemma 57)

3
1√
1
1
5. R3 = O(M12 log n + M14 n 4 δ + M2 n 2 ) (Lemma 58)
1

1

1

1

1

Substituting the value of M1 , M2 and M3 and using that δ . n− 3 , δ . α− 2 n− 3 and δ . α− 2 m− 4 ,
we have that
 1 1

1
1. ℓ1 = O α 2 m 4 δ + δ√log
n
√
√
2. R1 = O( n + α m)

1 √
√
√
4
3. R2 = O(n + α nm + αm +δ n log n )
√
√
3
3
1
4. R3 = O( n log n + α m log n + nδ + α 2 m 4 n 4 δ)

Now, we verify all the δ conditions for Theorem 30.
1
1
1
Using δ . n− 3 and δ . α− 2 m− 4 , we have that δ2 .
Using δ . n

− 13

− 12

and δ . α

m

− 41

1
R1 .

, we have that

1
1
1
δ5 R12 ℓ1 . δ5 (n + α2 m) · (α 2 m 4 δ + √
) ≤ ℓ0 .
δ log n

For the last condition, we note that
p
√
√
1
δ3 R2 + δ2 R3 .nδ3 + α nmδ3 + αm 4 δ2 + n log nδ2
p
p
3
1
3
+ n log nδ2 + α m log nδ2 + nδ3 + α 2 m 4 n 4 δ3
p
p
√
√
3
1
1
3
. αm 4 δ2 + n log nδ2 + α m log nδ2 + nδ3 + α nmδ3 + α 2 m 4 n 4 δ3
p
p
√
3
3
1
. n log nδ2 + α m log nδ2 + nδ3 + α nmδ3 + α 2 m 4 n 4 δ3

where we used that
Therefore, if

√
√
1
αm 4 . (1 + α m) at the end

 1

1
1
1
1
1
1
δ ≤ c · min n− 3 , α− 3 m− 6 n− 6 , α− 2 m− 4 n− 12

for small enough constant, then all the δ conditions for Theorem 30 are satisfied.
Acknowledgement. We thank Ben Cousins for helpful discussions. This work was supported in part
by NSF awards CCF-1563838, CCF-1717349 and CCF-1740551.

References
[1] D. Applegate and R. Kannan. Sampling and integration of near log-concave functions. In
STOC, pages 156–163, 1991.
[2] M. Betancourt. A Conceptual Introduction to Hamiltonian Monte Carlo. ArXiv e-prints,
January 2017.
[3] B. Cousins and S. Vempala. Bypassing KLS: Gaussian cooling and an O ∗ (n3 ) volume algorithm.
In STOC, pages 539–548, 2015.
[4] M. E. Dyer and A. M. Frieze. Computing the volume of a convex body: a case where randomness provably helps. In Proc. of AMS Symposium on Probabilistic Combinatorics and Its
Applications, pages 123–170, 1991.
56

[5] M. E. Dyer, A. M. Frieze, and R. Kannan. A random polynomial time algorithm for approximating the volume of convex bodies. In STOC, pages 375–381, 1989.
[6] M. Girolami, B. Calderhead, and S. A. Chin. Riemannian Manifold Hamiltonian Monte Carlo.
ArXiv e-prints, July 2009.
[7] Mark Girolami and Ben Calderhead. Riemann manifold langevin and hamiltonian monte carlo
methods. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 73(2):123–
214, 2011.
[8] A. T. Kalai and S. Vempala. Simulated annealing for convex optimization. Math. Oper. Res.,
31(2):253–266, 2006.
[9] R. Kannan, L. Lovász, and M. Simonovits. Isoperimetric problems for convex bodies and a
localization lemma. Discrete & Computational Geometry, 13:541–559, 1995.
[10] R. Kannan, L. Lovász, and M. Simonovits. Random walks and an O ∗ (n5 ) volume algorithm
for convex bodies. Random Structures and Algorithms, 11:1–50, 1997.
[11] R. Kannan and H. Narayanan. Random walks on polytopes and an affine interior point method
for linear programming. In STOC, pages 561–570, 2009.
[12] Yin Tat Lee and Aaron Sidford. Path finding methods for linear programming: Solving linear
programs in õ (vrank) iterations and faster algorithms for maximum flow. In Foundations
of Computer Science (FOCS), 2014 IEEE 55th Annual Symposium on, pages 424–433. IEEE,
2014.
[13] Yin Tat Lee and Santosh S. Vempala. Geodesic walks in polytopes. CoRR, abs/1606.04696,
2016.
[14] Yin Tat Lee and Santosh Srinivas Vempala. Eldan’s stochastic localization and the KLS hyperplane conjecture: An improved lower bound for expansion. CoRR, abs/1612.01507, 2016.
[15] Yin Tat Lee and Santosh Srinivas Vempala. Eldan’s stochastic localization and the KLS hyperplane conjecture: An improved lower bound for expansion. In Proc. of IEEE FOCS, 2017.
[16] Yin Tat Lee and Santosh Srinivas Vempala. Geodesic walks in polytopes. In Proceedings of the
49th Annual ACM SIGACT Symposium on Theory of Computing, STOC 2017, Montreal, QC,
Canada, June 19-23, 2017, pages 927–940, 2017.
[17] L. Lovász and M. Simonovits. Mixing rate of Markov chains, an isoperimetric inequality, and
computing the volume. In ROCS, pages 482–491, 1990.
[18] L. Lovász and M. Simonovits. On the randomized complexity of volume and diameter. In Proc.
33rd IEEE Annual Symp. on Found. of Comp. Sci., pages 482–491, 1992.
[19] L. Lovász and M. Simonovits. Random walks in a convex body and an improved volume
algorithm. In Random Structures and Alg., volume 4, pages 359–412, 1993.
[20] L. Lovász and S. Vempala. Fast algorithms for logconcave functions: sampling, rounding,
integration and optimization. In FOCS, pages 57–68, 2006.
[21] L. Lovász and S. Vempala. Hit-and-run from a corner. SIAM J. Computing, 35:985–1005, 2006.
57

[22] L. Lovász and S. Vempala. Simulated annealing in convex bodies and an O ∗ (n4 ) volume
algorithm. J. Comput. Syst. Sci., 72(2):392–417, 2006.
[23] L. Lovász and S. Vempala. The geometry of logconcave functions and sampling algorithms.
Random Struct. Algorithms, 30(3):307–358, 2007.
[24] Oren Mangoubi and Aaron Smith. Rapid mixing of hamiltonian monte carlo on strongly logconcave distributions. arXiv preprint arXiv:1708.07114, 2017.
[25] Radford M. Neal. MCMC using Hamiltonian dynamics. Handbook of Markov Chain Monte
Carlo, 54:113–162, 2010.
[26] R.M. Neal. Bayesian Learning for Neural Networks. Lecture Notes in Statistics. Springer New
York, 1996.
[27] Yurii Nesterov, Arkadii Nemirovskii, and Yinyu Ye. Interior-point polynomial algorithms in
convex programming, volume 13. SIAM, 1994.
[28] Burt Totaro. The curvature of a hessian metric.
15(04):369–391, 2004.

International Journal of Mathematics,

[29] S. Vempala. Geometric random walks: A survey. MSRI Combinatorial and Computational
Geometry, 52:573–612, 2005.

58

A

Matrix ODE

In this section, we prove Lemmas (20) and (21) for the solution of the ODE (3.4), restated below
for convenience.
d2
Ψ(t) = Φ(t)Ψ(t),
dt2
d
Ψ(0) = B,
dt
Ψ(0) = A.
Lemma 62. Consider the matrix ODE (3.4). Let λ = max0≤t≤ℓ kΦ(t)k2 . For any t ≥ 0, we have
that
√
√
kBk
kΨ(t)k2 ≤ kAk2 cosh( λt) + √ 2 sinh( λt).
λ
Proof. Note that
Z

′

t

Ψ(t) = Ψ(0) + tΨ (0) +
(t − s)Ψ′′ (s)ds
0
Z t
= A + tB +
(t − s)Φ(s)Ψ(s)ds.

(A.1)

0

Let a(t) = kΨ(t)k2 , then we have that
a(t) ≤ kAk2 + t kBk2 + λ

Z

t

(t − s)a(s)ds.

0

Let a(t) be the solution of the integral equation
a(t) = kAk2 + t kBk2 + λ

Z

t
0

(t − s)a(s)ds.

By induction, we have that a(t) ≤ a(t) for all t ≥ 0. By taking derivatives on both sides, we have
that
a′′ (t) = λa(t), a(0) = kAk2 , a′ (0) = kBk2 .
Solving these equations, we have
√
√
kBk
kΨ(t)k2 = a(t) ≤ a(t) = kAk2 cosh( λt) + √ 2 sinh( λt)
λ
for all t ≥ 0.
Lemma 63. Consider the matrix ODE (3.4). Let λ = max0≤t≤ℓ kΦ(t)kF . For any 0 ≤ t ≤
have that


t3
2
kΨ(t) − A − BtkF ≤ λ t kAk2 + kBk2 .
5
In particular, this shows that

Ψ(t) = A + Bt +

with kE(s)kF ≤ λ s2 kAk2 +

s3
5


kBk2 .

Z

t
0

(t − s)Φ(s)(A + Bs + E(s))ds

59

√1 ,
λ

we

Proof. Recall from (A.1) that
Ψ(t) = A + tB +

Z

t
0

(A.2)

(t − s)Φ(s)Ψ(s)ds.

Let E(t) = Ψ(t) − (A + tB). Using Lemma 20, we have that
kE(t)kF =

Z

≤λ

t
0

Z

Z

(t − s)Φ(s)Ψ(s)ds

F

t

(t − s) kΨ(s)k2 ds

0

√
√
kBk
(t − s) kAk2 cosh( λs) + √ 2 sinh( λs)ds
λ

0
√
√
√
kBk2
= λ kAk2 (cosh( λt) − 1) + √ (sinh( λt) − λt) .
λ
≤λ

t

√
√
√
Since 0 ≤ t ≤ √1λ , we have that cosh( λt) − 1 ≤ λt2 and sinh( λt) − λt ≤
the result.
The last equality follows again from (A.2)

λ3/2 t3
5 .

This gives

Next, we have an elementary lemma about the determinant.
Lemma 64. Suppose that E is a matrix (not necessarily symmetric) with kEk2 ≤ 41 , we have
|log det(I + E) − TrE| ≤ kEk2F .
Proof. Let f (t) = log det(I + tE). Then, by Jacobi’s formula, we have

f ′ (t) = Tr (I + tE)−1 E ,

f ′′ (t) = −Tr((I + tE)−1 E(I + tE)−1 E).

Since kEk2 ≤ 14 , we have that (I + tE)−1
f ′′ (t)

=
≤

2

≤

4
3

and hence

Tr((I + tE)−1 E(I + tE)−1 E)
T
Tr(E T (I + tE)−1 (I + tE)−1 E)

≤ 2 Tr(E T E) = 2 kEk2F .

The result follows from
Z

1

(1 − s)f ′′ (s)ds
f (1) = f (0) + f (0) +
0
Z 1
(1 − s)f ′′ (s)ds.
= Tr(E) +
′

0

60

B

Concentration

Lemma 65 ([13, Ver 3, Lemma 90]). For p ≥ 1, we have
p 


!1/p
 2
p+1 X
p/2
2 Γ( 2 )
t
p
p




√
+ kAk2→p t
Px∼N (0,I) kAxkp ≤
kai k2
.
≤ 1 − exp −
π
2
i

In particular, we have

and



Px∼N (0,I) kAxk44 ≤  3




Px∼N (0,I) kAxk22 ≤ 

C

X
i

X
i

kai k42

kai k22

!1/4

!1/2

4 

 2
t


+ kAk2→4 t
≤ 1 − exp −
2
2 

 2
t
.
+ kAk2→2 t  ≤ 1 − exp −
2

Calculus

Proof of Fact 48. Recall Definition 47 and write
dSγ−1
dAγ
=
A
dt
dt
dSγ −1
Sγ A
= −Sγ−1
dt 

d(Aγ − b)
−1
= −Sγ Diag
Aγ
dt

= −Diag Sγ−1 Aγ ′ Aγ
= −Diag(Aγ γ ′ )Aγ = −Sγ ′ Aγ .

For the second, using the first,
dAγ (ATγ Aγ )−1 ATγ
dPγ
=
dt
dt
d(ATγ Aγ )−1 T
dAγ T
dAγ
=
(Aγ Aγ )−1 ATγ + Aγ (ATγ Aγ )−1
+ Aγ
Aγ
dt
dt
dt
d(ATγ Aγ ) T
(Aγ Aγ )−1 ATγ
= −Sγ ′ Pγ − Pγ Sγ ′ − Aγ (ATγ Aγ )−1
dt

= −Sγ ′ Pγ − Pγ Sγ ′ + 2Aγ (ATγ Aγ )−1 ATγ Sγ ′ Aγ (ATγ Aγ )−1 ATγ
= −Sγ ′ Pγ − Pγ Sγ ′ + 2Pγ Sγ ′ Pγ .

And for the last,
dSγ ′
= Diag
dt



dAγ ′
γ + Aγ γ ′′
dt



= Diag(−Sγ ′ Aγ γ ′ + Aγ γ ′′ ) = −Sγ2′ + Sγ ′′ .

61

D

Basic definitions of Riemannian geometry

Here we recall basic notions of Riemannian geometry. One can think of a manifold M as a ndimensional “surface” in Rk for some k ≥ n.
1. Tangent space Tp M : For any point p, the tangent space Tp M of M at point p is a linear
subspace of Rk of dimension n. Intuitively, Tp M is the vector space of possible directions that
are tangential to the manifold at x. Equivalently, it can be thought as the first-order linear
d
c(t) is tangent
approximation of the manifold M at p. For any curve c on M , the direction dt
′
to M and hence lies in Tc(t) M . When it is clear from context, we define c (t) = dc
dt (t). For any
n
n
open subset M of R , we can identify Tp M with R because all directions can be realized by
derivatives of some curves in Rn .
2. Riemannian metric: For any v, u ∈ Tp M , the inner product (Riemannianq
metric) at p is given

by hv, uip and this allows us to define the norm of a vector kvkp = hv, vip . We call a
manifold a Riemannian manifold if it is equipped with a Riemannian metric. When it is clear
from context, we define hv, ui = hv, uip . In Rn , hv, uip is the usual ℓ2 inner product.

3. Differential (Pushforward) d: Given a function f from a manifold M to a manifold N , we
define df (x) as the linear map from Tx M to Tf (x) N such that
df (x)(c′ (0)) = (f ◦ c)′ (0)
for any curve c on M starting at x = c(0). When M and N are Euclidean spaces, df (x) is the
Jacobian of f at x. We can think of pushforward as a manifold Jacobian, i.e., the first-order
approximation of a map from a manifold to a manifold.
4. Hessian manifold: We call M a Hessian manifold (induced by φ) if M is an open subset of Rn
with the Riemannian metric at any point p ∈ M defined by
hv, uip = v T ∇2 φ(p)u
where v, u ∈ Tp M and φ is a smooth convex function on M .
5. Length: For any curve c : [0, 1] → M , we define its length by
L(c) =

Z

1

0

d
c(t)
dt

dt.
c(t)

6. Distance: For any x, y ∈ M , we define d(x, y) be the infimum of the lengths of all paths
connecting x and y. In Rn , d(x, y) = kx − yk2 .
7. Geodesic: We call a curve γ(t) : [a, b] → M a geodesic if it satisfies both of the following
conditions:
(a) The curve γ(t) is parameterized with constant speed. Namely,
for t ∈ [a, b].

d
dt γ(t) γ(t)

is constant

(b) The curve is the locally shortest length curve between γ(a) and γ(b). Namely, for any
family of curve c(t, s) with c(t, 0) = γ(t) and c(a, s) = γ(a) and c(b, s) = γ(b), we have
Rb d
d
c(t, s) c(t,s) dt = 0.
that ds
s=0 a dt
62

Note that, if γ(t) is a geodesic, then γ(αt) is a geodesic for any α. Intuitively, geodesics are
local shortest paths. In Rn , geodesics are straight lines.
8. Exponential map: The map expp : Tp M → M is defined as
expp (v) = γv (1)
where γv is the unique geodesic starting at p with initial velocity γv′ (0) equal to v. The
exponential map takes a straight line tv ∈ Tp M to a geodesic γtv (1) = γv (t) ∈ M . Note that
expp maps v and tv to points on the same geodesic. Intuitively, the exponential map can be
thought as point-vector addition in a manifold. In Rn , we have expp (v) = p + v.
9. Parallel transport: Given any geodesic c(t) and a vector v such that hv, c′ (0)ic(0) = 0, we
define the parallel transport Γ of v along c(t) by the following process: Take h to be infinitesimally small and v0 = v. For i = 1, 2, · · · , 1/h, we let vih be the vector orthogonal to c′ (ih)
that minimizes the distance on the manifold between expc(ih) (hvih ) and expc((i−1)h) (hv(i−1)h ).
Intuitively, the parallel transport finds the vectors on the curve such that their end points are
closest to the end points of v. For general vector v ∈ Tc′ (0) , we write v = αc′ (0) + w and we
define the parallel transport of v along c(t) is the sum of αc′ (t) and the parallel transport of
w along c(t). For a non-geodesic curve, see the definition in Fact 66.
10. Orthonormal frame: Given vector fields v1 , v2 , · · · , vn on a subset of M , we call {vi }ni=1 is an
orthonormal frame if hvi , vj ix = δij for all x. Given a curve c(t) and an orthonormal frame at
c(0), we can extend it on the whole curve by parallel transport and it remains orthonormal
on the whole curve.
11. Directional derivatives and the Levi-Civita connection: For a vector v ∈ Tp M and a vector
field u in a neighborhood of p, let γv be the unique geodesic starting at p with initial velocity
γv′ (0) = v. Define
u(h) − u(0)
∇v u = lim
h→0
h
where u(h) ∈ Tp M is the parallel transport of u(γ(h)) from γ(h) to γ(0). Intuitively, LeviCivita connection is the directional derivative of u along direction v, taking the metric into
d
u(x + tv). When u is defined on a curve
account. In particular, for Rn , we have ∇v u(x) = dt
d
d
n
c, we define Dt u = ∇c′ (t) u. In R , we have Dt u(γ(t)) = dt
u(γ(t)). We reserve dt
for the usual
derivative with Euclidean coordinates.
We list some basic facts about the definitions introduced above that are useful for computation and
intuition.
Fact 66. Given a manifold M , a curve c(t) ∈ M , a vector v and vector fields u, w on M , we have
the following:
1. (alternative definition of parallel transport) v(t) is the parallel transport of v along c(t) if and
only if ∇c′ (t) v(t) = 0.
2. (alternative definition of geodesic) c is a geodesic if and only if ∇c′ (t) c′ (t) = 0.
3. (linearity) ∇v (u + w) = ∇v u + ∇v w.
4. (product rule) For any scalar-valued function f, ∇v (f · u) =
63

∂f
∂v u

+ f · ∇v u.

5. (metric preserving)

d
dt

hu, wic(t) = hDt u, wic(t) + hu, Dt wic(t) .

∂c
6. (torsion free-ness) For any map c(t, s) from a subset of R2 to M , we have that Ds ∂c
∂t = Dt ∂s
where Ds = ∇ ∂c and Dt = ∇ ∂c .
∂s

∂t

7. (alternative definition of Levi-Civita connection) ∇v u is the unique linear mapping from the
product of vector and vector field to vector field that satisfies (3), (4), (5) and (6).

D.1

Curvature

Roughly speaking, curvature measures the amount by which a manifold deviates from Euclidean
space. Given vector u, v ∈ Tp M , in this section, we define uv be the point obtained from moving
from p along direction u with distance kukp (using geodesic), then moving along direction “v” with
distance kvkp where “v” is the parallel transport of v along the path u. In Rn , uv is exactly p + u + v
and hence uv = vu, namely, parallelograms close up. For a manifold, parallelograms almost close
up, namely, d(uv, vu) = o(kuk kvk). This property is called being torsion-free.
1. Riemann curvature tensor: Three-dimensional parallelepipeds might not close up, and the
curvature tensor measures how far they are from closing up. Given vector u, v, w ∈ Tp M , we
define uvw as the point obtained by moving from uv along direction “w” for distance kwkp
where “w” is the parallel transport of w along the path uv. In a manifold, parallelepipeds do
not close up and the Riemann curvature tensor how much uvw deviates from vuw. Formally,
for vector fields v, w, we define τv w be the parallel transport of w along the vector field v for
one unit of time. Given vector field v, w, u, we define the Riemann curvature tensor by
R(u, v)w =

d d −1 −1
τ τ τsu τtv w
ds dt su tv

.

(D.1)

t,s=0

Riemann curvature tensor is a tensor, namely, R(u, v)w at point p depends only on u(p), v(p)
and w(p).
2. Ricci curvature: Given a vector v ∈ Tp M , the Ricci curvature Ric(v) measures if the geodesics
starting around p in direction v converge together. Positive Ricci curvature indicates the
geodesics converge while negative curvature indicates they diverge. Let S(0) be a small shape
around p and S(t) be the set of point obtained by moving S(0) along geodesics in the direction
v for t units of time. Then,
volS(t) = volS(0)(1 −
Formally, we define
Ric(v) =

t2
Ric(v) + smaller terms).
2

X
ui

(D.2)

hR(v, ui )ui , vi

where ui is an orthonormal basis of Tp M . Equivalently, we have Ric(v) = Eu∼N (0,I) hR(v, u)u, vi.
kvk2 .
For Rn , Ric(v) = 0. For a sphere in n + 1 dimension with radius r, Ric(v) = n−1
r2
Fact 67 (Alternative definition of Riemann curvature tensor). Given any M -valued function c(t, s),
∂c
we have vector fields ∂c
∂t and ∂s on M . Then, for any vector field z,
R(

∂c ∂c
, )z = ∇ ∂c ∇ ∂c z − ∇ ∂c ∇ ∂c z.
∂t
∂s
∂s
∂t
∂t ∂s

Equivalently, we write R(∂t c, ∂s c)z = Dt Ds z − Ds Dt z.
64

Fact 68. Given vector fields v, u, w, z on M ,
hR(v, u)w, zi = hR(w, z)v, ui = − hR(u, v)w, zi = − hR(v, u)z, wi .

D.2

Hessian manifolds
2

Recall that a manifold is called Hessian if it is a subset of Rn and its metric is given by gij = ∂x∂i ∂xj φ
for some smooth convex function φ. We let gij be entries of the inverse matrix of gij . For example,
P
2
∂3
we have j gij gjk = δik . We use φij to denote ∂x∂i ∂xj φ and φijk to denote ∂xi ∂x
j ∂xk φ.
Since a Hessian manifold is a subset of Euclidean space, we identify tangent spaces Tp M by
Euclidean coordinates. The following lemma gives formulas for the Levi-Civita connection and
curvature under Euclidean coordinates.
Lemma 69 ([28]). Given a Hessian manifold M , vector fields v, u, w, z on M , we have the following:
1. (Levi-Civita connection) ∇v u =
and the Christoffel symbol

P

∂uk
ik vi ∂xi ek

Γkij =

+

P

k
ijk vi uj Γij ek

where ek are coordinate vectors

1 X kl
g φijl .
2
l

2. (Riemann curvature tensor) hR(u, v)w, zi =
Rklij =

3. (Ricci curvature) Ric(v) =

1
4

P

P

ijlk

Rklij ui vj wl zk where

1 X pq
g (φjkp φilq − φikp φjlq ) .
4 pq

ijlkpq g

pq g jl

(φjkpφilq − φikp φjlq ) vi vk .

65

