Accelerating Hamiltonian Monte Carlo via Chebyshev Integration Time
Authors: Jun-Kun Wang
Andre Wibisono
Venue: ICLR 2023 poster (2023)

Abstract: Hamiltonian Monte Carlo (HMC) is a popular method in sampling. While there are quite a few works of studying this method on various aspects, an interesting question is how to choose its integration time to achieve acceleration. In this work, we consider accelerating the process of sampling from a distribution $\pi(x) \propto \exp(-f(x))$ via HMC via time-varying integration time. When the potential $f$ is $L$-smooth and $m$-strongly convex, i.e. for sampling from a log-smooth and strongly log-concave target distribution $\pi$, it is known that under a constant integration time, the number of iterations that ideal HMC takes to get an $\epsilon$ Wasserstein-2 distance to the target $\pi$ is $O( \kappa \log \frac{1}{\epsilon} )$, where $\kappa := \frac{L}{m}$ is the condition number. We propose a scheme of time-varying integration time based on the roots of Chebyshev polynomials. We show that in the case of quadratic potential $f$, i.e. when the target $\pi$ is a Gaussian distribution, ideal HMC with this choice of integration time only takes $O( \sqrt{\kappa} \log \frac{1}{\epsilon} )$ number of iterations to reach Wasserstein-2 distance less than $\epsilon$; this improvement on the dependence on condition number is akin to acceleration in optimization. The design and analysis of HMC with the proposed integration time is built on the tools of Chebyshev polynomials. Experiments find the advantage of adopting our scheme of time-varying integration time even for sampling from distributions with smooth strongly convex potentials that are not quadratic. 


Published as a conference paper at ICLR 2023
ACCELERATING HAMILTONIAN MONTE CARLO VIA CHEBYSHEV INTEGRATION TIME
Jun-Kun Wang and Andre Wibisono
Department of Computer Science, Yale University (jun-kun. wang, andre. wiwi aosoooooolllleeee
ABSTRACT
Hamiltonian Monte Carlo (HMC) is a popular method in sampling. While there are quite a few works of studying this method on various aspects, an interesting question is how to choose its integration time to achieve acceleration. In this work, we consider accelerating the process of sampling from distribution T(x) exp(-f(z)) via HMC via time-varying integration time. When the potential f is L-smooth and m-strongly convex, i.c. for sampling from a log-smooth and strongly log-concave target distribution т, it is known that under a constant integration time, the number of iterations that ideal HMC takes to get an e Wasserstein-2 distance to the target is Ok log 2), where := L is the condition number. We propose scheme of time- varying integration time based on the roots of Chebyshev polynomials. We show that in the case of quadratic potential i.e. when the target t is a Gaussian distribution, ideal HMC with this choice of integration time only takes OlVR log !) number of iterations to reach Wasserstein-2 distance less than e; this improvement on the dependence on condition number is akin to acceleration optimization. The design and analysis of HMC with the proposed integration time is built on the tools of Chebyshev polynomials. Experiments find the advantage of adopting our scheme of varring integration time even for sampling from distributions with smooth strongly convex potentials that are not quadratic.
INTRODUCTION
Markov chain Monte Carlo (MCMC) algorithms are fundamental techniques for sampling from
probability distributions, which is a task that naturally arises in statistics (Duane et al.)
1987: Girolami
Calderhead T) optimization (Flaxman et al. [2055] Duchi et al. 2012 al. 2017). machine
learning and others (Wenzel et al., 2020] Salakhutdinoov & Mnih 2008) Koller & Friedman 2009
Welling & Tch 201 1, Among all the MCMC algorithms, the most popular ones perhaps are Langevin
methods (Li et al. 2022 Dalalyan] 2017 Durmus et al. 2019 Vempala & Wibisono 2019 Lee
et al. 202 Chewi et al. 2020) and Hamiltonian Monte Carlo (HMC) (Neal [012]] Betancourt
2017 Hoffman & Gelman 2014 Levy et al.) 2018). For the former, recently there have been
sequence of works leveraging some techniques in optimization to design Langevin methods, which
include borrowing the idea of momentum methods like Nesterov acceleration (Nesterov (013) to
design fast methods, e.g., (Ma et al. 2022 Dalalyan & Riou-Durand] (020).. Specifically, Ma et al. (202I) show that for sampling from distributions satisfying the log-Sobolev inequality, under-damped
Langevin improves the iteration complexity of over-damped Langevin from O( O(v 4) where is the dimension and is the error in KL divergence, though whether their result has an optimal dependency on the condition number is not clear. On the other hand, compared Langevin methods, the connection between HMCs and techniques in optimization seems rather loose. Moreover, to our knowledge, little known about how to accelerate HMCs with a provable acceleration guarantee for converging to a target distribution. Specifically, [Chen & Vempala (2019) show that for sampling from strongly log-concave distributions, the iteration complexity of ideal HMC is O log 4), and Vishnoi (2021) shows the same rate of ideal HMC when the potential is strongly convex quadratic in nice tutorial. In contrast, there are a few methods that exhibit acceleration when minimizing strongly convex quadratic functions in optimization. For example, while Heavy Ball (Polyak) 1964) does not have an accelerated linear rate globally for minimizing general smooth strongly convex functions, does show acceleration when minimizing strongly convex quadratic functions (Wang al. 2020

---

Published as a conference paper at ICLR 2023
Algorithm 1: IDEAL HMC
1: Require: an initial point 2b Ie number of iterations K, and a scheme of integration time {17 (K)
2: for k = 1 to do
3: Sample velocity È ~ N(0, Ta).
4: Set (Tk, vx)) = HMC do (T)-1,
5: end for
2021 2022). This observation makes us wonder whether one can get an accelerated linear rate of ideal HMC for sampling, i.e., Olv hi log 2), akin acceleration in optimization.
We answer this question affirmatively, at least in the Gaussian case. We propose a time-varying integration time for HMC, and show that ideal HMC with this time-varying integration time exhibits acceleration when the potential is a strongly convex quadratic (i.e. the target и is a Gaussian), compared to what is established in Chen & Vempala (2019) and Vishnoi (2021) for using constant integration time. Our proposed time-varying integration time at each iteration of HMC depends on the total number of iterations K, the current iteration index k. the strong convexity constant m, and the smoothness constant L of the potential; therefore, the integration time at each iteration is simple to compute and set before executing HMC. Our proposed integration time based on the roots of Chebysev polynomials, which we will describe in details in the next section. In optimization, Chebyshev polynomials have been used to help design accelerated algorithms for minimizing strongly convex quadratic functions, i.e., Chebyshev iteration (see e.g., Section 2.3 in d' A'spremnnt al. (2021)). Our result of accelerating HMC via using the proposed Chebyshev integration time can be viewed the sampling counterpart of acceleration from optimization. Interestingly. for minimizing strongly convex quadratic functions, acceleration of vanilla gradient descent can be achieved via a scheme of step sizes that is based our work is inspired by a nice blog article by Pedregosa (202T). Hence, our acceleration result of HMC can also be viewed as a counterpart in this sense. In addition to our theoretical findings, we conduct experiments of sampling from a Gaussian as well as sampling from distributions whose potentials are not quadratics, which include sampling from mixture two Gaussians, Bayesian logistic regression, and sampling from a hard distribution that was proposed in Lee et al. (2021a) for establishing some lower-bound results of certain Metropolized sampling methods. Experimental results show that our proposed time-varying integration time also leads to better performance compared to using the constant integration time of Chen & Vempala (2019) and Vishnoi (202T) l for sampling from the distributions whose potential functions not quadratic. We conjecture that our proposed timevaryiigg integration time also helps accelerate HMC for sampling from omosthh and strongly log-concave distributions, and leave the analysis of such cases for future work.
Chebyshev polynomial, see e.g., Agarwal al. (2021). and
2
PRELIMINARIES
2.1 HAMILTONIAN MONTE CARLO (HMC)
Suppose we want to sample from a target probability distribution v(z) x exp(-f(z)) on Rd, where f: Rd + IR is a continuous function which we refer to as the potential.
Denote E &d the position and E R the velocity of particle. In this paper, consider the
standard Hamiltonian of the particle (Chen & Vempala 2019 Neal 2012), which defined as
Hzz,v) f() + 111011
while we refer the readers to Girolami Calderhead (2011): Hirt al. (22 Brofos & Lederman (2021) and the references therein for other notions of the Hamiltonian. The Hamiltonian flow generated by H the flow of the particle which evolves according to the following differential equations:
dx
aH and dr aH
dt Dv For the standard Hamiltonian defined in T.. the Hamiltonian flow becomes
dt dr
dx = v and
du
dt = Vf(x).
(2)
dt
2

---

Published as a conference paper at ICLR 2023
We will write (zz, 0i) HMC (2o, vo) the position and the velocity of the Hamiltonian flow after integration time t starting from (Øo. vo). There are many important properties of the Hamiltonian flow including that the Hamiltonian is conserved along the flow, the vector field associated with the flow is divergence free, and the Hamiltonian dynamic is time reversible, see e.g. Section 3 in Vishnoi (202T).
The Ideal HMC algorithm (see Algorithm D proceeds as follows: in each iteration k, sample an initial velocity from the normal distribution, and then flow following the Hamiltonian flow with a pre-specified integration time ?)k- It is well-known that ideal HMC preserves the target density (a) x exp(-f(x))
see e.g., Theorem 5.1 in Vishnoi (202I). Furthermore, in each iteration, HMC brings the density of the iterates IĄk ~ Pk: closer to the target π. However, the Hamiltonian flow (M,, 00)) is in general difficult to simulate exactly, except for some special potentials. In practice, the Verlet integrator is commonly used to approximate the flow and Metropolis-Hasting filter applied to correct the induced bias arises from the use of the integrator (Tripuraneni et al.)
2017 Brofos & Lederman 2021 Hoffman et al.) 2021 Lee et al. 202 Chen et al. 2020)) In recent years, there have been some progress on showing some rigorous theoretical guarantees of HMCs for converging to target distribution, e.g.. Chen al. 2020); Durmus et al. (2017); (Bo--Rabee
& Eberle (2021); Mangoubi & Smith (2019) 2021); Mangoubs & Vishnoi (20 (2018). There are also
other variants of HMCS proposed in the literature, e-g.. Riou- Kiou-Durand & Vogrinc (2022); Bouu-Rbbee
& Sanz-Serna] y Zou & Gu (202T); Steeg & Galstyan (2021) y Hoffman & Gelman (014)):
Tripuraneni al. (2017); Chen al. (2014). 4. name just few.
Recall that the Wasserssnninn distance between probability distributions and My is
W2(v1, a.yET(PL,"a)
inf
E [||= y1
where [(21,12) represents the set of all couplings of 11 and V2.
2.2 ANALYSIS OF HMC IN QUADRATIC CASE WITH CONSTANT INTEGRATION TIME
In the following, we replicate the analysis of ideal HMC with a constant integration time for quadratic potentials (Vishnoi 2021), which provides the necessary ingredients for introducing method in the next section. Specifically, we consider the following quadratic potential:
= (x)==>; >ya3. where 0 < m S >3 S I.,
which means the target density is the Gaussian distribution π (O, where A the diagonal matrix whose juh diagonal entry is Ay. We note for general Gaussian target ((,,,) M) for some E R&d and Σ > O, we can shift and rotate the coordinates to make μ = 0 and Σ a diagonal matrix, and our analysis below applies. So without loss of generality, we may assume the quadratic potential is separable, as in 33).
In this quadratic case, the Hamiltonian flow (2) becomes a linear system of differential equations, and we have an exact solution given by sinusoidal functions, which are
xe [i] = cos (V2g xols]] + /2X; sin (V2Xzl) voli),
(4)
vr = - 2Xy 22 sin (22xyt) zoli] + cos (V2Xtt) voli).
In particular, we recall the following result on the deviation between two co-evolving particles with the same initial velocity.
Lemma 1. Vishnoi 2027) Let No. No E Rd. Consider the following coupling: (It, v.) HMC, (Zo. è) and (yt, u.) = HMC+(yb. {) for some & E Rdd Then for all 2 0 and for all j E [d], it holds that
2 1 ve [i] cos (22x)t) (xoli) yo(1))
Using Lemma I we can derive the convergence rate of ideal HMC for the quadratic potential as follows.
3

