Wasserstein mixing time of the unadjusted Langevin algorithm

Wasserstein mixing time of the unadjusted Langevin algorithm Francesco Pedrotti and Peter A. Whalley ETH Zürich August 4, 2026

arXiv:2608.02430v1 [stat.CO] 3 Aug 2026

Abstract We provide new estimates in Wasserstein distance for the asymptotic bias of the unadjusted Langevin algorithm, in the classical setting of log-smooth √ strongly log-concave measures. Our bound implies a Wasserstein mixing time of order κ d/ε, where κ is the condition number, √ d is the dimension, and ε is the target precision: this improves by a factor of d/ε over the previous state-of-the-art results.

1

Introduction

Sampling from a probability density π ∝ e−V in Rd is a fundamental problem in a variety of fields. A popular family of methods to produce these samples, known as MCMC algorithms, consists in designing a stochastic process which is ergodic with respect to π, and simulating it for a long enough time. One of the most fundamental choices is the Langevin dynamics (LD): given
d a (possibly) random initialization µ ∈ P2 R and a standard d-dimensional Brownian motion (Bt )t , it is given by the solution of the SDE √ (1) X0 ∼ µ, dXt = −∇V (Xt )dt + 2dBt . It is well-known that, under mild assumptions on V , this stochastic process is ergodic with respect to π, making it a natural candidate for an MCMC algorithm (see [Pav14]). Its practical implementation, however, requires considering a numerical scheme for (1), such as the Euler– Maruyama method, which leads to the unadjusted Langevin algorithm (ULA), described by the easily implemented recursion √
X̂0 ∼ µ, X̂k+1 = X̂k − h∇V (X̂k ) + 2 B(k+1)h − Bkh , (2) for a small step-size h > 0. Unfortunately, the above scheme is biased, which means that its limiting stationary distribution π̂ is in general different from π; this leads to a trade-off in the choice of the step-size, with smaller values of h leading to a smaller error in the bias, but to a higher required number of iterations of (2) in order to converge to π̂. Theoretical guarantees for ULA In recent years, considerable effort has been put into rigorously quantifying this trade-off [Che26], by proving error bounds in the approximation  µk := law X̂k ≈ π in different metrics, such as the Wasserstein distance, defined for distributions µ, ν ∈ P2 (Rd ) W2 (µ, ν) =

r h

inf (X,Y )∼(µ,ν)

i E |X − Y |2 ,

where the infimum runs over all Rd -valued random vectors X

and Y defined on the same probability space with law(X) = µ and law(Y ) = ν. Our main result belongs to this line of work, and holds under the following classical assumption on the target measure. 1

Assumption 1.1. The potential V ∈ C 2 Rd satisfies αId ≼ ∇2 V ≼ βId , for some 0 < α ≤ β. We denote the condition number by κ := αβ . Theorem 1. Under Assumption 1.1, if h ≤ β1 , we have p W2 (π, π̂) ≤ 6h κβd.

(3)

Theorem 1 immediately leads to convergence guarantees for the iteration complexity of ULA to achieve a desired accuracy in Wasserstein distance; more precisely, in this setting it is standard √ and more meaningful to work with the rescaled metric αW2 , in order to have scale-invariant results, cf. [Che26, Remark 4.1.4]. Under Assumption 1.1, for h ≤ β1 ULA converges exponentially fast to its biased limit, i.e. writing µk = law(X̂k ) we have W2 (µk , π̂) ≤ (1 − αh)k W2 (µ, π̂). Thus, √ the iteration complexity to achieve a desired accuracy in αW2 -distance is essentially determined by the asymptotic bias and the step-size that it imposed, and we deduce from Theorem 1 the following Corollary. Corollary 1.1. Under Assumption 1.1, if h ≤ β1 , we have p W2 (µk , π) ≤ (1 − αh)k W2 (µ, π̂) + 6h κβd.

(4)

In particular, for any ε > 0, choosing h=

1 β

 1∧

ε √ 12 d



it suffices to take

√ !   √ 12 d 2 αW2 (µ, π̂) k ≥κ 1∨ log 1 + ϵ ϵ √ iterations of ULA to ensure αW2 (µk , π) ≤ ε.

(5)

