SQRT(D) DIMENSION DEPENDENCE OF LANGEVIN MONTE CARLO
Year: 2022

Abstract:
This article considers the popular MCMC method of unadjusted Langevin Monte Carlo (LMC) and provides a non-asymptotic analysis of its sampling error in 2Wasserstein distance. The proof is based on a refinement of mean-square analysis in Li et al. (2019), and this refined framework automates the analysis of a large class of sampling algorithms based on discretizations of contractive SDEs. Using this framework, we establish an Õ (√ d/ ) mixing time bound for LMC, without warm start, under the common log-smooth and log-strongly-convex conditions, plus a growth condition on the 3rd-order derivative of the potential of target measures. This bound improves the best previously known Õ ( d/ ) result and is optimal (in terms of order) in both dimension d and accuracy tolerance for target measures satisfying the aforementioned assumptions. Our theoretical analysis is further validated by numerical experiments.

text: (√ d/ )
mixing time bound for LMC, without warm start, under the common log-smooth and log-strongly-convex conditions, plus a growth condition on the 3rd-order derivative of the potential of target measures. This bound improves the best previously known Õ ( d/ )
result and is optimal (in terms of order) in both dimension d and accuracy tolerance for target measures satisfying the aforementioned assumptions. Our theoretical analysis is further validated by numerical experiments.

heading: 1 INTRODUCTION
text: The problem of sampling statistical distributions has attracted considerable attention, not only in the fields of statistics and scientific computing, but also in machine learning (Robert & Casella, 2013; Andrieu et al., 2003; Liu, 2008); for example, how various sampling algorithms scale with the dimension of the target distribution is a popular recent topic in statistical deep learning (see discussions below for references). For samplers that can be viewed as discretizations of SDEs, the idea is to use an ergodic SDE whose equilibrium distribution agrees with the target distribution, and employ an appropriate numerical algorithm that discretizes (the time of) the SDE. The iterates of the numerical algorithm will approximately follow the target distribution when converged, and can be used for various downstream applications such as Bayesian inference and inverse problem (Dashti & Stuart, 2017). One notable example is the Langevin Monte Carlo algorithm (LMC), which corresponds to an Euler-Maruyama discretization of the overdamped Langevin equation. Its study dated back to at least the 90s (Roberts et al., 1996) but keeps on leading to important discoveries, for example, on non-asymptotics and dimension dependence, which are relevant to machine learning (e.g., Dalalyan (2017a;b); Cheng et al. (2018a); Durmus et al. (2019); Durmus & Moulines (2019); Vempala & Wibisono (2019); Dalalyan & Riou-Durand (2020); Li et al. (2019); Erdogdu & Hosseinzadeh (2021); Mou et al. (2019); Lehec (2021)). LMC is closely related to SGD too (e.g., Mandt et al. (2017)). Many other examples exist, based on alternative SDEs and/or different discretizations (e.g., Dalalyan & Riou-Durand (2020); Ma et al. (2021); Mou et al. (2021); Li et al. (2020); Roberts & Rosenthal (1998); Chewi et al. (2021); Shen & Lee (2019)).
Quantitatively characterizing the non-asymptotic sampling error of numerical algorithms is usually critical for choosing the appropriate algorithm for a specific downstream application, for providing practical guidance on hyperparameter selection and experiment design, and for designing improved samplers. A powerful tool that dates back to (Jordan et al., 1998) is a paradigm of non-asymptotic
error analysis, namely to view sampling as optimization in probability space, and it led to many important recent results (e.g., Liu & Wang (2016); Dalalyan (2017a); Wibisono (2018); Zhang et al. (2018); Frogner & Poggio (2020); Chizat & Bach (2018); Chen et al. (2018); Ma et al. (2021); Erdogdu & Hosseinzadeh (2021)). It works by choosing an objective functional, typically some statistical distances/divergences, and showing that the law of the iterates of sampling algorithms converges in that objective functional. However, the choice of the objective functional often needs to be customized for different sampling algorithms. For example, KL divergence works for LMC (Cheng & Bartlett, 2018), but a carefully hand-crafted cross term needs to be added to KL divergence for analyzing KLMC (Ma et al., 2021). Even for the same underlying SDE, different discretization schemes exist and lead to different sampling algorithms, and the analyses of them had usually been case by case (e.g., Cheng et al. (2018b); Dalalyan & Riou-Durand (2020); Shen & Lee (2019)). Therefore, it would be a desirable complement to have a unified, general framework to study the non-asymptotic error of SDE-based sampling algorithms. Toward this goal, an alternative approach to analysis has recently started attracting attention, namely to resort to the numerical analysis of SDE integrators (e.g., Milstein & Tretyakov (2013); Kloeden & Platen (1992)) and quantitatively connect the integration error to the sampling error. One remarkable work in this direction is Li et al. (2019), which will be discussed in greater details later on.
The main tool of analysis in this paper will be a strengthened version (in specific aspects that will be clarified soon) of the result in Li et al. (2019). Although this analysis framework is rather general and applicable to a broad family of numerical methods that discretize contractive1 SDEs, the main innovation focuses on a specific sampling algorithm, namely LMC, which is widely used in practice. Its stochastic gradient version is implemented in common machine learning systems, such as Tensorflow (Abadi et al., 2016), and is the off-the-shelf algorithm for large scale Bayesian inference. With the ever-growing size of parameter space, the non-asymptotic error of LMC is of central theoretical and practical interest, in particular, its dependence on the dimension of the sample space. The best current known upper bound of the mixing time in 2-Wasserstein distance for LMC is Õ ( d
) (Durmus & Moulines, 2019). Motivated by a recent result (Chewi et al., 2021) that shows better dimension dependence for a Metropolis-Adjusted improvement of LMC, we will investigate if the current bound for (unadjusted) LMC is tight, and if not, what is the optimal dimension dependence.

heading: Our contributions
text: The main contribution of this work is an improved Õ
(√ d ) mixing time upper bound for LMC in
2-Wasserstein distance, under reasonable regularity assumptions. More specifically, we study LMC for sampling from a Gibbs distribution dµ ∝ exp ( −f(x) ) dx. Under the standard smoothness and strong-convexity assumptions, plus an additional linear growth condition on the third-order derivative of the potential (which also shares connections to popular assumptions in the frontier literature), our bound improves upon the previously best known Õ ( d ) result (Durmus & Moulines, 2019)
in terms of dimension dependence. For a comparison, note it was known that discretized kinetic Langevin dynamics can lead to √ d dependence on dimension (Cheng & Bartlett, 2018; Dalalyan & Riou-Durand, 2020) and some believe that it is the introduction of momentum that improves the dimension dependence, but our result shows that discretized overdamped Langevin (no momentum) can also have mixing time scaling like √ d. In fact, it is important to mention that recently shown was that Metropolis-Adjusted Euler-Maruyama discretization of overdamped Langevin (i.e., MALA) has an optimal dimension dependence of Õ (√ d )
(Chewi et al., 2021), while what we analyze here is the unadjusted version (i.e., LMC), and it has the same dimension dependence (note however that our dependence is not as good as that for MALA; more discussion in Section 4). We also constructed an example which shows that the mixing time of LMC is at least Ω̃ (√ d ) . Hence, our
mixing time bound has the optimal dependence on both d and , in terms of order, for the family of target measures satisfying those regularity assumptions. Our theoretical analysis is further validated by empirical investigation of numerical examples.
A minor contribution of this work is the error analysis framework that we use. It is based on the classical mean-square analysis (Milstein & Tretyakov, 2013) in numerical SDE literature, however extended from finite time to infinite time. It is a minor contribution because this extension was
1possibly after a coordinate transformation
already pioneered in the milestone work of Li et al. (2019), although we will develop a refined version. Same as in classical mean-square analysis and in Li et al. (2019), the final (sampling in this case) error is only half order lower than the order of local strong integration error (p2). This will
lead to a Õ ( C 1 p2− 1 2 / 1 p2− 1 2 ) mixing time upper bound in 2-Wasserstein distance for the family of
algorithms, where C is a constant containing various information of the underlying problem, e.g., the dimension d. Nevertheless, the following two are new to this paper: (i) We weakened the requirement on local strong and weak errors. More precisely, Li et al. (2019) requires uniform bounds on local errors, but this could be a nontrivial requirement for SDE integrators; the improvement here only requires non-uniform bounds (although establishing the same result consequently needs notably more efforts, these are included in this paper too). (ii) The detailed expressions of our bounds are not the same as those in Li et al. (2019) (even if local errors could be uniformly bounded), and as we are interested in dimension-dependence of LMC, we work out constants and carefully track their dimension-dependences. Bounds and constants in Li et al. (2019) might not be specifically designed for tightly tracking dimension dependences, as the focus of their seminal paper was more on dependence; consequently, its general error bound only led to a Õ(d)-dependence in mixing time when applied to LMC (see Example 1 in Li et al. (2019)), whereas our result leads to Õ( √ d).

