Double Randomized Underdamped Langevin with Dimension-Independent Convergence Guarantee
Authors: Yuanshi Liu
Cong Fang
Tong Zhang
Venue: NeurIPS 2023 poster (2023)
Keywords: Langevin
Dimension dependence
Acceleration
TL;DR: We propose the double-randomized underdamped Langevin algorithm that achieves a dimension-independent convergence speed.

Abstract: This paper focuses on the high-dimensional sampling of log-concave distributions with composite structures: $p^*(\mathrm{d}x)\propto \exp(-g(x)-f(x))\mathrm{d}x$. We develop a double randomization technique, which leads to a fast underdamped Langevin algorithm with a dimension-independent convergence guarantee. We prove that the algorithm enjoys an overall $\tilde{\mathcal{O}}\left(\frac{\left(\mathrm{tr}(H)\right)^{1/3}}{\epsilon^{2/3}}\right)$ iteration complexity to reach an $\epsilon$-tolerated sample whose distribution $p$ admits $W_2(p,p^*)\leq \epsilon$.  Here,  $H$ is an upper bound of the Hessian matrices for $f$ and does not explicitly depend on dimension $d$. For the posterior sampling over linear models with normalized data, we show a clear superiority of convergence rate which is dimension-free and outperforms the previous best-known results by a $d^{1/3}$ factor. The analysis to achieve a faster convergence rate brings new insights into high-dimensional sampling.