We remark that in (3) the dependence on h and d can be shown to be optimal via explicit computations when the target is Gaussian; cf. [Wib18, Example 2]. Additionally, the convergence time of (2) towards its biased limit depends linearly on the condition number, already for anisotropic Gaussian targets. As a consequence, the linear dependence on the condition number in (5) is also optimal. Literature overview Assumption 1.1 is classical in the study of the quantitative guarantees of the Unadjusted Langevin Algorithm (see [Dal17b; DM17] and the book [Che26] for a more complete overview). Theoretical guarantees in Wasserstein distance under √ this assumption were first provided in [DM19; Dal17a], √ which establish the bound W2 (π, π̂) ≲ κ hd. This bound was later improved to W2 (π, π̂) ≲ κhd in [DMM19]. Linear dependence in h for the asymptotic bias was only established under additional assumptions on the third derivatives of V , but with a worse dependence on the dimension cf. [DM19; DK19]. Subsequently, the dimension dependence was improved, but only under even stronger assumptions on the third derivatives of V (see [LZT22; DE24; DK26]). Besides Wasserstein distance, convergence guarantees have also been provided in other metrics, such as total variation, cf. [Dal17b; DM17], and relative entropy, see [VW23; Che+25], which also relaxes the log-concavity assumption to appropriate functional inequalities. Additionally, recent work [Che+26] showed under suitable sparsity assumptions that dimension-free guarantees (up to logarithmic factors) can be obtained, but this crucially relies on replacing the 2-Wasserstein distance with an ℓ∞ -version. 2

2

Proof of Theorem 1

This section is devoted to the proof of Theorem 1. We first prove the theorem under the additional regularity Assumption 2.1 below, and remove it in Section 2.2 with an approximation argument. Assumption 2.1. The potential V is smooth and all of its partial derivatives of all orders have at most polynomial growth. This additional assumption is considered for simplicity, to make rigorous and more transparent the application of standard semigroup and spectral theory in our argument. For the convenience of the reader, we briefly recall some of these useful facts in the setting of Assumptions 1.1 and 2.1, and refer to the book [BGL14], the lecture notes [KL25] and the papers [KP23; Kla23] for details. First, the Laplace operator L associated to π is initially defined for f ∈ Cc∞ (Rd ) via Lf = ∆f − ⟨∇V, ∇f ⟩. From the definition and integration by parts, it follows that for f, g ∈ Cc∞ (Rd ) Z ⟨Lf, g⟩L2 (π) = − ⟨∇f, ∇g⟩ dπ, which shows that −L is a symmetric and monotone operator. The operator −L is essentially self-adjoint, and in what follows we denote by −L its minimal extension (with a mild abuse of notation), and by D(L) its domain. Under our assumptions, −L has a discrete spectrum, 0 = λ0 < λ1 ≤ . . . 2 and a corresponding orthonormal basis (φk )∞ k=0 in L (π), where φ0 = 1. The connection with the Langevin dynamics comes from the fact that L can also be identified with the L2 -generator of the associated Markov semigroup (Pt )t , which is initially defined probabilistically for bounded measurable functions f by

Pt f (x) = E[f (Xt ) | X0 = x], and then extended to a contraction Pt : L2 (π) → L2 (π) by Jensen’s inequality. Indeed, one can check that for f ∈ D(L) we have Pt f − f → Lf t

in L2 (π) as t ↓ 0,

and D(L) is precisely the set of functions f ∈ L2 (π) for which the left hand-side above admits a limit in L2 (π). This suggests the identity Pt f = etL f , which holds in fact for all f ∈ L2 (π). Here, etL : L2 (π) → L2 (π) is defined as an operator on L2 (π), even if the operator L is defined only on D(L) ⊂ L2 (π): more generally, for every smooth F : R → R bounded on (−∞, 0], we can define F (L) : L2 (π) → L2 (π) with spectral calculus, setting for f ∈ L2 (π) [F (L)]f =

∞ X

⟨f, φk ⟩L2 (π) F (−λk )φk .

k=0

3

2.1

Proof for regular distributions

In this subsection, we prove Theorem 1 under Assumption 2.1. Fix a step-size 0 < h ≤ β1 and consider synchronous coupling of (1) and (2) initialized at X0 = X̂0 ∼ π (i.e. use the same driving Brownian motion). Setting Zk = Xkh − X̂k , our aim is to prove that r h i p lim sup E |Zk |2 ≤ 6h βκd, (6) k→∞

from which the desired conclusion follows since Xt ∼ π for all t ≥ 0 and X̂k converges to π̂ in W2 . Observe that after subtracting (2) from (1), Zk satisfies the recursion Zk+1 = (Id − hAk )Zk + Dk ,

Z0 = 0, where

Z 1 Ak :=

∇2 V (X̂k + tZk ) dt,

