A Dynamical System View of Langevin-Based Non-Convex Sampling
Authors: Mohammad Reza Karimi Jaghargh
Ya-Ping Hsieh
Andreas Krause
Venue: NeurIPS 2023 spotlight (2023)
Keywords: Non-Convex Sampling
Langevin Dynamics
Dynamical Systems
TL;DR: We develop a novel framework that guarantees last-iterate convergence in Wasserstein distance for many advanced sampling algorithms at once.

Abstract: Non-convex sampling is a key challenge in machine learning, central to non-convex optimization in deep learning as well as to approximate probabilistic inference. Despite its significance, theoretically there remain some important challenges: Existing guarantees suffer from the drawback of lacking guarantees for the last-iterates, and little is known beyond the elementary schemes of stochastic gradient Langevin dynamics. To address these issues, we develop a novel framework that lifts the above issues by harnessing several tools from the theory of dynamical systems. Our key result is that, for a large class of state-of-the-art sampling schemes, their last-iterate convergence in Wasserstein distances can be reduced to the study of their continuous-time counterparts, which is much better understood. Coupled with standard assumptions of MCMC sampling, our theory immediately yields the last-iterate Wasserstein convergence of many advanced sampling schemes such as mirror Langevin, proximal, randomized mid-point, and Runge-Kutta methods.

A Dynamical System View of -anse Non-Convex Sampling
Mohammad Reza Karimi'
Ya-Ping Hsieh" ETH Zürich yaping hsieh@inf ethz. ch
Andreas Krause ETH Zürich
ETH Zürich mkarimi @inf. ethz. ch
krausea@ethz.co ch
Abstract
Non-convex sampling is a key challenge machine learning, central to non-convex optimization in deep learning well to approximate probabilistic inference. Despite its significance, theoretically there remain some important challenges: Existing guarantees suffer from the drawback of lacking guarantees for the last- iterates, and little known beyond the elementary schemes of stochastic gradient Langevin dynamics. To address these issues, we develop novel framework that lifts the above issues by harnessing several tools from the theory of dynamical systems. Our key result is that, for a large class of state-of-the-art sampling schemes, their last-iterate convergence in Wasserstein distances can be reduced to the study of their continuous-time counterparts, which much better understood. Coupled with standard assumptions of MCMC sampling, our theory immediately yields the last-iterate Wasserstein convergence many advanced sampling schemes such mirror Langevin, proximal, randomized mid-point, and Runge-Kutt Kutta methods.
1 Introduction
Many modern learning tasks involve sampling from a nigh-dimensional density æ or e- where f is non-convex potential representing, for instance, the loss function of deep neural network. To this end, an approach that has found wide success is to discretize the continuous-time Langevin diffusion
dL, (L,) dr v2 dB,
(LD)
where B, is a Brownian motion [57]. The idea behind this approach that, since is the stationary distribution of (LD), one can expect similar behavior for discretizations of (LD). Such framework has inspired numerous sampling schemes with per-iteration costs as cheap as stochastic gradient descent, which are particularly suitable for large-scale approximate probabilistic inference and Bayesian learning [2, 54, 57]. Moreover, several works have noticed that these Langevin-baat schemes provide deep insights about minimizing ff using stochastic oracles [22, 48], which serves as an important step toward explaining the empirical success of training deep neural networks.
The convergence of Langevin-basse based non-convex sampling has therefore attracted significant interest from both practitioners and theoreticians, whose intense study has led to plethora of new guarantees; see related work for details. Despite such impressive progress, several challenges remain for the fully non-conven setup:
The convergence is typically given on the averaged iterates instead of the more natural last iterates [4, 54]. This is especially problematic from the perspective of understanding the minimization of f, as in practice, the last iterates of an optimization algorithm play the most pivotal role for downstream tasks.
"Equal contribution.
37th Conference on Neural Information Processing Systems (NcurIPS 2023).

---

