Year: 2021

Abstract:
It is known (Shi et al., 2021) that Nesterov’s Accelerated Gradient (NAG) for optimization 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-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 W2 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 Õ( √ d/ ) iteration complexity in W2 distance, same as underdamped Langevin dynamics, but with a reduced constant. Empirical experiments are conducted to numerically verify our theoretical results.

text: √ d/ ) iteration complexity in W2 distance, same as underdamped
Langevin dynamics, but with a reduced constant. Empirical experiments are conducted to numerically verify our theoretical results.

heading: 1 INTRODUCTION
text: Optimization is a 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 a sampler based on this interplay.
In fact, approaches for sampling statistical distributions, such as gradient-based 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 optimizationbased methods. Although not entirely the same thing, optimization and sampling are 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 et al. (1998); Liu & Wang (2016); Dalalyan (2017a); Wibisono (2018); Zhang et al. (2018); Frogner & Poggio (2020); Chizat & Bach (2018); Chen et al. (2018a); Ma et al. (2021); Erdogdu & Hosseinzadeh (2021)). In fact, an unadjusted Euler-Maruyama discretization of overdamped Langevin dynamics (abbreviated as 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., 1996) and/or Langevin Monte Carlo (LMC). The convergence properties of the continuous dynamics of OLD, as well as asymptotic and non-asymptotic analyses of its discretizations have been extensively studied (e.g., Roberts et al. (1996); Villani (2008); Pavliotis (2014); Dalalyan (2017b); Durmus & Moulines (2016); Dalalyan (2017a); Durmus et al. (2019a;b); Vempala & Wibisono (2019); Cheng & Bartlett (2018); Dwivedi et al. (2019); Ma et al. (2019); Chewi et al. (2021); Erdogdu & Hosseinzadeh (2021)).
Meanwhile, the notion of acceleration is less quantified in sampling compared to that in optimization, although attention has been rapidly building up. Along this direction, one line is based on diffusion processes such as underdamped Langevin dynamics (ULD). For example, the convergence and nonasymptotics of discretized ULD have been studied in Cheng et al. (2018); Dalalyan & RiouDurand (2020); Ma et al. (2021), 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); Baudoin (2017); Eberle et al. (2019)). Another important line of research is related to accelerating particle-based approaches for optimization in probability spaces (Liu et al., 2019; Taghvaei & Mehta, 2019; Wang & Li, 2019), although we note there is no clear boundary between these two lines (e.g., Leimkuhler et al. (2018)). Additional interesting ideas also include Chen et al. (2018b); Deng et al. (2020). In general, it has been known that adding an irreversible part to the reversible dynamics of OLD1 accelerates its convergence (e.g., Hwang et al. (2005); Lelievre et al. (2013); Ohzeki & Ichiki (2015); Rey-Bellet & Spiliopoulos (2015); Duncan et al. (2016)), and this work can be viewed to be under this umbrella. Note, though, the discretization of an accelerated continuous process is also important, and it will be analyzed.
Specifically, we propose an accelerated gradient-based MCMC algorithm termed HFHR. It is motivated by a simple question: how to appropriately inject noise to NAG algorithm in discrete time, so that it is turned into an algorithm for momentum-accelerated sampling? Note we don’t add noise to the learning-rate→ 0 limit of NAG (this has been well studied in Ma et al. (2021)), because a finite-step-size discretization of this limiting ODE may not converge as fast as NAG with the same learning rate. However, we will still use continuous dynamics as intermediate steps.
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. (2021), a low-resolution limiting ODE (Su et al., 2014), albeit being a milestone leading to a new venue of research (e.g, Wibisono et al. (2016)), does not fully capture the acceleration enabled by NAG — for example, it can’t distinguish between NAG and another momentum method of heavy ball (Polyak, 1964). The main reason is, the lowresolution ODE describes the h → 0 limit of NAG, but in practice NAG uses a finite (nonzero) h. High-resolution ODE was thus proposed to include additionalO(h) terms to account for the finite h effect (Shi et al., 2021). The original form of high-resolution ODE involves Hessian of the objective function, which is computationally expensive to evaluate and store for high-dimensional problems, but this difficulty can be overcome using techniques introduced in, e.g., Alvarez et al., 2002; Attouch et al., 2020, which allows us to derive a High-Resolution and Hessian-Free limiting ODE for NAG.
Then we replace the high-resolution term’s coefficient in the HFHR ODE by a hyperparameter α ≥ 0, 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 is then discretized. We will see, both theoretically and empirically, that nonzero α 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 accuracy and even the speed of convergence (see e.g., Chen et al. (2015); Li et al. (2019); Shen & Lee (2019)), although such improvements often come with more computations per iteration. The discretization considered here is just a simple first-order scheme that uses one (full-)gradient evaluation per step, but it better utilizes the structure of HFHR dynamics than Euler-Maruyama.
Our presentation will be structured as follows. After detailing the construction of HFHR, we will analyze 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-concave 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.5.1. For our discretized HFHR algorithm, a non-asymptotic error bound will
1For irreversible-acceleration not from OLD, see e.g., Bierkens et al. (2019); Bouchard-Côté et al. (2018).
be obtained (Thm.5.2), 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.
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.5.2, Cor.5.4 & Rmk.5.5) and numerical experiments (Sec.6) are also provided to quantify the effectiveness of this idea.

heading: 2 BACKGROUND: LANGEVIN DYNAMICS FOR SAMPLING
text: Consider sampling from Gibbs measure µ whose density is dµ = 1∫ e−f(y)dy e−f(x)dx, where f : Rd 7→ 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 dynamics (abbreviated as ULD to comply with a convention of calling it underdamped Langevin). They are respectively given by
(OLD) dqt = −∇f(qt)dt+ √ 2dW t (ULD)
{ dqt = ptdt
dpt = −γptdt−∇f(qt)dt+ √ 2γdBt
where qt,pt ∈ Rd, W t,Bt are i.i.d. Wiener processes in Rd, and γ > 0 is a friction coefficient. Under mild conditions (e.g., Pavliotis (2014)), OLD converges to µ and ULD converges to dπ(q,p) = dµ(q)ν(p)dp, where ν(p) = (2π)− d 2 e−‖p‖ 2/2 , so its q marginal follows µ.
OLD and ULD are closely related. In fact, OLD is the γ → ∞ overdamping limit of ULD after time dilation (e.g., Pavliotis (2014)). However, OLD is a reversible Markov process but ULD is 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 (e.g., Dalalyan, 2017b, Lemma 1), synchronous coupling (Villani, 2008, p33-35)(Durmus et al., 2019b, Proposition 1), functional inequalities such as Poincaré’s inequality (Pavliotis, 2014, Theorem 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; Dolbeault et al., 2009; 2015; Roussel & Stoltz, 2018), coupling (Cheng et al., 2018, Theorem 5)(Dalalyan & Riou-Durand, 2020, Theorem 1)(Eberle et al., 2019, Theorem 2.3), LSI (Ma et al., 2021, Section 3.1), modified Poincaré’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 et al., 1994; Roberts et al., 1996). The non-asymptotic analysis of LMC discretization of OLD can be found in Dalalyan (2017b) and it shows the discretization achieves error, in TV distance, in Õ(d/ 2) steps. Subsequent results include Õ(d/ 2) in W2 (Durmus & Moulines, 2016), Õ(d/ ) in KL (Cheng & Bartlett, 2018), Õ(d/ ) in W2 under additional 3rd-order regularity (Durmus et al., 2019b), and Õ( √ d/ ) in W2 under additional 3rd-order regularity (Li et al., 2021). For discretized ULD, one has Õ( √ d/ ) iteration complexity in W2 (Cheng et al., 2018; Dalalyan & Riou-Durand, 2020) and Õ( √ d/ √ ) in KL (Ma et al., 2021). ULD is still generally conceived to be advantageous over OLD and sometimes understood as its momentum-accelerated version.

heading: 3 NOTATIONS AND CONDITIONS
text: We will use 2-Wasserstein distance to quantify convergence, i.e. W2(µ1, µ2) =( infπ∈Π(µ1,µ2) E(X,Y )∼π‖X − Y ‖ 2 ) 1 2 where Π(µ1, µ2) is the set of all couplings of µ1 and µ2.
Assume WLOG that 0 ∈ argminx∈Rd f(x). The following condition will also be frequently used.
A 1. (Standard Strong-Convexity and Smoothness Condition) Function f ∈ C1(Rd) : Rd 7→ R is m-stronly-convex and L-smooth, if there exist constants m,L > 0 such that ∀x,y ∈ Rd, we have
‖∇f(y)−∇f(x)‖ ≤ L‖y − x‖ and f(y) ≥ f(x) + 〈∇f(x),y − x〉+ m 2 ‖y − x‖2
For f ∈ C2, this condition is equivalent to mI ∇2f LI .