heading: 2 PRELIMINARIES
text: Notation Use symbol x to denote a d-dim. vector, and plain symbol x to denote a scalar variable. ‖x‖ denotes the Euclidean vector norm. A numerical algorithm is denoted by A and its k-th iterate by x̄k. Slightly abuse notation by identifying measures with their density function w.r.t. Lebesgue measure. Use the convention Õ (·) = O(·) logO(1)(·) and Ω̃ (·) = Ω(·) logO(1)(·), i.e., the Õ (·)/Ω̃(·) notation ignores the dependence on logarithmic factors. Use the notation Ω̃(·) similarly. Denote
2-Wasserstein distance by W2(µ1, µ2) = ( inf(X,Y )∼Π(µ1,µ2) E ‖X − Y ‖ 2 ) 1 2
, where Π(µ1, µ2) is the set of couplings, i.e. all joint measures with X and Y marginals being µ1 and µ2. Denote the target distribution by µ and the law of a random variableX by Law(X). Finally, denote the mixing time of an sampling algorithm A converging to its target distribution µ in 2-Wasserstein distance by τmix( ;W2;A) = inf{k ≥ 0|W2(Law(x̄k), µ) ≤ }.
SDE for Sampling Consider a general SDE
dxt = b(t,xt)dt+ σ(t,xt)dBt (1)
where b ∈ Rd is a drift term, σ ∈ Rd×l is a diffusion coefficient matrix and Bt is a l-dimensional Wiener process. Under mild condition (Pavliotis, 2014, Theorem 3.1), there exists a unique strong solution xt to Eq. (1). Some SDEs admit geometric ergodicity, so that their solutions converge exponentially fast to a unique invariant distribution, and examples include the classical overdamped and kinetic Langevin dynamics, but are not limited to those (e.g., Mou et al. (2021); Li et al. (2020)). Such SDE are desired for sampling purposes, because one can set the target distribution to be the invariant distribution by choosing an SDE with an appropriate potential, and then solve the solution xt of the SDE and push the time t to infinity, so that (approximate) samples of the target distribution can be obtained. Except for a few known cases, however, explicit solutions of Eq. (1) are elusive and we have to resort to numerical schemes to simulate/integrate SDE. Such example schemes include, but are not limited to Euler-Maruyama method, Milstein methods and Runge-Kutta method (e.g., Kloeden & Platen (1992); Milstein & Tretyakov (2013)). With constant stepsize h and at k-th iteration, a typical numerical algorithm takes a previous iterate x̄k−1 and outputs a new iterate x̄k as an approximation of the solution xt of Eq. (1) at time t = kh.
Langevin Monte Carlo Algorithm LMC algorithm is defined by the following update rule
x̄k = x̄k−1 − h∇f(x̄k−1) + √ 2hξk, k = 1, 2, · · · (2)
where {ξk}k∈Z>0 are i.i.d. standard d-dimensional Gaussian vectors. LMC corresponds to an EulerMaruyama discretization of the continuous overdamped Langevin dynamics dxt = −∇f(xt)dt+√
2dBt, which converges to an equilibrium distribution µ ∼ exp(−f(x)).
Dalalyan (2017b) provided a non-asymptotic analysis of LMC. An Õ ( d 2 ) mixing time bound in
W2 for log-smooth and log-strongly-convex target measures (Dalalyan, 2017a; Cheng et al., 2018a; Durmus et al., 2019) has been established. It was further improved to Õ ( d ) under an additional
Hessian Lipschitz condition (Durmus & Moulines, 2019). Mixing time bounds of LMC in other statistical distances/divergences have also been studied, including total variation distance (Dalalyan, 2017b; Durmus & Moulines, 2017) and KL divergence (Cheng & Bartlett, 2018).
Classical Mean-Square Analysis A powerful framework for quantifying the global discretization error of a numerical algorithm for Eq. (1), i.e., ek = { E ‖xkh − x̄k‖ } 1 2 , is mean-square analysis (e.g., Milstein & Tretyakov (2013)). Mean-square analysis studies how local integration error propagate and accumulate into global integration error; in particular, if one-step (local) weak error and strong error (both the exact and numerical solutions start from the same initial value x) satisfy
‖Exh − Ex̄1‖ ≤C1 ( 1 + E ‖x‖2 ) 1 2
hp1 , (local weak error)( E ‖xh − x̄1‖2 ) 1 2 ≤C2 ( 1 + E ‖x‖2 ) 1 2 hp2 , (local strong error) (3)
over a time interval [0,Kh] for some constants C1, C2 > 0, p2 ≥ 12 and p1 ≥ p2 + 1 2 , then the global error is bounded by ek ≤ C ( 1 + E ‖x0‖2 ) 1 2 hp2− 1 2 , k = 1, · · · ,K for some constant C > 0
dependent on Kh.
Although classical mean-square analysis is only concerned with numerical integration error, sampling error can be also inferred. However, there is a limitation that prevents its direct employment in analyzing sampling algorithms: the global error bound only holds in finite time because the constant C can grow exponentially as K increases, rendering the bound useless when K →∞.