αId ⪯ Ak ⪯ βId ,

(7) (8)

0

and

Z h {∇V (Xkh ) − ∇V (Xkh+r )} dr,

Dk :=

(9)

0

which we refer to as the local error at iteration k. Local error decomposition The key idea is to decompose each local error (9) into a telescoping increment and a martingale increment. At stationarity, sums of the local errors are expected 2 to √ exhibit central-limit-theorem scaling [Gor69; KV86; MT12]: their L -norm grows at rate k, rather than at the naive rate k, where k is the number of iterations (and summands). The decomposition makes this cancellation explicit, as the telescoping increments cancel across successive steps, whilst the conditionally centred martingale increments accumulate in L2 at rate √ k. We exploit both effects in the global error analysis below. We now derive the desired decomposition: denoting by Ft = σ((Xs )0≤s≤t ) the σ-algebra generated by the diffusion up to time t, we claim that there exist a vector field g ∈ L2 (π; Rd ) and random vectors Mk+1 such that Dk = g(Xkh ) − g(X(k+1)h ) + Mk+1 , and

∥g∥2L2 (π) ≤ h2 βd,

E[Mk+1 | Fkh ] = 0,

E|Mk+1 |2 ≤ 18 h3 β 2 d.

(10) (11)

To achieve this, the idea is to take conditional expectations with respect to Fkh in (10), which gives Z h (I − Ph )g(Xkh ) = E[Dk | Fkh ] = (I − Ps )∇V (Xkh ) ds. 0

We can then identify g by solving the discrete Poisson equation [MT12] Z h (I − Ph )g = (I − Ps )∇V ds. 0