An additional notable drawback of the current theory is its predominant focus on the basic Euler- Maruyama discretization of (LD) (see, e.g., [4, 20, 54]). As result, the convergence analysis of more advanced sampling schemes remains largely unexplored in the fully non-convex regime [1, 24, 25, 34, 53, 61].
ş Contributions and Approaches. To overcome the aforementioned challenges, our main contri- bution, from a high level, can be succinetly summarized as:
Under mild assumptions, we prove that the iterates of a broad range of Langevin-based (*)
sampling schemes converge to the continuous-fime (LD) in Wasserstein distance.
Combining (*) with classical results on Langevin diffusion [45] immediately yields the last- last-iterate convergence in Wasserstein distances for wide spectrum of sampling schemes, thus resolving all the challenges mentioned above. To illustrate this point, we state simple version our main result.
Theorem (Informal). Suppose we discretize (LD) as
Xk+1 Xx ( f (xx) + noise + bias) + 22 Ek+1
with step-sizes kEN and i.i.d. standard Gaussians (Ex}REN- Then, under easy-to-verify condition on the bias (see (5) in Assumption 3), {xe)KeN converges in Wasserstein distance to n. In addition, these conditions are satisfied by many advanced sampling schemes.
This result is achieved via a new dynamical perspective to study Langevin- based sampling. More specifically,
We introduce the Picard process, which the sampling analogue of Picard's method of successive approximations for solving ODEs [17]. Contrary to most existing analyses, the Picard process allows us to completely bypass the use of relative entropy, which is the culprit for the appearance of averaged iterates [20].
Using the Picard process, we will prove that the iterates of various Langevin- basedd schemes generate a so-called Wasserstein asymptotic pseudotrajectory (WAPT) for the continuous- time (LD). The main motivation for considering WAPT to connect Langevin- schemes to the dynamical system theory of Benaîm and Hirsch [7], which works for metric spaces and is last-iterate by design, and therefore particularly suitable for our purpose.
3. Finally, under standard stability assumptions in the literature [39, 51], we show how a tandem of WAPT result and dynamical system theory yields the desirable convergence of various existing schemes, as well as motivates more efficient algorithms that enjoy the same rigorous guarantees.
Related work. There vast literature on structured non-convex sampling, where one imposes extra assumptions on the target density. Under these conditions, one can derive non-asymptotic rates for Langevin- baased schemes [13, 15, 33, 35, 36, 41, 48, 56, 60, 63]. Our work is orthogonal to these works as we study generic non-convex sampling, an NP-hard problem whose convergence is asymptotic best.
Most relevant to our paper are the works [4, 8, 20. 29, 54], which study the asymptotic convergence of Langevin-based schemes under minimal regularity assumptions on f. Compared to their results, our guarantees either improve upon existing ones or are incomparable; see Section 5.4 for a more detailed comparison.
2 The Langerin-Robbin--Monn Template
We consider the following general template for sampling algorithms: Starting from initial point, the iterates {xx)een follow the recursion
X X ( Zu+13 <2Yk+1 (xxx) EL+1.
(LRM)
where yx 's are step sizes, is vector field, Zx's 's's are (random or deterministic) perturbations, or is the state-dependent diffusion matrix, and & 's are i.i.d. standard Gaussian random variables. In the sequel, we will further decompose the perturbation Z₃ Ux bg. where Ug is the (zero-mean)

---

noise and by is the bias. We call this recursion the -ieiii-RRobbbs---oon (LRM) template, as it is reminiscent of the Robbins-Monro template for stochastic approximation [50].
The generality of the LRM template allows us to capture many existing algorithms and suggests ways to design new ones. For illustration purposes, showcase instances of (LRM) with the following examples. Other examples (SGLD and proximal) are provided in Appendix A. In the first three examples, the vector field v in (LRM) is f and
Example 1. The Randomized Mid-Point Method [24, 53] is an alternative discretization scheme to Euler-Maruyama and has been proposed for both overdamped and underdamped Langevin diffusion. For the overdamped case, its iterates are
1/+1/2 = Xx YK+1 (xxx) + V2yk+1RR+1 E
(RMM)
Xk+1 = X2 Yk+1 V (x2+1112) + + EK+1,
where {(rx) are i.i.d. and uniformly distributed in [0, 1], Éx. are standard Gaussian random variables with cross-variance varg, and ff is noisy evaluation of f.. To cast (RMM) in the LRM template, we set UL+1 == fxxx+1122 (xx=1122 and PL+1 V / fxxx).
that it requires two gradient oracle calls at each iteration. Inspired by the optimistic gradient methods in optimization and online learning [16, 47, 49], we propose to "recycle" the past gradients:
Example 2. Inspecting the update rule of (RMM),
Xk+1/2 = Xx Yk+1 ax+1 a - (21t-104H151H
(ORMM)
X++1 = Xx Yk+1 (xx+112
22 5k+1
where Ex. and f are the same as (RMM). This is again an LRM scheme with U2+1 == (((++1/2) - V [x(++1/)2) and bx+1 = (*x++122 f(xx).
Notice that (ORMM) requires one gradient oracle, thereby reducing the per-iteration cost of (RMM) by 2. To our knowledge, the scheme (ORMM) is new.
Example In addition to the simple (stochastic) Fuler-Maruyama discretization in (SGLD), there exists a class of more sophisticated discretization methods of (LD) known as higher-order integrators. Runge-Kutta Kutta method [34] is an example of an order 1.5 integrator, with iterates
The Stochastic Runge-
h = xg + v 22kk (CIEx+ + c2fk.1)
h2 xx - 7k+1 V fxx) /2ykr (C3Ex+1+ c25
X
7 (h:) (fz))) 2y Ek+1,
where Ex+1 and E are independent standard Gaussian random variables, and c1. are suitably chosen integrator constants. This algorithm is LRM scheme with U2+1 == (( -V)())))(()))))) (V2) (h2)) and bx+1 = (hi) + f - fxx)).
Example 4. The Mirror Langevin algorithm [1, 25, 61], which is the sampling analogue of the celebrated mirror descent scheme in optimization [5, 43), uses strongly convex function & to adapt favorable local geometry. In the dual space (i.e., the image of Vd), its iterates follow
Xk+1 X Yk+1 V VVVVV'(VVVV + V2 V2 (xL) (1))) EL+1.
(ML.
where 8' is the Fenchel dual of ø [52]. In our framework, (ML) fits into (LRM) by taking f V and o = 2 1-122 Additionally, one can also consider a stochastic version of (ML) with noisy evaluations of f.
3 Technique Overview: A Dynamical System Perspective
The goal of our paper is to provide last-iterate guarantees for the general LRM schemes introduced in Section 2. There are two equivalent, commonly considered, ways of characterizing the dynamies of the iterates of an LRM scheme. The first one is to view the iterates (xx keN as a random trajectory in Rd, which is perhaps the most natural way of describing sampling algorithm. The second way to view the distributions (eK)keN of (xR a deterministic trajectory in the Wasserstein space.
3

---

With these two characterizations in mind, in this section, we will devise new framework based on the dynamical system theory and present its high-level ideas.
To understand our novelty, it is important to contrast our framework to the existing Wasserstein viewpoint towards Langevin-based sampling algorithms. Following the seminal work of Otto [44], one can view a sampling algorithm as the discretization of a class of well-studied dynamical systems-gradien flows. This viewpoint suggests using Lyapunov arguments, which has become the predominant approach in much prior work.
Despite its appealing nature, in the rest of this section, will argue that Lyapunov analysis of gradient flows is in fact not suited for studying generic non-convex sampling. In particular, we will show how our new framework is motivated to overcome the several important limitations of gradient flow analysis. Finally, we give a high-level overview of the techniques used in our paper.
$ Langevin Diffusion Gradient Flows. We denote by Pr the probability density of L, in (LD), and consider the continuous curve + p; in the Wasserstein space W2. In their seminal works, Jordan et al. [27] and Otto [44] discover that this curve is the (exact) gradient flow of the relative entropy functional; that is, defining the functional F p + (xI(pllle one has depr = grad F(e)), where "grad" is the gradient in the Wasserstein sense. This gradient flow viewpoint of (LDD thus provides a clear link between sampling in I and optimization in WW, Indeed, this suggests that the relative entropy is natural choice for the Lyapunov function of the discrete- time sampling algorithm. which is a prominent approach for analyzing sampling algorithms in recent years [4, 21, 58].
Although the gradient flow viewpoint has led to sequence of breakthroughs, it has few important shortcomings:
(a) The usual Lyapunov-type analysis for sampling algorithms focuses on bounding the change in relative entropy across iterations. This is extremely challenging when one considers more advanced sampling algorithms, as one has to understand the effect of the additive bias and noise of the algorithm on the change of relative entropy. Crucially, this makes the Lyapunov analysis applicable only to the simple Euler-Maruyama discretization of (LD), i.e., X&+1 X - Yk+1 V fx(x) v2yx+1 EK+1, and fails to capture more advanced and biased sampling schemes such as Examples 1-4. Even for the simple (SGLD), the presence of stochastic gradients significantly complicates the Lyapunov analysis and requires extra assumptions such as convexity [21 ] or uniform spectral gap [48].
(b) This gradient flow-based analysis often requires an extra averaging step to decrease the relative entropy (see, e.g., [4]). This is the main reason why many existing works provide guarantees only on the averaged iterates (px == & Σ Pt) instead of the last ones (px).
In this paper, we overcome these limitations by introducing a new perspective, whose two ingredients are as follows.
Wasserstein Asymptotic Pseudotrajectories. A notion that will play a pivotal role our analysis is the Wasserstein asymptotic pseudotrajectory (WAPT), which is a measure of "asymptotic closeness" in the Wasserstein sense, originally defined by Benaim and Hirsch [7] for metric spaces:
Definition 1 (Wasserstein asymptotic pseudotrajectorry). We say the stochastic process (X,) is Wasserstein asymptotic psendotrajectory (WAPT) of the SDE
dΦ, = u( (() ((D,) dB,
(SDE)
if, for all 7' O,
lim sup W2 l.:)) 0. OS8sT
Here, 0.0) is the solution of the SDE time initialized X,, and W2 is the 2-Wasserstein distance.
Despite the seemingly convoluted definition, WAPT can be intuitively understood as follows: Let {xx keN be the iterates of sampling scheme. Then, (1) simply posits that for sufficiently large
While the Lyapunov-type analysis has been applied to elementary (i.c., unbiased) discretization schemes for other SDEs, such as the under-damper (i.e., kinetic) Langevin dynamics [19], our primary focus in this paper remains centered on the overed Langevin diffusion and similar SDEs.

---

Dissigutivity (Thenrem
Dynamies via WAPT (Thorem I)
Dynamical Perspective
Stability (Eq. (10))
Dissipativity
(Theoremn 4)
Last Iterate Convergence in Wi (Theorem
Figure 1: High-level overview of the two components of the dynamical perspective.
m, one cannot distinguish between the "tail" iterates {xx kmm versus the SDE solution starting at Xm up to arbitrarily small error measured in terms of the Wasserstein distance. Since we are only interested in the asymptotic behavior of Xj, these controls on the tail iterates will suffice to conclude the last-iterate convergence."
Importantly, from the perspective of WAPT, the Langevin diffusion (L.D) (or more generally, ")")) is simply viewed as a generic dynamical system and not as gradient flow. In particular, relative entropy will play no role throughout our analysis, thereby resolving issue (b).
$ Langevin- Langevvn-Robbins--eee Re---Roobbbborrrr Schemes. We have seen that the LRM template in Section 2 is capable of capturing broad range of existing and new algorithms in unified way. To resolve the remaining issue (a), we will further rely on the LRM template: for proving that (LRM) generates WAPT of the corresponding SDE, we show that the key condition (1) WAPT can be reduced to checking casy-to-verify bound on the perturbation terms Zk.
To achieve this, the most important step in our proof, which distinguishes our analysis from all existing works in non-convex sampling, is the construction of the so-called Picard process, the natural generalization of the Picard's successive approximation method [17] from ordinary differential equations to stochastic differential equations. In the stochastic approximation literature, similar techniques have been successfully applied to study optimization and games various settings such as on Riemannian or primal-dual spaces [26. 28, 37J. The application to sampling has also been previously explored by Bubeck al. [11], Chau et [12] different contexts. What distinguishes our work from the existing literature is the advantage of generalizing the Picard process encompass a vastly wider class of algorithms, specifically the LRM schemes. Moreover, the integration of the Picard process with the theory of WAPT plays pivotal role in our analysis, and both of these aspects present original contributions.
$ Framework overview. To conclude, for proving last-iterate convergence, proceed as follows:
1. For a given LRM scheme (xx IReN, we first construct a continuous trajectory (X,)rz0 via interpolating the iterates (see (3)).
2. We prove that (X,) constitutes WAPT of the SDE (see Theorem 1). This step relies heavily the construction of the aforementioned Picard process.
3. By invoking the dynamical system theory of Benaim and Hirsch [7]. the convergence of LRM
schemes reduces simply checking the stability condition (Theorem 2). the Wasserstein space, this condition translates into boundedness of the second moments of the iterates (, for which there is plethora of approaches; we present two such methods in Section
Fig. 1 depicts a high-level overview of the ingredients needed in our framework, and their correspond- ing theorems.
The Dynamies of Langevin-Robinn-MMon Schemes
In this section, we view (LRM) noisy and biased discretization of (LD). To make this analogy precise, let (B,)a0 be a Brownian motion defined on a filtered probability space with filtration
Definition is phrased in terms of continuous-tine time stochastic process (X,),20- The discrete iterates (xR ken be converted to continuous -immoeuusttmeeeee process through a suitable interpolation; (3).
5

---

(Fi) and define тą -n= Yn to be the effective time that has elapsed at the iteration k. Using the Brownian motion, we can rewrite (LRM) as
X1+1 X Yk+1 (u(xx) ZR+1) (xx)) (Bzz.: BTT))
assuming that the filtration satisfies Zx F 4 The (continuous-time) interpolation (X))r20 of hee is then defined as the adapted process
x, = x (1 - (x(x)) + E[ 9:1) + (xx) (B, BB,), for E [Tk, Tk+1].
(3)
In addition, for fixed 1, consider the Brownian motion (B.")s )>0 where B(d) = Br+s B, and define the Langevin flow (T)") )=0 as the (strong) solution of (SDE) initialized at X, It is important to note that ((() and X are synchronously coupled by sharing the same Brownian motion.
4.1 Technical Assumptions and Requirements
We now introduce the basic technical assumptions and discuss their generality. Assumption 1. The vector field is L-Lipschitz, and satisfies (x, v(x)) s Co(1 + [|x]] for some C O. Moreover, or is L-Lipschitz and bounded
norm.
I Lipschitzness of is standard assumption and is also required to ensure the existence of a unique strong solution of (SDE). The second assumption on the vector field is exceedingly weak and when u - is satisfied even for distributions without moments. The assumptions on diffusion coefficient o are already satisfied when o = 1, and we show that it holds for practical schemes such as Example 4.
Assumption 2. The Bobbins-Monru summability conditions hold: E Yt oo and -k=1 y co.
Moreover, for some constant P to be defined in (20), we have
Yk+//k Pyky
Yk:
Vk.
The Robbins-Monro step size conditions are standard in the non-convex sampling literature [4, 20, 29, 30]. For (4), can be verified that condition is satisfied even for slowly-decreasing step sizes such as Yk (VI log k)- which hence is not restrictive.
Assumption 3. The noises (Ux)keN form martingale difference sequence, i.e., E[Ux+1 Uv] and have uniformly bounded second moments. In addition, the bias terms satisfy
E[]lbb+I|| 1| Fr. = 1(xx))llu(xx)11++++++ yxt)).
(5)
martingale difference sequence more general than an i.i.d. sequence, allowing the noise to be dete-dependent. The bias condition (5) simply states that the bias shall not overpower the signal v(xx), and, as we show later, is satisfied by all our examples.
4.2 From Discrete Continuous: LRM Schemes and WAPTs
We are now in position to state our main theorems. Our first result below establishes a precise link between the discrete-time (LRM) and the continuous-time (SDE).
Theorem 1. Under Assumptions 1-3, the interpolation (3) of an LRM scheme is a Wasserstein asymptotic pseudotrajectory of (SDE).
Sketch of the Proof for Theorem 1. The proof of this theorem is heavily based the notion of the Picard process and iterate moment bounds. The complete proof can be found in Appendix C.
Step The Picard Process. For a fixed > recall the construction of the interpolation (3) and the Langevin flow. Central to our analysis is the Picard process, defined as
X, + 1/" v(X(uu) du + ( (((++u) dB!").
(6)
*One augment the filtration of the Brownian motion the olaagerra of Zg times (Tk IeeN-

---

The Picard process is adapted and is (synchronously) coupled with the Langevin flow and the interpolation. We think of the Picard process as one step of the Picard iteration for successive approximations to solve ODEs. This means, intuitively, that its trajectory should be close to the original interpolation, as well to that of the Langevin flow, playing the role of "bridge".
Fix T > O. For E [], we decompose the distance between the interpolation X, in (3) and the Langevin flow as
|1 101 s |r!! 11 + || - 1!
We now bound each term of the decomposition. By synchronous coupling of the processes, Lipsehitz- ness of and Itô isometry, Lemma 3 bounds the first term
S [100"" X du.
| - - , s 2(T 1))?
(8)
This will be suitable for later use of Grönwall's lemma.
$ Step 2: Accumulated Noise and Bias. For the rest of the proof, we need some extra notation.
Define m(r) sup 0: T2 & f and the piecewise-constant process X == Xm(()) Going back to the second term of (7), observe that
J," -x
v(Xu) du + S'"
Xres Y
dB (z,
where (t, s) is the accumulated noise and bias from time I to time + s. expected that (I.z(I.)) ()]l eventually becomes negligible, since the step size becomes small. The next lemma confirms this intuition.
Lemma 1. Suppose Assumptions 1-3 hold. Then, for any fixed we have
lim sup Os5ST
(t,
Step 3: Gradient Moment Bounds. Based on (9) and Lemma 1, bounding the distance between the Picard process and the interpolation essentially reduces to bounding how much the discrete algorithm "moves" during one iteration in expectation. This, in turn, depends on how large the moments of [|{(x))] grow per iteration, which is controlled by the following lemma:
Lemma Let (xx keN be the iterates of (LRM) and suppose Assumptions 1-3 hold. Then, E||xx|
O(1/yk=1) This in turn implies Ellu(x))] (1//r++1 and El]bx+1l12 = (YKkt))
Using this lemma and Lemma we can obtain A, == suPo E||X,ss 11 - 0 as f - co, which shows that the Picard process gets arbitrarily close to the interpolation as f - co.
Step 4: Concluding the Proof. Let us go back to the dccomposition (7). Taking expectation and using (8) and Grönwall's lemma, we obtain E[]]X 1D/"113 s 4 A, exp( Thus,
lim sup sE[O.T]
E X ".."'' =.
As we coupled X,+x and ( in specific way (via synchronizing the Brownian motions), we directly get an upper bound on the Wasserstein distance.
5 Last-Iterate Convergence of Sampling Schemes
In this section we focus on last-iterate convergence of LRM schemes in Wasserstein space. We first explore the interplay between the convergence of WAPT' and stability. We then show that the existing stability results for simple Fuler-Maruyana discretization of the Langevin diffusion can be extended, with little to no extra assumptions, to the class of LRM schemes in Section This in turn readily implies the last-iterate convergence of wide class of LRM schemes.
7

---

5.1 From WAPT's to Convergence in W2
Since convergence of the distribution of xx to x in Wasserstein distance implies convergence of the second moments of xę to that of [3], convergence in the Wasserstein space should at least require:
sup
E|| co.
(10)
keN
turns out that, for WAPTs, the exceedingly weak necessary condition (10) is also sufficient:
Theorem 2. Let (X₂) be a Wasserstein asymptotic pseudotrajectory of the Langevin diffusion (LD)
generated by LRM scheme (xx) via (3). Then (2(x,, - 0 if and only if (10) holds.
Proof. The proof relies on the structure of compact sets in the Wasserstein space and limit-set theorems for dynamical systems [7]. Specifically, the closure of bounded subsets of W2 is compact [3], condition (10) implies that (law(X))) is pre-compact in W2- Moreover, Assumption 1 implies that the Langevin flow is globally integrable. Thus, (law(X))) x2))rz0 is a pre-compact WAPT of globally integrable flow, and we can apply the limit-set theorem for metric spaces [7, Theorem 0.1]
to conclude that the limit-set of (law(X) (X,), is an internally chain transitive (ICT) set.
Next, show that for the case of the Langevin flow, the only ICT set is (π), implying the desired convergence of our theorem. To see this, define V(-) = D x). It can be observed that V is a Lyapunov function for (LD), whose value is strictly decreasing along the flow (as the time derivative of V along the flow is negative of the relative Fisher information, which is strictly positive for all measures other than x). Thus, all requirements of [6. Prop. 6.4] satisfied, showing that the only point in the ICT set is π. This also shows the uniqueness of the stationary distribution of (LD).
Remark. From the proof of Theorem we observe that the supremum of the Wasserstein distance between [t,tT]] and (T!!)) [o.]] typically scales exponentially with T, which is common for weak approximation error in the literature, see [40]. Despite the exponential dependence on T, the
convergence of the last iterate is assured by Theorem 2 without need of uniform control in T. This is primarily attributed to the adoption of a dynamical system viewpoint and the application of corresponding tools, effectively harnessing the paradigm established by Benaïm and Hirsch.
Theorems 1-2 in tandem thus show that, as long as LRM scheme satisfies Assumptions 1-3 and the moment condition (10), the desirable last-iterate convergence in W/2 is immediately attained. Therefore, in the rest of this section, we turn our focus establishing (10) for LRM schemes.
5.2 Bounded Moments of LRM Schemes
There is long history of study on conditions that ensure (10) for iterative algorithms, which has culminated in the so-called dissipativity properties. We consider two such examples below.
Assumption (Dissipativityy. There exist constants a and B B20 such that
(x,v(x)) s - + B, Vx E Rd.
Under Assumption 4, it is classical that (10) holds for the simple Euler Maruyaaa discretization of (LD) with deterministic or stochastic gradient oracles [23, 29, 30, 38, 48, 51, 54]. These studies, however, cannot handle non-zero bias, which, as seen in Examples 1-3, is crucial for incorporating more advanced sampling schemes.
To this end, our next result shows that for wide class of LRM schemes, the stability (10) essentially comes for free under Assumption 4. The proof is provided in Appendix D.
Theorem 3. Let be vector field satisfying Assumptions I and 4 and or be diffusion coef- ficient satisfying Assumption I, and let (XL) be an LRM scheme. Assume that limg Yk = sup, 2| co, and the bias satisfies (5). Then, the stability condition (10) holds for {xz}.
weaker notion of dissipativity that has been studied in the literature is:
Assumption (Weak dissipativity). There exist constants a > 0, K E (0. 1], and B 2 such that
vx,) -a||x|| 1+x +B. Vx Vx.RS.
8

---

NOISE
BIAS
LAST-ITERATE
LAMBERTON AND PAGES [29]. LEMAIRE [30]
TEH ET AL. [54]
BENAÎM AL.
DURMUS AND MOULINES [20]
BALASUBRAMANIAN [4]
THIS WORK
Table I: Comparison to existing works on convergence of L.R schemes, All methods, except for [4], require bounded second moments of the iterates.
When K Assumption 5 is simply Assumption 4. As opposed to Assumption 4, which requires quadratic growth of f outside a compact set (when u = f),,, Assumption 5 only entails superlinear growth and therefore is considerably weaker.
For Suler-Maruyami discretization of (LD) with deterministic gradients, [20] prove that Assumption sufficient to guarantee bounded moments of the iterates. As for generic LRM scheme, consider the following general condition on the bias terms, which will suffice to cover all our examples in Section 2: For some constant c,
|| |b & c((r+, (1(x))1 12.1 Yk+1 1 yenillee-l 11+111x++111)))
(11)
where U/,1 is an extra noise term, and Ex,, is standard Gaussian independent of the noises and Ek- The price to pay with the weaker Assumption 5, however, is that we need to assume uub-Gaussianity of the noise. For a proof, see Appendix D.
Theorem 4. Let and let {xx} \ be an LRM scheme. Assume that lim Yk: = 0, the noises Us and Uú are sub- Gaussian, and the bias term of {xx} satisfies (11). Then, (10) holds for {xg] in when (i) 1, or (ii) f is Lipschitz and the LRM follows the Mirror Langevin algorithm (Example 4).
f be the targer distribution, where u f satisfies Assumptions I and 5,
5.3 Examples of Convergent LRM Schemes
We now illustrate the of Theorems 1-4 on examples in Section 2.
Proposition 1. Under Assumption and noise with uniformly bounded second moments, the following holds for Examples 1-6: (i) The bias has the form (11) and satisfies (5), (üi As result, under Assumptions 2 and 3, Examples 1-6 produce iterates that generate a WAPT of (SDE). (ii)) Under the additional conditions Theorem 3 or Theorem 4, Examples 1-6 enjoy last-iterate convergence the target distribution in Wasserstein distance.
5.4 Comparison to Existing Work
We now give more detailed comparison of our results existing literature; summary is given in Table 1, and additional comparison with prior works can be found in Appendix B.
$ Guarantees for LRM Schemes. Lamberton and Pages [29] and Lemaire [30] study the simple Mullr-Maruyama discretization of (LD) with deterministic gradients (i.e., U2 = bx = 0) and establish the weak convergence of the average iterates under moment condition that is slightly weaker than (10). Their analysis is further extended by [54] to incorporate stochastic gradients. Later, the last-iterate convergence of the simple Fuler-Maruyama discretization of (LD) is studied by [20], who prove the convergence in the total variation distance under Assumption 5. Another work on a similar setting as [20] is [8]. where the convergence criterion is given in an integral probability metric (IPM) [42] of the form dæ ( r) == se E E for certain class of test functions B that is known imply weak convergence, but not convergence in total variation or Wasserstein distances.
Compared to these results, our guarantees possess the following desirable features:
5 Although the condition in [29, 30) is stated in weaker form than (10), it is typically only verified on special case that is equivalent Assumption and thus implies (10). See e.g.. [29, Remark 3).

---

The convergence is always on the last iterates instead of the average iterates.
As we tolerate biased algorithms, the class of LRM schemes we consider is significantly more general than the ones in existing work.
Finally, we note that our results are incomparable to the recent work of Balasubramanian et al. [4], who derive the same result as in [29, 30J, i.c., average-iterate, weak convergence, deterministic Euler-Maruyama discretization. A remarkable feature of the analysis in [4] is that it does not require any bounded moments, and, in particular, their bounds can be applied to target distributions with unbounded variance. However, the downside of [4] is that, in the presence of U3 and bx, their analysis produces a bound that does not vanish as - oo; see [4, Theorem 15]. In contrast, our framework can tolerate quite general U2 and bx. gives stronger guarantees (W2 vs. weak convergence; last-iterate vs. average-iterate)
S On Analysis Techniques. While, to our knowledge, our framework is significantly different from previous works on sampling, we acknowledge that similar ideas of creating an auxiliary process in-between the iterates and the continuous-time flow is not entirely new and has been touched upon in the literature, e.g., [10, 12]. That being said, our specific approach in building the Picard process and its development into wider array of algorithms, i.e., Langevin-Robbee undoubtedly plays a pivotal role in our analysis. Moreover, the integration of the Picard process with offers dual benefits to study, and view these as
Monro schemes,
the theory of asymptotic eseudo-trajectore our unique contributions to this area of research.
Furthermore, the novel Picard process gives significant advantage in all of our results. The work of [8] also hinges on dynamical system theory-related ideas. Yet, missing the critical step of the Picard process has seemingly resulted in much weaker findings compared to our work. This observation not meant as a critique; rather, it merely highlights the potency of the unique method we have integrated into our study.
6 Concluding Remarks
In this paper, we provided a new, unified framework for analyzing a wide range of sampling schemes, thus laying the theoretical ground for using them in practice, well motivating efficient sampling algorithms that enjoy rigorous guarantees. We built on the ideas from dynamical system theory, and gave rather complete picture of the asymptotic behavior of many first-order sampling algorithms. In short, our results help with the following:
and more
Validating existing methods: Methods like mirror Langevin and randomized mid-point currently lack even asymptotic guarantees in fully non-convex scenarios, such as sampling from neural network-defined distributions. Our work fills this gap by offering the first rigorous justification for these schemes, supporting practitioners in utilizing these methods confidently.
Facilitating new algorithm design: Our work motivates novel sampling methods through a straightforward verification of Assumptions 1-3. An illustrative instance involves the randomized mid-point method and Runge-Kunn Kutta integrators, wherein a substantial 50% reduction in compu- tation per iteration can be achieved without compromising convergence by simply recycling past gradients, shown in Example 2. The balance between the benefits of saving gradient oracles and potential drawbacks remains open question, necessitating case-by-case practical evaluation. Nevertheless, our theory provides flexible algorithmic design template that extends beyond the current literature's 's scope.
While our WAPT result holds under very mild conditions, a severe limitation of our current framework that it only applies to Langevin-based algorithms, whereas there exist numerous practical sampling schemes, such as Metropolis-Hasting restriction arises as an artifact of our analysis, as the WAPT framework can principle applied equally well to any continuous-timg time dynamics. Lifting such constraint is an interesting future work.
that are not immediately linked to (LD). We believe that this
10

---

Acknowledgments and Disclosure of Funding
This work was supported by the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program grant agreement No 815943. YPH acknowledges funding through an ETH Foundations of Data Science (ETH-FDS) postdoctoral fellowship.
References
[1] Kwangjun Ahn and Sinho Chewi. Efficient constrained sampling via the mirror-langeve algorithm. Advances in Neural Information Processing Systems, 34:28405-2885 284188, 2021.
[2] Sungjin Ahn, Anoop Korattikara, and Max Welling. Bayesian posterior sampling via stochastic gradient fisher scoring. In Proceedings of the 29th International Coference on International Conference Machine Learning, pages 1771-1778, 2012.
[3] Luigi Ambrosio, Nicola Gigli, and Giuseppe Savaré. Gradient flows: the space of probability measures. Springer Science & Business Media, 2005.
metric spaces and in
[4] Krishna Balasubramanian, Sinho Chewi, Murat A Erdogdu, Adil Salim, and Shunshi Zhang. Towards theory of non-log-coneave sampling: first-order stationarity guarantees for langevin monte carlo. In Conference Learning Theory, pages 2896-2922 PMLR, 2022.
[5] Amir Beck and Mare Teboulle. Mirror descent and nonlinear projected subgradient methods for convex optimization. Operations Research Letters, 31(3): 2003.
[6] Michel Benaïm. Dynamics of stochastic approximation algorithms. In Seminaire de probabilites XXXIII, pages 1-68. Springer, 1999.
[7] Michel Benaim and Morris W Hirsch. Asymptotic pseudotrajectories and chain recurrent flows, with applications. Journal of Dynamics and Differential Equations, 8(1):141-176, 1996.
[8] Michel Benaïm, Florian Bouguet, and Bertrand Cloez. Ergodicity of inhomogeneous markov chains through asymptotic pseudotrajectories, The Annals of Applied Probability, 27(5):3004- 3049, 2017.
[9] Espen Bernton. Langevin monte carlo and jko splitting. In Conference On Learning Theory, pages 1777-1798. PMLR, 2018.
[10] Sébastien Bubeck, Ronen Eldan, and Joseph Lehec. Sampling from log-concave distribution with projected langevin monte carlo. arXiv preprint arXiv: 1507.0250 2015.
[11] Sébastien Bubeck, Ronen Eldan, and Joseph Lehec. Sampling from log concave distribution with projected langevin monte carlo. Discrete Computational Geometry, 59:757-783, 2018.
[12] Ngoc Huy Chau, éric Moulines, Miklos Rásonyi, Sotirios Sabanis, and Ying Zhang. On stochastic gradient langevin dynamics with dependent data streams: the fully non-convex case, 2021.
[13] Xiang Cheng, Niladri S Chatterji, Yasin Abbasi- Yadkori, Peter L Bartlett, and Michael I Jordan. Sharp convergence rates for langevin dynamics in the nonconvex setting. arXiv preprint 105 2018.
[14] Sinho Chewi. Log-coneave sampling, 2023.
[15] Sinho Chewi, Murat Erdogdu, Mufan Bill Li, Ruoqi Shen, and Matthew Zhang. Analysis of langevin monte carlo from poincar\' 'e to log-sobolev. arXiv preprint arXiv:2 2021.
[16] Chao-Kai Chiang, Tianbao Yang, Chia-Jung Lee, Mehrdad Mahdavi. Chi-Jen Lu, Rong Jin, and Shenghuo Zhu. Online optimization with gradual variations. In Conference on Learning Theory, pages 6-1. JMLR Workshop and Conference Proceedings, 2012.
[17] Earl A Coddington and Norman Levinson. Theory of ordinary differential equations. Tata McGraw-Hin Education, 1955.

---

[18] Arnak S Dalalyan and Avetik G Karagulyan. User-friendly guarantees for the langevin monte carlo with inaccurate gradient. arXiv preprint arXiv:1710.00095, 2017.
[19] Arnak S. Dalalyan and Lionel Riou-Durand. On sampling from a log- log-concave density using kinetic Langevin diffusions. Bernoulli, 26(3):1956 1988, 2020. doi: 10.3150/19 131150/99-BE1111 178. URL https: org/10. wttpg//....///..../// 3150/19- BEJ1178.
[20] Alain Durmus and Eric Moulines. Nonasymptotic convergence analysis for the unadjusted langevin algorithm. The Annals of Applied Probability, 27()):551-1587, 2017.
[21] 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 2711,
[22] Murat Erdogdu, Lester Mackey, and Ohad Shamir. Global non-convex optimization with discretized diffusions. Advances in Neural Information Processing Systems, 31, 2018.
[23] J Hale. Asymptotic behavior of dissipative systems. In Aerican Mathematical Society, 1988.
[24] Ye He, Krishnakumar Balasubramanian, and Murat A Erdogdu. On the ergodicity, bias and asymptotic normality of randomized midpoint sampling method. Advances in Neural Informa- tion Processing Systems, 13.7366-7376, 2020.
[25] Ya-Ping Hsieh, Ali Kavis, Paul Rolland, and Volkan Cevher. Mirrored langevin dynamics. Advances in Neural Information Processing Systems, 31, 2018.
[26] Ya-Ping Hsieh, Panayotis Mertikopoulos, and Volkan Cevher. The limits of min-max optimiza- tion algorithms: Convergence to spurious non-critical sets. In International Conference on Machine Learning, pages 4337-4348, PMLR, 2021.
[27] Richard Jordan, David Kinderlehrer, and Felix Outo. The variational formulation of the fokker- planck equation. SIAM journal on mathematical analysis, 29(1): 1998.
[28] Mohammad Reza Karimi, Ya-Ping Hsich, Panayotis Mertikopoulos, and Andreas Krause. The dynamics of riemannian robbins-monro algorithms. arXiv preprint arXiv:2206.0677 2022.
[29] Damien Lamberton and Gilles Pages. Recursive computation of the invariant distribution of diffusion. Bernoulli, pages 367-405, 2002.
[30] Vincent Lemaire. Estimation récursive de la mesure invariante d' un processus de diffusion. PhD thesis, Université de Marne la Vallée, 2005.
[31] Ruilin Li, Molei Tao, Santosh S. Vempala, and Andre Wibisono. The mirror langevin algorithm converges with vanishing bias, 2021.
[32] Ruilin Li, Hongyuan Zha, and Molei Tao. Sqrt(d) dimension dependence of langevin monte carlo. In The International Conference on Learning Representations, 2022.
[33] 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.
[34] Xuechen Li, Denny WI, Lester Mackey, and Murat A. Erdogdu. Stochastic Runge- Kutta Accelerates Langevin Monte Carlo and Beyond, February 2020.
[35] Yi-An Ma, Niladri S Chatterji, Xiang Cheng, Nicolas Flammarion, Peter I. Bartlett, and Michael I Jordan. Is there an analog of nesterov acceleration for gradient-based meme? Bernoulli, 1749)1942- 2021.
[36] Mateusz B Majka, Aleksandar Mijatović, and Lukasz Szpruch. Nonasymptotic bounds for sampling algorithms without log-concavity. The Annals of Applied Probability, 35(4):1534 1581, 2020.
[37] Panayotis Mertikopoulos, Ya-Ping Hsieh, and Volkan Cevher. Learning in games from stochastic approximation viewpoint. arXiv preprint 223212200.03392,, 2022.
12

---

[38] Sean P Meyn and Richard L Tweedie. Stability of markovian processes ii:: foster-lyapune criteria for continuous-tin processes. Advances in Applied Probability. 25(3):518-544 548, 1993.
[39] Sean P Meyn and Richard L Tweedie, Markov chains and stochastic stability. Springer Science & Business Media, 2012.
[40] Grigori N Milstein and Michael V Tretyakov. Stochastic numerics for mathematical physics, volume 39. Springer, 2004.
[41] Wenlong Mou, Nicolas Flammarion, Martin J Wainwright, and Peter I. Bartlett. Improved bounds for discretization of langevin diffusions: Near-optimal rates without convexity. Bernoulli, 187))1157711611 2022.
[42] Alfred Müller. Integral probability metrics and their generating classes of functions. Advances in Applied Probability, 29(2)42294433 1997.
[43] AS Nemirovsky and DB Yudin. Problem Complexity and Method Efficiency in Optimization. J. Wiley @ Sons, New York, 1983.
[44] Felix Otto. The geometry of dissipative evolution equations: the porous medium equation. Communications in Partial Differential Equations, 26, 2001.
[45] Grigorios A Pavliotis. Stochastic processes and applications: diffusion processes, the Fokker- Planck and Langevin equations, volume 60. Springer, 2014.
[46] Marcelo Pereyra. Proximal markov chain monte carlo algorithms. Statistics and Computing, 26 (4):745-760, 2016.
[47] Leonid Denisovich Popov. A modification of the arrow-hurwicz method for search of saddle points. Mathematical notes of the Academy of Sciences of the USSR, 28(5):845-88) 848, 1980.
[48] Maxim Raginsky, Alexander Rakhlin, and Matus Telgarsky. Non-convex learning via stochastic gradient langevin dynamics: a nonasymptotic analysis. In Conference on Learning Theory, pages 1674-1703. PMLR, 2017.
[49] Sasha Rakhlin and Karthik Sridharan. Optimization, learning, and games with predictable sequences. Advances in Neural Information Processing Systems, 26, 2013.
[50] Herbert Robbins and Sutton Monro. A Stochastic Approximation Method. The Annals of Mathematical Statistics, 22(3):400 - 407, 1951. doi: 1//11///mm//177729955 URL https:
//doi. 1214 /tt////os///////////// /.mmmm17772958.
[51] Gareth O Roberts and Richard L Tweedie. Exponential convergence of langevin distributions and their discrete approximations. Bernoulli, pages 341-363, 1996.
[52] Ralph Tyrell Rockafellar. Convex analysis. Princeton university press, 2015.
[53] Ruoqi Shen and Yin Tat Lee. The randomized midpoint method for log-concave sampling. Advances Neural Information Processing Systems, 32, 2019.
[54] Yee Whye Teh, Alexandre H Thiery, and Sebastian Vollmer. Consistency and fluctuations for stochastic gradient langevin dynamics. Journal Machine Learning Research, 17, 2016.
[55] Belinda Tzen, Anant Raj, Maxim Raginsky, and Francis Bach. Variational principles for mirror descent and mirror langevin dynamics, 2023.
[56] Santosh Vempala and Andre Wibisono. Rapid convergence of the unadjusted langevin algorithm: Isoperimetry suffices. Advances in neural information processing systems, 32, 2019.
[57] Max Welling and Yee Whye Teh. Bayesian Learning via Stochastic Gradient Langevin Dynam- ics. In ICML, page 8, 2011.
[58] Andre Wibisono. Sampling as optimization in the space of measures: The langevin dynamics as a composite optimization problem. In Conference on Learning Theory, pages 2093-3027. PMLR, 2018
13