heading: 3 MEAN-SQUARE ANALYSIS OF SAMPLERS BASED ON CONTRACTIVE SDE
text: We now review and present some improved results on how to use mean-square analysis of integration error to quantify sampling error. A seminal paper in this direction is Li et al. (2019). What is known / new will be clarified. In all cases, the first step is to lift the finite time limitation when the SDE being discretized has some decaying property so that local integration errors do not amplify with time.
The specific type of decaying property we will work with is contractivity (after coordinate transformation). It is a sufficient condition for the underlying SDE to converge to a statistical distribution. Definition 3.1. A stochastic differential equation is contractive if there exists a non-singular constant matrix A ∈ Rd×d, a constant β > 0, such that any pair of solutions of the SDE satisfy(
E ∥∥A (xt − yt)∥∥2) 12 ≤ (E∥∥A (x− y)∥∥2) 12 exp(−βt), (4)
where xt,yt are two solutions, driven by the same Brownian motion but evolved respectively from initial conditions x and y. Remark. As long as b and σ in (1) are not explicitly dependent on time, it suffices to find an arbitrarily small t0 > 0 and show (4) holds for all t < t0. Remark. Sometimes contraction is not easy to establish directly, but can be shown after an appropriate coordinate transformation, see (Dalalyan & Riou-Durand, 2020, Proposition 1) for such a treatment for kinetic Langevin dynamics. The introduction of A permits such transformations.
In particular, overdamped Langevin dynamics, of which LMC is a discretization, is contractive when f is strongly convex and smooth.
We now use contractivity to remove the finite time limitation. We first need a short time lemma. Lemma 3.2. (Milstein & Tretyakov, 2013, Lemma 1.3) Suppose b and σ in Eq.(1) are Lipschitz continuous. For two solutions xt,yt of Eq. (1) starting from x,y respectively, denote zt(x,y) := (xt − x)− (yt − y), then there exist C0 > 0 and h0 > 0 such that
E ∥∥zt(x,y)∥∥2 ≤ C0 ‖x− y‖2 t, ∀x,y, 0 < t ≤ h0. (5)
Then we have a sequence of results that connects statistical property with integration property. We will see that a non-asymptotic sampling error analysis only requires bounding the orders of local weak and strong integration errors (if the continuous dynamics can be shown contractive). Theorem 3.3. (Global Integration Error, Infinite Time Version) Suppose Eq.(1) is contractive with rate β and with respect to a non-singular matrix A ∈ Rd×d, with Lipschitz continuous b and σ, and there is a numerical algorithm A with step size h simulating the solution xt of the SDE, whose iterates are denoted by x̄k, k = 0, 1, · · · . Suppose there exists 0 < h0 ≤ 1, C1, C2 > 0, D1, D2 ≥ 0, p1 ≥ 1, 12 < p2 ≤ p1 − 1 2 such that for any 0 < h ≤ h0, the algorithm A has, respectively, local weak and strong error of order p1 and p2, defined as ∥∥E (xh − x̄1)∥∥ ≤ (C1 +D1√E ‖x‖2)hp1 ,( E ‖xh − x̄1‖2 ) 1 2 ≤ ( C22 +D 2 2E ‖x‖ 2 ) 1 2 hp2 , (6)
where xh solves Eq.(1) with any initial value x and x̄1 is the result of applying A to x for one step.
If the solution of SDE xt and algorithm A both start from x0, then for 0 < h ≤ h1 ,
min h0, 14β ,( √β4√2κAD2) 1p2− 12 , (
β
8 √
2κ2A(D1+C0D2)
) 1 p2− 1 2 , the global error ek is bounded as ek := ( E‖xkh − x̄k‖2 ) 1 2 ≤ Chp2− 12 , k = 0, 1, 2, · · · , where (7)
C = 2√ β κ2A
( C1 + C0C2 +
√ 2U(D1 + C0D2)√ β + C2 + √ 2D2U
) , (8)
C0 is from Eq. (5), κA is the condition number of matrix A and U2 , 4E ‖x0‖2 + 6Eµ ‖x‖2. Remark (what’s new). Thm.3.3 refines the seminal results in Li et al. (2019) in the sense that it only requires non-uniform bounds on the local error (6), whereas Li et al. (2019) requires uniform bounds, i.e., D1 = D2 = 0 in (6). Therefore, the refinement we present has wider applicability.
In general, local errors tend to depend on the current step’s value, i.e. D1 6= 0, D2 6= 0. Allowing local error bounds to be non-uniform enabled applications such as proving the vanishing bias of mirror Langevin algorithm (Li et al., 2021). For a simpler illustration, consider LMC for 1D standard Gaussian target distribution, then we have∥∥E (xh − x̄1)∥∥ = (e−h − 1 + h)E‖x‖ = (h2
2 + o(h2))E‖x‖.
One can see that the local error does depend on x and is not uniform. Meanwhile, our non-uniform condition still holds because ( h2 2 + o(h 2) ) E‖x‖ ≤ (h 2 2 + o(h 2)) √
E‖x‖2 (and thus p1 = 2). Note if the discretization does converge to a neighborhood of the target distribution, it is possible that E‖x‖2 and/or E‖x‖ become bounded near the convergence, and in this case the D1, D2 parts can be absorbed into C1 and C2; however, this ‘if’ clause is exactly what we’d like to prove.
Nevertheless, we state for rigor that the convention 1/0 =∞ is used when D1 = D2 = 0. Another remark is, even in this case, our bound has a different expression from the seminal results. We will carefully work out, track, and combine dimension-dependence of constants using our bound.
Following Theorem 3.3, we obtain the following non-asymptotic bound of the sampling error in W2:
Theorem 3.4. (Non-Asymptotic Sampling Error Bound: General Case) Under the same assumption and with the same notation of Theorem 3.3, we have
W2(Law(x̄k), µ) ≤ e−βkhW2(Law(x0), µ) + Chp2− 1 2 , ∀0 < h ≤ h1.
A corollary of Theorem 3.4 is a bound on the mixing time of the sampling algorithm:
Corollary 3.5. (Upper Bound of Mixing Time: General Case) Under the same assumption and with the same notation of Theorem 3.3, we have
τmix( ;W2;A) ≤ max  1βh1 , 1β ( 2C ) 1 p2− 1 2  log 2W2(Law(x0), µ) In particular, when high accuracy is needed, i.e., < 2Chp2− 1 2
1 , we have
τmix( ;W2;A) ≤ (2C)
1
p2− 1 2
β
1
1
p2− 1 2
log 2W2(Law(x0), µ) = Õ C 1p2− 12 β
1
1
p2− 1 2  . (9) Corollary 3.5 states how mixing time depends on the order of local (strong) error (i.e., p2) of a numerical algorithm. The larger p2 is, the shorter the mixing time of the algorithm is, in term of the dependence on accuracy tolerance parameter . It is important to note that for constant stepsize discretizations that are deterministic on the filtration of the driving Brownian motion and use only its increments, there is a strong order barrier, namely p2 ≤ 1.5 (Clark & Cameron, 1980; Rüemelin, 1982); however, methods involving multiple stochastic integrals (e.g., Kloeden & Platen (1992); Milstein & Tretyakov (2013); Rößler (2010)) can yield a larger p2, and randomization (e.g., Shen & Lee (2019)) can possibly break the barrier too.
The constant C defined in Eq. (7) typically contains rich information about the underlying SDE, e.g. dimension, Lipschitz constant of drift and noise diffusion, and the initial value x0 of the sampling algorithm. Through C, we can uncover the dependence of mixing time bound on various parameters, such as the dimension d. This will be detailed for Langevin Monte Carlo in the next section.
It is worth clarifying that once Thm.3.3 is proved, establishing Theorem 3.4 and Corollary 3.5 is relatively easy. In fact, analogous results have already been provided in Li et al. (2019), although they also required uniform local errors as consequences of their Thm.1. Nevertheless, we do not claim novelty in Theorem 3.4 and Corollary 3.5 and they are just presented for completeness. Our main refinement is just Thm.3.3 over Thm.1 in Li et al. (2019), and the non-triviality lies in its proof.