---

Published as a conference paper at ICLR 2023
Lemma 2. Vishnoi 2021) Let exp(-f) NO be the target distribution, where f(x) is
defined on (3). Let PK be the distribution of и generated by Algorithm a the final iteration K. Then for any Po and any K 1, we have
(2KK, T) I maxjeld I cos (V2XyTh
W2(po
We replicate the proof of Lemma and Lemma 2 in Appendix for the reader's convenience.
Vishnoi 2021
shows that by choosing
(Constant integration time)
(K)
(5)
2
22
one has that cos 2X
(K) < - (T) for all the iterations k E [K] and dimensions j [d].
Hence, by Lemma] 2 the distance satisfies
((1 (())) W W2(Po,
W2(PK, x)
after K iterations ideal HMC with the constant integration time. On the other hand, for general smooth strongly convex potentials f(-). Chen & Vemppla (2019) show the same convergence rate
1 (T) of HMC using a constant integration time ` [K] VT. where c > 0 is a universal constant. Therefore, under the constant integration time, HMC needs log !) iterations to reach error (2(Kk, TT < where = L is condition number. Furthermore, they also show that the relaxation time of ideal HMC with constant infegration time is Q(x) for the Gaussian case.
2.3
CHEBYSHEV POLYNOMIALS
We denote (-) the degree-K Chebyshev polynomial of the first kind, which is defined by:
cos(( arccos(z))
if x E [-1,1], if x > if x <
Φx(x) = cosh(K (K arccosh(z))
(6)
(-1)" cosh(K arccosh(z))
Our proposed integration time is built ocaled-and-shifter Chebyshev polynomial, defined as:
Φ (I(A))
4 k (A) =
k (h())))
where h(.) the mapping ((A)) == L+m-22 Observe that the mapping I--) maps all À E [m, L
L-
into the interval [-1, 1]. The roots of the degree- K scaled-and-shifted Chebyshev polynomial 4x (A) are
(Chebyshev roots)
Tk =
(K)
L
cos
(K-
2) m
(8)
K
where k = 1.2..... K, i.e., T5 (rk = 0. We now recall the following key result regarding the scaled-and-shifted Chebyshev polynomial & Lemma 3. (e.g., Section 2.3 in dd'Aspremont et al. (202)) For any positive integer K, we have
maxxe(m,LL c (A) <
=0 -
(9)
VT+
The proof of Lemma 3 is in Appendix B
CHEBYSHEV INTEGRATION TIME
We are now ready introduce our scheme of time-varying integration time. Let K be the pre-specified total number of iterations of HMC. Our proposed method will first permute the array [1, 2, K] before executing HMC for K iterations. Denote (k) the kn, element of the array 21, K] after an arbitrary permutation a. Then, we propose to set the integration time of HMC at iteration k. i.c.,
(K) T/s
as follows:
4

---

Published as a conference paper at ICLR 2023
10"
10-4
0.8
product
pi2 s4gtt)))
10
0.6
10
0.4
10-21
10
Chebyshev time Const time
0.0
20
40
60
80
100
120
140 160
180
20
40
60
80
100
Iteration
Figure 1:
Left: Set K
400,
1 and L
100.
The green solid line (Cheby-
shev integration time
(10) on the subfigure represents maxxg (m,m L... II
= COS
2Xm
IT
101000
v.s. & while blue dash (Constant integration time
represents
max (m. 1.........
1 cos (VeXn
- cos [
v.s. k. Since the cosine product con-
trols the convergence rate the W22 distance by l.mma 2 this confirms the acceleration via using the proposed scheme of Chebyshev integration over the constant integration time (Chen & Vempala 2019 Vishnoi 2021).
cos(
Right: (xx)
V.S. J.
(Chebyshev integration time)
(K) 1
1
(10)
(2
(M)
We note the usage of the permutation ø is not needed in our analysis below: however, it seems to help improve performance practice. Specifically, though the guarantees of HMC at the final iteration K provided in Theorem [] and Lemma below is the same regardless of the permutation, the progress of HMC varies under different permutations of the integration time, which is why we recommend an arbitrary permutation of the integration time in practice.
Our main result is the following improved convergence rate HMC under the Chebyshev integration time, for quadratic potentials.
Theorem 1. Denote the target distribution с cxp(-f(x)) = .N(0 where f(x) is defined on 3 and denote the condition number k =:= L. Let be the distribution of x generated by Algorithm Mat the final iteration K. Then, have
Wa(PK, m) s 1- (1-2
w2( =0 (0
W W2(Po,
VL
v
Consequently, the total number of iterations K such that the Wasserstein-2 distance satisfies We(pk.T) is Ses (VElog
Theorem shows an accelerated linear rate 1 - θ (vR using Chebyshev integration time, and
hence improves the previous result of 1 θ (2) as discussed above. The proof of Theorem relies on the following lemma, which upper-bounde the cosine products that appear in the bound of the W2 distance in Lemma 2 by the sedled-and-shifee Chebyshev polynomial & (A) (7.
Lemma 4. Denote [P
K
COS
Suppose À E [m. L). Then, we have for
bock)
any positive integer K,
(A)
xx (A)
(11)
The proof of Lemma maxxe|m.L] I cos cos
is available in Appendix
C Figure compares the cosine product
v V2
(K) in Lemma of using the proposed integration time and that

---

for 0, Sk do Published as a conference paper at ICLR 2023
Algorithm 2: HMC WITH CHEBYSHEV INTEGRATION TIME
and f() L-soooth m-sttrongly convex. 2: Require: number iterations and the step size the leapfrog steps
Given: a potential f(), where ((r) &
3: Define rk (KK L L Lm cos
(k- 1/== for k I,...,K. K.
K
4: Arbitrarily permute the array [1. 2........ ................................ K]. Denote (Ik) kn element of array after permutation.
5: for l 1.2...... .K do
6: Sample velocity &a ~ (0,, Ia).
7: Set integration time f) 2
(K)
(K)
(K)
8: Set the number of leapfrog steps Sa + (Ta, Bb) - (I'--1, (x)
9
9:
% Leapfrog steps
10:
11:
1 1 F. (F.e)
Jott F. 10 s
- [Vf(s.+1);
12:
end for
% Metropolis filter
13:
Compute the acceptance ratio ag min (1, expl H ((0
-xpi- 18,
expl
14:
Draw { Uniform(0, (0, 1).
If On then
I& 4 S.S.
17:
Else
18:
I6-1-
19:
end for
of using the constant integration time, which illustrates acceleration via the proposed Chebyshev integration time.
We now provide the proof of Theorem
Proof. (of Theorem T From Lemma
we have
W2(PK, T) I maxje[4] _ (eaccos(222n (v2Xm x (2(
(12)
We upper-hound the cosine product of any [d] as,
TTK cos (V2Xmm
(K)
(a)
II
cos
2
xx(())] S 2 1-
(0)
VT+./m
mm
(13)
where (a) is due to the of Chebyshev integration time 10 (b) is by Lemma and (c) is by
Lemma Combining (1) and (13) leads to the result.
HMC with Chebyshev Integration Time for General Distributions To sample from general strongly log-concave distributions, propose Algorithm 2 which adopts the Verlet integrator (a.ka. the leapfrog integrator) to simulate the Hamiltonian flow HMC, {) and uses Metropolis filter correct the bias. It is noted that the number of leapfrog steps S₃ in each iteration & is equal to (K) divided by the step size used in the leapfrog steps. More precisely, we the integration time 7
have Sg- in iteration A of HMC.
4 EXPERIMENTS
We now evaluate HMC with the proposed Chebyshev integration time (Algorithm 2p and HMC with the constant integration time (Algorithm 2 with line replaced by the constant integration time (5) in several tasks. For all the tasks in the experiments, the total number of iterations of HMCs is set to 10, 000, and hence we collect K = 10, 000 samples along the trajectory. For the step size in the leapfrog steps, let {0.001, 0.005, 0.01, 0.05). To evaluate the methods,

---

Published as a conference paper at ICLR 2023
Table 1: Ideal HMC with K 10,000 000 iterations for sampling from a PouussianNV 2). where μ
and
0 Here, Cheby. (W/) ideal HMC with a arbitrary permutation the Chebyshev integration
Σ 0 100 time, while Cheby. (W/O) ideal HMC without a permutation; and Const. refers using the constant integration time
Method
Mcan ESS
Min ESS
Cheby. (W7)
10.399008 11 +: 347. 25021
7172.503336 1: 257-21244
Cheby, (W/O) Coust.
10197.09964 = 276.94894
7043.55225 + 284.78037
7692.00382 196288
5533 26519
+ 213.31943
compute effective sample size (ESS), which common performance metric of HMCs (Girolami
& Calderhead 2011
Brofos & Lederman, [021]] Hirt et al. 2021 Riou-Durand & Vogrine [022]
Hoffman et al., 2021
Hoffman Gelman 2014] Steeg & Galstyan 2021). by using the toolkit ArViz
(Kumar et al. 2019). The ESS of sequence of N dependent samples is computed based on the
autocorrelations within the sequence different lags: ESS := N/(1 2 Σx 7(k)), where 7(k) (k) an
estimate of the autocorrelation at lag k. We consider 4 metrics, which are (1) Mean ESS: the average of ESS of all variables. That is, ESS is computed for each variable/dimension. and Mean ESS is the average of them. (2) Min ESS: the lowest value of ESS among the ESSs of all variables; (3) Mean FSSS/Sec.: Mean ESS normalized by the CPU time in seconds; (4) Min ESS/Sec.: Minimum ESS normalized by the CPU time in seconds. In the following tables, denote "Cheby." our proposed method, and "Const." as HMC with the the constant integration time (Vishnoi] 2021 CCen & Vempala 2019). Each of the configurations repeated 10 times, and we report the average and the standard deviation of the results. We also report the acceptance rate of the Metropolis filter (Acc. Prob) on the tables. Our implementation of the experiments is done by modifying a publicly available code of HMCs by Brofos & Lederman (2021). Code for our experiments can be found in the supplementary.
4.1 IDEAL HMC FLOW FOR SAMPLING FROM A GUSSIAN WITH DIAGONAL COVARIANCE
Before evaluating the empirical performance of Algorithm 2 in the following subsections, here discuss and compare the use of arbitrary permutation of the Chebyshev integration time and that without permutation (as well as that of using a constant integration time). We simulate ideal HMC
for sampling from a Gaussian N (p, E), where [0] and Σ =
0 It is noted that ideal 100
HMC flow for this case has a closed-form solution as shows. The result are reported on Table
From the table, the use of a Chebyshev integration time allows to obtain larger ESS than that from using a constant integration time, and a arbitrary permutation helps get better result. An explanation is that the ESS is quantity that computed along the trajectory of a chain, and therefore a permutation the integration time could make a difference. We remark that the observation here (a arbitrary permutation of time generates larger ESS) does not contradict to Theorem I since Theorem is about the guarantee in W2 distance at the last iteration K.
4.2 SAMPLING FROM A GAUSSIAN
We sample (H, )), where 0 and Σ
1
0.5
Therefore, the strong convexity constant 100
0.5
m is approximately 0.01 and the smoothness constant L is approximately 1. Table shows the results. HMC with Chebyshev integration time consistently outperforms that of using the constant integration time in terms of all the metries: Mean ESS, Min ESS, Mean ESS/See, and Min ESS/Sec.
We also plot two quantities throughout the iterations of HMCs on Figure Specifically, Sub-figure (a) on Figure plots the size of the difference between the targeted covariance and an estimated covariance Σg at each iteration & of HMC, where Σg is the sample covariance of 10, 000 samples collected from number of 10, 000 HMC chains at their ky, iteration. Sub-figure (b) plots a discrete TV distance that is computed as follows. We use a built-in function of Numpy to sample 10, 000 samples from the target distribution, while we also have 10, 000 samples collected from number
7

---

Published as a conference paper at ICLR 2023
Table 2: Sampling from Gaussian distribution. We report 4 metrics regarding ESS (the higher the better), please see the main text for their definitions.
Step Size
Method
Mean ESS
Min
ESS/Sec.
ESS/Sec.
Acc. Prob
0.001
Cheby.
5187.28 +: 261.13
307. 09.0090 = 21 $1.92
20.28 1 1.74
1.20 = 0.11
1.00 1 0.00
0.001
Const.
1912.76 72.10
39.87
13.77
15.87
0.33 0.11
1.00 1: 0.00
0,005
Cheby,
5146.71 = 257.65
304.126 = 19.09
97.84 + 9.23
5.79 0.68
1.00 : 0.00
0.005
Const.
1926.71 = 136.53
32.83 9.57
80.31 = 4.39
1.37 = 0.39
1.00 - 0.00
0.01
Chehy,
5127.90 +:
279.50 38.09
20.99
10.01
1.00 + 0.00
0.01
Const.
1832.87 + 77.47
35.71 = 11.74
147.53 + 12.69
2.85 0.95
1.00 = 0.00
0.05
Chehy,
5133.67 1: 195.07
316.87 1 36.27
871.72 1: 88.73
53.54 6.22
0.99 1 0.00
0.05
Const.
1849.15 92.75
34.98 14.70
615.73 +: 30.16
5.07 11.70
0.99 0.00
0.1
Cheby,
1918.46 = 144.03
281. = 11.79
1992,96
166.21
$1.39
13.04
0.99 = 0.00
0.1
Const.
1852.79 + 132.95
38.17 16.35
1033.54 82.34
21.44 9.31
0.99 + 0.00
distance the true coverience iteration
discrete TV discance iteration
100
Cheyyshev time
Chebyshev time
disance of conatianue
80
Consrant time
dizreere TV distacce
Constant time
60
100
250
0.0g
250
Iteration
Iteraiio
(a)
[
Σg F V.S. iteration &
(b) discrete (*t Pa) v.s. iteration
Sampling from a Gaussian distribution. Both lines correspond to HMCs with the same step size h = 0.05 used in the leapfrog steps (but with different schemes of the integration time). Please see the main text for the precise definitions of the quantities and the details of how compute them.
Figure 2:
of 10, 000 HMC chains at each iteration k. Using these two sets of samples, construct two histograms with 30 number of bins for each dimension, we denote them as स and Åv. The discrete TV (s)) iteration k is 0.5 times the sum of the absolute value of the difference between the number of counts of all the pairs of the bins divided by 10, 000, which serves as surrogate of the Wasserstein-2 distance between the true target π and Pk from HMC, since computing or estimating the true Wasserstein distance is challenging.
4.3 SAMPLING FROM A MIXTURE OF TWO GAUSSIANS
For a vector a E Rd and a positive definite matrix Σ E Rddd we consider sampling from a mixture of two Gaussians (aa, Σ) and (i-a, Σ) with equal weights. Denote b := Σ-'a a and = E-! The potential is f(x) = 1 - all - log(1 + exp(- -2x b)), and its gradient is ff(x) = Ar b+ 2b(1 + exp(-2: 22 bb) For each dimension E [d]. we set ald set the covariance > = diag, cica
A and la < 1, see e.g.. Riou-Durand & Vogrine (2022). We set d 10 in the experiment, and simply use the smallest and the largest eigenvalue of to approximate the strong convexity constant m and the smoothness constant I of the potential, which are m = 1 and Í. = 10 in this case. Table 3 shows that the proposed method generates a larger effective sample size than the baseline.
( ()) The potential is strongly convex if Ty
4.4 BAYESIAN LOGISTIC REGRESSION
We also consider Bayesian logistic regression to evaluate the methods. Given an observation (2,, W.). where 2; E d and y/ E (0, ), the likelihood function is modeled as p(ya [2, w) = T+exp(-9117 w))
Moreover, the prior the model parameter
assumed follow Gaussian distribution, p(w) (0, a TTa), where a > O is a parameter. The goal is to sample w E IRd from the posterior, (
p(w][z, 31 =
w). where n is the number of data points in dataset. The potential function f (w) can be written as
f(w) Ei= fá(w),, where fi(w)) log (1 + exp( y:w 21)) +
wl
(14)
2m
8

---

Published as a conference paper at ICLR 2023
Table 3: Sampling from a mixture two Gaussians
Step Size
Method
ESS/Sec.
0.001
Cheby.
2439.86 = 71.83
815.20 1 83.82
22.68 = 0.93
1.57 1 0.81
0.89 = 0.00
0.001
Const.
845.44 31.42
261.14
34.34
12.90
0.52
3.98 0.53
0.91 0.00
0.005
Cheby.
2399.50 100. 12
784.06 +: 82.07
105.97 + 8.78
34.58 4.12
0.89 0.00
0.005
Const.
876.61 = 25.62
277.72 + 30.62
63.80 4.67
20.22 2.62
0.91 0.00
0.01
Cheby.
2341.35 +: 89.99
794.27
1: 48.75
194.81
23.51
66.30 9.89
0.88
0.01
Const.
860.61 = 20.39
235.33
+ 33.73
110.62
14.09
30.40 6.34
0.91 - 0.00
0.05
Caeby,
22141.9 = 87.27
748.66
- 46.18
761.59 = 08.88
256.51
13.76
0.89 0.00
0.05
Const.
653.40 + 41.05
265.70
+ 37.41
376.54
67.83
116.45
22.25
0.91 + 0.00
0.1
Cheby.
206-1.42 = 67.44
657.45
+ 60:44
T162.67 = 84.
370.07
41.11
0.90 0.00
Coust.
182.88
43.91
0.00
We set a = in the experiments. We consider three datasets: Heart, Breast Cancer, and Diabetes binary classification datasets, which are all publicly available online. To approximate the strong convexity constant m and the smoothness constant L of the potential f(w), we compute the smallest eigenvalue and the largest eigenvalue of the Hessian V we use them as estimates of m and L respectively. We apply Newton's 's method to approximately find the maximizer of the posterior. The experimental results are reported on Table Qin Appendix [E.1 due to the space limit, which show that our method consistently outperforms the baseline.
the maximizer of the posterior, and
4.5 SAMPLING FROM A hard DISTRIBUTION
We also consider sampling from a step-size-dependent distribution =(x) x exp(- (x)),. where the potential fa is k-smooth and 1-strongly convex. The distribution is considered in Lee et al. (2021a) for showing a lower bound regarding certain Metropolized sampling methods using constant integration time and constant step size h of the leapfrog integrator. More concretely, the potential is
(1) fn(z): = =t= E (zx), where (xx)
2
2
i=1
(15)
sh
25isd.
In the experiment, we set & = 50 and 10. The results are reported on Table in Appendix E.2 The scheme of the Chebyshev integration time is still better than the constant integration time for this task.
DISCUSSION AND OUTLOOK
The Chebyshev integration time shows promising empirical results for sampling from a various of strongly log-concave distributions. On the other hand, the theoretical guarantee of acceleration that we provide this work is only for strongly convex quadratic potentials. Therefore, a direction left open by our work is establishing some provable accelcration guarantees for general strongly log-concave distributions. However, unlike quadratic potentials, the output (position, velocity) of a HMC flow does not have closed-form solution in general, which makes the analysis much more challenging. A starting point might be improving the analysis Chen & Vempala (2019), where contraction bound of two HMC chains under small integration time 7) ( is shown. Since
the scheme of the Chebyshev integration time requires large integration time 17 θ al some iterations of HMC, a natural question is whether a variant of the result of Chen Vempala (2019) can be extended
We state as an open question: can ideal HMC with a scheme of time-varying integration time achieve an accelerated rate O(v log ! for general smooth strongly log-concave distributions?
large integration time 7) = θ
The topic of accelerating HMC with provable guarantees is underexplored, and we hope our work can facilitate the progress in this field. After the preprint of this work was available arXiv, Tiang (2022) proposes randomized integration time with partial velocity refreshment and provably shows that ideal HMC with the proposed machinery has the accelerated rate for sampling from a Gaussian distribution. Exploring any connections between the scheme of Jiag (2022) and ours can be an interesting direction.

---

Published as a conference paper at ICLR 2023
ACKNOWLEDGMENTS
We thank the reviewers for constructive feedback, which helps improve the presentation of this paper.
REFERENCES
Naman Agarwal, Surbhi Goel, and Cyril Zhang. Acceleration via fractal learning rate schedules. ICML, 2021.
Michael Betancourt. A conceptual introduction Hamiltonian Monte Carlo. 1rXIv:1701.0245:. 2017.
Nawaf Bou-Rabee and Andreas Eberle. Mixing time guarantees for unadjusted Hamiltonian Monte Carlo. arXiv:2105.0088: 2021.
Nawaf Bou-Rabee and Jesus Maria Sanz-Serna. Randomized Hamiltonian Monte Carlo. Annals of Applied Probability, 2017.
James A. Brofos and Roy R. Lederman. Evaluating the implicit midpoint integrator for Riemannian manifold Hamiltonian Monte Carlo. ICML, 2021.
Tianqi Chen, Emily B. Fox, and Carlos Guestrin. Stochastic gradient Hamiltonian Monte Carlo. ICMI., 2014.
Yuansi Chen, Raaz Dwivedi, Martin Wainwright, and Bin Yu. Fast mixing of Metropolized Hamiltonian Monte Carlo: Benefits of multi-step gradients. JMLR, 2020.
Zongchen Chen and Santosh S Vempala. Optimal convergence rate of Hamiltonian Monte Carlo for strongly logconcave distributions. International Conference on Randomization and Computation (RANDOM), 2019.
Sinho Chewi, Thibaut Le Gouic, Chen Lu, Tyler Maunu, Philippe Rigollet, and Austin J. Stromme. Exponential ergodicity of mirror-langevin diffusions. NeurIPS, 2020.
Armak S. Dalalyan. Theoretical guarantees for approximate sampling from smooth and log-concave density. Journal of the Royal Statistical Society: Series B, 2017.
Arnak S. Dalalyan and Lionel Riou-Durand. On sampling from a log-concave density using kinetic Langevin diffusions. Bernoulli, 2020.
Alexandre d' Aspremont, Damien Scieur, and Adrien Taylor. Acceleration methods. Foundations and Trends in Optimization, 2021.
Simon Duane, A. D. Kennedy, Brian Pendleton, and Duncan Roweth. Hybrid monte carlo. Physics Letters B, 1987.
John C. Duchi, Peter L. Bartlett, and Martin Wainwright. Randomized smoothing for stochastic optimization. SIAM Journal on Optimization, 2012.
Alain Durmus, Eric Moulines, and Eero Saksman. On the convergence of Hamiltonian Monte Carlo. arXiv: 2017.
Alain Durmus, Szymon Majewski, and Blažej Miasojedow. Analysis of Langevin Monte Carlo via convex optimization. JMLR, 2019.
Abraham D. Flaxman, Adam Tauman Kalai, and H. Brendan McMahan. Online convex optimization in the bandit setting: gradient descent without a gradient. SODA, 2005.
Mark Girolami and Ben Calderhead. Riemann manifold Langevin and Hamiltonian Monte Carlo methods. Journal of the Royal Statistical Society, 2011.
Marcel Hirt, Michalis K. Titsias, and Petros Dellaportas. Entropy- adaptive Hamiltoniam Monte Carlo. NeurIPS, 2021.
10

---

Published as a conference paper at ICLR 2023
Matthew D. Hoffman and Andrew Gelman. The No-U-Tum sampler: Adaptively setting path lengths in Hamiltonian Monte Carlo. JMLR, 2014.
Matthew D. Hoffman, Alexey Radul, and Pavel Sountsov. An adaptive-MCCMC scheme for setting trajectory lengths in Hamiltonian Monte Carlo. AISTATS, 2021.
Qijia Jiang. On the dissipation of ideal hamiltonian monte carlo sampler. arXiv:2209.07438, 2022.
Chi Jin, Rong Ge, Pranceth Netrapalli, Sham M. Kakade, and Michael I. Jordan. How to escape saddle points efficiently. ICML., 2017.
Daphne Koller and Nir Friedman. Probabilistic graphical models: Principles and techniques. MIT Press, 2009.
Ravin Kumar, Colin Carroll, Ari Hartikainen, and Osvaldo Martin. Arviz unified library for exploratory analysis of bayesian models python. The Journal of Open Source Software, 2019.
Yin Tat Lee, Ruoqi Shen, and Kevin Tian. Lower bounds on Metropolized sampling methods for well-conditioned distributions. NeurIPS, 2021a.
Yin Tat Lee, Ruoqi Shen, and Kevin Tian. Structured logconcave sampling with restricted gaussian oracle. COLT, 2021b.
Daniel Levy, Matthew D. Hoffman, and Jascha Sohl-Dickstein. Generalizing Hamiltonian Monte Carlo with neural networks. ICLR, 2018.
Ruilin Li, Hongyuan Zha, and Molei Tao. Sqrt(d) dimension dependence of Langevin Monte Carlo. ICLR, 2022.
Yi-An Ma, Niladri Chatterji, Xiang Cheng, Nicolas Flammarion, Peter L. Bartlett, and Michael I. Jordan. Is there an analog of Nesterov acceleration for MCMC? Bernoulli, 2021.
Oren Mangoubi and Aaron Smith. Mixing of Hamiltonian Monte Carlo strongly logconcave distributions 2: Numerical integrators. AISTATS, 2019.
Oren Mangoubi and Aaron Smith. Mixing of Hamiltonian Monte Carlo on strongly logconcave distributions 1: Continuous dynamics. Annals of Applied Probability, 2021.
Oren Mangoubi and Nisheeth Vishnoi. Dimensionally tight bounds for second-order Hamiltonian Monte Carlo. NeurlPS, 2018.
Radford M. Neal. MCMC using Hamiltonian dynamics. arXiv: 2012.
Yurii Nesterov. Introductory lectures on convex optimization: basic course. Springer, 2013. 20133
Fabian Pedregosa. Acceleration without momentum, 2021. URL httt:://fa.banpp.cee
Dlog/2021/no-morree um/
B.T. Polyak. Some methods of speeding up the convergence of iteration methods. USSR Computa-
tional Mathematics and Mathematical Physics, 1964.
Lionel Riou-Durand and Jure Vogrine. Metropolis Adjusted Langevin trajectories: a robust alternative to Hamiltonian Monte Carlo. 1333 2022.
Ruslan Salakhutdinov and Andriy Mnih. Bayesian probabilistic matrix factorization using Markov chain Monte Carlo. ICML, 2008.
Greg Ver Steeg and Aram Galstyan. Hamiltonian dynamics with non-newtonian momentum for rapid sampling. NewrIPS, 2021.
Nilesh Tripuraneni, Mark Rowland, Zoubin Ghahramani, and Richard Turner. Magnetic Hamiltonian Monte Carlo. ICML, 2017.
Santosh S. Vempala and Andre Wibisono. Rapid convergence of the Unadjusted Langevin Algorithm: Isoperimetry suffices. NeurlPS, 2019.
11

---

Published as a conference paper at ICLR 2023
Nisheeth K. Vishnoi. An introduction Hamiltonian Monte Carlo method for sampling. arXiv:2108. 12107, 2021.
Jun-Kun Wang, momentum. ICLR, 2020.
Lin, and Jacob Abemethy. Escaping saddle points faster with stochastic
Jun-Kun Wang, Chi-Heng Lin, and Jacob Abernethy. A modular analysis of provable acceleration via Polyak's 's momentum: Training a wide ReLU network and deep linear network. ICML, 2021.
Jun-Kun Wang, Chi-Heng Lin, Andre Wibisono, and Bin Hu. Provable Acceleration of Heavy Ball beyond Quadratics for Class of Polyak-Lojasiewiez Functions when the Non-Convexity is Averaged-Out. ICML, 2022.
Max Welling and Yee Whye Teh. Bayesian learning via Stochastic Gradient Langevin dynamics. ICML, 2011.
Florian Wenzel, Kevin Roth, Bastiaan S. Veeling, Jakub Swiatkowski, Linh Tran, Stephan Mandt, Jasper Snoek, Tim Salimans, Rodolphe Jenatton, and Sebastian Nowozin. How good is the Bayes posterior in deep neural networks really. ICML, 2020.
Difan Zou and Quanquan Gu. On the convergence of Hamiltonian Monte Carlo with stochastic gradients. ICML, 2021.
A A CONNECTION BETWEEN OPTIMIZATION AND SAMPLING
To provide an intuition of why the technique of Chebyshev polynomials can help accelerate HMC for the case of the strongly convex quadratic potentials, we would like describe the work of gradient descent with the Chebyshev step sizes Agarwal et al. (2021) in more detail. because we are going to draw a connection between optimization and sampling to showcase the intuition. Agarwal al. (2021) provably show that gradient descent with scheme of step sizes based on the Chebyshev Polynomials has an accelerated rate for minimizing strongly convex quadratic functions compared to GD with constant step size, and their experiments show some promising results for minimizing smooth strongly convex functions beyond quadratics via the proposed scheme of step sizes. More precisely, define f(w) 1/w' Aw, where A Rdxd is positive definite matrix which has eigenvalues L:= \1 2 A2 Z ... Ad =: m. Agarwal et al. 2021) consider applying gradient descent
Wk+1 = Wi
Th f((k))
to minimize f()), where 7is is the step size of gradient descent at iteration Let be the unique global minimizer of (c). It is easy to show that the dynamic of the distance evolves as
Wk+1 w. = (Ia A) 7m 7
(Ia A)(1 - w.).
Hence, the size of the distance to at iteration K + 1 is bounded by
| w.ll S max jE(d]
II (1 mAj))llu w-
This shows that the convergence rate of GD is governed by maxje(4] II (1 n)d))). By setting ? as the inverse of the Chebyshev root rk definition), the polynomial
(KK) (K) or any permuted root ra(k) (see @ for the
I kk=1 (1 m)) is actually the K-degree scale-and-shifted polynomial, -I
(A) (see (7) for the definition). It is well-known in the literature of optimization and numerical linear algebra that the K K-degree scale-and-shifted polynomial satisfies
i.e., _ (1 ))
K =0 =0(11- (1-10(11 m
K
XE m..... [x((a) <2 1-2
2-122
/ m
)
max
VT. + vm
12

---

Published as a conference paper at ICLR 2023
which is restated in Lemma and its proof is replicated in Appendix B of our paper for the reader's 's convenience. Applying this result, one gets a simple proof of the accelerated linear rate of GD with the proposed scheme of step sizes for minimizing quadratic functions. A nice blog article by Pedregosa (2021) explains this in detail.
Now ready highlight its connection with HMC. In Lemma 1 of the paper, we restate a known result in HMC literature, where its proof is also replicated in Appendix B for the reader's convenience. The lemma indicates that the convergence rate of HMC is governed by maxje[d] I cos( 22 (x)) .) By way of comparison to that of GD for minimizing quadratic functions, i.e., max I us wonder if we could bound the former by the latter. We show in Lemma that cos 2 which holds for all x 2 0, and consequently,
m)dr)) it appears that they share similarity, which made
K
A
A
P & = II cos
II 1=1
1
= c (2)
(K)
(K)
(kk
(A)
The key lemma above implies that if we set the integration time as T) acceleration of HMC.
(K)
2 22
then get
B PROOFS OF LEMMAS IN SECTION 2
We restate the lemmas for the reader's convenience.
Lemma |Vishnoi 2021) Let Top Mo IR4. Consider the following coupling: (at, 0z) HMC, (MC, and (yt, u.) HMCo(yo. {) for some & Rd. Then for all 2 0 and for all j E [d]. it holds that
205 volil cos (V21jt) (xoli) yo(d1)-
ue) == (Vo, E), we have du: du, -vf(z)) + ff(v)) = 2A(v: 2²). Therefore, we have (= (aslli--(lll) veld)) (xy (i),, for all j E |d]. Because of the initial condition deglal duultl = E[i.. the differential equation de
Proof. Given (I,, = HMC, E) and (bt)
di
dea
dt implies that 2, yo |-] = cos (22X)t) x (zoly) yo (il).
is noted that the result also follows directly from the explicit solution
Lemma Vishnoi 2021) Let x exp(- (0, be the targel distribution, where (x) defined on (3). Let pR be the distribution of k generated by Algorithm at the final iteration K. Then for any pa and any K 2 1, we have
W2(pK ④) S maxjeld cos cos (V2A,m( W2(po.
Proof. Starting from so Po. draw an initial point yo such that (æo, Ho) has the optimal
We-coupling between Po and n. Consider the following coupling at each iteration k: (Ik- U%) =
HMC (x)- and (yk, ux) HMC (i) Ex) where & (0, I) is independent "h
Gaussian. We collect 1k and (ye ] from Algorithm 1 We know each yı n,, since is
13

---

Published as a conference paper at ICLR 2023
stationary distribution of the HMC Markov chain. Then by Lemma
we have
(pK, T) - E yx
yx11
E|E (Ix [i]
yK
[51)2]
c (v (22gyTh (x)) x (zoli) - 8ol61) ?
E JE(4]
(16)
(v2XmTk )) 2)
maxjel4]
EE (xols)) -601111
(v2X)1k*)) 2 W2 (po,
max/e[d]
IT c cos
Taking the square root on both sides leads to the result.
Lemma (e.g., Section 2.3 in d' d'Aspremont et al. (2021) ) For any positive integer K, we have
K = -0((1-0(D)))
maxxe[m,r1 (x $2 (1-2
vm VL+Vm
(17)
Proof. Observe that the numerator of & (A)
(b()
xx(0) satisfies | (h())) < 1, since
((A) E [] for X E [m, L] and that the Chebyshev polynomial satisfies [xx(-) S 1 when
its argument is [-1,1] by the definition. remains to bound the denominator, which Φx (h(()) cosh (K arccosh
L+m
Since
m.
arccosh
L+m L-m
log
L+t. L-m
Ltt L-m
-1
-log(0), where 0: = VT++.mm VI.-Vm
we have
$ xx(h(0)) = cosh (K'
arccosh
L+
exp(K log( +expp K log(O))
8 +9 k
2
gk
2
Combing the above inequalities, we obtain the desired result:
4 (h()))
2 OK
m K
max XE(m,L] AE(m.L]
(x<(A)|=
(
=
1-2 VL L + vm
=0 (1-1-0 0 (1 111
L
C PROOF OF LEMMA 4
Lemma Denote [PC (A)) = II cos
A
Suppose E [m, L). Then, we have for
2
I8T (ock)
any positive integer K,
PC (X)l (A)]-
(18)
Proof. We use the fact that the K Kddegree scaled-and-shifter Chebyshev Polynomial can be written
as.
xA) = IT (1
(19)
TKT rork)
14

---

Published as a conference paper at ICLR 2023
for any permutation a(),, since (R)) are its roots and xx (0) to
So inequality (18)
equivalent
(1
11
COS
11 Rm1
J.R (())
(20)
To show (20). let us analyze the mapping (x) == cos for æ 2 0, 0 1, with (1) =
4 by
continuity, and show that max [x(x))] < 1, 20 would be immediate. We have v'(x)) when
-- sin( Vz) + cos( vx) (I-z)Z Hence, ab
Lan(3vo)
4 x[l-2)
(21)
Denote an extreme point of (xx) as I, which satisfies
21
Then, using 21 we have [{)(
000
V
(12++ where we used cos( +
16 or
=(1-a)
(1
The denominator /162 7² 2 has the smallest value 0, which
v (1-
means that the largest value of [(b(1) happens at r 0, which is 1. The proof is now completed.
D A COMPARISON OF THE TOTAL INTEGRATION TIME (JIANG 2022)
Since the Chebyshev integration time are set to be some large values at some steps of HMC, it natural to ask the number of steps to get an 2- 2-Wasserstein distance a fair metric. In this (K)
section, we consider the total integration time k=1 7h to get an distance as another metric for the convergence. It is noted that the comparison between HMC with our integration time and HMC with the best constant integration time has been conducted by [iang (2022), and our previous version did not have such comparison. Below, we reproduce the comparision of Jiang (2022).
Recall the number of iterations to get an 2- W-ssseerstein distance to the target distribution is K = O (Vn log (t)) of HMC with the Chebyshev integration time (Theorem) in the paper). The average of the integration time is
K
K
K
(K)
π
π
K
K
2/2
(K)
K
2
(K)
1=
o(k)
k=1
where we recall that a permutation (-)) does not affect the average.
Then, if K is even, rewrite the averaged integration time as
K
K/2
K
K
K2
k=1
K
Otherwise, K is odd, and we can rewrite the averaged integration time as
K Σ
(K (K-1)/2
m
+
K
K22
K
Vr
(K)
Vre
V v((k)
We will show
1
(K)
(K)
(R)
(K) VKK-LK/2++1
rk
TK+1-k
for any l: = {1,
K [
KJ} soon. Given this, we can further upper-bound the averaged integration
time
[K]
1k
2
V
TLK/2]
(K)
1K-//111-1
(K)
(K)
15

---

Published as a conference paper at ICLR 2023
when K is even; when K is odd, upper-bound the averaged integration time as
K
1
1
K -
T
K
K
2/2
Vr
(K)
1=1
(K++1)/2
K/2]
V R [K/2]+1
Using the definition of the Chebyshev root, have
(K) K/2]
I+m
L -
m
(L 2)
L + m
cos
2
K
(1 J b)=
(1
where the approximation is because 0. Similarly, we can approximate
3h when K is large, and hence cos
(K) K [K/2]++
L + m
m 2
(K (K- K/2] 1
L +
m
cos
2
K
2
as
(K [K/2] +1 1
1 when K is large. Also, we can approximate (K)
(K+11//2
Lto
2
when K
odd and large for the same reason.
Combining the above, the total integration time of HMC with the Chebyshev scheme be approxi- mated as
number of iterations
average integration time
Vo ())
K
(K)
VKloo
())
1 VL + m
m
K
k=
When x == is large, the total integration time becomes
VK log (t)
1 0 (
(J log )t)
x
(22)
LL 1 +m
Now let us switch analyzing HMC with the best constant integration time 17 θ (+) (see e.g..
(5). Vishnoi (2021)), which has the non-accelerated rate. Specifically, it needs iterations to converge to the target distribution. Hence, the total integration time of HMC with the best constant integration time is
log (1)
number of iterations average integration time log (1) x
1 VL
VL m
=
log
(23)
see that the total integration time of HMC with the proposed scheme of Chebyshev integration time reduces by a factor Vk, compared with HMC with the best constant integration time.
By way of comparison (2 vs.
The remaining thing to show is the inequality
(24)
(R)
(K)
(K)
[K] V 1K+1-1K/21
rk
V
VKKII-K
V
ViK/2
for any l: {1, ...,
[K []]
16

---

Published as a conference paper at ICLR 2023
We have
1
(K)
V rK+1-k
x2x
(K=
(K-k )rr
+ m - (L m)cos
1+ + m - -LL m)cos
k
K
(k-
(k 2)
T. + (T. m)oos
I. + m + (T. m)ooss
(25)
Now let us define H(k) =
+
and
T+tm 1l-majoa (--b)r
l+m+ m)oos (k-
treat & as continuous variable.
The derivative of (k) is
(k 2
(k:)) 2K
(L m)sin
K
1
(k-
3/2
3)=
3/2
+ m (L m)cos
L+ m + (I m)cos
K
0.
(26)
That is, (k) is an increasing function of l when 1 S 1: S LK) which implies that the inequality 24h. Now have completed the analysis.
17

---

Published as a conference paper at ICLR 2023
E EXPERIMENTS
E.1
BAYESIAN LOGISTIC REGRESSION
Table Bayesian logistic regression
HEART dataset
2.59, 92.44)
Step Size
Method
Mean FSS
Min ESS
Mean FSS/Sec.
Min. ESS/Sec.
Acc. Proh
0.001
Cheby,
1693.71 +
520.43 :
18.54
5.69
0.00
0.001
Coust.
312.18 = 12,60 65
80.97 + 15.97
6.07 0.42
1.69 0.28
1.00 0.00
0.005
Cheby.
166-1.87 = 13.
181.76 + 19.00
82,90 = 16.51
24.08 4 5.72
0.99 4 0.00
Const.
329.48
75.78
31.87
7.40
0.99
0.01
Cheby.
16:18.25
17.50
508.69 = 49.81
157.09 = 26.70
18.15
0.99 = 0.00
0.01
Const.
307.52
8.77
82.85 +: 13.88
53,89
6.37
3.22
0.99 + 0.00
Cheby.
1421.21 = 54.03
439.85
+ 56.25
458.56
51.33
140.51
16.58
0.98 0.00
0.05
Const.
242.44 + 14.61
56.42 + 17.68
103.36 + 12.64
23.90 7.40
0.98 + 0.00
BREAST CANCER dataset (fi 1.81, 7.
69.28)
Step Size
Metbod
Mean ESS
Min ESS
Mean ESS/Sec.
Min.
ESS/Sec.
Acc. Prob
0.001
Cheby,
1037.98 + 34.46
676.72 + 41.14
9.40 0.31
5.21 0.331
1.00 1 0.00
0.001
Coust.
171.73 + 13.91
78.24 + 23.28
2.59 + 0.29
2.59 : 0.29
1.00 + 0.00
0.005
Cheby,
1010.49 = 24.15
571.03
+ 36.64
43.09
1.14
24.35
1.70
0.99 = 0.00
0.005
Const.
173.17 1: 11.40
79.76
1. 13.49
11.88
1.39
11.88
1.39
0.99 1. 0.00
0.01
Cheby.
1038.10 + 31.48
565.54
+ 30.51
82.82 = 3.51
15.14
4.41
0.99 + 0.00
Const.
162.64
18.92 1. 2.59
2.59
0.99
0.05
Cheby.
886.21 = 38.92
199.5/
+ 43.99
240.08 = 12.55
135.28
12.01
0.98 = 0.00
0.05
Const.
99.48 1: 10.10
44.70 +: 13.23
33.25 6.50
33.25
6.50
0.98 1: 0.00
BIABETES dataset (A
4.96, m 270.20)
Step Size
Method
Mean ESS
Min ESS
Mean ESS/Sec.
Min. ESS/Sec:
Ace. Prob
0.001
Cheby.
726.08 33.92
424.59 + 58.77
11.61 0.85
6,83 = 1.16
0.99 = 0.00
0.001
Const.
100.30 9.32
41.84
19.33
3.6 + 0.31
1.50 + 0.68
0.99 0.00
0.005
Cheby:
731.46 33.04
395.82 22 47.98
64.92 5.26
29.61 = 3.76
0.99
0.00
0.005
Const.
100.16 + 11.83
44.62 + 20.81
14.71 2.52
6.67 + 3.37
0.99 = 0.00
0.01
Cheby.
687.74
399.41
45.01
93.10 6.78
53.90 =: 5.38
0.98
0.00
0.01
Const.
83.04 9.36
36.39
12.43
20.87 3.331
9.09
3.25
0.98
0.00
0.05
Cheby.
546.80 + 37.40
330.09 + 31.31
206.07
17.76
125.07 : 18.87
0.96
0.00
4.34
E.2 SAMPLING FROM hard DISTRIBUTION
Table 5: Sampling from a distribution ((r) x exp(- (xx)) whose potential fa(-) is defined on
15
Step Size
Method
Mean ESS
Min ESS
Mean ESS/See.
Min. ESS/Scc.
Acc. Prob
sampling from T(x)
exp( fo (x]) Choby.
0.001
6222 21 + 88.90
453.03 + 30.35
114.71 7.59
8.36 = 0.83
1.00
# 0.00
0.001
Const.
2008.
46.56
63.53 + 15.00
82.31 5.29
2.50
1.00 +
sampling from (())) x expl - Jo oos (x))
0.005
Cheby,
6271.43 117.71
429.42
31.52
545. + 26.10
37.28
= 2.29
0.99
= 0.00
0.005
Const.
21.87 2125.36
67.42 16.51
361. +: 5.65
11.44 = 2.76
0.99 # 0.00
sampling m(x))
exp(=foo112]]
0.01
Cheby,
6523.21
95.65
459.48
38.83
1070.77 77 1 68.78
75.61
9.79
0.99 1. 0.00
0.01
Const.
2125 + 31.83
69.66 2 20.75
528..35 + 80.17
17.19 = 6.34
0.99 = 0.00
sampling from ((n)
exp(=fa.ns ox(r)))))
0.05
Cheby.
6457.21 4: 110.05
375.97 1 30.64
3319.51 + 134.92
193,066 = 14.49
0.97 =: 0.00
0.05
Const.
2796.41 56.89
62.33
13.26
1893.99 + 57.23
42.22
9.03
0.97 = 0.00
18


Reviews:
[{"review_id": "fGNJtqOyYZp", "reviewer": "Reviewer_XHWf", "summary": "", "questions": "", "limitations": "", "rating": 8, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "This paper considers the technique of sampling from a Gibbs distribution via the Ideal Hamiltonian Monte-Carlo (HMC) method. Specifically, the ideal HMC method is derived by iterative running the HMC from the current spatial position, but with a resampled velocity for a certain period of integration time. The integration time plays a key role in proving the mixing time of the algorithm: It is known that even for the Gaussian case, with a constant integration time, the mixing time is proportional to $\\Omega(\\kappa)$. To improve this, this paper proposes the Chebyshev integration time, which are derived from the roots of the Chebyshev polynomials. It is proved that using the proposed integrating time, the resulting ideal HMC method enjoys a mixing time of order $O(\\sqrt{\\kappa})$, for the Gaussian case. The authors also provide empirical evidence that the improved rate generalizes to non-Gaussian tasks.", "summary_of_the_review": "The Chebyshev integration time is novel, but the current analysis is limited to the Gaussian case. Intuition behind the success of this technique is missing.", "strength_and_weaknesses": "Strength:\nThis paper proposes the Chebyshev integration time which leads to a provably faster instance of the Ideal HMC method, in the Gaussian case. Empirical evidence shows that the improved rate generalizes to non-Gaussian tasks.\n\nWeakness:\nThe current theoretical result only applies to the Gaussian case. \n", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "3: Some of the paper’s claims have minor issues. A few statements are not well-supported, or require small changes to be made correct.", "technical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work.", "empirical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work."}, "extra_text": {}}, {"review_id": "YvDyheYzRpI", "reviewer": "Reviewer_dujW", "summary": "", "questions": "", "limitations": "", "rating": 6, "confidence": 2, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "This paper proposes a new time-varying integration duration for HMC using the roots of Chebyshev polynomials. \nIn a very restrictive setting, they paper proves an upper bound on the number of iterations required to reach Wasserstrain-2 distance less than an a specified error threshold. This setting is where exact HMC mechanism is applicable (i.e. the equations of motion can be solved in closed form). This closed-form solution only exists for Gaussian distributions (as far as I know).\nNonetheless, experiments shows that their scheme has benefits in case the target distribution is not Gaussian but is associated with a smooth strongly convex potential energy function.\n ", "summary_of_the_review": "This paper presents a new method to adapt the HMC's integration time. It is shown that on the experimental models, the proposed method can outperform the baseline. \nThey also find new bounds for a very restrictive case where HMC's state evolution can be solved in closed-form. As such, I am not sure if it can be called significant. Nonetheless, I think this finding is interesting and valuable.  It may also open doors to future theoretical works. \n \n\n", "strength_and_weaknesses": "The paper is well-written and the proofs seem correct. Though, I should highlight I did not check all the proof details closely. \nThe newly found bound is interesting though the studied theoretical setting, namely exact HMC (i.e. where the equations of motion are solved in closed-form) is only applicable to Gaussians and therefore, it is quite restrictive because clearly, for sampling from Gaussians there exist better algorithms and HMC is not used. \nNonetheless, this theoretical finding may open the way for future theoretical works that consider more general families of distributions. \nAlso, the presented empirical results show that the proposed algorithm can outperform the baseline on models other than Gaussian. \n  ", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "4: All of the claims and statements are well-supported and correct.", "technical_novelty_and_significance": "2: The contributions are only marginally significant or novel.", "empirical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work."}, "extra_text": {}}, {"review_id": "HwJJ4YcbuD", "reviewer": "Reviewer_GRbp", "summary": "", "questions": "", "limitations": "", "rating": 5, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "The paper considers the possibility of accelerating Hamiltonian Monte Carlo (HMC) methods for sampling from distributions $\\pi$. For a $L$-smooth and $m$-convex $f$, the complexity of sampling from $\\pi \\;\\alpha\\; e^{-f}$ (via an \"ideal method\" that does not discretize time) is of the order $\\kappa \\log(1/\\delta)$, where $\\delta$ is the desired precision and $\\kappa={L/m}$ is the condition number. The present paper shows that one can lower the complexity of the ideal method by a factor of $\\sqrt{\\kappa}$ when $f$ is quadratic (and so $\\pi$ is Gaussian). This requires knowledge of $L$ and $m$, which one uses to choose special step sizes defined via Chebyshev polynomials. Experiments suggest that the same choice of weights works beyond the case of quadratic $f$; in fact the method seems to \"beat the competition\" even for nonconvex likelihoods. ", "summary_of_the_review": "The paper has a nice idea, but, at it is present form, its interest is too limited from both theoretical and practical perspectives. \n\nUPDATE on Nov 30th\n\nI have increased my \"novelty\" and overall scores. I thank the authors for their responses. To comment on a few points. \n\na) I agree ESS is a standard surrogate for convergence. It would be nice, however, to see a problem with eg. a multimodal likelihood and show (empirically) that accelerated HMC converges faster than the standard method. \n\nb) I understand that the above point is (as the authors say) disjoint from their goals, as it'd go beyong the logconcave setting. However, given the lack of theoretical results for general logconcave f, I think the paper would greatly benefit from a more extensive empirical evaluation of the accelerated method in more practical settings. \n\nc) I thank the authors for the clarifiction regarding sensibility to the smoothness ans strong convexity parameters. ", "strength_and_weaknesses": "*Strengths*\n\n- This seems to be the first acceleration result for HMC.  \n- The use of Chebyshev polynomials is somewhat interesting.\n\n*Weaknesses*\n\n- Unfortunately, the theory does not cover any cases of practical interest.\n- The experimental success metrics are all defined in terms of effective sample size (ESS). This is a useful proxy of MCMC quality, but it does not prove convergence, nor does it quantify this convergence in harder problems. It would be nice to see other metrics as in: can one approximate the integrals of some \"interesting/hard functions\" from the MC trajectories of different methods?\n- In line with the previous remark, it'd be nice to look at some *hard* distributions coming from real-life Bayesian statistics. What do the credible regions look like? Does one reobtain known results at lower cost?\n\n", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "4: All of the claims and statements are well-supported and correct.", "technical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work.", "empirical_novelty_and_significance": "2: The contributions are only marginally significant or novel."}, "extra_text": {}}, {"review_id": "BQp2jqEzoJ", "reviewer": "Reviewer_Kxzm", "summary": "", "questions": "", "limitations": "", "rating": 5, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "This paper proposed a time-varying integration method for accelerating the mixing of Hamiltonian Monte Carlo. In particular, when the potential $f$ is quadratic function, i.e., the target distribution is Gaussian, the ideal HMC with the proposed time-varying integration enjoy a $O(\\sqrt{\\kappa}\\log(1/\\epsilon))$ iteration complexity in terms of Wasserstein-2 distance. The authors further develop a practical HMC algorithm using the proposed integration time. Experimental results show that the developed algorithm can also lead to better sampling performance for the general smooth and strongly convex potentials.", "summary_of_the_review": "My major concern is (1) the limitation of theory (i.e., only for quadratic potential) and (2) lack of thorough discussion and exploration of baseline method (e.g., constant integration time) in terms of both theory and experiment.\n", "strength_and_weaknesses": "Strength:\n\n* this paper proposes a novel time-varying integration method for HMC and proves a faster convergence compared to the constant-time integration method.\n* this paper incorporated the proposed integration into a practical HMC algorithm and demonstrate its better performance in experiments.\n\nweakness:\n* The theory is limited to quadratic cases, while prior works cover all smooth and strong convex cases.\n* The explanation of the acceleration effect of Chebyshev integration time is not clear.\n* The stepsize is not tuned for HMC with constant-time integration (as mentioned in Section 4, the stepsize for HMC with constant integration time is set by (5)).\n\nQuestions:\n* What’s the role of $\\bar \\Phi_K(\\lambda)$ in the algorithm or theory? Why can equation (11) illustrate the acceleration of the proposed Chebyshev integration time?\n* The authors mention that the permutation $\\sigma$ is not needed in the analysis, does it mean that Theorem 1 holds for any permutation?\n* (5) is only one choice of constant integration time, can we find a better choice to gain faster convergence? \n* It would be better to compare the total integration time of the proposed method and baselines.\n* Following the previous comment, can we fix the total integration time and directly optimize the integration time in each iteration according to Lemma 2?\n", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "4: All of the claims and statements are well-supported and correct.", "technical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work.", "empirical_novelty_and_significance": "2: The contributions are only marginally significant or novel."}, "extra_text": {}}]

Author rebuttal:

Decision: Accept: poster