---

[59] Andre Wibisono. Proximal langevin algorithm: Rapid convergence under isoperimetry. arXiv preprint 11Xiv:1911.01469, 2019.
[60] Pan Xu, Jinghui Chen, Difan Zou, and Quanquan Gu. Global convergence of langevin dynamics based algorithms for nonconvex optimization. Advances in Neural Information Processing Systems, 31, 2018.
[61] Kelvin Shuangjian Zhang, Gabriel Peyré, Jalal Fadili, and Marcelo Pereyra. Wasserstein control of mirror langevin monte carlo. In Conference Learning Theory, pages 3814-3841 PMLR, 2020.
[62] Kelvin Shuangjian Zhang, Gabriel Peyré, Jalal Fadili, and Marcelo Pereyra. Wasserstein control of mirror langevin monte carlo. In Jacob Abernethy and Shivani Agarwal, editors, Proceedings of Thirty Third Conference on Learning Theory, volume 125 of Proceedings of Machine Learning Research, pages 3814-3841. PMLR, 09-12 Jul 2020. URL https: / //proceedings.. mlr peoess/v22///e http/////////////// /hagg20a..22.... html.
[63] Difan Zou, Pan Xu, and Quanquan Gu. Sampling from non-log-concave distributions via variance-reduced gradient langevin dynamics. In The 22nd International Conference on Artificial Intelligence and Statistics, pages 2936-2945. PMLR, 2019.
14