heading: 4 NON-ASYMPTOTIC ANALYSIS OF LANGEVIN MONTE CARLO ALGORITHM
text: We now quantify how LMC samples from Gibbs target distribution µ ∼ exp ( −f(x) ) that has a
finite second moment, i.e., ∫ Rd ‖x‖ 2 dµ <∞. Assume without loss of generality that the origin is a local minimizer of f , i.e. ∇f(0) = 0; this is for notational convenience in the analysis and can be realized via a simple coordinate shift, and it is not needed in the practical implementation. In addition, we assume the following two conditions hold: A 1. (Smoothness and Strong Convexity) Assume f ∈ C2 and is L-smooth and m-strongly-convex, i.e. there exists 0 < m ≤ L such that mId 4 ∇2f(x) 4 LId, ∀x ∈ Rd.
Denote the condition number of f by κ , Lm . The smoothness and strong-convexity assumption is the standard assumption in the literature of analyzing LMC algorithm (Dalalyan, 2017a;b; Cheng & Bartlett, 2018; Durmus et al., 2019; Durmus & Moulines, 2019). A 2. (Linear Growth of the 3rd-order Derivative) Assume f ∈ C3 and the operator ∇(∆f) grows at most linearly, i.e., there exists a constant G > 0 such that
∥∥∇(∆f(x))∥∥ ≤ G (1 + ‖x‖). Remark. The linear growth (at infinity) condition on∇∆f is actually not as restrictive as it appears, and in some sense even weaker than some classical condition for the existence of solutions to SDE. For example, a standard condition for ensuring the existence and uniqueness of a global solution to SDE is at most a linear growth (at infinity) of the drift (Pavliotis, 2014, Theorem 3.1). If we consider monomial potentials, i.e., f(x) = xp, p ∈ N+, then the linear growth condition on∇∆f is met when p ≤ 4, whereas the classical condition for the existence of solutions holds only when p ≤ 2. Remark. Another additional assumption, namely Hessian Lipschitz condition, is commonly used in the literature (e.g., Durmus & Moulines (2019); Ma et al. (2021)). It requires the existence of a constant L̃, such that ‖∇2f(y)−∇2f(x)‖ ≤ L̃‖y − x‖. It can be shown that smoothness and Hessian Lipschitzness imply A2. Meanwhile, examples that satisfy A2 but are not Hessian Lipschitz exist, e.g., f(x) = x4, and thus A2 is not necessarily stronger than Hessian Lipschitzness.
Remark. Same as L and m in A1, we implicitly assume the constant G introduced in A2 to be independent of dimension. Meanwhile, it is important to note examples for which G depends on the dimension do exist, and this is also true for other regularity constants including not only L and m but also the Hessian Lipschitz constant L̃. This part of the assumption is a strong one. Remark. A2, together with Ito’s lemma, helps establish an order p1 = 2 of local weak error for LMC (see Lemma D.2), which enables us to obtain the √ d dependence.
To apply mean-square analysis to study LMC algorithm, we will need to ensure the underlying Langevin dynamics is contractive, which we verify in Section C and D in the appendix. In addition, we work out all required constants to determine the C in Eq. 7 explicitly in the appendix. With all these necessary ingredients, we now invoke Theorem 3.4 and obtain the following result: Theorem 4.1. (Non-Asymptotic Error Bound: LMC) Suppose Assumption 1 and 2 hold. LMC iteration x̄k+1 = x̄k − h∇f(x̄k) + √ 2hξk satisfies
W2(Law(x̄k), µ) ≤ e−mkhW2(Law(x0), µ) + CLMCh, 0 < h ≤ 1
4κL , k ∈ N (10)
where CLMC = 10(L2+G)
m 3 2
√ 2d+m ( E ‖x0‖2 + 1 ) = O( √ d).
Corollary 3.5 combined with the above result gives the following bound on the mixing time of LMC: Theorem 4.2. (Upper Bound of Mixing Time: LMC) Suppose Assumption 1 and 2 hold. If running LMC from x0, we then have
τmix( ;W2; LMC) ≤ max {
4κ2, 2CLMC m 1
} log 2W2(Law(x0), µ)
where CLMC is the same in Theorem 4.1. When high accuracy is needed, i.e., ≤ CLMC2mκ2 , we have
τmix( ;W2; LMC) ≤ 2CLMC m 1 log 2W2(Law(x0), µ) = Õ
(√ d ) .
The Õ (√ d ) mixing time bound in W2 distance improves upon the previous ones (Dalalyan, 2017a;
Cheng & Bartlett, 2018; Durmus & Moulines, 2019; Durmus et al., 2019) in the dependence of d and/or . If further assuming G = O(L2), we then have CLMC = O(κ2 √ m √ d) and Thm.4.2 shows the mixing time is Õ ( κ2√ m √ d ) , which also improves the κ dependence in some previous results (Dalalyan, 2017a; Cheng & Bartlett, 2018) in the m ≤ 1 regime. A brief summary is in Table 1. Remark (more comparison). The seminal work of Li et al. (2019) provided mean-square analysis (their Thm.1) and obtained a Õ ( d ) mixing time bound for LMC (their Example 1) under smoothness,
strong convexity and Hessian Lipschitz conditions, consistent with that in Durmus & Moulines (2019). By using our version (Thm.3.3) and tracking down constants’ dimension-dependence, we are able to tighten it to Õ (√ d ) . Worth clarifying is, dimension-dependence might not be the focus of Li et al.
(2019); instead, it considered -dependence and other discretizations, and showed for example that 1.5 SRK discretization has improved mixing time bound of Õ ( d 2/3 ) . The dimension dependence of
this discretization, for example, can possibly be improved by our results too. Optimality In fact, the Õ (√ d ) mixing time of LMC has the optimal scaling one can expect. This
is in terms of dependence on d and , over the class of all log-smooth and log-strongly-convex target measures. To illustrate this, consider the following Gaussian target distribution whose potential is
f(x) = m
2 d∑ i=1 x2i + L 2 2d∑ i=d+1 x2i , with m = 1, L ≥ 4m. (11)
We now establish a lower bound on the mixing time of LMC algorithm for this target measure. 1The dependence on κ is not readily available from Theorem 8 in Durmus & Moulines (2019). 2The G = O(L2) assumption is only for κ,m dependence. Removing it doesn’t affect d, dependence.
Theorem 4.3. (Lower Bound of Mixing Time) Suppose we run LMC for the target measure defined in Eq. (11) from x0 = 12d, then for any choice of step size h > 0 within stability limit, we have
τmix( ;W2; LMC) ≥ √ d
8 log
√ d
= Ω̃
(√ d ) .
Combining Theorem 4.2 and 4.3, we see that mean-square analysis provides a tight bound for LMC and Õ (√ d ) is the optimal scaling of LMC for target measures satisfying Assumptions 1 and 2.
Note that the above optimality results only partially, but not completely, close the gap between the upper and lower bounds of LMC over the entire family of log-smooth and log-strongly-convex target measures, because of one limitation of our result — A(ssumption)2 is, despite of its close relation to the Hessian Lipschitz condition frequently used in the literature, still an extra condition. We tend to believe that A2 may not be essential, but rather an artifact of our proof technique. However, at this moment we cannot eliminate the possibility that the best scaling one can get out of Assumption 1 only (no A2) is worse than Õ (√ d/ ) . We’d like to further investigate this in future work.
Discussion Besides Li et al. (2019) (see the previous Remark), let’s briefly discuss three more important sampling algorithms related to LMC. Two of them are Kinetic Langevin Monte Carlo (KLMC) and Randomized Midpoint Algorithm (RMA), both of which are discretizations of kinetic Langevin dynamics. The other is Metropolis-Adjusted Langevin Algorithm (MALA) which uses the one-step update of LMC as a proposal and then accepts/rejects it with Metropolis-Hastings.
The Õ (√ d ) mixing time in 2-Wasserstein distance of KLMC has been established for log-smooth
and log-strongly-convex target measures in existing literature (Cheng et al., 2018b; Dalalyan & RiouDurand, 2020) and that was a milestone. Due to its better dimension dependence over previously best known results of LMC, KLMC is understood to be the analog of Nesterov’s accelerated gradient method for sampling (Ma et al., 2021). Our findings show that LMC is able to achieve the same mixing time, although under an additional growth-at-infinity condition. However, this does not say anything about whether/how KLMC accelerates LMC, as the existing KLMC bound may still be not optimal. We also note KLMC has better condition number dependence than our current LMC result, although the κ dependence in our bound may not be tight.
RMA (Shen & Lee, 2019) is based on a brilliant randomized discretization of kinetic Langevin dynamics and shown to have further improved dimension dependence (and other pleasant properties). From the perspective of this work, we think it is because RMA is able to break the strong order barrier due to the randomization, and more investigations based on mean-square analysis should be possible.
For MALA, a recent breakthrough (Chewi et al., 2021) establishes a Õ (√ d )
mixing time in W2 distance with warm start, and the dimension dependence is shown to be optimal. We see that without the Metropolis adjustment, LMC (under additional assumptions such as A2) can also achieve the same dimension dependence as MALA. But unlike LMC, MALA only has logarithmic dependence on 1 . With warm-start, is it possible/how to improve the dependence of 1 for LMC, from polynomial to logarithmic? This question is beyond the scope of this paper but worth further investigation.