More precisely, consider the function G ∈ C ∞ (R) defined by ( G(x) =

h+(1−ehx )/x 1−ehx h 2

x ̸= 0 x = 0,

which satisfies h2 ≤ G(x) ≤ h for x ≤ 0. We then define g = G(L)∇V and observe that Z ∥g∥2L2 (π) ≤ h2 ∥∇V ∥2L2 (π) = h2 ∆V dπ ≤ h2 βd, 4

(12)

which is the first bound in (11). Next, notice that Z h E[Dk | Fkh ] = (∇V (Xkh ) − E[∇V (Xkh+s ) | Fkh ]) ds 0 Z h (I − Ps )∇V (Xkh ) ds = 0 Z h (I − esL )∇V (Xkh ) ds = 0 nh i o = (I − ehL )G(L) ∇V (Xkh ) = g(Xkh ) − Ph g(Xkh )   = E g(Xkh ) − g(X(k+1)h ) | Fkh .
Thus, defining Mk+1 := Dk − g(Xkh ) − g(X(k+1)h ) yields the desired decomposition in (10). It h i remains to bound E |Mk+1 |2 : to this end, we first bound h i D E 2 E g(X(k+1)h ) − g(Xkh ) = 2 g, (I − ehL )g 2 L (π) D E = 2 G(L)∇V, (I − ehL )(G(L)∇V ) 2 L (π) D E 2 hL ≤ 2h ∇V, (I − e )∇V 2 L (π) h i 2 = h2 E ∇V (X(k+1)h ) − ∇V (Xkh ) ≤ 6h3 β 2 d. where we have used that h i
E |∇V (Xkh+s ) − ∇V (Xkh )|2 ≤ β 2 d 2s2 β + 4s ≤ 6sβ 2 d, (13) h i for 0 < s ≤ β1 . Additionally using (13) we can also bound E |Dk |2 ≤ 3h3 β 2 d, and hence using h i triangle inequality and (10) we have that E |Mk+1 |2 ≤ 18 h3 β 2 d, which proves (11). Global error bound We can now insert decomposition (10) for Dk in the recursion (7) for Zk : after rearranging, we deduce that
Zk+1 + g X(k+1)h = (Id − hAk )(Zk + g(Xkh )) + hAk g(Xkh ) + Mk+1 (14) Taking the square and expectation on both sides, we obtain h h i h i h i
2i 2 2 2 E Zk+1 + g X(k+1)h ≤ (1 − αh)E |Zk + g(Xkh )| + hβ E |g(Xkh )| + E |Mk+1 | h i (15) ≤ (1 − αh)E |Zk + g(Xkh )|2 + h3 β 2 d + 18 h3 β 2 d, where we used (11) and that |(I − hAk )u + hAk v|2 ≤ (1 − αh)|u|2 + hv ⊤ Ak v for every u, v ∈ Rd since hβ ≤ 1. From (15) it follows that h i lim sup E |Zk + g(Xkh )|2 ≤ 19 h2 κβd. k→∞

Therefore, lim sup k→∞

r h

2

E |Zk |

i

≤ lim sup

r h

2

E |Zk + g(Xkh )|

k→∞

which is the desired bound (6). 5

i

+

r h

2

E |g(Xkh )|

! i

p ≤ 6h κβd,

2.2

Removing Assumption 2.1

We finally remove the extra regularity assumption 2.1, thus concluding the proof of Theorem 1. For any small ε > 0, we denote by γε := N (0, εId ) the density of a Gaussian random variable with mean 0 and covariance εId , and we set πε = π ∗ γε . Then, Vε := − log πε ∈ C ∞ (Rd ) and α Id ⪯ ∇2 Vε ⪯ βId , αε + 1 see [SW14] and equation (6) in [MS23] for the lower and upper bound respectively [KP23]. Moreover, Assumption 2.1 is satisfied too, cf. [KP23, Lemma 2.2]. For h ≤ β1 , we denote by π̂ε the biased limit of ULA with step-size h and targeting πε . By the triangle inequality, W2 (π, π̂) ≤ W2 (π, πε ) + W2 (πε , π̂ε ) + W2 (π̂ε , π̂). p An application of the previous yields W2 (πε , π̂ε ) ≤ 6h (1 + αε)κβd. Moreover, letting r step h i √ Z ∼ γε , we have W2 (π, πε ) ≤ E |Z|2 ≤ εd. Thus, to conclude the proof, it suffices to show
that W2 (π̂ε , π̂) → 0 as ε → 0. To see this, fix any 0 < ε < min 1, α−1 , let X̂0 = Ŷ0 ∼ π̂ and consider synchronous coupling for ULA targeting π, πε respectively, i.e. √
X̂k+1 = X̂k − h∇V (X̂k ) + 2 B(k+1)h − Bkh √
Ŷk+1 = Ŷk − h∇Vε (Ŷk ) + 2 B(k+1)h − Bkh ,     Defining Ẑi = X̂i − Ŷi , it follows by the bounds on ∇2 Vϵ that i∈N

i∈N

r h

E |Ẑk+1 |2

i

 r h i p αh ≤ 1− E |Ẑk |2 + h Eπ̂ [|∇Vε − ∇V |2 ]. 2

Taking the limit as k → ∞ this implies W2 (π̂ε , π̂) ≤

2p Eπ̂ [|∇Vε − ∇V |2 ]. α

Observe now that ∇Vε → ∇V pointwise as ε → 0. Moreover, for ε small enough and some constant C > 0 depending only on π, we have |∇Vε (x) − ∇V (x)| ≤ C(1 + |x|) ∈ L2 (π̂), since π̂ is subgaussian [VW23]. Thus, by dominated convergence, we have p Eπ̂ [|∇Vε − ∇V |2 ] → 0 as ε → 0, which implies that W2 (π̂ε , π̂) → 0 too, as desired. Acknowledgments We would like to thank Yuansi Chen for inspiring discussions during preparation of another work, and acknowledge the use of ChatGPT 5.6 for complementary proof checking and literature exploration.

References [BGL14]

Dominique Bakry, Ivan Gentil, and Michel Ledoux. Analysis and geometry of Markov diffusion operators. Vol. 348. Grundlehren der mathematischen Wissenschaften [Fundamental Principles of Mathematical Sciences]. Springer, Cham, 2014, pp. xx+552. 6

[Che+25]

Sinho Chewi, Murat A. Erdogdu, Mufan Li, Ruoqi Shen, and Matthew S. Zhang. “Analysis of Langevin Monte Carlo from Poincaré to log-Sobolev”. In: Found. Comput. Math. 25.4 (2025), pp. 1345–1395.

[Che+26]

Yifan Chen, Xiaoou Cheng, Jonathan Niles-Weed, and Jonathan Weare. “Convergence of unadjusted Langevin in high dimensions: delocalization of bias”. In: Comm. Pure Appl. Math. 79.6 (2026), pp. 1467–1491.

[Che26]

Sinho Chewi. Log-concave sampling. Available online at https://chewisinho.github. io/. Forthcoming, 2026.

[Dal17a]

Arnak Dalalyan. “Further and stronger analogy between sampling and optimization: Langevin Monte Carlo and gradient descent”. In: Conference on Learning Theory. PMLR. 2017, pp. 678–689.

[Dal17b]

Arnak S. Dalalyan. “Theoretical guarantees for approximate sampling from smooth and log-concave densities”. In: J. R. Stat. Soc. Ser. B. Stat. Methodol. 79.3 (2017), pp. 651–676.

[DE24]

Alain Durmus and Andreas Eberle. “Asymptotic bias of inexact Markov chain Monte Carlo methods in high dimension”. In: Ann. Appl. Probab. 34.4 (2024), pp. 3435–3468.

[DK19]

Arnak S. Dalalyan and Avetik Karagulyan. “User-friendly guarantees for the Langevin Monte Carlo with inaccurate gradient”. In: Stochastic Process. Appl. 129.12 (2019), pp. 5278–5311.

[DK26]

Arnak S Dalalyan and Avetik Karagulyan. “Improved Guarantees for Langevin Monte Carlo with Average Smoothness”. In: arXiv preprint arXiv:2605.31413 (2026).

[DM17]

Alain Durmus and Éric Moulines. “Nonasymptotic convergence analysis for the unadjusted Langevin algorithm”. In: Ann. Appl. Probab. 27.3 (2017), pp. 1551–1587.

[DM19]

Alain Durmus and Éric Moulines. “High-dimensional Bayesian inference via the unadjusted Langevin algorithm”. In: Bernoulli 25.4A (2019), pp. 2854–2882.

[DMM19] Alain Durmus, Szymon Majewski, and Błażej Miasojedow. “Analysis of Langevin Monte Carlo via convex optimization”. In: J. Mach. Learn. Res. 20 (2019), Paper No. 73, 46. [Gor69]

M. I. Gordin. “The central limit theorem for stationary processes”. In: Dokl. Akad. Nauk SSSR 188 (1969), pp. 739–741.

[KL25]

Boaz Klartag and Joseph Lehec. “Isoperimetric inequalities in high-dimensional convex sets”. In: Bull. Amer. Math. Soc. (N.S.) 62.4 (2025), pp. 575–642.

[Kla23]

Bo’az Klartag. “Logarithmic bounds for isoperimetry and slices of convex sets”. In: Ars Inven. Anal. (2023), Paper No. 4, 17.

[KP23]

Bo’az Klartag and Eli Putterman. “Spectral monotonicity under Gaussian convolution”. In: Ann. Fac. Sci. Toulouse Math. (6) 32.5 (2023), pp. 939–967.

[KV86]

C. Kipnis and S. R. S. Varadhan. “Central limit theorem for additive functionals of reversible Markov processes and applications to simple exclusions”. In: Comm. Math. Phys. 104.1 (1986), pp. 1–19.

[LZT22]

Ruilin Li, Hongyuan Zha, and Molei Tao. “Sqrt (d) dimension dependence of Langevin Monte Carlo”. In: The International Conference on Learning Representations. 2022.

[MS23]

Dan Mikulincer and Yair Shenfeld. “On the Lipschitz properties of transportation along heat flows”. In: Geometric aspects of functional analysis. Vol. 2327. Lecture Notes in Math. Springer, Cham, [2023] ©2023, pp. 269–290.

[MT12]

Sean P Meyn and Richard L Tweedie. Markov chains and stochastic stability. Springer Science & Business Media, 2012. 7

[Pav14]

Grigorios A. Pavliotis. Stochastic processes and applications. Vol. 60. Texts in Applied Mathematics. Diffusion processes, the Fokker-Planck and Langevin equations. Springer, New York, 2014, pp. xiv+339.

[SW14]

Adrien Saumard and Jon A. Wellner. “Log-concavity and strong log-concavity: a review”. In: Stat. Surv. 8 (2014), pp. 45–114.

[VW23]

Santosh S. Vempala and Andre Wibisono. “Rapid convergence of the unadjusted Langevin algorithm: isoperimetry suffices”. In: Lecture Notes in Math. 2327 ([2023] ©2023), pp. 381–438.

[Wib18]

Andre Wibisono. “Sampling as optimization in the space of measures: The Langevin dynamics as a composite optimization problem”. In: Conference on learning theory. PMLR. 2018, pp. 2093–3027.

8

Related documents

• Statistical finite elements via interacting particle Langevin dynamics#295964
• Rényi Tracking Bounds for Langevin Dynamics with Moving Targets#1011310
• Density functions for the overdamped generalized Langevin equation and its Euler--Maruyama method: smoothness and convergence#925711
• Fast Convergence for High-Order ODE Solvers in Diffusion Probabilistic Models#999662
• The high-dimensional asymptotics of first order methods with random data#152953