heading: 4 THE CONSTRUCTION OF HFHR DYNAMICS
text: HFHR is obtained by formulating NAG-SC as a Hessian free high-resolution ODE, lifting the highresolution term’s coefficient as a free parameter, and adding appropriate noises.
More precisely, let’s start with NAG-SC algorithm:
xk+1 = yk − s∇f(yk) (1) yk+1 = xk+1 + c(xk+1 − xk) (2)
where s is the learning rate (also known as step size), and c = 1− √ ms
1+ √ ms is a constant based on s and the strong convexity coefficient m of f ; the method also works for non-strongly-convex f though.
A high-resolution ODE description of Eq.(1) & (2) is obtained in (Shi et al., 2021, Section 2)
ÿ + √ s ( 2(1− c) s(1 + c) +∇2f(y) ) ẏ + 2 1 + c ∇f(y) = 0, (3)
which can better account for the effect of non-infinitesimal s than the s → 0 limit (note c depends on s). However, in this original form, Hessian of f is involved, which is expensive to compute and store especially for high-dimensional problems.
To obtain a Hessian-free high-resolution ODE description of Equation (1) and (2), we first turn the iteration into a ‘mechanical’ version by introducing position variable qk = yk and momentum variable pk = (yk−xk)/h. Replacing xk+1 in (1) and the first xk+1 in (2) by qk+1 and pk+1, the second xk+1 in (2) by qk − s∇f(qk), and the xk in (2) by qk and pk, we obtain{
qk+1 = qk + hpk+1 − s∇f(qk) pk+1 = cpk − c sh∇f(qk)
Now, choose γ, α and h as h = √ cs, γ = 1−ch , α = s h . It is easy to see that γ > 0, α > 0, then NAG-SC exactly rewrites as { qk+1 = qk + hpk+1 − hα∇f(qk) pk+1 = pk − hγpk − h∇f(qk) . (4)
Note the technique for bypassing the Hessian without introducing any approximation is already well studied in the literature (e.g., Alvarez et al. (2002); Attouch et al. (2020)).
So far, both h and α are actually determined by the hyperparameter s of NAG-SC. However, if we now consider α as an independent variable (i.e., ‘lift’ it) and let h → 0, we see (4) is a 1st-order discretization (with step size h) of the dynamics{
q̇ = p− α∇f(q) ṗ = −γp−∇f(q) . (5)
Note α, if inherited from NAG-SC, should be α = √ s/c = O(h), which, in a 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 α 6= O(h) can be advantageous. Before quantifying these advantages later on, we finish the construction by appropriately injecting Gaussian noises to Equation (5). This is just like how OLD can be obtained by adding noise to the
gradient flow. The right amount and structure of noise turn the ODE into a Markov process that can serve the purpose of sampling, and the detailed form of our noise is given by:{
dqt = (pt − α∇f(qt))dt+ √
2αdW t dpt = (−γpt −∇f(qt))dt+ √ 2γdBt . (6)
Here α ≥ 0, γ > 0 are constant parameters, and W t,Bt are independent standard Brownian motions in Rd. This irreversible process will be named as Hessian-Free High-Resolution(HFHR) dynamics. We write it as HFHR(α, γ) to emphasize the dependence on α and γ when needed.
Substitution into Fokker-Planck PDE shows HFHR dynamics is unbiased (proof in Appendix B.1): Theorem 4.1. π is the invariant distribution of HFHR described in Eq.(6), just like ULD.

heading: 5 THEORETICAL ANALYSIS OF THE HFHR DYNAMICS/ALGORITHM

heading: 5.1 HFHR DYNAMICS IN CONTINUOUS TIME
text: Let’s establish the exponential convergence of HFHR dynamics and its additional acceleration when compared to ULD, when the target measure has a strongly-convex and smooth potential.
Theorem 5.1. Assume Conditions A1 holds and further assume γ2 > L + m and α ≤ γ 2−L−m mγ . Denote the law of qt by µt. Then there exists κ ′ > 0 depending only on α and γ, such that
W2(µt, µ) ≤ κ′e−( m γ +mα)tW2(µ0, µ).
Detailed expression of κ′ can be found in Appendix A.
Thm. 5.1 state that HFHR converges to the target distribution exponentially fast in log-stronglyconcave and smooth setup. For ULD, Dalalyan & Riou-Durand (2020, Theorem 1) obtained exponential convergence result in 2-Wasserstein distance with rate √ m√
κ+ √ κ−1 using a simple and ele-
gant coupling approach, and showed this rate is optimal as it is achieved by the bivariate function f(x, y) = m2 x 2 + L2 y 2. In Thm 5.1, we use the same coupling approach to obtain an (asymptotically) equivalent rate √ m
2 √ κ . Since ULD is HFHR(0,γ) and our bound agrees with existing result when α = 0 and shows faster convergence for α > 0, the acceleration of HFHR in continuous time is evidenced. For example, if we set γ = 2 √ L and push α to the upper bound specified in Thm. 5.1, we obtain an O( √ L) rate in the log-strongly-concave setup. Compared with the rate in Dalalyan & Riou-Durand (2020), this is a speed-up of order κ.

heading: 5.2 HFHR ALGORITHM IN DISCRETE TIME
text: To obtain an implementable method, we now discretize the time of HFHR dynamics. As our main goal is to show the acceleration enabled by α won’t disappear after discretization (unlike a fake acceleration due to time rescaling), we’ll just analyze a 1st-order discretization (but a 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 a symmetric composition for
updating over each time interval [kh, (k+ 1)h]: xk+1 := φ h 2 ◦ψh ◦ φh2 (xk) where xk = [ qkh pkh ] , φ and ψ correspond to solution flows of split SDEs, respectively given by
φ :
{ dq = pdt
dp = −γpdt+ √ 2γdB ψ :
{ dq = −α∇f(q)dt+ √ 2αdW
dp = −∇f(q)dt ,
and φt(x0) and ψt(x0) mean x’s value after evolving φ and ψ for t time with initial condition x0.
Note that φ flow can be solved explicitly since the second equation is an Ornstein-Unlenbeck process and integrating the second equation followed by integrating the first one gives us an explicit solution
{ qt = q0 + 1−e−γt γ p0 + √ 2γ ∫ t 0 1−e−γ(t−s) γ dB(s),
pt = e −γtp0 +
√ 2γ ∫ t
0 e−γ(t−s)dB(s).
(7)
For an implementation of the stochastic integral part in Equation 7, denoting X =√ 2γ ∫ t
0 1−e−γ(t−s) γ dB(s) and Y = √
2γ ∫ t
0 e−γ(t−s)dB(s), and the covariance matrix of (X,Y )
is Cov(X,Y ) =  γh+4e−γ h2 −e−γh−3γ2 Id (1−e−γ h2 )2γ Id (1−e−γ h 2 )2
γ Id (1− e−γh)Id . As mean and covariance fully determine a Gaussian distribution, [ X Y ] = Mξ where M is the Cholesky decomposition of Cov(X,Y ), ξ is
a 2d standard Gaussian random vector, i.i.d. at each step, and φt can thus be exactly simulated.
However, ψ flow is generally not explicitly solvable unless f is a quadratic function in q. We simply choose to approximate ψh(x0) with one-step Euler-Maruyama integration ψh(x0) ≈
ψ̃h(x0) given by
{ qh = q0 − α∇f(q0)h+ √ 2αhη
ph = p0 −∇f(q0)h where η is a standard d-dimensional Gaus-
sian random vector, again i.i.d. each time ψ̃ is called.
Altogether, one step of an implementable Strang’s splitting of HFHR is hence φ h 2 ◦ ψ̃h ◦φh2 and we call this numerical scheme the HFHR algorithm, which is summarized in Algorithm 1.
Algorithm 1 HFHR Algorithm 1: Input: potential function f and its gradient ∇f , damping coefficients α and γ, step size h,
initial condition (q0,p0) 2: procedure DISCRETIZED HFHR(f,∇f, α, γ, h, q0,p0) 3: k = 0 and initialize [ q0 p0 ] 4: while not converge do 5: Generate independent standard Gaussian random vectors ηk+1 ∈ Rd, ξ 1 k+1, ξ 2 k+1 ∈ R2d
6: Run φ h 2 : [ q1 p1 ] = qkh + 1−e−γ h2γ pkh e−γ
h 2 pkh +Mξ1k+1 7: Run ψ̃h : [ q2 p2 ] = [ q1 − α∇f(q1)h+ √ 2αhηk+1 p1 −∇f(q1)h
] 8: Run φ h 2 : [ q3 p3 ] = q2 + 1−e−γ h2γ p2 e−γ
h 2 p2 +Mξ2k+1 9: [ q(k+1)h p(k+1)h ] ← [ q3 p3
] 10: k ← k + 1 11: end while 12: end procedure
Asψ in Strang splitting is replaced by a 1st-order approximation ψ̃, the method is of order 1, however with good constant. This is rigorously established by the following theorem (interested readers are referred to Appendix.D.5-D.7 and Li et al. (2021) for more technical details): Theorem 5.2. Under Assumption 1, we further assume γ − L+mγ ≥ mα and the function ∇∆f
satisfies a third-order growth condition, i.e., ∥∥∇∆f(q)∥∥ ≤ G√1 +‖q‖2,∀q ∈ Rd for someG > 0.
If (q0,p0) ∼ π0, then there exists h0, C > 0 such that when 0 < h < h0, we have
W2(µk, µ) ≤ √ 2κ′e−( m γ +mα)khW2(π0, π) + √ 2Ch (8)
where κ′ is a constant depending only on L,m, γ, α (see Appendix A for detail), µk is the law of the q marginal of the k-th iterate in Algorithm 1, and µ is the q marginal of the invariant distribution π.
In particular, C = O( √ d) and there exists b > 0, independent of α and is of order O( √ d), s.t.
C ≤ b m (α2 − α γ + 1 γ2 ). (9)
Remark 5.3. The linear growth (at infinity) condition on ∇∆f is actually not as restrictive as it appears. For example, for monomial potentials, i.e., f(x) = xp, p ∈ Z+, our linear growth condition is met when p ≤ 4, whereas a standard condition (Pavliotis, 2014, Theorem 3.1) for the existence of SDE solutions holds only when p ≤ 2. In addition, our condition is related to the Hessian Lipschitz condition commonly used in the literature (e.g., Durmus et al. (2019b); Ma et al. (2021)). Smoothness and Hessian Lipschitzness imply our condition. Meanwhile, examples that satisfy linear growth condition but are not Hessian Lipschitz exist, e.g., f(x) = x4, and thus linear growth condition is not necessarily stronger than Hessian Lipschitzness.
Inspecting the role of α in Equation (8), we see that α 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 of having a positive α > 0, at least for some α?, is reduced iteration complexity. Corollary 5.4. Consider the same assumption as in Thm. 5.2. If (q0,p0) ∼ π0, then there exists h0, C > 0 (same as that in Theorem 5.2; recallC = O( √ d)) such that for any target error tolerance
> 0, if we choose h = h? , min{h0, 2√2C }, then for < 2 √ 2Ch0, after
k? = 2 √ 2 C
m γ +mα
1
log
2 √ 2κ′W2(π0, π)
= Õ
(√ d ) . (10)
steps, we have W2(µk, µ) ≤ . Remark 5.5. Recall from Thm.5.2 that C ≤ bm (α 2− αγ + 1 γ2 ), so if we consider the minimizer α ? of
an upper bound of Cm γ +mα , α? = argminα≥0 b m2
α2−αγ + 1 γ2
1 γ+α
= √
3−1 γ . This suggests that by choosing
an optimal α > 0, one could effectively reduce iteration complexity. Note, however, that this α? may not be the true optimal one as bounds may not be tight. If they were, k?α? = (2 √ 3 − 3)k?α=0 ≈ 0.46k?α=0; i.e., steps needed by ULD (discretized by Alg.1 with α = 0) can be halved by HFHR.
Rmk.5.5 shows HFHR algorithm can lead to a similar bound on iteration complexity as ULD algorithm but with an improved constant, and thus is 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 to Appendix E, where we show, for 2D Gaussian target with condition number κ 1, the convergences of Euler discretization of ULD under optimal parameters and HFHR under suboptimal parameters are, respectively, given by (1 − 1/κ + o(1/κ))n and (1 − 2/κ + o(1/κ))n, where n is the number of iterations. The latter (HFHR) is faster despite that its hyperparameters may not be optimal.

heading: 6 NUMERICAL EXPERIMENTS
text: We now empirically study the acceleration enabled by α 6= 0 by comparing HFHR algorithm and the popular KLMC discretization of ULD (Dalalyan & Riou-Durand, 2020). For fairness, discretizations of the same order and cost are compared (Appendix F has an additional comparison).

heading: 6.1 VERIFICATION OF THEORETICAL RESULTS IN SECTION 5.2
text: This subsection numerically verifies the d 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 O( √ d) dependence would just be a consequence of using 2-Wasserstein distance for quantifying statistical accuracy. To inspect a more interesting example, we consider a potential that is no longer additive across different dimensions, namely f(x) = log (ex1 + · · ·+ exd) + 12‖x‖ 2 .
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 a 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 as a surrogate because∥∥Eµkq − Eµq∥∥ ≤W2(µk, µ) 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 √ d. To numerically verify the linear dependence on h, we work with d = 2 and ran ULD algorithm for sufficiently long time with tiny step size (h = 0.0005) to obtain 108 independent realizations and use them (as benchmark) to empirically estimate Eµq. We then set γ = 2, α = 1 and use Algorithm 1, with h ∈ {2k| − 7 ≤ k ≤ 0}. For each h, we run Th (with T = 50) iterations in Algorithm 1 to ensure the Markov chains are well-mixed and the contribution to final error from exponential decay is order-of-magnitude smaller than discretization error. The results are plotted in Figure 1a. The observed linear dependence on the step size h is consistent with Theorem 5.2.
To numerically verify the O( √ d) dependence, we extensively experiment with d ∈ {1, 2, 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, α = 1, h = 0.1, T = 10 and for each d, we run 1,000 independent realizations of HFHR algorithm for Th = 100 iterations. Experiment results are plotted in Figure 1b and the linear trend demonstrates that the bound in Thm. 5.2 is tight in terms of d-dependence.
The final experiment compares Algorithm 1 with ULD algorithm in terms of iteration complexity. 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 accordingly. To do so, we push both ULD and HFHR to their respective largest h 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.5 suggests that with appropriately chosen α, HFHR algorithm effectively reduces the constant factor of the iteration complexity, implying accelerated sampling. To just provide one empirical verification of this improvement, we choose d = 10 and use the error of mean
∥∥Eµkq − Eµq∥∥ to measure sampling accuracy. The benchmark, i.e., Eµq, 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 wellmixed. The initial measure is chosen as a Dirac measure at (100 × 1d,0d), where 1d,0d are d-dimensional vectors filled with 1 and 0 respectively. We pick threshold = 0.1, and for each α ∈ {0, 0.001, 0.002, 0.005, 0.01, 0.02, 0.05, 0.1, 0.2, 0.5, 1, 2, 5, 10, 20, 50, 100}, we try all combinations of (γ, h) ∈ {0.1, 0.2, 0.5, 1, 2, 5, 10, 20, 50, 100} × { 0.1× [50] } for Algorithm 1 (we also run ULD algorithm when α = 0), and empirically find the best combination that requires the fewest iterations to meet
∥∥Eµkq − Eµq∥∥ ≤ . We find that h = 5 already surpasses the stability limit of ULD algorithm, hence the range of step size covers the largest step size that are practically useable for ULD algorithm. Experiments are repeated with 100 different seeds to further reduce variance.
The results are shown in Figure 2. When α > 0, HFHR algorithm consistently outperforms ULD algorithm (note it also does so when α = 0 because HFHR uses a efficiency-wise-comparable but more refined discretization than ULD algorithm). In particular, when α = 0.5 and 1, which are
empirically the best values we found for this experiment, HFHR algorithm achieves the specified -closeness nearly 6× times faster than ULD algorithm, and its decreased mixing time (compared to α = 0 for the same algorithm) is consistent with the≈ 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.

heading: 6.2 BAYESIAN NEURAL NETWORK
text: 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 training, which amounts to a practical, high-dim., multi-modal example of sampling tasks. It no longer satisfies the conditions of our analysis, and our goal is to show HFHR still accelerates. We use fullyconnected 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-parameter for Algorithm 1 and ULD algorithm are systematically investigated. For each pair (γ, α) ∈ {0.1, 0.5, 1, 5, 10, 50, 100}2, 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 γ, we plot the negative log likelihood of HFHR algorithm (with different α choices) and ULD algorithm on training and test data in Figure 3. Cases where α 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-concave assumption required in Theorem 5.2 does not hold for 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. It showcases the applicability and effectiveness of HFHR as a general sampling algorithm.

heading: 7 CONCLUSION AND DISCUSSION
text: This paper proposes HFHR, an accelerated gradient-based MCMC method for sampling. To demonstrate 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 is often already considered as an acceleration of Overdamped Langevin Dynamics. As HFHR is based on a new perspective, which is to turn NAG-SC optimizer 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 performance of SG-HFHR? Can the generalization ability of HFHR trained learning models (e.g., BNN) be quantified, and how does it compare with that by LMC and/or KLMC? These will be future work.

heading: A ADDITIONAL NOTATIONS
text: We introduce a few notations that are used in the main text as well as some proof. When ∇f is L-Lipschitz, the drift term [ p− α∇f(q) −γp−∇f(q) ] in HFHR dynamics is also L′-Lipschitz, as proved in
Lemma D.3, where
L′ = √ 2 max {√ 1 + α2 max { 1√ 2 , L } , √ 1 + γ2 } .
We show in Lemma D.5 that a linear-transformed HFHR dynamics satisfies the nice contraction property, the linear transformation P we use is defined as
P = [ γI I 0 √ 1 + αγI ] ∈ R2d×2d.
Denote the largest and the smallest singular value of P by
σmax =
√ αγ
2 + γ2 2 +
√ α2γ2 − 2αγ3 + 4αγ + γ4 + 4
2 + 1,
σmin =s
√ αγ
2 + γ2 2 − √ α2γ2 − 2αγ3 + 4αγ + γ4 + 4 2 + 1
and its condition number by
κ′ = σmax σmin =
√√√√√ αγ2 + γ22 + √ α2γ2−2αγ3+4αγ+γ4+4 2 + 1
αγ 2 +
γ2 2 − √ α2γ2−2αγ3+4αγ+γ4+4 2 + 1 .
The rate λ′ of exponential convergence of transformed HFHR dynamics is characterized in Lemma D.5 and is defined as
λ′ = min
{ m
γ + αm, γ2 − L γ } given that γ2 > L.

heading: B PROOFS FOR THE CONTINUOUS DYNAMICS
text: Notations and definitions can be found in Sec.3.

heading: B.1 PROOF OF THEOREM 4.1
text: Proof. The Fokker-Plank equation of HFHR is given by
∂tρt = −∇x ·
([ p
−∇f(q)
] ρt ) + α ( ∇q · (∇f(q)ρt) + ∆qρt ) + γ ( ∇p · (pρt) + ∆pρt
) where∇x = (∇q,∇p). For π ∝ e−f(q)− 1 2‖p‖ 2 , we have
∇x ·
([ p
−∇f(q)
] π ) = 〈 [
p −∇f(q)
] ,∇xπ〉 = 0,
∆qπ = −∇q · (π∇f(q)) ∆pπ = −∇p · (πp)
Therefore ∂tπ = 0 and hence π is the invariant distribution of HFHR.

heading: B.2 PROOF OF THEOREM 5.1
text: Proof. Consider two copies of HFHR that are driven by the same Brownian motion{ dqt = (pt − α∇f(qt))dt+ √ 2αdB1t
dpt = (−γpt −∇f(qt))dt+ √ 2γdB2t ,
{ dq̃t = (p̃t − α∇f(q̃t))dt+ √ 2αdB1t
dp̃t = (−γp̃t −∇f(q̃t))dt+ √ 2γdB2t ,
where we set (q̃0, p̃0) ∼ π, p0 = p̃0 and q0 such that W 22 (µ0, µ) = E [ ‖q0 − q̃0‖22 ] , q0 ∼ µ0
Denote [ φt ψt ] = P [ qt − q̃t pt − p̃t ] where P is defined in Appendix A. By Lemma D.5 and the assumption on α, γ, we have ∥∥∥∥∥ [ φt ψt ]∥∥∥∥∥ 2 ≤ e−2( m γ +mα)t ∥∥∥∥∥ [ φ0 ψ0 ]∥∥∥∥∥ 2 .
Therefore we obtain
W 22 (µt, µ) = inf (qt,q̃t)∼Π(µt,µ) E‖qt − q̃t‖ 2
≤ inf (qt,q̃t)∼Π(µt,µ),(pt,p̃t)∼Π(νt,ν) E ∥∥∥∥∥ [ qt − q̃t pt − p̃t ]∥∥∥∥∥ 2
≤E‖P−1‖22 ∥∥∥∥∥ [ φt ψt ]∥∥∥∥∥ 2
≤E‖P−1‖22e −2(mγ +mα)t ∥∥∥∥∥ [ φ0 ψ0 ]∥∥∥∥∥ 2
≤(κ′)2e−2( m γ +mα)t ∥∥∥∥∥ [ q0 − q̃0 p0 − p̃0 ]∥∥∥∥∥ 2
=(κ′)2e−2( m γ +mα)tW 22 (µ0, µ)
Taking square root yields the desired result.

heading: C ARBITRARY LONG TIME DISCRETIZATION ERROR OF ALGORITHM 1
text: Theorem C.1. Under Conditions A1 and further assume the function∇∆f grows at most linearly, i.e., ∥∥∇∆f(q)∥∥ ≤ G√1 +‖q‖2,∀q ∈ Rd. Also suppose γ in HFHR dynamics satisfy γ2 > L. Then there exist C, h0 > 0, such that for 0 < h ≤ h0, we have(
E‖xk − x̄k‖2 ) 1 2 ≤ Ch
where x̄k is the k-th iterate of Algorithm 1 with step size h starting from x0, xk is the solution of HFHR dynamics at time kh, starting from x0. This result holds uniformly for all k ≥ 0 and k can go to∞. In particular, C = O( √ d) and if γ − L+mγ ≥ mα, then there exists b > 0, independent of α and is of order O( √ d), such that
C ≤ b m (α2 − α γ + 1 γ2 ). (11)
Proof. Denote tk = kh, the solution of the HFHR dynamics at time t by x0,x0(t), the k-th iterates of the Strang’s splitting method of HFHR dynamics by x̄0,x0(kh). Both x0,x0(t) and x̄0,x0(kh)
start from the same initial value x0. The linear transformation P defined in Appendix A, transforms the solution of HFHR dynamics into y0,Px0(t) = Px0,x0(t) and the Strang’s splitting discretization of HFHR into ȳ0,Px0(t) = P x̄0,x0(t).
For the ease of notation, we write y0,y0(tk) as yk and ȳ0,y0(tk) as ȳk. We have the following identity
E ∥∥yk+1 − ȳk+1∥∥2 =E∥∥∥ytk,yk(h)− ȳtk,ȳk(h)∥∥∥2
=E ∥∥∥ytk,yk(h)− ytk,ȳk(h) + ytk,ȳk(h)− ȳtk,ȳk(h)∥∥∥2
=E ∥∥∥ytk,yk(h)− ytk,ȳk(h)∥∥∥2︸ ︷︷ ︸
1
+E ∥∥∥ytk,ȳk(h)− ȳtk,ȳk(h)∥∥∥2︸ ︷︷ ︸
2 +2E 〈 ytk,yk(h)− ytk,ȳk(h),ytk,ȳk(h)− ȳtk,ȳk(h) 〉 ︸ ︷︷ ︸
3
By Lemma D.5, when 0 < h < 12λ′ , term 1 can be upper bounded as
E ∥∥∥ytk,yk(h)− ytk,ȳk(h)∥∥∥2 ≤e−2λ′hE‖yk − ȳk‖2
≤ ( 1− 2λ′h+ 2(λ′)2h2 ) E‖yk − ȳk‖ 2
≤ ( 1− λ′h ) E‖yk − ȳk‖ 2
where the second inequality is due to e−x ≤ 1− x+ x 2
2 ,∀x > 0.
For term 2 , we have by Lemma D.8 that
E ∥∥∥ytk,ȳk(h)− ȳtk,ȳk(h)∥∥∥2 ≤ σ2max E∥∥xtk,x̄k(h)− x̄tk,x̄k(h)∥∥2 ≤ σ2max C22h3
where σmax is the largest singular value of matrix P .
For term 3 , we have by Lemma D.1 that
2E 〈 ytk,yk(h)− ytk,ȳk(h),ytk,ȳk(h)− ȳtk,ȳk(h) 〉 =2E 〈 yk − ȳk + z,ytk,ȳk(h)− ȳtk,ȳk(h)
〉 = 2E 〈 yk − ȳk,ytk,ȳk(h)− ȳtk,ȳk(h) 〉 ︸ ︷︷ ︸
3a
+ 2E 〈 z,ytk,ȳk(h)− ȳtk,ȳk(h) 〉 ︸ ︷︷ ︸
3b
For term 3a , by the tower property of conditional expectation, we have
2E 〈 yk − ȳk,ytk,ȳk(h)− ȳtk,ȳk(h) 〉 =2E E[〈yk − ȳk,ytk,ȳk(h)− ȳtk,ȳk(h)〉 ∣∣∣∣Fk ]
=2E 〈 yk − ȳk,E [ ytk,ȳk(h)− ȳtk,ȳk(h) ∣∣∣∣Fk ]〉
≤2 √ E‖yk − ȳk‖ 2 √√√√√E ∥∥∥∥∥∥E [ ytk,ȳk(h)− ȳtk,ȳk(h) ∣∣∣∣Fk ]∥∥∥∥∥∥ 2 ≤2 √ E‖yk − ȳk‖ 2 √√√√√σ2maxE ∥∥∥∥∥∥E [ xtk,x̄k(h)− x̄tk,x̄k(h) ∣∣∣∣Fk ]∥∥∥∥∥∥ 2
≤2 √ E‖yk − ȳk‖ 2 √ σ2maxC 2 1h 4
≤2σmaxC1 √ E‖yk − ȳk‖ 2 h2.
For term 3b , when 0 < h < 14L′′we have by Lemma D.1 and Lemma D.8
2E 〈 z,ytk,ȳk(h)− ȳtk,ȳk(h) 〉 ≤2 √ E‖z‖2 √ E ∥∥∥ytk,ȳk(h)− ȳtk,ȳk(h)∥∥∥2
=2 √ E‖z‖2 √√√√√E E[∥∥∥ytk,ȳk(h)− ȳtk,ȳk(h)∥∥∥2 ∣∣∣∣Fk ]
=2 √ E‖z‖2 √√√√√σ2maxE E[∥∥xtk,x̄k(h)− x̄tk,x̄k(h)∥∥2 ∣∣∣∣Fk ] ≤2σmax √ C̃E‖yk − ȳk‖ 2 h2 √ C22h 3
≤2σmaxC2 √ C̃ √ E‖yk − ȳk‖ 2 h 5 2
where C̃ = 2 ( L′′ )2 = 2(κ′)2 ( L′ )2 is from Lemma D.1 and Lemma D.3.
Recall bothC1 andC2 depend on‖xk‖ and we would like to upper bound this term. To this end, consider x̃(t), a solution of HFHR dynamics with initial value x̃0 that follows the invariant distribution x̃0 ∼ π and realizes W2(π0, π), i.e., E‖x̃0 − x0‖2 = W 22 (π0, π).
Denote x̃k = x̃(kh) and ek = ( E‖yk − ȳk‖ 2 ) 1 2 , we then have
E‖x̄k‖2 =E‖xk + x̄k − xk‖2
≤2E‖xk‖2 + 2E‖x̄k − xk‖2
≤4E‖x̃k‖2 + 4E‖x̃k − xk‖2 + 2E‖x̄k − xk‖2 =4E‖x̃k‖2 + 4E ∥∥∥P−1P (x̃k − xk)∥∥∥2 + 2E∥∥∥P−1P (x̄k − xk)∥∥∥2
≤4 (∫
Rd ‖q‖2 dµ+ d
) + 4 σ2min E ∥∥P (x̃k − xk)∥∥2 + 2 σ2min E‖ȳk − yk‖ 2
(i) ≤4 (∫
Rd ‖q‖2 dµ+ d
) + 4
σ2min e−2λ
′khE ∥∥P (x̃0 − x0)∥∥2 + 2
σ2min e2k ≤4 (∫
Rd ‖q‖2 dµ+ d
) + 4κ2W 22 (π0, π) + 2
σ2min e2k
,Fe2k +G
where (i) is due to Lemma D.5. Recall from Lemma D.8, we have
C1 ≤ A1 √ E‖x̄k‖2 +B1 ≤ A1 √ Fek + (A1 √ G+B1) , U1ek + V1
C2 ≤ A2 √ E‖x̄k‖2 +B2 ≤ A2 √ Fek + (A2 √ G+B2) , U2ek + V2
where
A1 =(L+G) max{α+ 1.25, γ + 1}(1.74 + 0.71α) B1 =(L+G) max{α+ 1.25, γ + 1} [ 0.5α+ (1.26 √ α+ 1.14α √ α+ 2.32 √ γ) √ hd ] A2 =Lmax{α+ 1.25, γ + 1}(1.92 + 2.30αL) √ h B2 =Lmax{α+ 1.25, γ + 1}(2.60 √ α+ 3.34 √ γh) √ d
Combine the above and bounds for terms 1 , 2 , 3a and 3b , we then obtain e2k+1 ≤(1− λ′h)e2k + σ2max C22h3 + 2σmaxC1ekh2 + 2σmaxC2 √ C̃ekh 5 2
≤(1− λ′h)e2k + σ2max2(U22 e2k + V 22 )h3 + 2σmax(U1ek + V1)ekh2 + 2σmax(U2ek + V2) √ C̃ekh 5 2
= ( 1− λ′h+ 2σ2maxU22h3 + 2σmaxU1h2 + 2σmaxU2 √ C̃h 5 2 ) e2k
+ ( 2σmaxV1 + 2σmaxV2 √ C̃h ) ekh 2 + 2σ2maxV 2 2 h 3
≤ ( 1− λ′h+ 2σ2maxU22h3 + 2σmaxU1h2 + 2σmaxU2 √ C̃h 5 2 ) e2k + λ′
8 he2k
+ 2 ( 2σmaxV1 + 2σmaxV2 √ C̃h )2
λ′ h3 + 2σ2maxV 2 2 h 3
= ( 1− 7
8 λ′h+ 2σ2maxU 2 2h 3 + 2σmaxU1h 2 + 2σmaxU2
√ C̃h 5 2 ) e2k
+ 2 ( 2σmaxV1 + 2σmaxV2 √ C̃h )2
λ′ + 2σ2maxV 2 2
h3
(i) ≤(1− 1 2 λ′h)e2k + 2 ( 2σmaxV1 + 2σmaxV2 √ C̃h )2 λ′ + 2σ2maxV 2 2 h3 ,(1− 1
2 λ′h)e2k +Kh 3
where (i) is due to h < min{h1, h2, h3} and
h1 =
√ λ′
4 √ 2κ′Lmax{α+ 1.25, γ + 1}(1.92 + 2.30αL) ,
h2 = λ′
16 √ 2κ′(L+G) max{α+ 1.25, γ + 1}(1.74 + 0.71α) ,
h3 = λ′
8κ′Lmax{α+ 1.25, γ + 1}(1.92 + 2.30αL) .
Unfolding the above inequality, we arrive at
e2k ≤ ( 1− λ ′
2 h
)k e20 + ( 1 + (1− λ ′
2 h) + · · ·+ (1− λ
′
2 h)k−1
) Kh3
(i) ≤Kh3 ∞∑ i=0 ( 1− λ ′ 2 h )i = 2K
λ′ h2
where (i) is due to ek = 0. Therefore
( E‖xk − x̄k‖2 ) 1 2 = ( E ∥∥∥P−1(yk − ȳk)∥∥∥2) 12 ≤ 1σmin ek ≤ 1σmin √ 2K λ′ h
Collecting all the constants and we have
1
σmin
√ 2K
λ′ ≤8κ
′
λ′ (L+G) max{α+ 1.25, γ + 1}(1.74 + 0.71α) √∫ Rd ‖q‖2 dµ+ d+ κ′W2(π0, π)  + 4κ′
λ′ (L+G) max{α+ 1.25, γ + 1}
( 0.5α+ (1.26 √ α+ 1.14α √ α+ 2.32 √ γ) √ d )
+ 8κ′√ λ′ (√ κ′L′√ λ′ + 1 ) Lmax{α+ 1.25, γ + 1}(1.92 + 2.30αL) √∫ Rd ‖q‖2 dµ+ d+ κ′W2(π0, π)  +
4κ′√ λ′ (√ κ′L′√ λ′ + 1 ) Lmax{α+ 1.25, γ + 1}(2.60 √ α+ 3.34 √ γ) √ d
, C
It is clear that in terms of the dependence on dimension d, we have C = O( √ d). In the regime where γ 2−L γ ≥ m γ +mα, then λ ′ = mγ +mα. Recall the definition of κ ′ and there exist A′, B′ > 0 such that κ′ ≤ A′ √ α+B′. It follows that
C ≤ a1α 3 + a2α 5 2 + a3α 2 + a4α 3 2 + a5α+ a6α 1 2 + a7
λ′ ≤ b
α3 + 1γ3
λ′ = b α3 + 1γ3 m γ +mα = b m (α2−α γ α+ 1 γ2 )
for some positive constants a1, a2, a3, a4, a5, a6, a7, b > 0 and independent of α, in particular, we have b = O( √ d).

heading: C.1 PROOF OF THEOREM 5.2
text: Proof. Denote the k-th iterate of the Strang’s splitting method of HFHR by x̄k with time step h, the solution of HFHR dynamics at time hk by xk. Both x̄k and xk start from x0 = [ q0 p0 ] . Also denote
the solution of HFHR dynamics starting from x̃0 at time kh by x̃k where x̃0 = [ q̃0 p0 ] , (q̃0, p̃0) ∼ π
and E ∥∥∥∥∥ [ q0 − q̃0 p0 − p̃0 ]∥∥∥∥∥ 2 = W 22 (π0, π). Since π is the invariant distribution of HFHR dynamics, it
follows that x̃k ∼ π. By Lemma D.5 and Theorem C.1, we have
W 22 (µk, µ) = inf ξ∈Π(µk,µ) E(q1,q2)∼ξ‖q1 − q2‖ 2
≤ inf ξ∈Π(πk,π) E(x1,x2)∼ξ‖x1 − x2‖ 2
≤E‖x̄k − x̃k‖2 ≤2C2h2 + 2E ∥∥∥P−1P (xk − x̃k)∥∥∥2
≤2C2h2 + 2‖P−1‖22E ∥∥P (xk − x̃k)∥∥2 ≤2C2h2 + 2‖P−1‖22e−2λ ′khE
∥∥P (x0 − x̃0)∥∥2 ≤2C2h2 + 2(κ′)2e−2λ ′khW 22 (π0, π)
Take square root on both sides and apply √ a2 + b2 ≤ a+ b, we obtain
W2(µk, µ) ≤ √ 2Ch+ √ 2κ′e−λ ′khW2(π0, π).

heading: C.2 PROOF OF COROLLARY 5.4
text: Proof. By Theorem 5.2, we have
W2(µk, µ) ≤ √ 2Ch+ √ 2κ′e−λ ′khW2(π0, π).
Given any target accuracy > 0, if we run the Strang’s splitting method of HFHR with h? = min{h0, 2√2C }, then after k ? = 1λ′ max{ 1 h0 , 2 √ 2C } log 2 √ 2κ′W2(π0,π) , we have
W2(µk? , µ) ≤ √ 2Ch+ √ 2κ′e−λ ′khW2(µ0, µ) ≤
2 + 2 = .
Recall C = O( √ d), when high accuracy is needed, e.g. < 2 √ 2Ch0, the iteration complexity to reach -accuracy under 2-Wasserstein distance is k? = O( √ d log 1 ) = 2 √
2 Cλ′ 1 log
2 √
2κ′W2(π0,π) = Õ( √ d ). Recall from Theorem C.1, C ≤
b m (α 2 − αγ + 1 γ2 ), we
have
C λ′ ≤ b m2
α2 − αγ + 1 γ2
1 γ + α
Denote g(α) = bm2 α2−αγ + 1 γ2
1 γ+α
, simple calculation shows that α? = argminα≥0 g(α) = √ 3−1 γ =
O( 1γ ).

heading: D TECHNICAL/AUXILIARY LEMMAS AND THEIR PROOFS

heading: D.1 DEPENDENCE OF ERROR OF SDE ON INITIAL VALUES
text: Lemma D.1. Consider the following two SDE with different initial condition{ dxt = a(xt)dt+ σdW t,
x(0) = x0
{ dyt = a(yt)dt+ σdW t,
y(0) = y0
where a(u) ∈ Rd is L-Lipschitz, and σ ∈ Rn×n is a constant matrix. For 0 < h < 14L , we have the following representation
xh − yh = x0 − y0 + z
with
E‖z‖2 ≤ 2L2‖x0 − y0‖ 2 h2
Proof. Let z = (xh − yh)− (x0 − y0) = ∫ h 0 a(xs)− a(ys)ds. Ito’s lemma readily implies that
E‖xh − yh‖ 2 =‖x0 − y0‖ 2
+ 2E ∫ h
0
〈xs − ys,a(xs)− a(ys)〉ds
≤‖x0 − y0‖ 2 + 2L ∫ h 0 E‖xs − ys‖ 2 ds
By Gronwall’s inequality, it follows that
E‖xh − yh‖ 2 ≤‖x0 − y0‖ 2 e2Lh ≤ 2‖x0 − y0‖ 2 , for 0 < h <
1
4L
and
E‖z‖2 = ∥∥∥∥∥∥E [∫ h 0 a(xs)− a(ys)ds ]∥∥∥∥∥∥ 2 ≤ (∫ h 0 ∥∥∥E [a(xs)− a(ys)]∥∥∥ ds )2
≤ ∫ h
0
12ds ∫ h 0 ∥∥∥E [a(xs)− a(ys)]∥∥∥2 ds ≤h ∫ h
0
E ∥∥a(xs)− a(ys)∥∥2 ds
≤L2h ∫ h
0
E‖xs − ys‖ 2 ds
≤2L2‖x0 − y0‖ 2 h2

heading: D.2 GROWTH BOUND OF SDE WITH ADDITIVE NOISE
text: Lemma D.2. Consider the following SDE with constant diffusion{ dxt = a(xt)dt+ σdW t,
x(0) = x0
where a(x) ∈ Rd is L-smooth, i.e., |a(y) − a(x)| ≤ L|y − x|, a(0) = 0 and σ ∈ Rd×d is a constant matrix independent of time t and xt. Then for 0 < h < 14L , we have
E‖xh − x0‖2 ≤ 2.57 ( ‖σ‖2F + 2hL2‖x0‖ 2 ) h.
Proof. We have
E‖xh − x0‖2 =E ∥∥∥∥∥ ∫ h
0
a(xt)dt+ ∫ h 0 σdW t ∥∥∥∥∥ 2
≤2E ∥∥∥∥∥ ∫ h
0
a(xt)dt ∥∥∥∥∥ 2 + 2E ∥∥∥∥∥ ∫ h
0
σdW t ∥∥∥∥∥ 2
(i) =2E ∥∥∥∥∥ ∫ h
0
a(xt)dt ∥∥∥∥∥ 2 + 2 ∫ h 0 ‖σ‖2F dt
≤2E (∫ h 0 ∥∥a(xt)∥∥ dt)2 + 2h‖σ‖2F
≤2E (∫ h 0 ∥∥a(xt)− a(x0)∥∥ dt+ ∫ h 0 ∥∥a(x0)∥∥ dt)2 + 2h‖σ‖2F
≤2E (L∫ h 0 ‖xt − x0‖ dt+ h ∥∥a(x0)∥∥)2 + 2h‖σ‖2F ≤4E L2(∫ h 0 ‖xt − x0‖ dt )2 + h2 ∥∥a(x0)∥∥2 + 2h‖σ‖2F
(ii) ≤ 2h‖σ‖2F + 4h2 ∥∥a(x0)∥∥2 + 4L2h∫ h
0
E‖xt − x0‖2 dt
where (i) is due to Ito’s isometry, (ii) is due to Cauchy-Schwarz inequality and ‖σ‖F is the Frobenius norm of σ. By Gronwall’s inequality, we obtain
E‖xh − x0‖2 ≤ ( 2h‖σ‖2F + 4h2 ∥∥a(x0)∥∥2) exp{4L2h2} .
Since ∥∥a(x0)∥∥ =∥∥a(x0)− a(0)∥∥ ≤ L‖x0‖, when 0 < h < 14L , we finally reach at
E‖xh − x0‖2 ≤ 2 ( ‖σ‖2F + 2hL2‖x0‖ 2 ) e 1 4h ≤ 2.57 ( ‖σ‖2F + 2hL2‖x0‖ 2 ) h.

heading: D.3 LIPSCHITZ CONTINUITY OF THE DRIFT OF HFHR DYNAMICS
text: Lemma D.3. Assume ∇f is L-Lipschitz, i.e. ∥∥∇f(x)−∇f(y)∥∥ ≤ L‖x− y‖, then the drift term of HFHR dynamics [ p− α∇f(q) −γp−∇f(q)
] is L′-Lipschitz, where L′ , √ 2 max{ √ 1 + α2 max{ 1√
2 , L},
√ 1 + γ2}. Let P be defined in Ap-
pendix A and [ φ ψ ] = P [ q p ] , then [ φ ψ ] satisfies the following SDE
[ dφ dψ ] = P [ p(φ,ψ)− α∇f(q(φ,ψ)) −γp(φ,ψ)−∇f(q(φ,ψ)) ] dt+ P [√ 2αI 0 0 √ 2γI ] [ dW dB ] and the drift term
P [ p(φ,ψ)− α∇f(q(φ,ψ)) −γp(φ,ψ)−∇f(q(φ,ψ)) ] is L′′-Lipschitz, where L′′ = κ′L′ and κ′ is the condition number of P .
Proof. By direct computation and Cauchy-Schwarz inequality, we have∥∥∥∥∥ [ p1 − α∇f(q1) −γp1 −∇f(q1) ] − [ p2 − α∇f(q2) −γp2 −∇f(q2) ]∥∥∥∥∥ =
√∥∥∥−α (∇f(q1)−∇f(q2))+ (p1 − p2)∥∥∥2 +∥∥∥− (∇f(q1)−∇f(q2))− γ(p1 − p2)∥∥∥2 ≤ √ 2α2 ∥∥∇f(q1)−∇f(q2)∥∥+ 2‖p1 − p2‖2 + 2∥∥∇f(q1)−∇f(q2)∥∥+ 2γ2‖p1 − p2‖2
≤ √
(2α2L2 + 2L2)‖q1 − q2‖+ (2 + 2γ2)‖p1 − p2‖ 2
≤ √ 2 max{L √ 1 + α2, √ 1 + γ2} ∥∥∥∥∥ [ q1 − q2 p1 − p2 ]∥∥∥∥∥ ≤ √ 2 max{ √
1 + α2 max{ 1√ 2 , L},
√ 1 + γ2} ∥∥∥∥∥ [ q1 − q2 p1 − p2 ]∥∥∥∥∥ ,L′ ∥∥∥∥∥ [ q1 − q2 p1 − p2
]∥∥∥∥∥ By Ito’s lemma, we have[
dφ dψ
] = P [ p(φ,ψ)− α∇f(q(φ,ψ)) −γp(φ,ψ)−∇f(q(φ,ψ)) ] dt+ P [√ 2αI 0 0 √ 2γI ] [ dW dB ]
Using the Lipschitz constant obtained for the drift of HFHR, we further have∥∥∥∥∥P [ p(φ1,ψ1)− α∇f(q(φ1,ψ1)) −γp(φ1,ψ1)−∇f(q(φ1,ψ1)) ] − P [ p(φ2,ψ2)− α∇f(q(φ2,ψ2)) −γp(φ2,ψ2)−∇f(q(φ2,ψ2)) ]∥∥∥∥∥ ≤σmax ∥∥∥∥∥ [ p1 − α∇f(q1) −γp1 −∇f(q1) ] − [ p2 − α∇f(q2) −γp2 −∇f(q2)
]∥∥∥∥∥ ≤σmaxL′ ∥∥∥∥∥ [ q1 − q2 p1 − p2
]∥∥∥∥∥ ≤σmaxL′ ∥∥∥∥∥P−1 [ φ1 − φ2 ψ1 −ψ2
]∥∥∥∥∥ ≤σmaxL′ 1
σmin ∥∥∥∥∥ [ φ1 − φ2 ψ1 −ψ2 ]∥∥∥∥∥ =κ′L′ ∥∥∥∥∥ [ φ1 − φ2 ψ1 −ψ2
]∥∥∥∥∥ where σmax, σmin and κ′ are the largest, smallest singular values and the condition number (w.r.t. 2-norm) of matrix P .
Remark D.4. The following inequalities associated with L′ will turn out to be useful in many proofs
L′ ≥ 1, L′ ≥ √ 2γ, L′ ≥ √ 2α,L ≥ √ 2L and L′ ≥ √ 2αL.

heading: D.4 CONTRACTION OF (TRANSFORMED) HFHR DYNAMICS
text: Lemma D.5. Suppose f is L-smooth, m-strongly convex and γ2 > L. Consider two copies of HFHR dynamics [ qt pt ] , [ q̃t p̃t ] (driven by the same Brownian motion) with initialization [ q0 p0 ] , [ q̃0 p̃0 ] respectively, then we have ∥∥∥∥∥P [ qt − q̃t pt − p̃t ]∥∥∥∥∥ ≤ e−λ′t ∥∥∥∥∥P [ q0 − q̃0 p0 − p̃0
]∥∥∥∥∥ where P = [ γI I 0 √ 1 + αγI ] and λ′ = min{mγ + αm, γ2−L γ }.
Proof. Consider two copies of HFHR that are driven by the same Brownian motion{ dqt = (pt − α∇f(qt))dt+ √ 2αdB1t
dpt = (−γpt −∇f(qt))dt+ √ 2γdB2t ,
{ dq̃t = (p̃t − α∇f(q̃t))dt+ √ 2αdB1t
dp̃t = (−γp̃t −∇f(q̃t))dt+ √ 2γdB2t .
Based on Taylor’s expansion, the difference of the two copies is expressed as
d
dt [ qt − q̃t pt − p̃t ] =− [ αHt −I Ht γI ] [ qt − q̃t pt − p̃t ] , −A [ qt − q̃t pt − p̃t ] where Ht = ∫ 1 0 ∇2f(q̃t + s(q − q̃t))ds. Denote the eigenvalues of Ht by ηi, 1 ≤ i ≤ d, by strong convexity and smoothness assumption on f , we have m ≤ ηi ≤ L, 1 ≤ i ≤ d.
Denote [ φt ψt ] = P [ qt − q̃t pt − p̃t ] and consider Lt = 12 ∥∥∥∥∥ [ φt ψt ]∥∥∥∥∥ 2 , we have
d dt Lt =− [ φt ψt ]T PAP−1 [ φt ψt ] =− [ φt ψt ]T 1 2 (PAP−1 + (P−1)TATPT ) [ φt ψt
] =− [ φt ψt ]T 1 γ [ (1 + αγ)Ht 0d×d 0d×d γ 2I −Ht ] [ φt ψt
] ,− [ φt ψt ]T B(α) [ φt ψt
] It is easy to see that
λmin(B(α)) = min i=1,2,··· ,d {min{ηi γ + αηi, γ − ηi γ }} ≥ min{m γ + αm, γ2 − L γ } , λ′.
Therefore we have ddtLt ≤ −2λminB(α)Lt ≤ −2λ ′Lt. By Gronwall’s inequality, we obtain∥∥∥∥∥ [ φt ψt ]∥∥∥∥∥ 2 ≤ e−2λ ′t ∥∥∥∥∥ [ φ0 ψ0 ]∥∥∥∥∥ 2 .
and the desired inequality follows by taking square root.

heading: D.5 LOCAL ERROR BETWEEN THE EXACT STRANG’S SPLITTING METHOD AND HFHR DYNAMICS
text: Lemma D.6. Assume f is L-smooth and 0 ∈ argminx∈Rd f(x), i.e. ∇f(0) = 0. If 0 < h ≤ 14L′ , then compared with the HFHR dynamics, the exact Strang’s splitting method has local mathematical expectation of deviation of order p1 = 2 and local mean-squared error of order p2 = 2, i.e. there exist constants Ĉ1, Ĉ2 > 0 such that∥∥Ex(h)− Ex̂(h)∥∥ ≤ Ĉ1hp1(
E [∥∥x(h)− x̂(h)∥∥2]) 12 ≤ Ĉ2hp2
where x(h) = [ q(h) p(h) ] is the solution of the HFHR dynamics with initial value x0 = [ q0 p0 ] and
x̂(h) = [ q̂(h) p̂(h) ] is the solution of the implementable Strang’s splitting with initial value x0 = [ q0 p0 ] , p1 = 2 and p2 = 2. More concretely, we have
Ĉ1 = Lmax{α+ 1.25, γ + 1} ( 1.74‖x0‖+ (1.26 √ α+ 2.84 √ γ) √ hd ) ,
Ĉ2 = Lmax{α+ 1.25, γ + 1} ( 1.92‖x0‖+ (1.30 √ α+ 3.22 √ γ) √ hd ) .
Proof. The exact Strang’s splitting integrator with step size h reads as φ h 2 ◦ ψh ◦ φh2 where
φ :
{ dq = pdt
dp = −γpdt+ √ 2γdB ψ :
{ dq = −α∇f(q)dt+ √ 2αdW
dp = −∇f(q)dt .
The φ flow can be explicitly solved and the solution is{ q(t) = q0 + 1−e−γt γ p0 + √ 2γ ∫ t 0 1−e−γ(t−s) γ dB(s)
p(t) = e−γtp0 + √ 2γ ∫ t 0 e−γ(t−s)dB(s)
.
The ψ flow can be written as{ q(t) = q0 − ∫ t 0 α∇f(q(s))ds+ √ 2α ∫ t 0 dW (s)
p(t) = p0 − ∫ t 0 ∇f(q(s))ds
.
The solution of one-step exact Strang’s splitting integrator with step size h can be written as q3 = q2(h) + 1−e−γ h 2 γ p2(h) + √ 2γ ∫ h h 2 1−e−γ(h−s) γ dB(s) p3 = e −γ h2 p2(h) + √ 2γ ∫ h h 2 e−γ(h−s)dB(s) q2(r) = q1 − ∫ r 0 α∇f(q2(s))ds+ √ 2α ∫ r 0 dW (s) (0 ≤ r ≤ h) p2(r) = p1 − ∫ r 0 ∇f(q2(s))ds q1 = q0 + 1−e−γ h 2 γ p0 + √ 2γ ∫ h 2 0 1−e−γ( h 2 −s) γ dB(s)
p1 = e −γ h2 p0 +
√ 2γ ∫ h 2
0 e−γ( h 2−s)dB(s)
Therefore, we have q̂(h) = q3, p̂(h) = p3 and q̂(h) = √ 2γ ∫ h h 2 1− e−γ(h−s) γ dB(s) + q1 − ∫ h 0 α∇f(q2(s))ds+ √ 2α ∫ h 0
dW (s)︸ ︷︷ ︸ q2(h)
+ 1− e−γ h2
γ
p1 − ∫ h 0
∇f(q2(s))ds︸ ︷︷ ︸ p2(h)  = √ 2γ
∫ h h 2 1− e−γ(h−s) γ dB(s)− ∫ h 0 α∇f(q2(s))ds+ √ 2α ∫ h 0 dW (s)− 1− e −γ h2 γ ∫ h 0 ∇f(q2(s))ds
+ q0 + 1− e−γ h2
γ p0 +
√ 2γ
∫ h 2
0
1− e−γ(h2−s)
γ dB(s)︸ ︷︷ ︸
q1
+ 1− e−γ h2
γ
e−γ h2 p0 +√2γ ∫ h 2
0 e−γ( h 2−s)dB(s)︸ ︷︷ ︸
p1  =q0 + 1− e−γh
γ p0 −
( α+
1− e−γ h2 γ )∫ h 0 ∇f(q2(s))ds
+ √ 2α ∫ h 0 dW (s) + √ 2γ ∫ h h 2 1− e−γ(h−s) γ dB(s) + √ 2γ ∫ h 2 0 1− e−γ(h2−s) γ dB(s)
+ 1− e−γ h2
γ
√ 2γ
∫ h 2
0
e−γ( h 2−s)dB(s)
p̂(h) =e−γ h 2 p1 − ∫ h 0
∇f(q2(s))ds︸ ︷︷ ︸ p2(h)
+ √ 2γ ∫ h h 2 e−γ(h−s)dB(s)
=e−γ h 2 e−γ h2 p0 +√2γ ∫ h 2
0 e−γ( h 2−s)dB(s)︸ ︷︷ ︸
p1
− e−γ h2 ∫ h
0
∇f(q2(s))ds+ √ 2γ ∫ h h 2 e−γ(h−s)dB(s)
=e−γhp0 − e−γ h 2 ∫ h 0 ∇f(q2(s))ds+ e−γ h 2 √ 2γ ∫ h 2 0 e−γ( h 2−s)dB(s) + √ 2γ ∫ h h 2 e−γ(h−s)dB(s)
It is clear that q̂(h), p̂(h) should be compared with the exact solution of HFHR at time h, which can be written as
q(h) =q0 + 1− e−γh
γ p0 − ∫ h 0 ( 1− e−γ(h−s) γ + α ) ∇f(q(s))ds+ √ 2α ∫ h 0 dW s + √ 2γ ∫ h 0 1− e−γ(h−s) γ dBs
p(h) =e−γhp0 − ∫ h
0
e−γ(h−s)∇f(q(s))ds+ √ 2γ ∫ h 0 e−γ(h−s)dB(s)
Subtracting q(h),p(h) from q̂(h), p̂(h) respectively, we obtain
q̂(h)− q(h) =− ( α+
1− e−γ h2 γ )∫ h 0 ∇f(q2(s))−∇f(q(s))ds
+ ∫ h 0 ( 1− e−γ(h−s) γ − 1− e −γ h2 γ ) ∇f(q(s))ds
p̂(h)− p(h) =− e−γ h2 ∫ h
0
∇f(q2(s))−∇f(q(s))ds+ ∫ h
0
( e−γ(h−s) − e−γ h2 ) ∇f(q(s))ds
It should be clear now that we will need to bound the term∇f(q2)−∇f(q) and ∇f(q). Since
q2(r) =q0 + 1− e−γ h2
γ p0 +
√ 2γ
∫ h 2
0
1− e−γ(h2−s)
γ dB(s)− α ∫ r 0 ∇f(q2(s))ds+ √ 2α ∫ r 0 dW (s)
q(r) =q0 + 1− e−γr
γ p0 − ∫ r 0 ( 1− e−γ(r−s) γ + α ) ∇f(q(s))ds+ √ 2α ∫ r 0 dW (s)
+ √ 2γ ∫ r 0 1− e−γ(r−s) γ dB(s),
we then have
q2(r)− q(r) = e−γr − e−γ h2
γ p0 − α ∫ r 0 ∇f(q2(s))−∇f(q(s))ds+ ∫ r 0 1− e−γ(r−s) γ ∇f(q(s))ds
+ √ 2γ
∫ h 2
0
1− e−γ(h2−s)
γ dB(s)−
√ 2γ ∫ r 0 1− e−γ(r−s) γ dB(s)
By Lemma D.3 and D.2, when 0 < h < 14L′ , we have the following for the solution of HFHR dynamics
E[ ∥∥x0,x0(h)− x0∥∥2] ≤ Ĉ0h
where Ĉ0 = 5.14 { (α+ γ)d+ h ( L′ )2‖x0‖2} and hence
E [∫ r
0 ∥∥∇f(q(s))∥∥2 ds] ≤E [2 ∫ r 0 ∥∥∇f(q(0))∥∥2 ds+ 2 ∫ r 0 ∥∥∇f(q(s))−∇f(q(0))∥∥2 ds] ≤E [ 2L2r
∥∥q(0)∥∥2 + 2L2 ∫ r 0 ∥∥q(s)− q(0)∥∥2 ds] ≤2L2r‖x0‖2 + 2L2E [∫ r 0
∥∥q(s)− q(0)∥∥2 ds] ≤2L2r‖x0‖2 + 2L2Ĉ0 ∫ r 0 sds
≤L2r ( 2‖x0‖2 + hĈ0 )
≤L2r ( 2.33‖x0‖2 + 5.14(α+ γ)dh )
(12)
Now E [ ‖q2 − q‖ 2 ] can be bounded as follow
E [∥∥q2(r)− q(r)∥∥2]
≤5  ( e−γr − e−γ h2 γ )2 ‖p0‖ 2 + α2E ∥∥∥∥∫ r 0 ∇f(q2(s))−∇f(q(s))ds ∥∥∥∥2 + E ∥∥∥∥∥ ∫ r 0 1− e−γ(r−s) γ ∇f(q(s))ds ∥∥∥∥∥ 2 
+ 5 2γE ∥∥∥∥∥ ∫ h 2
0
1− e−γ(h2−s)
γ dB(s) ∥∥∥∥∥ 2 + 2γE ∥∥∥∥∥ ∫ r
0
1− e−γ(r−s)
γ dB(s) ∥∥∥∥∥ 2  (Cauchy-Schwartz Inequality)
≤5 h24 ‖x0‖2 + α2L2r ∫ r
0
E ∥∥q2(s)− q(s)∥∥2 ds+ ∫ r
0
( 1− e−γ(r−s)
γ
)2 ds ∫ r 0 E ∥∥∇f(q(s))∥∥2 ds  + 5 { γdh3
12 +
2γd
3 r3
}
≤5
{ h2
4 ‖x0‖2 + α2L2r ∫ r 0 E ∥∥q2(s)− q(s)∥∥2 ds+ h33 E [∫ r 0 ∥∥∇f(q(s))∥∥2]+ 3γd 4 h3 }
≤5
{ h2
4 ‖x0‖2 +
3γd
4 h3 +
h3 3 L2 ( 2.33‖x0‖2 + 5.14(α+ γ)dh ) r + α2L2r ∫ r 0 E ∥∥q2(s)− q(s)∥∥2 ds }
≤5h2 { 1
4 ‖x0‖2 +
3γd
4 h+
h2 3 L2 ( 2.33‖x0‖2 + 5.14(α+ γ)dh )} + 5α2L2h ∫ r 0 E ∥∥q2(s)− q(s)∥∥2 ds
By Gronwall’s inequality and 0 < h ≤ 14L′ , we have
E [∥∥q2(r)− q(r)∥∥2] ≤5h2
{ 1
4 ‖x0‖2 +
3γd
4 h+
h2 3 L2 ( 2.33‖x0‖2 + 5.14(α+ γ)dh )} exp{5α2L2h2}
≤5h2 { 1
4 ‖x0‖2 +
3γd
4 h+
h2 3 L2 ( 2.33‖x0‖2 + 5.14(α+ γ)dh )} e 5 32
≤5.85h2 { 0.28‖x0‖2 + (0.06α+ 0.81γ)hd }
≤h2 { 1.64‖x0‖2 + (0.36α+ 4.74γ)hd } . (13)
With bounds in Equation (12) and (13), we are now ready to show p1 and p2. For p1, i.e. the order of the mathematical expectation of deviation, we have
∥∥∥∥∥∥E [[ q̂(h) p̂(h) ] − [ q(h) p(h) ]]∥∥∥∥∥∥ ≤ ∥∥∥E [q̂(h)− q(h)]∥∥∥+∥∥∥E [p̂(h)− p(h)]∥∥∥ ≤ ( α+
1− e−γ h2 γ )∥∥∥∥∥ ∫ h
0
E [ ∇f(q2(s))−∇f(q(s)) ] ds ∥∥∥∥∥+ ∥∥∥∥∥∥ ∫ h 0 ( 1− e−γ(h−s) γ − 1− e −γ h2 γ ) E [ ∇f(q(s)) ] ds ∥∥∥∥∥∥ + e−γ h 2 ∥∥∥∥∥ ∫ h 2
0
E [ ∇f(q2(s))−∇f(q(s)) ] ds ∥∥∥∥∥+ ∥∥∥∥∥ ∫ h
0
( e−γ(h−s) − e−γ h2 ) E [ ∇f(q(s)) ] ds ∥∥∥∥∥ ≤ ( α+ 1 + h
2
) L ∫ h 0 E ∥∥q2(s)− q(s)∥∥ ds
+ ∫ h 0 ∣∣∣∣∣1− e−γ(h−s)γ − 1− e−γ h 2 γ ∣∣∣∣∣+ ∣∣∣e−γ(h−s) − e−γ h2 ∣∣∣ ∥∥∥E [∇f(q(s))]∥∥∥ ds
≤L ( α+ 1 + h
2 )∫ h 0 E ∥∥q2(s)− q(s)∥∥ ds
+  ∫ h 0 ∣∣∣∣∣1− e−γ(h−s)γ − 1− e−γ h 2 γ ∣∣∣∣∣ 2 ds  12 +(∫ h 0 ∣∣∣e−γ(h−s) − e−γ h2 ∣∣∣2 ds) 12  (∫ h 0 ∥∥∥E [∇f(q(s))]∥∥∥2 ds) 12
≤L ( α+ 1 + h
2 )∫ h 0 ( E ∥∥q2(s)− q(s)∥∥2) 12 ds+ 1 + γ 2 √ 3 h 3 2 ( E ∫ h 0 ∥∥∥[∇f(q(s))]∥∥∥2 ds) 12 ≤L ( α+ 1 + h
2
) h2 { 1.64‖x0‖2 + (0.36α+ 4.74γ)hd } 1 2 + 1 + γ
2 √ 3 h2L
( 2.33‖x0‖2 + 5.14(α+ γ)dh ) 1 2
≤L (α+ 1.25)h2 ( 1.29‖x0‖+ √ 0.36α+ 4.74γ √ hd ) + (1 + γ)h2L ( 0.45‖x0‖+ √ 0.43α+ 0.43γ √ dh )
≤Lh2 max{α+ 1.25, γ + 1} ( 1.74‖x0‖+ (1.26 √ α+ 2.84 √ γ) √ hd )
The above derivation proves p1 = 2 with
Ĉ1 = Lmax{α+ 1.25, γ + 1} ( 1.74‖x0‖+ (1.26 √ α+ 2.84 √ γ) √ hd ) .
We now proceed with p2, i.e. mean-square error
E ∥∥∥∥∥ [ q̂(h) p̂(h) ] − [ q(h) p(h) ]∥∥∥∥∥ 2
≤2 ( α+ h
2
)2 E ∥∥∥∥∥ ∫ h
0
∇f(q2(s))−∇f(q(s))ds ∥∥∥∥∥ 2 + 2E ∥∥∥∥∥∥ ∫ h 0 ( 1− e−γ(h−s) γ − 1− e −γ h2 γ ) ∇f(q(s))ds ∥∥∥∥∥∥ 2
+ 2E ∥∥∥∥∥ ∫ h
0
∇f(q2(s))−∇f(q(s))ds ∥∥∥∥∥ 2 + 2E ∥∥∥∥∥ ∫ h
0
( e−γ(h−s) − e−γ h2 ) ∇f(q(s))ds ∥∥∥∥∥ 2
≤2 ( (α+ h
2 )2 + 1
) L2E (∫ h 0 |q2(s)− q(s)|ds )2 + 2 ∫ h 0 ∣∣∣∣∣1− e−γ(h−s)γ − 1− e−γ h 2 γ ∣∣∣∣∣ 2 ds ∫ h 0 E ∥∥∇f(q(s))∥∥2 ds
+ 2 ∫ h 0 ∣∣∣e−γ(h−s) − e−γ h2 ∣∣∣2 ds ∫ h 0 E ∥∥∇f(q(s))∥∥2 ds
≤2 ( (α+ h
2 )2 + 1
) L2h ∫ h 0 E|q2(s)− q(s)|2ds+ 1 + γ2 6 h3 ∫ h 0 E|∇f(q(s))|2ds
≤2 ( (α+ h
2 )2 + 1
) L2 { 1.64‖x0‖2 + (0.36α+ 4.74γ)hd } h4 + 1 + γ2 6 L2 { 2.33‖x0‖2 + 5.14(α+ γ)hd } h4
≤L2 max{(α+ 1.25)2, 1 + γ2} ( 3.67‖x0‖2 + (1.68α+ 10.34γ)hd ) h4
The above derivation implies p2 = 2 with Ĉ2 = Lmax{α+ 1.25, 1 + γ} ( 1.92‖x0‖+ (1.30 √ α+ 3.22 √ γ) √ hd ) .

heading: D.6 LOCAL ERROR BETWEEN ALGORITHM 1 AND THE EXACT STRANG’S SPLITTING METHOD
text: Lemma D.7. Assume f is L-smooth, 0 ∈ argminx∈Rd f(x), i.e. ∇f(0) = 0 and the operator ∇∆f grows at most linearly, i.e. ∥∥∇∆f(q)∥∥ ≤ G√1 +‖q‖2. If 0 < h ≤ 14L′ , 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 p1 = 2 and local mean-squared error of order p2 = 1.5, i.e. there exist constants C̄1, C̄2 > 0 such that∥∥Ex̂(h)− Ex̄(h)∥∥ ≤ C̄1hp1(
E [∥∥x̂(h)− x̄(h)∥∥2]) 12 ≤ C̄2hp2
where x̂(h) = [ q̂(h) p̂(h) ] is the solution of the exact Strang’s splitting method for HFHR with initial
value x0 = [ q0 p0 ] and x̄(h) = [ q̄(h) p̄(h) ] is the one-step result of Algorithm 1 with initial value
x0 = [ q0 p0 ] , p1 = 2 and p2 = 1.5. More concretely, we have
C̄1 = α(α+ 1.125)(L+G) [ 0.5 + 0.71‖x0‖+ (1.14 √ α+ 0.21 √ γh) √ hd ]
and C̄2 = L(α+ 0.73) ( 2.30 √ hαL‖x0‖+ (2.27 √ α+ 0.12 √ γh) √ d ) .
Proof. The solution of one-step exact Strang’s splitting integrator with step size h can be written as q3 = q2(h) + 1−e−γ h 2 γ p2(h) + √ 2γ ∫ h h 2 1−e−γ(h−s) γ dB(s) p3 = e −γ h2 p2(h) + √ 2γ ∫ h h 2 e−γ(h−s)dB(s) q2(r) = q1 − ∫ r 0 α∇f(q2(s))ds+ √ 2α ∫ r 0 dW (s) (0 ≤ r ≤ h) p2(r) = p1 − ∫ r 0 ∇f(q2(s))ds q1 = q0 + 1−e−γ h 2 γ p0 + √ 2γ ∫ h 2 0 1−e−γ( h 2 −s) γ dB(s)
p1 = e −γ h2 p0 +
√ 2γ ∫ h 2
0 e−γ( h 2−s)dB(s)
and the solution of one-step implementable Strang’s splitting integrator with step size h can be written as  q̄3 = q̄2(h) + 1−e−γ h 2 γ p̄2(h) + √ 2γ ∫ h 2 0 1−e−γ( h 2 −s) γ dB( h 2 + s) p̄3 = e −γ h2 p̄2(h) + √ 2γ ∫ h 2 0 e−γ( h 2−s)dB(h2 + s) q̄2(r) = q1 − ∫ r 0 α∇f(q1)ds+ √ 2α ∫ r 0 dW (s) (0 ≤ r ≤ h) p̄2(r) = p1 − ∫ r 0 ∇f(q1)ds q1 = q0 + 1−e−γ h 2 γ p0 + √ 2γ ∫ h 2 0 1−e−γ( h 2 −s) γ dB(s)
p1 = e −γ h2 p0 +
√ 2γ ∫ h 2
0 e−γ( h 2−s)dB(s)
Note that in the implementable Strang’s splitting method, φ flow can be explicitly integrated and hence q1,p1 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 q2(h) − q̄2(h) and p2(h)− p̄2(h). We have{
q2(h)− q̄2(h) = −α ∫ h 0 ∇f(q2(s))−∇f(q1)ds
p2(h)− p̄2(h) = − ∫ h 0 ∇f(q2(s))−∇f(q1)ds
(14)
Square both sides of the first equation in (14) and take expectation, we obtain
E ∥∥q2(h)− q̄2(h)∥∥2 =α2E ∥∥∥∥∥ ∫ h
0
∇f(q2(s))−∇f(q1)ds ∥∥∥∥∥ 2
≤α2E (∫ h 0 ∥∥∇f(q2(s))−∇f(q1)∥∥ ds )2
≤α2L2E (∫ h 0 ∥∥q2(s)− q1∥∥ ds )2
≤α2L2h ∫ h
0
E ∥∥q2(s)− q1∥∥2 ds
Note that q2 is the solution of a rescaled overdamped Langevin dynamics whose drift vector field is αL-Lipschitz, by conditional expectation version of Lemma D.2, for 0 < h < 14L′ < 1 4αL , we have E ∥∥q2(h)− q1∥∥2 ≤ C̄0h with C̄0 = 5.14{αd+ h(αL)2E‖q1‖2} and it follows that{
E ∥∥q2(h)− q̄2(h)∥∥2 ≤ α2L2C̄0h3
E ∥∥p2(h)− p̄2(h)∥∥2 ≤ L2C̄0h3.
Now consider p1, i.e., the deviation of mathematical expectation. By Ito’s lemma, we have q2(h)− q̄2(h)
=− α ∫ h
0
∇f(q2(s))−∇f(q1)ds
=− α ∫ h
0 [∫ s 0 −α∇2f(q2(r))∇f(q2(r))dr + α ∫ s 0 ∇∆f(q2(r))dr + ρ ] ds (15)
where ρ is a stochastic integral term. Take expectation and norm for Equation (15), we have∥∥∥E [q2(h)− q̄2(h)]∥∥∥ =α2 ∥∥∥∥∥ ∫ h
0
E [∫ s
0
∇2f(q2(r))∇f(q2(r))dr − ∫ s
0
∇∆f(q2(r))dr ] ds ∥∥∥∥∥ ≤α2
∫ h 0 E [∫ s 0 ‖∇2f(q2(r))‖2 ∥∥∇f(q2(r))∥∥ dr + ∫ s 0 ∥∥∇∆f(q2(r))∥∥ dr] ds ≤α2
∫ h 0 E [ L ∫ s 0 ∥∥q2(r)∥∥ dr + ∫ s 0 G(1 + ∥∥q2(r)∥∥)dr] ds
=α2(L+G) ∫ h 0 ∫ s 0 E ∥∥q2(r)∥∥ dr + α2Gh22
≤α2(L+G) ∫ h
0 ∫ s 0 E ∥∥q2(r)− q1∥∥+ E‖q1‖ dr + α2Gh22
≤α2(L+G) ∫ h
0 ∫ s 0 √ E ∥∥q2(r)− q1∥∥2 + E‖q1‖ dr + α2Gh22
≤α2(L+G) √ C̄0h h2
2 + α2(L+G)
h2
2 E‖q1‖+ α2G
h2
2 ≤α2 {√
C̄0h+ E‖q1‖ 2 (L+G) + G 2
} h2
≤1 2 α2(L+G)
{√ C̄0h+ E‖q1‖+ 1 } h2
Similarly, we have ∥∥∥E [p2(h)− p̄2(h)]∥∥∥ ≤ 12α(L+G){√C̄0h+ E‖q1‖+ 1}h2.
For p2, i.e., mean-square error, we have
E ∥∥q2(h)− q̄2(h)∥∥2 ≤α2E {∫ h 0 ∥∥∇f(q2(s))−∇f(q1)∥∥ ds }2
≤α2E {∫ h 0 1ds ∫ h 0 ∥∥∇f(q2(s))−∇f(q1)∥∥2 ds }
≤α2L2h ∫ h
0
E ∥∥q2(s)− q1∥∥2 ds
≤α 2L2C̄0
2 h3
Similarly we obtain E ∥∥p2(h)− p̄2(h)∥∥2 ≤ L2C̄02 h3. Recallq3 − q̄3 = q2(h)− q̄2(h) + 1−e −γ h 2
γ (p2(h)− p̄2(h)) p3 − p̄3 = e−γ h 2 (p2(h)− p̄2(h)) .
and it follows that when 0 < h ≤ 14L′ < 1∥∥∥∥∥E [ q3 − q̄3 p3 − p̄3 ]∥∥∥∥∥ ≤α(α+ 1 + h2 )(L+G) √ C̄0h+ E‖q1‖+ 1 2 h2 (16)
E ∥∥∥∥∥ [ q3 − q̄3 p3 − p̄3 ]∥∥∥∥∥ 2 ≤L2C̄0 ( α2 + 1 2 + h2 4 ) h3. (17)
Finally we need to bound E‖q1‖ 2 by E‖x0‖2, to this end, we have
E‖q1‖ 2 =E ∥∥∥∥∥q0 + 1− e−γ h 2 γ p0 + √ 2γ ∫ h 2
0
1− e−γ(h2−s)
γ dB(s) ∥∥∥∥∥ 2
≤(1 + h 2
4 )E‖q0‖ 2 + (1 +
h2
4 )E‖p0‖ 2 + 2γd
∫ h 2
0
( 1− e−γ(h2−s)
γ
)2 ds
≤(1 + h 2
4 )E‖x0‖2 +
γd 12 h3 (18)
=(1 + h2
4 )‖x0‖2 +
γd 12 h3 (19)
Collecting all pieces together, including (16), (17), (19), the definition of C̄0 and 0 < h < 14L′ , it is not difficult to obtain the following ∥∥∥∥∥E [ q3 − q̄3 p3 − p̄3
]∥∥∥∥∥ ≤C̄1h2E∥∥∥∥∥ [ q3 − q̄3 p3 − p̄3 ]∥∥∥∥∥ 2  12 ≤C̄2h 32
with
C̄1 = α(α+ 1.125)(L+G) [ 0.5 + 0.71‖x0‖+ (1.14 √ α+ 0.21 √ γh) √ hd ]
and
C̄2 = L(α+ 0.73) ( 2.30 √ hαL‖x0‖+ (2.27 √ α+ 0.12 √ γh) √ d )

heading: D.7 LOCAL ERROR BETWEEN ALGORITHM 1 AND HFHR DYNAMICS
text: Lemma D.8. Assume f is L-smooth, 0 ∈ argminx∈Rd f(x), i.e. ∇f(0) = 0 and the operator ∇∆f grows at most linearly, i.e. ∥∥∇∆f(q)∥∥ ≤ G√1 +‖q‖2. If 0 < h ≤ 14L′ , then compared with the HFHR dynamics, the implementable Strang’s splitting method has local weak error of order p1 = 2 and local mean-squared error of order p2 = 1.5, i.e. there exist constants C1, C2 > 0 such that ∥∥Ex(h)− Ex̄(h)∥∥ ≤ C1hp1(
E [∥∥x(h)− x̄(h)∥∥2]) 12 ≤ C2hp2
where x(h) = [ q(h) p(h) ] is the solution of HFHR with initial value x0 = [ q0 p0 ] and x̄(h) = [ q̄(h) p̄(h) ] is the solution of the implementable Strang’s splitting with initial value x0 = [ q0 p0 ] , p1 = 2 and p2 = 1.5. More concretely, we have
C1 = (L+G) max{α+1.25, γ+1} [ 0.5α+ (1.74 + 0.71α)‖x0‖+ ( 1.26 √ α+ 1.14α √ α+ 2.32 √ γ )√ hd ]
and C2 = Lmax{α+ 1.25, γ + 1} [ (1.92 + 2.30αL) √ h‖x0‖+ (2.60 √ α+ 3.34 √ γh) √ d ]
Proof. Denote by x̂(h) = [ q̂(h) p̂(h) ] the solution of the exact Strang’s splitting method with initial
value x0 = [ q0 p0 ] . By triangle inequality and Minkowski’s inequality, we have∥∥Ex(h)− Ex̄(h)∥∥ ≤∥∥Ex(h)− Ex̂(h)∥∥+∥∥Ex̂(h)− Ex̄(h)∥∥ ,(
E ∥∥x(h)− x̄(h)∥∥2) 12 ≤(E∥∥x(h)− x̂(h)∥∥2) 12 + (E∥∥x̂(h)− x̄(h)∥∥2) 12 .
By Lemma D.6 and D.7, we have∥∥Ex(h)− Ex̂(h)∥∥ ≤ Ĉ1h2, ∥∥Ex̂(h)− Ex̄(h)∥∥ ≤ C̄1h2( E ∥∥x(h)− x̂(h)∥∥2) 12 ≤ Ĉ2h 32 , (E∥∥x̂(h)− x̄(h)∥∥2) 12 ≤ C̄2h 32
and hence ∥∥Ex(h)− Ex̄(h)∥∥ ≤(Ĉ1 + C̄1)h2( E ∥∥x(h)− x̄(h)∥∥2) 12 ≤(Ĉ2 + C̄2)h 32
with
Ĉ1 + C̄1 ≤C1 ,(L+G) max{α+ 1.25, γ + 1} [ 0.5α+ (1.74 + 0.71α)‖x0‖+ ( 1.26 √ α+ 1.14α √ α+ 2.32 √ γ )√ hd ]
Ĉ2 + C̄2 ≤C2 , Lmax{α+ 1.25, γ + 1} [ (1.92 + 2.30αL) √ h‖x0‖+ (2.60 √ α+ 3.34 √ γh) √ d ]
E α DOES CREATE ACCELERATION EVEN AFTER DISCRETIZATION: AN ANALYTICAL DEMONSTRATION
If α → ∞ while γ remains fixed, then dq = −α∇f(q) + √
2αdW is the dominant part of the dynamics, and in this case the role of α could be intuitively understood as to simply rescale the time of gradient flow, which does not create any algorithmic advantage, as the timestep of discretization has to scale like 1/α in this case. However, finite α no longer corresponds to solely a time-scaling, but closely couples with the dynamics and creates acceleration. This is true even after the continuous dynamics is discretized by an algorithm .
We will analytically illustrate this point by considering quadratic f . In this case, the diffusion process remains Gaussian, and it 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 1 (1D, γ given; infinite acceleration). Consider f(x) = x2/2, γ fixed. The mean process is {
q̇ = p− αq ṗ = −q − γp
Consider, for simplicity, an Euler-Maruyama discretization of the HFHR dynamics, which coressponds to a Forward Euler discretization of the mean process (other numerical methods can be analyzed analogously):[
qk+1 pk+1
] = A [ qk pk ] , A = [ 1− αh h −h 1− γh ] .
We will show that, unless γ = 2, an appropriately chosen α will converge infinitely faster than the case with α = 0, if both cases use the optimal h.
To do so, let us compute A’s eigenvalues, which are
1
2
( 2− (α+ γ)h± h √ −4 + (α− γ)2 ) Consider the case where |α − γ| ≤ 2, then the eigenvalues are a pair of complex conjugates. Their modulus determines the speed of convergence, and it can be computed to be
1
2
√ (2− (α+ γ)h)2 + h2(4− (α− γ)2) = √ 1− (α+ γ)h+ (1 + αγ)h2
Minimizing the quadratic function gives the optimal h that ensures the fastest speed of convergence, and the optimal h is
h = α+ γ
2(1 + αγ)
and the optimal spectral radius is √ 1− (α+ γ) 2
4(1 + αγ) .
When one uses low-resolution ODE, in which α = 0, the optimal rate is 1 − γ2/4 (note it is not surprising that the critically damped case, i.e., γ = 2, will give the fastest convergence).
If γ 6= 2, the additional introduction of α can accelerate the convergence by reducing the spectral radius. For instance, if α = γ + 2, upon choosing the optimal h = 11+γ , the optimal spectral radius
is 0 (note in this case A actually has Jordan canonical form of [ 0 1 0 0 ] and thus the discretization
converges in 2 steps instead of 1, irrespective of the initial condition).
Demonstration 2 (multi-dim, γ, α and h all to be chosen; acceleration quantified in terms of condition number). Consider quadratic f with positive definite Hessian, whose eigenvalues are 1 = λ1 < · · · < λn = −1 for some 0 < 1. Assume without loss of generality that f = q21/2 +
−1q22/2. Similar to Demonstration 1, the forward Euler discretization of the mean process is
q1,k+1 p1,k+1 q2,k+1 p2,k+1
 = [A1 00 A2 ]
q1,k p1,k q2,k p2,k
 , A1 = [1− αh h−h 1− γh ] , A2 = [ 1− α −1h h − −1h 1− γh
] (20)
We will (i) find h and γ that lead to fastest convergence of the ULD discretization, i.e. the above iteration with α = 0, and then (ii) constructively show the existence of h, γ and α that lead to faster convergence than the optimal one in (i) — note these may not even be the optimal choices for HFHR, but they already lead to significant acceleration. More specifically,
(i) In a ULD setup, α = 0. It can be computed that the eigenvalues of A1 and A2 are respectively
1
2
( 2− hγ ± h √ −4 + γ2 ) and
1
2
( 2− hγ ± h √ −4 −1 + γ2 ) We now seek γ > 0, h > 0 to minimize the maximum of their norms for obtaining the optimal convergence rate. This is done in cases.
Case (i1) When γ ≤ 2, both A1 and A2 eigenvalues are complex conjugate pairs. To minimize the maximum of their norms, let’s first see if their norms could be made equal.
A1 eigenvalue’s norm squared ×4 is
(2− hγ)2 − h2(−4 + γ2) = 4(h− γ/2)2 + 4− γ2 (21)
A2 eigenvalue’s norm squared ×4 is
(2− hγ)2 − h2(−4 −1 + γ2) = 4 −1(h− γ/2)2 + 4− γ2 (22)
It can be seen that for (21) is always strictly smaller than (22) for any h > 0. Therefore, the max of the two is minimized when h = γ/2, and the corresponding max value is 4− γ2. γ that minimizes this max value is γ = 2. Corresponding rate of convergence is
√ 1− .
Case (i2) When γ ≥ 2 −1/2, both A1 and A2 eigenvalues are real. Since 1, we can order them×2 as
2−hγ−h √ −4 + γ2 < 2−hγ−h √ −4 −1 + γ2 < 2−hγ+h √ −4 −1 + γ2 < 2−hγ+h √ −4 + γ2 < 2.
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− hγ − h √ −4 + γ2) = 2− hγ + h √ −4 + γ2.
This gives h = 2/γ (which does verify the assumption that the smallest of four is negative). Corresponding max of their norms is thus √ 1− 4/γ2. γ that minimizes this max value is γ = 2 −1/2, which gives rate of convergence of √ 1− .
Case (i3) When 2 ≤ γ ≤ 2 −1/2, A1 eigenvalues are real and A2 eigenvalues are complex conjugates. Again, the max of their norms is minimized if the norms can be made all equal.
Note A1 eigenvalues cannot be of the same sign, because otherwise 2 − hγ − h √ −4 + γ2 =
2−hγ+h √ −4 + γ2, which means either h = 0 or γ = 2, but if γ = 2 then 2−hγ+h √ −4 + γ2
being equal to 2*norm of A2 eigenvalue, which is √
4 −1(h− γ/2)2 + 4− γ2, leads to h = 0 again.
Therefore, the equality of norms of A1, A2 eigenvalues means −(2− hγ − h √ −4 + γ2) = 2− hγ + h √ −4 + γ2 = √ 4 −1(h− γ/2)2 + 4− γ2.
The first equality gives hγ = 2, which, together with the second equality, gives h = ± √
2 1+ . Selecting the positive value of optimal h, we also obtain optimal γ = √
2(1 + ) −1/2, which is ≤ 2 −1/2 and thus satisfying our assumption (2 ≤ γ ≤ 2 −1/2). The corresponding rate of convergence is thus
1
2
( 2− hγ + h √ −4 + γ2 ) = √ 1− 1 + .
Summary of (i) Since √
1− 1+ <
√ 1− , the ULD Euler-Maruyama discretization converges the
fastest when
h =
√ 2
1 + , γ =
√ 2(1 + ) −1/2,
and the corresponding discount factor of convergence (i.e. base of exponential convergence) is√ 1− 1 + , where = 1/κ with κ being Hessian’s condition number. (23)
(ii) Now consider the HFHR setup. Let’s first state a result: when
γ =
√ 4c2 4 + 8c2 3 + 4c2 2 + 2 − 2 + 1 + + 3
2c 2 + 2c > 0, (24)
α = − √
4c2 4 + 8c2 3 + 4c2 2 + 2 − 2 + 1 + 3 + 1 2c 2 + 2c > 0, h = c (25)
for any c > 0 independent of , the iteration (20) converges with discount factor
1√ 2(1 + )
√ (1− ) ( 1− + √ 4c2 4 + 8c2 3 + (4c2 + 1) 2 − 2 + 1 ) . (26)
While the exact expression is lengthy, it can proved that the HFHR non-optimal discount factor (26) is strictly smaller than the ULD optimal discount factor (23) for not only small but also large ’s.
For some quantitative intuition, discount factors respectively have the following Taylor expansions in :
HFHR non-optimal: 1− 2 +
( c2
2 + 2
) 2 +O ( 3 )
(27)
ULD optimal: 1− + 2
2 +O
( 3 )
(28)
The exact expressions of discount factors are also plotted in Fig.4 (c = 1 was arbitrarily chosen) and one can see acceleration for any (not necessarily small) .
(ii details) How were values in (25) chosen? Following the idea detailed in (i), we consider a case where A1 eigenvalues are both real, A2 eigenvalues are complex conjugates, and all 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 at least the following equations{ trA1 = 0 detA1 + detA2 = 0 (29)
One can solve this system of equations to obtain α and γ as functions of h. Following the idea of choosing h small enough to resolve the stiffness of the ODE{
q̇2 = p2 − α −1q2 ṗ2 = − −1q2 − γp2 ,
pick h = c . Then (29) gives
γ =
√ 4c2 4 + 8c2 3 + 4c2 2 + 2 − 2 + 1 + + 3
2c 2 + 2c
α = − √
4c2 4 + 8c2 3 + 4c2 2 + 2 − 2 + 1 + 3 + 1 2c 2 + 2c
or
γ = − √
4c2 4 + 8c2 3 + 4c2 2 + 2 − 2 + 1 + + 3 2c 2 + 2c
α =
√ 4c2 4 + 8c2 3 + 4c2 2 + 2 − 2 + 1 + 3 + 1
2c 2 + 2c
The former is our choice (25) because it can be checked that the latter leads to detA1 > 0 which violates the assumption of a pair of plus and minus real eigenvalues.
It is possible to find optimal α, γ, h for HFHR for the Gaussian cases. One has to minimize detA2 under the constraint detA2 > 0 in addition to (29). And then do similar calculations for the other 3 cases, and then finally the best among the 4 cases. Doing so however does not give enough insights to determine optimal hyperparameters for sampling general distributions.

heading: F RANDOMIZED MIDPOINT DISCRETIZATION OF HFHR

heading: F.1 THE ALGORITHM
text: HFHR is based on a continuous dynamics that adds HFHR corrections to the Underdamped Langevin Dynamics (ULD). It can be turned into a sampling algorithm via either a low-order time discretization (e.g., HFHR Algorithm 1) or a more accurate one. To complement the main text, this section demonstrates the latter, based on a powerful recent progress in discretizing ULD, known as 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 betterO(d 13 ) dimension dependence of mixing time than first-order discretization of ULD, e.g., 1st-order KLMC (Dalalyan & Riou-Durand, 2020). Although RMA is originally designed specifically for ULD only, it is a general idea and already adapted to overdamped Langevin (He et al., 2020). Here we show it can be easily adapted to HFHR as well, as illustrated by the following Algorithm 2. Red highlights algorithmic changes we made to account for the HFHR corrections of ULD.
Algorithm 2 Randomized Midpoint Algorithm from Shen & Lee (2019), adapted for HFHR 1: Input: potential function f and its gradient ∇f , damping coefficients α and γ, step size h,
initial condition (q0,p0) 2: procedure RMA-HFHR(f,∇f, α, γ, h, q0,p0) 3: k = 0 and initialize [ q0 p0 ] 4: while not converged do 5: Generate an independent uniform random variable θk ∼ U(0, 1) 6: Generate Gaussian random vectors ( W 1k+1,W 2 k+1,W 3 k+1 ) ∈ R3d as in (Shen & Lee,
2019, Appendix A) 7: Generate Gaussian random vectorsB1k+1,B 2 k+1 ∈ Rd as described by (31)
8: qk+ 12 = qk + 1 γ (1 − e −γθkh)pk − 1γ ( θkh− 1γ (1− e −γθkh) ) ∇f(qk) + W 1 k+1
−αθkh∇f(qk) + √
2αB1k+1 9: qk+1 = qk + 1 γ (1 − e −γh)pk − 1γh(1 − e −γ(h−θkh))∇f(qk+ 12 ) + W 2 k+1 −
αh∇f(qk+ 12 ) + √ 2α(B1k+1 +B 2 k+1)
10: pk+1 = pke −γh − he−γ(h−θkh)∇f(qk+ 12 ) + 2W 3 k+1 11: k ← k + 1 12: end while 13: end procedure
The red parts basically correspond to two Euler-Maruyama time-steppings of an auxiliary dynamics that contains only the HFHR correction terms
dq = −α∇f(q)dt+ √ 2αdBt, (30)
first over a θkh timestep, and then over an h timestep. These two steps originate from an operator splitting treatment of the full HFHR dynamics (eq.6), which is split into ULD and (30). Therefore, it is natural to see that
B1k+1 = ∫ h(k+θk) hk dBt, B 2 k+1 = ∫ h(k+1) h(k+θk) dBt,
and therefore B1k+1 and B 2 k+1 are, when conditioned on θk, centered Gaussian vectors independent from each other and theW ’s, each being d-dimensional with i.i.d. entries, and they can be generated via
B1k+1 = √ θkhξ 1 k+1, B 2 k+1 = √ h− θkhξ2k+1, (31)
where ξ1k+1 and ξ 2 k+1 are i.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 α. However, since we have already used this letter for the HFHR correction coefficient, we use instead θ to denote this uniform random variable.
Remark F.2. From the red text, it is easy to see that if α = 0, Algorithm 2 degenerates to RMA for ULD. Nevertheless, Algorithm 2 is again just one RMA discretization of HFHR but not the only one.

heading: F.2 NUMERICAL RESULTS: HFHR AGAIN ACCELERATES
text: To numerically compare the RMA discretization of HFHR dynamics and ULD dynamics (note we don’t compare 1st-order HFHR Algorithm 1 with RMA-ULD as we’d like to compare apple with apple), we conduct an experiment very similar to that in Sec.6.1, with the same nonlinear potential function. We run both RMA for ULD and RMA for HFHR with dimension d = 10, initial value (100 × 1d,0d), h = 1 (chosen to be near the stability limit of RMA-ULD), a family of γ ∈ {0.1, 0.2, 0.5, 1, 2, 5, 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, 5, 10, 20, 50, 100}. For each algorithm and each set of parameter values, we run 1,000 independent realizations to compute statistics and estimate the mean time of reaching ε = 0.1 neighborhood of the target distribution. Then, for each α (including α = 0, which is the original RMA), we optimize over γ choices to get the best results. To further reduce variance, we also repeat the experiment with 100 different random seeds.
Too large α values with which Algorithm 2 fails to reach -neighborhood are not plotted and the final results are shown in Figure 5. It clearly suggests that with appropriated chosen α (α = 0.5 in our case), RMA discretized HFHR dynamics requires fewer iterations than RMA discretized ULD, which suggests a better iteration complexity.

Review summaries:
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). 

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.

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.

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.