heading: 5 NUMERICAL EXAMPLES
text: This section numerically verifies our theoretical findings for LMC in Section 4, with a particular focus on the dependence of the discretization error in Theorem 4.1 on dimension d and step size h. To this end, we consider two target measures specified by the following two potentials:
f1(x) = 1
2 ‖x‖2 + log  d∑ i=1 exi  and f2(x) = 1 2 ‖x‖2 − 1 2d 1 2 d∑ i=1 cos ( d 1 4xi ) . (12)
It is not hard to see f1 is 2-smooth and 1-strongly convex, f2 is 32 -smooth and 1-strongly-convex, and both satisfy Assumption 2. f2 is also used in (Chewi et al., 2021) to illustrate the optimal dimension dependence of MALA. Explicit expression of 2-Wasserstein distance between non-Gaussian distributions is typically not available, instead, we use the Euclidean norm of the mean error as a surrogate because
∥∥Ex̄k − Eµx∥∥ ≤ W2(Law(x̄k), µ) due to Jensen’s inequality. To obtain an accurate estimate of the ground truth, we run 108 independent LMC realizations using a tiny step size (h = 0.001), each till a fixed, long enough time, and use the empirical average to approximate Eµx.
To study the dimension dependence of sampling error, we fix step size h = 0.1, and for each d ∈ {1, 2, 5, 10, 20, 50, 100, 200, 500, 1000}, we simulate 104 independent Markov chains using LMC algorithm for 100 iterations, which is long enough for the chain to be well-mixed. The mean and the standard deviation of the sampling error corresponding to the last 10 iterates are recorded.
To study step size dependence of sampling error, we fix d = 10 and experiment with step size h ∈ {1, 2, 3, 4, 5, 6, 7, 8, 9, 10}×10−1. We run LMC till T = 20, i.e., dTh e iterations for each h. The procedure is repeated 104 times with different random seeds to obtain independent samples. When the corresponding continuous time t = kh > 10, we see from Eq. (10) that LMC is well converged and the sampling error is saturated by the discretization error. Therefore, for each h, we take the last d 10h e iterates and record the mean and standard deviation of their sampling error.
Results shown in Fig.1 are consistent with our theoretical analysis of the sampling error. Both linear dependence on √ d and h can be supported by the empirical evidence. Note results with smaller h are less accurate because one starts to see the error of empirical approximation due to finite samples.

heading: 6 CONCLUSION
text: Via a refined mean-square analysis of Langevin Monte Carlo algorithm, we obtain an improved and optimal Õ (√ d/ )
bound on its mixing time, which was previously thought to be obtainable only with the addition of momentum. This was under the standard smoothness and strongly-convexity assumption, plus an addition linear growth condition on the third-order derivative of the potential function, similar to Hessian Lipschitz condition already popularized in the literature.
Here are some possible directions worth further investigations. (i) Combine mean-square analysis with stochastic gradient analysis to study SDE-based stochastic gradient MCMC methods; (ii) Is it still possible to obtain √ d-dependence without A2, i.e., only under log-smooth and log-stronglyconvex conditions? (iii) Applications of mean-square analysis to other SDEs and/or discretizations; (iv) Motivated by Chewi et al. (2021), it would be interesting to know whether the dependence on 1 can be improved to logarithmic, for example if LMC is initialized at a warm start.

heading: ACKNOWLEDGMENTS
text: The authors thank Andre Wibisono, Chenxu Pang, anonymous reviewers and area chair for suggestions that significantly improved the quality of this paper. MT was partially supported by NSF DMS1847802 and ECCS-1936776. This work was initiated when HZ was a professor and RL was a PhD student at Georgia Institute of Technology.

heading: A PROOF OF RESULTS IN SECTION 3