---

A Further Examples of LRM Schemes
Example 5. The classic Stochastic Gradient Langevin Dynamics [57] iterates as
Xk+1 = X Yk+1 Vxxx) + /2yR.1 Ek+1.
(SGLD)
where f is the gradient of the negative log-likelihood of random batch of the data. (SGLD) fits the LRM template by setting Uk+1 = fxxx)) - V V and bx+1 == O.
Example The Proximal Langevin Algoritmm [9, 46, 59] is defined via
Xe+1 X Y V (x (2 EL+1.
(PLA)
This algorithm is implicit, and it is assumed that one can solve (PLA) for XK+1- By setting bx+1 == ff(++)) fxxx) and Uk+1 == O, that this algorithm also follows the LRM template.
B
Additional Related Work
Our paper studies the behavior of a wide range of Langevin-based sampling algorithms proposed in the literature in the asymptotic setting under minimal assumptions. This allows us to give last-iterate guarantees in Wasserstein distance. As stressed Section I, our goal is not to provide non-asymptotic rates in this general setting as the problem is inherently NP-Hard. However, given more assumptions and structures on the potential f. there is plethora of works which prove convergence rates for the last iterates in Wasserstein distance. In this appendix, provide additional backgraound for these works and the methods used in the literature.
A powerful framework for quantifying the global discretization error of a numerical algorithm is the mean- meean-quaree analysis framework [40]. This framework furnishes a general recipe for controlling short and long-term integration errors. For sampling, this framework has been applied to prove convergence rates for Langevin Monte-Carlo (the Buler-Maruyano discretization of (LD)) in the strongly-conven, setting [32, 34]. Similar to our work, the convergence obtained in these works is last-iterate and in Wasserstein distance. One of the essential ingredients the latter work the contraction property of the SDE, which is ensured by the strong convexity assumption. This, in turn, implies strong non non-asymptotic convergence guarantees.
It is an interesting future direction to study the combination of the Mean-Squared analysis together with the Picard process and its applicability more sophisticated algorithms (such LRM schemes with bias and noise), as well as non-convex potentials.
As explained in Section 3, one of the main themes in proving enror bounds for sampling is the natural relation between sampling and optimization in the Wasserstein space. This point of view, when applied strongly-convex potentials, has produced numerous non-asymptotic guarantees; see [14, 18] for a recent account and the references therein. Note that strong convexity is crucial for the analysis used in the aforementioned work. Moreover, the error bounds for biased and noisy discretizations do not decrease with the step-size or iteration count; see [18, Theorem 4, Eqn. (14)). This means that while the bound is non-asymptotic, it does automatically result an asymptotic convergence. Finally, we stress that these approaches are orthogonal to our techniques: We view a sampling algorithm as a (noisy and biased) discretization of a dynamical system (and not necessarily gradient flow), and use tools from dynamical system theory to provide asymptotic convergence results.
C Proofs for Section 4
C.1 Proof of Theorem
In this appendix, bring the detailed proof of Theorem 1. Recall that we interpolate the iterates of the LRM scheme {xx} as
E[Z [1} + (xxx (B, B.z). X, = Xx + ( (1(x2)
(3)
Moreover, for a fixed > 0, we considered the Brownian motion B," = Br+s B, and constructed two important processes: the Langevin flow defined via
d v( v( ds (x) dB
4 i1 X,,
(12)
15

---

and the Picard process (6) constructed as
y X, + S ul v(X+tu) du + S-" o(Xetu) dB
(6)
Let fix T > O, and for Langevin flow as
[0, 7] decompose the distance between the interpolation and the
5 1 - - 11 x ||r!" -6!"" 11x,s -Y 11
(7)
where we have used ||a + b11 21|a|| + 2||b|l We now bound each term of this decomposition. Notice that due to the synchronous coupling of the processes, the Brownian motion cancels out in the differences.
The first term controls how close the Picard process is to the Langevin flow, and is bounded in the following lemma.
Lemma 3. For fixed 1,T > 0 and 0 s s s T, the distance of the Picard process and the Langevin flow is bounded as
1 -0, 1 2(T 1)22 [ 110 Xxtull du.
Proof. By the auxiliary Lemma below, Lipschitzness of e, or Itô isometry (see, e.g., [62]) and s - T, we have
11 =E [/6 (() o( ((x+t) du f( o( o o(( dB!!
EIlY!
! Elucco!)) v u(Xitu)) || du 2E I ( du $ 2(T 1)L2 [ Ell0" ( - Xy+ull du.
2s
For the rest of the proof, we need to define the continuous- time piecewise-constant processes Zk+1. and Z(Tk + s) [E[ Frats, for OS < Yk+1- Also, let m1) = sup 2 0: Tk $ r) so that Tm(r) S Tm(r)+1-
X(Tr s) XR. 7(Th + s) Yk+1, Z(TA
To bound the second term in (7), we have seen that
1+s
v(X(u)) du S( v(Xreu) du
X14s ((1) =
i+s
o(X(u)) dB, - S." ((X(+u) dB(()
+ st,
where (z) plays the role of accumulated noise and bias from time to s, and is defined as
Az(1,s) = Zi+1 (1+ T2 [() Frtt] (t T,) E( F.).
(13)
with & = m(1+ s) and n= = m(). We therefore have
X - 11 s
o(Xu) - u(X())) du
( 148
o(Xx) ((((e) dB, 12
3E
(1.
El(((x)) - ((x(u)) du
$ 3s
(e(X)))((X((x))) o(X(u))] du 3E,]) 33(s 13($+1)12 1)1.3 ( EIJX, x21)112 3E,IAZ(1,s)1' $
+3E
(14)
16

---

For bounding the term inside the integral, we have
E| - x()lll = Ell( Tm(la) (o(X(u)) Z(u)) +(((()) (Bu, - B BTm )11
4y( 4
2y(w) Tu(o((())TTxxx
We have used the fact that
El())(X((u) (Bu B Brmii))ll'
E( )o((() ((((u)) (Bu, B Brimis)
= =Etr(o(X((u) (X(()) (Bu - B Brmia) (Bu B Brimin)) )7)
= =E(E((((((((u)) To (X(1)) (Bu - B B (Bu - B B )1) & Fiminl)
= -u- (T- ((()())))))))))))))))
Notice that since conditional expectation is a projection in L.2, we have ElIZ(u)]1 ElIZ(011 Using this fact, along with boundedness of (+)) by Co. and Lemma we get
E X (1 47 (T(x)(E
27(4) Ex((((())[([X]XXXX))) (x)) (x]
2CC,((u) C7(u),
for some constant C > 0. Plugging this estimate into (14) after taking expectation yields
E[IIX++> - y, 11 s 3(s+ 1)L2c J,
7(u) du A1,
$ 3(s+ + l)sL'C
sup
7(u) + 3E (z.)Z((
we[r.rts]
$ 3(T + 1122.22
sup we[r,s+1]
7(u) + sup (I,)
ue[0.T]
Taking supremum over s E [0,7] and noticing that the right-hand-side is independent of s and Yx O, together with Lemma yields
A, = sup
1| 110
(15)
055ST
3(' 1
sup 7(u) +3 sup
1SUSH4
OsusT
co,
showing that the Picard process gets arbitrary close to the original interpolation, as - co. Let us return to the decomposition (7). By taking expectation and using (8) and (15) we obtain
|| - 2 |2]
2(T 1)L2
E|1IX 1/ du 2A,
2A, exp s7 1)1.2)
$ 2A, exp((T 1)²
where in the last line we have used the Grönwall lemma. Thus,
lim 1-000
sup sE(O.T]
|| ! |2 =
Recall that the Wasserstein distance between X,+x and ( is the infimum over all possible couplings between them, having the correct marginals. As ( has the same marginal as the Langevin diffusion started from X, time s, and the synchronous coupling of the interpolation and the Langevin flow produces specific coupling between them, we directly get
W2( (X14 Φ (1) SE || - ① 11
which implies
lim sup W2 ((())s).(!)))) =0. 1-00 sc(0,7]
as desired.
17

---

C.2
Auxiliary Lemmas
Lemma 1. Suppose Assumptions 1-3 hold. Then, for any fixed we have
lim sup 053ET
(t)
Proof. Define As and Au the same way as (13). By Cauchy-Schward we have
(1,
2
R-1
Yiti |bb+il| (I++s Tx)]|E[b<+]
Fitlll + (1 [T)]][[bnn)) [9,]ll
n2yn=1 +
Σ Weillbull ( Fit +1 Tn) T,) 1,
where the last inequality comes from Σk-1 Y/41 S 5, TK S Yk+l, - Tn S Yntl, and yk+1 $ Yntl-
Noticing that conditional expectation is contraction 2 and letting k' m(1 7), we get
sup (1, s) s (2+T) 0sssT
Y:t F+ b Σ
7j+1Ellb 112 Ynti EE| b 112
sup nsjsk'+]
Now, invoking Lemma yields
K--1
sup ElIl = 0sxsT
Σ Vi41
sup
2/t1 Vnti
nsisk441
K-11
s C(2+7)
Σ
Vi+1
2y
s
sup
7(1 s).
0s:5T
As f - the last quantity vanishes, since yn - 0.
For the noise we have
k-1
| (t, ()1 2
Σ
Vi+l U++1
+ +4|1 - T(UL( FF+,0111 +4][(1 E,U, [F]1
1=
s
Σ
Vi+1
U++1
+ 14 1/
|=
Taking expectations and then sup, we get
sup (t. 2 US
sup n+1 sk skm(teT)
E Σ
Yin Ui+1
+4y 44 0² + 4y2 02
Since (U:) is martingale difference sequence, we have that {xR-- Pios Uit1 Fesn is martingale. Thus, by the boundedness of the second moments of U;, we get
k-1
Yi4
U1+1
Σ Ion
24
E+|
2 s
Σ 1m
yiu-
Hence,
lim sup Ell
>
Vist 12 in < k m(T, +7)
s lim 0²
Σ
Yi+1 = 0.
i=n
n-000
i=
18

---

Lemma 2. Let (x& KeN be the iterates of (LRM) and suppose Assumptions 1-3 hold. Then, E||xx|1²= 0(1/yk=)). This in turn implies El| (11/ and =
Proof. Without loss of generality, suppose has stationary point O. We repeatedly use the fact that [l[ x L Moreover, by Assumption we have (v(x),x) Co(Ilx] + 1), and
(x) Co.
Define ak E| We have
ax+1 ak v Ello(xx)) ZK+1 (1: 2y v(xx) + Zk+1)
+ k+1 E(xx, o (xx) 27 + ZR-1.
$ 22 ak + 22 21 yk+1Co 2yk+1 Cyk= 1) + 2yk+1 vax v x
2 k+1 CC v /
(16)
By Assumption there is some C > 0 such that El|b++1|1 Cb(yk) at Yk-1), and we have
E| s 2 1| 21E| s 2Co(( + Yk+1) + 20
(17)
Moreover, vp + 9 vp- vg. have
1I s (yC4(vk vyk+1) + V20.
(18)
Plugging the bounds from (17) and (18) into (16) gives
ax+1 dk s 21. vk v ax + 4CDyk- ax + 4Cbbkk1 47 k-1 x² Co+ 2yk+1 Co Van 2y
3/2 val + 22220 Vak 2c a + 22
(19)
+2 512 k+1 Vax C 2y v2
Py ax @Yk+1 Vax RyR+1,
where
P P= + +4cby2 + +2/2Cb
Q = 2C, + 2C, Vy 2V2or 2V + 2 3/2
= 4Cby2.1 + 4 Co + 2C, + 2 2c + 2y 2C
1/2
The exact values of P, Q. and R are irrelevant, and we only need upper bounds for them. Assuming that Yk+1 < for all k, we replace the three quantities by
P = 2L + 4C, +2/2C, Q = 2C, 2v + 2V20r 22/CC + +2/2C6C.
(20)
R = 4C, 2 C. + 2C, 2V2C,C + +2V2Cu00
Now, define h2 = v2 1 ak. The recursion (19) in terms of h4 becomes
F222 + vhelr + Ryk-1 hx s (1 + Py
TE+1
We now prove that there exists some M O so that h2 MyL+1 by induction. Suppose it is the case for k, and we prove for + 1. Using the induction hypothesis we get
T22 ha41 & Myx+1 Py
MMY.. RyK-112.2 R2-
22
=M(1 Py2
+ +VM @ + RyL+1Vi42
Yk+1
19

---

For the last to be less than yx+2 we have to verify
M(1 + Pyk-1)
Yk+2
+ VMQ + Ryk+1yk-2 $ M
Yk+1
equivalently,
M Yk+2 Pyk+11142211 VMQVILHIYK=> + RYk+17142 $ O.
yk+1
This a quadratic equation VM, and for this inequality hold, we prove that the leading coefficient is negative, and the largest root is bounded above by some constant not depending on
Negativity of the leading coefficient is equivalent
Yk+2 Pyk+1 < 1,
Yk+1
which is implied by our assumption on the step size.
The larger root of the equation
(-4y k42 PR Fk+1 (Y2:2 + 4RY 1/2 + +
2(1 Yk+l ( Yko2/
V Q + RY
Tk+1 7k42//k+1)
VKk+17k+1 VRykk
(1 Yk+1YR=2P 7k+2/7k.) 7k+1)
By our assumption on the step size that
7k+2
PY < yk+1,
Yk+1
we get that the larger root is smaller than
V VRyk+
V VR < C+
Yk+1
Letting M == Q + VR gives the desired result.
The second argument of the lemma follows from Assumption 3 and the first result of the lemma. Lemma 4, For a vector valued function g (RR; Rd), one has 1..0 g(u) du 1
(f, [lo(w))l du) S s l1l(w)ll2 du.
D Proofs for Section
D.1 Proof of Theorem 3
For brevity, let us write Fx instead of F Opening up || ||x (x=x()(xx)) Zx+13 (x4))r(x))) Ex+ and ignoring every term that is zero-mean under E(. Fx), we get
E[]|xen|| Fx] E| [|xx|]] + 2yk+1 ((x,v(()) Zent)
v(+1llu(xx) Z++1ll (x=lllo (o baa)) Fi
||x 2 + 2yk+ ((x,, C.r (2
E
Fs
|| 22 + 2yk= (NA.. (VN.D(N2)) ((x) C.,
k+1 C., / + + 11xx))1
(21)
+ +E(2> [1z 11
+ +
2E(yx (xx, bxr])
20

---

Recalling (5) in Assumption 3, we have for some C > 0
E|| 20 + (y1 +x-11)
(22)
Without loss of generality, assume yk z and || 2 I (so that (E]]xx&|13) 2 EllxxI19) for all k. Then, [ I L? |2||x I1 together with Assumption 4 and the Cauchy- Bouchy-Schwartz inequality on the last term of (21), implies
1 E||xa+1 21? s E||xx|| (arx++Ellxxx 2yk+1 B Co + 2/ 2 Co
+ /
+ 2
20 2C L Ellxx|1 + xK+1))
+ v C
(L Ellx|| + yka1)
El|xx11
+ + VC
1.2 (E]1)x2133* y
s (1 -1 C2Y C37k+
for some constants C1, C2. C3 depending on L, L.C, B, and d. Since Yk there exist â. and kn such that, for all k 2 ko.
E|| 11 s E||xx|I ā + BYR+I, 1- Ty++ > 0.
simple induction yields
B T Elj.xap
sup [F]]xx11 s max
which concludes the proof.
D.2 Proof of Theorem 4 for Constant Diffusion
Before proceeding, we need a lemma which can be distilled from [20, Proposition 8]:
Lemma 5. Suppose f is L-Lipschitz. Fix E Rd and y O, let = x yVfxx) 22y.. Then
exp 1 (V(W)) xx) [11 xlp) s (11-
(23)
Let Re+1 Yk+1 Vx(x)) /2yxr1 Ex+1 so that x2+1 -R+l Yk+i (Ust bat1).
Conditioned on xR. UKII E and using the L-LLpschitzness of VS, we get
Ext (x211)
$ Eexp
2 (Vf(xx), - xx) + 1 |1 xxll3)
(24)
s Eexp
1 (Vx(xx)) Ex+1 -xx) / fxxx,, kkm
(25)
コ (V
L + " [1 -xx12 + LrinnllUee + Lyhallbbe+11
(26)
Let e (0, 1). Since
(V (V((x)) s x yε., 11..
2 Yribe4)) $ y [VV [| +1|bx+11?
21

---

we have
els (se)))
(27)
$ Eexp (V)(c))) - xx) z ||a - xull
(28)
+y2in + + + (LRET v |1 ( 1)1beanlP
(29)
Invoking (11) denoting c' (LY + 1) we get
é (f Eel S sea Eexp ] (V - -xx) + 1lik+t - xxxll c
,
(30)
where,
A₁ 2 (vrri
2-5 r2es + +cyrk.)) (
(Lrkn / ||Uk-11?
(31)
+ + |U++1 +
Recalling that (2 = Eatl + Yk+1 (x(x),, we have FR s || xxll2 2 VV f(xx) and thus
Eexp [ (V)( + (2 (xal suarr. -
e Eet se Eexp
(32)
where 11 A2 c'r fxx)) Lemma 5 then implies
e - (x)) (xeal) se (1-))+1L')
(33)
where = Yerk [V{(x)|
We now take the expectation over UK+1, Uku E41 (in other words, we are now only conditioning on xx). Set E (1 - O. Since Uk+1, Ee1 are sub-Gaussian and since yk - 0, for k sufficiently large we have
EAK s
exp (- Yk+1
+ Ykal c'ybn + +ccrtn)))IVVCC
(34)
4
$ (1+ -- Thrl (x))
(35)
To summarize, we have shown that, conditioned on Xg.
ts (xx) (xx) Eetf (xea) < (1 - L') "hd |VV (xi)l²
(36)
simple induction la [20, Lemma & Proposition 8] then concludes the proof.
D.3 Proof of Theorem for Mirror Langevin
Here, we bring the proof of Theorem 4 for the case of Example 4 and without noise. The proof for the noisy case the same as in Appendix D.2.
Define
x V (V
where a standard Gaussian random variable. Let U(x) = f(Vo' (x)). For a fixed x, we have
Ee U
(2m)d/12 / exp U
1 -()
11E11
dg
2
22

---

Notice that we have
(v2 (xx)) x VV
v2y
which implies
dg = (V2y) /det (x) dx*
Thus, the integral, after the change of variable from & to to.x* becomes
1 ?/ ʃ exp -2x 4y
1
o dx' (37)
with C = (4xy) det (20' Now we use the smoothness of f:
= -
(VV' -xx) + f1xx - x112
On the other hand, we have
1Iv 60 - VV (x))12
= IVV - x)1 (
2y 00 (V) (x)xx*
Notice that in (37), the colored terms cancel out, and what we are left with is
Ee -U0)0003000
1 $ C f exp 11xx x11 4y
(
(V - x)1 (1
As, by our assumption, 22 is bounded from above and below, we get the exact form as in Lemma 5. The rest of the proof the same as Appendix D.2.
D.4 Proof of Proposition
In this section, we prove that Examples 1-6 satisfy our bias conditions, which, as we have seen in Section 5, implies Proposition For brevity, we write F2 for FTL.
$ Proof for Example For randomized mid-point method, by replacing (((x)) and A (x2+112 with fx(x) U/ and V f(x=1122 + UK+1 respectively, we have
xk-1/2 = x 7 (x)(xx) Uk+r1 } + 22
XK+1 X Yk+1 (((()))))))) + U2+13 22
where }ax} are i.i.d. and uniformly distributed in [0, 1], {Ux} and {0,} are noises in evaluating f the corresponding points, and 5% are independent standard Gaussians.
Notice that the Lipschitzness of V,, and the fact that a& S implies that the bias term bet1 == V (xx+1122 (xxx) satisfies
E[ |12 F < L'El|INuH1/22x11111 Fi
$ 1.2 (22, [[IIV]((x)) U/,1 1 1x]++ 2Yk+1d)
$ 21.2 2 21. dy
+
23

---

$ Proof for Example 2. Recall that the new algorithm Optimistic Randomized Mid-Point Method has the iterates
Xk+1/2 = X YK+104+1 (xxxx)))+ + 27k-104+1E
Ak+1 =x Yk+1 (xx+11/2 22 Ek+1,
where (axt, E. and ff are the same as in (RMM), and the noise and bias are Uk+1 fx(xx1122) - V (x(+1122 and bx+1 = (f((x=1122= f(xx). We have
E[[|b++1] 112 [. = -[I==E((((<>>>>>>> - - (xxx)] [x
1.-E[[1=R-112-20 xx 1 [
LL'Ell- 7k++ V fx 22
- E f Fx1 + +2L YK- o² + dyk+1.
Similar to the proof for Example 6, notice that | (x
2 / vs(x))1
2 As yk 0, can assume that 2 k < 1. and we get
E|||b++1 | = IV
0² 8L.2 dy 1V rx.)),
as desired.
$ Proof for Example 3. The iterates of stochastic Runge- Kutta Langevin algorithm is follows:
h1 = Xk 2yk++1 (1/2+ + 1/V6) Ek+1 + E /V1
h/2 = Xt: Yk-1 (V=
2y (1/2 - 1/V6) Ek+1 + E /V12
Yk+1
(Vf(h) U + 22 EK+1,
X1+1 = X₁
2
where Ex+1 and 5k+1 are independent standard Gaussian random variables independent of xx, and Uk+1 and U/1 are noise in the evaluation of f-
Observe that
2 1u+1 1 (Vf(h) - (V/hh)
We have
2 1/12) =0
and
- 2 2 1 22d(1/4 1/6+ 1/ 1/12)yk*1)
= = +
We thus have
E[Ilb*+1 111 } - (
= (V 7k+1),
as desired.
$ Proof for Example Suppose a Legendre function [52] for Rdd and consider the iterates
Xk+1 Xx - Yk+1 VVVo' (x+)'(xx)) + 22yx., (xx) -1)1/2 EL+1.
where 6² is the Fenchel dual of d. that is, 6" (x) sup suppoRd (x. (()). Also recall that [52]
V xx, (v
Vx E Rd
Let u= -f V and 1120")1122 First, mention what assumptions imply on f:

---

The Lipschitzness of corresponds to similar condition in [31, A2]: [V - V ff(il s 4IIV((x) vb(y))l
The Lipschitzness of or in Frobenius norm corresponds to modified self-concordance in [31, A1]:
[V' (y) !/2 || s LIV vb(y))l.
Boundedness of or in Hilbert-Schmidt norm implies
v (xx()) s Cor
Dissipativity and weak-dissipativity of corresponds the conditions below, respectively:
(v)(x), vf(x)) 2 al]Vo(x)]] B. (Vé(x), vf((x) 2 a]lVo(x)]l B.
If and & satisfy the conditions above, then the mirror Langevin algorithm Example 4 fits into the (LRM) scheme.
Remark. Note that this version of Mirror Langevin cannot handle the case where e-x supported on compact domain; in that case, the Hessian of & has to blow up the boundary, and will not satisfy our boundedness assumption. The version of mirror Langevin we consider in this paper, though, be thought setting has also been studied in the literature, see [55].
adaptive conditioning method for densities supported Rd This
$ Proof for Example 6. The iterates of (PLA) follow
X1+1 = X4 Yk+1 V fx(x+1) + v 22 Ek+1-
(PLA)
We mentioned that the bias term is bk+1 = ( (xx)). Now remains to prove that satisfies the conditions (5) and (11). We have
E[llb=ill| Fx]
Vf(xx)11 Fx
< L |
= -. V + 2 Ek+I [1A]]
$ 21.3y2 1 (x]] ||xx+ +
Now, notice that [V(++)) z 2|[ V f21 As 0, one can
assume that 2 3, and we get
E[[|b++1|] 2| | $
Elllbx+1 1 1| [VV[<x)1 41.2
which implies
E[[bb=-1|||3]] s [Vf(x))1 +8L + =1 yk+1),
as desired.
25


Reviews:
[{"review_id": "4r7e0tMqIP", "reviewer": "Reviewer_hN2a", "summary": "This paper gives a unified asymptotic analysis of a broad class of stochastic algorithms that encompasses several variants of the Langevin algorithm. In particular, it can handle issues of inexact gradients, bias, noise, and problems beyond gradient-based algorithms. The key technique is the introduction of an intermediate process, termed the *Picard process*, which fits between the iterates and the continuous-time ", "questions": "* Can you give examples beyond Langevin-type algorithms for which you can apply the  method?\n* Can you get some more quantitative guarantees beyond asymptotic convergence?", "limitations": "These are adequately addressed.", "rating": 7, "confidence": 4, "soundness": 4, "presentation": 3, "contribution": 3, "strengths": "The main strength of this paper is the unified nature of the convergence results. The general framework encompasses a fairly large variety of Langevin-type algorithms, and likely a decent class of problems beyond Langevin algorithms. It also gives a clean, unified analysis. ", "weaknesses": "My only major criticism is that the paper seems to overstate the novelty of the analytic methods. In particular, I have not seen the specific Picard process defined here, but several works introduce related processes that fit between the iterates and the continuous time process. This enables similar triangle-inequality based convergence proofs. For example (but not limited to):\n\nChau, Ngoc Huy, et al. \"On stochastic gradient langevin dynamics with dependent data streams: The fully nonconvex case.\" SIAM Journal on Mathematics of Data Science 3.3 (2021): 959-986.\n\nBubeck, Sébastien, Ronen Eldan, and Joseph Lehec. \"Sampling from a log-concave distribution with projected Langevin Monte Carlo.\" Discrete & Computational Geometry 59 (2018): 757-783.\n\nAdditionally, last iterate convergence guarantees are not particularly rare. Both works cited above give last-iterate bounds from corresponding stationary distributions. Many of the works citing these papers do as well. \n\nOn a minor note, there are some confusing notations.\n* $b_k$ is used for the bias, but then gets re-defined in the proof of Lemma 2.\n* Using $\\sigma$ for both the diffusion matrix and the variance bound on the gradient noise is mildly confusing.\n* At the end of the proof of Lemma 1, it should be the limit as $n\\to\\infty$, instead of $k\\to \\infty$", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We are sincerely grateful for pointing out the missing references and remarks. We reply to your questions below, and we will revise our manuscript accordingly in the upcoming revision.\n\n> My only major criticism is that the paper seems to overstate the novelty of the analytic methods. In particular, I have not seen the specific Picard process defined here, but several works introduce related processes that fit between the iterates and the continuous time process. This enables similar triangle-inequality based convergence proofs. For example (but not limited to): Chau, Ngoc Huy, et al. \"On stochastic gradient langevin dynamics with dependent data streams: The fully nonconvex case.\" SIAM Journal on Mathematics of Data Science 3.3 (2021): 959-986. Bubeck, Sébastien, Ronen Eldan, and Joseph Lehec. \"Sampling from a log-concave distribution with projected Langevin Monte Carlo.\" Discrete & Computational Geometry 59 (2018): 757-783.\n\nWe agree with Reviewer hN2a's assessment and sincerely appreciate the references they have provided. We will duly update the paper and highlight our contributions, as well as acknowledge that similar ideas have been explored in prior works.\n\nWhat distinguishes our work from the existing literature is the advantage of generalizing the Picard process to encompass a vastly wider class of algorithms, specifically the Langevin-Robbins-Monro schemes. Moreover, the integration of the Picard process with the theory of asymptotic pseudo-trajectories plays a pivotal role in our analysis, and both of these aspects present original contributions.\n\n> Additionally, last iterate convergence guarantees are not particularly rare. Both works cited above give last-iterate bounds from corresponding stationary distributions. Many of the works citing these papers do as well.\n\nWe agree that there are numerous results on last iterates for different settings and algorithms. What we intended to express is that little is known in the **generic non-log-concave setup** (without convexity or functional inequalities) and scenarios involving **biased discretization**. For example, Bubeck et al. (2018) is for log-concave target distributions, and  Chau et al. (2021) is for unbiased gradient estimates (see Eqn. (7) in their paper).\n\nTo alleviate any misunderstandings, we are committed to revising our exposition on related work.\n\n> On a minor note, there are some confusing notations. b_k is used for the bias, but then gets re-defined in the proof of Lemma 2. Using sigma for both the diffusion matrix and the variance bound on the gradient noise is mildly confusing. At the end of the proof of Lemma 1, it should be the limit as n -> infty , instead of k -> infty.\n\nThanks a lot for pointing this out. The end of the proof of Lemma 1 should be $n \\to \\infty$ as you mentioned. We will make the notation more succinct in the final version.\n\n> Can you give examples beyond Langevin-type algorithms for which you can apply the method?\n\nEssentially, our framework can be applied across a spectrum of continuous-time dynamics, such as Underdamped/Higher-order Langevin, Hamiltonian Monte Carlo, Neural SDEs, and more. In addition, the scope extends to diverse discretization methods for these dynamics. Examples of such include Euler-Maruyama, leap-frog, and symplectic integrators.\n\n\n\n> Can you get some more quantitative guarantees beyond asymptotic convergence?\n\nAt the expense of stronger assumptions, yes: We are able to obtain results that are close in spirit to non-asymptotic results. Assuming step-sizes $\\gamma_n = \\Theta(n^{-1})$ and exponential convergence of continuous-time dynamics (which holds under, e.g., LSI), we can show that an LRM scheme converges in $\\mathcal{W}^2_2$ at a rate of $O(n^{-1})$ under a warm-start condition, **under the same noise and bias conditions**. While we cannot estimate the duration of the burn-in phase, our argument holds for a wider range of stochastic and biased algorithms than in the current literature. However, as our primary focus lies in the generic setting where such assumptions are unavailable, we have chosen to defer these studies to future work.\n\n\n---\n\nWe hope that the above addresses your questions - but please let us know if any of the above is not sufficiently clear.\n\nThank you again for your input and positive evaluation,\n\nThe authors"}, {"review_id": "yJVqvQTxfv", "reviewer": "Reviewer_vdDX", "summary": "The work studies when a discretized Langevin dynamics under the Robbins-Monro-type stepsizes can converge to the Gibbs distribution. The paper obtains asymptotic results with very mild assumptions, and the framework not only includes Euler discretization, but many other sampling schemes as well, such as mirror Langevin, proximal, randomized mid-point and Runge-Kutta methods. The analysis builds upon constructing a continuous-time trajectory via interpolating the iterates, Wasserstein asymptotic pseudotrajectory and checking the stability condition by invoking the dynamical system theory.", "questions": "(1)\tPage 4. ``The usual Lyapunov-type analysis for sampling algorithms focuses on bounding the change in relative entropy across iterations…” ``this makes the Lyapunov analysis applicable only to the simple Euler-Maruyama discreteization of (LD)’’ I am not too sure whether these two statements are accurate. Lyapunov functions are often used in analyis of Langevin algorithms, to show uniform bounds on the moments, e.g. Raginsky et al. (2017), in the coupling methods, e.g. Dalalyan and Riou-Durand (2020). It is definitely applicable beyond the Euler-Maruyama scheme, e.g. Dalalyan and Riou-Durand (2020) uses the disretization proposed in Cheng et al. (2018) to analyze kinetic Langevin dynamics.\n\n(2)\tOne technical point I would like to see more discussions is that in equation (1) in your Definition 1, it is for fixed $T>0$. Actually the dependence on $T$ can be exponential in $T$, which is quite common for weak approximation error in the literature. However, in order for the Langevin algorithm to converge to the Gibbs distribution, one often needs uniform-in-time guarantees, and would you need $T\\rightarrow\\infty$ in order to obtain Theorem 2?\n\n(3)\tTheorem 2 is a very nice and clean result. But I am surprised that you only need assumption (10) which is an assumption on the discretized dynamics only. The reason I am asking is that it seems to me that Assumptions 1-3 alone do not guarantee that the continuous-time Langevin SDE has a unique stationary distribution. If Theorem 2 holds, that means assumption (10) can imply that the continuous-time Langevin SDE has a unique stationary distribution? The existence of $\\pi$ is necessary for Theorem 2 to hold.", "limitations": "I did not see such discussions about limitations.", "rating": 5, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "The paper is built upon solid mathematical analysis and it is a nice contribution to the vast literature of Langevin algorithms in machine learning. The builds provides a unified framework for the asymptotic guarantees under the Robbins-Monro scheme.", "weaknesses": "Even though the mathematical theory is nice, the results may not have too much practical importance. The author(s) emphasized that this work is about asymptotic analysis instead of the non-asymptotic analysis. However, because the results are of asymptotic nature, it not clear to me what insights the results can provide concerning the schemes such as mirror Langevin, proximal, randomized mid-point and Runge-Kutta methods, because you only have asymptotic guarantees, it is impossible to use these results to compare these algorithms with more classic and basic Euler discretization of Langevin algorithms. \n\nAlso, I find the discussion of existing literature less satisfactory. Some of the claims and statements may not be that accurate. For example, on page 1, ``Existing guarantees suffer from the drawback of lacking guarantees for the last-iterates’’  ``the convergence is typically given on the averaged iterates instead of the more natural last iterates’’. To the best of my knowledge, there are numerous works on Langevin algorithms in the past decade and most of them are about last iterates guarantees in Wasserstein, KL or other distances; see e.g. Dalalyan and Karagulyan (2019), Dalalyan and Riou-Durand (2020), Ma et al. (2019), Ma et al. (2021), Raginsky et al. (2017), Gao et al. (2022). ``and little is known beyond the elementary schemes of stochastic gradient Langevin dynamics.’’ This is not accurate either. There have been many studies in the literature about underdamped Langevin, high-order Langevin, non-reversible Langevin and other variants of SGLD such as Dalalyan and Riou-Durand (2020), Hu et al. (2020), Ma et al. (2021), Mou et al. (2021), Gao et al. (2022).", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We sincerely appreciate Reviewer vdDX's thoughtful criticisms. After a thorough reading, we believe that these critiques primarily stem from presentational issues, which we fully acknowledge exist and will commit to improving them following your suggestions.\n\nIn light of this, we sincerely ask for a re-evaluation of our work based on the points addressed in our rebuttal. We are open to further discussions and eager to address any additional questions or concerns they may have.\n\n> Because the results are of asymptotic nature, it not clear to me what insights the results can provide such as mirror Langevin, proximal, randomized mid-point and Runge-Kutta methods.\n\nIn our revised version, we intend to highlight the following two insights that our results offer to practitioners:\n\n1. **Validating existing methods:** We first note that methods like mirror Langevin and randomized mid-point currently lack even asymptotic guarantees in fully non-convex scenarios, such as sampling from neural network-defined distributions. Our work fills this gap by offering the first solid justification for these schemes, supporting practitioners in utilizing these methods confidently.\n\n2. **Facilitating new algorithm design:** Our work motivates novel sampling methods through a straightforward verification of **Assumptions 3**. An illustrative instance involves the randomized mid-point method and Runge-Kutta integrators, wherein a substantial 50% reduction in computation per iteration can be achieved without compromising convergence by simply recycling past gradients; see **Example 2**.\n\n    We do acknowledge the inherent limitations of our approach: The balance between the benefits of saving gradient oracles and potential drawbacks remains an open question, necessitating case-by-case practical evaluation. Nevertheless, our theory provides a flexible algorithmic design template that extends beyond the current literature's scope.\n\n> There are numerous works on Langevin algorithms in the past decade and most of them are about last iterates guarantees...\n\nWe apologize for the confusion and would like to clarify our intended meaning concerning the absence of guarantees for **fully non-log-concave setups**. This refers to scenarios lacking assumptions of convexity or functional inequalities like LSIs or Poincaré. The references cited by the Reviewer fall within the latter classification: Dalalyan and Karagulyan (2019), Dalalyan and Riou-Durand (2020), and Ma et al. (2021) are centered around the (strongly-)log-concave context. Meanwhile, assumptions 1-3 in Raginsky et al. (2017) and assumption 1 in Gao et al. (2022) imply Poincaré inequalities.\n\nWe express our gratitude to the Reviewer for bringing this point of confusion to our attention; we will make the necessary modifications accordingly.\n\n> On ``little is known beyond the elementary schemes''\n\nWe would like to clarify that our reference to \"elementary schemes\" pertains specifically to discretizing the **Langevin diffusion** mentioned in equation (LD). While Reviewer vdDX has rightfully pointed out the existence of elementary discretization schemes for **other SDEs**, such as the underdamped or higher-order Langevin, our primary focus in this paper remains centered on the Langevin diffusion. [It is worth noting that our framework can be adapted to analyze these schemes.]\n\nAs for the existing SGLD variants, which encompass variance reduction methods, it is important to highlight that their core emphasis continues to revolve around unbiased gradients.\n\n### Q1\nWe thank the Reviewer for raising this question, which again stems from our presentational issue: Our primary focus is the impact of **bias** in these analyses. As the Reviewer rightfully pointed out, handling noisy gradients' impact on entropy is possible, but bias effects are more complex for current Lyapunov analyses. Consequently, these analyses are unsuitable for application to schemes like the mirror Langevin approach. It is within this context that our framework offers an alternative proof methodology, bypassing the necessity to track any Lyapunov function.\n\n\n### Q2\nThis is a very nice question. Indeed, the dependence on $T$ is exponential (see the last equation in line 547). However, the defining property of being a WAPT is that for **any** fixed $T > 0$, as the beginning of the time window $[t, t+T]$ goes to infinity, the sup distance goes to zero; see **Definition 1**. Our analysis then shows that it suffices to have a **finite** (but arbitrary) $T$ in equation (1), and no uniform control in $T$ is required.\n\n\n### Q3\nThe uniqueness of the stationary distribution of the continuous-time Langevin diffusion can be established in various ways under our assumptions. One possible route, in line with our dynamical system analysis, is as follows.\n\n1. Assumption 1 ensures that the continuous-time Langevin diffusion has strong solutions; this is standard. Denote the distribution at time $t$ by $\\mu_t$.\n\n2. Following the notation of M. Benaïm (1999), \"Dynamics of stochastic approximation algorithms\" section 6.2, let $M$ be the set of absolutely continuous probability measures in $W_2$, $\\Lambda =$ {$\\pi$} and define $V(\\mu) = D_{\\mathrm{KL}}(\\mu \\| \\pi)$. Then, it is clear that $V$ is a Lyapunov function of the Langevin dynamics, whose value is strictly decreasing along the flow (as the time derivative of $V$ along the flow is negative of the relative Fisher information, which is strictly positive for all measures other than $\\pi$).\n\n3. Thus, all requirements of Proposition 6.4 are satisfied, showing that every limit set of $\\{\\mu_t\\}_{t\\geq 0}$ is contained in $\\Lambda$. In other words, the only point possible limit point is $\\pi$.\n\n---\nWe hope that the above addresses your questions - but please let us know if any of the above is not sufficiently clear.\n\nThank you again for your input and thoughtful criticisms,\n\nThe authors"}, {"review_id": "Y1xKUdRgdN", "reviewer": "Reviewer_Vq75", "summary": "The authors offer theoretical guarantees on the convergence of the last iterate of a very generic class of sampling methods to the stationary  distribution for a very large classe of sampling schemes in non-convex settings. This is achieved by showing that a large class of discrete sampling schemes can be mapped to continuous time ones and then to a Wasserstein asymptotic pseudo-trajectory. The distribution at large times of this process then converges to the stationary one under quite general assumptions on the Langevin noise and an opportune annealing scheme for the learning rate.", "questions": "1. Within this framework, is it possible to extract how the typical time to convergence scales with the size of the system? (for example in the sense of obtaining a lower bound on the number of steps T(epsilon) to obtain a distance between the distribution after T steps and the stationary one less than epsilon).\n2. Related to the previous point: can you say something quantitative on the performance of ORMM beyond requiring less gradient calls?\n3. Can you comment on the case where sigma(x_k)=1/beta_k in the definition of LRM, where beta_k is positively diverging as k increase?", "limitations": "This paper contains no applications or experiments, as it's expected by a paper of this kind.", "rating": 8, "confidence": 3, "soundness": 3, "presentation": 4, "contribution": 3, "strengths": "The paper addresses the interesting and actively researched problem of sampling non-convex, high dimensional densities. It explores a general framework that can be specialised to a wide range of sampling schemes, thereby providing a valuable guarantee for practical sampling scenarios. \n\nNotably, the authors employ a novel proof technique, which, to the best of my knowledge, is unique in the context of sampling problems. This innovative approach contributes to the originality of the paper and distinguishes it from existing research in the field, which mostly study these problems as gradient flows. \n\nFurthermore, the paper exhibits exceptional writing quality, being highly readable, well-structured, and accessible to a broad audience. These qualities enhance its overall impact and make it easier for readers to grasp the main idea.\n\nConsidering these strengths, I strongly believe that this paper is of high caliber. It offers significant contributions to the field, with its practical relevance, original proof technique, and excellent presentation.", "weaknesses": "I think this is a solid paper without major weaknesses, so I just have a minor remark.\nThe main result of this paper is showing that at very large times the sampling schemes converge to the desired distribution. It would be interesting to relax this condition and obtain bounds on the Wasserstein distance for a large (but finite) number of steps.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your input and remarks. We reply to your questions below, and we will revise our manuscript accordingly in the upcoming revision.\n\n\n\n> The main result of this paper is showing that at very large times the sampling schemes converge to the desired distribution. It would be interesting to relax this condition and obtain bounds on the Wasserstein distance for a large (but finite) number of steps.\n\nFrom our analysis, a doubly exponential bound can be deduced: $\\mathcal{W}_2$ error = $O\\left(\\frac{1}{\\log\\log n}\\right)$. However, we believe that this loose bound offers no significant advantage over asymptotic convergence. Consequently, we have chosen to omit it. It is important to note that obtaining polynomial bounds is infeasible due to NP-hardness.\n\nWe are also able to obtain results that are close in spirit to non-asymptotic results: Assuming step-sizes $\\gamma_n = \\Theta(n^{-1})$ and exponential convergence of continuous-time dynamics (which holds under, e.g., LSI), we can show that an LRM scheme converges in $\\mathcal{W}^2_2$ at a rate of $O(n^{-1})$ under a warm-start condition, **under the same noise and bias conditions**. While we cannot estimate the duration of the burn-in phase, our argument holds for a wider range of stochastic and biased algorithms than in the current literature. \n\n> Within this framework, is it possible to extract how the typical time to convergence scales with the size of the system? (for example in the sense of obtaining a lower bound on the number of steps T(epsilon) to obtain a distance between the distribution after T steps and the stationary one less than epsilon).\n\nLower bounds on the time complexity for sampling problems is a very hard problem. There has been some recent works and progress is being made for specific algorithms and simpler problem classes, see, e.g., Chatterji et al. (2021) \"Oracle Lower Bounds for Stochastic Gradient Sampling Algorithms\" and Chewi et al. (2023) \"Query lower bounds for log-concave sampling\". \n\nIn our setup, the class of target distributions is rich enough to encompass NP-Hard problems, implying that the most difficult problems within this class require at least exponential time to solve.\n\n> Related to the previous point: can you say something quantitative on the performance of ORMM beyond requiring less gradient calls?\n\nUnfortunately, we are unable to provide such comparisons in theory, although we point out that the lack of such quantitative bounds for ORMM is not specific to our framework: It is present already in the context of purely **deterministic convex-concave optimization**, where optimistic methods originate, since their convergence rates do not surpass those of their non-optimistic counterparts. In most cases, the advantages of these methods are largely empirical. \n\nIn this regard, our primary contribution lies in offering flexibility: By recycling past gradients, one can preserve the asymptotic convergence of randomized mid-point, demonstrating the potential to maintain favorable convergence behavior while saving 50% of the per-iteration costs.\n\n\n> Can you comment on the case where sigma(x_k)=1/beta_k in the definition of LRM, where beta_k is positively diverging as k increase?\n\nThis is a very nice question. Diffusions with vanishing diffusion coefficient usually arise from simulated annealing processes. The problem with these time-inhomogeneous systems is that the flow (the $\\Phi$ in the paper) is not going to be well-defined. We can make two possible comments here: \n\n1. One can take $\\Phi$ to be the flow of the resulting *ODE* when $\\beta_k \\to \\infty$ (rather than the SDE), and the LRM scheme with the solutions of the ODE. This is a similar approach as in our reference [6] for \"asymptotically autonomous with limit equation\" (see section 2 of [6]). Our analysis goes through without any changes, and the result will be \"$(X_t)$ will almost surely converge to an ICT set of the ODE; in the case of annealed Langevin dynamics as the main process, this means that $(X_t)$ converges to a *stationary point* of the potential.\" This result is, however, unsatisfactory, as one usually performs simulated annealing to converge to a **global optimum**.\n\n2. One can convert the time-inhomogeneous system to a time-homogeneous system (by, e.g., introducing a new variable for tracking time). This dynamics will have a well-defined flow. However, this flow does not have any equilibria (as the time variable goes to infinity).\n\nWe will definitely explore this direction in our future works.\n\n\n---\nWe hope that the above addresses your questions - but please let us know if any of the above is not sufficiently clear.\n\nThank you again for your input and positive evaluation,\n\nThe authors"}, {"review_id": "3qAkt6fTnQ", "reviewer": "Reviewer_3YTY", "summary": "This paper presents a general framework for studying the convergence of last-iterate, noisy, and possibly biased Langevin-like discrete time approximations to the continuous Langevin flow for sampling from a distribution. ", "questions": "What is the theorem of [6] being invoked here? Can the authors explain its conditions and why they are met.", "limitations": "This work is purely asymptotic, and provides no rates of convergence. ", "rating": 7, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "This paper is (with one small quibble which I will explain below) very well written and well explained. Using the machinery of asymptotic pseudotrajectories, it establishes asymptotic convergence of a wide range of discrete time sampling schemes to the Boltzman distribution $e^{-f}$. The advantage of the pseudotrajectory framework is that it replaces standard \"Cauchy-type\" convergence (which requires checking that an infinitely long tail of iterates converge) with a weaker notion of convergence that only requires \"finite length\" tails. \n\nThe upshot is that, due to cited source [6], this weaker sense of convergence is sufficient for convergence to a limit point of the original flow $\\Phi$, provided that $\\Phi$ has a unique fixed point. This allows for much more freedom in the analysis, and hence encompasses a wider range of conditions. The technique also has advantages of discretization-based approaches as explained by the author.", "weaknesses": "My biggest complain with the paper is that perhaps the most interesting part of the framework is the 2-line proof of Theorem 2 relegated to the appendix. This argument relies on  the limit-set characterization of [6] to show that asymptotic pseudotrajectories are sufficient to ensure convergence to the fixed point of the Langevin flow, namely the correct Boltzman distribution. I think the authors should explain this point in the body, and explain which theorem of [6] is being invoked (I suspect Thm 0.1), and how it's conditions are met. That way, the reader understands not simply how the proof framework of the paper works, but why we should believe that is somehow the \"correct\" one. ", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your input and remarks. We reply to your questions below, and we will revise our manuscript accordingly in the upcoming revision.\n\n\n\n> My biggest complain with the paper is that perhaps the most interesting part of the framework is the 2-line proof of Theorem 2 relegated to the appendix. This argument relies on the limit-set characterization of [6] to show that asymptotic pseudotrajectories are sufficient to ensure convergence to the fixed point of the Langevin flow, namely the correct Boltzman distribution. I think the authors should explain this point in the body, and explain which theorem of [6] is being invoked (I suspect Thm 0.1), and how it's conditions are met. That way, the reader understands not simply how the proof framework of the paper works, but why we should believe that is somehow the \"correct\" one.\n\nWe agree with the reviewer. We only made this choice due to the page limit. We will expand the proof of the theorem to include more details. Indeed, Theorem 0.1 (see Theorem 5.7 (i) in M. Benaïm (2006), \"Dynamics of stochastic approximation algorithms\" for an extended version) is used. \n\n> What is the theorem of [6] being invoked here? Can the authors explain its conditions and why they are met.\n\nTheorem 0.1 is used here. The assumptions are:\n\n- $(\\mathrm{law}(X_t))_t$ is an APT in the corresponding metric space (which we prove in Theorem 1 in the paper)\n- $(\\mathrm{law}(X_t))_t$ is precompact (which we prove in Theorem 2 that it is equivalent to Assumption (10))\n- $\\Phi$ is integrable (which is implied by Assumption 1 in our paper)\n\nThese imply that the limit-set of $(\\mathrm{law}(X_t))_t$ is an *internally chain transitive* set. \n\nThese conditions are met becuase Assumption (10) implies that the trajectory $(X_t)$ is precompact in the Wasserstein space. If the trajectory is also a Wasserstein APT, then it is guaranteed to converge to an ICT set of the flow corresponding to the SDE. For the case of the Langevin SDE, we can show that the ICT set is the singleton $\\{ \\pi \\}$. Following the notation of M. Benaïm (2006), \"Dynamics of stochastic approximation algorithms\" section 6.2, let $M$ to be the set of absolutely continuous probability measures in $W_2$, $\\Lambda = \\{\\pi\\}$ and define $V(\\mu) = D_{\\mathrm{KL}}(\\mu \\mid \\pi)$. Then, it is clear that $V$ is a Lyapunov function of the Langevin dynamics, whose value is strictly decreasing along the flow (as the time derivative of V along the flow is negative of the relative Fisher information, which is strictly positive for all measures other than $\\pi$). Thus, all requirements of Proposition 6.4 are satisfied, showing that every ICT set is contained in $\\Lambda$. In another words, the only point in the ICT set is $\\pi$.\n\n---\n\nWe hope that the above addresses your questions - but please let us know if any of the above is not sufficiently clear.\n\nThank you again for your input and positive evaluation,\n\nThe authors"}]

Author rebuttal:
Dear AC and dear reviewers,

We wish to express our sincere gratitude for your dedicated efforts. Your insightful critiques and favorable evaluation have been acknowledged, and we have responded to all your inquiries in a detailed point-by-point manner, presented below.

After thorough consideration of your remarks, we wish to affirm the validity of the concerns highlighted by the reviewers. We firmly believe that these concerns can be aptly addressed through minor revisions, as detailed in our individual response. We appreciate your invaluable contributions towards refining our manuscript.


With utmost appreciation,

The authors


Meta-review: This paper considers the convergence of discretized Langevin diffusion using Robbins-Monro step sizes schedules towards the target Gibbs measure. Their framework can cover a variety of numerical discretization schemes, and provides asymptotic analysis under mild assumptions. Authors accomplish this by taking a standard interpolation argument and modifying it to construct a continuous version of the trajectory and verifying the stability conditions from the dynamical systems theory. The scope covers LMC, mirror LMC, proximal LMC, Randomized Mid-Point, and Runge-Kutta methods.

This paper is reviewed by 4 expert reviewers and after the discussion period the paper received the following Rating/Confidence scores: 5/3, 8/3, 7/3, 7/4. All four of the reviewers agree that the paper makes good contributions and it is above the bar. Majority of the reviewers found the contributions significant. AC also agrees with the reviewers and recommends including this paper to conference program.
Decision: Accept (spotlight)