Hessian-Free High-Resolution Nesterov Acceleration for Sampling
Authors: Ruilin Li
Hongyuan Zha
Molei Tao
Venue: ICLR 2022 Submitted (2022)
Keywords: Markov Chain Monte Carlo
Nesterov Accelerated Gradient
accelerated sampling

Abstract: It is known (Shi et al., 2021) that Nesterov's Accelerated Gradient (NAG) for optimization starts to differ from its continuous time limit (noiseless kinetic Langevin) when its stepsize becomes finite. This work explores the sampling counterpart of this phenonemon and proposes an accelerated-gradient-based MCMC method, based on the optimizer of NAG for strongly convex functions (NAG-SC): we reformulate NAG-SC as a Hessian-Free High-Resolution ODE, change its high-resolution coefficient to a hyperparameter, inject appropriate noise, and discretize the resulting diffusion process. Accelerated sampling enabled by the new hyperparameter is quantified and it is not a false acceleration created by time-rescaling. At continuous-time level, additional acceleration over underdamped Langevin in $W_2$ distance is proved. At discrete algorithm level, a dedicated discretization is proposed to simulate the Hessian-Free High-Resolution SDE in a cost-efficient manner. For log-strong-concave-and-smooth target measures, the proposed algorithm achieves $\tilde{\mathcal{O}}(\sqrt{d}/\epsilon)$ iteration complexity in $W_2$ distance, same as underdamped Langevin dynamics, but with a reduced constant. Empirical experiments are conducted to numerically verify our theoretical results.

Under review as a conference paper at ICLR 2022
HESSIAN-FREE HIGH-RESOLUTION NESTEROV AC- CELERATION FOR SAMPLING
Anonymous authors Paper under double-blind review
ABSTRACT
is known (Shi et al. 2021) that Nesterov's Accelerated Gradient (NAG) for opti- mization differs from its continuous time limit (noiseless kinetic Langevin) when its stepsize becomes finite. This work explores the sampling counterpart of this phenonemon and proposes an accelerated seteerrttddgrrddeet-bbasss MCMC method, based the optimizer of NAG for strongly convex functions (NAG-SC): reformu- late NAG-SC as a Hessian-Free High-Reson Resolution ODE, change its high-resolutin coefficient hyperparameter, inject appropriate noise, and discretize the result- ing diffusion process. Accelerated sampling enabled by the new hyperparame- ter is quantified and not a false acceleration created by time-rescaling. At continuous-time level, additional acceleration over underdamped Langevin in W2 distance is proved. At discrete algorithm level, a dedicated discretization is pro- posed to simulate the Hessian-Free High-Reae Regh-Resolution SDE in a cost-efficient man- ner. For stronttcooccaee aadeeee--ndd-iee and-smooth target measures, the proposed algorithm achieves Ó(Vd/c) iteration complexity in W2 distance, same as underdamped Langevin dynamics, but with reduced constant. Empirical experiments con- ducted to numerically verify our theoretical results.
INTRODUCTION
Optimization is major machinery that drives both the theory and practice of machine learning in recent years. Since the seminal work of Nesterov (1983). acceleration has played a key role in gradient-based optimization methods. A notable example is Nesterov's Accelerated Gradient (NAG), which is an instance of a more general family of "momentum methods". NAG in fact consists of multiple methods, including NAG-C and NAG-SC, respectively for convex and strongly convex functions. Both of them provably converge faster than vanilla gradient descent (GD) in their corresponding setups (Nesterov] 1983] 2013). Newer perspectives of acceleration continue to be revealed, e.g.. [Su et al. (2014); Wibisono et al. (2016); Wilson et al. 2021); ); Hu & Lessard (2017); Attouch et al. (2018); Shi et al, (2021). many based on the interplay between continuous and discrete times. This work aims at turning NAG-SC into sampler based this interplay.
In fact, approaches for sampling statistical distributions, such as gradient- Markov Chain Monte Carlo (MCMC) methods, are also of great importance in machine learning, for example due to their links to statistical inference and abilities to represent uncertainties lacking in optimization- based methods. Although not entirely the same thing, optimization and sampling closely related: besides seeing a large class of sampling dynamics as optimization dynamics with additional noise, viewing sampling as optimization in probability space is another profound perspective that led to fruitful discoveries (e.g.. Jordan er al. (1998): ; Liu & Wang (2016); Dalalyan (2017a); Wibisono (2018) ; Zhang et (2018); Frogner & Poggio (200):: Chizat & Bach (2018): Chen al. (2018a); Ma et al. (202)): Erdogdu & Hosseinzadeh (2021). In fact, an unadjusted Euler-Maruyama dis- cretization of overdamped Langevin dynamics (abbreviated OLD here) is commonly considered as the analog of GD in sampling (although many other discretizations are also possible), and often referred to as Unadjusted Langevin Algorithm (ULA) (Roberts et al. (996) and/or Langevin Monte Carlo (LMC). The convergence properties of the continuous dynamics of OLD, as well as asymptotic and non- -symtoti analyses of its discretizations have been extensively studied (c.g., Roberts et al.
(1996); Villani (2008); Pavliotis (2014); Dalalyan (2017b); Durmus & Moulines (2016); Dalalyan
(2017a); Durmus et al. (2019a b): Vempala & Wibisono (2019); 2019):[Chen, Cheng Bartlett (2018) Dwivedi
et al. (2019): Ma et al. 2019); Chewi al. (202)
Erdogdu & Hosseinzadeh (2011)

---

Under review as a conference paper at ICLR 2022
Meanwhile, the notion of acceleration is less quantified in sampling compared to that in optimiza- tion, although attention has been rapidly building up. Along this direction, one line is based on diffusion processes such underdamped Langevin dynamics (ULD). For example, the convergence and nonasymptotics of discretized ULD have been studied in Cheng et al. (2018); y: Dalalyan & Riou- Durand (2020); Ma et al.] (202T), and were demonstrated provably faster than discretized OLD in suitable setups. These are not only great progresses but also forming perspectives complementary to the extensive studies of the convergence of continuous ULD in the mathematical community (e.g, Mattingly et al. (2002] Cao et al. (2019); Dolbeault et al. (2009) 2015); Villani (2009): |Eckmann & Hairer (2003); Baudoi (2017): Eberle et (2019). Another important line of research related to accelerating paticle-based approaches for optimization in probability spaces (Liu et al. 2019.
TTghvaei & Mehta 2019; Wang & Lii 2019), although we note there no clear boundary between these two lines (e.g.. Leimkuhler et al.] (20T8)) Additional interesting ideas also include (Cen et al. (2018)); Deng et al. (2020) In general, it has been known that adding an irreversible part to the reversible dynamics of OLD accelerates its convergence (c.g., Hwang et al. (2005): Lelievre et al. [213): Ohzeki & Ichiki (2015); RRey-Blllet & Spiliopoulos ( Duncan et al. (2016)), and this work can be viewed to be under this umbrella. Note, though, the discretization of accelerated continuous process is also important, and it will be analyzed.
Specifically, we propose an accelerated gradient-based MCMC algorithm termed HFHR. It is mo- tivated by a simple question: how to appropriately inject noise to NAG algorithm in discrete time, so that it is turned into an algorithm for coomentuu-accelereen to the learning-rate-> 0 limit of NAG (this has been well studied in Ma et al. (202I)) because a finite-step-size discretization of this limiting ODE may not converge as last as NAG with the same learning rate. However, we will still use continuous dynamics as intermediate steps.
sampling? Note we don't 't add noise
More precisely, our first step is to combine existing tools to prepare a non-asymptotic formulation for the later steps. The goal is to better account for NAG's behavior when a finite (not infinitesimal) learning rate is used. As pointed out in (Shi et al. (202T), a low-resolution limiting ODE (Su et al.) 2014). albeit being a milestone leading to new venue of research (e.g. Wiibisono et al. (2016)) does not fully capture the acceleration enabled by NAG - for example, it can't I distinguish between NAG and another momentum method of heavy ball (Polyak] 1964). The main reason is, the low- resolution ODE describes the h - O limit of NAG, but in practice NAG uses a finite (nonzero) h. High-resolution ODE was thus proposed to include additional O(h) terms to account for the finite h effect (Shi et al. 202Tb. The original form of high-resolution ODE involves Hessian of the objective function, which is computationally expensive evaluate and store for high-dimensional problems, but this difficulty can be overcome using techniques introduced in, e.g., Alvarez et al.) 2002 Attouch ct al. 2020) which allows us to derive a Reighheesottton and Hessian- limiting ODE for NAG.
Then we replace the high-resolution term's coefficient in the HFHR ODE by a hyperparameter 2 and then add noise to the resulting ODE in a specific way, which turns it into an SDE suitable for the sampling purpose. This SDE will be termed as HFHR dynamics.
To obtain an actual algorithm, the HFHR SDE then discretized. We will see. both theoretically and empirically, that nonzero a can lead to accelerated convergence of the sampling algorithm; this acceleration is not an artificial consequence of time-rescaling, which would not give acceleration after discretization with an appropriate step size. Meanwhile, note our discretization is just one of the many possible schemes. It was known that high-order discretizations can improve statistical ac- curacy and even the speed of convergence (see e.g., Chen et al. (2015): 5): Li et al. (2019); [Shen & Lee) 2019)), although such improvements often come with more computations per iteration. The dis- cretization considered here is just simple first-order scheme that uses one (full))gradiem evaluation per step, but it better utilizes the structure of HFHR dynamics than Euler- Martermaaaaaa....
e.g.,Chen
Our presentation will be structured as follows. After detailing the construction of HFHR, we will an- alyze its convergence, at both the continuous level (HFHR dynamics) and the discrete level (HFHR algorithm). For precise theoretical results, we will consider the setup of log-strongly-comee
target
distributions, which are commonly considered in the literature (Kim et al., 2016 Bubeck et al. 2018 Dalalyan| 2017b Dalalyan & Riou-Durand 2020 Dwivedi et al., 2019 Shen & Lee| 2019) The additional acceleration of HFHR when compared to ULD in continuous time will be demonstrated
explicitly in Thm For our discretized HFHR algorithm, a nom-syymio error bound will
For imreeesible-aaceeeraaion not from OLD, see e.g., Bierkens 2019 Bouchard-Côté al. 2018
2

---