heading: A.1 PROOF OF THEOREM 3.3 (GLOBAL INTEGRATION ERROR, INFINITE TIME VERSION)
text: Proof. We write the solution of an SDE by xt0,xt0 (t0 + t) when the dependence on initialization needs highlight. Denote tk = kh and xtk = xk for better readability.
We will first make an easy observation that contraction and bounded 2nd-moment of the invariant distribution lead to bounded 2nd-moment of the SDE solution for all time: let y0 be a random variable following the invariant distribution of Eq. (1), i.e., y0 ∼ µ, then yt ∼ µ and
E ‖xt‖2 ≤2E ‖xt − yt‖ 2 + 2E ‖yt‖ 2
≤2E ‖x0 − y0‖ 2 exp(−2βt) + 2E ‖yt‖ 2
≤4E(‖x0‖2 + ‖y0‖ 2 ) exp(−2βt) + 2E ‖yt‖ 2 =4E ‖x0‖2 exp(−2βt) + ( 2 + 4 exp(−2βt) ) Ey∼µ ‖y‖2
≤4E ‖x0‖2 + 6 ∫ Rd ‖y‖2 dµ , U2
and then it follows that
E ‖x̄k‖2 ≤ 2E ‖x̄k − xk‖2 + 2E ‖xk‖2 ≤ 2e2k + 2U2. (13)
Denote 〈x,y〉A = 〈Ax, Ay〉, ‖x‖A = ‖Ax‖ and
fk = { E ‖xk − x̄k‖2A } 1 2
(14)
where A is the non-singular matrix from Equation (4). Also denote that largest and smallest singular values of A by σmax and σmin, respectively, and the condition number of A by κA = σmaxσmin . Recall ek = E ‖xk − x̄k‖, it is easy to see that
σminek ≤ fk ≤ σmaxek. (15)
Further, we have the following decomposition
f2k+1 =E ‖xk+1 − x̄k+1‖ 2 A =E ∥∥∥xtk,xtk (tk+1)− xtk,x̄k(tk+1) + xtk,x̄k(tk+1)− x̄k+1∥∥∥2A
=E ∥∥∥xtk,xtk (tk+1)− xtk,x̄k(tk+1)∥∥∥2A︸ ︷︷ ︸
1
+E ∥∥xtk,x̄k(tk+1)− x̄k+1∥∥2A︸ ︷︷ ︸
2
(16)
+ 2E〈A ( xtk,xtk (tk+1)− xtk,x̄k(tk+1) ) , A ( xtk,x̄k(tk+1)− x̄k+1 ) 〉︸ ︷︷ ︸
3
.
Term 1 is taken care of the contraction property
E ∥∥∥xtk,xtk (tk+1)− xtk,x̄k(tk+1)∥∥∥2A ≤ f2k exp(−2βh). (17)
Term 2 is dealt with by the bound on local strong error
E ∥∥xtk,x̄k(tk+1)− x̄k+1∥∥2A ≤ σ2max (C22 +D22E ‖x̄k‖2)h2p2 . (18)
Term 3 requires more efforts to cope with, and by the decomposition in Eq. (5) we have
E〈(xtk,xtk (tk+1)− xtk,x̄k(tk+1),xtk,x̄k(tk+1)− x̄k+1〉A =E〈xk − x̄k,xtk,x̄k(tk+1)− x̄k+1〉A + E〈zh(xk, x̄k),xtk,x̄k(tk+1)− x̄k+1〉A (i) =E〈xk − x̄k,E[xtk,x̄k(tk+1)− x̄k+1|Fk]〉A + E〈zh(xk, x̄k),xtk,x̄k(tk+1)− x̄k+1〉A
(ii) ≤ fk ( E ∥∥E[xtk,x̄k(tk+1)− x̄k+1|Fk]∥∥2A) 12 + (E∥∥zh(xk, x̄k)∥∥2A) 12 (E∥∥xtk,x̄k(tk+1)− x̄k+1∥∥2A) 12
(iii) ≤ σmaxfk ( E ∥∥E[xtk,x̄k(tk+1)− x̄k+1|Fk]∥∥2) 12 + σ2max (E∥∥zh(xk, x̄k)∥∥2) 12 (E∥∥xtk,x̄k(tk+1)− x̄k+1∥∥2) 12
(iv) ≤ σmaxfk ( C1 +D1 √ E ‖x̄k‖2 ) hp1 + κAσmaxC0fk √ h ( C2 +D2 √ E ‖x̄k‖2 ) hp2
(v) ≤κAσmax(C1 + C0C2)ekhp2+ 1 2 + κAσmax(D1 + C0D2) √ E ‖x̄k‖2fkhp2+ 1 2 (19)
where (i) uses the tower property of conditional expectation and Fk is the filtration at k-th iteration, (ii) uses Cauchy-Schwarz inequality, (iii) is due to the relationship between ek and fk, (iv) is due to local weak error, local strong error and Eq. (5), and (v) is due to p1 ≥ p2 + 12 and 0 < h ≤ h0 ≤ 1. Now plug Eq. (17), (18) and (19) in Eq. (16), we obtain f2k+1 ≤f2k exp(−2βh) + σ2max ( C22 +D 2 2E ‖x̄k‖ 2 ) h2p2 + κAσmax(C1 + C0C2)fkh p2+ 1 2
+ κAσmax(D1 + C0D2) √ E ‖x̄k‖2fkhp2+ 1 2
(i) ≤ (
1− 7 8 βh
) f2k + σ 2 max ( C22 +D 2 2E ‖x̄k‖ 2 ) h2p2 + κAσmax(C1 + C0C2)fkh p2+ 1 2
+ κAσmax(D1 + C0D2) √ E ‖x̄k‖2fkhp2+ 1 2
(ii) ≤ (
1− 7 8 βh
) f2k + κAσmax ( C1 + C0C2 + √ 2U(D1 + C0D2) ) fkh p2+ 1 2 + 2κ2AD 2 2f 2 kh 2p2
+ √
2κ2A(D1 + C0D2)f 2 kh p2+ 1 2 + σ2max ( C22 + 2D 2 2U 2 ) h2p2
(iii) ≤ (
1− 7 8 βh
) f2k + κAσmax ( C1 + C0C2 + √ 2U(D1 + C0D2) ) fkh p2+ 1 2 + 3β
8 f2kh
+ σ2max ( C22 + 2D 2 2U 2 ) h2p2
= ( 1− 1
2 βh
) f2k + κAσmax ( C1 + C0C2 + √ 2U(D1 + C0D2) ) fkh p2+ 1 2
+ σ2max ( C22 + 2D 2 2U 2 ) h2p2
(iv) ≤ (
1− 1 2 βh
) f2k + β
4 f2kh+
κ2Aσ 2 max ( C1 + C0C2 + √ 2U(D1 + C0D2) )2 β h2p2
+ σ2max ( C22 + 2D 2 2U 2 ) h2p2
= ( 1− 1
4 βh
) f2k + κ 2 Aσ 2 max  ( C1 + C0C2 + √ 2U(D1 + C0D2) )2 β + C22 + 2D 2 2U 2 h2p2
where (i) is due to the assumption 0 < h ≤ 14β and e −x ≤ 1 − x + x
2
2 for 0 < x < 1,
(ii) is due to the upper bound on E ‖x̄k‖2 in Eq. (13), (iii) holds provided when h ≤
min ( √β4√2κAD2) 1p2− 12 , (
β
8 √
2κ2A(D1+C0D2)
) 1 p2− 1 2  and (iv) is due to Cauchy-Schwarz inequality.
Unfolding the above inequality gives us
f2k+1 ≤ 4
β κ2Aσ 2 max
 ( C1 + C0C2 + √ 2U(D1 + C0D2) )2 β + C22 + 2D 2 2U 2 h2p2−1. Taking square root on both sides and using √ a2 + b2 + c2 ≤ a+ b+ c,∀a, b, c ≥ 0 yields
fk+1 ≤ 2√ β κAσmax
( C1 + C0C2 +
√ 2U(D1 + C0D2)√ β + C2 + √ 2D2U
) hp2− 1 2 .
Finally using the relationship between ek and fk, we obtain
ek ≤ 2√ β κ2A
( C1 + C0C2 +
√ 2U(D1 + C0D2)√ β + C2 + √ 2D2U
) hp2− 1 2 .

heading: A.2 PROOF OF THEOREM 3.4 (NON-ASYMPTOTIC SAMPLING ERROR BOUND: GENERAL CASE)
text: Proof. Let y0 ∼ µ and (x0,y0) are coupled such that E ‖x0 − y0‖ 2 = W 22 (Law(x0), µ). Denote the solution of Eq. (1) starting from x0,y0 by xt,yt respectively, and tk = kh. We have
W2(Law(x̄k), µ) ≤W2(Law(x̄k),Law(xtk)) +W2(Law(xtk), µ) ≤ √ E ∥∥x̄k − xtk∥∥2 +√E∥∥∥xtk − ytk∥∥∥2
(i) ≤ek + √ E ‖x0 − y0‖ 2 exp (−2βtk)
=ek + exp (−βtk)W2(Law(x0), µ)
where (i) is due to the contraction assumption on Eq. (1). Invoking the conclusion of Theorem 3.3 completes the proof.

heading: A.3 PROOF OF COROLLARY 3.5 (UPPER BOUND OF MIXING TIME: GENERAL CASE)
text: Proof. Given any tolerance > 0, we know from Theorem 3.4 that if k is large enough and h is small enough such that
exp (−βkh)W2(Law(x0), µ) ≤ 2 . (20)
Chp2− 1 2 ≤
2 (21)
we then have W2(Law(x̄k), µ) ≤ . Solving Inequality (20) yields
k ≥ 1 βh log 2W2(Law(x0), µ) , k? (22)
To minimize the lower bound, we want pick step size h as large as possible. Besides h ≤ h1, Eq. (21) poses further constraint on h, hence we have
h ≤ min h1, ( 2C ) 1 p2− 1 2  .
Plug the upper bound of h in Eq. (22), we have
k? = max  1βh1 , 1β ( 2C ) 1 p2− 1 2  log 2W2(Law(x0), µ) . When high accuracy is needed, i.e., < 2Chp2− 1 2
1 , we have
k? = (2C)
1
p2− 1 2
β
1
1
p2− 1 2
log 2W2(Law(x0), µ) = Õ C 1p2− 12 β
1
1
p2− 1 2
 .

heading: B PROOF OF RESULTS IN SECTION 4

heading: B.1 PROOF OF THEOREM 4.1 (NON-ASYMPTOTIC ERROR BOUND: LMC)
text: Proof. From Lemma C.1 we know that Langevin dynamics is a member of the family of contractive SDE, and with a contraction rate of strong-convexity coefficient β = m (w.r.t. identity matrix Id×d).
Next, we will need to work out the constants C0, C1, D1, D2, C2 needed in Theorem 3.3. We have C0 = √ m 2 , implied from Lemma C.3.
The local strong error and local weak error are bounded in Lemma D.1 and D.2 respectively. Note that the coefficient C̃1/C̃2 in the bound for local strong/weak error depends on initial value, which changes from iteration to iteration. Combined with Lemma D.3, we would obtain C1 and C2, namely
C̃1 ≤ 2(L2 +G) ( d
4κL + E ‖x0‖2 +
8d 7m + 1
) 1 2 ≤ 2(L2 +G) √ 2d
m + E ‖x0‖2 + 1 , C1
and
C̃2 ≤ 2L ( d+ m
2
( E ‖x0‖2 + 8d
7m
)) 12 ≤ 2L √ m √ 2d
m + E ‖x0‖2 + 1 , C2.
We collect all constants here in the proof for easier reference
A = Id×d, κA = 1, β = m, h0 = 1
4κL , C0 =
√ m
2 ,
C1 = 2(L 2 +G)
√ 2d
m + E ‖x0‖2 + 1, D1 = 0
C2 = 2L √ m
√ 2d
m + E ‖x0‖2 + 1, D2 = 0.
Then the constant in Theorem 3.3 for LMC algorithm simplifies to
C = 2√ β
( C1 + C0C2√
β + C2
) ,
≤ 10(L 2 +G)
m 3 2
√ 2d+m ( E ‖x0‖2 + 1 ) , CLMC.
Assuming L,m,G are all constants and independent of d, then clearly CLMC = O( √ d). Then applying Theorem 3.4 to LMC, we have
W2(Law(x̄k), µ) ≤ e−mkhW2(Law(x0), µ) + CLMCh (23)
for 0 < h ≤ 14κL .

heading: B.2 PROOF OF THEOREM 4.3 (LOWER BOUND OF MIXING TIME)
text: Proof. If we start from x0 = 12d and run LMC for the potential function in Eq. (11), we then have
(x̄k)i =
{ (1−mh)k(x0)i + √ 2h ∑k l=1(1−mh)k−l(ξl)i, 1 ≤ i ≤ d
(1− Lh)k(x0)i + √ 2h ∑k l=1(1− Lh)k−l(ξl)i, d+ 1 ≤ i ≤ 2d
and hence
(x̄k)i ∼ N ( (1−mh)k, 2m(2−mh) ( 1− (1−mh)2k )) , 1 ≤ i ≤ d N ( (1− Lh)k, 2L(2−Lh) ( 1− (1− Lh)2k )) , d+ 1 ≤ i ≤ 2d
Clearly, stability requires h < 2L .
The squared 2-Wasserstein distance between the law of the k-th iterate of LMC and target distribution is
W 22 (Law(x̄k), µ) =d(1−mh)2k + d
m
(√ 2
2−mh
√ 1− (1−mh)2k − 1
)2
+d(1− Lh)2k + d L
(√ 2
2− Lh
√ 1− (1− Lh)2k − 1 )2 .
Suppose W2(Law(x̄k), µ) ≤ , we then must have
d(1−mh)2k ≤ 2 (24)
d
m
(√ 2
2−mh
√ 1− (1−mh)2k − 1 )2 ≤ 2. (25)
A necessary condition of Eq. (25) is that
1 + √ m√ d ≥
√ 2
2−mh
√ 1− (1−mh)2k (i) ≥ √ 2
2−mh
√ 1− 2
d (26)
where (i) is due to Eq. (24). It follows from Eq. (26) and m = 1 that
h ≤ 4 1 + √
d
√ d ≤ 4 √ d . (27)
Revisiting Eq. (24) yields
2 ≥ d(1−mh)2k (i) ≥ d ( 1− 2mh+ (2mh) 2
2
)2k (ii)
≥ de−4mkh
⇐⇒ k ≥ 1 2hm log
√ d
(28)
where (i) is due to mh < 2κ < 1 2 and (ii) is due to e
−x ≤ 1− x+ x 2
2 , 0 < x < 1.
Combine Eq. (27) and (28), we then obtain a lower bound of the mixing time
k ≥ √ d
8m log
√ d
=
√ d
8 log
√ d
= Ω̃
(√ d ) .

heading: C SOME PROPERTIES OF LANGEVIN DYNAMICS

heading: C.1 CONTRACTION OF LANGEVIN DYNAMICS
text: Lemma C.1. Suppose Assumption 1 holds. Then two copies of overdamped Langevin dynamics have the following contraction property
{ E ‖yt − xt‖ 2 } 1 2 ≤ { E ‖y − x‖2 } 1 2 exp(−mt)
where x,y are the initial values of xt,yt.
Proof. First assume x,y are deterministic. Suppose xt,yt are respectively the solutions to
dxt =−∇f(xt)dt+ √ 2dBt dyt =−∇f(yt)dt+ √ 2dBt
whereBt is a standard d-dimensional Brownian motion. Denote Lt = 12E ‖yt − xt‖ 2 and take time derivative, we obtain
d dt Lt = −E〈yt − xt,∇f(yt)−∇f(xt)〉
(i) ≤ −mE ‖yt − xt‖ 2 = −2mLt
where (i) is due to the strong-convexity assumption made on f . We then obtain Lt ≤ L0 exp(−2mt) and it follows by Gronwall’s inequality that
{ E ‖yt − xt‖ 2 } 1 2 ≤ ‖y − x‖ exp(−mt).
When x,y are random, by the conditioning version of the above inequality and Jensen’s inequality, we have
E [ E ‖yt − xt‖ 2 ∣∣∣∣x,y ] 1 2 ≤ { E ‖y − x‖2 exp(−2mt) } 1 2 = { E ‖y − x‖2 } 1 2 exp(−mt).

heading: C.2 GROWTH BOUND OF LANGEVIN DYNAMICS
text: Lemma C.2. Suppose Assumption 1 holds, then when 0 ≤ h ≤ 14κL , the solution of overdamped Langevin dynamics xt satisfies
E ‖xh − x‖2 ≤ 6 ( d+ m
2 E ‖x‖2
) h
where x is the initial value at t = 0.
Proof. We have
E ‖xh − x‖2 =E ∥∥∥∥∥− ∫ h
0
∇f(xt)dt+ √ 2 ∫ h 0 dBt ∥∥∥∥∥ 2
≤2E ∥∥∥∥∥ ∫ h
0
∇f(xt)dt ∥∥∥∥∥ 2 + 4E ∥∥∥∥∥ ∫ h
0
dBt ∥∥∥∥∥ 2
(i) =2E ∥∥∥∥∥ ∫ h
0
∇f(xt)dt ∥∥∥∥∥ 2 + 4hd
≤2E (∫ h 0 ∥∥∇f(xt)−∇f(x)∥∥ dt+ ∫ h 0 ∥∥∇f(x)∥∥ dt)2 + 4hd
≤2E (L∫ h 0 ‖xt − x‖ dt+ h ∥∥∇f(x)∥∥)2 + 4hd ≤4E L2(∫ h 0 ‖xt − x‖ dt )2 + h2 ∥∥∇f(x)∥∥2 + 4hd
(ii) ≤ 4hd+ 4h2E ∥∥∇f(x)∥∥2 + 4L2h∫ h
0
E ‖xt − x‖2 dt
where (i) is due to Ito’s isometry, (ii) is due to Cauchy-Schwarz inequality. By Gronwall’s inequality, we obtain
E ‖xh − x‖2 ≤ 4h ( d+ hE ∥∥∇f(x)∥∥2) exp{4L2h2} . Since
∥∥∇f(x)∥∥ = ∥∥∇f(x)−∇f(0)∥∥ ≤ L ‖x‖, when 0 < h ≤ 14κL , we finally reach at E ‖xh − x‖2 ≤ 4e 1 4 ( d+ 2hL2E ‖x‖2 ) h ≤ 6 ( d+ m
2 E ‖x‖2
) h.

heading: C.3 BOUND ON EVOLVED DEVIATION
text: Lemma C.3. Suppose Assumption 1 holds. Let xt,yt be two solutions of overdamped Langevin dynamics starting from x,y respectively, for 0 < h ≤ 14κL , we have the following representation
xh − yh = x− y + z
with
E ‖z‖2 ≤ m 4 E ‖x− y‖2 h.
Proof. Let z = (xh − yh)− (x− y) = − ∫ h 0 ∇f(xs)−∇f(ys)ds. Ito’s lemma readily implies that
E ‖xh − yh‖ 2
=E ‖x− y‖2 − 2E ∫ h
0
〈xs − ys,∇f(xs)−∇f(ys)〉ds
(i) ≤E ‖x− y‖2 − 2m ∫ h
0
E ‖xs − ys‖ 2 ds
≤E ‖x− y‖2
where (i) is due to strong-convexity of f . We then have that
E ‖z‖2 = ∥∥∥∥∥∥E [∫ h 0 ∇f(xs)−∇f(ys)ds ]∥∥∥∥∥∥ 2
≤ (∫ h 0 ∥∥∥E [∇f(xs)−∇f(ys)]∥∥∥ ds )2
≤ ∫ h
0
12ds ∫ h 0 ∥∥∥E [∇f(xs)−∇f(ys)]∥∥∥2 ds ≤h ∫ h
0
E ∥∥∇f(xs)−∇f(ys)∥∥2 ds
≤L2h ∫ h
0
E ‖xs − ys‖ 2 ds
≤L2E ‖x− y‖2 h2
(i) ≤m 4 E ‖x− y‖2 h
where (i) is due to h ≤ 14κL .

heading: D SOME PROPERTIES OF LMC ALGORITHM

heading: D.1 LOCAL STRONG ERROR
text: Lemma D.1. Suppose Assumption 1 holds. Denote the one-step iteration of LMC algorithm with step size h by x̄1 and the solution of overdamped Langevin dynamics at time t = h by xh. Both the discrete algorithm and the continuous dynamics start from the same initial value x. If 0 ≤ h ≤ 14κL , then the local strong error of LMC algorithm satisfies{
E ‖x̄1 − xh‖2 } 1 2 ≤ C̃2h 3 2
with C̃2 = 2L ( d+ m2 E ‖x‖ 2 ) 1 2 .
Proof. We have for 0 ≤ h ≤ 14κL ,
E ‖x̄1 − xh‖2 =E ∥∥∥∥∥ ∫ h
0
∇f(xs)−∇f(x)ds ∥∥∥∥∥ 2
≤E (∫ h 0 ∥∥∇f(xs)−∇f(x)∥∥ ds)2
≤L2E (∫ h 0 ‖xs − x‖ ds )2 (i)
≤L2h ∫ h
0
E ‖xs − x‖2 ds
(ii) ≤ 3L2 ( d+ m
2 E ‖x‖2
) h3
where (i) is due to Cauchy-Schwartz inequality and (ii) is due to Lemma C.2. Taking square roots on both side completes the proof.

heading: D.2 LOCAL WEAK ERROR
text: Lemma D.2. Suppose Assumption 1 and 2 hold. Denote the one-step iteration of LMC algorithm with step size h by x̄1 and the solution of overdamped Langevin dynamics at time t = h by xh. Both the discrete algorithm and the continuous dynamics start from the same initial value x. If 0 ≤ h ≤ 14κL , then the local weak error of LMC algorithm satisfies
‖Ex̄1 − Exh‖ ≤ C̃1h2 with C̃1 = 2(L2 +G) ( d 4κL + E ‖x‖ 2 + 1 ) 1 2 .
Proof. By Ito’s lemma, we have
d∇f(xt) = −∇2f(xt)∇f(xt)dt+∇(∆f(xt))dt+ √ 2 ∫ t 0 ∇2f(xt)dBt.
It follows that
‖Ex̄1 − Exh‖ = ∥∥∥∥∥E ∫ h
0
∇f(xs)−∇f(x)ds ∥∥∥∥∥ = ∥∥∥∥∥∥E {∫ h 0 ∫ s 0 −∇2f(xr)∇f(xr) +∇(∆f(xr))drds+ √ 2 ∫ h 0 ∫ s 0 ∇2f(xr)dBrds
}∥∥∥∥∥∥ = ∥∥∥∥∥∥E {∫ h 0 ∫ s 0 −∇2f(xr)∇f(xr) +∇(∆f(xr))drds
}∥∥∥∥∥∥ ≤ ∫ h
0 ∫ s 0 E ∥∥∥∇2f(xr)∇f(xr)∥∥∥ drds+ ∫ h 0 ∫ s 0 E ∥∥∇(∆f(xr))∥∥ drds
≤L ∫ h
0 ∫ s 0 E ∥∥∇f(xr)∥∥ drds+ ∫ h 0 ∫ s 0 E ∥∥∇(∆f(xr))∥∥ drds
(i) ≤(L2 +G) ∫ h
0 ∫ s 0 E ‖xr‖ drds+ G 2 h2
≤(L2 +G) (∫ h 0 ∫ s 0 E ‖xr − x‖ drds+ h2 2 E ‖x‖ ) + G 2 h2
(ii)
≤ (L2 +G) (∫ h 0 ∫ s 0 √ E ‖xr − x‖2drds+ h2 2 E ‖x‖ ) + G 2 h2
(iii)
≤ (L2 +G) ∫ h 0 ∫ s 0 √ 6 ( d+ m 2 E ‖x‖2 ) rdrds+ h2 2 E ‖x‖ + G 2 h2
=(L2 +G) 4√6 15 √( d+ m 2 E ‖x‖2 ) h+ 1 2 E ‖x‖ h2 + G 2 h2
(iv) ≤ (L2 +G)h2 √( d+ m
2 E ‖x‖2
) h+ 1
2 E ‖x‖2 + G 2 h2
(v) ≤ (L2 +G)h2 √ d
4κL + E ‖x‖2 + G 2 h2
≤(L2 +G)
(√ d
4κL + E ‖x‖2 + 1
) h2
≤2(L2 +G) ( d
4κL + E ‖x‖2 + 1
) 1 2
h2
where (i) is due to Assumption 2, (ii) is due to Jensen’s inequality, (iii) is due to Lemma C.2, (iv) is due to √ a+ √ b ≤ √ 2 √ a2 + b2 and (v) is due to h ≤ 14κL . It is worth noting in the third equation that the Ito’s correction term ∇∆f can also be written as ∆∇f as the two operators commute for C3 functions.

heading: D.3 BOUNDEDNESS OF LMC ALGORITHM
text: Lemma D.3. Suppose Assumption 1 holds. Denote the iterates of LMC by x̄k. If 0 ≤ h ≤ 14κL we then have the iterates of LMC algorithm are uniformly upper bounded by
E ‖x̄k‖2 ≤ E ‖x0‖2 + 8d
7m , ∀k ≥ 0
Proof. We have E ‖x̄k+1‖2 =E ∥∥∥x̄k − h∇f(x̄k) +√2hξk+1∥∥∥2
(i) =E ‖x̄k‖2 + h2E ∥∥∇f(x̄k)∥∥2 + 2hd− 2hE〈x̄k,∇f(x̄k)〉 =E ‖x̄k‖2 + h2E
∥∥∇f(x̄k)−∇f(0)∥∥2 + 2hd− 2hE〈x̄k,∇f(x̄k)〉 (ii)
≤ E ‖x̄k‖2 + h2L2E ‖x̄k‖2 + 2hd− 2hE〈x̄k,∇f(x̄k)〉 (iii)
≤ E ‖x̄k‖2 + h2L2E ‖x̄k‖2 + 2hd− 2mhE ‖x̄k‖2
(iv) ≤ (
1− 7 4 mh
) E ‖x̄k‖2 + 2hd
where (i) is due to the independence between ξk+1 and x̄k, (ii) is due to Assumption 1, (iii) is due to the property of m-strongly-convex functions, 〈∇f(y)−∇f(x),y−x〉 ≥ m ‖y − x‖2 ∀x,y ∈ Rd, and (iv) uses the assumption h ≤ 14κL . Unfolding the inequality, we obtain
E ‖x̄k‖2 ≤ (1− 7
4 mh)kE ‖x̄0‖2 + 2hd
( 1 + 7
4 mh+ · · ·+ (7 4 mh)k−1
) ≤ E ‖x0‖2 + 8d
7m

Review summaries:
The manuscript considers the unadjusted Langevin Monte Carlo (LMC) algorithm, and performs non-asymptotic analysis of its convergence with respect to the 2 Wasserstein distance. The main contribution is a mixing time bound of O(d^0.5/\epilson), which improves upon the existing O(d/\epsilon) bound of Durmus and Moulines. The O(d^0.5/\epilson) is also shown to be optimal for the classes of probability distributions considered. 

The paper is concerned with the non-asymptotic analysis of SDE-based sampling algorithms and has two main contributions. Firstly, it improves upon the general framework of [Li et al. (2019)](https://proceedings.neurips.cc/paper/2019/hash/7d265aa7147bd3913fb84c7963a209d1-Abstract.html) and, in particular, does not require uniform boundedness assumptions for the local strong/weak errors. This means their first result, Theorem 3.3, is easier to apply to numerical schemes for contractive SDEs and can already be seen to unify previous works (such as the analysis of KLMC in [Dalalyan et al. (2020)](https://projecteuclid.org/journals/bernoulli/volume-26/issue-3/On-sampling-from-a-log-concave-density-using-kinetic-Langevin/10.3150/19-BEJ1178.short)). Secondly, the authors consider the Unadjusted Langevin Algorithm (ULA) and derive a 2-Wasserstein bound of $\mathcal{O}(\sqrt{d} h)$, which has optimal dependence on both the dimension $d$ and step size $h$. To achieve this, they impose an additional smoothness condition on $f$ that is less restrictive than having the Hessian $\nabla^2 f$ Lipschitz continuous. The authors support this theory with numerical examples, where they demonstrate that a lower bound on the 2-Wasserstein error scales as $\mathcal{O}(\sqrt{d})$ and $\mathcal{O}(h)$.