Double Randomized Underdamped Langevin with Dimension-Inddpeedent Convergence Guarantee
Yuanshi Liu", Cong Fang"*, Tong Zhang?
School of Intelligence Science and Technology, Peking University
Department of Computer Science & Engineering, the Hong Kong University of Science and Technology
liuu_yuashh, fangcong @ pku. edu.ee tongzhang @ woonzzanng-ml.org
Abstract
This paper focuses on the high-dimensional sampling of log-concave distributions with composite structures: p" (dxx exp(-g(x) f(x))xx. We develop a double randomization technique, which leads to a fast underdamped Langevin algorithm with a dimension iaarsiioonnddeeddeent convergence guarantee. We prove that
the algorithm enjoys overall Ö (ur(H))* iteration complexity reach <2/3
c-tolerated sample whose distribution p admits W2(P. p²) S Here, H is upper bound of the Hessian matrices for f and does not explicitly depend on dimension d. For the posterior sampling over linear models with normalized data, show a clear superiority of convergence rate which is dimension frre and outperforms the previous best- results by a d/33 factor. The analysis to achieve a faster convergence rate brings new insights into high-dimensional sampling.
Introduction
Sampling from a ligh-dimenssonal distribution serves as one of the key components in statistics, machine learning, and scientific computing, and constitutes the foundation of the fields including
Bayesian statistics and generative models (Liu and Liu 2001 Brooks et al. 2011 Song et al.
2020 Recently, there an emerging trend in designing provably faster Markov Chain Monte Carlo
(MCMC) algorithms using techniques from first-order optimization [Dalalyan 2017
Durmus et al.
2019 Cheng and Bartlett 2018 Vempala and Wibisono 2019 Chewi et al. 2021 One typical
MCMC algorithm that allows following the idea is the Langevin-type algorithms, which are of the central interest of this paper.
Langevin-type algorithms originated in statistical physics discretize a stochastic differential equation with stationary distribution corresponding to the target. For Gibbs distribution p(x) xx eU((x) the standard overdamped Langevin algorithm uses Fuler Maruyama scheme discretize the diffusion process dx, (x,(x,) + V2dB, The algorithm iteratively performs the updates as for n = 0,1,
X X, (x,( Vhen.
(1.1)
where en (O, follows the normal distribution and h 0 the step size.
This paper focuses on the dimension dependence of the convergence behavior of Langevin-type algorithms. Specifically, consider sampling problems over potentially high-dimensional and strongly log-concave distributions with composite structure:
e e-U(x) e -((x))((x) where
g(x) = 2 |x| Regarding p as posterior distribution, g=((x) corresponds to Gaussian prior and e f(x) corresponds to the likelihood. This structure also includes general m-strongly convex
37th Conference on Neural Information Processing Systems (NcurIPS 2023).

---

potential functions (xx), which can be split into a strongly convex term g(x) convex one f(x) (x) z ||
2 ||x|| and a weakly
The analysis of Langevin sampling from an optimization viewpoint may date back to Jordan et al. [1998]. 'he viewpoint has led to a surge of works that establish quantitative convergence guarantees for overdamped Langevin algorithms
sample log-concave distributions [Dalalyan 2017 Durmus and Moulines 2016]. It should be noted that these convergence guarantees always involve an extra d dimension dependence due to the injection of non-negligible Gaussian noise, whereas the convergence rates of first-order optimization are often dimension free [Nesterov] 2003]. Various accelerated methods have been proposed that can mitigate the dimension dependence and also achieve faster convergence. First, as the momentum acceleration of Langevin algorithms (Ma et al. 2021]. underdamped Langevin Monte Carlo has a more stable trajectory and is known to exhibit faster convergence rate [Cheng et al.] 2018 Ma et al., [02TT] Zhang et al.] [203]] Apart from using different trajectory, more stable discretization schemes of the same diffusion process also lead to better convergence (Shen and Lee, 2019. Wibisono 2019, He et al.) 2020 Li et al. 2019]. However, the dimension dependence remains. To the best of our knowledge, the fastest randomized midpoint
4//8
method for underdamped Langevin under Wasserstein distance achieves a O
convergence
rate and still has a d=/3 dimension dependence [Shen and Lee 2019].
Recently another thread of works studies the convergence with weaker dimension dependence. Some researchers explore further general assumptions on the target distribution. The idea by observing that some smoothness conditions can average out the dimension-dependent errors brought by the noise using Ito's formula. For example, Li et al. [2021) achieves a v d dimension dependence with a 3-rd order growth condition. Another approach investigates the curvature of the target distribution and illustrates that dimension often does not determine the complexity of the sampling problem. Vono jet al. [022] propose an ADMM-type splitting algorithm with a dimension-free convergence rate when the likelihood is separable. Along the same thread, Freund et al. [202] study the convergence rate of Langevin algorithms with no explicit dependence on dimension and propose to characterize the convergence rate by the upper bound of Hessian matrices of f. Specifically, by letting H be the upper bound of the Hessian matrices of f. they show a variant of the overdamped Langevin algorithm
achieves convergence rate of O ar(M1) in KI. divergence. This result improves the rate under wide
conditions because many high-ding sigh-dimensional sampling problems are intrinsically -owdiiensional in the sense that tr tr(H)= o(d), which frequently appears in machine learning. For example, when the potential function has ridge separable structure with mild conditions, tr( H) can be dimension-free (see Section 3.3 for more details).
Though these works succeed in obtaining a convergence rate of the Langevin algorithm with a weak dependence on the dimension, some important questions remain:
How to design provably faster algorithms with weak dependencies on dimension even for the log-concave sampling?"
Such a question is significant for understanding high-dimensional sampling and already includes lots of applications in real practice.
In this paper, we follow the regime of Freund et al. [002 to design provably faster Langevin algorithms, which answers the above question affirmatively. We propose a double-randomized algorithm showing that variants of underdamped algorithms inherit the lower dimensional dependence in overdamped Langevin and a special design of random stepsize can still keep the property when using more stable discretization.
Specifically, we develop a double randomization technique and obtain the Guable-Randomized Underdamped Langevin (DRUL) algorithm. We consider an averaged contraction effect to avoid dimension dependence. We design two distributions to achieve the mid-point acceleration with a
randomized stepsize. DRUL is proven to enjoy an overall Ö (tr(#))))/ iteration and gradient
ca/a
complexity. For the posterior sampling over linear models, we show a clear superiority of DRUL,
which achieves a dimension-free Ô I convergence rate. The novel perspective of DRUL introducing the random step size at each update. Such a random step size combined with a random midpoint point proposed by [Shen and Lee [2019) reduces the discretization error and provides new insights into designing provably faster sampling algorithms.

---

Table 1: Comparison of the convergence rates for most related works. We consider the convergence in Wasserstein distance to a strongly log-concave and log-Lipschitz smooth target distribution (we summarize below the methods if additional assumptions required). EU stands for Euler- MaallerMaruyamaa discretization and RMM stands for Randomized Midpoint discretization. 'Ridge-Separable Case" refers to the convergence rate when f admits ridge-separable structure (see (3.2) ). Note that some
results are established the KI. divergence. We convert these convergence rates into the ones in Wasserstein distance using Talagrand's inequality.
Method
Convergence Rate
Ridge-Separable Case
Overdamped as Composite OptimizCRation
Durmus et al. 2019
O (
(Proved in KL divergence)
Overdamped with EU
[Li et al. 2021]
(With an additional Imear growth condition)
Underdamped with RMM
d=/3 2273
4=/3 2473
[henn and Lee 2019
Overdamped as Composite Optimization
Frreund et al. 2022
o
[((H+| 11
0
(Proved in KL divergence)
DRUL (Ours)
[Ix, 113)1/3
1 2373
c473
In summary, the contributions of the paper are listed below:
(A) We propose the Donble-Randomized technique and design the DRUL algorithm.
(B) We show the DRUL converges to the target distribution in Wasserstein distance in
+|x-)) iterations. For posterior sampling over generalized linear models, 24/3
dimension-fre Ö 27/3
complexity can be achieved.
2 Related works
There has been surge of works investigating the asymptotic guarantees for the Langevin-type algorithms [Roberts and Tweedie| 1996] Mattingly et al.) 2002). And a series of recent works establish the non-asymptotic quantitative analysis framework of the Langevin-type algorithms. For performance on strongly log- log-concave distributions, early works [Dalalyan 2017 Durmus and Moulines] (2016) establish the non-asymptotic convergence rate in TV distance with Lipschitz gradients assumption for overdamped dynamics. Along the same setting, similar convergence rates are achieved in KL divergence using the Wasserstein gradient flow [Cheng and Bartlett 2018, Durmus al. 2019 or Wasserstein distance using the contractive property for overdamped Langevin. Underdamped Langevin accelerates the vanilla Langevin algorithms (Ma et al. [202T). With the same setting mentioned above, faster convergence rate can be established [Cheng al. 2018, Dalalyan and Riou-Durand, 2020 Shen and Lee) 2019 Zhang et al. 2023) for underdamped Langevin algorithms. setting, other examples les also show better results can be achieved, such as Durmus and Moulines [2019), Li et al. [2019) for overdamped Langevin algorithms. Previous works also investigate the convergence rate without strongly log-convex assumptions, such as the log-Sobolev inequality (LSI) condition which is similar to the Polyak-Lojasiewir condition in optimization. With target distributions satisfying LSI and log-Lipschitz-smoole the results can be extended for both overdamped Langevin /Vempala and Wibisono 2019 and underdamped Langevin
Beyond the log Lipschitz-smoot
[Ma et al. 2021 Zhang et al. 2023].
Another thread of works explores the dimension dependence of Langevin-type algorithms. Despite the great similarity between the analysis in Langevin algorithms and optimization, the convergence of Langevin-type algorithms depends on dimension the log-coneave setting whereas convex optimization algorithms can often achieve dimennsion-nndeeendeee results. With general strongly log-
3

---

concave (or LSI) and log-Ling Leg-LLpsshizzsmooole condition, previous works establish O (1² convergence Durmus et al. 2019] and O Ld
convergence rate for under- damped algorithms [Cheng et al. 2018 to ensure finding a solution 2n, with (x(Law(xxn)))) <
rate for overdamped algorithms
or K L (xn), S
With an additional linear growth condition on third-order derivative, Li
et al. 2021 achieve O Vd
convergence rate for overdamped Langevin. And for underdamped
Langevin algorithms, Shen and Lee [2019) improve the dependence of dimension and obtains a
convergence rate. Recent work of Freund al. [2022] characterizes the dimensional dependence of convergence rate by the upper bound of the Hessian matrix of f. As discussed in Section Freund et al. [2022] establish a O (t(())xxx111" convergence rate in KL divergence,
O
da/3
which implies a O br(
convergence rate in Wasserstein distance. And when f admits a segae-sep formula, Freund et al. [2022] obtain dimension-free O (2²) convergence rate with some mild assumptions. We summarize the comparison of these most related works Table [
Preliminary and problem setup
3.1 Notations
We use the convention O (.) and Q (-) to denote lower and upper bounds with a universal constant. (-) ignores the polylogarithmic dependence. And use f S to denote f Ogg). Use W2(H, v) to denote the 2- W-ssseerttein distance of distribution μ and v. Use Law(X) to denote the distribution of the random variable X. The Frobenius norm is denoted by ||, while || ||2 stands for operator 2-norm for matrices.
3.2 Sampling problem
We consider the distributions with a composite structure:
dp(x) x exp(-U(x))dxx = exp(-g(x) fxx) dx,
(3.1)
where g(x) 3 ||| is a quadratic function. In the context of posterior sampling, the associated
task is sampling from a distribution with a Gaussian prior. 'he composite structure also includes the general m-strongly convex function, which be divided into g(x) 3 |||| and the weakly
convex function fxx Ux) z [|x]] We make the following assumption on
Assumption 3.1. f E C2 is convex and has I.-Lipschitz continuous gradients, i.e. 0 V LI.
corresponds to making m-strongly convex and I. + m-LLischitz smooth assumptions on the potential function U, which is basic setting and widely studied in the Langevin sampling literature (see e.g. Dalalyan 2017] Cheng and Bartlett [2018). CCeng et al. (2018) Shhn and Lee (2019). )
In addition to the previous assumption, we follow Freund et al. [2022 to characterize the convergence rate by a new factor H to avoid explicit dimension dependence.
Definition 3.2. Let H be upper bound of the Hessian matrices of i.e. H Z fx)
Note that the Lipschitz smooth condition in Assumption 3.1 provides a loose bound for H since we always have H < LI. This implies that tr(H) will reach dL in the worst case, whereas this quantity can be much smaller than dL under wide conditions. One typical example when f admits so-called ridge separable structure shown below.
3.3 Example: ridge separable functions
Definition 3.3. f is said admit the ridge separable form if
((x) Σ or(aaf(
(3.2)
where o, are all univariate functions, and a, given vectors in R4

---

Ridge separable functions contain many applications in machine learning, such as regression or classification over generalized linear models well (deep) neural networks in the neural tangent kernel regime [Gelman and Hill, 2006 Jacot et al. 2018). We follow the argument of Freund et al. [2022]. showing realizable conditions that ensure tr(H) to be dimension-free. Assumption 3.4. The function 0; C² has bounded second derivative, i.e. of Lo for all i [n].
Assumption 3.5. For all i [n]. then norm of is bounded by R, i.e. ||a:|| 2 S R2.
When Assumptions
3.4 and 3.5 hold, the Hessian has the upper bound <²ʃ(x) =
E n (a? a)
La Σ a.a? Let
La -i= a.a and
tr(H) = tr Σ
Lo a.a
Lo Σ |a:l & Lolt",
thus illustrates that (r( has a dimension-free R² upper bound. Meanwhile the Lipschitz constant of V can be bounded by
Lo n Σ
(i) & ||
a T
LoRt"
It indicates that worst-case upper bounds of ir( (r(H) and Lip(V are the same.
Note that Assumptions 3.4 and 3.5 are easy to achieve practice. For example, consider using posterior sampling to compute the Bayes estimator over linear model, whose advantages against maximum posterior estimator have been discussed broadly (see e.g. [Audibert] [2009]. Here, is associated with the data. So Assumption 3.5 can be realized by simply normalizing the data. Moreover, o; corresponds to the loss function and is only required to have a bounded second derivative by Assumption 3.4 Note that sampling from a ridge separable potential functions extensively studied in related literature (see [Mou et al. 2021 Vono et al. 2022
Lee et al. 2018) as examples).
3.4 Preliminary
Overdamped Langevin. The overdamped Langevin dynamics for target distribution 3.1 is diffusion process that evolves along the following SDE
dx(()= - V ()(((()) + V2dB,
(3.3)
where B, is the standard Brownian motion. The overdamped Langevin algorithms simulate and discretize the SDE (.3). Different algorithms vary mainly by different discretization methods.
Underdamped Langevin Underdamped Langevin algorithms accelerate the convergence rate of overdamped Langevin algorithms and instead discretize the following dynamics
dx(t) v(t)dt, dv(t) (t) -u-V -tt
(3.4)
(((()))dttuVoaaa (V)(( (v(d)dt+++++++ 2 /udB
The diffusion process has the stationary distribution P(x.v) X exp Zu || Ux)).
Contractive Property One notable property that yields the convergence of Langevin dynamics is the contractive property [Cheng and Bartlett 2018 Cheng et al. 2018 given as follows.
Definition 3.6. A stochastic differential equation has contractive property if there exists positive constant m and
Bl|x(()
(0)) = exp(-m(),
(3.5)
for any pair of solutions x(t) and yt) driven by the same Brownian motion.
Remark 3.7. ((t)) - (+) is an upper bound of squared Wasserstein distance of Law (aw(x(t)) and Law(y()). and thus 3.5) implies a geometric convergence of the Wasserstein distance.
Contractive property can be established for both underdamped and overdamped I Langevin dynamics under suitable conditions. The convergence analysis for our proposed algorithms follows a similar argument the contractive property.
5

---

Algorithm Double-Randomized Underdamped Langevin (DRUL)
Require: Iteration N, target function U' - ||| + f(x), initial point (Xo, vo). max step size
I and K Ltm
u m+L p(dl) ( hx dt with support [0, At]-
set (d()) x t-h
1. dt with support [0, h).
for n - 1,2,
do
Sample Ot p' and Bn ~.p p. Obtain the covariance matrix E(xm: 3r).
Sample random vector (G, H, W) from distribution (V,, E).
8, + (1 + (12(Onnv. Vn + u f A12(s - an (((x,) s + 2/uH.
Update Xn+1 (n, + A12(Bn)v, + S0" A12(s B,,) (x,, + 2 Update Vn+1 A21 (1,)), + A22(Bn) + +u u A22(s V (2,), +
end for
4 Double-randomized underdamped Langevin algorithm
In this section, we introduce our double-randomized sampling method and present our main result. We will illustrate our intuition in Section
The proposed algorithm is built upon the following discretization scheme given a fixed point Tn
dx, (t) vv, (t)dt, dt,, = -Vg - VV (v + 2/ddB,.
(4.1)
The purpose of splitting the strongly convex part of (xx is to avoid the mn dimension dependence of Hr) Although m is reasonably small in practice, one cannot obtain fully dimensional-free convergence rate if the term with respect to is discretized.
Denote the solution of process (4.1) at time t given starting point (x,,, Vn), Brownian motion [B. step size I and point In by [B, Xn; (Brlosego, (Xn, vn)). The double-randomizee algorithm performs the following one-step update:
Xk+1 = [[,, Xni (BrJosesa: (Xn, vn)); (BrJosess, (X.s vn))..
xn;
(4.2)
The algorithm introduces random step size 3 p deterministic h. And let a follows the distribution p' x (e ch
defined on [0, A] other than
Tox
on 00, h).
The analysis focuses on a Gaussian prior setting. Under a Gaussian prior, the following Lemma points out that J is linear in starting point (X,, Vn,) and has decoupled Brownian motion, and thus the update (4.2 can be easily implemented.
Lemma 4.1. Assume g(x) = m ||x If xx, follows the discretized Langevin diffusion process 4.1 with starting point (Xn. Vn) and a given Xn, for any t >> 0, it satisfies the integral equation
x, (t) Jt; (Brlssse vn)) An (t),, A12(t)v..
+ u L' A1z(s 1 (x,) 2/ū 1' A12(s - -tdB
(4.3)
where
VI um
VI um
An1(1)
6 um
(-1+vI-umj
2/1 11
wm
2/1 um
A12(t) 2
1
um
+
1
(-1+/I-um))
um
2 umni

---

are deterministic functions of t. Moreover, if x,, (t) follows the exact Langevin process (3.4) with starting point (X,, Vn), then for any
x (n1())n A1z(O)n +u (xV())))V(((((()))
(4.4)
+ +2// / -
Given the integral formula (4.3), the algorithm can be summarized as Algorithm ! In Algorithm I A21, A22 deterministic sealar functions, and Σ R2² - R3dx3d is also deterministic, For explicit formula, please refer to Appendix A The distribution of (H, G, W) is induced by the coupling between (BrJostga and (B:Josigg
(.2
4.1 Main theorem
The convergence guarantee of DRUL can be stated as the following theorem.
Theorem 4.2 Main theorem, convergence DRUL). For any tolerance e E (0, 1), denote the minimizer of (x) by and set the step size
2/3
S min
12C2k'
(24Crmi m[L+m] tr (f) + ||x
where C2 2 is universal constant.
With initial point (Xo. vul. define Zo (x.(.v (11x0- - x|2 + ||xo + Vo x1|²).
Then under Assumptions 3.1 when
Sek
2EL
log
h
2
Algorithm outputs x,, such that W2 WLaw(x(X, (x,): p) S
Theorem 4.2 shows an overall δ (lx 11)
iteration and gradient complexity to find an
e-approximate sample in 2- 2-sssserstein distance. By pre-finding using convex optimization algorithm and linearly drifting the coordinates, we can assume x. without loss of generality.
DRUL admits the stepsize reaching O
22/a and the overall complexity Ö (t)(H))
Now we consider a realizable case where our result achieves better complexity on d. The concrete example is that admits separable structure discussed in Section 3.3 Corollary 4.3 below shows that when f admits a ridge separable structure, Algorithm enjoos dimension-free iteration and gradient complexity.
Corollary 4.3. Follow the notations in Theorem 4.2 Further, if f(x) admits a ridge-separable form and satisfies Assumptions 3.4 3.5) Then if we set the step size
2/3
min
12C2x' (24C2k( + 11)) 1/3
With initial point (Xo, vo) and Qo Ex-p.v-N(0..)) (1xxo x12 + 11x0 + vo - - --x13.. Then when
log Sek
?
Algorithm outputs x,, such that W2(Taw(X.). p) 6.
Theorem 4.2 establishes the convergence rate of DRUL with weak dimension depen- dence. We shall note that for lots of high-dimensional sampling problems, the trace of the Hessian matrices for the potential function is much smaller than times the largest eigenvalue because the eigenvalues often drop rapidly. In practice, this situation has been considered in lots of topics. For
Discussion.

---

Egensaluue spectrum the Hessian
10k
1000
130
30 hrdden units
hidden
120
100
$20
692
700
(a.
(b)
A demonstration of the eigenvalues of the Hessian matrix. (a) The eigenvalues of the Gram matrix of MNIST data. (b) Eigenvalues of three-layer neural network from Sagun al. 2016].
Figure I:
example, one defines the effective rank (e.g. Hsu al. [2012] to represent the acting dimension of the data on linear models. In distributed machine learning, one can show fewer bits are needed to be transmitted between machines when the eigenvalues decrease fast (e.g- Hanzely et al.) (2018). As shown in Figure T(a]] the Gram matrix of the MNIST dataset, which is also the Hessian of Bayesian ridge linear regression, has rapidly decreasing eigenvalues. Similar empirical results are observed on deep neural network models by Sagun al. (2016), in Figure [[b) In fact, beyond the ridge separable case, there are many problems admitting tr tr(H) = of o(d), such as the concentrated posterior distributions with bounded gradients and neural networks with regularization. However, notice that the result of Theorem 4.2 is based on uniform upper bound the Hessian matrices. We think this is the basic case to show that algorithm achieves convergence with an intrinsically low-dimension dependence. It is possible to relax the condition to the upper bound the traces for local Hessian matrices with additional Hessian smoothness assumption. Please more discussions in Appendix
may
Intuition
In this section, plicity of the notation, let (x)
provide the intuitions of the algorithm design. For the sim-
[8, Xn; (B, Jishs (Xn, vx)) and xx) Ts, (a,x;; Xn; (BrJostso; (Xn, vn)); ;); (BrJosiso) (Xn: Va)). Then given a and in Algorithm (s) = [s, (1,, (a);; (Be)) (B.ooeegs, (Xn: vn)) and x,, (3) X,+1-
Our analysis tracks the dynamics of the distance to the stationary distribution. Specifically, analysis the dynamics of (t) = xx, - xx (f) + (t,) + +E] (t, x (f) and EΩ, EΩ,, (0), which upper bound the squared 2- Wassersstnin distance W₂ (LaL (t)), p) and
Wy (La respectively. Here, (x", (f), v" (t)) evolving along the exact Langevin diffusion
3.4 follows the stationary distribution and is synchronously coupled with (x,((), v,, (t)). One can
find its formal definition in Appendix B Via tracking the flow and using the contraction property of the process, the following Lemma characterizes the one-step discretization.
Lemma 5.1. Let X,, be the n-step output of Algorithm and x", (t) be defined as above. Set
any h L+m Under Assumptions 3.1 given that X,, and x", are coupled synchronously, we have for
E (B < Egnae EE2, +
(5.1)
where is
È 2uEEarp E Egpp 1 2 (x) (3()) vn((s) v.) Vf(x, (a))
Lemma indicates at each step, EO,, contracts with local discretization error. Telescoping
5.1 and upper bounding the term E yields the Theorem 4.2 To obtain the convergence rate in the

---

main theorem, we will show that Algorithm| I guarantees an upper bound of & which (1) achieves the state-of-the-art dependence on stepsize and (2) satisfies that the dimension dependence can be controlled by the contraction.
Improved discretization error. We accomplish the first goal by matching the expectation, which is detailed by Lemma 5.2]
Lemma 5.2. Let p be probability measure defined on [0, k] satisfying p(dt) oc ( knd dt, and the probability measure p' defined on [0, h] satisfies p'dtt) (e t- - dt on [0. hl. Then for measurable function F(t), p and p' satisfy that
(A) There exists positive constant C1 such that Ec-p Soe " F(t)ds C2hEe F(t).
(B) There exists positive constant C2 such that Ex-p/ [F(I CRepp
Claim (A) in Lemma 5.2 states that by choosing a random step size B p', we can leverage the low discretization error of the randomized midpoint method. Denote the random weight w., a)
Xn (8) - x (+) + vn (8) (s). One can split w,, (8, a)) into (wn (8, a) w., (0.00) 0)) w. (., 0). The former term can be bounded using the one-step move, and by claim (A), the dominating latter one is
2u EuEBarp E Ea-p l - ((xx,)) ([(R(())))
-2uChh EE E (V](xx(+)) - (+))))
which attains low discretization given that Bgge V (8),(()) is low biased approximation Eg-p/ (x),,,))).
Averaged contraction can control the weight variance. Then we consider the variation of the weight un as) The time difference of wn(s a) writes
w(s, a) (xn +V Vn v% / (vn(r) (.(r) um((x xx,
(5.2)
(r) v,,(r) -V xVf +
(5.2) indicates the variance Var Varssp (w(s, a)) is dimension dependent for nondegenerated distribu- tions on [0, h), since it will introduce the (r) whose difference to (3) or vn (0) is dimension dependent. We control the dimension dependence via the averaged contraction. And the randomized step size makes it possible to consider the averaged effect. We have the following Lemma to bound the variation of the weight.
Lemma 5.3. Let x, (t), (t)(,,, vi and vi (1) be defined as above and u L+m Under Assump-
tion 3.1 for any 1, a S h, we have
Xx (L) X x (1) + x + vnn(1) vn - (1( + vlll?
<10 EE 2n (1) + E 2 L (r(H) + 11111.11
(5.3)
m
Now the dimension dependence of Var w(s a:) is contained in Ex (t), which can be controlled using the averaged contraction under the stochastic step size.
Conclusion
This paper proposes
technique and designs the DRUL algorithm. We prove that with strongly convex and Lipschitz smooth assumptions potentials, the algorithm converges the target distribution in Wasserstein distance in Ō 1|x iterations. The result illustrates
that many sampling tasks in machine learning can achieve a dimension-inndeprndenn complexity. The proposed DRUL algorithm can be potentially much faster than existing algorithms for high- dimensional problems. As a concrete example, when the negative log-likelihood function admits a
ridge-separable structure, under mild conditions, dimension-free Ô can obtained by DRUI.. We hope our technique brings new insights for designing dimension- independent algorithms for high-dimensional sampling.
(273 iteration complexities

---

7 Acknowledgement
C. Fang was supported by National Key R&D Program of China 00022D0114900) the NSF China (No. 62376008) and Wudao Foundation. T. Zhang was supported by the General Research Fund (GRF) of Hong Kong (No. 16310222).
10

---

References
Jun Liu and Jun S Liu. Monte Carlo strategies in scientific computing, volume 75. Springer, 2001.
Steve Brooks, Andrew Gelman, Galin Jones, and Xiao-Li Meng. Handbook of markov chain monte carlo. CRC press, 2011.
Yang Song. Jascha Sohl-Dickstein, Diederik P Kingma, Abhishek Kumar, Stefano Ermon, and Ben Poole. Score-based generative modeling through stochastic differential equations. arXiv preprint arXiv:201 13456, 2020.
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, 2017.
Alain Durmus, Szymon Majewski, and Blažej Miasojedow. Analysis of langevin monte carlo via convex optimization. The Journal of Machine Learning Research, 20(1):2666-277 2711, 2019.
Xiang Cheng and Peter Bartlett. Convergence of langevin meme in kl-divergence. In Algorithmic Learning Theory, pages 186-211. PMLR, 2018.
Santosh Vempala and Andre Wibisono. Rapid convergence of the unadjusted langevin algorithm:
Isoperimetry suffices. Advances in neural information processing systems, 32, 2019.
Sinho Chewi, Murat A Erdogdu, Mufan Bill Li, Ruoqi Shen, and Matthew Zhang. Analysis of langevin monte carlo from poincar', log-sobolev. arXiv preprint arXiv:2 126. 2021.
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.
Alain Durmus and Eric Moulines. Sampling from strongly log log-concave distribution with the unadjusted langevin algorithm. 2016.
Yurii Nesterov. Introductory lectures on convex optimization: A basic course, volume 87. Springer Science & Business Media, 2003.
Yi-An Ma, Niladri S Chatterji, Xiang Cheng, Nicolas Flammarion, Peter L Bartlett, and Michael I Jordan. Is there analog of nesterov acceleration for gradient-based meme? Bernoulli, 27(3): 1942-1992, 2021.
Xiang Cheng, Niladri S Chatterji, Peter I. Bartlett, and Michael Jordan. Underdamped langevin meme: A non-asymptotic analysis. In Conference on learning theory, pages 300-323. PMLR, 2018.
Matthew Zhang, Sinho Chewi, Mufan Bill Li, Krishnakumar Balasubramanian, and Murat Er- dogdu. Improved discretization analysis for underdamped langevin monte carlo. arXiv preprint arXiv:2302.08044 2023.
Ruoqi Shen and Yin Tat Lee. The randomized midpoint method for log-concave sampling. Advances in Neural Information Processing Systems, 32, 2019.
Andre Wibisono. Proximal langevin algorithm: Rapid convergence under isoperimetry. arXiv preprint arXiv:1911.01469, 2019.
Ye He, Krishnakumar Balasubramanian, and Murat A Erdogdu. On the ergodicity, bias and asymptotic normality of randomized midpoint sampling method. Advances in Neural Information Processing Systems, 33:7366-7376 2020.
Xuechen Li, Yi Wu, Lester Mackey, and Murat A Erdogdu. Stochastic runge-kutta accelerates langevin monte carlo and beyond. Advances in neural information processing systems, 32, 2019.
Ruilin Li, Hongyuan Zha, and Molei Tao. Sqrt (d) dimension dependence of langevin monte carlo. arXiv preprint ar arXiv:2109.0388 2021.

---

Maxime Vono, Daniel Paulin, and Arnaud Doucet. Efficient meme sampling with dimension- convergence rate using admm- type splitting. The Journal of Machine Learning Research, 23(1): 1100-1168, 2022.
Yoav Freund, Yi-An Ma, and Tong Zhang. When is the convergence time of langevin algorithms dimension independent? a composite optimization viewpoint. Journal of Machine Learning Research, 23(214):1-32, 2022.
Gareth O Roberts and Richard L Tweedie. Exponential convergence of langevin distributions and their discrete approximations. Bernoulli, pages 341-363, 1996.
Jonathan C Mattingly, Andrew M Stuart, and Desmond J Higham. Ergodicity for sdes and approx- imations: locally lipschitz vector fields and degenerate noise. Stochastic processes and their applications, 101(2):185 185 232, 2002.
Arnak $ Dalalyan and Lionel Riou-Durand. On sampling from log-concave density using kinetic langevin diffusions. Bernoulli, 36(3):1956-11) 2020.
Alain Durmus and Eric Moulines. nigh-dimensional bayesian inference via the unadjusted langevin algorithm. Bernoulli, 25(4A))2854-2888 2019.
Andrew Gelman and Jennifer Hill. Data analysis using regression and miliieee//hhiierrhhheed models. Cambridge university press, 2006.
Arthur Jacot, Franck Gabriel, and Clément Hongler. Neural tangent kernel: Convergence and generalization in neural networks. Advances in neural information processing systems, 31, 2018.
lean-Yves Audibert. Fast learning rates statistical inference through aggregation. 2009.
Wenlong Mou, Yi-An Ma, Martin Wainwright, Peter L Bartlett, and Michael Jordan. High-order langevin diffusion yields accelerated meme algorithm. The Journal Machine Learning Research, 22(1):1919-199) 1959, 2021.
Yin Tat Lee, Zhao Song, and Santosh Vempala. Algorithmic theory of odes and sampling from well-conditioned logconcave densities. arXiv preprint arXiv: frXiv:18/2.06243. 2018.
Levent Sagun, Leon Bottou, and Yann LeCun. Eigenvalues of the hessian in deep learning: Singularity and beyond. arXiv preprint arXiv:1611.07476, 2016.
Daniel Hsu, Sham M Kakade, and Tong Zhang. Random design analysis of ridge regression. In Conference on learning theory, pages 9-1. JMLR Workshop and Conference Proceedings, 2012.
Filip Hanzely, Konstantin Mishchenko, Peter Richtárik, and A. Sega: Variance reduction via gradient sketching. Advances in Neural Information Processing Systems, 31, 2018.
Gilles Hargé. A convex/log-concave correlation inequality for gaussian measure and application to abstract wiener spaces. Probability theory and related fields, 130(3):415 2004.
Sanjib Basu and Anirban DasGupta. The mean, median, and mode of unimodal distributions: characterization. Theory of Probability & Its Applications, 41(2):210-223 1997.
12

---

A Simulation of the discretization process
A.1 Integral form of the discretized process
Proof of Lemma 4.1
Proof. Define the aligned vector Zn (s) (x), (s)) and write the diffusion process as
dz (s) (s) uu((,, ds 2/dBdB,..
(A.1)
where
I
=
umI -2 a0=
(xx) B
0 I
And B, a Brownian Motion Red, By a variation of constant method, the solution of A.1 1.1pr is
(1) -0 (0) (f. elo--t))ss a 2/ / eallttt)
To obtain the explicit expression of x.) and (t),),, respectively, it suffices to present block-wise form
An()) I A21(()
At
(2)
of the exponential At And direct decomposition yields that
/22 Aum - 2 (-2- 1-2-v22-4um)t//
22 4um
(n1
6 (-22v22-4uum)e
2/ 4um
2/22 twmn
A12(t) V2 4um
el ((2- -u-m//-uun
el- (-22 (1-2+222 2" 2umn4mm)t//
V22 4um
(A.2)
(21 2 4um
um
(-2 22 4umjt/2 +
um
/22 4um
A22(1) 2 v 22 4wm
(-2 1um
2 4wm
2
2² 4um
2
4um
Since zn (t) = (t), vn (t)), we obtain the close-form solution of the discretized Langevin process of x, (f)
(t) (t)x)), + A12(t)v, + S/'' A12(s A 2/u J A12(s - t)dB.
similar proof can be applied to the exact diffusion process and thus (t) satisfies
x (t) (n)()),, Ara(t)vn + +u L A1z(s - (x) + 2/u f' A12(s t)dB..
Thus we can obtain the result in Lemma
A.2 Controls of the integral forms
Lemma A.1. The functions An (1) and A12(t) Lemma 4.1 satisfy
[A12(t)] St, An -
fort E [0, h..
13

---

Proof. In the proof, we will repeatedly use the inequality e¹ - = O(t). For notational simplicity. denote vI- umr by k. For the first claim,
(12 1
kt
1
1+kt
=e-
1 e
ekt - -1 ) O(1).
24:
2k
2k:
2k
As for the second claim, we have
A12(t) - k - 1
k + 1
t-
+
e
-t+kit -
2k
2k:
Ek
E
+
tekt 2
+
1
e-
ekt -
) Okl) + O(1..
=
2
2k
2k:
Note that h /1 I 1 and h 1 1. When S h our claims hold true.
A.3 Simulation of the Brownian motion
Lemma A.2. Let B, be a standard Brownian Motion. Algorithm 1 we have 16 A12(s () H = A addB, and W So A22(s - l)dB.. Then conditioned on the randomness of a and B, (G, H, W) follows a zeron Gaussian distribution with covariance
(min((,aa)
[ RG) [H)
A12(s A)A12) adds I,
[B(( (G) (G)G)"))))
A12(s - d) I,
B[(H EH) EH)"1
A12(s - aa)2ds I,
(A.3)
(%, A22(s - ds 1,
B[(W - SW - EW)"
(min(a..)) 112 - A22 - t)ds) I
[[(W (H - RH)")
( A12(s A) t)ds) I.
B[(W- EW) ]G)"
Proof. By the properties of Brownian motion, (G, H. is zero-mean Gaussian variable. Its variance is given by
B[(G EH) - -E [
A12(s l)dB, /' A12(s aa)dB, minn(t,))
A12(s ()A12(s -a)ds a)ds a)dd)1 I.
The rest of the claims can be proved similarly.
B Proof of the main theorem
In this section, we present the proof of Theorem 4.2
B.1 Contractive lemma
begin by defining the n-step exact process (x,, vit) and (x),,,, vū(()). Denote dy" x -U())- l|/ dxdv. Let (x,, vo) ~ p" and
dx = (()d(, (l) = (x) + V2udBr, (x, v0) = (x (B,-1),x;_1 1. (3n-1));
(B.1)
x"
14

---

For the simplicity of notation, we denote (x,, v,;) = (0), v. (0)). And we have for any and t, (x,, (t), V, (t)) ~ p². With slight abuse of notation, let x, (t) and (t) denote inaccurate starting point discretization process and midpoint process in Algorithm respectively. And recall that in Algorithm| we have
(1) = (,)(,), + (12(v), + A12(5 ()V-tt))(xx)))) + 2 / A12(s ()dB, (B.2)
and
(t)(t) =An(tt)x, + A12()v,, + / A1z(s ()),s) + +2/u /. A12(s t)dB..
(B.3)
And we assume (t) is driven by the same Brownian Motion as (t) and (t) at the n-th step.
Define 2,, (t) x (1)l (t)
(() v (1) and Ω ||x, x 12 |xx,, x + Vn - v 2 The proof concentrates on the
x dynamics the Ωn. Note that Wasserstein
distance minimizes all the couplings and 22 (t) 2 ||x, (t) x (+) Hence EΩ2, upper bounds of the Wasserstein distance (VaLL((X,)))) p). Given synchronously coupled assumption, the dynamics of (t) can be characterized as follow.
B.1.1 Proof of Lemma 5.1
Proof. Let the random process B given B,, be
B. B.-B. d
0 t
dB;
a-t dB
where B; is an independent Brownian motion. Note that B, is Brownian bridge and B B, Then when (a) is given and x,, (t), x, (t) are coupled synchronously, given Ba, we have dx (1) = v,, (1)dl,
(t) -u (t) -V
(t 2t) 2 2/ddB;
and
(x,,() (t)dt,
dt,)(t) d)) -vv + 2V
For simplicity, denote (t) (t) x and (t) (t) (t). Then by Taylor's Theorem
+ VU(X, (t)) (U(xx = / (1 (1- (t))dszn(
(B.4)
Define H4 = fo V-U(sx, (t) + (1 s)xx (f))ds, then
dl,, (t)/dt =2 (t) + (t), vn u(g(xx - g(x)(())) - -u "VV((E((a)) - vf(xi(())) -V(((((())))(((((()))) 78n(())
+ 2 2n (t), 2/ (t)
= - 2 (z, (1) + 4/ (1), (7 - (1)) + uHi,,,
(2 (z),
+ 2u(zu (t) + 2h (1), Vf(x, (a)))
where in the second equality use B.4] and divide f((,((a)) into V fx(x (t)) and the difference vf((((a)) - (x)x)())). Note that A can be converted to quadratic form
A = nm (t) (t) 2n (1)
1
1 juH, 1
]
Zn (t) + (t) 2, (t)
JuH4 -
15

---

whose eigenvalues X; are given by characteristic function (1- ^)² (ud:(HH)) + um - -2)² 0. Given the strong-convexity and Lipschitz smoothness of the potential (xx,, when is set to be T. the eigenvalue of greater than 2x' Therefore,
(1)/l) S - (1) + 2m(zn (1) + (4),), (V)-)((().))))))
K
Then by multiplying e# and taking the integral form 0 to B, we obtain
(3) Se (0) + 2u !-²
e -
(xn(s) (s) + (x) vi(),, ((xx(())) V)(RR
In Algorithm the random step size B ~ p and thus we reach that
EΩ,-1 S Egnpe 2 EQ, ZuEBBag 1oe -- (=((xn()) (xx()))x;(()+++ + Vv (b), VS(xx())) - -V((s()))VV((((((()))))) (a()aa))))d..
B.2
Proof of Theorem
We devote the rest of this section to the proof of the main theorem. Lemma 5.1 shows that the dynamics of the EQ, is driven by a shrinkage term and a local discretization error. To achieve the claimed convergence rate, the proof aims to establish an upper bound of local error which is: (1) dimensional-indee -nneepieeenndddeeeeeeenn,, (2) optimal in the sense of max step size h.
Proof. Given that EΩ,, is an upper bound of the squared 2- W-sssersteinn distance, it suffices to prove EΩ2, S e By Lemma 5.1 and Lemma 5.2
we have
(9) S B= EN
1o
+ ZuBE gmpp
(x, - (x) (1))+vn())) -10), - (a))
= ESSEE
En
2uhC EE
B (3) (3) -Z))) VS(x, (3/
(o))
A
(B.5)
Next, we bound the error term A by
A =2= 2CwuhEEB3... 3mpp (B) - X,, xx) +x) + (B) v, - (3) +v,, V/(xx, - (3)(((()))))((
+ (CCCuhEExs, (x, x + Vn Vn Eanp' Esnp (Vf(xx(8)) (a())())))
wh uh
E x,, (3) Xn x (B) + x + (B) Vn (3) + vill?
+ u 22 V
V/( (a)))1 LuhBBli - X vn v
whh Ea B (x (8)) - vf(xx
h h 1EE,,,2,,113 (2) En, + 2h (H) + h 111120 112)
m
+ 22 1 + h2L tr( + 1'3 ||-, 1122) + 1E
m
+ 22 I1LL"ED + 4
tr(H H + 1 | 112
m
2 2 I
S (3) + h E2, +
tr( H + ntt xx.1
16

---

where in inequality S. we apply Lemma D.I Lemma 5.3 and Lemma D.5 inequality follows by the observation L constant C2 such that
L+m S Then by plugging the upper bound A into (B.5 there exists
T tr(H) 12 11x1111)
EQ, (B) B EQ, + C2
I + + m
(B.6)
B.6 can simplified as
(3) Baupe + O2h2
O2h
( + ||x,1 11x.11²)
Enn
(B.7)
1 - C2h
C21
m
Note that when ~~ (3) = EΩ And EΩ,, is an upper bound of Wasserstein distance.
Equation [B.7) characterizes the one-step discretization of the Langevin process, where the first indicates a linear convergence if without discretization error, and the second term is the one-step
discretization error.
Without loss of generality, let C2 1. Since h K, we have Ex-pe = e 1
Since S 120 we have h I 2C2r and S VC2
and therefore C2h2 2 0, 2C2h² I 6x
Hence the first term of B.7 can be bounded by
E Egnpe + C2he2 EQ,
1
O2h2
3x
EP
O2h²
1 C2h2
h
+ 202h2
+1 E2
C2l12
(B.8)
h
h
3x
6K
+1
EPn
Czh2
h ERns
As for error term
B.7 tr H) |x. IP) and therefore
since h S 120'x S V2C,'
we have
Caha I-C2h² S 2C2ha. Denote
6x
and & 2C2h
m
EΩN S rEEQN-1 + E < NED + E(1 + +
N11
NE + 5
(B.9)
1
When N >
6 log
2B52 2²
we have
2
Nh N NBO S e eon
EDO S
(B.10)
2
And when h S
2/3
1/3, we
(24C2k( 22 (lx))
E r 2
<
(B.11)
2
Plugging (B.0)) and B.1 into B.9 yields
EΩN < 2
and therefore completes the proof.
C
Supporting lemmas
We assume that Assumption 3.1 holds in this section and Section D
17

---

C.1 Proof of Lemma 5.2
We would like to point out that how match the random step size a and p is critical for obtaining low-bias estimator. And the control for local discretization error in Lemma [.1 involves a challenging multiple integral. These considerations inspire our choice of p and p'. And in the following lemma, we summarize the main properties of p and p'.
Lemma C.1. Let be probability measure defined on [o, h] satisfying p(dt) x (h. hnm dt, and (e xx" - :)) dt on h). Then for the probability measure p' defined on [0, h] satisfies p' (dt)
measurable function F(t), p and p' satisfy that
(A) There exists positive constant C1 such that Ee-p ee F(1) = C.hBeap F()).
(B) There exists positive constant C2 such that Eenp' [F(I) S C.e)alF(())
Proof. We begin with the first claim. By the definition of and integration by parts, we obtain
E [1 (' - F(s)ds
h I
1 -
E=
F(s)dsct
1
l e F(s)ds) 1:a P
1. h F(t) 1 [e F(())ss)))))))) de
h
k
F(t)dt - F(t)dt+ 1/' ho L c F(s)dsdt, F()ddsd,
(C.1)
which indicates that
(1 [
I == F(s)dadt = = =1% (ctth F())dt.
hi
has a normalizing constant which is of order @(1). Denote the normalizing constant of ds by scalar function Zh)h). To obtain the first claim, it suffices to show Z'(h) is of order (h). Then given that e - 2 1 ht Zh) satisfies
x h hu Lh
Jo e -
1) z'(h) 1." h (e (5= - di (1-e
c-h
h
h
h
<
h
2
2
which proves the first claim.
Then we prove the second claim of the proposition. Note that when h S K
Z'(h) = K. (1 11-e-
11-
1
1 e
2
=
-
2
2
2
Thus
1
1 1
e
e
S
(C.2)
z'(h)
Z'(h)
h
As for denote the normalizing constant of p by 7(h)
f
1
t hv
Since h = o(1) and k > 1,
h
h
Z(h) = = (1))
2x
Then we have
1
Z h
(C.3)
((h)
h
jk
h
where the last inequality follows by that S 011) and i. 1. Combining (C.2 and C Claim (B).
yields
18

---

C.2
Convergence lemma
Lemma C.2. Let x, be the minimizer of the potential function, that is, x, arg min(g(x) + f(x))
and A is positive definite matrix that ml N A = LI. Let (t) and x,, (t) be defined as (B.3) and Assume that L+m Then we have the following bound
BB
(0)) tr(A) + 12 113 E| Ax (t) (t)
m
AI (1) (4)
LA x (1) + v.) (()( - v1, tr( x1,
and
L ttr(A?) + L' 1|.111
18 (x))) &
m
Proof.
We begin with the first bound. By Young's inequality and Lipschitz smoothness of f.
E Ax (x (1) - ()) xx.0 x.) 2 |x
(C.4)
< xm, (t) ()) xx)]|
+
Since (-) is m strongly convex and || ||||| is convex function with respect to x, Theorem 1.1 of
Hargé [0044] implies that
1
E (x (t) Ex²) 11 S Bx # |x
ttr(A)).
(C.5)
m7
On the other hand, by Theorem 7 Basu and DasGupta 1997].
(1) x. ). (Ex; (1) x.) S 3,
where Σ is the covariance of p and therefore bounded by 1. Hence
|A(Ex; - x.)1|² < ||-||| I - tA(()))
1
(C.6)
m
Combining equation C.5 and (C.6) we obtain that
AI|A(x; x.) || < E A(x; Ex;))1 + + x.) - Ir(P))
(C.7)
m
Combining (C.4 and C and then we achieve the first claim.
To bound E Axx, + V., ( we have
(I),(x,,t) + (t) E A(x, (t) + v,, - (t) (1)
<L' E + (t,)))+)) x, 20 + (()((X) x.)11
+ - (t) v0]l + tr((²)+ +
1 tr(AA) + 1111
+ (()
m
m
1 <L'EIx,(1) + (x) x (t) - 2 tr
+ L |x, 1 112
where the second inequality follows by the Lipschitz smoothness of f(x) and Young's 's inequality; in the third inequality
C.7 and ~N(0, )
To prove the third claim, we have
(IIAVf(x, (+))
A - (/1
(x.)((x))
<L E X (t) + + A(xx A(x))
xx)Il 2 +
17
ST. (1)
m (r(?)) 1.2m?11x.1
where in last inequality, apply C.7 and ff(x,) + mx. 0.
19

---

D Technical bounds
In this section, we present some controls that are useful in the proof of Theorem 4.2
D.I Bound of (x)))(x.()) (=(
Lemma D.1. Let x,, (1), (t) defined as above and L+m For any fixed 1, a S h.
(i) When Assumption 3.1 hold and H is a uniform upper bound i.c. v2f(x) = for all x e Rd, x (f)and x, (a) satisfies
h? L3 |
E| V (8)
+
m
iii) When Assumptions 3.1
and
hold, we have
12L
E [VV (a))l $ 1 I, +
M
1,2 x.
123221
m
Proof.
For any x, y E R, we have
[Vʃ(x) - ()x)) 11 4
f(1 v²
k)x + ky) y)dk
(D.I)
k(y x))dk and H depends on and y. Substitute x, x.. (,, and then obtain | fx,, (t) (a) [[H((x(()) (a)a)) Then by plugging in the close-form solution into |H, (t) - (a))
Denote H - -
- - (a))l1
Hl
An(a)) An(()) - A12(t) + A1((0))x, (A12(1) A12(a))(xx, +n)+
- A12(s- - (a))s -u 1/' A1z(s - a)(ss))Vf((x)d 20 /' A 2
A12(s
2 / - acdss
(tAn((t) AAi(a) A12(t) (12(a))211Hxx,, 1|2 (A12(t) A12(0))3117((,, + vv,)l
+ 2² S A12(s (to--ttWW/( L-d A12(s - (x,,)((x))ds
+ //" A1z(s
A12(s- - t)HlB Aly - f A1z(s HdB
(D.2)
where in a we use Lemma 4.1
follows by the Jensen's inequality. Besides Lemma A. implies that
(An(()) (a) A12 (l S 12ll Hx, and (A12(1) A (x, +
vn)ll2 1 H*,, 2
20

---

For the expectation of term ① we have
A12(s - ()) -.o aa A12(s-a)) a- (x,
E
% A12(s () - vs((x))) 1 +
2
! A12(s - (al(s-a)AV/((s// + (A12(s - -12((s (l))Vf((x)))
E
(D.3)
<Et
A12(s ( Vf(xx))]12
+ E (t) /- 4 [x,)(>))1
EEt / (A12(s - t) A12(s d))
S (VF(,,())) 1f
where in the last inequality we use that a < h and A12(t) O(t) Let H = (a) (a) )d.)... By the definition of (x) EIiH (V((x((a)))V(((())) 112 can be upper
bounded by
E H. (V vf(x,))
<L'E|
(V)(&,(())) (x(x,)) 11
(((((nn(a) 1)x, A12(a)), +u (-" A12(t (oltta)Fffaal (xxx + +2v0 f" A1z(t B12(taa)BB)
=L"E
<LLEAA1((a) ! + 3
A12(z)
1) 2) x,,''
+ B2x 1 1211A12(1 a) H' F + WEI!2 1 A12 - a)H'd
2h2L² 17'' (Xen
+ 2 H "
// A12(1- - a)' a)'BBB
+ uL u
where in the last inequality we apply Lemma A.1 Note that x,, x,, = X,, and v,,(0) = vn- Then by Lemma
E
- vf((x)))
<I. L (a)) - vf(x,)))
2 Eir(H -x (|(xx11))) + (12)))+1((1111111)1 11 xx,1
2/2 (L (L* + -x-- - -
m
+ 1 Etr(H(?)
I.2
1211x, 112) u² ( (L Xn x,,12
Etr(H(?) I I'm² m2111x,11?) (D.4)
m
m
10-
+ u
Aiz(t F
121,2 5he L L'EO, + Etr(F(?) + h² 4 ||x.|| + u 1/ Arz(t Anz(t-a)f''a a)
where in the last inequality we use - < min { }. Plugging D.4 into D.3 and then we have upper bound of the expectation of term @ Plugging the upper bound of term ! into (D2)).
21

---

and thus we have the following inequality
E VIV v
< 7x || +
1 + + (xx))]
1B + 1 Etr (FL')) + 1| 11 uEL? 1/ A12(t - a)H dB l')
+ h²h²
+ wE /
A12(s -t)HdB, / A
S1aT2 + h? E (tr ((F())) + tr (FE) + 1.27 || xx +12 12,,,,2
L A12(t -a)
m
+ WE 1 A12(8 - t)HdB - l0" A12(s o.)ddB
where in the last inequality we use Lemma and the definition of 9n
When assume V [x) has uniform upper bound H, have
||
wE
A12(s t)HdB,
A12(s - H)HdB.)
2
(D.5)
h3
h
u.E
hHdB,
tr(H²)
tr H
m
m
and
hh LZB
1/ A12(t - a1)''B))''B, S Su
S hat (r(())) =
htL
(r(H).
m
Besides we have ((H)))) (722 s Ltr (H') Ltr (H) s Ltr(H). These bounds yield the control
12T. tr + 12
V (t) (a)il L +
(D.6)
And when is -e-Hessi smooth, by Lemmas D.2 and D.3
we have
V (t)
54 L + h2 E tr (H²)*) + tr (HP) 3 LM + 12,2111x.11
m
+ L3 d2 + L2 2 10d2
(D.7)
S,,, + I. -. (tr (H') + tr (H)) + 1
112
M + 1.,,1|x.111 1|+ 13
m
m
1,R,, + LhnM
xx 1
1.
m
Lemma D.2. Let H. be defined as above. When Assumptions E.1 and E.2 hold, we have
E % A1z(s
t)HdB. - o A12(s aa).
< LL + 12gh 82
22

---

Proof. First, we have
A12(8- - 1)FAdB, - S-¹ A ||"
H (L/ A12(s - ()dB, - S/¹ -A12(0 - A11((-a)ddB)
H E (H])F)) (%. A1a(s-)))lB
1)dB, - 1/- A12(s- - a)dB,
S
(D.8)
E [M (o A12(5- A12(s-1)dB, -jdB, - 1/¹ A12(8 A12(8-
+
72 E [m] (1." (o((s-))dB, - J A12(s (olssa))BB) )1)
S
E H|)T A12(s - -)ddB - 11" A12(5 A12(5 a))B )l²
where S follows by dividing H into E' [F]] and H - E [[]F)] In the last inequality, EE [H] Fr.] is independent of B,, and thus
A1a(s-()dBB ( A12(s (oo(s-6ddB f)dB, 1088--oo A12((- ajdB' ||
E E
H|
L/- h
(D.9)
S1²B
E [M]]FM]]
cB,
ShºBtr ( (H(
As for the firs term, let [B'Jocset be a standard Brownian motion which is indepen- dent of [B.Jossh. Define & (a)
J(a, Xni (B (X,, vn)) and x,(1)
Haa (a); (B Jusase (Xn; vn)).. Then we have
ΣE [ || (t)x (t) (a),
E B - E |1
+ ( dk
LB /
h
- |11--)) (x))))xx(tt)) x((()) - (a))
$1.20P1 PB
c
1'' A12(8 (V VI(2( ds
1//0 A12(s -a)dB, - / A12(s a)dB'
+ L2 u
1./2 A12(s - -tdB, - //" A12(s-t)4 -
+ L3 L3 WE
(D.
where 2 follows by the definition of H. and Jensen's 's inequality; in S we use the Hessian smoothness
of S follows by the definition of (t) x", (t), (a) and & (a). For the first term of the last inequality, by the definition of xn, (a) and X, (a), we have
ST.71Au 1./0 A12(s (hass-addB, - 1." A1z(s - a)dB'
23

---

And thus we have
([/ A12(s - l)dB, f A12(s -a)dB. )1²)
|| - [H] 1/],,11)))
1// A12(s t)dB.
EL3u %// A12(s - t)dB'
5 <LguB
(D.11)
1100 A12(s a)dB,
+ L3 1,00 A12(s a)d
+ LLzuE
LLzgh&t de
Plugging D.11
and D.9 in to control (.8 yields
1/¹ A12(s - (P)dBB)
E
A12(s t)HdB
Sh Etr (IB(H(),1)
dd
<h LM L
where the last inequality follows by Assumption E.2 and the observation that Etr (EE()) - LE (H) [F]] S LM.
Lemma D.3. Let H and H' be defined as above. When Assumptions E.2 and
hold, we have
uL'E S Ala(t -ha(t-a)Hadp S h M +
Proof. Similar the proof of Lemma
D.2 have
uL E
A12(t a))'d
H'' /
=uL?B PE
- ad)d
(H'' H [7]) 1'' A1z(t - a)dB, aa)(BB
2
E [7(] /-' A12(t - a)dB,
SuL"E
+ uL'E
| 1171'- E H |7' I| 1// Arz(t -
he(t-o)dBB) o)dB, oo/dB.))) + ut 1 (B (F(
SuL'E
(D.12)
We follow the notation of BB, and X (a) in D.10). Then
| [ 1| 1/00 Ara(t a)dB. I'
E
(1- k)) 18, (a) & (a)ll dk 1/ Ala(t -ofttaalo (ale-aadBB) (")
<E
L3
(1/% A1z(s a)dB, /0 A12(s - hhassB))BBB) | 1/1// A1z(1 -a)da) ()
(D.13)
SLIE
uu
10"0
LLLuB
A1z(s A addB,
SIzuh&d'
24

---

Combining D.12) and D.13 and using tr Frl) I LM yields
uL E / A12(l-
S SuL + 13
D.2 Bound of (t) X, -x) +v(() Vn (1) + +1
To prove Lemma 5.3 consider the following Lemma which also includes the Hessian smooth case.
Lemma D.4. Let x., (t), (t),),, (t) and (t) be defined as above and t L+m
ii) Under Assumption
3.1 for any t, h, we have
B x, X, (1) + x + vn(()) - Vn vv) + vill?
(D.14)
S p On (1) + 1.2 + 2 ir(H) + 111111.1??
ii Under Assumptions 3.1 E.1
for any t, a h, we have
El (1) (x) X vl)
v2 + vill
Sh Emp + En2 + u h4 LM
+ h² 1x11>.1' 112 L3
m
Proof. Since the discretized process and the continuous process are synchronously coupled,
x,(()= X, x(, + x + (() V, -,-v,
v1(() + vi
=/6 (V)(s) (s) um -s))-)((vv(((
+ (UV((())()))+++((((()()())))))
(D.15)
Hence the square norm satisfies
(t) X - xx(()) + + (t) v, - (t(() + v
5 31 / [|)(() - (8) +
xx, (s) (x) + 3(1 um) 1/ (xn(s) - (s)
+ +uu2 ( [VT((x (s(x
33t I (n(ds + 3u²t / [V)({x(()) (s)
(D.16)
2 where 2 follows by equation [D.55 and Jensen's 's inequality; in S we use 11-um| < and substitute square norm by Pr (x). bound the second term in the last inequality, we have
- vf(x)
1 [V)(x.(s)) (f(x (s)) /
(D.17)
2 S (n(ss)ss+ [VV((R(())))) - [[(x,(s)]
25

---

where in the second inequality, we use the Lipschitz smoothness of f and (t), xx) - (t).
By combining bound (D.16) and (D.17). we arrive at that
B Xn X, (1) + x + v.) Vn. - vv + vl
(s)ds u²t / (( (s)
St
(D.18)
<h²
Th EE (1)dl + w²11 + (IV - v (+)]l
h
By C.3
we have h? So h (t) h h) EΩ2 (t) S h (tn
When the Hessian of f is upper bounded by H. by D.6). we have
E X. (t) Xn, - x (1) + x + VR (t) + v
h
Eh'EamEElln (t) + w² E +
tr ()) + L
Sha 2n + h-t EΩ,, + u h
( + h 12?
m
which completes the first claim.
When the is Le-Hessian smooth, by D.7 ,, we have
E (t) Xx (t)() x (t)(t) vn vz((t + vi.l?
+ u21² (RL2EER
12LM
|x= -22
L3
m
Eh"EEEr + 1 E + u2 1.T.MMM
h 1=1|x 1|2 + 12
m
which establishes the second claim.
D.3
Bound of B| Banip' (V (x, (8)) (((x())))
Lemma D.5. Let (f) and x2 (f) be defined as above and 1+m Then,
When Assumption 3.1 holds, and H be defined as Definition 3.2 Then
L
BB E Baop (VS(x, - ((((. > 1 + -
+ + 11,
iii When Assumptions 3.1 F.I and E.2 hold, we have
E8 El|E Ea-g' (x (8)) - (a))
L M + 1 1 10
LEE, +
Proof. Note that x,, (t) depends on a. First, we have
E
E Ea Earp V ((((x, - (a)))
E Eanp ESS
[ (3)) Eanp (a))
Baw of Eanp'
(V((x(()) - (+)((
< E.Eap Egn
V (x(x((s)) - (s)((
= -Bang Eyn E
| (S(6
26

---

Define 71" Lor + (( (,)))) and recall that H' - So (xx,(a)))) + k(x,
{a(a)))dkk And then E [|V]f(x(()) (s)]1 a, can be bounded by
V (x,, (s)) (s)
|| (x (s) - [x))))) a.
UFH / A12(t - s) - (x(x))) a,
E
Su Su² 1.2h41,m -V...)la))
where _ is by plugging in the close-form solution (s) and (x); S follows by A12(1 - 8) S and the Lipschitz continuous of V..
pi (f(xx,)))) - (a))
-
SuZn" 122,480, +
L Ftr 113 1. |x,11 w1/78 / A12(1 A o)F''B,
m
h6
Ett( h 1, + L 1.00 A12(t- - a) 1?
L, +
m
where X follows by
D.4 When 2 has uniform upper bound H, have
E B| Banp' V (3) VV
21 RR,,
+ 1, uh tr(?)
L m
< I,
Ir Ir(H) 1
where in the last inequality, we use that ( - LEtr (F') s Ltr(H) and tr(HH) s Ltr(H).
And when Assumptions and E.2 hold, by Lemma D.3
E Banp Hanp (x (.
I1 + h Eir(H? + 11 x. x2. uh IM + 13321
m
L M + xx 19
21 + m
E Local Hessian bound results
As discussed in Section 4.1, here show that the uniform Hessian upper bound can be relaxed to local trace control with an additional Hessian smooth assumption. Specifically, we make the following assumptions on
Assumption E.1. has a Le-Lipschitz Hessian.
Assumption E. a standard assumption that frequently appears the literature of optimization and sampling Ma et al. 2021 Dalalyan and Riou-Durand 2020 Durmus and Moulines 2019).
Assumption E.2. There exist M > such that for all x E R, tr ( f(x) < M.
Assumption E.2 relax the uniform Hessian bound V² [x) < for x E Rd under which we obtain convergence rate with tr(H). Assumption E.2 only requires that the trace of Hessian can be controlled locally. Under Assumptions E.T and 2 the convergence of DRUL only has a weak
dimension dependency, stated in Theorem E.3
27

---

Theorem E.3. For any tolerance E (0, 1), denote the minimizer of (xx by x.. Assume that
Assumptions 3.1 and E.2 hold. Set the step size
2/3
2/7
h S min 12C2n' (48C2k 1 M [1x, 1/3' (28C212uld2)) 1/7 m
where C2 2 1 is universal constant.
With initial point (xo. vol.
define Go when
Ex-p.v~((0..) (1xx - x ||xo + Vo
V x1(2). -xll Then under Assumptions 3.1
n 2 SEA
2ES
log
h
e²
Algorithm Cariimmoouttuuuut outputs x,, such that W2((La(Xn): p) S
Without loss of generality, assume x. = 0. Then Theorem E.3 indicates convergence rate of Ô (M1 + 2 which has lower dimension dependency both in the sense of the order of d and the relevant coefficient. When the Hessian of any E Rd has dimension-free trace, we arrive at convergence rate of O (E-2/3 + 6 2/T).
F. Proof of Theorem
In this section, prove the convergence guarantee using the upper bound of the local Hessian trace.
Proof. We follow the notations in the proof of Theorem in Appendix B.2 Similarly, denote the error in (3.5 by
A = 2uhC ELE (3) -
+ +x) - (B), V.f(xx - (a),
Then
A <2C]uhee (8) - Xn - x2, + + v. Vn - v, + -/5
+ 2C, uhlE (x, x, + Vn V Bang' Esnup (Vf(x, (a())
E x,, (3) Xn (8) + x + Vn (B) Vn (3) + vill? wh
wh
+ " 11 E
(x (a)))1 LuhBlix, - x vn v
whlB Larp B3 (x (3)) - (a(xl1)))
he Pn MED, u² h LM + + + 1o
m
+ u2 11212ER + 12L M + 1,21,11x.) 112 + 13221002) + 2
m
+ 2 11 + L M + 1/ 11 1| 12322100 d2
1n + + EE, + 2 L
M + 11 x 2| + L3
where 2 follows by Lemmas the above control with B.5 we have
D.I D.4 D.5
We shall note that E8 (3) Combining
Egnpe + Czh
C2h1
Iu2L
El
EΩ,
M + ||xx,12 2 + L3u h4d2
C2h2
1 Ozh2
m
h
L M + 2 133 02
1
EΩ2 + 2C2h2
6k
m
28

---

where the last incquality follows by the control (88)) and the choice of h which satisfies 1 1- C2h2
h and E¹ = 2C2hh [ M + x 112 + L2 (d/d)) and therefore
Denote
E
EΩN S -E + E r. ESo + E(1 ++ +r S
EΩo + --r
When N 2
6 log 2852 we have
2²
No BB S en IES20 S
(E.1)
2
And when
2/3
2/7
S
min
u (48C2x
M +
|xx. [12))
1/3'
(48C2L3u*d))
m
we get I-r
2 Together with the (E. we obtain the convergence rate Theorem
E.3
Finally, we would say that the trace of Hessian should be a more realistic factor to characterize the complexity of Langevin algorithms. Beyond the ridge separable functions described in Section 3.3 we illustrate the Bayesian sampling over the two-layer neural network model and show that the neural network also has bounded trace under suitable conditions. Note that our convergence analysis cannot be applied to this task since neural networks are non-convex models, and we leave the extensions of analysis to the general non-log-concave case as future work.
T x). where a is the activation function. When ||x|| S FI, ||| < and (x))) S a, we have ur (V2 ((,, w)) S ar1r2.
Proposition E.4. Define f(W, w) = w (
Proof. By direct computation, have
0ʃ
0 (d'(W (d'(Wx)) @ Øx,
a(W x),
Ow
oW
0 f
0 Diaglo" "( x) x)
Ow2 0,
OW2
Therefore,
tr 2 f(W, w))) = ||x| tr (Diag(o"(W) x) (WTx) @))
< s (WW x)
(E.2)
< 00122
29


Reviews:
[{"review_id": "fUL13zWKaw", "reviewer": "Reviewer_TP8b", "summary": "In this paper, the authors propose a Langevin-type algorithm for sampling a strongly log-concave distribution with a composite structure. Their method can be viewed as a variant of the randomized midpoint method, with two key modifications: (i) they only discretize the smooth convex part of the negative log likelihood but retain the strongly convex part; (ii) they draw both step sizes in the algorithm randomly according to carefully crafted distributions. It is shown that the algorithm achieves an accelerated rate without explicit dependence on the dimension.", "questions": "- While this work focuses on the dimension dependence, the condition number $\\kappa = L/\\mu$ can also impact the convergence rate greatly. How is the result in this paper compared with the existing works in terms of the dependence on $\\kappa$?\n- In Lemma 5.1, it is unclear to me what it means that \"$x_n$ and $x_n^*$ are coupled synchronously\". Do you mean that they are driven by the same Brownian process?\n- Page 9, Lines 277-281: I am confused by this paragraph. By definition, isn't the random weight $w_n(s,\\alpha)$ a $d$-dimension vector? If so, how can you apply Claim (A) in Lemma 5.2?\n- Page 9, Lines 282-292: It is also unclear to me why \"the randomized step size make it possible to consider the averaged effect\".\n- Page 12, (A.1): I am not sure why the authors introduce the extra parameter $\\gamma$. As far as I can see, $\\gamma$ is fixed as 2 in the rest of the proof.  \n- Page 14, Section B.1.1: I don't see why the random process $B_t^{\\alpha}$ is a Brownian bridge. Is it supposed to be the random process $B_t$ conditioned on $B_{\\alpha}$? And why do we need to introduce this process in the first place?\n- Page 14, the equations under Line 431: Here the authors exchange the order of differentiation and expectation, which is not justified. Indeed, it is not even clear if $\\Omega(t)$ is differentiable since it involves the solution trajectories of SDEs. \n\n\nTypos in the paper:\n\n- The convergence rates reported in the introduction are inconsistent with the ones in Table 1. Specifically, the rate by [Shen and Lee, 2019] should be $\\tilde{O}(\\frac{d^{1/3}}{\\epsilon^{2/3}})$ (Page 2, Line 47), and the rate by [Freund et al., 2022] should scale linearly with $O(\\frac{1}{\\epsilon})$ (Table 1, Row 5). \n- Page 1, Line 49: \"convergence dependence\" -> \"dimension dependence\"\n- Page 5, Definition 3.6: $A$ is undefined. \n- Page 9, Line 275: the integral should be $\\int_{0}^t e^{\\frac{s-t}{\\kappa}} F(s) ds$.", "limitations": "Yes, the authors addressed the limitations of their work. ", "rating": 6, "confidence": 2, "soundness": 3, "presentation": 2, "contribution": 3, "strengths": "- The result is interesting and noteworthy. To sample a strongly log-concave distribution, the best-known iteration complexity bound either has an undesirable dimension dependence, such as $\\tilde{O}(\\frac{d^{1/3}}{\\epsilon^{2/3}})$ in [Shen and Lee, 2019], or a worse dependence on $\\epsilon$, such as $\\tilde{O}(\\frac{\\mathrm{tr}(H)}{\\epsilon})$ in [Freund et al., 2022]. In this work, the authors manage to achieve the best of both worlds and prove a complexity bound of $\\tilde{O}(\\frac{(\\mathrm{tr}(H))^{1/3}}{\\epsilon^{2/3}})$. \n- The authors introduce the double randomized technique to reduce the discretization error, which seems novel to me. ", "weaknesses": "I think the presentation of the paper can be improved.\n\n- In particular, the explanation in Section 5 is not very helpful, and it remains unclear to me why the randomized step size helps reduce the discretization error, and why the authors choose the specific distribution in Lemma 5.2. Moreover, it would be helpful if the authors can compare their analysis with the one in [Shen and Lee, 2019] to better explain how they remove the dimension dependence.  \n\n-  Also, there are numerous typos in the main text and the proofs in the appendix, which sometimes make it hard to understand. Please see the \"Questions\" section for more details. ", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": " Page 9, Lines 282-292: It is also unclear to me why \"the randomized step size makes it possible to consider the averaged effect\".\n● In particular, the explanation in Section 5 is not very helpful [...] better explain how they remove the dimension dependence.\n\nIn short, the introduction of the random step size is to bound the error such as $\\|\\int_0^t x_n(s) - x_n^*(t)\\mathrm{d} s\\|^2$ in the discretization analysis. The $\\|\\int_0^t x_n(s) - x_n^*(s)\\mathrm{d} s\\|^2$ may be dimension-dependent as a standard bound will lead to a $dt$ additional discretization error. We seek to bound this kind of error using the descent on $\\Omega_n(t)$. Under the stochastic step size, $\\Omega_{n+1}$ will be $ \\mathbb{E} \\Omega_n(t) $,  and thus we can bound the dimension-dependent discretization error as $\\|\\int_0^t x_n(s) - x_n^*(s)\\mathrm{d} s\\|^2 \\lesssim \\mathbb{E}\\Omega_{n}(t)$. And it remains a careful analysis to obtain the dimension-free error and delicate device of $\\rho, \\rho'$.\n\n\n● While this work focuses on dimension dependence, the condition number $\\kappa=L/μ$ can also impact the convergence rate greatly. How is the result in this paper compared with the existing works in terms of the dependence on $\\kappa$?\n\nThanks for your suggestions. This paper mainly considers dimension dependence. It is very interesting to study the dependence on $\\kappa$. However, we currently do not take a very careful analysis of it. We have achieved a $\\kappa^{4/3}$ dependency. We think it is possible to improve this dependence using the technique in Cao et al. (2019) (On explicit $ L^ 2$-convergence rate estimate for underdamped Langevin dynamics). And we will work on it in the near future.  \n\n● In Lemma 5.1, it is unclear to me what it means that \"$x_n$ and $x_n^∗$ are coupled synchronously\". Do you mean that they are driven by the same Brownian process?\nIt exactly means that they are driven by the same Brownian motion.\n\n● Page 9, Lines 277-281: I am confused by this paragraph. By definition, isn't the random weight $w_n(s,\\alpha)$ a $d$-dimension vector? If so, how can you apply Claim (A) in Lemma 5.2?\n\nThanks for pointing out the ambiguity. Line 277 illustrates how our choice of $\\rho$ and $\\rho’$ inherits the midpoint discretization. And Claim (A) shows that the match of expectation is critical for obtaining this acceleration. In rigorous proof, we apply (A) of Lemma 5.2 with $F(s) = \\langle w_n(0, \\alpha), \\nabla f(x_n(s)) - \\nabla f(\\hat{x_n}(\\alpha)) \\rangle$. And this analysis neglects a discretization error on $\\bar{w}$ and $w_n(s, alpha)$ as discussed in line 280. And it is supposed to be an inner product in line 279. We will make it clear in a further revised version.\n\n\n\n● Page 12, (A.1): I am not sure why the authors introduce the extra parameter $\\gamma$. As far as I can see, $\\gamma$ is fixed as $2$ in the rest of the proof.\n\nThanks for raising this issue with us. We indeed fix $\\gamma = 2$ in the rest of the proof. We will specify that $gamma = 2$ in a further revised version.\n\n● Page 14, Section B.1.1: I don't see why the random process $B_t^{\\alpha}$ is a Brownian bridge. Is it supposed to be the random process $B_t$ conditioned on $B_{\\alpha}$? And why do we need to introduce this process in the first place?\n\n$B_t^s$ is a Browning ridge which is conditioned on $B_{\\alpha}$. And the process is the SDE of the Brownian bridge. We introduce $B_t^s$ for technical considerations. The diffusion process in line 428 is an adapted process conditioned on $B_{\\alpha}$. And process $(x_n(t), v_n(t))$ driven by $B_t^{\\alpha}$ is not adapted to the filtration when $t < \\alpha$.\n\n\n● Page 14, the equations under Line 431: Here the authors exchange the order of differentiation and expectation, which is not justified. Indeed, it is not even clear if $\\Omega(t)$ is differentiable since it involves the solution trajectories of SDEs.\n\nThanks for raising this to us. The exchange of differential and expectation requires further validation. And a rigorous proof will be considering the evolution of the process not in expectation but only $\\mathrm{d} \\Omega(t)/\\mathrm{d} t$ and then taking the expectation in line 437. And in line 437, we take the differential of $\\Omega(t)$ using Ito’s formula.\n\n\n● The convergence rates reported in the introduction are inconsistent with the ones in Table 1. Specifically, the rate by [Shen and Lee, 2019] should be $\\tilde{O}\\left( \\frac{d^{1/3}}{\\epsilon^{2/3}} \\right)$ (Page 2, Line 47).\n\nThank you for pointing out this issue. It is indeed $\\epsilon^{2/3}$ in the denominator. \n\n● And the rate by [Freund et al., 2022] should scale linearly with $O\\left( \\frac{1}{\\epsilon} \\right)$ (Table 1, Row 5).\n\nWe compare the convergence rate in Wasserstein distance whereas the $\\frac{1}{\\epsilon}$ rate is proven in KL distance, which implies a $\\frac{1}{\\epsilon^2}$ rate in Wasserstein distance.\n"}, {"review_id": "t2LJknG8Tu", "reviewer": "Reviewer_SZ7b", "summary": "The paper adapts the Randomized Midpoint Method to the composite optimization context considered in Freund et al, and consequently improves the dependence from $O(tr(H)/\\epsilon)$ to $O((tr(H)/\\epsilon)^{1/3})$.", "questions": "It seems inappropriate to compare this to the original randomized midpoint work/other work for the standard Langevin Monte Carlo, since they do not assume the composite structure of the problem. The primary highlighted comparison is with respect to Freund et al. (2022), in which case the gain is more like $tr(H)^{1/3}/\\epsilon^{4/3}$. If $tr(H)$ is $O(1)$ then this is only a gain in epsilon, which is usually smaller than $d$.\n\nWhat is the previous proof referred to in L. 253?\n\nTerminology of “acceleration” should probably be avoided since this is classically used only to refer to sqrt(kappa) rates. A better term might be “improved discretization error”.\n\nTypos:\n\nL. 140 What is being made strongly convex?\n\nL. 185 squred -> squared\n\nL. 197 denote solution -> denote the solution\n", "limitations": "I have outlined my concerns already in the previous sections.", "rating": 5, "confidence": 4, "soundness": 4, "presentation": 4, "contribution": 2, "strengths": "The application of randomized midpoint in this composite sampling is novel, and there is genuine improvement in the rate estimate when compared to Freund et al.\n\nThe technique of double randomization is new and requires some novel analysis when contrasted with prior works.\n\nThe authors do a decent job of illustrating that the trace of the Hessian is $o(d)$ through some figures and discussion.\n", "weaknesses": "The primary contributions of this paper are not particularly original and mostly stem from combining the framework in Freund et al. with the known analysis for randomized midpoint in Shen and Lee. This in my view is the primary weakness of the paper.\n\nIn general, claims about the “dimension-free” nature of the convergence guarantees need to be careful since the composite structure assumption is quite strong, although the authors in general do a good job of qualifying their claims.\n\nOverall, while this work contains some novel claims and results, and is a bona fide improvement on prior work. However, the technical novelty is not significant, and I am borderline on this paper as a result.\n", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "● The primary contributions of this paper are not particularly original and mostly stem from combining the framework in Freund et al. with the known analysis for the randomized midpoint in Shen and Lee. This in my view is the primary weakness of the paper.\n\nWe would like to emphasize that our method is not a simple combination of the previous work.\n\nFirst, the two-stage optimization in Freund et al. cannot be directly extended to the accelerated process or complicated discretization scheme due to its form. Meanwhile, the midpoint discretization analysis introduces dimension-dependent errors such as in line 447, which does not exist in the analysis in the EU discretization of the overdamped Langevin analysis and underdamped analysis. This requires refined discretization schemes to eliminate the dimension-dependent error and thus there are underlying difficulties to extend to complicated discretization schemes. We design a discretization scheme and prove the scheme can achieve a convergence rate with $\\mathrm{tr}(H)$ dependence. And this can be achieved for methods based on both accelerated processes and complicated discretizations. This implies that a wide range of Langevin algorithms can be adapted to achieve low dimension dependence. \n\n\n● In general, claims about the “dimension-free” nature of the convergence guarantees need to be careful since the composite structure assumption is quite strong, although the authors in general do a good job of qualifying their claims.\n\nSorry for the ambiguity of our statement. We would like to clarify that our analysis is not restricted to composite structure and applies to general strongly convex target distributions. As discussed in section 3.2, the composite structure also includes the general $m$-strongly convex function, which can be divided into $g(x) = \\frac{m}{2} \\|x\\|^2$ and the weakly convex function $f (x) = U (x) −\\frac{m}{2} \\|x\\|^2$. This leads to the same $\\frac{\\mathrm{tr}(H)^{1/3}}{\\epsilon^{2/3}}$ convergence rate. One of the purposes of introducing composite structure in the proof is to alleviate the dimension dependence in the Hessian trace. The composite structure does not serve as the prerequisite of the analysis. Otherwise, the $\\mathrm{tr}(H)$ will have an extra $md$ term. The scaling of $md$ is problem specific and may be dimension dependent when $\\frac{1}{m} = o(d)$. The natural composite structure is ubiquitous in sampling tasks. \n\n● It seems inappropriate to compare this to the original randomized midpoint work/other work for the standard Langevin Monte Carlo, since they do not assume the composite structure of the problem. \n\n● L. 140 What is being made strongly convex?\n\nAs discussed above, the algorithm applies to the general strongly convex distributions. $U$ is strongly convex. We will discuss the composite structure more in a further revised version.\n\n● What is the previous proof referred to in L. 253?\n\nThe proof in Shen and Lee (2019). We follow the seminal contraction-based method by Shen and Lee (2019) and Cheng et al. (2018). This is to illustrate the difference of the proof to one familiar with the proof by Shen and Lee. Our analysis tracks the difference between the implemented process and the exact process with the initial distribution corresponding target distribution, and considers to bound the local error in the process. It is consistent with the proof based on the gradient flow.\n\n\n● Terminology of “acceleration” should probably be avoided since this is classically used only to refer to sqrt(kappa) rates. \n\nThanks for pointing out this issue. Acceleration is indeed improper. \n"}, {"review_id": "O0qGFA4zDk", "reviewer": "Reviewer_5mo9", "summary": "The paper suggests a novel version of the Unadjusted Langevin Algorithm with sample complexity in Wasserstein-2 distance scaling with the effective dimension of the problem (trace of the potential's Hessian) instead of the ambient space dimension in case of strongly log-concave distributions. This result completes and generalizes the results of [Shen and Lee, 2019].", "questions": "1. Is it possible to add any numerical findings illustrating the superiority of the doubly randomized ULA (with random step size) against the one with constant or decreasing step size? If one could trace the precise dependence upon the $\\operatorname{trace}{H}$ even in the toyish setup, I would lean towards increasing my score.\n \n2. Are there any novel technical contributions developed to prove the result of Theorem 4.2? If yes, please add the corresponding discussion to the main text.", "limitations": "The paper is theoretical and no negative societal impact is expected.", "rating": 6, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "The research direction towards studying the sample complexity rates in terms of the effective dimension is interesting and potentially allows to explain the successful behaviour of the ULA-type algorithms in the high-dimensional problems.", "weaknesses": "First of all, the suggested Algorithm 1 (DRUL) does not seem to be really an implementable one, since the discretization error in $x_{n+1}$ and $v_{n+1}$ is to appear in the practical implementations. It is not clear if the control of this discretization error would not yield an explicit dimension dependence in the stepsize $h$ in Theorem 4.2. \n\nSecond, the sample complexity scaling as $\\varepsilon^{-2/3}$ is not completely convincing. For example, for the ridge separable potentials, which are the main motivating example towards the paper, the Hamiltonian Monte-Carlo method is known to obtain a sample complexity of order $(d/\\varepsilon)^{1/4}$, see e.g. [Mangoubi et al, 2017]. Thus there is a natural question if the $\\varepsilon^{-2/3}$ complexity optimal for Langevin-type algorithms? Again it seems that the particular rate can degrade after the intergral discretization in Alg. 1 taken into account.\n\nReferences:\nMangoubi, O., & Smith, A. (2017). Rapid mixing of Hamiltonian Monte Carlo on strongly log-concave distributions. arXiv preprint arXiv:1708.07114.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "● First of all, the suggested Algorithm 1 (DRUL) does not seem to be really an implementable one, […] It is not clear if the control of this discretization error would not yield an explicit dimension dependence in the stepsize $h$ in Theorem 4.2. \n\nWe apologize for the possible misunderstanding in the presentation of ALgorithm 1. We would like to clarify that Algorithm 1 is implementable and there is no dimension dependence error or necessity of further discretization. The discretized process (4.1) is linear in $x_n(t)$ and $v_n(t)$ and thus closed form solution is given by a Gaussian distribution. The exact solution is given by (4.3). Also note that the integral in Algorithm 1 is tractable, e.g. $\\int_0^{\\alpha_n} A_{12}(s-\\alpha_n)\\nabla f(x_n)\\mathrm{d}s = \\int_0^{\\alpha_n} A_{12}(s-\\alpha_n)\\mathrm{d}s \\nabla f(x_n)$ and $A_{12}$ defined in Lemma 4.1 has closed-from integral. To implement Algorithm 1, one only needs to sample from a Gaussian distribution whose mean is tractable and the variance matrix is given by (A.3). We will specify the solution of the integral in Algorithm 1 and the covariance in (A.3) in a later revised version. \n\n● Second, the sample complexity scaling as $\\epsilon^{-2/3}$ is not completely convincing. […] can degrade after the integral discretization in Alg. 1 is taken into account.\n\nWe would like to clarify that the ridge separable case is an illustration that $\\mathrm{tr}(H)$ can be dimension independent. The analysis is not restricted to the ridge separable case and reduces the dimension dependence for a wide range of problems that have rapidly dropping Hessian eigenvalues. One can refer to our discussion at the end of section 4.1.\n\nAnd Mangoubi et al (2017) do not assume the ridge separable structure. They assume in assumption 1.7 that the potential $U(x)$ can be separated according to the blocked state space, that is $U(x) = \\sum_{i=1}^{d/m}U(x_{i})$ where $x_i\\in\\mathbb{R}^{m}$. Although Mangoubi et al (2017) prove that HMC achieves $d^{1/4}$ or $\\epsilon^{1/4}$ rate, the convergence guarantees require higher-order discretization schemes. These higher-order ODE solvers require higher-order smoothness which is not necessary in our algorithm. Such as the leap-frog scheme requires third-order smoothness to obtain the guarantee.\n\n$\\left(\\frac{d}{\\epsilon}\\right)^4$ can be achieved via higher-order Langevin schemes by Mou et al. (2021) (High-order Langevin diffusion yields an accelerated MCMC algorithm) for ridge separable case. However, it is based on a different oracle $\\Delta U$ other than the gradient oracle considered in our algorithms. And the convergence has a much worse dependence on conditioning number. The difference oracle $\\Delta U$ is usually intractable. For general target distributions, they require numerical integration and the convergence rate degrades.\n\nAs for the optimality, the discretization scheme with convergence rate $\\epsilon^{-2/3}$ is known to be optimal for underdamped Langevin algorithms [Cao et al. (2020)] (Complexity of randomized algorithms for underdamped Langevin dynamics). It is a discretization lower bound.\n\n● Is it possible to add any numerical findings illustrating the superiority of the doubly randomized ULA (with random step size) against the one with constant or decreasing step size?\n\nWe have conducted some synthetic studies to show that our DRUL (with random step size) can achieve better mixing distribution compared to the randomized midpoint with a constant step size. Here we illustrate the dimension dependence using synthetic data. We consider Bayesian ridge regression with $g(w) = 0.1\\|w\\|^2/2$ and $f(w) =\\frac{1}{2n} \\|Xw - Y\\|^2$, where rows of $X$ have unit norm. $f(x)$ has dimension-independent trace. We evaluate the result by a projection to one dimension and calculate the empirical Wasserstein distance (since it is non-trivial to obtain precise high dimensional distance using empirical data). We fix the step size and compare the mixing result. We take the average of 3 independent evaluations. Since the algorithms achieve high precision, we require a large number of samples (500000) to get an accurate estimate of Wasserstein distance even in one dimension. It is a huge computation cost in high-dimension regime. Thus we fixed a large step $h=1$ (under which the mixing step is low) size and $T = 300$ steps to guarantee mixing. Below is the results, where we compare the result of the randomized midpoint method (RMM) and our DRUL method when the dimension $d = 5,10,100,1000$.\n\n| Dimension | 5      | 10     | 100    | 1000   |\n|-----------|--------|--------|--------|--------|\n| RMM       | 0.0084 | 0.0083 | 0.0180 | 0.0416 |\n| DRUL      | 0.0085 | 0.0071 | 0.0050 | 0.0077 |\n\n\nFrom the table, DURL achieves more accurate mixing distributions in high-dimension space, thus showing superiority as the dimension grows. Note that by directly sampling from the target distribution 500000 times, the empirically estimated distance from the target distribution is around 0.001-0.005. DRUL method may achieve better results than the ones shown in the table. \n\nAs a reminder, our work is not a stochastic step size Midpoint method. The stochastic step size Midpoint method will induce a differential $\\rho’$ that does not satisfies claim (B) in Lemma 5.2. Besides, the decreasing stepsize also applies to our method.\n\n● Are there any novel technical contributions developed to prove the result of Theorem 4.2? If yes, please add the corresponding discussion to the main text.\n\nWe develop a new discretization scheme and bound the error-dependent term by an averaged effect as discussed in Section 5. The analysis requires delicate control and devising on the choice of $\\rho$ and $\\rho'$. And our analysis bounds the local error in the process, which leads to our choice of $\\rho$ and $\\rho'$. We will make a more clear and thorough discussion in a revised version.\n"}, {"review_id": "cfNlJeuqTk", "reviewer": "Reviewer_RFtN", "summary": "This paper considers the problem of sampling from a Gibbs distribution $p(x) \\propto e^{-U(x)}$ using discretized Langevin dynamics. Since the approximation error of such methods usually depends on $\\mathrm{Tr}(\\nabla^2 U)$, the proposed algorithm splits $U$ into a quadratic part $g(x) = \\frac m2 ||x||^2$ and a remainder $f$, and only discretizes the dynamics according to $f$.\n\nSuch a split was already considered in [Freund et al. '21]; however, the main novelty of this article is that the discretization on $f$ is performed through a two-step method with random step-sizes, instead of a simple gradient update. This scheme shaves off a factor of $\\epsilon^{-1/3}$ from the required time complexity to reach a Wasserstein error of $\\epsilon$.\n\nThe proof is based on an approximate contraction bound for a quantity $\\Omega_n$ that bounds the desired Wasserstein distance, followed by a precise analysis and optimization of the error in the aforementioned bound. ", "questions": "- How important is the choice of $\\rho$ ? I understand from Lemma 5.2 and its implications that for a given $\\rho$, the choice of $\\rho'$ is important to ensure these conditions, but why can't I, for example, choose a uniform prior for $\\beta$?\n- In section 5.2, it seems like you take $\\bar w_n = x_n - x_n^* + v_n - v_n^*$ as the expectation of $w_n(s, \\alpha)$; why is it the case ?", "limitations": "N/A", "rating": 7, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "This paper presents a novel algorithm for Langevin simulation, that achieves state-of-the-art performance for the dependency on both the precision requirement $\\epsilon$ and the ambient dimension $d$. The proposed algorithm is fairly simple (at least for quadratic $g$) and easy to implement, and the main ideas behind it are clearly explained. Overall, the paper is fairly well-written, with only a few typos here and there.", "weaknesses": "Clarity/soundness: some of the proofs are very hard to parse due to the amount of simplifications made at once from one line to the next. This is especially felt in B.2, where the first two inequalities below l.447 contain around 5 sub-inequalities to check each. The specifications $\\alpha \\sim \\rho'$ and $\\beta \\sim \\rho$ are also used inconsistently, which makes it hard to understand with respect to what quantities each expectation is taken.\n\nNovelty/significance: The relationship to [Shen and Lee '19] and [Freund et al. '22] should probably be expanded. From what I understand, this paper unites the randomized midpoint method of the former with the combined optimization viewpoint of the latter, but it is unclear if there are challenges other than computational to this endeavor. Namely, neither of those methods require such a complicated step-size scheme, which seems to be the main novelty of the paper, but the need for it is unclear.\n\nMinor remarks:\n- eq. (3.5): what is A?\n- l.185: \"squre\"\n- the equation below l.279 should be a scalar product.\n- in (B.4), shouldn't $z_n(t)$ be $\\hat x_n(\\alpha) - x_n^*(t)$ instead ? This change does ripple through the proof of Lemma 5.1, so I'm not actually sure of how minor it is.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "● Neither of those methods require such a complicated step-size scheme, which seems to be the main novelty of the paper, but the need for it is unclear.\n\n\nWe introduce the random step size to bound the error such as $\\|\\int_0^t x_n(s) - x_n^*(t)\\mathrm{d} s\\|^2$ in the discretization analysis. The $\\|\\int_0^t x_n(s) - x_n^*(t)\\mathrm{d} s\\|^2$ may be dimension-dependent as a standard bound will lead to a $dt$ additional discretization error. We seek to bound this kind of error using the descent on $\\Omega_n(t)$. Under the stochastic step size, $\\Omega_{n+1}$ will be $\\mathbb{E} \\Omega_n(t)$, and thus we can bound the dimension-dependent discretization error as $\\|\\int_0^t x_n(s) - x_n^*(t)\\mathrm{d} s\\|^2\\lesssim \\mathbb{E}\\Omega_{n}(t)$. And then there reminds additional conditions on $\\rho$ and $\\rho'$ to \n\n● The specifications $\\alpha\\sim\\rho'$ and $\\beta\\sim\\rho$ are also used inconsistently, which makes it hard to understand with respect to what quantities each expectation is taken.\n\nThanks for pointing out the clearness issue in the presentation. We frequently use Claim (B) in our proof which indicates that we can replace the expectation under $\\rho’$ with the one under $\\rho$ up to a constant level. Mostly, its purpose is to control the error term with the expectation of $\\rho'$ using the contraction which is an expected value over $\\rho$. And we use Claim (A) in Lemma D.3. We will make more explanations for our proof especially on the specification of the $\\rho$ and $\\rho'$ in a further revised version.\n\n\n● How important is the choice of $\\rho$? \n\nWe believe that $\\rho$ is not restrictive to the proposed one. However, this choice of $\\rho$ has to satisfy several conditions. Our choice of $\\rho$ induces a simple $\\rho'$ and is a natural choice. We require its expectation is of $\\Theta(h)$ order which guarantees sufficient descent and the average descent can control the dimension-dependent error as in (D.18). And the corresponding $\\rho’$ should satisfy claim (B) in Lemma 5.2. Uniform $\\rho$ is possible given that corresponding $\\rho’$ is chosen accordingly by Lemma 5.2. But this requires further validations and induces a more complicated $\\rho'$. As a reminder, uniform $\\rho$ does not lead to the random step size version of the randomized midpoint method. The distribution of $\\rho’$ of the stochastic step size randomized midpoint method is different and it does not satisfy claim (B) in Lemma 5.2.\n\n\n\n● in (B.4), shouldn’t $z_n(t)$ be $\\hat{x}_n(\\alpha) - x_n^*(t)$.\n\nWe appreciate the careful revision and pointing out this issue. This is our typo and does not influence the following analysis. This typo can be fixed only by changing the equation in line 430 and line 431. $z_n(t)$ is $x_n(t) - x_n^*(t)$. But Equation (B.4) should be $\\nabla U(x_n(t)) - U(x_n^*(t)) = \\int_0^1 \\nabla^2 U(s\\hat{x}_n(\\alpha) + (1-s)x_n^*(t)) \\mathrm{d} s z_n(t)$. And the following analysis remains unchanged. The last inequality follows by dividing $\\nabla f(\\hat{x}_n(\\alpha)) - \\nabla f(x_n^*(t))$ into $\\nabla f(\\hat{x}_n(\\alpha)) - \\nabla f(x_n(t))$ and $ f(\\hat{x}_n(\\alpha)) - \\nabla f(x_n^*(t))$.\n\n\n● In section 5.2, it seems like you take $\\bar{w}_n= x_n - x_n^* + v_n - v_n^*$ as the expectation of $w_n(s,\\alpha)$.\n\nThanks for pointing out the unclearness in our presentation. We do not take $\\bar{w}$ as the expectation of $w_n(s,\\alpha)$. There is a higher-order discretization error gap between the approximation $\\bar{w}$ and $w_n(s,\\alpha)$. Here we illustrate that claim (A) of 5.2 leads to a randomized mid-point discretization. One can refer to the first equation in line 447 for the rigorous proof. And we will specify the error term and our choice of $\\bar{w}$ in the main text in a later refined version.\n\n\n"}]

Author rebuttal:

Meta-review: The paper proposes a Langevin Monte Carlo variant to sample from a strongly log-concave target distribution with a composite structure. Authors are inspired by Randomized Midpoint Method and establish sample complexity in Wasserstein-2 distance. Authors achieve this by only discretizing the smooth convex part of the negative log likelihood. Resulting algorithm enjoys faster rates.

The paper is reviewed by 4 reviewers with the following Rating/Confidence scores: 5/4, 7/3, 6/4, 6/2. The methodology is an interesting extension of Randomized Midpoint method and improvements are significant. I recommend including this paper to conference proceedings.
Decision: Accept (poster)