Under review as a conference paper at ICLR 2022
which confirms that the additional acceleration in continuous time carries through to the discrete territory. Finally, numerical experiments are provided, verifying the validity and tightness of our theoretical results, and empirically showing HFHR remains advantageous for the nonconvex and high-dimensional example of Bayesian Neural Networks.
be obtained (Thm 5.2)
The main contribution of this article is the idea of turning NAG-SC optimization algorithm into a sampler, which also introduces a new dynamics that is neither overdamped or underdamped Langevin. Nevertheless, theoretical analyses (e.g.. Thm Cor & Rmk(5.5) and numerical experiments 6 are also provided to quantify the effectiveness of this idea.
2 BACKGROUND: LANGEVIN DYNAMICS FOR SAMPLING
Consider sampling from Gibbs measure μ whose density is du Je- Te-T (w)) da, where f
Fd + R will be called the potential function. Two diffusion processes popular for sampling (and modeling important physical processes too) are named after Langevin. One is overdamped Langevin dynamics (OLD), and the other is kinetic Langevin dynamies (abbreviated as ULD to comply with convention of calling it underdamped Langevin). They are respectively given by
dq. p.dt
(OLD) dq. -V
(ULD)
dp: -p,tt fl(g;) /2ydB.
where q:, P. E RA, W, B, are i.d Wiener processes in F., and y > is a friction coeffi- cient. Under mild conditions (e.g.. Pavliotis (2014), OLD converges μ and ULD converges to dr(q. = (p) where vpp) (2ir)
/|pll// ---llell so its q marginal follows μ.
OLD and ULD are closely related. In fact, OLD is the y + oo overdamping limit of ULD after time dilation (e.g.. Paaliotti (2014)). However, OLD is a reversible Markov process but UL.D irreversible, and thus both their equilibrium and non-equilibrium statistical mechanics are different, although closely related too. We will only focus on the convergence to statistical equilibrium (see e.g., Souza & Tao (2019) for non-equilibrium aspects).
Many celebrated approaches exist for establishing the exponential convergence (a.k.a. geometric ergodicity) of OLD, including the seminal work of Roberts et al. (1996). the ones using spectral gap (c.g.. Dalalyan 2017b] Lemma 1), synchronous coupling Viillani 22008] p33- ( et al.]
2019h Proposition 1), functional inequalities such as Poincaré's inequality (Pavliotis) [014] The- orem 4.4) and log Sobolev inequality (Vempala & Wibisono 2019 Theorem 1). There are also fruitful results for ULD, including the ones leveraging Lyapunov function (Mattingly et al. 2002 Theorem 3.2), hypocoercivity (Villani 2009 Dolbcault et al. 2009 2015 Roussel & Stoltz 2018
coupling Cheng et al. 2018 Theorem (Dalllya & Riou-Durand 2020 Theorem () al.
2019 Theorem 2.3), LSI (Ma et al. [021 Section 3.1), modified Poincaré's inequality (Cao et al.,
2019
Theorem 1), and spectral analysis (Kozlov) 1989: Eckmann & Hairer 2003).
The study of asymptotic convergence of discretized OLD dates back to at least the 1990s (Meyn af. 1994: Roberts et al. (996). The non-asymptotic analysis of LMC discretization of OLD can be found in Dalalyan (201 7b) and it shows the discretization achieves error, in TV distance, in Ō(422) steps. Subsequent results include ((d/2) in W2 (Dumuss & Moulines 2016). O(d/e) in KI. [Chen & |Bartlet 2018), O(d/e) in W₂ under additional 3rd-order regularity (Durmus et al. 2019b), and (VV4/c) in W2 under additional 3rd-order regularity (Li al. 2021). For discretized
ULD. one has Ö(Vife) iteration complexity in W2 (Cheng et al. 2018 Dalalyan & Riou-Durand,
2020 and Ö(va/vz) in KL (Ma et 2021). ULD is still generally conceived to be advantagcous over OLD and sometimes understood as its momentum--cceleraali version.
3 NOTATIONS AND CONDITIONS
We will use 2-Wasserstein distance to quantify convergence, i.c. W2(11,
(infgen(an)k)) X r11))
& where (2) the set of couplings of μi and jı2-
Assume WLOG that 0 aargmin f(x) The following condition will also be frequently used.
3

---

Under review as a conference paper at ICLR 2022
A 1. (Standard Strong-Convexity and Smoothness Condition) Function f (RR) :R R is m-strronlycoovvx and L-smooth, if there exist constants m, such that vz, y E RR, we have
[V V\(>)) v fly - æ:l and f 2 (() + (V)())) - x) + m
x:
For f C², this condition is equivalent to ml 3 V LI.
4 THE CONSTRUCTION OF HFHR DYNAMICS
HFHR is obtained by formulating NAG-SC Hessian free high-resolution ODE, lifting the high- resolution term's coefficient free parameter, and adding appropriate noises.
More precisely, let's start with NAG-SC algorithm:
0s++ y/x sVy((
yk+1 m1++1 + c(#x+1 - 13k))
where s is the learning rate (also known as step size), and 1-y/ms is a constant based on and the strong convexity coefficient m of f: the method also works for non-strongly-convex f though.
I+yms
high-resolution ODE description of Eq. I & (2 is obtained in (Shi et al. 2021 Section 2)
ÿ + V 2(1-c) + v2 f(v( V + I+c V f(y) = 0, c((+c)
(3)
which can better account for the effect of non-infinitesimal than the - 0 limit (note depends involved, which is expensive to compute and
on 8). However, this original form, Hessian
store especially for high- nigh-dimensional problems.
To obtain a Hessian-free high-resolution ODE description of Equation II and (2. we first turn the iteration into 'mechanical' version by introducing position variable qs y/ and momentum
variable Р (8r, -w/hh. Replacing 0x+1 in ( and the first 0++1 in (2) by qx+1 and Ps+1> the
second x++1 in (2) by qs fVg())) and the 2y in (2) by 9s and Px. we obtain
4k+1 a: + hpa+1 sVf(gn)
Ph+1 = cp, c V)(
Now, choose 7, and h h Ves, h a f.. It is easy to that 0, > 0, then NAG-SC exactly rewrites as
a++1 4s + 11px+1 ha ((q))
(4)
Pk+1
Pk hyPp hV
Note the technique for bypassing the Hessian without introducing any approximation is already well studied in the literature Alvarez al.) (2002); Attoouh et al. (2020)).
So far, both and a are actually determined by the byperparameter s of NAG-SC, However, if we now consider a as an independent variable (i.e., 'lift' it) and let + 0, we see (4) is a 1st-order discretization (with step size h) of the dynamics
= aV p= - V.((g)
(5)
Note a, if inherited from NAG-SC, should be a= -s/c= Ohh), which, in low-resolution ODE will be discarded, and this eventually leads to ULD rather than HFHR. However, we now allow it to be a free parameter and will see that # O(h) can be advantageous.
Before quantifying these advantages later on, we finish the construction by appropriately injecting Gaussian noises Equation (5) This is just like how OLD be obtained by adding noise the
4

---

Under review as a conference paper at ICLR 2022
gradient flow. The right amount and structure of noise turn the ODE into Markov process that can serve the purpose of sampling, and the detailed form of our noise is given by:
dq. (p. + 2adW: dp, (- VS(g,)) )d + v2ndB,
(6)
Here a: 2 0 > O are constant parameters, and W. B. are independent standard Brownian motions in R. This irreversible process will be named as Hessian-Free Higgh-Reolluton (FFHR) dynamics. We write it as HFHR(a, 7) emphasize the dependence
and 7 when needed.
Substitution into Fokker-Planck PDE shows HFHR dynamics unbiased (proof in Appendix Theorem 4.1. I is the invariant distribution of HFHR described in Eq.(6) just like ULD.
5
THEORETICAL ANALYSIS OF THE HFHR DYNAMICS/ DYNAMICS/ALGORRTII ALGORITHM
5.1 HFHR DYNAMICS IN CONTINUOUS TIME
Let's establish the exponential convergence of HFHR dynamics and its additional acceleration when compared to ULD, when the target measure has strongly conveex and smooth potential.
Theorem 5.1. Assume Conditions Al holds and further assume 72 > + m and a I -
Denote the law of by pt. Then there exists K. > 0 depending only on a and % such that
W2(Us, (1) S K'e - +ma)t
W W2(Ho.H1)
Detailed expression of k can be found in Appendix
Thm. 5.1 state that HFHR converges to the target distribution exponentially fast in log-strongly- concave and smooth setup. For ULD, Dalalyan & Riou- Burand 2020 Theorem 1) obtained ex- ponential convergence result in 2- Wasseesteei distance with rate VK+VE= using a simple and ele- gant coupling approach, and showed this rate is optimal as it is achieved by the bivariate function f(z v) = - + y. In Thm 5.1 we use the same coupling approach to obtain an (asymptot- ically) equivalent rate 2yr when a 0 and shows faster convergence for > O, the acceleration of HFHR in continuous time is evidenced. For example, if we set y = 2vL and push to the upper bound specified in Thm. 5.1 we obtain an O(VL) rate in the log-strongly-conee Riou-Durand 20200), this is speed-up of order K.
mm
ym Since ULD is HFHR(0,7) and our bound agrees with existing result
setup. Compared with the rate in Dalalyan &
5.2 HFHR ALGORITHM IN DISCRETE TIME
To obtain implementable method, we now discretize the time of HFHR dynamics. As our main goal is to show the acceleration enabled by a won't disappear after discretization (unlike a fake acceleration due time rescaling), we'll just analyze a Ist-order discretization (but high-accuracy discretization adapted from RMA (Shen & Lee 2019) will also be provided, in Appendix F.
For simplicity, we'll work with constant step size h. Inspired by Strang splitting for differential equations (Strang 1968 McLachlan & Quispel 2002). we consider symmetric composition for
updating over each time interval [kh, (k + 1)hl: #h++1 := Φ² thth o & (æx) where 2 and di correspond to solution flows of split SDEs, respectively given by
4sh Pun.
à: dq pdl
lg (QV(( 22rdW
v:
dp "ppdt
2
dp dt
and o² (xx) and (i'(() (xx) mean z's value after evolving and vb for t time with initial condition Zo- Note that & flow can be solved explicitly since the second equation is an -nnttinnUnnleeeckk process and integrating the second equation followed by integrating the first one gives explicit solution
5

---

Under review as a conference paper at ICLR 2022
4: 90 + 1-e Po 22 So
1- dB(s).
e "P + 22 So 7(t
(8)).
For an implementation of the stochastic integral part Equation denoting
v lo 1-e 1--- ds(s) and Y 22 e (->(t-=) B (s), and the covariance matrix of (X Y)
7
2h446 - 11
(ov(X,X) Y)
As mean and covariance fully determine
(1-6 7
(1 (1-e-th
X Y
a Gaussian distribution,
.M where M is the Cholesky decomposition of Cov( (X, Y). is
2d standard Gaussian random vector, i..d. at each step, and of can thus be exactly simulated.
However, & flow is generally not explicitly solvable unless f is quadratic function in We simply choose approximate w/a (3ro) with one-step Euler-Maruyama integration v/th (xo) R3
qn do aVV 22hn
A (#o) given by
where is standard d-dimensional Gaus-
Pn. - Po ((go)
sian random vector, again i.i.d. each time called.
Altogether, one step of an implementable Strang's splitting of HFHR hence of T o σ² and we call this numerical scheme the HFHR algorithm, which is summarized in Algorithm
Algorithm HFHR Algorithm
Input:
potential function and its gradient Vj, damping coefficients a and 7, step size h,
initial condition (g) Po)
2: procedure DISCRETIZED HFHR(f, Vf, Po)
and initialize
qo
Po
4: while not converge do
5: Generate independent standard Gaussian random vectors 7641 Rdd 1 Rd
1-0-24
6: Run ob
91
4sh
7
PRR
+ +ME++1
P1
1. Pkh
7:
Run Th
92 P2
=
91 aV aV )h 221hm411
P1 - vf(gg))
+ 1-c-73
8:
Run of
[3
92
7 P2
+ M
P3
e P2
(1=+1)A
433
9:
P(K+1)A
P3
+ k
11:
end while
12:
end procedure
As v in Strang splitting is replaced by stt-rrder approximation vi the method is of order 1, however with good constant. This is rigorously established by the following theorem (interested readers referred to Appendix D.7 and Li et al. (2021) for more technical details):
Theorem 5.2. Under Assumption we further assume L+m > and the function V
satisfies a third-order growth condition, i.e.. iee.| G 1+ Vq E R&d for some G 0.
If (go- Po) ~ To. then there exists ba, C > 0 such that when < ho. we have
Wa(ue, n) S -2k''-1" +ma)kh tma (2 π) + /2Ch.
where re' a constant depending only on L,m, 7, a (see Appendix f detail), bhg is the law of the q marginal the kithh iterate in Algorithm and p is the marginal of the invariant distribution m.

---

Under review as a conference paper at ICLR 2022
In particular, C' O(vd) and there exists b > 0, independent of a and is of order O
b (a² m
a
C
7?
Remark 5.3. The linear growth (at infinity) condition VDf is actually not restrictive as appears. For example, for monomial potentials, i.e., f(x) = x², E Z,, our linear growth condition is met when p S 4. whereas standard condition (Paviionis 2014 Theorem f.r) the existence of SDE solutions holds only when S 2. In addition, our condition is related to the Hessian Lipschitz condition commonly used in the literature (e.g.. Durmus al. (2019b); |a et al. (202T). Smoothness and Hessian Lipschitzness imply our condition. Meanwhile, examples that satisfy linear growth condition but are not Hessian Lipschitz exist, e.g. f(x) x and thus linear growth condition is not necessarily stronger than Hessian Lipschitzness.
Inspecting the role of a in Equation 8 we see that a clearly increases the rate of exponential decay, but at the same time it can also increase the discretization error (see (9); assuming h is fixed). However, as the following Corollary 5.4 and its remark will show, the net effect having a positive 0, at least for some a", is reduced iteration complexity.
If (qu:Pa) ~ To. then there exists ha, C > 0 (same as that in Theorem 5.2 recall C = O d.) such that for any target error tolerance
Corollary 5.4. Consider the same assumption as in Thm. 5.2
> 0, if we choose h = ht min ho, 2v25 }. then < 2/2Cho, after
C log
2V2n' W2( 2/2k'W2(TO..
vd
k* == 2
Ô
(10)
+ max
steps, we have W2(URR 11) S.
Remark 5.S. Recall from Thm 5.2 that C b (a² 77 so if we consider the minimizer a' of an upper bound of m m+maa a* argmin m +
a² -+ 23-1 This suggests that by choosing
C
an optimal a > 0, one could effectively reduce iteration complexity. Note, however, that this a may not be the true optimal one as bounds may not be tight. If they were, k;, (2v33 3)k² R 0.46k iee. steps needed by ULD (discretized by Alg 7ith a 0) can be halved by HFHR.
Rmk 5.5 shows HFHR algorithm can lead to similar bound iteration complexity as ULD algo- rithm but with an improved constant, and thus a more efficient algorithm. This improvement also shows that the acceleration of HFHR can be carried through from continuous time to discrete time. The same conclusion has been consistently observed in numerical experiments too.
Remark 5.6. Readers interested in more explicit condition number dependence are referred Ap- pendix E where we show, for 2D Gaussian target with condition number k. >>> 1, the convergences of Euler discretization of ULD under optimal parameters and HFHR under suboptimal parameters are, respectively, given by (1 - 1/k + o(1/x))) and (1 - 2/N o(1/k))", where is the number of iterations. The latter (HFHR) is faster despite that its hyperparameters may not be optimal.
NUMERICAL EXPERIMENTS
We now empirically study the acceleration enabled by a * 0 by comparing HFHR algorithm and the popular KLMC discretization of ULD (Dalalyan & Riou-Durand 2020). For fairness, discretiza- tions of the same order and cost are compared (Appendix has an additional comparison).
6.1 VERIFICATION OF THEORETICAL RESULTS IN SECTION 5.2
This subsection numerically verifies the and h dependence of HFHR algorithm as well as the genuine acceleration enabled by For this purpose, we will not use Gaussian targets, because otherwise HFHR will decouple across different (orthogonal) dimensions, and hence its discretization error having a ((vd) dependence would just be a consequence using 2- 2-Wasserstein distance for quantifying statistical accuracy. To inspect a more interesting example, consider a potential that no longer additive across different dimensions, namely f(x) log (e +... edd) + |||=||?
7

---

Under review as a conference paper at ICLR 2022
It's not hard to see that all dimensions will be coupled in the HFHR dynamics (and ULD too). Moreover, the new potential f is still a strongly convex function and satisfies the assumption in Theorem 5.2 When the target measure is non-Gaussian, we no longer have closed form expression for 2- Wasserstein distance and it is computationally expensive to approximate the 2- Wasserstein distance by samples. Therefore, we use the error of mean instead surrogate because 9 q < (2 μ1) and hence the bound in Equation (8) also applies to the error in mean, and so does the iteration complexity bound in Equation (10).
Theorem 5.2 says the final sampling error is upper bounded the discretization error that is linear in h and vd. To numerically verify the linear dependence on h. we work with d' = 2 and ran ULD 0.0005) to obtain 10 independent realizations and use them (as benchmark) to empirically estimate E. We then set y 2,
algorithm for sufficiently long time with tiny step size
and use Algorithm| [ with h E (2*1 - < k < 0}. For each h, we run F (with T = 50) iterations in Algorithm| [ to ensure the Markov chains are well-mixed and the contribution to final error from smaller than discretization error. The results are plotted in
exponential decay is
Figure
la
The observed linear dependence on the step size h is consistent with Theorem 5.2
Euxili
HFHR Mgonttin
Euxila
NETHR Algeritmmm
of ftors to mach e-dmssnnsss
Agoethm Agorthm
Step size
Dimension
(a) Linear dependence of (h) Linear dependence of discretization of Al- discretization error of Al- gorithm on
gorithm on vd
Figure I: Illustration of the consistency between Figure 2: Improvement of Algorithm | the theoretical bound in Theorem 5.2 and experi- over ULD algorithm in iteration com- ment results.
plexity. (vertical bar 1 std.)
To numerically verify the O(vi) dependence, we extensively experiment with E 21, 5, 10, 20, 50, 100, 200, 500, 1000]. For each d, we run 1,000 independent realizations of ULD algorithm until well converged with tiny step size (h- = 0.005) and use their empirical average as the 'true' mean. Then we fix 2, a = 1, h 0.1, T 10 and for each d, we run 1,000 independent realizations of HFHR algorithm for T. = 100 iterations. Experiment results are plotted in Figure Ib and the linear trend demonstrates that the bound in Thm. 5.2 is tight in terms of d-dependence.
The final experiment compares Algorithm I with ULD algorithm in terms of iteration complex- ity. The goal is to demonstrate the genuine acceleration of HFHR is not an artifact due to time rescaling, which would disappear after discretization as the stability limit changes accord- ingly. To do so, we push both ULD and HFHR to their respective largest values that still allow monotonic convergence at a large scale, and compare their mixing times. For general nonlinear problems like the one here, Remark [..5 suggests that with appropriately chosen a, HFHR algorithm effectively reduces the constant factor of the iteration complexity, implying ac- celerated sampling. To just provide one empirical verification of this improvement, we choose d = 10 and use the error of mean || - E,gg to measure sampling accuracy. The bench- mark, i.e., E, is again obtained from 1,000 independent realizations of ULD algorithm with tiny step size (h = 0.005), ran for long enough to ensure the corresponding Markov chain is well- mixed. The initial measure is chosen as a Dirac measure at (100 x 1,, Oa), where 1,, O₂ are d-dimensional vectors filled with 1 and 0 respectively. We pick threshold a E {0, 0.001, 0.002, 0.005, 0.01, 0.02, 0.05, 0.1, 0.2, 0.5, 1, 2, 5, 10, 20, 30, 1003. we try all combi- nations of (7. E 1.. 0.5, 1, 2, 10, 20, 50, 100) x {0.1 x [50]] for Algorithm (we also 0). and empirically find the best combination that requires the fewest
0.1, and for each
run ULD algorithm when iterations to meet | E, g| < e. We find that = 5 already surpasses the stability limit of ULD algorithm, hence the range of step size covers the largest step size that practically useable for ULD algorithm. Experiments are repeated with 100 different seeds to further reduce variance.
The results are shown in Figure When algorithm (note it also does so when a = because HFHR uses a siieency-wis--compparr but more refined discretization than ULD algorithm). In particular, when 0.5 and 1, which
0, HFHR algorithm consistently outperforms ULD
8

---

Under review as a conference paper at ICLR 2022
empirically the best values found for this experiment, HFHR algorithm achieves the specified e-closeness nearly 6x times faster than ULD algorithm, and its decreased mixing time (compared to 0 for the same algorithm) is consistent with the se 0.46 factor in Rmk 5.5). This empirical study corroborates that the acceleration HFHR dynamics creates also carries through to its discretization, and the acceleration of HFHR algorithm over ULD algorithm can be significant.
6.2 BAYESIAN NEURAL NETWORK
Now consider Bayesian neural network (BNN) which is a compelling learning model (Wilson, 2020]; however, the focus won't be on its learning capability, and instead we just consider its train- ing, which amounts to a practical, high-dim., multi- example of sampling tasks. It no longer satisfies the conditions of our analysis, and our goal is to show HFHR still accelerates. We use fully- connected network with [22, 10, 2] neurons, ReLU, standard Gaussian prior for all parameters, and compare ULD and HFHR on data set Parkinson from UCI repository (Dua & Graff 2017).
Choices of hyper- syper-parameter for Algorithm I and ULD algorithm systematically investigated. For each pair (7, a) E {0.1, 0.5, 1,5, 10, 50, 100]?, we empirically tune the step size to the stability limit of ULD algorithm, simulate 10,000 independent realizations, and use the ensemble to conduct Bayesian posterior prediction. HFHR will then use the same step size. For each y, we plot the negative log likelihood of HFHR algorithm (with different a choices) and ULD algorithm on training and test data in Figure 3 Cases where a is too large for numerical stability are not drawn.
From Figure 3 we find that HFHR converges significantly faster than ULD in a wide range of setups. In general, the log-strongly-coneee multimodal target distributions. However, this numerical result shows that HFHR still accelerates ULD for highly complex models such as BNN, even when there is no obvious theoretical guarantee. showcases the applicability and effectiveness of HFHR as a general sampling algorithm.
assumption required in Theorem 5.2 does not hold for
12
Medi
Moalla Lat
7 0.1 (h 0.005)
(b) 7 = (h- 0.01)
(c) 7 = 10 (h- = 0.05)
(d) 7 = 100 (h = 0.1)
Mocile ie 6.1
wilile 0.2
Mar Dess
suvie of Paase
(c) 0.1 (A- 0.01)
(h- 0.02)
gg)
10 0.1)
(h) 100 (h- = 0.2)
Figure 3: Training Negative Log-Likelihood (NLL) for various Row I: step sizes are close to the stability limit of ULD algorithm; Row 2: further increased step size exceeds that stability limit.
7 CONCLUSION AND DISCUSSION
This paper proposes HFHR, an accelerated gradient- based MCMC method for sampling. To demon- strate the acceleration enabled by HFHR, the geometric ergodicity of HFHR (both the continuous and the discretized versions) is quantified, and its convergence is provably faster than Underdamped Langevin Dynamics, which by itself often already considered as an acceleration of Overdamped Langevin Dynamics. As HFHR is based on new perspective, which is to turn NAG-SC opti- mizer with finite learning rate into a sampler, there are a number of interesting directions in which this work can be extended. Besides further theoretical investigations that aim at refining the error bounds, examples also include the followings: to scale HFHR up to large data sets, full gradient may be replaced by stochastic gradient (SG)- how to quantify, and hence optimize the perfor- mance of SG- SG-HFHR? Can the generalization ability of HFHR trained learning models (c.g.. BNN) be quantified, and how does it compare with that by LMC and/or KLMC? These will be future work.

---

Under review as a conference paper at ICLR 2022
REFERENCES
Felipe Alvarez, Hedy Attouch, Jérôme Bolte, and Patrick Redont. A second-order gradient-like dissipative dynamical system with hessian-driven damping.: Application to optimization and me- chanics. Journal de mathématiques pures et appliquées, 81(8):747-779, 2002.
Hedy Attouch, Zaki Chbani, Juan Peypouquet, and Patrick Redont. Fast convergence of inertial dynamics and algorithms with asymptotic vanishing viscosity. Mathematical Programming, 168 112)::13--17,, 2018.
Hedy Attouch, Zaki Chbani, Jalal Fadili, and Hassan Riahi. First-order optimization algorithms via inertial systems with hessian driven damping. Mathematical Programming, pp. 1-43, 2020.
Fabrice Baudoin. Bakry-emery meet villani. Journal of Functional Analysis, 2017.
GNJC Bierkens, Paul Fearnhead, and Gareth Roberts. The zig-zag process and super-efficient sam- pling for bayesian analysis of big data. Annals of Statistics, 47(3), 2019.
Alexandre Bouchard-Côté, Sebastian J Vollmer, and Arnaud Doucet. The bouncy particle sam- pler: A nonreversible rejection-free frre markov chain monte carlo method. Journal of the American Statistical Association, 113(522):855-885 2018.
Sébastien Bubeck, Ronen Eldan, and Joseph Lehec. Sampling from a log-concave distribution with projected langevin monte carlo. Discrete Computational Geometry, 59(4):757-773, 2018.
Yu Cao, Jianfeng Lu, and Lihan Wang. On explicit 12 converegeee rate estimate for underdamped langevin dynamics. arXiv preprint arXiv:1908. 0444:16908.04766. 2019.
Changyou Chen, Nan Ding, and Lawrence Carin. On the convergence of stochastic gradient MCMC algorithms with high-order integrators. NIPS, 2015.
Changyou Chen, Ruiyi Zhang, Wenlin Wang, Bai Li, and Liqun Chen. A unified particle- optimization framework for scalable bayesian sampling. In The Conference on Uncertainty in Artificial Intelligence, 2018a.
Yi Chen, Jinglin Chen, Jing Dong, Jian Peng, and Zhaoran Wang. Accelerating nonconvex learning via replica exchange langevin diffusion. In International Conference Learning Representa- tions, 2018b.
Xiang Cheng and Peter L Bartlett. Convergence of langevin meme in kl-divergence. PMLR 83, (83):186-21)) 2018
Xiang Cheng, Niladri Chatterji, Peter I. Bartlett, and Michael Jordan. Underdamped langevin meme: A non-asymptotic analysis. Proceedings of the 31st Conference On Learning Theory, PMLR, 2018.
Sinho Chewi, Chen Lu, Kwangjun Ahn, Xiang Cheng, Thibaut Le Gouic, and Philippe Rigollet. Optimal dimension dependence the metropolis-adjusted langevin algorithm. COLT, 2021.
Lenaic Chizat and Francis Bach. On the global convergence of gradient descent for over- parameterized models using optimal transport. In Advances in neural information processing systems, pp. 3036-300 3046, 2018.
Arnak Dalalyan. Further and stronger analogy between sampling and optimization: Langevin monte carlo and gradient descent. Conference on Learning Theory, pp. 678-689. PMLR, 2017a.
Arnak S Dalalyan. Theoretical guarantees for approximate sampling from smooth and log-concave densities. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 79(3): 651-676, 2017b.
Arnak S Dalalyan and Lionel Riou-Durande On sampling from log-concave density using kinetic Langevin diffusions. Bernoulli, (6(3):1956-1988 1988, 2020.
10

---

Under review as a conference paper at ICLR 2022
Wei Deng, Qi Feng, Liyao Gao, Faming Liang, and Guang Lin. Non-convex learning via replica exchange stochastic gradient memc. In International Conference on Machine Learning, pp. 2474 2483. PMLR. 2020.
Jean Dolbeault, Clément Mouhot, and Christian Schmeiser. Hypocoercivity for kinetic equations with linear relaxation terms. Comptes Rendus Mathematique, 347(9-0)):51-1511 2009.
Jean Dolbeault, Clément Mouhot, and Christian Schmeiser. Hypocoercivity for linear kinetic equa- tions conserving mass. Transactions of the American Mathematical Society, 367(6):3877-3888 2015.
Dheeru Dua and Casey Graff. UCI machine learning repository, 2017. URL http://arccive. ics. uci Ledu/mi /ml
Andrew B Duncan, Tony Lelievre, and GA Pavliotis. Variance reduction using nonreversible langevin samplers. Journal of statistical physics, 163(3):457-44) 491, 2016.
Alain Durmus and Eric Moulines. Sampling from strongly log-concave distributions with the unad- justed langevin algorithm. arXiv preprint arXiv:1605.01559. 5, 2016.
Alain Durmus, Szymon Majewski, and Blažej Miasojedow. Analysis of langevin monte carlo via convex optimization. The Journal of Machine Learning Research, (M(1):2666-271) 2711, 2019a.
Alain Durmus, Eric Moulines, et al. High-dimensional bayesian inference via the unadjusted langevin algorithm. Bernoulli, 25(4A):2854 2882, 2019b.
Raaz Dwivedi, Yuansi Chen, Martin J Wainwright, and Bin Yu. Log-concave sampling: Metropolis- hastings algorithms are fast. Journal of Machine Learning Research, 20(183); 2019.
Andreas Eberle, Arnaud Guillin, Raphael Zimmer, et Couplings and quantitative contraction rates for langevin dynamics. The Annals of Probability, 47(4): 2010, 2019.
J-P Eckmann and Martin Hairer. Spectral properties of hypoelliptic operators. Communications in mathematical physics, (35(2):233-255) 2003.
Murat A Erdogdu and Rasa Hosseinzadeh. On the convergence of langevin monte carlo: The inter- play between tail growth and smoothness. COLT, 2021.
Charlie Frogner and Tomaso Poggio. Approximate inference with wasserstein gradient flows. In International Conference on Artificial Intelligence and Statistics, 2020.
Ye He, Krishnakumar Balasubramanian, and Murat A Erdogdu. On the ergodicity, bias and asymp- totic normality of randomized midpoint sampling method. Advances Neural Information Pro- cessing Systems, 33, 2020.
Bin Hu and Laurent Lessard. Dissipativity theory for nesterov's accelerated method. In Pro- ceedings of the 34th International Conference Machine Learning-Volie Volume 70, pp. 1549-1557. JMLR.org, 2017.
Chii-Ruey Hwang, Shu-Yin Hwang-Ma, Shuenn-Jyi Sheu, et al. Accelerating diffusions. Annals of Applied Probability, 15(2):1433-1444, 2005.
Richard Jordan, David Kinderlehrer, and Felix Otto. The variational formulation of the fokker- planck equation. SIAM journal on mathematical analysis, 29(1):1-17, 1998.
Arlene KH Kim, Richard Samworth, et al. Global rates of convergence in log-concave density estimation. The Annals of Statistics, 446)))2756-2779 2016.
S. M. Kozlov. Effective diffusion the fokker-planck equation. Mathematical notes of the Academy of Sciences of the USSR, 45:360-368, 1989.
Benedict Leimkuhler, Charles Matthews, and Jonathan Weare. Ensemble preconditioning for markov chain monte carlo simulation. Statistics and Computing, 2872),277-290, 2018.

---

Under review as a conference paper at ICLR 2022
Tony Lelievre, Francis Nier, and Grigorios A Pavliotis. Optimal non-reversible linear drift for the convergence to equilibrium of a diffusion. Journal of Statistical Physics, 152 2013.
Ruilin Li, Hongyuan Zha, and Molei Tao. Sqrt (d) dimension dependence of langevin monte carlo.
arXiv preprint wixv:2209.03839, 2021.
Xuechen Li, Denny Wu, Lester Mackey, and Murat A Erdogdu. Stochastic Runge-Kutta accelerates Langevin Monte Carlo and beyond. NeurIPS, 2019.
Chang Liu, Jingwei Zhuo Pengyu Cheng, Ruiyi Zhang, and Jun Zhu. Understanding and acceler- ating particle-ban baseed variational inference. International Conference on Machine Learning, pp- 4082-4092, 2019.
Qiang Liu and Dilin Wang. Stein variational gradient descent: A general purpose bayesian inference algorithm. Advances in neural information processing systems, pp. 3378-2386, 2016.
Yi-An Ma, Yuansi Chen, Chi Jin, Nicolas Flammarion, and Michael Jordan, Jordan. Sampling can be faster than optimization. Proceedings of the National Academy of Sciences, 116(42):20811-2088 2019.
Yi-An Ma, Niladri Chatterji, Xiang Cheng, Nicolas Flammarion, Peter Bartlett, and Michael I Jor- dan. Is there an analog of nesterov acceleration for memc? Bernoulli, 2021.
Jonathan C Mattingly, Andrew M Stuart, and Desmond Higham. Ergodicity for sdes and ap- proximations: locally lipschitz vector fields and degenerate noise. Stochastie processes and their applications, 101(2):18: ((12):85-2333 232, 2002.
Robert McLachlan and G Reinout Quispel. Splitting methods. Acta Numerica, 11:341, 2002.
Sean P Meyn, Robert L Tweedie, et al. Computable bounds for geometric convergence rates of markov chains. The Annals of Applied Probability, 4(4):981-1011, 1994.
Yurii Nesterov. A method for unconstrained convex minimization problem with the rate of conver- gence o (1/k2 2). In Doklady AN USSR, volume 269, pp. 543-547, 1983.
Yurii Nesterov. Introductory lectures on convex optimization: A basic course, volume 87, Springer Science & Business Media, 2013.
Masayuki Ohzeki and Akihisa Ichiki. Langevin dynamics neglecting detailed balance condition. Physical Review E, 92(1):012105, 2015.
Grigorios A Pavliotis. Stochastic processes and applications: diffusion processes, the Fokker-Plann Planck and Langevin equattons, volume 60. Springer, 2014.
Boris Polyak. Some methods of speeding up the convergence of iteration methods. USSR Com- putational Mathematics and Mathematical Physics, 4(5):1-17, 1964.
Luc Rey-Bellet and Konstantinos Spiliopoulos. Inreversible langevin samplers and variance reduc- tion: a large deviations approach. Nonlinearity, 28(7):2081, 2015.
Gareth O Roberts, Richard L Tweedie, et al. Exponential convergence of langevin distributions and their discrete approximations. Bernoulli, 2(4):341-363, 1996.
Julien Roussel and Gabriel Stoltz. Spectral methods for langevin dynamics and associated error estimates. ESAIM: Mathematical Modelling and Numerical Analysis, 52(3):1051-1033 2018.
Ruoqi Shen and Yin Tat Lee. The randomized midpoint method for log-concave sampling. In Advances in Neural Information Processing Systems, pp. 1098-2109, 2019.
Bin Shi, Simon S Du, Michael Jordan, and Weijie J Su. Understanding the acceleration phe- nomenon via high-resolution differential equations. Mathematical Programming, pp. 1-70, 2021.
Andre N Souza and Molei Tao. Metastable transitions in inertial langevin systems: What can be different from the overdamped case? European Journal of Applied Mathematics, 30(5):830-852, 2019.
12

---

Under review as a conference paper at ICLR 2022
Gilbert Strang, On the construction and comparison of difference schemes. SIAM journal merical analysis, 533):506-5111 517, 1968.
Weijie Su, Stephen Boyd, and Emmanuel Candes. A differential equation for modeling nesterov's accelerated gradient method: Theory and insights. In Advances Neural Information Processing Systems, pp. 2510-2518, 2014.
Amirhossein Taghvaei and Prashant Mehta. Accelerated flow for probability distributions. In Ka- malika Chaudhuri and Ruslan Salakhutdinov (eds.), Proceedings of the 36th International Con- ference Machine Learning, volume 97 of Proceedings of Machine Learning Research, pp- 6076-6000 6085- Long Beach, California, USA, 09-15 Jun 2019, PMLR.
Santosh Vempala and Andre Wibisono. Rapid convergence of the unadjusted langevin algorithm: Isoperimetry suffices. In Advances in Neural Information Processing Systems, pp. 10092-1044 2019.
Cédric Villani. Optimal transport: old and new, volume 338. Springer Science & Business Media, 2008.
Cédric Villani. Hypocoercivityy Memoirs of the American Mathematical Sociery, 202(950), 2009.
Yifei Wang and Wuchen Li. arXiv:1909.02102. 2019.
Accelerated information gradient flow.
arXiv preprint
Andre Wibisono. Sampling as optimization in the space of measures: The langevin dynamies composite optimization problem. In Conference On Learning Theory, pp. 2093-3027, 2018.
Andre Wibisono, Ashia C Wilson, and Michael Jordan. A variational perspective accelerated methods in optimization. Proceedings of the National Academy of Sciences, 113(47):E7351- E7358, 2016.
Andrew Gordon Wilson. The case for bayesian deep learning. arXiv preprint arX19:2001.10995, 2020.
Ashia C Wilson, Ben Recht, and Michael Jordan. A Iyapunov analysis of accelerated methods in optimization. Journal of Machine Learning Research, 22 13):1-34. 2021.
Ruiyi Zhang, Changyou Chen, Chunyuan Li, and Lawrence Carin. Policy optimization as wasser- stein gradient flows. In International Conference on Machine Learning, pp. 5737-5746, 2018.
13

---

Under review as a conference paper at ICLR 2022
A
ADDITIONAL NOTATIONS
We introduce a few notations that are used in the main text as well as some proof. When Vf is HFHR dynamics also L'-Lipschitz, as proved in
a (q)) P f((g)
L-Lipschitz, the drift term
Lemma D.3
where
1 + max (/2 ]. vI v
L' max
We show in Lemma D.5 that a finear-transformed HFHR dynamics satisfies the nice contraction property, the linear transformation P we defined as
7 0
P
VI+ayl
ER 2dx2d
Denote the largest and the smallest singular value of P by
a7
72
Va222 227 + 4ary + 7p
mmax
1,
2
ay
72
Va2g² 207 + 4ay + grt +4
0 min
+
+1
2
2
and its condition number by
Vazyz 22 +1ay
K
0max
+
miin
07
Va2ga 4
+
The rate X' of exponential convergence of transformed HFHR dynamics is characterized Lemma D.5 and is defined as
am,, -2- L
X min
7
given that 72 > L.
B PROOFS FOR THE CONTINUOUS DYNAMICS
Notations and definitions can be found in Sec
B.1 PROOF OF THEOREM
Proof. The Fokker-Plank equation of HFHR is given by
aept = Vv (1-P/10) P: +a (Vf(g)p) + Agpi) (Vf())(V
+ (Vp (pp:) + Appt)
where Va = (Va: Vp). For x e -(()))))- 1(g)- bllel? we have
V (1
- (V( ([-P/() x) - - Van) =0,
p
=
A = V, (TVS(g))
APT -V (TD)
Therefore ब, 0 and hence is the invariant distribution of HFHR.
14

---

Under review as a conference paper at ICLR 2022
B.2 PROOF OF THEOREM
Proof. Consider two copies of HFHR that are driven by the same Brownian motion
[V.() -P:==(PpP(( dq. -
+ + V
dq. (p: VV 2 dp. = (->P, Vf(g)) v2jdB
dp; = (-MP: VV(q:))) 2
where we (āo: Po) n. Po - Pu and qo such that
W (mo, p)) EE Ilao 8o1l3 qu po
Φ 4/ on a, 7, have
9: 9: PP: P.
Denote
P
where is defined Appendix A
By Lemma D.5 and the assumption
Φ
2( +maet
Φ
Therefore obtain
W2 (Pr: (4)
inf
Elqg: - q.l12
2
inf (9. 1.)
E
q: - q1
p: P:
<B| P- 112
<B111-113 2 +majr
1
[440
-0 <(6'))e-2115) +ma)
90 90
Po
Po
20 - +ma) W2(HH) 11)
Taking square root yields the desired result.
C ARBITRARY LONG TIME DISCRETIZATION ERROR OF ALGORITHM
Theorem C.I. Under Conditions A further assume the function V grows at most linearly, iee, [ S / + Vq E Rdd, Also suppose y in HFHR dynamics satisfy 722 L. Then there exist C, ho O, such that for 0 < h S bo. we have
E|| =x1?) Ch
where 2y is the k-th iterate of Algorithm 7 with step size h starting from To. I& is the solution of HFHR dynamics at time khh starting from Po. This result holds uniformly for all k and k can go to sc. In particular, C= O(Vd) and if Ltm 2 ma, then there exists b > 0, independent of
and is of order O (Vd), such that
b
a
1
C
(11)
Proof. Denote t4 = khh, the solution of the HFHR dynamics at time t by Zp (t), the k-th iterates of the Strang's splitting method of HFHR dynamics by To (kh).. Both mo (1) and Ta (ksh)
15

---

Under review as a conference paper at ICLR 2022
start from the same initial value zu. The linear transformation P defined in Appendix A transforms the solution HFHR dynamics into Pa (t) Pro, (t) and the Strang's splitting discretization of HFHR into Par (t) = P (t).
For the ease of notation, we write yo.u. (4x) as Uv and yo.u, (1x) as Dr- We have the following identity
y 12. (h) -y.. 32 (ti]l* =B
9 ykt1
Dx+1
9) Vv (hi Yt.. v2 (1t) + yl, De. (h) Vi.. (hill'
2 + IE
Vie.ve (l)) Vt Rk (1t)
Vce.8e (h) yes Mx (h))
+2B E Uix.tr (h)
Des. (h) Uon.D (h), Ve. (1)
By Lemma D.5
when 0 h < 2x term 1
1
can be upper bounded
Vta.ur (h) - Ve (n))
< 2X' + =(22)22²)
< (1 X'h) Blyya 5z
where the second inequality is due to --
< x 2 Vx > O.
For term we have by Lemma D.8 that
B Vor.Ds (h) I1. de w/|l
S oi Bs (h) Bry.os (A)l < a Ch
where Tmax is the largest singular value of matrix P.
For term
have by Lemma D.1 that
2F Vinunr (h) Vcc ur (h); Yea Nr. (h) Dee vr (n)
-2k V Ms z, Dic. 8z (h) V Ven.D. (11)
2E
Uı Uu Uu.. (h) Vo.. (h) 2E 2, F...u.. (h)
100--.r..00 Ver. (h)
3a
3b
16

---

Under review as a conference paper at ICLR 2022
For term
3a by the tower property of conditional expectation, we have
21E
yk
yk. Mr (1) - Dr (1i) -2
R
Uv
Uu Vr. (1))
Der (hi) FR
=2
yk
Dc R
Vc..8e (?)
Der.Dr (1t) Fk
2
22
1E
11ll
E
y 112.81 (h)
Dea 20.00 Br (h)
(h) FR
22 19 1 22ll
2 o max
EF
E
Besida (1))
De (1) Fx
52 12 - Bull sil/voo Vonan GP
220 Ci CiV BBlly U/ 12?.
For term 3b
when
h 17.7 we have by Lemma D.I and Lemma D.8
2E 2, Dr. (h)
De.. (1)) 22 B| 2|
Vta.Wn (1)) - Dc. (n))
2
-2 =|
E E
Bec.Du (h)
Ven.in (h)
-2V 1|
E
Deu (h)
De 22 Fi
<20max
ČBllya Uxll" 12 v C3ss
$20 C2 C.VOVTI El
where ő = 21,) = 2(x))²
(I.)² is from Lemma D.1 and Lemma D.3
Recall both C₁ and C2 depend on || and we would like to upper bound this term. To this end, con- sider ie)), solution of HFHR dynamics with initial value To that follows the invariant distribution æt and realizes W2(To, n), i.e., Eliao 20ol² (TD,
W3 T).
17

---

Under review as a conference paper at ICLR 2022
(Blyy>- FR 1) 2
Denote žt = =(kh)) and ek
we then have
E xx|1
=E p + 3t. æx|
S2B xx 2 + 21 By 2x 1|
<4E|
+ 4E|
bt xxl 21 + 2E 3tK xx 11
-45| mxe|"
4E
p- 1--1((0x xx))
+ +2E p- P(ox xzx)
du + +d)
+
E
P(#x - 2x>i|? +
2
E E 3 11
min
24
(LL
4 22
2
||ql? dh +
+
el
min
0
(fL
2 ox
<4
||q1| du + d)
Wx2 (To, T)
72
Rd
min
Fe +G
where (i) is due to Lemma
D.5 Recall from Lemma
D.8
we have
C S A. V B| 2 + +B, S AIVFE + (AIVG+ Bi) _ Uies + Vi
C2 S A2V || 2 + B2 s A2V + (A2VG + BB) 2 Uzes + V2
where
A - G) max 1.25, + )
Bi = + L 1| 0.5a + (1.26/ L.14aaa + 27) Vhd
42 +
B, max(a + + 3.34
18

---

Under review as a conference paper at ICLR 2022
Combine the above and bounds for terms 2 2.2
3a and 3b
we then obtain
224 5(1 Xhjee + mas C2 20 m Cregh + 20 mar C2
<(1 Xhhee + max 2(U3et + V/2 n² 3 + 20x (UIEH + ViJech² 20x +
Xhh + 20 U3h2 + max U2h + 20 U2 The ez
max 1/2 + max 122 V2Vch Ch exli 2 + 20 max 1/21² +
Xh + 20 U2h + u2h + 20 U2 V Che ee
X
het
8
2 m2x V1 20 max 1/2 v2 2
+ 20 max Vi2n3 3
+
X
1 7 N'h + 20 max U3h + 20 max + 20 1/2 vchă ex
8
2 max V1 + 20 max V22 VVhhh) 2
+
X
+ 20 V22
max 2
h
2 12x V1 + 20 max V22 Ch
(
1
5(1 Xh)ez
20 V/2 3
2
X
(1 Xh)e + Kh
2
where (i) is due to h min 22 ha) and
VX
h1 2/n' L max + 1.25, 7 + + 2.30aL)
X'
h2 16 (6 maxfa + 1.25, +
X
ha 8k' max + 1.25, 7 + ()(922 +.30al)
Unfolding the above inequality, we arrive at
X
k
X h) + + (1
X
"n 1) ) Kh3
2
2
en +
1 + (1
2
X
Kh
1
h
==
2K
1
X
where (i) due to eg Therefore
[Blaz 2213) =
1
2K
E p
(yx
s
cx
min
min
X
19

---

Under review as a conference paper at ICLR 2022
Collecting all the constants and have
(W |lal| du + d + n x'Wa(mo:T))
2K
St
[L - G) max + 1.25, 7 + () (1.7744..777)))
X
X
min
4x'
+ X
(I.
G) maxfa + 1.25, +1} (0.5a + (1.26v 1.140a++22
Sk!
'R'L'
(Win lu1l du d + K (2((T..
+1
mmxxfa + 1.25, 7 + (1.92 2.30aL)
VN
VX
4k
VETT
+1 +1
L max + 12 + +334
VX
VX
C
It is clear that in terms of the dependence on dimension d, we have C = where - I + ma, then X' m + ma. Recall the definition of n' and there exist A', B' such that K² < A' Va + B'. It follows that
O(va). In the regime
7
C a + agat aza² + agat + asa + agat + az <
a² +
+ X
a
a
(a a a+ 7
X
max
for some positive constants 01 02, 43, a4,, as, ac, a7, and independent of a, particular, we have b = O(vd).
C.I PROOF OF THEOREM 5.2
Proof. Denote the k-th iterate of the Strang's splitting method of HFHR by 3g with time step h, the
solution of HFHR dynamics at time hk by æg. Both T and ay start from zo
90
Also denote
Po
40
the solution of HFHR dynamics starting from To time kh by 36 where ão-
(ğo- Po)
Po
go
and E
T). Since is the invariant distribution of HFHR dynamics,
Pp Po
follows that i π.
By Lemma D.5 and Theorem C.I we have
W (uk, p))
inf (εT((A.))))
Elgn 42ll
< inf
E E0 x211
(EII(TT..)
<E Txll
<20 + 2E p- P(xx - 2xx
$20 + +2 P(x 22)
220 - 22 P( 3o) 11
<20 =20 2((k 2X kh W (To x)
Take square root on both sides and apply v < + b. we obtain
W2(me, v) S V2Ch- V2x' X XX'A T).
20

---

Under review as a conference paper at ICLR 2022
C.2
PROOF OF COROLLARY 5.4
Proof. By Theorem 5.2 we have
W2(Me, S /2Ch V2k'e - (2( (2
Given any target accuracy 6 0, if we run the Strang's splitting method of HFHR with h
mhn 2/20 J. then after 1s 1 max 1. 2v2C 2v W)
log
we have
Wa(me, n)) & v2Ch + T 1)) S 2 =.
Recall C O(Vi), when high accuracy is needed, c.g. < 2v2Cho. the itera-
tion complexity reach e-accuracy under 2- Wasserssnein distance is k* 0 0(xx log !)
2v log 2v2x Wa have
فقا)) Recall from Theorem C.I C m (a² 7
C b a² 2 m
+
Denote g(a)
-+ +
simple calculation shows that a argmin g(a) .3-1
m²
OG)
D TECHNICAL/AUXILIARY LEMMAS AND THEIR PROOFS
D.1 DEPENDENCE ERROR OF SDE ON INITIAL VALUES
Lemma D.1. Consider the following two SDE with different initial condition
dr. a(x:)dt + odWi x(0) ID
dy: = a adW y yo
where auu) E TRd is I--iipscitzz and a E Raxx is constant matrix. For 0 < < T. we have the following representation
Kh Un Yo +
with
E E| < 2L2 1|200 20ll22 h2
Proof. Let z = (2% yn) - -(x0 - vo0 = - a(xx) -a(y.)ds. Ito's lemma readily implies that
2E 1" (x, V. a(x, aly.))ds
E E3 ynll² =1|200 yoll
S|100 - voll 21 / El y.1l ds
By Gronwall's inequality, it follows that
< 211100 - yoll. for h 1
Zy Vn.l < |1000
yoll 2hh
4L
21

---

Under review as a conference paper at ICLR 2022
and
h a(xr)) a(y,)ds
( [ [a(xx) a(v))) ds
2
||
E
ds (| [a( - a(v,)) 1|
<h / Blla(1r)) a(y,)ll ds
<1²h E I3 v.ll ds
$21.2
yol 1 1:2
D.2 GROWTH BOUND OF SDE WITH ADDITIVE NOISE
Lemma D.2. Consider the following SDE with constant diffusion
dx; = a(x; + odW..
x00) Zu
where a(x)) Red is L-smooth, ie., la(y) a axx) S fly - ar. a(0) = and E Rdxd
constant matrix independent of time and Tt. Then for 0 < h < aI. we have
El|aa 2oll < 2.57 |10117 + 2h (1121120013)) h.
Proof. We have
=I /' alx;)dt //h odW,
E 1,
xol12
<2R
a Fz dtt
+ +
odW,
1
h
@ 2E
a a(x;) ((;x)dd
+
| dt
27
<2E
||(aa)) |dt)t
+ 22||or1
0 ||a(x.) a(zo)) dt + 1'" ||(z))]| at)
<2E
+ 2h o 113
27
L //"" || =xoll di + h |
<2B
+ 2/1
2 (1/)) || xoll dt dt)
<4E
+ n1?1a((20)1?
+ +2hllo
(66)
|02 4h2 [2(2)))) 4L2h S' Blla,, xoll2
S 2h
22

---

Under review as a conference paper at ICLR 2022
where (i) is due to Ito's isometry, (ti) is due to Cauchy-Schward inequality and ||o|] is the Frobe- nius norm a. By Gronwall's inequality, we obtain
E|a æoll² _ 2 | + a4ha))(()))11) exp /4L
Since a = -[(((2))) a(O) -| when 0 < म., we finally reach at
Bllan - sall² s 2 (1101 2 1 ean = 2.57 (11011 + 22 h.
D.3 LIPSCHITZ CONTINUITY OF THE DRIFT OF HFHR DYNAMICS
Lemma D.3. Assume S is I--lipsciiz, i.e. of HFHR dynamics
Ví(x) vs(D))] & L| yll, then the drift term
p p- aV)())
-YP Vf(4)
Lippishchitz where & V2max vI az max + L), +1))) Let P be defined in Ap-
pendix and
P
q
then
satisfies the following SDE +/b
Φ
r/
p
ddb)
p(Φ. 2b)
aVVlq((..... b))
22xT
0
dW
di +
drb
-xplo. vf(g(o, 2)))
1
dB
and the drift term
p( -Vol)))aV((q(() b))
- - Vf(a(( +1)
is L." Lipschitz, where L," = x L2 and w' is the condition number of P.
Proof. By direct computation and Cauchy-Schwarz inequality, we have
Pi- Pi-aVf(q))
P2 P2- - ((g))) -yp2 Vf(92)
>P: -
a (V/(g)) v)((g)) + + -p))
p2)
v\(()) 7pp -P2)
5 20 v)(1)) < +21lp: - -221" + 2(( - v + -221"
SV (20212 + 2T2) 42ll + (2+2+)))pp; - 1
[q 12
<
ma a2 v
PPi P2
max V 1 max V2 Lt,, v1+
92 Pi P2
L
91 42 Pi
P2
By Ito's lemma, we have
dØ dup
p( v) /Vfll((. 4b))
V2al
0
dd
P
dt +
1pp( V/(g(á, 2b)
0
V27I
dB
23

---

Under review as a conference paper at ICLR 2022
Using the Lipschitz constant obtained for the drift of HFHR, we further have
p 1( - - -V +11)
p(pg, 1/2) aV/(q(0), 4b)) - 422) 2f(g(02.222)
P
> 21) Vf(g(01, +1))
Pi a VV
P2 - [V)(())
<a
-PP vf(g))
-7p2 Vf((a)
< I!
91 q2
P1 P2
Φ1 Φ
<0 max L'
p-1
1/1 2
1
Φ1 22 2/2
< max I'
min
141-1
1 $2
=k'L'
1
12
where max, min and are the largest, smallest singular values and the condition number (w.r.t 2 2-norm) of matrix P.
Remark D.4. The following inequalities associated with will turn out be useful many proofs
221,''222,,,''' V2a, L V2L and L' V2aL.
D.4 CONTRACTION OF (TRANSFORMED) HFHR DYNAMICS
Lemma D.S. Suppose HFHR dynamics
is
L-smooth, m--stronly convex and 722 L. Consider two copies of
qe
qt Pt.
(driven by the same Brownian motion) with initialization
go
pj
Po Po.
respectively, then we have
q, q.] Pt Pc
S -ett
90 9 Po Po.
where P 71 1
0 VI+a7l
and X' = min(" m am,
-- }.
Proof. Consider two copies HFHR that are driven by the same Brownian motion
dq: = (p--aV/(o (V:))V/()))d Bl dp, Vf(g)) 22,dB,
[dg. (p. a VV + 22d dp. (->p, ddt + v2ydB3
Based on Taylor's expansion, the difference of the two copies is expressed as
d 4 d. Tatt -1 Ho
q
d. Pt
A
4: q: Pt Pt
dt P: P.
71
P.
where H, - v2 s - ds.. Denote the eigenvalues of H, by n/. convexity and smoothness assumption on f, we have m n L, < i d.
S Sis d, by strong
24

---

Under review as a conference paper at ICLR 2022
Denote
Φ
P
9: q: P: P:
and consider L2
1
have
111
2
d
Φ
T
[6.
L1
PAP-1
dt
1/.
[4/1]
0
T
Φ
(PAP-1 + +(P-))=A'P"
1/1
2
4
T
(1
(a))
0axd
]
Φ
1.
7
Oaxd
22 H₂.
1
Φ
Φ
B(x)
It is easy to see that
Th }} min( + am, 22 2²44 L a X.
Amin (B(nBB(a)) = min
{min 11,2,
Thi
+ +07,,
7
Therefore we have &. & S 22 2 ((a) - 22 By Gronwall's incquality, we obtain
2
2X's
Se
[//]
(40
and the desired inequality follows by taking square root.
D.5 LOCAL ERROR BETWEEN THE EXACT STRANG'S SPLITTING METHOD AND HFHR DYNAMICS
Lemma D.6. Assume f is L-smooth and 0 argmin f(x), i.e. Vf(0) = 0. If0 h S 4 then compared with the HFHR dynamics, the exact Strang's 's splitting method has local mathematical expectation of deviation of order P1 2 and local mean-squared error of order p2 2, L.e. there
exist constants C₁ C2 >> 0 such that
[[xx(h) [E)(())] Chhp
[1)) =(())1T)
} < CChhz
where (h) =
[(h)
is the solution of the HFHR dynamics with initial value xpp =
90
and
p(h)
Po.
((h)= = [q(h)]
90
is the solution of the implementable Strang's 's splitting with initial value
[p(h)
Po
P1 and P2 More concretely, we have
δ L max(a + 1.25, 7 (1.74 + 22266
O2 I. maxfa 1.25, + 1) (1.92| (1.30/ā 3.22
Proof. The exact Strang's splitting integrator with step size h reads as o wh o o where
à: dg = pdl
lq = a (g)dt + 22ndW
dp ypdt + 2
dp v (()ddt
The ś flow be explicitly solved and the solution is
(t))) = 4o + 1-e-se
Po + 2 I 1-e ((( (( dB(s)
7 p(t) -nit ep"p Po + / So (t-((t))) [1
25

---

Under review as a conference paper at ICLR 2022
The t flow can be written as
a(() Qo Jo (V(((())) 22x So dW(s) (s) (t()) - Po - Jo Vf(g(s))ds
The solution of one-step exact Strang's s splitting integrator with step size h can be written
93 (2(h)) + 1-e 1-e 2
(h) + v J 1-c 7hh-))
dB())
P3 e ph) + v 2 1 1
(2(r) q1 -. av (qg (8)]ds + /2x So ds(ss
(05rsh)
(2(r) Pı -Jo V2 )s
91 = do + 1-e 4
1-c - dB(s)
Po + 27 l0
Pi - Po + 2 1/2 E - dB (9)
Therefore, we have qh) = 93: p(h) = P3 and
1 ())-s)
dB ()) q1 /'" a @VF(g)(s))ds + V2a 1-" ds(s)
)(h)=v21
42(h)
1- 1-E-74 74
h
P1 (s)g
7
P2(h)
1 ()
aV [V[(gg()))ds + V2a / h ds(s)
1 1- 74
/ V[(gg(s))ds
= =/22
dB(s)=
γ
- e 7 h 7
1 ( -)
(B(s) 1
1² E -
+ 90
Po +
22
Po + V2
dB
7
41
Pi
e
7h
1-
Vf(gg(s))ds
-40
Po
7
7
11 1 1- ((h-s)
1 e -- -s)
+
20
dW (8) + +/22
dB(s) + /27
dB
7
0
7
+
27
e e- (BB(s)
h V
((h) -
/27 e-(hh-))
P1
(h)
Po + 24 10 -- -5) BBss) - // V + 2 1" (h BB())
= 11 iho c 7 // V)(gg(s))ds - /
- dB(s) /27 e
26

---

Under review as a conference paper at ICLR 2022
It is clear that (()), p(h) should be compared with the exact solution of HFHR time h, which can be written as
q(h) -00 + 1 2
1 (1- 11---(hh-0) (1-0000000 +
(1 + /27 ("" e-1(h-s) ds(ss)
V[(g(s))ds + V2a /-" dW, /27 c 6 ((h--)
Po
7 dB
ph) =e se-1po 1-¹ h
e
Subtracting q(),, p(h) from @(h), ((h) respectively, we obtain
1 - 1--
! Vf((g(s)) vf(g(s))ds
@(h) ((h)
7
+ /-¹
h
(
1 - e1- 1-
v[((()))ds
7
7
p(h) p(h) -hh) f V[((2())) ((((()))))
(
1/¹ (6-100-0)
v V[(q(s))ds
It should be clear now that we will need to bound the term v5(g2) -V(()) and Vf(g). Since
-
1/2 1 -( dB(s) - / V[(gg()))ds + V2a L dw(s)
(r( Q0 ?
Po + /2y
.0
- ~m
Pu '-"
1---
Vf(g(s))ds v2a (" dw(s)
((r) =0
+a
7
r ((r-s)
+ v 2
(s)),.
0
we then have
1T -
Po a //" vl(g2(>)) - VI((()))ds + { (x(-)
(2(1) - (())
v[((()))ds
?
- e 716 dB(s) /27 1'' 1- dB(s)
27
By Lemma
D.3 and
D.2
when 0 < h
417. we have the following for the solution of HFHR
1
dynamics
E 00 (h) 20 Coh
where Co = 5.14 (a +a (L')² [00 and hence
/ [[ ds
[[](g(0))|| ds +2 f V Vf(q())) ds 5
E
<E
2T.2 [" [|(((s) -=()]] ds 1212r|10011
R2
2T. r
+12E
[lq(s) q(0))|² ds
+2L2O ["" sds
521
SL'r (21100112 + +hCO)
S (.333 20 2 + 5.14(a + 7)dh)
(12)
27

---

Under review as a conference paper at ICLR 2022
Now ||q2 - q|² can be bounded as follow
E
(r) - [(()]]²]
-0 -- 2
/' V[(g((s)) Vi(g(s))
2
((--s)
<5
1pol +
foll2+00 a
+ E
6 (3
2
1 - ((-
+
2yE
dB(s)
+ +22B
dB(s))
(Cauchy- Seuuch--Schwatz Inequality)
7
/'' Bllgz(s) q(>)]| ds /"
(1-
-1--
<5
|
+ a222²r
ds
E B)
ds
7
Ydh
27d
+
+
3 p
12
1/2
/ E
3yd h
<5
+ a 22r
[2)))()) qq())l² ds +
-E
E
+ 4
3
3yd
3
[2.33 5.14(a + z)(h) r a2 / EE [l q(s)]l ds
<5
[0ll
+
113
+
L
4
3
37d
h²
1 # +5.14(a + m)hh } 55500 f El(g2(s) - a(>) ds
S5h²
Fo
+
T?
2.33 2o
4
3
By Gronwall's inequality and 0 < h S
L. we have
|000| + 3yd h
1i2 1 (2.33) æa 2 5..44 + t)dh)
E (1))
<(r))/]]
<5h²
exp(5a2 1.?2
| + 3yd h
h 2 (2.33)
<5h²
æull +5.14(a + t)dh)
e
0..28| <5.85h2
2 (0.06a + 0.817)
2h²
(1.44|
(0.36a + 4.74
(13)
28

---

Under review as a conference paper at ICLR 2022
With bounds in Equation (12) and (3), ready show P1 and p2. For P1, ie. the order of the mathematical expectation of deviation, we have
((h) (phh)
(h)
E
p(h)
(h)
+ [p(h)
[p(h) - p[h)]
)110 [V((gg(()) - [ ds + $1//" ( 1 (1--1(h
(1- 1-
1
a+
[f(((())]] ds
1 78 / E (s)) v((g(s))] do d +|/ ( - E [V((g(s))] ds
a 1 +
L (l) q(5)) ds
+ 1."
(---)()-)) ((- 1
+ + - -
|=
7
γ
2)
E 1/"
SI
a +++ 2
(2() q(()]| ds
2 1/2
: }
s( 1 e-
+ (f% (h-
78 - ds)
(1. [VV[[((
| ds
ds
7
1
h
1/" (((42(s)-((())))) (l)a)())((<s]()))) ds 2
1 1+7,3
|| d.
ST.
+ 1
h
E.
2/3
<2 a + 1+ h
122 [1.41|1011 (0.36a + 4.747) 2
1 1+ n22 (2.3 æoll *+ 5.14(a + 7)dhh)
2/3
<L (a <L 1.25) /22 (1.29|/zool + 00.36a + 4.747 hd)
+ (1 + 7)h2 L (0.45 |2ol|
V0.3a 0.43yv
STh² max(a + 1.25, + 1} (1.74||zol| + (1.26 + /784 hd
The above derivation proves Pl = with
Ō = L max + 1.25,7+1] (1.74]|zoll + 22
29

---

Under review as a conference paper at ICLR 2022
We now proceed with P2, i.. mean-square error
[(h]
(h) [p(h)
E
p(1)]
a + h 2
//-0'00 Vf(()(()) - -
+ +2E 1/"
(1- (1-e-
<2
v vf(((s))dds
2
?
Vf(q2(5)) vf(g(s))ds + /// (e-T0h -
+ 2E
V Vf(g(s))ds
h y +1 L2 (2((8) ((s) +2 ph
7(h-a)
1 1
ds 1 (|] ds
22 aa +
/
e 1 ds ( E [ ds
+
(h-5))
11 2
72 3 L E EV
(a
+1 h Elq2(s) q(s) ds +
6
<2
(a +
y +1 1,2 { 1.64||20o) + (0.36a + 4.747) htt +
h1 1/2
1 I
+5146a 2.33 2 - +7)hd)
6
<L² ma + 1.25) + > o.07 + (1.68a + 10.347)ld) 11
The above derivation implies p2 = 2 with
E2 Lmxxfa + 1.25, 7) (1.92| + (1.30/a++22 3.22.7222222 ()
D.6 LOCAL ERROR BETWEEN ALGORITHM 1 AND THE EXACT STRANG'S SPLITTING METHOD
Lemma D.7. Assume f is L-smooth, 0 argmin f(x), i.e. Vf(0) = 0 and the operator
ar. then compared with the exact Strang's splitting method of HFHR dynamics, the implementable Strang's splitting method has local mathematical expectation of deviation of order = 2 and local mean- mean-squared error of
Af grows most linearly, i.e. VDf(a)) GV1++lg9 If0
order P2 1.5, i.e. there exist constants C1, O2 > 0 such that
[E - [b)]l s O.hm
E |1((h)) [|0(()-=(<)]f() =(M))1)) 3 < C2hPa
[((h)
where drh) =
is the solution of the exact Strang's splitting method for HFHR with initial
p(h)]
[q(h]]
value
4o
and at(h)
is the one-step result of Algorithm I with initial value pphh]
Po
æo
90 Po
P1 2 and p2 1.5. More concretely, we have
C = a(a aaa G) 0..5 + 0.71| + (..14V
and
C2 = L(a + 0.73) (2.00 +2270
30

---

Under review as a conference paper at ICLR 2022
Proof. The solution of one-step exact Strang's splitting integrator with step size h can be written as
93 = 92(h) + 1-e ~
P2 + +2 th 1-e vCh dB(s)
P3 = e 2 p + v e dB(s)
92(r) 91 - -J' [V + 22 ds((s)
(05rsM)
(2(r) Pı -S Vf(92 )ds
91 do + 1-e75 7
Po + 27 1- -- dB(s)
e - p0 + 22 Sot e - (s())
and the solution of one-step implementable Strang's splitting integrator with step size h can be written as
93 = 92(h) + l- -4 (h) + /2 10 1-c 7 "
dB
( -s)
P3 e e- p - /2 fo² - dB +
92 91 -S0 [VV[(qg)ds + 2a So dW(ss)
(05rsh)
(r( = =P1- -Jo V )ds
do + 1-c
Po +
27 Jo
1-c - dB(s)
Pi - Po 22
e (s)))
Note that in the implementable Strang's ' splitting method, & flow can be explicitly integrated and hence Pi are the same as that in the exact Strang's splitting method.
First, we will bound the deviation of mathematical expectation and mean squared error of 42(h) (h) and p. P2 (h).. We have
(h) (2(h) -a ( Vf(g))
(14)
(2(h) P2(h) - - vf(g2(s)) V (f
Square both sides of the first equation in 14 and take expectation, we obtain
(h) - (2))1
=
V)((g(+)) VS(g,)
Z ( [ vf(g
L2E
[19x(s) - -qilds
<aa L² h
(2 9, 11
ds
Note that qa is the solution of a rescaled overdamped Langevin dynamics whose drift vector field is a--ipschiit, by conditional expectation version of Lemma D.2 for / aaL. we have E (1)) - qull 11 S Coh with C. 5.14 { ad + M(..))
|| and follows that
E (h) 92(1)]1 a²12Cgh2
E (2(h)) - p2(h1)] L2Coh3.
Now consider Pt. i.c., the deviation of mathematical expectation. By Ito's lemma, we have (h) (h)
=-
V V[(gg(s))- -
=- a 1/") [/"" -oV> +a VAS(qa()))dr +p/ ds (15)
31

---

Under review as a conference paper at ICLR 2022
where p is stochastic integral term. Take expectation and norm for Equation 15 92 (h) - (b))
we have
( ds
v2 ((
h
l V/(a ([IVA/(((())))) ds
Za² 2
E
||
L G(L G(1+|1a>())]) ds
5a²
E
! []|ga(T)|| dr +a²G h²
= 6)
L. Bllg(()) 9 E dr + a²G 1t
Sa G
2
L. V (1)) 9. 1| +B dr + a
h2 2
(O((L+))
12
E 91 + a²c h?
<a + GJV Cal + a (L
+
2
2
VChh + 91
G
sa²
(L + G) 2
h
2
< (T. G) Coh + +B| + +1] +1} 122
[22 (h)) S { G Elq, |||++11 +1] +1]12. 12.
Similarly, we have
For P2, i.e., mean-square error, we have
[}([g(s)) - V/laall 2
E (1)
o
<I²FR { Ids /"11 -V/ ds
<a² 1 /"' Ell q, 1212ds
a² L.200 113
2
Similarly we obtain E P2(h) - p>(h))l?
L'Ca 1², Recall
2
9a 43 (h) (h) + 1-c-74 (p2(1) P2))))
Ps P3 = e - 2 (p)(h) (2))
and it follows that when 0 < < 37 <1
VCah +Ellqil ll + +1,²
43 4a
<a(a + + () [ G)
(16)
P3 Ps
2
E
43 q3
<I'C. a² +
12
h3
(17)
Pa Ps
32

---

Under review as a conference paper at ICLR 2022
Finally need to bound Bl
11 by Bl|æo|l to this end, we have
1
1-
1 1
11 2
-E
40
Po + v27
dB(s)
h
BBl1o 12 1 +1+ h
4
1--
1
-s)
<1 +
Ellpol + +2yd 2nd
ds
7
7d
<(1 )B|12001 112
12 133
(18)
42
yp
-(1 ||z0|||
+ 113
(19)
12
Collecting all pieces together, including not difficult obtain the following
16
19).
the definition of C. and 0 < h <
T it
E
Q3
Q3
<O1h²
P3
P3
93 43
CChh aa
E
Ps
Ps
with
C = + (1125 G) [0.5 0.71|/z0l + 4..14/ Vhd
and
C2= = L 0.73) (2.30VhaLL(z)) +
+
0.12
D.7 LOCAL ERROR BETWEEN ALGORITHM 1 AND HFHR DYNAMICS
Lemma D.S. Assume f L-smooth, 0 argmin f(z), i.e. Vf(0) and the operator
V grows at most linearly, i.e. [af(g))) GV +|lall +Vllll If0 < h S 37/ then compared with the HFHR dynamics, the implementable Strang's splitting method has local weak error of order
P1 = 2 and local mean-sque seaan-squard ervor of order P2 1.5, iee there exist constants C1. C2 >> 0 such
that
[[Ex(h) - (x) s c
E [[(((h) ((e(h)==(())(()) < C2hP2
where wh) = [qhh] [(h)]
[9o]
[q(h]
is the solution of HFHR with initial value =
and æ(h) =
Po
[ph)]
[q]]
is the solution of the implementable Strang's splitting with initial value xo
Pi 2 and
Po
P2 1.5. More concretely, we have
C
7+1} (0.5a + (1.74 0.71a) (1.26 1.14a/a + 2.32 Vhd
and
C₂ = L max + 1.25, + 1} (1.92 2.30a +2.00
33

---

Under review as a conference paper at ICLR 2022
[q(h] [P(h]
Proof. Denote by ((h)
the solution of the exact Strang's splitting method with initial
40 Po
value vo
By triangle inequality and Minkowski's inequality, we have
(h)
Ex(h)) (x))) [=)( (a)()) - (x))))
[B))=(hh) =0()))]))) < (Blla0(1) =(())l²) (Bll )
By Lemma
D.6 and
D.7
we have
Ex(hh)) Ex(h)] s Cih?,
[[ Ex(A)] & h2 (B|@(h)) (aochh->((h((()) =({))(²)* 2 C2h
(Bllz(hh) =(m))i)) - s C2h
and hence
(h)(h) (t)()) || s(Ōi + +i)A² [B] =(h))l²)* <(Ō² + C2)
with
C SC
+ G) max + 1.25, 1}
0.5a + O.70))l)) +.1 Vhd
Ö2 + C2 <C₂ e Lmxx(a + 1.25,
(1.92 + 2.30aL) +2.00
E a DOES CREATE ACCELERATION EVEN AFTER DISCRETIZATION: AN ANALYTICAL DEMONSTRATION
a + oo while remains fixed, then dg = -VV((g + 22dW is the dominant part of the dynamics, and in this case the role of could be intuitively understood simply rescale the time of gradient flow, which does not create any algorithmic advantage, as the timestep of discretization has to scale like 1/a in this case. However, finite a longer corresponds solely a time-scaling, but closely couples with the dynamics and creates accelcration. This is true even after the continuous dynamics discretized by an algorithm.
We will analytically illustrate this point by considering quadratic f. In this case, the diffusion process remains Gaussian, and suffices to quantify the convergence of its mean and covariance. In fact, it can be shown that both have the same speed of convergence, and therefore for simplicity we will only consider the mean process. Two demonstrations (with different focuses) will be provided.
Demonstration (1D, y given; infinite acceleration).
Consider f(x) = x² /2, y fixed. The mean
process is
p aq
i = -q 7p
Consider, for simplicity, Mnler-Maruyama. discretization of the HFHR dynamics, which coressponds to Forward Euler discretization of the mean process (other numerical methods can be analyzed analogously):
9k-1
A
g Pk
A
1 ah h
Pk+1
1 1-7h
We will show that, unless = 2, an appropriately chosen a will converge infinitely faster than the case with 0, if both cases use the optimal h.
34

---

Under review as a conference paper at ICLR 2022
To do let compute A's eigenvalues, which are
1 2 (a ( -1- (a-
2
Consider the case where la - > - 2, then the eigenvalues are a pair of complex conjugates. Their modulus determines the speed of convergence, and it can be computed to be
(2 (a (a+)))+)))--- 4 7)))- = V (a + +)4 ar)
Minimizing the quadratic function gives the optimal h that ensures the fastest speed of convergence, and the optimal h is
+ h = 2(1 + ar)
and the optimal spectral radius is
(a + 7)² 4(1 + +ay)
When one uses low-resolution ODE, in which a 0, the optimal rate is 1 72/4 (note it is not surprising that the critically damped case, i.e., 2, will give the fastest convergence).
If y *
the additional introduction of a can accelerate the convergence by reducing the spectral radius. For instance, if a y + 2, upon choosing the optimal h I+y" the optimal spectral radius
(note in this case A actually has Jordan canonical form of 0 and thus the discretization 0
converges in 2 steps instead of 1, irrespective of the initial condition).
Demonstration 2 (multi-dim, 7, and h all to be chosen; acceleration quantified in terms of condition number). Consider quadratic f with positive definite Hessian, whose eigenvalues are 1 À = /2 + E 9///2. Similar to Demonstration 1, the forward Euler discretization of the mean
< An = for some 0 < < 1. Assume without loss of generality that
process is
91,k-1
Q1,k
P1,4+1
A
Pi,k
A
1 ah h
Aa
I ae 1/2 - -1-
h
92,
0
A2
92.k
1 1-79h -h
1 1-ch
P2,641
P2,1
(20)
We will (i) find h and y that lead fastest convergence of the ULD discretization, i.e. the above iteration with a = O, and then (ii) constructively show the existence of h, 7 and a that lead to faster convergence than the optimal one (i) note these may even be the optimal choices for HFHR, but they already lead to significant acceleration. More specifically,
(i) In ULD setup, a 0. It be computed that the eigenvalues of A₁ and A₂ respectively
2 2 - by = h/ 4+ 4+72
and 2 (2-Ma) h- + ht -4 + +
We now seek 7 0, h > 0 to minimize the maximum of their norms for obtaining the optimal convergence rate. This done cases.
Case (il) When < 2, both A and Ag eigenvalues are complex conjugate pairs. To minimize the maximum of their norms, let's first see if their norms could be made equal.
A₁ eigenvalue's norm squared x4 is
(2 hr7 - 1 4(h- + +4- 7-72
(21)
Ag eigenvalue's norm squared x4 is
(2 17)² 1/2 (-46 -6
7²) 4e - er/2)² +4 632
(22)
35

---

Under review as a conference paper at ICLR 2022
It can be seen that for (1 is always strictly smaller than 22 for any h 0. Therefore, the max of
the two is minimized when = 2,/2, and the corresponding max value is - eey. y that minimizes this max value is y 2. Corresponding rate of convergence is
-
Case (12) When y 2 2c 2/-1 both A and A2 eigenvalues are real. Since e < 1, we can order themx2 as
2 2- -4 2 2-1 -4- < 2- 2-1
< 22 -4+ <.
To minimize the max of their norms, consider cases in which the smallest of four is negative, in which case at optimum one should have
-(2 - bry 4 2 by hy
This gives 2/y (which does verify the assumption that the smallest of four negative). Corre- sponding max of their norms is thus /I - 4/72. y that minimizes this max value is y = 2e -1/2 which gives rate of convergence of
v1
Case (i3) When S Y 2c 1/2 A1 eigenvalues are real and A₂ eigenvalues are complex conju- gates. Again, the max of their norms is minimized if the norms be made all equal.
Note A eigenvalues cannot be of the same sign, because otherwise 2 hy h 1VV-441 4
hu + / / -44 + 72, which means either h or 7 2, but if 2 then hy + hi - 72 e7 leads to h being equal to n of A2 eigenvalue, which is 4e 4e- (h + again.
Therefore, the equality of norms of A1, A2 eigenvalues means
-(2 - buy li 22) 2 by + /v -4+ / 4e 46 - + 672
The first equality gives by = 2, which, together with the second equality, gives h = + = 2r
Ite'
Selecting the positive value of optimal h, we also obtain optimal γ 2 e)e-1/2 which < 2e -1/2 and thus satisfying our assumption (2 I S 21/2). The corresponding rate of convergence is thus
(2 mr + hv -4+
1-c Ite
vI- the ULD Euler-Maruyama discretization converges the Summary of (i) Since life
fastest when
2c
? 2 e)ee--
Itte
and the corresponding discount factor of convergence (i.c. base of exponential convergence) is c
where = 1/x with being Hessian's 's condition number.
(23)
1 I+te
(ii) Now consider the HFHR setup. Let's first state a result: when
V4cee + 8 4c + e² 2e 1
> 0,
(24)
206 + 20x
-v 80283 4 4 2c + 36
0, = ce
(25)
2ce2 + 2cE
for any c > 0 independent the iteration 20p converges with discount factor
v2(1+e) (11 ()
(1-c)(1-e_)) (1- 4c2e 8 + (4c2 1) 2 2e +
(26)
36

---

Under review as a conference paper at ICLR 2022
While the exact expression is lengthy, can proved that the HFHR non-optimal discount factor (26)
strictly smaller than the ULD optimal discount factor (23) for not only small but also large e's.
For some quantitative intuition, discount factors respectively have the following Taylor expansions
in E:
c²
è 2+0 20)(2))
HFHR non-optimal:
2e
+
2
(27)
2 + +
ULD optimal:
1 c 2
(28)
acceleration enabled by HFHR
0.8
f factor
0.6
0.4
0.2
HF HR non-optimal hyper-parameters
ULD optimal pareeessssmeters
0
0.2
0.4
0.6
0.8
Figure 4: Acceleration of HFHR algorithm over ULD algorithm (despite of an additional constraint may place on h) for multi-dimensiona. quadratic objectives. 1/c is the condition number.
The exact expressions of discount factors are also plotted in Fig (c 1 was arbitrarily chosen) and one can see acceleration for any (not necessarily small)
(ii details) How were values 25 chosen? Following the idea detailed in (i), we consider a case where A₁ eigenvalues are both real, A2 eigenvalues are complex conjugates, and their norms are equal. Note there are 3 more cases, namely real/real, complex/real, and complex/complex, but we do not optimize over all cases for simplicity the real/complex case is enough for outperforming the optimal ULD.
This case leads to least the following equations
trA1
(29)
det A1 + det A2 0
One can solve this system of equations to obtain a and functions of h. Following the idea of choosing h small enough to resolve the stiffness the ODE
q2
ae- = - -1 122 7P2
P2
pick h Then 29 gives
2 + 8 2 + 2 2e + + +
2ce² + 2cc
-4 + 8 + 4c + 2 2e + + +
a
200 + 2ct
or
V4cee 8 + 4 + t 2 + + 33
=
20e + 2cc
V 44cee + 20 + 4622 + 2 2e + 1 3E +1
2002 + 2cc
37

---

Under review as a conference paper at ICLR 2022
The former is choice (25) because it can be checked that the latter leads to det A1 0 which violates the assumption of a pair of plus and minus real eigenvalues.
It is possible to find optimal a, y, h for HFHR for the Gaussian cases. One has to minimize det A
under the constraint det A2 > 0 in addition to (29). And then do similar calculations for the other cases, and then finally the best among the 4 cases. Doing however does not give enough insights to determine optimal hyperparameters for sampling general distributions.
F RANDOMIZED MIDPOINT DISCRETIZATION OF HFHR
F.1 THE ALGORITHM
HFHR is based on continuous dynamics that adds HFHR corrections to the Underdamped Langevin Dynamics (ULD). It can be turned into a sampling algorithm via either low-order time discretization (e.g., HFHR Algorithm D more accurate one. To complement the main text, this section demonstrates the latter, based on a powerful recent progress in discretizing ULD, known Randomized Midpoint Algorithm (RMA) (Shen & Lee 2019), and shows that the acceleration created by the HFHR correction terms persists.
More specifically, RMA is a high-order discretization scheme for ULD that achieved a better O(dš) dimension dependence of mixing time than first-order discretization of ULD, e.g., Ist-order KLMC (Dalalyan & Biou-Durang 2020). Although RMA is originally designed specifically for ULD only, is general idea and already adapted overdamped Langevin (He et al. 2020). Here we show can be easily adapted to HFHR as well, as illustrated by the following Algorithm Red highlights algorithmic changes we made to account for the HFHR corrections of ULD.
Algorithm 2 Randomized Midpoint Algorithm from [hen & Lee (2019), adapted for HFHR
1: Input: potential function f and its gradient V., damping coefficients a and 7, step size h, initial condition (qn- Pn) (f. V ar 7, 4o 4o.
2: procedure
3: = and initialize
4o
Po
4: while not converged do
5: Generate an independent uniform random variable 02 ~ U(,
6:
Generate Gaussian random vectors [W] &+11 k+
W
W k+1
E 3d as in Shen & Lee
2019
Appendix A)
7:
Generate Gaussian random vectors B
k411 B k+1 Rd as described by B1
8:
au+# Q: + (1 e
--
Oah (1 E -Bab)) Vi(ga) + W k+1
a0.h (qt)) + 12 k+1
9:
9s+1 = 9s + 1(1 <>>>)pp t/1l1 e-1(h- - W
ah V f + 22a 1 + +Bb+1)
Pk+1 = e 7h he-1hhh 7(h-0ath)
v +2
11:
k+
12:
end while end procedure
13:
The red parts basically correspond to two Fuler-Maruyama time-steppings of an auxiliary dynamics that contains only the HFHR correction terms
dq
VV + 20d.B
(30)
first over f&h timestep, and then over an h timestep. These steps originate from operator splitting treatment of the full HFHR dynamics (eq(o), which is split into ULD and (30). Therefore,
it is natural to see that
h(k+0%)
h(k+1)
BL+
dB..
B
/ dB,,
Akk
38

---

Under review as a conference paper at ICLR 2022
and therefore + and Bi++1 are, when conditioned on 8k, centered Gaussian vectors independent from each other and the W's, ss each being d-dimensional with i.d. entries, and they can be generated via
B = Vohhk++1
B = v 0
(31)
where Sk+1 and Ex+1 are i..d. standard d-dimensional Gaussian vectors.
Remark F.1. In the original RMA (Shen & Lee 2019 Algorithm 1), the uniform random variable for the midpoint's proportional location was denoted by a. However, since we have already used this letter for the HFHR correction coefficient, we use instead
denote this uniform random variable.
Remark F22 From the red text, it is easy to see that if O, Algorithm 2 degenerates to RMA for ULD. Nevertheless, Algorithm i again just one RMA discretization of HIFR but not the only one.
F.2 NUMERICAL RESULTS: HFHR AGAIN ACCELERATES
To numerically compare the RMA discretization of HFHR dynamics and ULD dynamics (note we don't compare Ist-order HFHR Algorithm [ with RMA-ULD as we'd like to com- pare apple with apple), we conduct an experiment very similar that in Sec with the same nonlinear potential function. We run both RMA for ULD and RMA for HFHR with dimension
1 (chosen to be near the stabil- ity limit of RMA-ULD), a family of 7 E {0.1, 0.2, 1, 10, 20, 50, 100) and {0, 0.001, 0.002, 0.005, 0.01, 0.02, 0.05, 0.1, 0.2, 0.5, 0.55, 0.6, 0.65, 0.7, 0.75, 0.8, 0.85, 0.9, 0.95, 1, 2, 10, 20, 50, 100}. For each algorithm and each set of parameter values, we run 1,000 indepen- dent realizations to compute statistics and estimate the mean time of reaching e 0.1 neighborhood of the target distribution. Then, for each (including a 0, which is the original RMA), optimize over "y choices to get the best results. To further reduce variance, we also repeat the experiment with 100 different random seeds.
10, initial value (100 x 1a, 0a). h
Too large a values with which Algorithm 2 fails to reach e-neighborhood are not plotted and the final results shown in Figure 5 It clearly suggests that with appropriated chosen (a: 0.5 in our case), RMA discretized HFHR dynamics requires fewer iterations than RMA discretized UL.D, which suggests better iteration complexity.
RMA
# of iters to reach E-closeness
HFHR RMA
16
12
10
α
Figure 5: Improvement of RMA for HFHR (Algorithm 2 over the original RMA (for ULD) in iteration complexity. (vertical bar 1 std.)
39


Reviews:
[{"review_id": "uAD32DXgZ6M", "reviewer": "Reviewer_4H3p", "summary": "", "questions": "", "limitations": "", "rating": 3, "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": "This paper has solid theoretical analysis. \nHowever the author missed an important comparison with randomized midpoint method (Shen & Lee, 2019). \nI think the author proved an iteration complexity that is worse to the method in that paper. This invalidates author's claim on constant factor speed-up over underdamped Langevin dynamics.\n\nSpecific Comments:\n1) In Assumption A1, $m \\lVert y-x \\rVert \\leq \\lVert \\nabla f(y)-\\nabla f(x) \\rVert \\leq L \\lVert y-x \\rVert$ seems insufficient to establish m-stronly-convex and L-smooth.\n\n2) Does this paper consider \"log-concave\" case? In the introduction, author said \" will consider the setup of log-concave / log-strongly- concave target distributions\". However, it seems Assumption A1 require $m > 0$. \n\n3) On page 9, author wrote \"(strongly) log-concave assumptions required in Theorem 5.2\", why there is a bracket? Does it imply that Theorem 5.2 could be applied in log-concave case?\n\n4) In Remark 5.5, at the last line, the author claimed \"steps needed by ULD can be halved by HFHR\". Which algorithm is the author referring to for \"ULD\". Is it algorithm in (Cheng et al. 2018) or HFHR with $\\alpha=0$? It seems these two algorithms are actually different.\n\n5) In section 6.1, author proposed to use error of mean as a surrogate of 2-Wasserstein distance. Although error of mean is a lower bound, how tight it is? In figure 2, the mean value at $\\alpha=0.5$ in the Figure 2 seems close to 1, which means HFHR achieves $\\varepsilon\\leq0.1$ accuracy in just one iteration. Clearly one iteration is not enough for burn-in. Does it means error of mean is a bad indicator of mixing?\n\n6) In figure 1, could the author also add lines for ULD and randomized midpoint method?", "summary_of_the_paper": "This paper introduces Hessian-Free High-Resolution SDE inspired by Accelerated Gradient (NAG). The author shows that continuous solution achieves an acceleration over the underdamped Langevin. A discrete algorithm also has speed-up with a constant factor.", "summary_of_the_review": "I thinks this paper could be improved substantially, but currently would like to recommend reject. ", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "2: Several of the paper’s claims are incorrect or not well-supported.", "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": "kPmp4QP_xNr", "reviewer": "Reviewer_97Mf", "summary": "", "questions": "", "limitations": "", "rating": 1, "confidence": 1, "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": "I do not think the problem addressed by this paper does not make sense.  The gradient-based algorithm added by the noise in the convex case does not play the essentially different role with that without the noise.  ", "summary_of_the_paper": " This paper discuss the accelerated-gradient-based MCMC method and propose some rigorous proof.", "summary_of_the_review": "Please the authors refer to \n\nBin Shi, Weijie Su and Michael I. Jordan,  On Learning Rates and Schrödinger Operators, https://arxiv.org/abs/2004.06977\n\nBin Shi, On the Hyperparameters in Stochastic Gradient Descent with Momentum, https://arxiv.org/abs/2108.03947\n\n", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "2: Several of the paper’s claims are incorrect or not well-supported.", "technical_novelty_and_significance": "1: The contributions are neither significant nor novel.", "empirical_novelty_and_significance": "1: The contributions are neither significant nor novel."}, "extra_text": {}}, {"review_id": "f2I-Usp47nM", "reviewer": "Reviewer_awRv", "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": "The paper is very technical and not easy to follow. The main contribution of the paper is to provide an accelerated first order diffusion process in continuous time, yielding a novel sampling method after discretization. I have a few questions/concerns to be clarified regarding the claims in the paper:\n\n1 Acceleration in continuous time \nComparing to the rate obtained in Dalalyan 2020, the proposed algorithm improves by a factor of $\\kappa$. However, I am confused as no improvement is achieved for the function $f(x,y) = mx^2+Ly^2$, this seems contradictory to the claim, do I misunderstand anything? Is it related to the time scaling? \n\n2 Discretization \nWhen discretizing, the paper applies a second-order symmetric composition. What is the motivation behind this step? What would be the benefit compared to applying forward Euler discretization? \n\n3 Acceleration in discrete time \nUnless we show that the constant term is improved by a non-trivial factor (L, or $\\kappa$ for example), I wouldn't call it an acceleration as those are only upper bounds in the analysis. \n\n4 ULD = HFHR(0,$\\gamma$)\nI understand that the continuous ODE are equivalent by taking $\\alpha = 0$, is the discretized algorithm still follows the same equivalence? (for instance, different discretization method may be applied)\n\n5 Comments on Experiments \na. In figure 2, it seems surprising to me that only 2 iterations is enough to reach $\\epsilon$ closeness, which may suggest the problem is too easy\nb. In Figure 3, why does some figure only has $\\alpha=0.1$ and the others has $\\alpha = 0.1, 0.5, 1$\n", "summary_of_the_paper": "The paper proposes an accelerated MCMC method for sampling, motivated by Nesterov's Accelerated Gradient (NAG) method. Starting from the high resolution ODE of NAG obtained in Shi et al, the paper applies a two stage mechanism to remove the Hessian-dependency. The obtained first order ODE system serves as the backbone of the proposed diffusion process, Hessian-Free High-Resolution (HFHR) dynamics. Discretization the HFHR dynamics leads to the proposed sampling method. Theoretical convergence are provided for both the continuous and the discretized variant, showing acceleration to existing underdamped Langevin method (ULD). ", "summary_of_the_review": "Overall, I believe the theoretical result of the paper has merit but the presentation need to be improved ", "strength_and_weaknesses": "", "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": "aRzhDwGyR8J", "reviewer": "Reviewer_r34U", "summary": "", "questions": "", "limitations": "", "rating": 6, "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": "This paper builds a NAG-based MCMC sampler, which is named as Hessian-Free High-Resolution (HFHR($\\alpha, \\gamma$)) and generalizes the underdamped Langevin (also can be written as HFHR($0, \\gamma$)) with an additional gradient drift and Brownian motion for the update of the position variable. \n\n**Pros**: The most prominent feature of the extended algorithm yields much faster exponential convergence with **rate $O(L)$**  in continuous time than the standard underdamped Langevin algorithm  **(Cheng'18, rate $O(\\frac{m}{L})$, Dalalyan'20, rate $O(\\frac{m}{\\sqrt{L}})$)**. Although the acceleration becomes less significant in numerical algorithms due to an increase of discretization error, some speed-up by a constant factor can be still achieved.\n\n**Cons**:  the writing is not very clear, which affects the readability experience and made me hard to check the details. For example 1) the derivation of Formula (6) is slightly ad-hoc and more related work or motivations are suggested; 2) at the end of section B, some link on the definitions of notations may be suggested; 3) page 24, $H_t$ is derived based on Taylor expansion instead of saying nothing is better for presenting to the readers.\n\nRemark: Theorem 5.1 is checked carefully; Theorem 5.2 is not; the logic seems to be reasonable.\n\nMinor issues: \n\n1. Experiments conducted on section 6.2 is non-convex and doesn't match the theory.\n\n2. Missing part in related works: Replica-exchange (a.k.a parallel tempering): \n\n[1]. Accelerating Nonconvex Learning via Replica Exchange Langevin Diffusion. ICLR'19\n\n[2]. Non-convex Learning via Replica Exchange Stochastic Gradient MCMC. ICML'20\n\n\nCheng'18: Underdamped Langevin MCMC: A non-asymptotic analysis. COLT'18.\n\nDalalyan'20: A. S. Dalalyan and L. Riou-Durand. On sampling from a log-concave density using kinetic Langevin diffusions. Bernoulli\n\n", "summary_of_the_paper": "The author proposed an accelerated-gradient-based MCMC method based on Nesterov's accelerated gradient (NAG). In continuous time, the algorithm is able to achieve a tremendous acceleration over the underdamped Langevin algorithm. Numerical schemes can propose a speed-up with a constant factor.\n", "summary_of_the_review": "The NAG-based MCMC sampler extends the underdamped Langevin with a drift term and Brownian motion and is proved faster than the alternative. I would recommend this paper to be accepted.", "strength_and_weaknesses": "", "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": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work."}, "extra_text": {}}, {"review_id": "Pjwy4-erZSK", "reviewer": "Reviewer_xkFf", "summary": "", "questions": "", "limitations": "", "rating": 3, "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": "Strengths:\nThe HFHR process seems novel and has interesting connections with NAG. The paper is mostly clear and easy to follow.\n\nWeakness:\n1. HFHR has noise injected into position directly\n\nThis property is actually harmful for discretization because directly coupling Brownian motion into $q_t$ makes the trajectory of $q_t$ extremely non smooth and make it hard to estimate $\\nabla f(q_t)$.\n\nThe ULD avoid directly inject noise into position to achieve acceleration compared to the overdamped LD.\nHFHR, unlike from ULD, has noise term $dW_t$ in position. However it still achieve same asymptotic convergence as certain discretized ULD at the price of an extra condition on higher order derivative.\n\n2. Rely on additional third-order growth condition\n\nAlthough this condition is weaker than Hessian Lipschitz and has been used in prior art to accelerate the convergence of LD [Li et al., 2021], it still makes the convergence guarantee weaker than many discretized ULD, which have same or better asymptotic convergence speed and don't require third-order growth condition.\n\n3. The condition number $\\kappa$ dependence and third-order growth constant $G$ dependence in iteration complexity is unclear\n\nThe main paper doesn't specify the asymptotic dependence of $\\kappa$ and $G$ but hide it in the constant $C$. Including that information in the main paper would be appreciated.\n\nAccording the appendix page 20, where C is defined, I guess that the $\\kappa$ dependence is $O(\\kappa^{3/2})$ and $G$ dependence is $O(G)$. Could the author help to verify the above statements?\n\n4. Iteration complexity is worse than certain discretization of ULD\n\nAlthough the authors cited [Shen & Lee, 2019], they didn't discuss methods in this paper at all.\n\nThe randomized midpoint method is a special discretization of ULD, which achieves $\\widetilde{O}(d^{1/3}/\\varepsilon^{2/3})$ iteration complexity and doesn't require third-order growth condition.\n\n5. The method is not well motivated\n\nAccording to the paper, HFHR is motivated by one question \"how to appropriately inject noise to NAG algorithm in discrete time\". We should note that, although there are underlying connections between optimization and sampling, injecting noise into a good optimization method doesn't necessarily yield a good sampling algorithm.\nUnfortunately, this paper seems like this is the case. More specifically, discrete HFHR with $\\alpha \\neq 0$ have worse iteration complexity compared to certain discretization of ULD.\n\n6. Points in Figure 2 is biased estimation\n\nIf I understand correctly, the y-axis value in Figure 2 is smallest number $k$ such that $\\lVert \\mathbb{E} {\\mu_k} q-\\mathbb{E} \\mu q \\rVert \\leq \\varepsilon$.\nIn order to estimate $\\mathbb{E} {\\mu_k} q$, $\\frac{1}{n} \\sum_i^{n} q_k^{(i)}$ is calculated for $n=1000$ samples (according to line 3 in generate_config.py) which are generated and run in parallel for $k$ iterations.\n\n$\\frac{1}{n} \\sum_i^{n} q_k^{(i)}$ is indeed unbiased estimation of $\\mathbb{E} {\\mu_k} q$. However \n$\\min k$ $s.t. \\lVert \\frac{1}{n} \\sum_i^{n} q_k^{(i)}-\\mathbb{E} \\mu q \\rVert \\leq \\varepsilon $\n is a biased estimation of \n$\\min k$ $s.t. \\lVert \\mathbb{E} {\\mu_k} q-\\mathbb{E} \\mu q \\rVert \\leq \\varepsilon$.\n\nTypo:\n\ndt is missing in eq 6.\nBracket error on 3rd-to-last line on page 8. ", "summary_of_the_paper": "This paper introduces a stochastic process called HFHR. The SDEs of HFHR is derived by rewriting Nesterov’s Accelerated Gradient (NAG) into phase-space representation, formulating it as ODEs, and then injecting noise into both position and momentum variables. HFHR can be used for sampling because it has stationary distribution as the target distribution. A discretization of HFHR is given by operator splitting and Euler integration. A mixing time bound of $\\widetilde{O}(\\sqrt{d}/\\varepsilon)$ is obtained for sampling log-strong-concave-and-smooth target distribution with extra third-order growth condition.", "summary_of_the_review": "The method is novel in the sense that it is derived based on the idea \"inject noise into Nesterov’s Accelerated Gradient\".\nHowever the derived theory is not significant as it rely on extra assumptions and has worse iteration complexity.", "strength_and_weaknesses": "", "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": "PXcxfITFRj3", "reviewer": "Reviewer_22pJ", "summary": "", "questions": "", "limitations": "", "rating": 6, "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": "This paper introduces a Hessian-free high-resolution Nesterov’s Acceleration approach which proves faster convergence than both Underdamped and Overdamped Langevin Dynamics. In addition, it also proves that the acceleration cannot be achieved by time-rescaling. Experiments results support the argument in optimization acceleration as well as the theoretical results. \n\nSome questions the reviewer minds:\n\n-The theorem 5.1 depends on the assumption on \\gamma and \\alpha. If such an assumption often holds in real applications.\n\n-It is will be better to demonstrate the acceleration superiority in with a more complex network and a large-scale dataset.", "summary_of_the_paper": "This stuyd focused on developing a Hessian-free high-resolution Nesterov’s Acceleration approach which proves converging faster and being  non-trivial (not a time-rescaling trick).", "summary_of_the_review": "The reviewer thinks this paper contributes a new and effective gradient based optimization solver. Overall, the reviewer would like to recommend acceptance.", "strength_and_weaknesses": "", "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": "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": {}}]

Author rebuttal:

Meta-review: The paper considers the high resolution continuous limit of Nesterov's Accelerated Gradient (NAG) algorithm and its connections to sampling (MCMC methods). The paper develops a Hessian-Free High Resolution (HFHR) ODE and injects noise into it to obtain an accelerated sampling algorithm. Further, the paper provides a discrete-time variant of the algorithm by appropriately discretizing HFHR using simple discretization schemes. For strongly log-concave potential functions (log-densities), the paper proves convergence of the order $\tilde{O}(\sqrt{d}/\epsilon)$ in Wasserstein-2 distance. In the asymptotic sense, the result matches the convergence of the underdamped Langevin algorithm; however, the paper argues that the constants in the proposed algorithm are smaller and empirically shows that the proposed algorithm is faster in practice. The main contributions of the paper are theoretical; however, the theoretical results are supplemented by numerical experiments.

Overall, the reviewers found the contributions interesting and the theoretical contributions of the paper technically sound. The main concerns that were not completely addressed were related to the presentation of the results and reproducibility of some of the numerical experiments. While both seem minor and possible to address, ultimately there was not enough support to recommend acceptance. However, the paper is solid and merits acceptance after suitable revisions. Thus, the authors are encouraged to revise the paper and resubmit it to one of the conferences in the equivalence class of ICLR.
Decision: Reject