Advancing Wasserstein Convergence Analysis of Score-Based Models: Insights from Discretization and Second-Order Acceleration
Authors: Yifeng Yu
Lu Yu
Venue: NeurIPS 2025 poster ()
Keywords: Diffusion model
Wasserstein convergence
Convergence Analysis
Primary area: theory

Abstract: Score-based diffusion models have emerged as powerful tools in generative modeling, yet their theoretical foundations remain underexplored. In this work, we focus on the Wasserstein convergence analysis of score-based diffusion models. Specifically, we investigate the impact of various discretization schemes, including Euler discretization, exponential integrators, and midpoint randomization methods. Our analysis provides the first quantitative comparison of these discrete approximations, emphasizing their influence on convergence behavior. Furthermore, we explore scenarios where Hessian information is available and propose an accelerated sampler based on the local linearization method. We establish the first Wasserstein convergence analysis for such a Hessian-based method, showing that it achieves an improved convergence rate of order $\widetilde{\mathcal{O}}\left(\frac{\sqrt{d}}{\varepsilon}\right)$, which significantly outperforms the standard rate $\widetilde{\mathcal{O}}\left(\frac{d}{\varepsilon^2}\right)$ of vanilla diffusion models. Numerical experiments on synthetic data and the MNIST dataset validate our theoretical insights.

Advancing Wasserstein Convergence Analysis of Score-Based Models: Insights from Discretization and Second-Order Acceleration
Yifeng Yu Tsinghua University yyf220mails ttsinguua edw.. edu. cn
Lu Yu City University of Hong Kong lu. yu@cityu..ueeee edu hk
Abstract
Score-based diffusion models have emerged as powerful tools in generative model- ing, yet their theoretical foundations remain underexplored. In this work, we focus on the Wasserstein convergence analysis of score- based diffusion models. Specifi- cally, we investigate the impact of various discretization schemes, including Euler discretization, exponential integrators, and midpoint randomization methods. Our analysis provides the first quantitative comparison of these discrete approximations, emphasizing their influence on convergence behavior. Furthermore, we explore scenarios where Hessian information available and propose an accelerated sam- pler based on the local linearization method. We establish the first Wasserstein convergence analysis for such a Hessian- based method, showing that it achieves an improved convergence rate of order ((Vd/e), which significantly outperforms the standard rate (dde²² =²) of vanilla diffusion models. Numerical experiments on synthetic data and the MNIST dataset validate our theoretical insights.
Introduction
Diffusion models have become pivotal framework in modern generative modeling, achieving notable success across fields such as image generation [TT [7 31 34), natural language processing (30), and computational biology [ 42). These models add noise data via forward process and learn to reverse it, reconstructing data from noise. This approach enables them to capture the underlying structure of complex, high-dimensional data distributions. For detailed review of diffusion models, we refer the readers to
39 43)
A widely adopted formulation of diffusion models is the score-based generative model (SGM). implemented using stochastic differential equations (SDEs) [36]. Broadly speaking, SGMs rely on two key stochastic processes: forward process that gradually transforms data samples into pure noise, and backward process that reverses this transformation, recovering the target data distribution from noise.
Despite the remarkable empirical success of diffusion models across various applications, their theoretical understanding remains limited. In recent years, there has been rapidly expanding body of research on the convergence theory of diffusion models. Generally, these contributions can be divided into two main approaches. each focusing on different metrics and divergences. The first category investigates convergence bounds based on a-divergence, including the Kullback-Leibler (KL) divergence and the total variation (TV) distance (see e.g., 7,77 8 22 26 41).. Among these works, several explore acceleration techniques that leverage higher-order information about log The second category focuses on convergence bounds in Wasserstein
density (see e.g., [9 24 27
"Corresponinng author
39th Conference on Neural Information Processing Systems (NcurIPS 2025).

---

distance, which is often considered more practical and informative for estimation tasks. One line of work within this category assumes strong con-coniyv of the data distribution and access to accurate estimates of the score function 3 13
37
38 Another line of work focuses on specific structural
assumptions of the data distribution 10
14 29
Much of the existing literature on the convergence theory of diffusion models relies on the Euler discretization method. Notably, [5] compare the behavior of Euler discretization and exponential integrators [8 in terms of KL divergence. Additionally, [T0] provide a comparative analysis of these two schemes, though without formal theoretical guarantees. A comprehensive and systematic understanding of how different discretization schemes influence convergence performance in diffusion models remains underexplored. Furthermore, while convergence analyses of accelerated diffusion models primarily focus on TV or KL distances, studies investigating Wasserstein convergence for these accelerations remain lacking
In this work, we address these challenges by analyzing the Wasserstein convergence of score-based dif- fusion models when the data distribution has smooth and strongly log-concave density. Specifically, we investigate the impact of different discretization schemes on convergence behavior. Beyond the widely used Euler method and exponential integrator, we explore the midpoint randomization method. This method was initially introduced in [32) for discretizing kinetic Langevin diffusion 9 and then has been extensively studied in log-concave sampling complexity theory [66 21
46 was later
applied to diffusion models [I5] 223. showing improved KL and TV convergence performance over vanilla models and offering easy parallelization.
We also consider scenarios where accurate estimates of the Hessian of the log density are accessible. Inspired by [3), we propose a novel sampler based on the local linearization method. Our analysis shows that this approach significantly accelerates convergence in Wasserstein distance.
Our contribution can be summarized as follows.
We establish convergence guarantees for SGMs in the Wasserstein-2 distance under various discretization methods, including the Euler method, exponential integrators, the midpoint randomization method, and a hybrid approach combining the latter two.
We introduce novel Hessian-based accelerated sampler for the stochastic diffusion process, leveraging the local linearization method. We then establish its Wasserstein convergence analysis in Theorem achieving state-of-the-ar order of Ö(vd/E).
Section 5 compares the performance of SGMs under four discretization schemes and the proposed Hessian bease method on both synthetic data and the MNIST dataset. The results align with our theory and highlight the acceleration of the proposed second-order method.
In summary, our analysis provides a quantitative comparison of different discrete approximations, offering practical guidance for choosing discretization. Moreover, we present the first Wasserstein convergence analysis of an accelerated sampler that leverages the second-order information about log-densities. This accelerated sampler achieves a faster convergence rate Ö(Vd/E), compared to the standard rate Ö(d/=?) of vanilla diffusion models. These results contribute the understanding of Wasserstein convergence in score-based models, shedding light on aspects that have not been extensively explored before.
Notation. Let RRd be the d-dimensional Euclidean space and Ia the identity matrix. The gradient and the Hessian of a function f Rdd - RR. are denoted by Vf and 22f. Given any pair of measures and v, the W distance between μ and y is defined as
1/2
inf / We(u,v) (eET(a))) RR
| - vll de(z,))
where the infimum is taken over all joint distributions e that have μ and v as marginals. For two symmetric d x d matrices A and B, we use A & B or B & A to denote the relation that B- A is positive semi-definite. For any random vector X, let L(X) denote its law, and define || == [P] where || || denotes the Fuclidean norm. For any random matrix A, define |Al|L : E A where || llz denotes the matrix 2-norm.
'Due to space constraints, we defer the discussion of related theoretical advancements in accelerating samplers for diffusion models to Appendix|A

---

2 Background and Our Setting
Framework. We consider the forward process
dx; ((x,, dt + t)
(I)
where the initial point Xo ~ Po follows the data distribution, and B, denotes the standard d- d-dimensional Brownian motion. Here, the drift f Rd R 7 R. and the function g : R x R - Rddd are diffusion parameters. Some conditions are necessary to ensure that the SDE T well-defined. In practice, various choices for the pair (J.g) are employed, depending on the specific needs of the model; for a detailed survey, we refer to [39]. For clarity, we adopt the simplest possible choice in this work by setting ((x,, the -uaieinn process, which is described by the following SDE:
X,/2 and g (Xx, 1) 1. This results in
dx, =
2 X, dt + dB,
(2)
The forward process is run until sufficiently large time T 0, at which point the corrupted marginal distribution of Xr, denoted by pr. is expected to approximate the standard Gaussian distribution. Then, diffusion models generate new data by reversing the SDE which leads to the following backward SDE
ax; = 2x, (X; + 2V Po dt dW.,
(3)
where X5- ~ ~PT, and the term V log Pt, referred to the score function for pt. is represented by the gradient of the log density function of pt. Additionally, W, denotes another standard Brownian motion independent of B,. Under mild conditions, when initialized X Pr, the backward process {X, JostEr has the same distribution as the forward process {Xx-tJugter [2.4 As a result, running the reverse diffusion X, from = 0 to T will generate sample from the target data distribution po- Note that the density pr is unknown; we approximate it using the distribution
as proposed in Therefore, we derive reverse diffusion process defined by
1 dY, = 2 IV dt + aW,, Yo pr.
Score Matching, Another challenge in working with (3) is that the score function V log is unknown, the distribution Pt not explicitly available. In practice, rather than using the exact score function v log Pr- approximate estimates for it are learned from the data by training neural networks on a score-matching objective (20 35 40) This objective is given by
minimize #EB
E s6
X,) V log Pt (X.))l
where {sp : 0 e {} is a sufficiently rich function class, such as that neural network. Substituting the learned score estimate s. into the backward process 33, obtain the following practical continuous-time backward SDE,
dxt = ] (x 2T. -t, dt dW.
(5)
Since this continuous backward SDE cannot be simulated exactly, is typically approximated using discretization methods.
Discretization Schemes. We outline the four discretization methods considered this work for solving the practical reverse SDE (5). Let h > 0 be the step size. Without loss of generality, assume T = Nh, where N a positive integer. For simplicity, we denote 1 X + (T X,) by (T - t, x;)),
), and define
1+h Ar == With - We AWW, = - / dw.
EULER-MARUYAMA SCHEME:
Given the step size ht the following approximation holds
Xtth = Xt + $1
T[ fnn((
vt, dt AnW, ( t, x,x A.W.
3

---

We derive the following discretized process for = ,,.... , 1:
JEM (1 + h/2)0 hs. nh, ( + v Vhgn.
where 55th pr and Sn NN Ia)-
EXPONENTIAL INTEGRATOR: Inspired by the work 18 47 propose a more refined discretiza-
tion method which solves the backward SDE (3) explicitly, yielding the following approximation
Xith = etxx + S/" e s.
s. - - X,,) A,W,
t xx X + 2(e (t, x() + Ar
We derive the following discretized process for n = ,,..... - 1;
0 e 0 + 2(e 1)s nh, ) + Veh 15,
with the initial point o0 ~ pr and En ~ Ia).
Unlike the Euler method, the midpoint randomization method evaluates the function >(T - X;=) a random point within the time interval [0. rather than at the start. Let U be random variable uniformly distributed in [0, 1) and independent of the Brownian motion W,. The randomized midpoint method exploits the approximation
VANILLA MIDPOINT RANDOMIZATION:
X X + / -
- t- X du + AnW,
(6)
x + hT --U, XUU A,W,
The idea behind the randomized midpoint method is to introduce an U in [0, 1], making hy(T (U,,(tt)) an estimator for integral f >(T - -t v, x de.
Furthermore, the intermediate term X is generated by employing the Euler method. We then derive the following discretized process for n = 0, 0...., 1:
Step 1
Generate S ~ NO, Ia) and Un Unif [0, 1). Set S, VUnsh + V1 Un
Step 2 With the initialization @REM ~ pr, define
REM REM
hU,>(T nh, RREM) + hhun.. n
REM n
+ hyy(T - (n + U, UREM + VhÊ..
EXPONENTIAL INTEGRATOR WITH MIDPOINT RANDOMIZATION: Combining midpoint
randomization with the exponential integrator approach, we propose the following new discretization process for .,...,
Step 1 Generate En El! NO, Ia) and Un Unif (0, 1]. Set &n Pu + VI el with
Pn 41446ll (1 (1-ce hhUi) [[chUn 1)(c -1)C= -1) 1)] 1/2
Step 2 With the initialization OREI pr, define
a+U hU /2,n/2 + 2 hUn/2 (T nh, ,REI) + VehUn. - 16 REI 4/2,RRI + he (1 3n)//2, (T - (n + Un)h, (RtT) ) + Veht - 16n
The resulting discrete process then solved generate new samples that approximately follow the data distribution Po.
3 Wasserstein Convergence Analysis under Various Discretization Schemes
In this section, we study the convergence diffusion models under EM, EI, REM and REI dis- cretizations of the continuous backward SDE 3. Specifically, we establish the upper bounds on

---

the Wasserstein-2 distance between the distribution of the N-th output of the SGMs under these discretization schemes and the target distribution
()(()))
E (EM, EI, REM, REI}.
Additionally, we analyze the number of iterations required for the Wasserstein distance to achieve pre-specified level under each discretization scheme. For clarity of presentation, omit constants in the main text, retaining only the key components that affect convergence rates. Full bounds with constants are provided in the proofs.
To establish the convergence analysis, we require the following assumption on target density Po- Assumption The target density po is mo-strongly log-concave, and the score function V log Po LopLischiiz.
Under Assumption 1 Pt is m({)-sttongly log-concave, V Po Pt is (())-Lipschitz. Moreover, m() is lower bounded by Munin min(1, ma), and T.()) is upper bounded by Tamax = + Lo, as summarized in Lemma and (see Appendix B We also assume that the score function V Pl (c (1)
exhibits linear growth in ||=]]] stated below.
N - 1, it holds that Assumption 2. There exists constant M1 > such that for n. = 1,.,...,
sup (n
VlogPr-s [VogPr-1<z)
< M1h(1 + |=]])), Vz.
The above condition a relaxation of the standard Lipschitz condition on the score function. More- require the following assumption on the score matching approximation at each point 0n. Assumption 3. Given a small E O, the score estimator satisfies
sup log - s. [ - nh, 11 & Esc
N
Assumption [ and 3 are standard in the Wasserstein convergence analysis of the score-based diffusion model. These assumptions were previously adopted in [12 13 and can be easily verified in the Gaussian case.
3.1 Euler-Maruyama Method and Exponential Integrator
In the following theorem, we quantify the Wasserstein distance between the distribution of &x E (EM, El) and the target distribution po-
Theorem 1. Suppose that Assumptions' 7 2 and 3 hold, it holds that
W2(C(U)) PO) S Momin Xolll C7 ddh + Ey Eacs
{EM, El),
(7)
Liuxx + 1/2
where
EM $
and E5 Mmii 1/2
1/2
Linax
6E
and 6 Mmii 1/2
Minin - 1/2
with mmin min(1, mb) and Lmax + Lo
Before comparing the above bound with existing ones, we state direct consequence.
Corollary 2. Given a small & > 0 and Esc = O the Wasserstein distance satisfies
W2(C(0%) ((0)).)) E, E (EM,El) after N o log (v)) iterations, provided that T O log( and = 0
We now explain the upper bound in Theorem which decomposes the total error into three parts: initialization, discretization, and score-matching errors. The first term on the right-hand side of display 7 bounds the error arising from initializing the reverse idealized continuous-time SDE (4) pr instead of pr The second term in captures the discretization error from the discretization scheme, while the third term reflects the score matching error.
5

---

The term -1/2 in Mmin - 1/2 arises from the drift term of the forward SDE 2 demonstrated in Lemma in the appendix. We refer interested readers the appendix for further details. More generally, if the forward SDE takes the form dX, = -3X, di + dB, then the bound becomes m(t) - 3, reflecting dependence on the coefficient 3. Although the current theory still assumes that the data distribution Po is strongly log-concave, the condition mg > 1/2 naturally generalizes to mo > B for any 3 > 0.
The convergence rates of EM method obtained in Theorem 1 and Corollary 2 align with Theorem 2 and Proposition in [3]. Moreover, we note that the convergence rate of EM and El schemes are comparable, which is consistent with the error bounds for these two schemes in KL divergence established in Theorem 1 of [5], where N (2(T22 in their setting.
3.2 Randomized Midpoint Method
Since the randomized midpoint method involves the i.i.d. uniformly distributed random variable Un. the resulting score matching function ST 1, can be evaluated any E [0, 77. To proceed, we impose a regularity condition on the deviation of the estimated score from the true score at these points. For this, we introduce the following auxiliary stochastic processes. Given Un u, define the conditional realization of the random vector GREM via
n-+U7
GREM = (1
wh) REM + uhs, (T - nh, (REM) + Vuhgn.
Similarly, we define the following conditional realization of @REI n+U
G == 22/22,20 RREI + 2(e 1)s, (T - nh, RRI)
+ V Veuh 1%
We impose the following assumption score estimates, which extends Assumption Assumption 4. There exists constant E 0 such that for any u [0, 1] and O,
V log Pr (r-(n))) (U ST - (n + u/h, o Wlz, & Escr a a (REM, REI].
In what follows, we provide the upper bound for the Wasserstein-2 distance between the law of &% E {REM, REIJ and the data distribution Pu-
Theorem 3. Suppose that Assumptions 1 and 4 hold, then for a E (REM, REI),
W2(C(O%) Po) S Mainn |olll (d) 62 Eec
where @REM(()) (d) = d/3Lmax + +1/2/3
3
and gREM
REM
Mymin 1/2
mmmin 1/2
V d/3Lmaxx and '
3
REI (d)
(mmin - 1/2)
Mmin 1/2
with Temax and Mmin as defined in Theorem
The convergence rate of these two schemes is generally consistent with EM, El methods, differing only the coefficients. As shown in display (6). the key idea behind the randomized midpoint method is to introduce a uniformly distributed random variable Un to evaluate the term (T t. x,) at random point within the time interval [0, h). In the proofs provided Appendix B.3 and B.4 (corresponding to a = REM and a = REI. respectively), we demonstrate that Ev, p n+ yields an accurate approximation to the randomization also introduces variance term ||00+ - Eu, [0% 1ll1,, of order ((33/2) ((h)sse. which obscures the benefits of the improved estimation. As a result, the midpoint randomization offers improvement in convergence over EM and El methods.
distribution, with error of order (15/22) (h) However, this
Although midpoint randomization itself does not improve convergence our setting, enables paral- lel computation [ 23), significantly reducing computational complexity and enhancing efficiency.

---

4 Second-order Acceleration
In this section, we propose accelerated sampler that leverages Hessian estimation. The idea behind the acceleration is the Local Linearization Method, introduced in [33]. which approximates the drift term of an SDE using its Itô expansion over small time intervals. To illustrate this, we begin with general framework for the backward process.
dr y(T - 1, + a dW.,
(8)
where Ø > 0 and W2 is the d-dimensional Brownian motion. We assume that "(t, 2) E C1,,.R xR*) and approximate it by a linear function in both state and time within each discretization step. Applying Itö's formula to ((T 1, 2),, we derive the following approximation for 7(T' (x)))) - r(T [ Ts)
0² 83 (T- --- s. at 2 a.
0% (T (rs-sx.) s.z,))
(t s)=s)=r + Oy (T - s. (xt - -ee)
dz
Here and henceforth, we abbreviate the partial derivative express 7(T - x.) in the following form
8 g(s) 02"
as
0g (zo). This allows us to 02
T(T tz) = T(T - -5, + Le(exe - 2x,) Mo(t 8)
(9)
with
L, 07 -T Dr
and M, 0² 02 Or. (T 8, 2.,)
07 (T- x.). at
Here, Lg is the first-order spatial derivative of 7, capturing its local variation with respect to position F. The term M, represents the temporal evolution of y, incorporating information about how its shape changes over both space and time. Substituting display ( into the original SDE (8). we obtain
dx, = (T 8, x.) + [,, xx.) + M, - s)] di + + dW,
This formulation ensures that the discretized process preserves the essential structure of the original dynamics while remaining computationally tractable. Unlike discretization schemes in Section 2 which rely on direct numerical integration, this transformed SDE allows for analytical solutions within each small time interval.
Setting t, = x/2 v log o = 1, let At [0, h]] and nh. By Itō's formula, we derive the analytical expression for z₁. In the resulting expression, we denote by 0so. Then, we obtain that for any L [nh, (n + 1)h), it holds that
1" nh
( gsso + V log PogPr-nhh mh (JS)) + To (Tu (05) + Ma(u Ma(u-nhh)
nhl) du + f dW. nb (10)
T2 = 050 +
2
where
d In マ 2 logPr-nn (050), Mn = Σ j=1
82 Ox3
0 Pr (05))
V MogPr-nhi (15
a
Thus, Ln contains the Hessian of the score function. The term M,, measures the difference between the spatial and temporal changes in the score function, reflecting the balance between curvature effects and temporal adaptation in the diffusion process. Notice that ((++1) is the point we aim to approximate. For this. we need estimate the score function and its higher-order derivatives to obtain accurate estimates of Ln, and Mn, denoted by (L) and s (M) respectively.
By the work of [28] higher-order derivatives of log px( with respect to can be accurately estimated. Moreover, we show in Appendix C that a V log (x can be expressed as a combination of up to second-orde partial derivatives of log Pr (z) with respect to æ. Thus, estimating O V (o))(()) requires no additional assumptions beyond those for accurately estimating 02V log Pz(x). We also require the following assumption.
Assumption 5. For some constants score function satisfies that
(L)) ,6 > 0, the estimate for high-order derivatives of the
(L) 80 (T (SD) - I.n
sup DE -1
< $0(0)
1
(M) (T- 050) Mn
sup -
11
7

---

This assumption has been adopted in prior works 24 and 27 Substituting these estimates into display (TO) then gives
-550
[T(T nh,, 50) + 80 (T - -nh, 100) ) - -00)
(15)
3t
nh
8x (M)) (T nh, ) nh)) du + f" dWw.
nh
Let ust1 denote ((n+1)h- The second discretization scheme is given by
x+1 -550 (T) (T nh, 0 ( (T 1 Ls) 5 + s+ sT nh, 1200
+ sx
(T- (T
nh, -000 (T-h, 5: (T ( nh, $50h Is) s(M)(T nh, 50
(a+1)h 1)h
(T ) M 1)h 1jh dw.
nhh We assume additional smoothness on the score function.
Assumption 6. Let || |r denote the Frobenius norm. There exists positive constant LF such that
V2 log V log < LF ||
yl1
Vz, y E RR
As shown in Theorems 4 and [27]. this condition plays crucial role in bounding the Wasserstein distance for Hessian estimates and can be easily verified in the Gaussian case.
Assumption 7. There exists constant M2 0 such that, for any .....,N N 1 and [nh, (n + 1)h,, it holds that
V log pr-(())( 2 log (x)] < M2h(1 + [0=I),, Vo E Rd
We now quantify the W2 distance between the generated distribution ((05) and the target po- Theorem 4. Suppose that Assumptions
and hold, then
W2(C(U5))) ),P))
Mimin [ + (g)((d)) 55
(5ar (Emtt 3
2 Vide!L) 2 + he
1 (11)
where (050(d (d) el 1/2)h. vd(Lin2 3/2 + v2Lr/4)
clhma 11/23h with Lmax and
mmin - 1/2 and 6500
mmin - 1/2
Miin as defined in Theorem
Before diseussing the results of this theorem, let us state its direct consequence.
Corollary 5. For a given E the Wasserstein distance satisfies W2(C(d po) < after
= 0
Ld log Ld
iterations, provided that T O (log( (oog(()) and h 0 =))
In the above, we present the first Wasserstein convergence analysis an accelerated sampler that utilizes accurate score function estimation and second-order information about log-densities. The total error arises from the same resources as in Theorem T The first term on the right-hand side of display (TT) captures the initialization error. The second term reflects the benefit of second-order acceleration, where improved discretization reduces the error in approximating the reverse SDE. The third term accounts for errors in estimating both the spatial and temporal components of the score function, as well as the higher order terms L,, and M,.
The accelerated convergence of this method is driven by two key innovations: approximating the drift term through its Itô expansion rather than endpoint evaluations, and deriving a closed-fom solution to the integral equation (IO) using the Itô formula, akin to Exponential Integrator techniques. Compared to the four schemes described in Section 3 the proposed second-order algorithm offers a clear computational ad advantage: it requires only O(1/s) iterations (v.s. O(1/=²)), and permits
Due space limitations, refer readers Appendix for complete derivation.

---

EM
EI
REM
REI
S0
2
3
3
3
4
4
4
/
4
5
5
5
5
S
6
6
6
12
7
?
9
7
1.0
8
8
8
7
9
9
9
Sisp Sies
(a)
(b)
Figure 1: (a): Errors of SGMs under EM, EI, REM, REI and SO with different choices of h. (b);
Samples generated by five different algorithms on the MNIST dataset.
larger step size (e)) enabling faster progression in each iteration instead of OE As result, the second-order method achieves the same accuracy with fewer iterations and a larger step size, making a more efficient method for approximating the target distribution.
We note that the convergence rate ((1/E) matches that of the accelerated Denoising Diffusion Proba- bilistic Models (DDPM) sampler [17] proposed in [27], which achieves this rate in KL. divergence and relies on Hessian estimation. Additionally, this rate aligns with the iteration complexity results in TV distance presented in [99 when setting 1 in their framework.
5 Numerical Studies
In this section, we compare the performance of the SGMs under EM, EI, REM, REI discretization schemes described in Section as well as the second-order acceleration method (SO) proposed in Section 4 evaluated on both synthetic data and the MNIST dataset.
5.1 Experiments on Synthetic Data
We apply the five algorithms to the posterior of penalized logistic regression, defined by (o) oc ex)(-f(0)), with the potential function f defined via
|1012 f(0) = 2
1
Σ 1=1
log (1 exp(-v<x, 0)),
Tidata
where > 0 denotes the tuning parameter. The data (1i, yr Mer, composed of binary labels
itd. (O, 100). Set y: E {-1.1} and features 2; = (I..1. .............). :T.) E Idd generated from æi,
= 100 with d 20 and Pldata 100. Appendix F provides additional results and implementation details. Figure Ta presents the W2 distance measured along the first dimension between the empirical distributions of the samples from the five algorithms and the target distribution under different choices of h. These results support our theoretical findings: all discretization schemes from Section exhibbi similar convergence behavior, while the proposed Hessian-based sampler consistently achieves superior performance.
5.2 Real Data Analysis
apply the four SGM discretization schemes and the second-order algorithm to the MNIST dataset. To accelerate the SO algorithm, we use Hessian-vector products (HVPs) instead of explicitly computing the Hessian. The results in Figure Tb demonstrate that SO outperforms the others. Additional results and implementation details provided in Appendix G

---

Discussion
The W distance, used this paper, serves natural and practical metric for measuring emrors in diffusion models. However, recent work on the convergence theory of diffusion models has also explored alternative metrics such as total variation distance and KL divergence. A promising direction for future research is to establish convergence guarantees with respect to these alternative distances.
Moreover, while this work makes progress in provably accelerating SDE-based diffusion sampling in Wasserstein distance, it would also be valuable explore deterministic samplers based on probability flow ODEs.
Finally, for clarity and simplicity, we foeus specific choice of drift functions and g in the forward process, corresponding to the more general framework with broader choices of g)) is an interesting avenue for future research.
process. Extending this analysis to a
Acknowledgments and Disclosure of Funding
Lu Yu is supported by the City University of Hong Kong Startup Fund and Hong Kong RGC Grant 21306325.
References
[1] Namrata Anand and Tudor Achim. Protein structure and sequence generation with equivariant denoising diffusion probabilistic models. arXiv preprint arXiv:2205.15011 2022.
[2] Brian DO Anderson. Reverse-ting time diffusion equation models. Stochastic Processes and their Applications, 12(3):313-326 1982.
[3] Stefano Bruno, Ying Zhang, Dong-Young Young Lim, Ömer Deniz Akyildiz, and Sotirios Sabanis. On diffusion-based generative models and their error bounds: The log-concave case with full convergence estimates. arXiv preprint arXiv:2311.3353
2023.
[4] Patrick Cattiaux, Giovanni Conforti, Ivan Gentil, and Christian Léonard. Time reversal of diffusion processes under a finite entropy condition. In Annales de I' l'Institut Henri Poincaré (B) Probabilités et Statistiques, volume 59, pages 1844-1881. Institut Henri Poincaré, 2023.
[5] Hongrui Chen, Holden Lee, and Jianfeng Lu. Improved analysis of score-based generative modeling: User-friendly bounds under minimal smoothness assumptions. In International Conference on Machine Learning, pages 4735-4763. PMLR, 2023.
[6] Minshuo Chen, Song Mei, Jianqing Fan, and Mengdi Wang. An overview of diffusion models: Applications, guided generation, statistical rates and optimization. arXiv preprint arXiv:2404.07771 2024.
[7] Sitan Chen, Sinho Chewi, Jerry Li, Yuanzhi Li, Adil Salim, and Anru R Zhang. Sampling is as easy learning the theory for diffusion models with minimal data assumptions. arXiv preprint arXiv:2209. 1215, 2022.
[8] Sitan Chen, Sinho Chewi, Holden Lee, Yuanzhi Li, Jianfeng Lu, and Adil Salim. The probability flow ode is provably fast. Advances in Neural Information Processing Systems, 36, 2024.
[9] Xiang Cheng, Niladri S Chatterji, Peter L Bartlett, and Michael I Jordan. Underdamped langevin meme: A non-asymptotic analysis. In Conference on learning theory, pages 300-323. PMLR, 2018.
[10] Valentin De Bortoli. Convergence of denoising diffusion models under the manifold hypothesis. arXiv preprint arXiv:2208.05314, 2022.
[11] Prafulla Dhariwal and Alexander Nichol. Diffusion models beat gans on image synthesis.
Advances in neural information processing systems, 34:8780-8794, 2021.
10

---

[12] Xuefeng Gao and Lingjiong Zhu. Convergence analysis for general probability flow odes of diffusion models in wasserstein distances. arXiv preprint arXiv:2401.17956 2024.
[13] Xuefeng Gao, Hoang M Nguyen, and Lingjiong Zhu. Wasserstein convergence guarantees for general class of score-based generative models. arXiv preprint arXiv:2311.1100) 2023.
[14] Marta Gentiloni-Silee and Antonio Ocello. Beyond log-concavity and score regularity: Improved convergence bounds for score-based generative models in w2-distance. arXiv preprint arXiv:2501.02298, 2025.
[15] Shivam Gupta, Linda Cai, and Sitan Chen. Faster diffusion- baasd sampling with randomized midpoints: Sequential and parallel. arXiv preprint arXiv:2406.00924, 2024.
[16] 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, 33:7366-7376, 2020.
[17] Jonathan Ho, Ajay Jain, and Picter Abbeel. Denoising diffusion probabilistic models. Advances in neural information processing systems, 33:6840-6851, 2020.
[18] Marlis Hochbruck and Alexander Ostermann. Exponential integrators. Acta Numerica, 19: 209-288 209-286,22 2010.
[19] Daniel Zhengyu Huang, Jiaoyang Huang, and Zhengjiang Lin. Convergence analysis of probability flow ode for score-based generative models. arXiv preprint arXiv:2444.09730, 2024.
[20] Aapo Hyvärinen and Peter Dayan. Estimation of non-normalized statistical models by score matching. Journal of Machine Learning Research, 6(4), 2005.
[21] Saravanan Kandasamy and Dheeraj Nagaraj. The poisson midpoint method for langevin dynam- ies: Provably efficient discretization for diffusion models. arXiv preprint 24Xiv:2405. /7068, 2024.
[22] Gen Li and Changxiao Cai. Provable acceleration for diffusion models under minimal assump- tions, arXiv preprint arXiv:2410.23285,, 2024.
[23] Gen Li and Yuchen Jiao. Improved convergence rate for diffusion probabilistic models. arXiv preprint 1331:22410.13388, 2024.
[24] Gen Li, Yu Huang, Timofey Efimov, Yuting Wei, Yuejie Chi, and Yuxin Chen. Accelerating convergence of score-based diffusion models, provably. arXiv preprint arXiv:2403.03852, 2024.
[25] Gen Li, Yuting Wei, Yuejie Chi, and Yuxin Chen. A sharp convergence theory for the probability flow odes of diffusion models. arXiv preprint arXiv:2408.02320, 2024.
[26] Jiadong Liang, Zhihan Huang, and Yuxin Chen. Low-dimensional adaptation of diffusion models: Convergence in total variation. arXiv preprint arXiv:2501.12982. 2025.
[27] Liang, Ju, Y Liang, and N Shroff. Broadening target distributions for accelerated diffusion models via a novel analysis approach. arXiv preprint arXiv:2402.13900 2024.
[28] Chenlin Meng, Yang Song, Wenzhe Li, and Stefano Ermon. Estimating high order gradients of the data distribution by denoising. Advances in Neural Information Processing Systems, 34: 25359-253699 2021.
[29] Nikiforos Mimikos- Sttiioos-ttmmttoouooosss Benjamin Zhang, and Markos Katsoulakis. Score-based generative models are provably robust: an uncertainty quantification perspective. arXiv preprint arXiv:2405. 15754, 2024.
[30] Vadim Popov, Ivan Vovk, Vladimir Gogoryan, Tasnima Sadekova, and Mikhail Kudinov. Grad- tts: diffusion probabilistic model for text-to-speech. International Conference on Machine Learning, pages 8599 8608. PMLR, 2021.

---

[31] Aditya Ramesh, Prafulla Dhariwal, Alex Nichol, Casey Chu, and Mark Chen. Hierarchical text-conditional image generation with clip latents. arXiv preprint 07X1522004.06/25, 1(2):3 2022.
[32] Ruoqi Shen and Yin Tat Lee. The randomized midpoint method for log-concave sampling. Advances in Neural Information Processing Systems, 32, 2019.
[33] Isao Shoji. Approximation of continuous time stochastic processes by a local linearization method. Math. Comput., 67:287-298, 1998. URL http://////pp.sssmttiisschhiiiie org/
Corpus ID: 13430357
[34] Jiaming Song, Chenlin Meng, and Stefano Ermon. Denoising diffusion implicit models. arXiv preprint arXiv:2010.02550 02502, 2020.
[35] Yang Song, Sahaj Garg, Jiaxin Shi, and Stefano Ermon. Sliced score matching: scalable approach to density and score estimation. In Uncertainty in Artificial Intelligence, pages 574-584 584, PMLR, 2020.
[36] Yang Song, Jascha Sohl-Dickstein, Diederik P Kingma, Abhishek Kumar, Stefano Ermon, and Ben Poole. Score-based generative modeling through stochastic differential equations. arXiv preprint arXiv:2011.13456 2020.
[37] Stanislas Strasman, Antonio Ocello, Claire Boyer, Sylvain Le Corff, and Vincent Lemaire. An analysis of the noise schedule for score-based generative models. arXiv preprini arXiv:2402.046500 2024.
[38] Wenpin Tang and Hanyang Zhao. Contractive diffusion probabilistic models. arXiv preprint orXiv:2401.13115 2024.
[39] Wenpin Tang and Hanyang Zhao. Score-based diffusion models via stochastic differential equations-a technical tutorial. arXiv preprint arXiv:2402.07487. 2024.
[40] Pascal Vincent. A connection between score matching and denoising autoencoders. Neural computation, 3377):661-16774 2011.
[41] Yuchen Wu, Yuxin Chen, and Yuting Wei. Stochastic runge-kutta methods: Provable accelera- tion of diffusion models. arXiv preprint arXiv:2410.0476/ 2024.
[42] Minkai Xu, Lantao Yu, Yang Song, Chence Shi, Stefano Ermon, and Jian Tang. Geodiff: geo- metric diffusion model for molecular conformation seneraaatiiggnooiiooonn arXiv preprint arXiv:2203.02923. 2022.
[43] Ling Yang, Zhilong Zhang, Yang Song. Shenda Hong, Runsheng Xu, Yue Zhao, Wentao Zhang, Bin Cui, and Ming-Hsuan Yang. Diffusion models: A comprehensive survey of methods and applications. ACM Computing Surveys, 56(4):1-39, 2023.
[44] Lu Yu and Arnak Dalalyan. Parallelized midpoint randomization for langevin monte carlo. arXiv preprint arXiv:2402. 14434, 2024.
[45] Lu Yu, Avetik Karagulyan, and Arnak Dalalyan. Langevin monte carlo for strongly log-coneave distributions: Randomized midpoint revisited. arXiv preprint arXiv:2306.08494, 2023.
[46] Yifeng Yu and Lu Yu. Randomized midpoint method for log- log-concave sampling under constraints. arXiv preprint arXiv:2405. /5379, 2024.
[47] Qinsheng Zhang and Yongxin Chen. Fast sampling of diffusion models with exponential integrator. arXiv preprint 23X4v:2204.13002, 2022.
12

---

NeurIPS Paper Checklist
1. Claims
Question: Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope?
Answer: [Yes]
The main claims in the abstract and introduction are substantiated by the detailed theoretical analyses presented in the paper. These sections clearly outline our contributions and scope, which are consistently reflected throughout the findings.
Justification:
Guidelines:
The answer NA means that the abstract and introduction do not include the claims made in the paper.
The abstract and/or introduction should clearly state the claims made, including the contributions made the paper and important assumptions and limitations. A No or
NA answer to this question will not be perceived well by the reviewers.
The claims made should match theoretical and experimental results, and reflect how much the results can be expected to generalize to other settings.
It is fine to include aspirational goals as motivation long is clear that these goals are not attained by the paper.
2. Limitations
Question: Does the paper discuss the limitations the work performed by the authors?
Answer: [Yes]
Justification: In the last section, some strong assumptions made in this work are highlighted, along with potential directions for relaxing these assumptions and generalizing the current findings.
Guidelines:
The answer NA means that the paper has limitation while the answer No means that the paper has limitations, but those are not discussed in the paper.
The authors are encouraged to create separate "Limitations" section in their paper.
The paper should point out any strong assumptions and how robust the results are to violations of these assumptions (e.g.. independence assumptions, noiseless settings, model well-specification, asymptotic approximations only holding locally). The authors should reflect on how these assumptions might be violated in practice and what the implications would be.
The authors should reflect on the scope of the claims made, e.g., if the approach was only tested few datasets or with few runs. In general, empirical results often depend on implicit assumptions, which should articulated.
The authors should reflect on the factors that influence the performance of the approach. For example, facial recognition algorithm may perform poorly when image resolution is low or images are taken in low lighting. Or speech-to-text system might not be used reliably provide closed captions for online lectures because it fails handle technical jargon.
The authors should discuss the computational efficiency of the proposed algorithms and how they scale with dataset size.
If applicable, the authors should discuss possible limitations of their approach to address problems of privacy and fairesss.
While the authors might fear that complete honesty about limitations might be used by reviewers as grounds for rejection, worse outcome might be that reviewers discover limitations that aren't acknowledged
the paper. The authors should use their best judgment and recognize that individual actions in favor of transparency play impor- tant role in developing norms that preserve the integrity of the community. Reviewers will be specifically instructed to not penalize honesty concerning limitations.
3. Theory assumptions and proofs
13

---

Question: For each theoretical result, does the paper provide the full set of assumptions and complete (and correct) proof?
Answer: [Yes]
Justification: For the theoretical results, we provide the full set of assumptions and complete proofs for them.
Guidelines:
The answer NA means that the paper does include theoretical results. All the theorems, formulas, and proofs in the paper should be numbered and cross- referenced.
All assumptions should be clearly stated or referenced in the statement of any theorems.
The proofs can either appear in the main paper or the supplemental material, but if they appear in the supplemental material, the authors are encouraged provide short proof sketch to provide intuition.
Inversely, any informal proof provided in the core of the paper should be complemented by formal proofs provided in appendix or supplemental material.
Theorems and Lemmas that the proof relies upon should be properly referenced.
4. Experimental result reproducibility
Does the paper fully disclose all information needed to reproduce the main ex- perimental results of the paper to the extent that it affects the main claims and/or conclusions of the paper (regardless whether the code and data are provided or not)?
Question:
Answer: [Yes]
Justification: All information necessary to reproduce the main experimental results is provided in the main text and appendix.
Guidelines:
The answer NA means that the paper does not include experiments.
the paper includes experiments, a No answer to this question will not perceived well by the reviewers: Making the paper reproducible is important, regardless of whether the code and data are provided or not.
If the contribution is dataset and/or model, the authors should describe the steps taken to make their results reproducible or verifiable.
Depending on the contribution, reproducibility can be accomplished in various ways. For example, if the contribution a novel architecture, describing the architecture fully might suffice, or if the contribution is specific model and empirical evaluation, it may be necessary to either make it possible for others to replicate the model with the same dataset, provide access to the model. In general. releasing code and data is often one good way to accomplish this, but reproducibility can also be provided via detailed instructions for how replicate the results, access hosted model (e.g.. in the of a large language model), releasing of a model checkpoint, or other means that are appropriate to the research performed.
While NeurIPS does not require releasing code, the conference does require all submis- sions to provide some reasonable avenue for reproducibility, which may depend on the nature of the contribution. For example
(a) If the contribution is primarily new algorithm, the paper should make it clear how to reproduce that algorithm.
(b) the contribution is primarily a new model architecture, the paper should describe the architecture clearly and fully.
(c) the contribution is a new model (e.g., a large language model), then there should either be way access this model for reproducing the results way to reproduce the model (e.g.. with an open-source dataset or instructions for how to construct the dataset).
(d) We recognize that reproducibility may be tricky in some cases, in which case authors are welcome describe the particular way they provide for reproducibility. In the case of closed-source models, it may be that access to the model is limited in some way (c.g., to registered users), but it should possible for other researchers to have some path to reproducing verifying the results.
14

---

5. Open access to data and code
Question: Does the paper provide open access to the data and code, with sufficient instruc- tions to faithfully reproduce the main experimental results, described in supplemental material?
Answer: [Yes]
Justification: The datasets used in this work either synthetic publicly available, and the code for reproducing the results provided in the supplementary materials.
Guidelines:
The answer NA means that paper does not include experiments requiring code.
Please see the NeurlPS code and data submission guidelines (https://nips.com
for more details.
While we encourage the release of code and data, we understand that this might not be possible, "No" an acceptable answer. Papers cannot be rejected simply for including code, unless this is central to the contribution (e.g., for new open-source benchmark).
The instructions should contain the exact command and environment needed to run to reproduce the results. See the NeurlPS code and data submission guidelines (https
CodeSuomi ssininoiscsyyyy for more details.
The authors should provide instructions on data access and preparation, including how to access the raw data, preprocessed data, intermediate data, and generated data, etc.
The authors should provide scripts to reproduce all experimental results for the new proposed method and baselines. only subset of experiments are reproducible, they should state which ones omitted from the script and why.
At submission time, preserve anonymity, the authors should release anonymized versions (if applicable).
Providing as much information as possible in supplemental material (appended to the paper) recommended, but including URLs to data and code permitted.
6. Experimental setting/details
Does the paper specify all the training and test details (c.g., data splits, hyper- parameters, how they were chosen, type of optimizer, etc.) necessary understand the results?
Question:
Answer: [Yes]
Justification: Experimental settings and details are provided.
Guidelines:
The answer NA means that the paper does not include experiments.
The experimental setting should be presented in the core of the paper to a level of detail that is necessary appreciate the results and make sense of them.
The full details can be provided either with the code, appendix, or as supplemental material.
7. Experiment statistical significance
Question: Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?
Answer: [No]
Justification: Statistical significance the experiments is not applicable in this context.
Guidelines:
The answer NA means that the paper does not include experiments.
The authors should answer "Yes" if the results are accompanied by error bars, confi- dence intervals, or statistical significance tests, least for the experiments that support the main claims of the paper.
The factors of variability that the error bars are capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, or overall with given experimental conditions).
15

---

The method for calculating the error bars should be explained (closed form formula, call a library function, bootstrap, etc.)
The assumptions made should be given (e.g., Normally distributed errors).
It should be clear whether the error bar is the standard deviation the standard error of the mean.
It is OK to report 1-sigma error bars, but one should state it. The authors should preferably report a 2-sigma error bar than state that they have 96% CI, if the hypothesis of Normality of errors is not verified.
For asymmetric distributions, the authors should be careful figures symmetric error bars that would yield results that are out of range (e.g. negative rates).
show in tables or
If error bars are reported in tables or plots, The authors should explain the text how they were calculated and reference the corresponding figures tables in the text.
8. Experiments compute resources
For each experiment, does the paper provide sufficient information on the puter resources (type of compute workers, memory, time of execution) needed to reproduce the experiments?
Question:
Answer: [Yes]
Justification: Sufficient information about the computational resources used is provided in appendix.
Guidelines:
The answer NA means that the paper does not include experiments.
The paper should indicate the type of compute workers CPU or GPU, internal cluster, cloud provider, including relevant memory and storage.
The paper should provide the amount of compute required for each of the individual experimental runs as well as estimate the total compute.
The paper should disclose whether the full research project required more compute than the experiments reported in the paper (e.g., preliminary or failed experiments that didn't make into the paper).
9. Code of ethics
Question: Does the research conducted the paper conform, in every respect, with the NeurlPS Code of Ethics https::/neurrig:: httet
Answer: [Yes]
Justification: This work fully complies with the NeurIPS Code of Ethics.
Guidelines:
The answer NA means that the authors have reviewed the NeurIPS Code of Ethics.
If the authors answer No, they should explain the special circumstances that require deviation from the Code of Ethics.
The authors should make sure to preserve anonymity (e.g., there a special consid- eration due to laws or regulations in their jurisdiction).
10. Broader impacts
Question: Does the paper discuss both potential positive socictal impacts and negative societal impacts of the work performed?
Answer: [NA]
Justification: This is theoretical work without addressing or having any societal impact.
Guidelines:
The answer NA means that there is no societal impact of the work performed.
If the authors answer NA or No, they should explain why their work has no societal impact why the paper does not address societal impact.
16

---

Examples of negative societal impacts include potential malicious unintended uses (e.g., disinformation, generating fake profiles, surveillance), faireess considerations (e.g., deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations.
The conference expects that many papers will be foundational research and not tied to particular applications, let alone deployments. However, if there is direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster.
The authors should consider possible harms that could arise when the technology is being used intended and functioning correctly, harms that could arise when the technology is being used as intended but gives incorrect results, and harms following from (intentional unintentional) misuse of the technology.
If there are negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition to attacks, mechanisms for monitoring misuse, mechanisms monitor how system learns from feedback over time, improving the efficiency and accessibility of ML).
11. Safeguards
Does the paper deseribe safeguards that have been put in place for responsible release of data or models that have a high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?
Question:
Answer: [NA]
Justification:
This work poses no such risks.
Guidelines:
The answer NA means that the paper poses no such risks.
Released models that have high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere to usage guidelines or restrictions to access the model or implementing safety filters.
Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images.
We recognize that providing effective safeguards is challenging, and many papers do not require this, but we encourage authors to take this into account and make best faith effort.
12. Licenses for existing assets
Are the creators or original owners of assets (e.g., code, data, models), used in the paper, properly credited and are the license and terms of use explicitly mentioned and properly respected?
Question:
Answer: [NA]
Justification: This work does not use existing assets.
Guidelines:
The answer NA means that the paper does use existing assets. The authors should cite the original paper that produced the code package dataset.
The authors should state which version of the asset used and, possible, include URL.
The name of the license (e.g., CC-BY 4.0) should be included for each asset.
For scraped data from a particular source (e.g., website). the copyright and terms of service of that source should be provided.
If assets are released, the license, copyright information, and terms use in package should be provided. For popular datasets, paperswithcode com/ datasets has curated licenses for some datasets. Their licensing guide can help determine the license of dataset.
17

---

For existing datasets that are re-packaged, both the original license and the license of the derived asset (if has changed) should be provided.
If this information is not available online, the authors are encouraged to reach out to the asset's s creators.
13. New assets
Question: Are new assets introduced in the paper well documented and is the documentation provided alongside the assets?
Answer: [NA]
Justification:
This work does not release new assets.
Guidelines:
The answer NA means that the paper does not release new assets.
Researchers should communicate the details of the dataset/code/moddl as part of their submissions via structured templates. This includes details about training, license, limitations, etc.
The paper should discuss whether and how consent was obtained from people whose asset used.
At submission time, remember to anonymize your assets (if applicable). You can either create anonymized URL or include anonymized zip file.
14. Crowdsourcing and research with human subjects
Question: For crowdsourcing experiments and research with human subjects, does the paper include the full text of instructions given to participants and screenshots, if applicable, as well as details about compensation (if any)?
Answer: [NA]
Justification: This work does not involve crowdsourcing research with human subjects.
Guidelines:
The answer NA means that the paper does not involve crowdsourcing nor research with human subjects.
Including this information in the supplemental material is fine, but if the main contribu- tion of the paper involves human subjects, then as much detail as possible should be included in the main paper.
According the NeurIPS Code of Bthies, workers involved in data collection, curation, or other labor should be paid at least the minimum wage in the country of the data collector.
15. Institutional review board (IRB) approvals or equivalent for research with human subjects
Question: Does the paper describe potential risks incurred by study participants, whether such risks were disclosed to the subjects, and whether Institutional Review Board (IRB) approvals (or equivalent approval/review based the requirements of your country or institution) were obtained?
Answer: [NA]
Justification: This work does not involve crowdsourcing nor research with human subjects.
Guidelines:
The answer NA means that the paper does not involve crowdsourcing nor research with human subjects.
Depending on the country in which research is conducted, IRB approval (or equivalent) may be required for any human subjects research. If you obtained IRB approval, you should clearly state this in the paper.
We recognize that the procedures for this may vary significantly between institutions and locations, and we expect authors to adhere to the NeurlPS Code of Ethics and the guidelines for their institution.
For initial submissions, not include any information that would break anonymity (if applicable), such the institution conducting the review.
18

---

16. Declaration of LLM usage
Question: Does the paper describe the usage of LLMs if it is an important, original, or non-standard component of the core methods in this research? Note that if the LLM is used only for writing, editing, or formatting purposes and does not impact the core methodology. scientific rigorousness, originality of the research, declaration is not required.
Answer: [NA]
Justification: This work does not involve LLMs as any important, original, or non- non-standard components.
Guidelines:
The answer NA means that the core method development in this research does not involve LLMs as any important, original, or non-standard components.
Please refer to our LLM policy ftteps://neuriie /t/CCooffeeencss//2222 for what should or should be described.
19

---

A Discussion of Theoretical Advances in Accelerating Samplers for Diffusion Models
Score-based diffusion models can be formulated using either SDEs or their deterministic counterparts, known as probability flow ODEs [3).. While SDE-based samplers generate samples through stochas- tic simulation, ODE-based samplers provide deterministic alternative. Theoretical advancements in accelerating these samplers have emerged only recently. A significant step toward designing provably accelerated, training-fie free methods made by [24]. who propose and analyze acceleration for both ODE- and SDE-based samplers. Their accelerated SDE sampler leverages higher-order expansions of the conditional density to enhance efficiency. This was followed by the work of [25], which provided convergence guarantees for probability flow ODEs. Furthermore, [19] studies the conver- gence properties of deterministic samplers based on probability flow ODEs, using the Runge- Kutta integrator: [1] propose and analyze training-free free acceleration algorithm for SDE-based samplers, based on the stochastic Runge- method. [27] proposes a novel accelerated SDE-based sampler when Hessian information is available. Another line of work involves the midpoint randomized method. In particular,
15] explore ODE acceleration by incorporating a randomized midpoint method, leveraging its advantages in parallel computation. A more recent work by [33 improved upon the ODE sampler proposed by [5], achieving the state-of-the-ar convergence rate.
We note that all of these works provide convergence analysis in terms of either KL divergence or TV distance. Among these, [27] accelerates the stochastic DDPM sampler by leveraging precise score and Hessian estimations of the log density, even for possibly non-smooth target distributions. This achieved through a novel Bayesian approach based on tilting factor representation and Tweedie's formula. [19] accelerates the ODE sampler by utilizing p-th (p > 1) order information of the score and employing early stopping. These two works are the most similar to our proposed accelerated sampler in that they all rely on the Hessian information of the log density. However, their settings differ from ours, and their convergence analyses are neither directly applicable to our framework nor precisely expressed in terms of Wasserstein distance
function, with the target distribution supported on compact
B Proof of Section 3
We define several stochastic processes associated with the backward process X and the sample path O,, First, recall that X;- is described by the following SDE:
dxt= = (5xt + (Voggp((x)) + dW,, X ~ ~TT,
and {0%, 0 & & N) satisfies the iterative law:
d n+ 9% (Wihnsss((+)))))
where a E (EM, EI, REM, REI, SO].
Based on X7_, we define the following two processes, (Y,, 0 < t & T} and (Ýr, < & h). Y, satisfies the SDE
dY: = (2y + Vlogpr-- dt + dW, Yo ~ pr.
Y, actually relies on X on the time interval [nh, (n + 1)4] for each n. However, we only need this notation in the proof of one-step discretization error, then we allow for some slight abuse of notation by omitting since it will not lead to any confusion. Therefore, (Y,,0 0 < l < h} satisfies
dÝ, (2r++ Vlogpr- dt + dw, Yo D,.
(12)
'When the target distribution is compactly supported, Pinsker's 's inequality allows translating TV or KL divergence into Wasserstein distance. However, this often yields loose bounds, especially in high dimensions. where the actual Wasserstein distance may be much smaller.
20

---

Recall that we have defined two processes nna and Pet in Section 3.2 that is, for any a E [0, and 0, N
(1+ uh R + uTs nh, (REM) + v Vuh5.:
I REM
1)s, @ 1 -1) nh, (RE) + Veuh 16
This section is devoted to proving the convergence rate of the diffusion model under various dis- cretization schemes. To this end, we need the following auxiliary lemmas.
Lemma (Lemma 9 in 13 Suppose that Assumption holds. Then, V logge((>) is L(t)-Lipschin, where [tt) is given by
Le
if log(1 1
I(t) = min - e t eet)))
Lg if > log(1 + 1) Lg
(1-
Therefore.
I(t) & Lu +
Lemma 7 (Proof of Proposition 7 m(()-strongly log-concave, where m(t given by
13 Suppose that Assumption 7 holds. Then, V log lo is
1 m(l) e-t/mo + (1 - e--)
Therefore,
m(t) > 1inn(,mo))
Combining these two lemmas, conclude that the Hessian matrix of log Pa satisfies the following condition
L(t)La & l' & -m
We will frequently use Grönwall's inequality in the proof. Below, we present a specialized form tailored the relevant processes.
Lemma 8. Suppose the Assumption (holds, consider stochastic processes Hę and G2 defined on the time interval [1,, 12). if they satisfy the same SDE, especially motivated by the same Brownian motion, which means that
1 { He + v log po dt dW,,
dH,
dG, 2
G, + Vlogrr-((G)) + aw,
then for each t E [t1, t2].
[H, - -Gll,, < (m(T=s)-))) HH,, Geilla
Applying Lemma g different processes and time intervals, we derive the following inequalities essential for our proof.
\ Yi|l,, Ee (m(t((((TT))) (2) ds Ynn - Yo 12 Vt [0,
(13)
||Y, - xx- 111, & -e- (6 ds 1 Xo 1
Vt [0,7].
(14)
which follow from the fact that same SDE as Lemma 8 by applying a time-shifting operator to y,-
T). (Ý. h). (Xt,0
50 t T) all satisfy the
21

---

B.1 Proof of Theorem I Part
In this part, we provide the first part of the proof of Theorem I with respect to Mareyama method. To achieve this, we first prove the one-step discretization error in the following proposition.
Proposition 9. Suppose that Assumption and 3 are satisfied. Then, the following claims hold.
(1) Firstly, it holds that
Ph
9EM
< (C)(((())++ + Vnn -0EM
L2
+ h (1n)) (((m)Ca(m) + c/ + Ca(m)) + Mr (1 (2))) (1++)))))+
+ 1/32 (nC)
+ h5
where
1 1 (n+1)h
((n)) 2 h Jnh
+
L(T - t) dt.
C2(n) = e -0 (m(T-0)-2) 1|o XTll,,
1 (n+1)h
Ca(n) T. J
(dL(T 12)
C4 sup ||eelaa.
(2) As result,
(Mn+1Ja - DeMille,
11,, < <TEM
EM ||Yn, JEMIle. 1I1o + 122CEM + /dC (n) + he
where
EM - pin 11 (1)((( - Mi). n
CEM (i ( + Mi (1 C2(n) C.)
Proof. We prove the two claims sequentially.
Proof of Claim (1) Rewrite display (12) in the integral form,
Vn % % //" (5ri V VlogpT-nh--(f. di
(x+1)e
dW.
Juh
For Euler-Maruyama method, we can write JEM in integral form as follows
(n+1jA
(x+1)
JEM
-EM +
o EM
(T - nh, (EM) dt + Juh dW,
2+1
2
nh
22

---

Note that Yo- = DEM, then, it holds that
1 h
(Y, ) dt + / A (VlogpT--h-- ST - nh, )raM) )) dt
|| 1rh
DEM
L1
2 0
12
ch (Y,- 8)) de + (((( (Vlogr--n-- PT Pn db
+
(V pr Poggr-nh--(000 (1 V log pr nn (DEM) dt
(V log PT nh (D - (T nh, (EM) dl
(15)
Here, we decompose the term individually.
into sum of three terms and then control each term n++
For the term of inequality 15 by Assumption 1 and Lipschitzness of V log
Pt. we obtain
2 - -Y)) + S-" (VlogpT-mn- - V Po (Y)) dt
=
ch
//h ph
| 8all, di + jo
L(T nh nhhnh - 1) 10|1/ Pr - $0l1, dt
2 Ju
(n+
h +
- dt dt
sup oggi
| - Yo L.
nb
[ro 8o
We then proceed to derive the upper bound for the term suPostec
Lemma 10. When po satisfies Assumption 7 holds that
(n+1)b
sup 0g
V.
Y
h
+
L(T- -t)dt || - -EM 1.2
L1
2
nhh
(n+1je
+
2
h
+
Juh
L(T - - t) dt
e (m((())))))) I - Xx 111
1
((+t)h
2
h
sup
|
Juh
(dI - -1) dt vah.
Notice that we have no initial limit on the Y, in Lemma 10 which means that can use this lemma to any discretization scheme.
For the term II of (15), first rely on Assumption to derive
11
V log (l0 V log Pr
(DEM) dt
lis
V log V log pT -nd (2 (DEM) dt
1/2 M1 +
||EM l1:))
23

---

Using the triangle inequality and 14 we obtain
MaMil 11
Yuh
JEM I1,, + nu Xxhll + X I Ii
Ynk
Dem IL + (m(--t) de 11o6 XTlh12 sup ogrET
(16)
For the term III of follows from Assumption that
III
V log PT (n sT nh, )EM)) dt
1.2
h
V log PT Iog (D - (T nh, ( Ille dl
0
<
hE
Combining these terms above, we obtain that
DEM
11,
(²(
Mi) |/nn -DEMI 11
(i(m) (Ci(m)C((m) + -c + Cs(m)) + Mı (1(1+(()))+ (17)
h
3/ (dC(()))
he
where
(i(n) 1 (n+1)h
L(T dt,
h nnh
(2(n) = e (m
1)- dt 3|oo XTll
((++))A
Ca(n) = h Jnh
(dL - t)) dt,
C sup ||
0E3ET
This completes the proof of Claim (1)
Proof of Claim (2) By the triangle inequality, we have
JEM oetille, Ile, < M(n+iya - En 111 + +115A 3EM
Y (I++i
(18)
L2
Applying 113 to the first term
18). we obtain that
Mm+1ia Yo < ((n-1) (2 (2m(T-1)-1)) Ynnh 8o1l
(19)
Notice that y, = JEM it then follows that
Yoo+Daa - 8%//1 I <
pin (m) -) Yun
JEM n
11
Claim (2) follows directly from our previous results and Claim (1) Since this step is independent of the discretization method, it applies to all the schemes discussed in this section. In the following analysis, we omit this step and proceed directly with the proof of the first claim.
We now proceed to derive the upper bound the Wasserstein distance between the sample distribution generated after N iterations and the target distribution Po, based on the one-step discretization error bound given by Proposition
First, note that
W2(C(0 ).Po)
Noll1, < VYA
DEM| I1, + +|NN Xolle
24

---

Invoking Proposition 7 of [13], we have
YYn Xolla, e So m(t) dt || |
(20)
According to Proposition by induction, we obtain that
YNK MM
11,
EM N-1
| (Nj
DEM
11,
(NiCRM, + 13/3 VdCC 1) hesc)
EM ]
\
11,
N-1
Σ
N-1
II
EM
)
(NºCEM +/3/22 1/3/2 /dCi + hEsc)
k=0
j=R+1
N-1
N-1
EM r
(1CCE + h/3/2 (d/2)))((())) + he2c)
k=
j=++1
(21)
where we define
N-1 EE 1= ]
1. Notice that
N-1
N-1
EM
II
1 (m(T-r)-2)dr (k) Mi))
j=k+1
j=b+1
N-1
M(Mamn - = hifmii -NN
j=b+1
Therefore, we have
N
YNN JEMIL
Mmmiin ()N--k-1) (RºCEM + 1/22 (dC) hEsc)
112
max OEkEN-1
CEM + 1/2/222
max
Ci(()) + hs
M(meiiu
0
Vdh
max UEKEN-1
Ci(k) + Esc
Mmin
1/2
(22)
Recall the definition of (1)) and the upper bound of L(t), follows that
max (1k) & + Lmax
-
2
and thus we obtain that
Limax + +1/2
|YNA EM IL, S vdh-
E.c
Mmmn 1/2
Mmii 1/2
Plugging this back into the previous display then we have
W Po) S - S3" mitt de [Xolle, + vdh. Imaa + 1/2
+ Esc
Mmii - 1/2
Mmii 1/2
which completes the first part of the proof for Theorem I
B.2 Proof of Theorem Part II
This part aims to prove the Wasserstein convergence result for the Exponential Integrator (EI) scheme. We will prove this theorem using the same method as in Theorem [ Following this approach, we first establish the one-step discretization error in the proposition below.
Proposition 11. Suppose that Assumption and 2 hold, then one-step discretization error for Exponential Integrator scheme is obtained from the following two bounds.
25

---

(I) It holds that
2(ch/2 Nnn D su
1% yE1
< h h²
Cs(n)Ci(m)++ + MI
a+1
h
+ h² Ct(n)) (Ci())a(()) 2/ + 2(04/2 +
h
1/32 (dCs(m)
+ h 2(eh/22 - -1)
Exc:
h
where
1 f
(((+ tt = Ci(n) 1
Cs(n)
h Juh
2"
(2) Therefore, we have the bound for one-step discretization error
2(e 1)
o + c, + 3//2 +13/2 VdCo(n)+ h h- (nntijh - 0 untille, II1, < < </
6s0
h
where
v = (m(T +)- (1) 122 (c)(m) + M1 2(24//221
Ce = Cs(n) (Ci(m)(((m++ + 1 ( +( + 2(c -1)
Mi(1 C2(m)) + CC).
h
Here, the constants C,, = 1, 1,2,3,4 4 are as defined in Proposition
Proof. We prove two claims in succession.
Proof of Claim (1), Consider the process defined in 12 which satisfies the SDE
1 Y, log Pr nh (vi) dt + dW,,
2
Instead of integrating both sides of the SDE, we use Ito's formula to e Y., then we have
d(e 1 11) = 2ah 3y, + e dύ, = - (V (Voogrr-n-= dt + dw.). and we notice that write it in integral form.
== 1/3% fod dd (t log (Y.)
{ fulltt e3((n+1)) dW.. nh
Then we obtain that
/' c log log - s
Yx El n+
nh, (E!)) dt.
We make decomposition the same the in 15 that is
V (X)) - (T nh, 5E = V PogPT-wee log
V log -hh (Yo)
+ V vogpr-nhh-(EE) V log PT n/ (DE))
+ V logPT-nM( su nh, )El)
then follows that
h
Yn.
e (1--)
4 PoggT- (T - nh, العد dl
1.2
0
e (1--)
loggr-nh- (Yi) V log pT nh n (1o) dt
e (h 1)
ni 1o () V logpr nk (, ( de
+
+
e (h--)
V log (v
(T- - nhh ,5 dt.
26

---

Note that apart from the exponential term (d(h)t), the derivation of the remaining parts is completely consistent with that of [T]) until we encounter the term involving @e! at which point we obtain
nh -() ( -hh
Yr
Yo
n+1
sup og
112
1/
(%/² c3(0--0) at Mih (1 (1+ || (111,,))))
h + % ct(h-e) dt
Esc
By Lemma
10
we can bound the first term on the right-hand side of the previous display. More- || L,, can be bounded similarly. Substituting all coefficients with (2(n) from
over, from Proposition
16
we obtain
< h²a Cb(n) Ci(n) (|Yn I.l C((())))()) J/c, ca(m))
2
12. 2(e4/2- M1 1 + |Yn JEMiL, +C((n) + C.
h
1/2 VdC5(m)
+ h. 2(ch/22 - Esc
= (ol())((m)) + Mi 2(eh/2-11
|nn -DEM! 1
h
ht2 (C(m) (C(())C2(n)) 1 _ +c((m))
2(eh/2 -
M1 C2(n)+Oa)) +Ca)
h
+ 3/2 (dCs(n)
h. 2/eh/22 - 1)
Escr
where
C2(n) = h nh
(n+1)) (
ed(m+110-0 L(T -) dt x0 Ci(n) 2
1
Proof of Claim 2
The proof is omitted for brevity, it merely requires incorporating
Ym++1h 13h Y
into the conclusion of Claim (1) following a similar argument as in the proof of
Claim (2) in Proposition
For the second part of the proof for Theorem [ recall that in the first part, the three key steps 220 and (22) lead to the desired result. We now revisit these steps within the framework of other discretization schemes.
Since 20 is independent of the discretization scheme, we can directly apply it throughout the proofs of Theorems Bnd For (2I), we note that the 1²² term in r is neglected, which results in the
same upper bound for N 1i k+ rf across all discretization schemes.
Given the consistency of these two steps, for the remaining discretization schemes, we can directly derive an analogue of (22) from Claim (2) Therefore, in the subsequent proofs of these theorems, after establishing the corresponding proposition, we proceed directly from an expression similar to
22
27

---

For this theorem, we begin the proof with the following inequality
|
UN
1
hC + 1/22Va
max Cs(k) + Esc
Miii 1/2
N
1
Vdh
max OEREN-1
Cs(k) + Esc
Mmmin - 1/2
Limax
vdh.
Esc
Mmin - 1/2
Momi 1/2
Combining this with the bound of |Xoo Ynn l1.. we obtain
Lruax
W2(C(UN)) S mmin T | + vdh.
+ Esc
Mmin 1/2
mmin 1/2
as desired.
B.3 Proof of Theorem I
In this section, we prove the Wasserstein distance between the generated distribution L((REM) and the target distribution. The following proposition established for the one-step discretization error.
Proposition 12. Suppose that Assumptions 7 2 and A are satisfied, the following two claims hold.
(1) It holds that
Y/h
0
REM a+1
L2
1/2
Iu U|[(T (n w)h) ( ( L(T (nhi- + +Mi 1²
h22
du du
1
L(T mh) + 8v3
|| -OREM 1ro
4
/ [w (w- v) LLT
(1 C2(n) #2 -/ca (dL(T- (aL(T-nh))11) L(T (n w)h)
L(Th) -nh))
2
1/2
+ M1 + C2(n) + (ca)
dudu
(Ca(n) + C4) + 1
[(T L(T-nh)C2(1) + (dL(T- - n2)))
8/3
4/3
/ L(T (n + u)h) ol du
1/2
+ 13/2
+
2/3
+ 2hE
(2) Moreover, holds that
[[( -REMI || < T
REM
Yaa 0ulli, + N²CRRM + 1/3 33/2CCRMM
n.2 + 3he
where
.REM
c p/m 1)th (m)))=))=+)))))
n
1/' [ ITT - (n + u)h) 11 + I.T mh)) - -M )2
1/2
+ 12
dudv
1
L(T - nh) +
4v3
8/3
28

---

11.1 (t (w--w) w) [( L(T nuh)) C2(n) + 1
CREM
+22 c + (dL(T -di)) L(T (n + u)h)
n,1
+ M1(1+ C2(n) + Cul 2
1/2
dudu
1 (C2(n) Ca) 1
[[(T -2h)Ca(n)) + (dT(T-
(dh))-nh))
8V3
4v3
1/2
CREM Va
I(T (n + u) 1)) v| du de
2v3
Proof of Proposition 12
Proof of Claim (1) We make the following decomposition of one-step
discretization error
REM n+1
<
Yx Bu, [ORMM] 11.2 + |uu, [PRT1]]
REM n++
(23)
1,2
We first derive the upper bound for the term P/ - Eu,, [
By the definitions of RREM and
Yn, we have
Yn Eu,, [OREM
12
h dt + S.¹
VlogPr dt hBu, (DREN) - hEU, (sa (n + U. @ 11
Notice that
l Y. h.u, n
ch V log PT nh IEE VoggPrnnnUUsM (Yr,,))).
Plugging this back into the previous display then gives
Vn EU, [
| WEB [Yrr.n) + NBU, A log Pr- (Yu,n)) 1 ,,, (DREM) -EEv, (s- (n Un) BREM ))
FREM 1 h Eu, (Yo.n att
h Eu. lot (Yu..) s,, ST (n + Un) 2 u+u
By the definition of Yu,n and RREM we have
Eu. (Yu.h - BREM.
11
Unh
Unh (VlogpT-nh--(P)) S. nhh, gRRM )) dt
Run
(V, -F,- (REM) di +
11
Unh
Unh VlogPr-nh-o(X (T - nh, 1 REEM
EU.B
1/
RREM
dt +
I||
df
2 J
BU..
1
REEM
di
V Vlogpr-nt-s(9) (T nhh, gREM
di
2
GREM
dt + Pog (T - nh, REM dt
The second inequality arises because the integrand is non-negative, the last inequality follows from the fact that the random variables inside the inner expectation Ey,, are independent of Un, and thus the inner expectation can be ignored. Then using the same argument in the proof of Proposition g
29

---

especially adopting the same procedure as the one following [5]. we can apply the conclusion of Proposition to the term above, then we obtain that
(n+1))
Buu, (oin - v
REM
h
-
sup Ogtis
1/ - Yo |12
2
1.2
nh
(m+
log 1og ((UREM 8, (T' nh, (REM) IL1 df
(C) M1) Vinh (REM 11a,
122
+ 112
(i(n) (G(m))2(m)) + +22 C4 Ca(m)) + +M C2(m) + C2(n)+
+M((1+ Ca)
+ 1//2
+ hE
2 1i2ri |Yar
-REMM l1, + hhr2 + 1//3//2// /2// (m) hesc:
where
(n))) + M1,
(((())2(M)+ + + M1 Ct).
C1 (n)
We now derive the upper bound the second term in Note that
Evs
v loggr-(ntU.i)) (m+U.))) (T - (n + Un) REM
1.60 (Vlog PT Progrr-m+uu) OY ((n+ujn) 8. - (n + u)h R
du
L2
| VlogPr- (Y ST (n + w/h, REMM
du
//' (IV ( (M(ta)))) V Pog REEM)
(24)
L2
| log Pog s. (T (n + u/h, RREM
du
L2
< L(T L(T (n + w)h) P(mtaja
gREM
du + Esc
the second inequality follows from the triangle inequality, and the last inequality depends on Assump- tion and 4 By (17). changing the value of wh, we have
Mintuya
REM
(C, + M1) Pnnn -REM
L1
+ (uh)² Cin((u) ( (C)n(()C((m))) + f. + Can(()) M1 ((())
(25)
+ VdCi,,()
+ whesc:
where Ci, (u) and Cz, (u) is the wh-version of C₁ (n) and Cz(n). respectively, that is
1
(n+a)h
CI,(u)= = +
L(T t) dt,
2
wh nh
(n-ujh wjh
1
(n) uh. nh
(dL 12) dt,
30

---

Plugging the previous display 25 expression, yields
back into display 24
then rearranging and simplifying the
Eu, v 1o
REE - S. -n -n Un)h,
U.
< 1/2 (o T(T - (n (() + Mi) Mi)du) du) l|Ynn -DEM llt.
h {oo
L(T - (n u)) C2,(()) (2)(( ] Cs Csm((w) MM1 +ca) du
(fo L(T -(n ++)u))/>² du (n)
1/12
L(T (n +u))uddu (f)
+
+ Exc
- rs Ynuh~ -DeMll 11, + hhrr + 1//3 rss + hroEsc Escr
where
[ L(T - (n + u))) (CI,(u)) M1)
33
[. L(T (n + u²) (r,((w)) (I,(())((()) + I -/cs + Can(u)) du Mi + C2(n) + C4),
76 ( L(T (n+ w) du) VdCi(n),
98 I(T (n + u)))udu
From the bounds we have obtained for two terms, it follows that
Vn Be, (DRM))
1.2
1 h/ (55 + r3) Ynn DPMI Il, + h (r +22++ra)+000 + + (dCi (VCC))) +rs) (arrs=h+ =h²( + ro)ecc + he
(26)
Considering the second term of one-step discretization error
@REM EU., [OREM]
7++1
[BREM
Ev, RREEM [] + [S. (n + Un)h, gREM Eu; (T. (n + Un Un),, RREM
h(Un OORM + hhUUn - )5. (T - nh, gREM)
2d
1
Juh | ((nnttr) dW, - S.' (flleta
(n+uan
dW,
du
nh
+ (T (n + ), RRMM)
Eu., (T (n + Un)h, 0
(27)
The second equality follows from the fact that
(n+Un>h
Ey,, [OREM] = .REM + hEE, [U,]UREM + +hEr, [Un]s. nh, @REM + Ev
dw,
2
nh
2 hs. (T- - nh, gREM +f, (Jr
(intaja
REM + REM
dW,
du,
since U, independent PREM
We proceed to bound each term in
27 For the first term, still notice that the independence between
31

---

Un and @REM, then we find that
J [2]
1/2
UUn
gREM
EU,
(Un - 2
1/2
EUs [1Un -100|RRMMM1 |1²
1/2
@REM
||
12
||OREM
2/3
L2
The bounding of another part of the first term follows in a similar manner, we obtain that
(Un - 1 ( nh, -h,, 0 @REM)
1
(T nh, RREM I|اa
2v3
1/
For the second term of (27). notice that due to Itô's isometry formula, for any well-defined stochastic process X₂ and its Itō stochastic integral (X) = So Xa dMu,, we have
[X) { x2 d2
(28)
then we can establish lemma.
Lemma 13. Suppose Wt is d-dim standard Brownian motion, then
(ntUw)h dW, - aWa-Jo (1, ((xtu))
dWe
du
<
For the third term of (27). we get
s. (T (n + h,) (REN) - -EU, [s. (T (n + Un h,
REM 6
8+ (T - (n + Un),, (REN) Se (T - (n + u h v REM du
La
2
1/2
E J
ST - (n + u)h, (REM) (T (n + w) ())
do
du
/ [T - un+ (REM) (T - (n + 1 hh RREM) du du
1/2
E J
/. [T - (n + h/ REEM) (T (n + v)h, RREM) 11
1/2
du de
Then by the triangle inequality and Assumption we have
s. (T' - (n + h//, (REM) s. (T (n + wh REEM
L2
(T (n + w)h, RREM) log (n+ (1REMM
12
||s. (T (n 1) )REM) V log PT Po
//
REM
(29)
+ |/ log logpr-(utu)) RRM)
V log PT-PT (REM))
L2
< 25 + log (m+a)) ( (RREM) V log Po (n+ (v REM
L,
Combining the three terms of (27 together, we have
||
Ev,, REM1
2 ||OREM
IL, + 4v3
//2
1
122 I|s nn,, REM
+
1//2
8/3
2/3"
h S. J.
1/2
V log loggo- (n+u)) (0 (OREM) V PlggPr-(nta)) (1 (BREM)) 1|=
du do
+ 2/Eac
32

---

By applying the same technique used in the proofs of Proposition 9 a Proposition 11 the upper bounds for | L and [|.. nh, @ III, follows readily. Thus, the proposition follows
immediately from the bound on the second last term. We now consider the case for u > v, due to Assumptions T and
V logpr- Po' wn) (DRFM) V log Pr- ( (DRTM) IL2.
< L(T (n w//h || MIOREM ORTMILL. + M [1 1
Since
| 12 + uuh 1|s. 1.,,,.RRM Mlro + / dW:
|| R+u
19 REM
12
2
12
2 (u w)
(1Ynh - @REMIl 11, + C2((n) + +ca)
Esc + L(T- - -n)) (1Yan - BREMIl,, - -c((m) + (dI(T--hi))
+ (u t)
Vuu
v)h
(u- v)/
+ L(T L(T'-hh)
||nn -0 @REM
2
L2
(u v)h
[ + L(T -nh)))2((m) + jot + (dL(T-nh))'
1
(u v)dh - (u - w)hsse.
The second incquality follows from [6). Assumptions and Lemma 18
Similarly,
|| R+v
< ||
|REMM
l|r, + -| l1.o
n+o
1
< vh
+ L(T -nhh) ||Ynn -YEM-
OREM
2
L2
( + L(T - (2)))(())) + _ + (dL(T- 12h))
1
+ vh
+ Vodh + vhs
+ Inn GREM II2, + Ca(nn) Ca.
Therefore, we obtain that
V log P (UREM)
V log (n (aae))(((RR0)) (BREM) Le
< h
(u
w)
+ [TT - mh)]
L(T - (n + uih) + M1 11
RREM I1Le
2
+ M1M
+ L(T - mbh)]
|| - -
2
+NPoM1
(
+ +(T - nli) C2(m)= + 2
f/c,+ + (dL(T-
2
+ -/3/2M/Ved
1 + h (u - 1) [2= +1+L L--LT--N) Ca(n) / c, + (UL(T-A)))V?
L(T (n + u)))
M1 + +g((n) + +cii
+14/2 - (y ())) vld
+ -)] - (n u)h)Eac + 122 M2
33

---

We claim that we only consider the lowest order of each part, which means the relative higher order term with the combination of d and will be ignored. Then take the supremum with respect to v,
{ L.
1/2
V log (r+to)) (ato) V (r-tta) ( I du cu
h !/" L--
||- u|T(T- u)++
1
+ L(T -nhi)
+ Mi
du dudv
1/2
||nn
gREM
11
2
[-¹ ! T(T (n +-)))/ rl du du
1/2
+ 1/0 L' (u- o)- (n + du du
1/2
E0c
(30)
Combining the above,
10 RREEM- EU,, (BRMM1 L2
h2a
(|Y'nh - (REMM lle, + C2(n))+ Cc)
8/3
+
1
li²
[Ean + L(T -nh) (1Kmn - -RREMI lا, + +Ca(n)) + (dL(T nh))2/2)
4 4/3
1/222
2/3
[ L(T -nhi) L(T--nt) + MM] 1²
1
1/2
12
[u - v/ L(T (n + u)hi
du du
1|Yrnh - -0 -RREM]
[ (n u) ( + I(T -nh)
C2(n) 1 C4 + (dI(T nh))
+ 122
I(T uu))
2
+ +1 + (2 + +ca)p ?
1/2
du
1 L(T (n + un)) r/ du du
1/2
3//2
122 [/ 1/' (u L)T (n w)n)² ())))dudr) + 2he
1/2
Esc
3
/- [u -u- LLL (n + u)h)
1
L(T -nh))
+ Mi
2
du du
1/2
1
+ 4/3
L(T nh) +
Ynn @REM 11
8v3
222
1'' { (u - 0)
[62 L(T -LT-- nh))
C2(m) + t + (aL(T (d)((-nh))112 L(T (n+ u)h)
2
+M (M((1C2ggm) + +cus } du
2 1/2
1 (C2(m) C.) 1
(L(T nh)C2(()
- (dL(T-
8/3
4/3
Va [/' / L(T - (n + u)h)2 v| du
1/2
1
+ 3/2
+
2v3
+ 2he
(31)
34

---

Compared to the term Yn EU,, [OPTT we can focus on the lower-order terms, ignoring the score matching error. Therefore, we have
Mootaa Dtt1
-UREM
1.
{ 9166- Lo Iu L|L - (n w)h) + L(T -ahh) + Mil ]
1/2
< 122
du du
|| nnn - @REM IL.
L(T - -nhi + 1
+ 4v3
8v3
1'" { (u - () ((: + L(T nhh) C2(n) + 1 c. + (dI(T - (an)(/n L(T (m + u)h)
122
1/2
[M + ci) dudv
1 (C2(n) C.)++ + 1
(L(T (hh)c))) (dL(T- m2)))
8v3
4/3
+ 3/2 Vd J/ [ L(T - (n + w)h) - wa du
1/2
1
+
2v3
+ 3h
Returning to the proof of Theorem 3 by the conclusion of Proposition 12
we have
YNn
N
11, 2
Mmin 1/2 h max
CREM
+ 1/2
REM k,2
+ 35
OEKEN-1
OEKEN
Va/3Tmax + 2v3
3
S
Vh.
Ean
Mmin
1/2
Mmin 1/2
This completes the first part of proof for Theorem
B.4 Proof of Theorem 3 Part II
We begin with the following proposition. Proposition 14. Suppose that Assumptions 7 and|aaa are satisfied, the followin two claims hold
(I) It holds that
Yn
REI
L2
0 1 1u v|L(T (n u))) ( L(T (((r-nhh) -nti) + M1
< 2
dudu 11/2
1/2
+ 2 [u ia ITT - (n + v) ()) (w)
||
0REI
1
1 (1-1)) 4)4 (1- wh L(T - (n + w)h)
+ 2/2
(2(n) + -
nh)CC(()) - (dL(T nh))
h
+ (1 Mi [ + toh - +dL (dL(T' (h))441))) n1)))/2) + C2(m)
(ež(1-w))
et (1-00)
[L(T (L(T- (n+ (dL un))))1)))
+
h
35

---

+ +10/2 Vo S/' / L(T (n u v/ de
1/2
+ 30
(2) Furthermore, it holds that
Mn++IJ DRT1 1nt1ll1 lll
REI < &n ||Y:nn
oRell 11. + CR, + 1 n, + +3hssa
where
.EEI n
e Jun 11n (m((T-)--)dd
122 1/ 1/' Iu o L(T - (n (n+ unn L(T -hh)
+ M1 + マタ v|L(T (n + (v)) 12
du du 1/2
e (1- vih 2
CREI
e I(T - (n + w)h C2(m) + -T) (dL(() (dh(T-nh))
(dL(T
h
+ et (1-)A M2 1 + ub
(L(T - + ( + C2(n) (c())+ca)) +C
+ le um- -
- v()) + (n
(n+)))))))
0, Val I/ [' L(T -
1/2
L(T (n + u)h)? uu - o| du du
Proof of Proposition 14 This proposition can be proven following the same approach in the proof of Proposition [2] with the only difference being the inclusion of the exponential coefficient term. However, this term does not significantly affect the overall proof. Similarly, we make a decomposition as
vh
0
REE
< P/--- Ee, [UREI
+ Eu, [OREL lttlllaa
(32)
12
12
Note that
Yn Eu,
[RREI
//¹ e (1 VogPr--Woo
hBu, [c s- nh Unh, GREI U..)
e (1-wah (Voggrr HlogpT- nh ub un (Yun) - (T 17--m nh - -uh, )REW)) du
h e 2 (1-w/h ( VVogPr (Yah) - V Vlog pr uh ( (URtL)) du
h 1' (2(1-w)h (o so nh uh, D du.
Then, we obtain
Yn - Ev, (OREL
< h /-/ [e -no (Yah) - -s, (T nh - wh, URE >) du
cl--w) T(T- - nh -h) /ann OREI l|., du du+ h L''' (1 duese
< h
(1(1-)) L(T nh uh)u² (G,n + M1 2(ouh/8-11 du |Ynn (RRE!
13
uh
11g
+ 13 ff (i(-- L(T nh uh)u² Cam(u) (Lm(())2(m))) + c/ + Can(ul)
36

---

2(euh -1)
M2(1+ +1(1))(()) + Ca) du
wh
+
- nh - uh)uu//ca du vd
e e2(1-u)h I(T nh ul) 2(c 2(e
du Ecc 2(e - 1) Excr
(33)
h
where
(n-wjk L(T
Cu, uh nh
In the third inequality, can directly bound
Yuh GREI
as it special case of Proposition
where the step size replaced by uh. For the second term of (32), we have
n+1 Rv,, [0,
[E x+1]
-he-(1-U.)h S. -hh O n+t) h.EE
[100U. nh Unh, ontu) )
= =h [ [ddlu-Uu) s. - nh - Unnh OEEI
ontte,)) e2(1-0)1 (T nh - wh, ) ) de.
Similar to display
29
we then obtain
UREE EU,
REI
11
2
1/2
E 0
h
edll-u)/ss ST ( u)h, n+u) - (1-00) ST (n + v)h, DREI a4p) de
du
/''' I.'
1/2
h
c [ (T (n + uhh, ( (1-1i u/ s. (T (n + h)h, )
du du
/ /.'
|| |dd [ed0--i) (3(1-u)h with v log (o (URE)) (1) v log Pr (n+e)) (1 (
2 du dv
1/2
+ 2h (f( cll-w)h du
1/2
Eac:
Using the same strategy in display
30
we arrive at
( log (RREI)
(1 wit V log (n+vi RREE
(1-1h V log (n) (U n+u) - V Po u)+ REI
+ (1 log (DREL) m+t) V Vlo PT (n bDj (OREL 1
+ - (1- - (1 V VlogPr-(ntojd
(OREL (OREL))lE II
(1) L(T (n w)h) |OREL gREI
+ + t(1-a) M1h(1 + R IIL.)
= - (1-0m [(T (n + v)h) (IlYatoja - oRttella + Ca(mi) + (dI (n 2/2)
The second inequality follows from Assumptions and
We bound the term
V VlogPr-oo-o> (ØREI
l1, by decomposing it as follows
ntt
log (n--r-(p-o) (D (DREL
V log (JREL) V ( (Mnte))lles
L2
v log Pr PT- ((tt))))= V Voo loggT PT (atoja (X (irttnnaa)l
+
V (X
37

---

Without loss of generality, we consider the case where u > the other case follows similarly.
|
AAtelLL
-eeth - ežin)) || L. + wh
L.m et' dt
dt (T 107--h, nh, (REI) 111,
+ / (Z((++())-)) dw, plater e(((+t))=o) dW,
nh
(ei(u-o)i 1) (1Yau GREI 1l1, +C(nn)+Ca)
+ 2(e2((u-n)) 1)e2mh @ec L(T -5 nh) (1Y,u - - llea + Ca(n)) + (dI(T- - n2
+ [ -1) +
2(e""" -1 2((
-1)
E " 1/2 Vd.
(34)
Here, we apply the formula in (28) to bound the last term.
(ntu)h
- eb((xto))-0)) ((
+ (o dw.
J (lattsh
(atu))
1/2
Vd
( dt +
(ntu))h- dt
(m+D)a
vd
(ceh
-1)
2(e "
" 1/2
We then bound the term
|| lL, following display (34 above. To this end, let 0, we then have
|| |
AaRtillee 111
< (ctth 1)((Inn
@REII 1|اo +C2(n) C4)
+ 2(etvh 1) Esc + L(T - nh) (|Ynhh OREIL 112. -2((n)) + (dL(TT nh))
+ V - 11/1/2
+ Yun @REI 1|1 + +C2(m) C4.
Additionally, we can bound || [((ntjja @REI L error under the Exponential Integrator scheme, where the step size is replaced by vh. Specifically, we have
is a special case of the one-step discretization
|| (ituja - OREI l1I,
< (v) || - OREILL, + 122 + 3/22 VdCo.n(1) 12 + wh 2(e 1) vh Esc
where
ia (v) = Sun plntwih 11)- (1) dt + 12 (1) + M2 2(e2vh- -1)
wh
I c + Csm(v)) + 2(e2 - -1) M1 + C2(n) + C4). (el( (o, (Cin(V)Ca(n)+ +
wh
Then, we obtain
11 (1 Dogpr-(ntu)) (DRtL) 3 (1- log PT (m+uj) (OREL
< 2() et (ll-a))lettohh L(T (n + w)h) 1 L(T - mih)
+ + (1)h M₂h [eetth - -1) + 2 mh I)L(T nh) +
38

---

+ - (((- L(T - (n+ ())r)
Ynn
+H-LTT - (n + wlhh) | ef(1-0}/ (v)
h
let(1-ujh e2(1-1ij)
+ + (n + b)11
u 3/2 (U)
h
+ + et uh -1
(C2(1) +C1)
h
+ +///2 Mh u/- Vod eeth -1) -1) 1/2
wh
(2(1--0) ed(l-uhh h
h
+ @ wh L(T (n u (0-+ Ch))((())))))) + (dI(T nh))
-dill M1 [1 + nh (L(T (aL
Cc+ +Gp
1e2(1-
())-
[I(T (n+ + (dI(T - (m+
h
(euh - 1)+ (e - 1) 2( - -c
+ L - (n (1)))
Vd
h
let(1-ajh- ef(1-njk) 2(e -1)
+ h (n + v)ti)
h
eduh edch
+ hE 20- (n (n u)h)
+ Mietoh
Ignoring the oideer-order terms, take the supremum with respect and substitute it back into the original expression, yielding
1
REI 1+1
EU. n+1
11,
1/2
e ( with sx (T u)h, DREL - s1 (n v)h, Ø 2
h
du du du
L
//- (1-vv[((2) - - uw)) 2 L(T whh)
<
+ +M,
1/2
+ lu - - - u)krrr du du
|| - &REI
+ 1/222 { {o L' dL(T (n + uu)) of dauu
1/2
+ 2hEac
(35)
This completes the proof.
Now, we have
REI YNN
2
1
h max OEREN-1
CREI
Vh max
CRREI 35
11g
Mmin 1/2
OgkEN-1
Lmax
27
dh
+ Exc.
V3(mmin - 1/2)
Mmnin 1/2
The desired result follows readily.
39

---

C
The proof of the upper bound of error of the second-order acceleration scheme
This section is dedicated to proving the Wasserstein convergence result for second-order acceleration. To this end, we first establish the following proposition. Proposition 15. Suppose that Assumptions
are satisfied, the following results hold.
(1) First, we have an upper bound
Yn:
as follows.
L1
v/ -55 Ir, & An.,e(L(uu)-b))) 2 ||an 1 111, + An, (
+ (neset
+ 2 van (L) +
1 n22c(0))
cll(TTnb))
where
1 / (1/0 [ L(T- u) nh -u) An,1 nhgte(n+1)W
sup
12
(L+ nh)] du) ds,
1 S ( [( L(T (h - UT nh - u)
An,2 = nn
sup
12
+(1+ L(T mh))
L( L(T nh)] du) ds C2(n)
va 1.
L(T nh w) L(T nh u)
+ 32 I(T - nh) I(T -h) n/) du ds
1.' (1,1
L-T-LLTTTi (I(T nh - u) + - mh)))) nh)) ds Ca
V2
+ VdLF
4
(2) Furthermore, it holds that
500 ([(11))
50
Ynth - -050 I11, 1:2050
+ hE +
2
Van3/22000
3/2
(L)
+ sec
1 (M))
(L(T-nh)) 2)h
3
2
where
Jun (1 (m(T-c)- + 122 (((-nN)))))))
50 e
CSO
A (4T-
Proof. Recall the expression in display 44 which states that
-50 + /² Ft
yso + (0500) + La(x, - (50) + Mn -nh)) do + S" aw,
nh
with
Ln Ia + 2
log Pr (1 Raxdd
2
8²
M, Σ
V log Pnh (050)
V log pr (1 E ERd.
dx
Øt
j=l
40

---

Plugging the estimates of V loggrr Pr-h (0x) Ln and Mn into the previous display yields the follow- ing process for 50 2'7
2t
50 5050
+
50 (T nh, 100 (2 50
+ ( nh, 500 ds
2
nh
(1)
00 8x (T nh, 1 ( nh) ds { dW..
(M)
nh
nh
Then, we obtain
Jt It (t- (cdnn)(ss((1 (T. nh, 50 v
50
+ (1) (T - -nth (55) Ln) [ (1
(50 do Iu I
(1.50 xz)
n/
nh
(s (T nh, 500) Mn) -1 nl))
2
Notice that
T la 2 lt (050)
1 < T.(T nh)
(36)
2
Combining this with Assumptions B and then provides with
1 Iz IellL.
[" 11:50 100l,, l1,,2 ds + ||TnllL, '" 11::50 - -xa|l ds + 2 (ch)=())) =nh)²((2))
< (t nh)esc + +o66) (L)
nh
< [(T - nh) [nh
11
[ 15
2 zzI ds + nh)E Vddtt n2)3/22210 E00 + 2(( nh)²
1
2
nh
2
We omit the constant of term
Esc in the last inequality. To handle the resulting integral inequality, we invoke the following Grönwall-type inequality.
Vd((t nh)3/22
3
Lemma 16. Let 2(t) > to satisfy the following inequality:
z(t) < ait) /' ((s)z(8) t> to:
where B(s) is non-negative, and to is the initial time. Then, the solution z(t) satisfies the following bound:
2(t) a(t) + [ ((s))(s) exp (()))))))) ds. t t> to
Additionally, af) is non-decreasing function, then
x(1) < a(l) exp (()))))))) 126.
exp
Let
=(t) 10:50 xeel1a
(tt) = (t - -nn)=ec
+
2
Vd(t - nh) (2(0)2 + 2 n/) ((M)
3
B(t) L(T nh)
and set nh. By Lemma 16 have
||050 Mm+1jallz, & hesc +
3/2
+ (M)
2 (L(T-nh))
(37)
Ese
41

---

The original SDE can be rewritten as follows
Y = J' (6 V PT nh- -Y- ds
fm dW..
nh
nnh
Combining this with the definition of ot in we then have
Ye- Ino+e = L/' Yr, + V Po (Y.)
-0o V Pnnh - Ln (Inhts 0000) ) - -Mns ds
2
2
I. Ln Inh+t)
ds !
1 y, + v log PT Pogpr-nh-s(f.))
1 U00 Vlogpr-nh ds 2"n
1/" (Jo" Im dr du ds S:d (Jo Mn du ds
S: Ln(Y, - ds I'" (To d
( +V VlogpT-nn-coo Plog ds
1.' Ln di du ds L'¹
Ln
(No Mn du ds.
1 Yu + V log Poggr--nb-u(f.) Recall the definitions
We then apply the Itô formula to the term
(
of Tn and Ma and after rearranging the expression,
obtain
Y Inh+t [/¹ LÝ, - Inh+t) ds + S
42 log PI -nb-u - V² log pr nn (%)) dř,
ds
II
0²
log -nh- Oe log
pr -n uln)
0x 2
82
log PT Pn (%)) 8, log Pn (8o) du ds
2
dmy
In what follows, derive the upper bounds for each term the right-hand side of the previous display.
Upper bound for term 1: The upper bound of the term follows directly from the fact that
| < L(T - nh) 2
Upper bound for term II: To derive the upper bound for the second term, we expand the term dÝ,. yielding
2 logpr-nh-uuY, v log (Yo)
L2
1/ ( pr po v2 log md (Fo)) (
(r + V Voo logpr-nh-u du
2 log PT (Y,) 2 log (Yo)) dw,,
L2
/ 2 logpT-nt-u 2 logpr- (1))
Y, + V VlogPr-nh-of du
(1/1 ( PlogPT-nh-UXX 2 log (2o)
1/2
du
42

---

The second inequality follows from display 28
We note that by Assumptions 6 and 7 it holds that
v² logpr-nh-u(F,
V2 loggr -n (Yo)
2 logpr-nt-u(f. V2 log (Fr)
11, +
2 log Prghr-nt(Pa)
Volog Pr P nn (Fo)
& (2h(11 + Yu ) + LF 8/- -Volls
< M2h + (LF + M2li)
\%x- - -Volla + M2he 1
I LpVdu,
Combining this with the previous display provides us with
v log pr -
V nn
La
L VEvdu- du
2 Ya + -logprr-mh-(kk
du + (1/
Lj
1/2
Hence, we obtain
-
11. (1, 22 log ProgPr-nh-iie Pog (oo ds
La
11/2
2 log PT -nb-a
V2 pog (Po) dř
ds
L1
1.
[ LFV du
1 Yv + V Vlo logpr-nh-u(fa) du (%o Lid du
1/27
ds
v2
4 LEVat.
(38)
Upper bound for term III:
In Section it is claimed that the partial derivative of V log with respect to f can be estimated without requiring additional assumptions. This is achieved by transforming the t-derivative into x--errvaaivve via the Fokker-Plande equation, as detailed below.
d +
1 T Vp((x) 2t
1 +
(x)
(xp((x) pp((
(39)
2
dx²
We need the following auxiliary lemma.
Lemma 17. Let Pt be the probability density function of Xx, then
82 p2)
Tr (o=lgg((x))) + [[Vogg((>))
Ox?
P((x)
1=1
8 (x)(
V logT(())gg<(>)))) V
dm?
(=(x))
(oV' + | Vogg (z) P.))..
We begin by taking the gradient log Pt. and then compute the partial derivative of log p: with respect to f. This results in
(Vp((t)
8, V.(r))
Vp((a) Oxp((x) px(x)
a, lo) =
pc(x))
Pc(x)
pz(x)
43

---

Under certain regularity conditions, we can interchange the operators 0 and V in the term d Vp((z), and substitute (.p,(z) by (39). it follows that
(x)p(())
(d 1 T V Vlogp((a) 2
1 d
8 px)
(z)
logp((r) Pe(x:)
2"
+
2
Ox²
Pc(x)
1=1
and
V0(p((z)
d
1 2" Vpi(x))
(x)
V
22((2) +
+
Σ
pr(x)
Pt(x)
dx²
1=1
d Vp((a)
1 Vp((z)
ppp(
1 v
Σ
pp((z)
1
2 (x)
2 p((
Pc(x)
2
dx?
pe(x)
d + llggp (x)
1 \2p((<)
1 v
(x)
1
Σ
2 (x(x)
Ou?
Pr(x))
Therefore, we obtain that
V.V logpe(0) d+1 V (o)
V'p((e)
d
8 (z(2)
+ 2 (+(x)
+ 2
2
Ox²
pe(x))
1=1
V Voggt(zz)
1 -
T V logpe(z)
+
1
pe
1
2
2
dx7
pr(x)
1=1
1
('p((z)
+ log
(1
2
(r(x))
d 8
8 p((x) d
1
log (x)
dx?
pe(x)
ax?
(x
By Lemma 17 the last two terms above can be calculated. Additionally, it holds that (p)
V'p((z)
V log
Pr(x)
pr(r)²
Thus, 0 V log (z) can be simplified to 0, V (x) マ (og + 1
V([V)) [(H
v
Tr(V2 log [Vgg((())) -
v (z)
Tr (o
+ |
2
v2 logp(())= + (oT((\)og)()))) + _ Volgp)))))))
log
+
2
2
2
Notice that
d
82
V log (x))(()) = V (x))
Ox
1=1
Then, it follows that
d
82
(ogg:(x) - (z)
2
dx
1=1
V (og + 1 v logg:(<)) ((
VV (o) + (zg + log v Vlogp(()).
2
44

---

Therefore we have
log (og
V log Log
dxy 3
1
1
V log PT nh
42 log Ps (Yu)
Yu
2
L1
L1
2 log PT (u)
< log pr (v))
V log PT
(8
+ 1 L(T nh u) Yu
+ L(T - nh u) V -nn- -logpr-nt--uaal
uul (Pal Ie
The second inequality
follows from Assumption
The bounds for
and Yx
1.2
log PT -nb can be derived according to the proof of Lemma 10 We then find
||,
1/2
8 V log PT loggr-nh-u(Yu)) V log log
2
V log PT MogPT-nh-uuYY)) (Ya)
+ + nh -w)
+ L(T nh w) V log PT P - (PW)
2
1,
L2
+ L(T nh -1 u))
L(T nh - u)) 10:0
11, + C2(nn) + (dL(T nh - 1/2]
+ 2 L(T I(T - nh wu) (1) 055 11, + C2(nn) = +Ca).
Therefore, we obtain
8
log PP nh- 0, log pr log
du
2
j=1
82
+
v log JogPr-nin logPr-nt 0 PT Pnb (8))
du ds
2
dz3
(11 + L(T - nh
(T)nn)))((TT nh)] du ||Y.h y50 l1a,
1.¹ 1 11+ L(T - n/ -d (1 L(T -T)) -nh]] ds (2
( LLT - - nh 2/2 + +1= L(T n()) n2/i2/22 du ds
L. I (L(T (L(T nh -
u) nh-
L(T nh)) du ds C4.
(40)
For simplicity, focus on the lowest-order term. Recall equations (36). 38 and 440). which lead to the following expression
1/, xnntee Inn <(L(T- (Lhh--nh)-
nh) - -h)- 2) 1/ [111- Inh+s ||,, ds + (An.1 || osolle + An.) -t?
where
1 //' [' [[1 + L(T nh w)) nh u) (1+ L(T -Th))((T du 12 J
An,1 nhgrg(at1)
sup
45

---

and
sup 1 (2(m)) J/" [ [(1 L(T nh w)) nh u) + (1 L(T nh)) nh)] du ds g
An,
Vd
[ L(T nh L)) nh - u) + [ L(T - -h)) -h) 12 ds
(L(T nh - u) + I(T -nh)) dus ds Co
2
V2
+
VdLF.
4
Using Lemma with
=(1) = 8, - Inh+t 11.2"
a(I) (An,1 ||nnh
- + An,2) 12,
((1) = I(T - -nh) 1
2'
set to nhh, we then obtain
Yn s(ntaaa I1. & (An,1 || Ynh - O20ll, + An,2) 4n,2))?²xp exp (L((T-nh) n/) - -
-(h)
1e (LCT-nh)- || 100 11, An })n
Invoking display (37). we arrive at
Y/ 550
Yn I(n+1)/ 12 a+1 - I(n+1)all
+ +111051
L2
A AnzelL(r nd) - || - 00 Ili
(L
(L
+ hErr 2 //3/
- 11 60
(M)
e(l(T-nh)--)
Furthermore, we can bound the coefficients An, and An, as follows,
1 An,1 & 8² ('' (.00 2(1 + Txmna) du ds = (1 Temax) ) Immax
An, & (1 + LInxx (2(x) + vāl Lmmax 11/2 + 1
2 22 + V2
vdLr
V22
$ VVl I.max)) L 4
Vdl.p.
Collecting all the pieces then gives
||Ynn
+50 1ol < r. |nnn 150 lll,
+
+050// + hEsc + 2 3/2 + z (
(L((T-nh)
where
50 - (m d) (
c500 An.e (LT-
From the result above, we finally obtain
YNA 050|| 11 - Mmin - 1/2 h max CSO + (Eeo + 2 221/2
+ 1 = el 1316
O
Sh
Vil( 3/2
2 4)elLmax $1h
Exc + 2 Vhdelll
+ 2 ghe cllman f3A
Maiin
1/2
This completes the proof of Theorem 4
46

---

D Proof of Auxiliary Lemma
D.I Proof of Lemma
We have
d H - Gell
d(H, - Gt)
= He - GG:
dt
dt
= 2 H: G.. (H G1) + (Vlogpr-(HH) - Vlogpr- (bog----(()))Vbbooxx-(Gi))
(41)
||1, G.ll 2(11
G., V PoggTT logpr-
< (1 2m 1)) -
The last inequality follows from Lemma 6
Then,
take the derivative of
Sú (2m(T- 1) ||, Gell
d (Li (2m(T-o)-1)ds H Gtll?)
dt
(T - 1)e (m)( |H; G,ll2 + t (m d Grll 00
dt
Therefore, obtain that
el (2m(T= || - G.ll < |H, - Gr, 11.
taking the expectation of both sides and then applying the square root yields the desired result.
D.2 Proof of Lemma 10
By the definition of Yn. we have
1 80ll, = 1
mu+te
( Ye + V log pr pr-nn-e ds
dW,
nh
L
1
Y. di 1/1|/loorT-n (8.)
nh+t
ds
+
dW,
2
nh
To bound the first term, observe that for any [0, h), the following holds
Y Wnn+ 1 + 1| - Ynntolle
Yntts X o lini+alle,, + 1l Ynkts
L2
lo (m( di Ilo X + (nth-(n + + (m f) | 9EM
< 1 (m(T } dt \ Xrlha, + sup
| + +1|Yn
D llt
0
Here, the second inequality follows from the Grönwall inequality applied on Yants xo
and V, - Y nitts and the fact that [Xellı | To bound the second term, we need the following lemma.
Lemma 18. If the target distribution po satisfies Assumption it holds that
|| loogp((X,))l. (2
47

---

According to Lemma 18 it follows that
v log PT uh
log PT nh- (Y,) - v nng-r-nh-s(A phh 1/22
+ log -x
L(T nh
Y, X
+dLL(T nh anh-s))1/2
1
< L(T - nh
Yo 1
+ (dL(T- - -nh - =s)](/2
< L(T nh 8) (1%0 Yn
+ +|\Ym Xinllea) + (dL(T nh s))
112
L, + +e (m((-)))) ||Yo XT + (dL nh s) 1/2
< L(T nh 8) (1 D
Here, we use the fact that Yo- JEM and Grönwall inequality are used in the third inequality and the last one. This completes the proof.
D.3 Proof of Lemma 13
For the stochastic integral of process X, we have
B(L(X))² FL X2
Then, we obtain
Unth
dW, S.'
(fr (latush dW, du
(n+1JA
E
-1 ((0.sa))(000000))sssssro + (or dW, du
nn
p(ntIyM
E
-
+ LU, (o ((0-))kr dW, du
nh
(n+1)
(1tt,))tttt + 1 ( ((ota)hrtttt+000.)) dt du
nh
(E (L(UnEup(u Un)h + LeUn w)hi) du
(u2 u + 1 du
h Jo
3 h.
D.4 Proof of Lemma 16
Define the function w(s) via
w(s) exp (-1,- B(r)dr)) ' 3(r))()) Vs lo
Differentiating this function gives
(()) = (2(s) - S. (())((r) dr) B(s)) ((5) exp (- 3(r) dr < <a(s))(s) exp (-f²² 3 dr)
to
48

---

Note that w(to) = 0. Integrating the function w from to to t yields
w(t) < L a(s))(s) exp (- ((r) dr) ds.
By the definition of w(s), we also have
[ ((s)z()) = (L B(r) dr) w(t).
Combining the previous two displays provides us with
L ((s))()) I a()) (f. B(r) dr ds.
By substituting this estimate into the inequality, a is non-decreasing, then for any s < t, holds that a(s) < a(t). This leads to
obtain the first desired result. Furthermore, if
=(1) & a() a(l) [ ((5) exp (f: ((r)ar)) ds.
which can be simplified to
=(1) & a(1) exp (f" (())rr) 13 l0
This completes the proof.
D.5
Proof of Lemma 17
Notice that
V log = pr(z)) Vp(())
Vep((x)
pr(x)
(logg(())Voog((e)) Pe(x)
V'p((x),
which indicates
d
8
(x)
P
V2 (2p((tt)
Tr
2
Ox?
(r(
2
Pr(x)
1=1
1
Tr
222 V2 log + v V
2
1
Tr
log
|oog
2
Additionally, we have
82 log pc()))
pp)( On
(Qpi(x)
dx?
(x(x)
Or: Pz(x)
(pl((z)
ppp(x))
(logp(())- v ((Slogm(a)) 2)
Ør; )
dx?
(x())
Ox? px(x)
Then, we obtain
Pn((a) Ox?
(x(
(Tr(V2 log (z))) + Tr(Vloogp( + V log
Vlogp((>) V(IV (x)])
49

---

D.6 Proof of Lemma 18
Note that
V po)1 (t(x) dr
lim R-+00 ((O.R)
(Vlogpr(z), V (1)
lim R+000 B(0.R)
(V Volggt()), (pl())) dr,
where (0, R) denotes the Euclidean ball with radius centered at the origin. Using integration by parts, then obtain
2 Olog dS
lim R B(O..))
-ogg))
JOB(O.N) p((x)
Ori
= S p((x) (-A dr < dL((),
where af V 77 represents the directional derivative along the normal vector iž and dS denotes Oi
the surface integral over the spherical surface. Here we use the fact that (z converges to at an exponential rate ||> approaches infinity, and the fact that
log
Tr( log (x)) E [0, dl(()),
which follows from Lemma 6
E Details for second-order acceleration
In this section, we present complete derivation of the second-order acceleration scheme, detailing the implementation of Itô-Taylor expansions and Itô's formula. Building upon the general backward process framework
dæ, = (( - t, 21) dt + dW,
(42)
where a > 0 and W, is the d-dimensional Brownian motion. We apply Itô's formula to >(T t.)..
This procedure generates an approximated structure of SDE
da, = |T(T 8, xx) + Lo xx)) + M. - s)] dt + dW..
(43)
with
07 (T - s, and M, dr
σ² 829 (T s, s,
0n (T 8, x.)
L,
2 O.
Ot
which serves as the foundation for our subsequent second-order discretization.in Section we demonstrate that this approximation preserves the core dynamical structure the original SDE (42) while admitting closed-fore solution. This is achieved by replacing the intractable drift term rT - f, x.) with its Itô-expanded counterpart, which remains analytical tractable through explicit integration.
Applying Itô's formula to e F. yields
d(e (et = =
- (yL s.x.) - L,x, + M,(t s)) dt te L.d L4. dW,
For fixed s, both sides of the equation permit closed-form integration. While the Brownian integral I gst eLat dW, formally appears non-analytic, it equivalently manifests as a Gaussian random variable with explicitly computable variance. This enables full analytical representation for Zs+Ar
50

---

when integrating over (s, 8 At).
stAr
Is+ = claAt xx +
clulstA(-1) d 8,3,)) lIgtas)
+ yotar clalotar-t) - s) M, to rotii elo(04At-1) dW,
Zg Lg (o () s, L - - 11)- L, LM,,,,
rt
+ clalstAt- dW.
Having established the general framework, ization: set y(7 - t,
specialize core through the parameter- 2 VlogPr-(()) a 1, At [0, h) and nh..
In the resulting expression, we denote I, by $ Then for any [nak, (n + 1)h), the solution admits the semi-analytic representation
-55
+
050
V logg Po
(44)
nh
2
Mn nh)) du + f dW
Ln(a 020)
nh
where
In = Ta +
72 Po PT - (0.5)
2
1
82
① V log (2-
M, 2
V log 1
Ot
Though Ln and M, are theoretically defined through exact derivatives in SDE 43 their practical evaluation requires approximations due to the score function's 's computational intractability. We imple- ment these approximations via numerical methods neural networks, with concrete techniques for Ln and M,, estimation provided separately in Appendix and Substituting these approximations into the SDE 44 yields
''' (7(T nh, (05) + 1T - -nh, on (50)) (50)
050
nh
[" dw..
s. MM nh, 1550 (u - nh)) du +
nu
Crucially, this substitution preserves the closed-for form integrability of the original framework. Adopting the same exponential integration strategy above, we derive a closed form of xt. Let 5so, denote 1(++1)/ the second-order discretization scheme given by
0so
000
( nh, 050)
(0 (T nnh - -Ia) iss + sT nh, (20))
1 (T nh, 0500 -2 (e 200m (10) (T nh, 0000 )h 1a) s!M (T nh,, 050)
iiatijn e M dW.
nh Implementation specifics for handling the matrix exponentials and stochastic integral are addressed in Appendix G
Numerical Studies on Synthetic Data
We apply the five schemes to the posterior density of penalized logistic regression, defined by po( x exp(-f(0) with the potential function
> f(0)
2 1011
1 ndata
log(1 + exp(-yar, 0)),
2
51

---

144
6.5
033
0.0
3.59
0.3
Step
Figure 2: Error of various discretization schemes and second-order sampler with different choice of step size.
where > 0 denotes the tuning parameter. The data 3e composed of binary labels iåd NO, 100).
V/ E { 111,1} and features 2; = (Fi,1," ,3,a) Rd generated from Fi.
F. Implementation Details
In the numerical studies, we set T 10, and the number of Monte Carlo iterations is chosen as the floor of T//, where h varies according the step size indicated in the figure. Figure 2 shows the Wasserstein distance measured along the first dimension between the empirical distributions of the N-th outputs from SGMs and the target distribution, with different choices of the step size h. In this simulation, we use the Konte-Carlo method to estimate the score function and the Hessian matrix.
Calculation
In this part, we derive explicit formulas for each coefficient term we need, First, the score function can be computed as
Yiu; exp( (-pt 0)
log)o(0) = 1O +
Tidata
1 exp( - yo. 0)
i=1
Thdata
413s
X0
+
Tddata
+ exp( 0)
i=1
For simplicity, we denote the logistic sigmoid function a(u) = 1 e- then
Tidara
(o)) X0 +
yi (-)xr,)))))))))))
Addtt
Since d')) = a(u)[1 - [(u)]. we have
log =
Pa +
(- 0) 00- 0)] xtar?
Rdata
1=1
XLa
a- (T ()) [1 - o1-oo-yia0 (;x 0)] n.p.
iddta
As 1. > 0. 2 2 log (t) z- -Ad. We also have that a(1 y: () (0, 1), then
2 (0) > Ald
4ridata
i=1
(A +
Amnaxe
24 Ia
Tddata
52

---

Therefore,
mo X, = =+
Amax
Tidata
Recall that the transition probability Pr (r.((0))))) (00;; M4, Ee), where Me = = Ex e LLa, and @00, ue Σ) denotes the probability density function of (u, E), then we have
/
Pr(O,) = Jpc Peio( (0: (006 d0
1
exp(
200. (00 pe)) do
1
J (2
See
exp(
2(1 - -e-1) 110 -e ()ool)))((00)) d0o
[2
[2 e-)))a/2 Edg
exp(
2(1 e-t) 110, c -3388611)
1
Hence,
1
1
||0 - --- -++86613)) )
Vpr(O,)
E Ego-po (2 e- -1]]1/2
cxp(
2(1 - e c-r)
1
2(1 e e- 100 - -(0, -= 1- e-
E8 EBorgo [2x(1 e ()]]d//
exp(
(p))
1
Ed Eoorpo
exp(
2(1 e-t) 110 - ((001l3)
1
[2x(1 =-)))//2
(0, e- +++000 (0+ - = 1-c-t Ta
(1 e)-
We can approximate pt (0:), Vp.(O.) and (O,) even higher order derivative tensor of p (0.) by Monte Carlo method, therefore, we can compute score function and its high order derivative by
Vp.(O.) (0,)
V (og = ('Dr(0.) (P.(O.) (0.)
log (O,;)
(r,) po(0.)²
G Real Data Analysis
G.1 Implementation Details
We set the step size h = 0.2 and N = 2/h. We conducted experiments on NVIDIA RTX 4060 GPU (16GB VRAM). The training process required GPU hours over 100 epochs with a batch size of 32, using CUDA 12.4, PyTorch 2.4, and torchvision 0.20.0. Figure 3 shows the digits generated by five algorithms, using the same score functions. The execution times for the algorithms are as follows: EM method 2 hour 12 min 49 El method 2 hour 12 min 50s, REM method 2 hour 13 min s0s, REI method 2 hour 13 min 47 s, SO method hour 14 min 05s. s.
G.2 Score matching function for second order acceleration
For the MNIST dataset, we have demonstrated in the proof of Proposition 15 that computing third- order derivatives is unnecessary. Unlike existing high-order methods for estimating second-order scores (28) which require the joint training of score functions and Hessian matrices and consequently incur substantial computational overhead, our second-order algorithm avoids explicit computation of the Jacobian matrix. Furthermore, by employing Hessian-vector products (HVPs), we efficiently capture higher-order information, enabling our second-order acceleration method to achieve improved performance with reduced iteration complexity and manageable computational cost.
More specifically, in the experiments of the MNIST dataset, we construct U-Net architecture incor- porating time and label embeddings to train the score function, where the time embedding operates
53

---

EM
EI
REM
REI
50
D
D
0
O
0
Q
O
D
O
O
0
1
1 1
/ 2 2 3 3 3
1
J
L
1
1
/ 2 3 3 4 4
/
1 2
2
2 22
2
2 2 2
2
2 B 4 5
3
2
3 3 3
3 - 3
3
3
3
5 7 4
4/ 4 4
4 144 4 5 5 5 3 6 6
4
9
4 5 6 7 8
5 -
5 5 5 6 666 6
5
5 S 6 7 8
6 6 6 7
6 7
6
6
7 977
7 7 4 8 8 8 9
7 7 7
7
7
7
9 8 8
8 7 3
8
r 8 9
8
9
9
7
9
9
9
9
7
9
9
9
9
Figure 3: Comparative visualization of generated MNIST digits under various discretization schemes.
on the temporal variable t of the score function, while the label embedding leverages MNIST's cate- gorical digit labels. This conditional formulation expresses the score function V 1og enabling per-class score estimation through discriminative embedding propagation.
Recall iteration rule of the SO algorithm, we have
n+1 -0500 (L) (T nhh 050
22 ( ( (Th -11) 1a) pgso + +s nh, 1200
-la)
o (1) (T nh, (55) - nT 5t (T nh, (5O0 ) 1a) s. (T nh, )50)
(nt1jn !0 (T-nh ((x dW.
nh
Note that
(Ih
20 (T ((n+)))) dWW dw, N
1 (4) 7 (T- nh, 050 (222) (T (T-nh.020 -la))
Jnh
Let (.) = (T - nh,,)) denote the score matching function at time T - nh. Although the approximation of the Hessian matrix 2 log (n)((()) will not explicitly appear in the algorithmic implementation, we formally designate it as H.(.)) (.. for notational clarity. Consequently, the estimators of Ln, and Mn are chosen be
(1) (T - nh, 550 = la + H (USO)
2
MM) (T nh, م55 = 2 (n H (U50) yso + (o(0))
Employing the Taylor expansion, we have
8
11 (
(o! nol (c - -Ta)
k!
k=1
8
ht (1
(s(2))-2 (chos") -hs(") -Ia)
2=
54

---

1 (o 2
(2hhs!")
1/2
vh
as aeshs
where (ao: 01, d3, 04, d5, "....) (1, 2: 21 IG 5760 1280. 193036" ..... ) We thus reformulate all operators in the discretization scheme using matrix multiplications, which is crucial step that avoids the explicit computation and storage of the full Hessian matrix H,- This is achieved by leveraging Hessian-vector products (HVPs) via automatic differentiation, which reduces memory complexity to ((d) while retaining second-order curvature information. Specifically, given that In corresponds the neural network's output and Hi, represents its Jacobian matrix, compute H for any vector v through Py PyTorhh's 's reverse- reverse-mode differentiation (torch. autograd, grad). By iteratively applying this HVP procedure & times, efficiently construct H for any k > Through Taylor series expansion, these HVP-powerd computations enable precise evaluation of each term the discretization scheme.
79
55


Reviews:
[{"review_id": "ApMP0nEpVt", "reviewer": "Reviewer_hQEd", "summary": "The authors analyze the effects of the Euler discretization, exponential integrator, midpoint randomization, and a novel exponential integrator with midpoint randomization scheme in the convergence of score-based diffusion models within the Wasserstein distance. They show that midpoint randomization introduces a variance term that cancels out the improved estimation, offering no benefit in convergence. The authors also propose an accelerated sampler that uses second-order information, achieving an improved convergence rate of $O(1/\\epsilon)$ iterations. A numerical study is offered, showing that the accelerated sampler allows for smaller error than the four previous discretizations (among which the randomized midpoint and Euler schemes perform the best).", "questions": "Can this requirement of exact second-order information be relaxed for a $O(1/\\epsilon)$ rate, as is the case for convergence results under the TV distance as in Li et al. (2024) and Li et al. (2025)?\n \nThis is a good paper. I am on the fence between a score of 4 and a score of 5. My concerns regarding whether the accelerated sampler is truly accelerated or not lead me to lean towards a score of 4, as this heavily affects the significancy of the paper. I will be happy to raise my score if the authors alleviate my concerns.", "limitations": "Yes", "rating": 4, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 3, "significance": 2, "originality": 3, "strengths_and_weaknesses": "## Strengths\n**Analysis of various discretization schemes.** \nThe authors study four different discretization schemes, and qualify the differences in performance therein from a theoretical standpoint. This analysis is novel, and welcome to the community, representing a step forward in understanding the behavior of discretization schemes within diffusion model sampling.  \nA nice numerical experiment (for a theory paper) is provided, allowing for an empirical comparison of these discretization schemes in practice that complements the theory provided. \n\n**Well-written.** The paper is very well-written. I commend the authors with respect to the clear amount of effort taken to present the results in a clear and coherent manner. Intuition with respect to the results and methodology is readily provided, and overall enhances the quality of the paper. \n \n## Weaknesses\n**Acceleration requires second-order acceleration.** Previous work suggests that this may not be necessary — Li et al. (2024) achieve a speedup for a deterministic sample in https://arxiv.org/pdf/2403.03852, though their guarantees are in the TV distance.\nIn fact, the guarantee provided by the “accelerated” stochastic sampler of that paper mirrors what I suspect might be the case here. The rate for the “accelerated” stochastic sampler within Li et al. (2024) was shown to be achieved by the vanilla DDPM sampler by another paper written by the same first author (Li et al. (2025)) a year later. It may be possible that the first-order sampler may be able to achieve the $O(1/\\epsilon)$ rate, though as Gao et al. (2025) remark, this may require additional assumptions.\nThe authors may have missed this remark, as seen in their remarks within Appendix A. The authors state that the “accelerated SDE sampler” of Li et al. (2024) “leverages higher-order expansions of the conditional density to enhance efficiency”. Yet, the rate achieved by the stochastic sampler of Li et al. (2024) is worse than the DDPM rate achieved by Li et al. (2025).\nThe authors do not quantify the dimensional dependence within the sample complexity explicitly. This hurts readability and clarity.\n \n## References\nLi et al. (2024), Accelerating Convergence of Score-Based Diffusion Models, Provably\n\nLi et al. (2025), O(d/T) Convergence Theory for Diffusion Probabilistic Models under Minimal Assumptions\n\nGao et al. (2025), Wasserstein Convergence Guarantees for a General Class of Score-Based Generative Models", "was_revised": false, "final_justification": "", "rebuttal": "We appreciate the reviewer's feedback. Please see our detailed responses to each comment below.\n\n1. >`Acceleration requires second-order acceleration. Previous work suggests that this may not be necessary — Li et al. (2024) achieve a speedup for a deterministic sample in https://arxiv.org/pdf/2403.03852, though their guarantees are in the TV distance. In fact, the guarantee provided by the “accelerated” stochastic sampler of that paper mirrors what I suspect might be the case here. The rate for the “accelerated” stochastic sampler within Li et al. (2024) was shown to be achieved by the vanilla DDPM sampler by another paper written by the same first author (Li et al. (2025)) a year later. It may be possible that the first-order sampler may be able to achieve the $O(1/\\epsilon)$\n rate, though as Gao et al. (2025) remark, this may require additional assumptions...Can this requirement of exact second-order information be relaxed for a $O(1/\\epsilon)$ rate?`\n\nWe first consider the result of Li et al. (2025), which establishes a sharper convergence rate for the vanilla DDPM sampler. Their improvement arises from a novel proof technique: instead of relying on the standard inequality ${TV}^2(p_{X_1},p_{Y_1})\\leqslant \\frac{1}{2}{KL}(p_{X_1}\\| p_{Y_1})$, they directly analyze the total variation (TV) distance. This direct control of TV allows for a tighter bound. Specifically, they introduce three auxiliary reverse processes to decompose the total error into estimation error and discretization error, with the improved rate stemming from tighter control of the latter. However, the analysis does not consider second-order samplers under the same proof framework, so it cannot be used to conclude that second-order information is ineffective. In fact, whether similar improvements can be achieved for second-order methods under this framework remains an open and interesting question.\n\nIn our work, within the framework of Wasserstein distance and score-based diffusion models, both second-order information and smoothness assumptions contribute to improving the convergence rate. Our analysis adopts a unified proof skeleton for both first-order and second-order methods, introducing auxiliary processes to decompose and control the discretization error. The incorporation of second-order information plays a clear and essential role in reducing the discretization error in the backward dynamics. This improvement is not an artifact of the proof technique but reflects the numerical behavior of the sampler. Therefore, regardless of how the proof is structured, the benefit of second-order information remains fundamental and cannot be dismissed. Meanwhile, the smoothness assumptions are necessary to control the behavior of the time-dependent score function $s_t(x)$ for $t\\in[0,T]$, which plays a central role in the backward dynamics and differs from the DDPM setting. While the optimal convergence rate in Wasserstein distance remains unknown, our results show that second-order acceleration is indeed effective in reducing the error introduced by discretizing the drift term in the backward process.\n\n\n2. >`The authors may have missed this remark, as seen in their remarks within Appendix A. The authors state that the “accelerated SDE sampler” of Li et al. (2024) “leverages higher-order expansions of the conditional density to enhance efficiency”. Yet, the rate achieved by the stochastic sampler of Li et al. (2024) is worse than the DDPM rate achieved by Li et al. (2025). ` \n\nWe thank the reviewer for pointing out this important clarification. In Appendix A, our use of the term \"accelerated SDE sampler\" for Li et al. (2024) was intended to reflect the authors' design motivation, which involves higher-order expansions of the conditional density, rather than to claim superior convergence guarantees. We agree that, as noted in Li et al. (2025), the stochastic sampler proposed in Li et al. (2024) does not outperform the DDPM rate in terms of theoretical convergence. To avoid confusion, we will revise the wording in Appendix A to more accurately reflect this distinction and acknowledge the comparison raised by the reviewer.\n\n3. `The authors do not quantify the dimensional dependence within the sample complexity explicitly. This hurts readability and clarity.`\n \n We thank the reviewer for raising this point. We have explicitly quantified the dimensional dependence in the sample complexity for Theorems 1, 3, and 4 in our revision. For clarity, we summarize the results below:\n\n| | EM / EI / REM / REI | SO |\n|----| -------- | -------- | \n|Sample Complexity| $\\tilde{\\mathcal{O}}(d/\\varepsilon^2)$ | $\\tilde{\\mathcal O}(d/\\varepsilon)$  |"}, {"review_id": "a7CD37Vw4V", "reviewer": "Reviewer_PXDo", "summary": "This paper established a Wasserstein convergence theory for the score-based diffusion model. It provided upper error bounds for several popular discretization methods: Euler-Maruyama, exponential integrator, and randomized midpoint methods. All discretization methods offer the same rate asymptotically and require $\\tilde{O}(\\epsilon^{-2})$ steps to reach $\\epsilon$ accuracy in $W_2$ distance, though midpoint randomization enables parallel computation and reduces computation time. The author also provides an acceleration method when the second-order information is available, and slashes the total sample steps to $\\tilde{O}(\\epsilon^{-1})$. The theoretical contributions were confirmed by the experiments on synthetic data and the MNIST dataset.", "questions": "1. As you mentioned, some other works have also focused on the convergence of diffusion in $W_2$ distance. Can you make a brief comparison of your result and theirs?\n\n2. As mentioned in Weakness 1, is there a way to drop the Lipshitz continuity of the score function when analyzing convergence in $W_2$ distance?\n\n3. Can you comment on the dimension dependence of the convergence result? Is it $\\widetilde{O}(d/\\epsilon^2)$ for non-accelerated version and $\\widetilde{O}(d/\\epsilon)$ for accelerated version?", "limitations": "Yes", "rating": 4, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 3, "significance": 3, "originality": 3, "strengths_and_weaknesses": "Strengths:\n\n1. The paper delivers a thorough Wasserstein-2 convergence analysis across multiple discretisation schemes, filling a gap left by prior work that centred on KL/TV metrics and relied almost exclusively on Euler steps.\n\n2. The author proposed a second-order acceleration method to improve the convergence rate. This is the first acceleration that only requires $\\tilde{O}(\\epsilon^{-1})$ sample steps to achieve $\\epsilon$ accuracy in $W_2$ distance.\n\n3. Overall, the paper is well written and easy to follow, and the result seems correct.\n\nWeaknesses:\n\n1. The theoretical guarantees hinge on strong assumptions—especially Assumptions 2 and 6—that may not hold for many non-Gaussian or multimodal data distributions. Recent TV-distance analyses (e.g., Li et al., 2024) achieve convergence without such restrictive conditions.\n\n2. The accelerated sampler needs accurate Hessian–vector products of $\\nabla^2\\log p_t(x)$ at every diffusion step. Computing these in high dimensions (e.g., $64\\times64$ ImageNet) may be impossible. Alternative accelerations that approximate Hessians via finite differences of score functions (again Li et al., 2024) to avoid this bottleneck. As a result, I think the assumptions for this paper can be further weakened if this strategy is applied (though Li et al., 2024, focus on the convergence of TV distance).\n\n3. Although the authors mention existing works focus on analyzing the convergence of $W_2$ distance, no comparison was given. It would be better to include a brief comparison of the setting/convergence rate between this paper and previous work.\n\nTypos:\n\n1. References [13] and [14] are the same.\n2. Page 16, Line 493, missing one right bracket \")\" in the first line.\n\nReferences:\n1. Li et al., 2024, \"Accelerating Convergence of Score-Based Diffusion Models, Provably\".", "was_revised": true, "final_justification": "The paper is comprehensive and solid, but there are a few weaknesses, especially in acceleration part, that limits the implement of the method. Therefore. I decide to maintain my original positive score.", "rebuttal": "We thank the reviewer for the feedback. Please find our point by point response below.\n\n1. >`The theoretical guarantees hinge on strong assumptions—especially Assumptions 2 and 6—that may not hold for many non-Gaussian or multimodal data distributions. Recent TV-distance analyses (e.g., Li et al., 2024) achieve convergence without such restrictive conditions...As mentioned in Weakness 1, is there a way to drop the Lipshitz continuity of the score function when analyzing convergence in $W_2$ distance?`\n\nWe thank the reviewer for this insightful comment. We agree that the Lipschitz continuity of the score function (Assumption 2) is a strong condition. However, this assumption is currently necessary to control the discretization error of the continuous-time backward process when analyzing convergence in the $W_2$ distance, as also indicated in [13]. \n\nSimilarly, for the SO scheme, we require Assumption 6, which imposes Lipschitz continuity of the Hessian. The only relaxation we are aware of is given in [28], where the Hessian Lipschitz condition is assumed only at the initial distribution $p_0$, not along the entire path. However, [28] focuses on accelerating DDPM under access to higher-order derivatives and analyzes convergence in KL divergence, relying on tools fundamentally different from those used in our Wasserstein-based analysis.\n\nIndeed, unlike analyses in the TV or KL divergence, convergence in the $W_2$ distance is particularly sensitive to the regularity of the score function due to its intrinsic dependence on the geometry of the space. Nonetheless, we acknowledge the importance of relaxing these assumptions and will add a discussion of this issue and related directions in our revision.\n\n2. >`The accelerated sampler needs accurate Hessian–vector products of $\\nabla^2 \\log p_t(x)$ at every diffusion step. Computing these in high dimensions (e.g., $64\\times 64$ ImageNet) may be impossible. Alternative accelerations that approximate Hessians via finite differences of score functions (again Li et al., 2024) to avoid this bottleneck. As a result, I think the assumptions for this paper can be further weakened if this strategy is applied (though Li et al., 2024, focus on the convergence of TV distance).`\n\nWe thank the reviewer for raising this important point. We agree that computing accurate Hessian–vector products at every diffusion step can be challenging in high-dimensional settings. In our current implementation, we use PyTorch's automatic differentiation to compute Hessian–vector products efficiently without explicitly forming the full Hessian, which makes the approach feasible in moderate dimensions. We also appreciate the reviewer's reference to the finite-difference approximation strategy discussed in Li et al. (2024). That work focuses on the acceleration of DDPM and DDIM samplers and establishes convergence guarantees under the total variation distance by approximating the Jacobian of the score function through finite differences. While their framework and target metric differ from ours, extending such approximation techniques to our stochastic setting and analyzing their effect under Wasserstein-2 convergence is indeed a compelling direction. We view this as promising future work that could potentially help relax some of our current assumptions while maintaining theoretical guarantees.\n\n3. >`Although the authors mention existing works focus on analyzing the convergence of $W_2$ distance, no comparison was given. It would be better to include a brief comparison of the setting/convergence rate between this paper and previous work...As you mentioned, some other works have also focused on the convergence of diffusion in $W_2$ distance. Can you make a brief comparison of your result and theirs?`\n\nWe thank the reviewer for the helpful suggestion. Our work is indeed motivated by the limited understanding of Wasserstein convergence for score-based diffusion models. To our knowledge, the only prior comparable theoretical results in $W_2$ distance are for the Euler–Maruyama (EM) scheme, as presented in [14,38]. This comparison is briefly discussed in line 176 of the main text. While other works also study the Wassertein convergence for EM scheme, they typically do so under assumptions that differ substantially from ours. In contrast, no existing theoretical results are available under settings similar to ours for the EI, REM, REI, or SO schemes, making our work the first to provide such an analysis. We will revise the manuscript to clarify this point and include a more detailed comparison with existing EM-based results to highlight the distinctions in assumptions and convergence behavior.\n\n4. >`Can you comment on the dimension dependence of the convergence result? Is it $\\tilde O(d/\\epsilon^2)$ for non-accelerated version and $O(d/\\epsilon)$ for accelerated version?`\n\nWe have explicitly quantified the dimensional dependence in the sample complexity for Theorems 1, 3, and 4 in our revision. For clarity, we summarize the results below:\n\n| | EM / EI / REM / REI | SO |\n|----| -------- | -------- | \n|Sample Complexity| $\\tilde{\\mathcal{O}}(d/\\varepsilon^2)$ | $\\tilde{\\mathcal O}(d/\\varepsilon)$  |\n\n\n5. Typos: We have corrected the typos in the revised version."}, {"review_id": "3XmLCqfdyR", "reviewer": "Reviewer_cN6S", "summary": "This paper provides a theoretical convergence analysis of diffusion models in the Wasserstein-2 distance.\nSpecifically, the authors first conduct a comprehensive comparison of different discretization schemes, including EM method, exponential integrator (EI) and randomized midpoint methods. \nThese schemes are shown to have a similar convergence rate of $\\mathcal{O}(\\frac{1}{\\epsilon^2} \\log (\\frac{1}{\\epsilon}))$ in $W_2$ distance.\nIn addition to these first first-order discretization results, the authors introduce an accelerated sampler that leverages the Hessian estimation. \nTheoretically, they prove that this second-order discretization scheme achieves a significantly improved convergence rate of $\\mathcal{O}((\\frac{1}{\\epsilon} \\log (\\frac{1}{\\epsilon}))$.\n\nIn the experiments, the proposed second-order accelerated method is compared with the above first-order discretization schemes on both a synthetic posterior of a penalized logistic regression and the MNIST data distribution.\nThe numerical results in Figure 1 confirm that the proposed second-order method converges faster, which aligns with the theoretical analysis.", "questions": "- My major question is about the relationship between the convergence rate and the dimension of the target distribution.\n- In the MNIST experiments, the running times for all five methods are reported to be close (within a 2-minute range for a >2 hour run) in Appendix G.1. \nCould the authors explain why the additional HVP computations in the second-order method did not lead to a noticeable increase in running time.", "limitations": "yes", "rating": 4, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 2, "clarity": 4, "significance": 2, "originality": 3, "strengths_and_weaknesses": "### Strength\n- The paper provides a systematic comparison of multiple popular discretization schemes (EM, EI, REM, REI) in terms of $W_2$ distance.\n- The use of Hessian-vector products in the real-data experiment like MNIST makes the second-order approach computationally feasible.\n\n### Weaknesses\nThe theoretical results are comprehensive, but further analysis could help clarify their properties and underlying assumptions.\n- In contrast to prior work [1], which relies on the relatively mild assumption of a bounded second moment for total variation (TV) divergence analysis, Assumption 1 imposes a stronger condition—strongly log-concavity. This assumption implies the Log-Sobolev inequality, an assumption that is generally not adopted in most recent theoretical results of diffusion models. \n- For the data generation task, the paper would benefit from the inclusion of quantitative metrics. \nTo better align the experiments with the theoretical results, the authors could report the FID or $W_2$ distance on toy 2D tasks, adopting an approach similar to that shown in Table 1 and Table 2 of paper [2]. \n\n[1] Fan, Ying, and Kangwook Lee. \"Optimizing ddpm sampling with shortcut fine-tuning.\" ICML 2023.\n\n[2] Gupta, Shivam, Linda Cai, and Sitan Chen. \"Faster Diffusion Sampling with Randomized Midpoints: Sequential and Parallel.\" ICLR 2025.", "was_revised": true, "final_justification": "I agree with the authors' response regarding the stronger assumptions. Overall, they have also addressed my other concerns. I recommend acceptance of the manuscript.", "rebuttal": "We thank the reviewer for the constructive remarks. Our responses to each point are provided below.\n\n1. >`In contrast to prior work [1], which relies on the relatively mild assumption of a bounded second moment for total variation (TV) divergence analysis, Assumption 1 imposes a stronger condition—strongly log-concavity. This assumption implies the Log-Sobolev inequality, an assumption that is generally not adopted in most recent theoretical results of diffusion models.`\n\nWe thank the reviewer for this thoughtful observation. It is true that Assumption 1 (strong log-concavity) is stronger than the bounded second-moment assumption used in prior TV-based analyses such as [1]. However, our analysis focuses on Wasserstein convergence for score-based diffusion models on non-compact domains, which requires stronger assumptions to ensure stability.\n\nIn particular, our analysis relies on the Wasserstein contraction properties of the reverse SDE to control discretization and score-matching errors. These contraction results critically depend on time-uniform strong convexity and smoothness of the score function $\\nabla \\log p_t$, which we derive from Assumption 1, as also demonstrated in [13].\n\nSome existing works avoid strong convexity by assuming compactly supported data and introducing early stopping with projection onto the support, e.g. [5]. In contrast, our setting allows for unbounded data distributions and thus necessitates a different technical framework.\n\nWe agree that relaxing Assumption 1 is an important and challenging direction for future work. We will revise the manuscript to clarify the role of this assumption and contrast it more explicitly with assumptions used in TV-based analyses.\n\n2. >`For the data generation task, the paper would benefit from the inclusion of quantitative metrics. To better align the experiments with the theoretical results, the authors could report the FID or $W_2$\n distance on toy 2D tasks, adopting an approach similar to that shown in Table 1 and Table 2 of paper [2].`\n\nWe thank the reviewer for the helpful suggestion. Figures 1 and 2 in our paper already report the $W_2$ distance for each scheme across different step sizes, aligning closely with our theoretical analysis. In the revised version, we will include the FID for the MNIST data generation task to provide a more comprehensive evaluation of sample quality.\n\n3. >`My major question is about the relationship between the convergence rate and the dimension of the target distribution.`\n \nWe have explicitly quantified the dimensional dependence in the sample complexity for Theorems 1, 3, and 4 in our revision. For clarity, we summarize the results below:\n\n| | EM / EI / REM / REI | SO |\n|----| -------- | -------- | \n|Sample Complexity| $\\tilde{\\mathcal{O}}(d/\\varepsilon^2)$ | $\\tilde{\\mathcal O}(d/\\varepsilon)$  |\n\n4. >`In the MNIST experiments, the running times for all five methods are reported to be close (within a 2-minute range for a >2 hour run) in Appendix G.1. Could the authors explain why the additional HVP computations in the second-order method did not lead to a noticeable increase in running time.`\n\nWhile second-order methods do introduce additional computational overhead due to HVP computations, in our implementation this overhead remains moderate. Specifically, we use PyTorch’s automatic differentiation to compute HVPs, which requires an additional backward pass to apply the Hessian to a given vector. As a result, the total cost per step (forward + backward) is indeed higher than that of first-order methods, which is both reasonable and expected. Moreover, since our experiments did not use GPU acceleration, the observed overhead may appear more pronounced than it would in a GPU-accelerated setting."}, {"review_id": "RIBexSSi2A", "reviewer": "Reviewer_mdah", "summary": "This work first studies the convergence of various discretization methods for the diffusion backward process, measured by the Wasserstein-2 distance. This work also derives a new discretization method based on the higher-order derivatives.", "questions": "1. Continuing Weakness 2b, What is special about the Wasserstein-2 distance?\n2. In all theoretical results, a lot of coefficients are related to $m\\_{\\min}-0.5$. What is special about the $-0.5$ term? Does it pose an implicit restriction on $m_{\\min}$ that $m_{\\min}>0.5$?", "limitations": "As is discussed in \"Weaknesses\", the experiments do not seem to be adequate. The paper does not seem to have any negative societal impact.", "rating": 4, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 2, "clarity": 3, "significance": 2, "originality": 2, "strengths_and_weaknesses": "## Strengths\n\n1. The theoretical results are presented in good logic, and the new discretization method is also derived with sound reasoning.\n2. The assumptions are clearly stated and moderate.\n\n## Weaknesses\n\n1. As a whole, this work looks like a simple combination of two parts: one about the analysis of Wasserstein-2 distance, and the derivation of the new discretization. These two parts seem to be only loosely connected.\n2. For the first part, the work seems to lack **detailed** comparison against previous works.\n- 2a. Comparison against previous works on Wasserstein-2 distance: In the last paragraph of Section 3.1, results in [14] are mentioned. How do the results in this work compare against those in [14]?\n- 2b. The comparison against results based on KL and TV are mentioned but not explained in detail. It is unclear whether the results about Wasserstein-2 distance imply different insights.\n- 2c. The authors could also discuss the results in Li et al. (2024).\n\nLi, Runjia, Di, Qiwei, and Gu, Quanquan. Unified Convergence Analysis for Score-Based Diffusion Models with Deterministic Samplers. 2024\n \n3. For the second part, the experiments do not seem to be adequate.\n- 3a. The experiments are only conducted on MNIST, which is a relatively easy task. I would suggest that the authors include more challenging tasks so that the possible instability in the estimation of higher-order derivatives can be evaluated.\n- 3b. I would suggest that the authors include baselines based on ODE sampling, especially the Runge-Kutta discretization.\n- 3c. A possible drawback of the algorithm is that the estimation of high-order derivatives may cause higher computational complexity. However, the results does", "was_revised": true, "final_justification": "I find the contributions of this work above the threshold of acceptance.\n\n1. The first part about the characterization of the sampling error measured by W2 distance is new. Although the assumptions are quite strong, it is a minor drawback given that there have been few other works studying the W2 distance.\n2. The proposed new method is an effective method that leverages the information of the hessian of $\\log p$. The hessian-vector product avoids introducing significant computational overhead, and the empirical results are good.\n3. The main drawback of the current version is that some discussions that I find important, including why the study of W2 distance is important, comparison with previous theoretical works (Weakness 2a), and comparison with other methods that use second-order information is inadequate. I am looking forward to improvements in the revision.\n\nGiven the strengths that overweigh weaknesses, I would like to recommend accepting this paper.", "rebuttal": "We thank the reviewer for the feedback. Our responses to each point are provided below.\n\n1. >`As a whole, this work looks like a simple combination of two parts: one about the analysis of Wasserstein-2 distance, and the derivation of the new discretization. These two parts seem to be only loosely connected.`\n\nWhile the paper may appear to have two separate components, both sections are closely connected through the shared goal of analyzing score-based diffusion models under the Wasserstein-2 ($W_2$) distance. The first part studies first-order discretization schemes and provides convergence guarantees in $W_2$. The second part builds on this by analyzing second-order methods, which to our knowledge have not been well explored in the context of Wasserstein convergence. The development of the second-order scheme leverages techniques and insights from the first-order analysis. Together, both parts contribute to a unified theoretical understanding of the Wasserstein behavior of discretized diffusion samplers.\n\n2. >`Comparison against previous works on Wasserstein-2 distance: In the last paragraph of Section 3.1, results in [14] are mentioned. How do the results in this work compare against those in [14]?`\n\nThe results in [14] focus on score-based diffusion models using the Euler discretization scheme and provide $W_2$ convergence guarantees under that setting. In our work, we compare our results with those of [14] immediately following Corollary 2 (line 176), highlighting the consistency between the results. \n\n3. >`The comparison against results based on KL and TV are mentioned but not explained in detail. It is unclear whether the results about Wasserstein-2 distance imply different insights...What is special about the Wasserstein-2 distance?`\n\nOur work specifically focuses on convergence analysis under the $W_2$ distance, which is fundamentally different from KL and TV divergence in both its geometric structure and the techniques required for analysis. While KL and TV are widely used in prior work, they are not directly comparable to $W_2$, and convergence in one does not generally imply convergence in the others. For this reason, our analysis is tailored to the Wasserstein setting, and a direct comparison with KL or TV-based results may not be appropriate due to the fundamental differences between these metrics.\n\n4. >`The authors could also discuss the results in Li et al. (2024). Li, Runjia, Di, Qiwei, and Gu, Quanquan. Unified Convergence Analysis for Score-Based Diffusion Models with Deterministic Samplers. 2024`\n\nWe thank the reviewer for pointing out the recent work by Li et al. (2024), which provides a unified convergence analysis for deterministic samplers under the total variation distance. In contrast, our work studies stochastic samplers based on discretized SDEs and provides convergence guarantees under the Wasserstein-2 metric. The theoretical settings, assumptions, and proof techniques are thus quite different. While both works share the broader goal of understanding diffusion-based generative models, a direct comparison may not be appropriate due to these fundamental differences. We will add a brief mention of this work in the related literature to clarify the distinction.\n\n5. >`The experiments are only conducted on MNIST, which is a relatively easy task. I would suggest that the authors include more challenging tasks so that the possible instability in the estimation of higher-order derivatives can be evaluated.`\n\nThe numerical studies and MNIST-based experiments in our current submission are designed to illustrate the theoretical findings and validate the convergence behavior provided by our analysis. Since this work is submitted under the theory track, our primary goal is to establish and support theoretical guarantees rather than to optimize empirical performance on complex datasets. While we acknowledge that evaluating higher-order methods on more challenging datasets may be of practical interest, we believe such extensions are better suited for follow-up work focused on empirical exploration. \n\n6. >`I would suggest that the authors include baselines based on ODE sampling, especially the Runge-Kutta discretization.`\n\nOur work focuses on stochastic samplers derived from SDE discretizations, with theoretical guarantees established in the Wasserstein-2 distance. In contrast, ODE-based methods correspond to deterministic samplers and are typically analyzed under different metrics, such as total variation or KL divergence. Since their sampling dynamics and theoretical frameworks differ fundamentally from ours, including them as baselines would fall outside the scope of our current analysis.\n\n\n7. >`A possible drawback of the algorithm is that the estimation of high-order derivatives may cause higher computational complexity.`\n\nWe thank the reviewer for raising this point. While high-order derivative computations can be computationally intensive in general, our implementation alleviates this issue by using Hessian–vector products (HVPs), which avoid explicitly forming or storing the full Hessian matrix. This strategy substantially reduces both time and memory costs, making the method practical in our experimental setting. Nevertheless, developing more efficient and accurate techniques for estimating higher-order derivatives of the score function remains an interesting and valuable direction for future work.\n\n\n8. >`In all theoretical results, a lot of coefficients are related to $m_{min}-0.5$. What is special about the $-0.5$ term? Does it pose an implicit restriction on $m_{min}$ that $m_{min}>0.5$?`\n\nThe $-0.5$ term arises from the drift term of the forward SDE, as demonstrated in Lemma 8. Specifically, the expression $m(T - s) - \\frac{1}{2}$ follows from Line 723, where we use the fact that\n\n$$\n\\begin{aligned}\n\\left\\langle X_t, \\frac{{\\rm d} X_t}{{\\rm d} t} \\right\\rangle = \\left\\langle X_t, \\frac{1}{2} X_t + \\nabla \\log p_{T - t}(X_t) \\right\\rangle \n= \\frac{1}{2} \\\\|X_t\\\\|^2 + \\left\\langle X_t, \\nabla \\log p_{T - t}(X_t) \\right\\rangle \n\\leq \\left(\\frac{1}{2} - m(T - t)\\right) \\\\|X_t\\\\|^2,\n\\end{aligned}\n$$\n\nwhich reflects the contribution of the drift term $-\\frac{1}{2} X_t$ in the forward SDE.\n\nMore generally, if the forward SDE takes the form\n\n$$\n{\\rm d} X_t = -\\beta X_t  {\\rm d} t + {\\rm d} B_t,\n$$\n\nthen the bound becomes $m(t) - \\beta$, reflecting a dependence on the coefficient $\\beta$. While this setup still requires the data distribution $p_{\\textsf{data}}$ to be strongly log-concave, the condition $m_0 > \\frac{1}{2}$ can be generalized to $m_0 > \\beta$ for any $\\beta > 0$."}, {"review_id": "CXeVkGqLww", "reviewer": "Reviewer_mBSW", "summary": "The paper studies the Wasserstein convergence of four different schemes for score-based diffusion models: Euler-Maruyama, exponential integrator, vanilla midpoint randomization, and exponential integrator with midpoint randomization. The paper observes that those two randomized midpoint methods do not give an improved convergence result because of the variance. The paper also proposed an accelerated sampler based on the Local Linearization Method. An approximation to the drift is designed according to Ito's formula, and Hessian estimation of the log density is employed. This leads to an improved dependence on the discretization $h$, so that only $O(\\frac{1}{\\varepsilon}\\log(\\frac{1}{\\varepsilon}))$ is required to be $\\varepsilon$ close to the target distribution.", "questions": "Questions:\n- The theorems involve $m_{\\min} - 1/2$, which is supposed to be positive number. I understand that $m_{\\min} > 1/2$ is a simple condition for one-step contraction in the coupling argument. But I want to ask what if $m_{\\min} > 1/2$ is not satisfied. Do you need to consider other reverse SDEs by changing the coefficient in front of the score in eq.(3) so that the drift in (3) will lead to contraction during coupling? Currently there is no remark on the fact that $m_{\\min} > 1/2$ and I think maybe comments on this in the paper will be helpful to readers.\n- Currently the complexity result (eg. Corollary 2 and Corollary 5) only shows the dependence on $\\varepsilon$. Maybe the dependence on the dimension $d$ is also important especially for high dimensions. Ideally dependence on the condition number $L/m$ (or other counterparts such as $L(t)/m(t)$) is also interesting to people.\n- For $\\lVert X_0\\rVert_{\\mathbb{L}_2}$ in the bounds, if we assume $\\log p_0$ is minimized at $x=0$, then it has the upper bound $\\sqrt{d/m_0}$ (see for example, Lemma 4.0.1 in Sinho Chewi's book draft https://chewisinho.github.io/main.pdf). I want to ask why the authors left the term in the theorem rather than gave it a bound. I am also not sure whether such considerations could simplify/improve some bounds in the proof. (will the minimizer stays at $x=0$ for $p_t(x)$?)\n- How tight are the bounds? For example, the authors can illustrate this when the target is a Gaussian.\n- In Appendix G.2, implementation details are illustrated for the proposed scheme. How do you choose $k$ for the truncation of matrix exponentials?\n\n\nTypos:\n- page 8 footnote 2: the derivation is in Appendix E rather than F. Also, at the end of Appendix E, the matrix exponentials are addressed in Appendix G rather than F.\n- Appendix D.1 is the proof of Lemma 8 rather than Lemma 7.\n- References [13] and [14] are the same.", "limitations": "yes", "rating": 5, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 2, "significance": 3, "originality": 2, "strengths_and_weaknesses": "Strengths:\n1. The paper uses a clear coupling framework to study the convergence of four different practical schemes in a Wasserstein distance. \n2. An improved scheme using Hessian information is proposed, and the experiment on MNIST shows the effectiveness of the new scheme.\n\nWeaknesses:\n1. The analysis is kind of standard: using coupling technique and analyzing the contraction rate and the one-step discretization error. It is kind of incremental, but such a paper to illustrate and compare different discretization schemes could be helpful for people to refer to.\n2. I am not sure how practical the newly proposed scheme is. For example, there is no experiment on CIFAR-10. I understand that this is a theory-focused paper, but I just want to know whether there is some implementation difficulty on CIFAR-10.", "was_revised": true, "final_justification": "The authors have clearly answered all my questions. And the connection to existing works mentioned in other reviews and rebuttals helps me to better understand the current situation for this research direction. So I think the paper can be interesting in this field for future research.", "rebuttal": "We're grateful for the reviewer's thoughtful comments. Please find our detailed responses below.\n\n\n1. >`The analysis is kind of standard: using coupling technique and analyzing the contraction rate and the one-step discretization error. It is kind of incremental, but such a paper to illustrate and compare different discretization schemes could be helpful for people to refer to.`\n\nWe agree that the coupling-based proof technique used in our analysis of first-order samplers follows a well established structure. However, to our knowledge, a systematic comparison of different discretization schemes under the Wasserstein metric in constrained settings has not been previously studied. Applying these tools in this context provides new insights and clarifies the trade-offs between different methods.\n\nMoreover, the second-order sampler we propose, which is based on a local linearization approach, is new within the diffusion model literature. Although some elements of the analysis share similarities with first-order methods, the algorithm itself and the convergence guarantees we establish are original.\n\n2. >`I am not sure how practical the newly proposed scheme is. For example, there is no experiment on CIFAR-10. I understand that this is a theory-focused paper, but I just want to know whether there is some implementation difficulty on CIFAR-10.`\n\n The proposed scheme does not present any fundamental implementation difficulty on datasets such as CIFAR-10. The main challenge lies in computational overhead, primarily due to the increased cost of evaluating second-order information. However, with GPU acceleration and efficient techniques such as Hessian vector products, we believe the method remains feasible and can be implemented within a reasonable runtime. Exploring such high dimensional applications is an interesting direction for future empirical validation, building on the current theoretical foundation.\n\n3. >`The theorems involve $m_min-1/2$, which is supposed to be positive number. I understand that \n is a simple condition for one-step contraction in the coupling argument. But I want to ask what if $m_min>1/2$ is not satisfied. Do you need to consider other reverse SDEs by changing the coefficient in front of the score in eq.(3) so that the drift in (3) will lead to contraction during coupling? Currently there is no remark on the fact that $m_min>1/2$ and I think maybe comments on this in the paper will be helpful to readers.`\n\nWe thank the reviewer for the helpful suggestion. The $-0.5$ term arises from the drift of the forward SDE, as shown in Lemma 8. In particular, the expression $m(T - s) - \\frac{1}{2}$ follows from Line 723, where we use the identity\n\n$$\n\\begin{aligned}\n\\left\\langle X_t, \\frac{{\\rm d} X_t}{{\\rm d} t} \\right\\rangle = \\left\\langle X_t, \\frac{1}{2} X_t + \\nabla \\log p_{T - t}(X_t) \\right\\rangle \n= \\frac{1}{2} \\\\|X_t\\\\|^2 + \\left\\langle X_t, \\nabla \\log p_{T - t}(X_t) \\right\\rangle \n\\leq \\left(\\frac{1}{2} - m(T - t)\\right) \\\\|X_t\\\\|^2,\n\\end{aligned}\n$$\n\nwhich reflects the influence of the drift term $-\\frac{1}{2} X_t$ in the forward process.\n\nMore generally, if the forward SDE is written as\n\n$$\n{\\rm d} X_t = -\\beta X_t {\\rm d} t + {\\rm d} B_t,\n$$\n\nthen the bound becomes $m(t) - \\beta$, showing explicit dependence on the coefficient $\\beta$. While this still requires $p_{\\textsf{data}}$ to be strongly log-concave, the condition $m_0 > \\frac{1}{2}$ can be generalized to $m_0 > \\beta$ for any $\\beta > 0$.\n\nIf the condition $m_{\\min} > 1/2$ is not met, one can adjust the drift accordingly (e.g., using a smaller $\\beta$), which does not affect the main results of the paper. We will add a remark clarifying this point in the revised manuscript.\n\n4. >`Currently the complexity result (eg. Corollary 2 and Corollary 5) only shows the dependence on $\\epsilon$. Maybe the dependence on the dimension $d$\n is also important especially for high dimensions. Ideally dependence on the condition number $L/M$ (or other counterparts such as $L(t)/m(t)$) is also interesting to people.`\n\nWe thank the reviewer for raising this point. We have explicitly quantified the dimensional dependence in the sample complexity for Theorems 1, 3, and 4 in our revision. For clarity, we summarize the results below:\n\n| | EM / EI / REM / REI | SO |\n|----| -------- | -------- | \n|Sample Complexity| $\\tilde{\\mathcal{O}}(d/\\varepsilon^2)$ | $\\tilde{\\mathcal O}(d/\\varepsilon)$  |\n\nRegarding the condition number, Theorems 1, 3, and 4 involve constants that depend on $L_{\\max}$ and $m_{\\min}$. In particular, the term $\\frac{L_{\\max} + 1/2}{m_{\\min} - 1/2}$ can be interpreted as an analogue of the condition number when $m_{\\min}$ is not too close to $1/2$. We will clarify this interpretation in the revised manuscript.\n\n5. >`For $\\|X_0\\|_{\\mathbb{L}_2}$ in the bounds, if we assume $\\log p_0$ is minimized at $x=0$, then it has the upper bound $\\sqrt{d/m_0}$. I want to ask why the authors left the term in the theorem rather than gave it a bound. I am also not sure whether such considerations could simplify/improve some bounds in the proof. (will the minimizer stays at $x=0$ for $p_t(x)$?)`\n\nWe thank the reviewer for the thoughtful question. The term $\\lVert X_0 \\rVert_{\\mathbb{L}_2}$ appears in our bounds as it reflects the dependence on the initialization and arises naturally in the standard decomposition of the Wasserstein error (see, for example, Theorem 5 in [13] and Theorem 2 in [12]). \nWe chose to keep the term explicit in our theorems to maintain generality and avoid imposing potentially unnecessary assumptions on the form of $p_0$.\n\nAlthough assuming that $\\log p_0$ is minimized at $x = 0$ enables the upper bound $\\sqrt{d/m_0}$, this condition does not generally hold for $\\log p_t(x)$ as time $t$ evolves. The forward SDE can shift the distribution, and without strong symmetry or centering assumptions, the mode of $p_t$ may move away from the origin. As a result, relying on such assumptions would reduce the generality of our analysis and may not yield significant simplifications in the proof. Therefore, we opt to make the dependence on $\\lVert X_0 \\rVert_{\\mathbb{L}_2}$ explicit in our bounds.\n\n\n6. >`How tight are the bounds? For example, the authors can illustrate this when the target is a Gaussian.`\n\nWe thank the reviewer for the insightful question regarding the tightness of our bounds. To our knowledge, the only existing lower bound result in the Wasserstein distance is provided in [14], where the target distribution is assumed to be Gaussian. In this setting, the lower bound for the iteration complexity is shown to be of order $\\Omega(\\sqrt{d}/\\varepsilon)$. As noted in [14], the same EM scheme is analyzed, and an upper bound with a convergence rate comparable to ours is established. This bound is also believed to be tight under the same assumptions (see [14] for further discussion). This suggests that the current upper bound may not be further improved without imposing additional structural conditions on the data distribution. Investigating whether such assumptions on $p_0$ could lead to improved upper bounds that match the lower bound in [14] remains an interesting direction for future research.\n\n7. >`In Appendix G.2, implementation details are illustrated for the proposed scheme. How do you choose $k$ for the truncation of matrix exponentials?`\n\nWe thank the reviewer for this helpful question. The choice of $k$ reflects a trade-off between computational cost and approximation accuracy. In our implementation, we use the Taylor expansion to approximate functions of $s_*^{(L)}$, but selecting an optimal $k$ is nontrivial, as we have limited knowledge of the error introduced by estimating $\\nabla^2 \\log p_{T - t}(x)$ using PyTorch’s autograd.\n\nHowever, we can still provide some guidance based on the Taylor expansion residual:\n\n$$\n\\\\|R_n(A)\\\\| \\leq \\frac{\\|A\\|^{n+1}}{(n+1)!} \\sup_{t \\in [0,1]} \\\\|f^{(n+1)}(A)\\\\|.\n$$\n\nFor sufficiently smooth functions $f$, the residual decays roughly as $1 / (n+1)!$, indicating that higher-order terms contribute less.\n\nIn practice, however, the overall score matching error $\\varepsilon_{\\text{sc}}$ is not arbitrarily small. It is lower bounded by the expressiveness of the neural network. In our U-Net training, the DSM loss plateaus around 0.01 after 50 epochs. This suggests that the total approximation error is dominated by training error rather than the truncation error from the matrix exponential.\n\nTo balance efficiency and performance, we tested three values of $k$: 3, 5, and 10. We ultimately chose $k = 5$, as it performed similarly to $k = 10$ but significantly better than $k = 3$, offering a trade-off between accuracy and runtime.\n\n8. Typos: We have corrected the typos in the revised version."}]

Author rebuttal:

Meta-review: This paper advances the theoretical understanding of score-based diffusion models by providing the first systematic Wasserstein convergence analysis across multiple discretization schemes (Euler, exponential integrator, randomized midpoint variants). A key contribution is the introduction of a second-order accelerated sampler leveraging Hessian information, shown to improve the convergence rate. Reviewers found the analysis rigorous and well written, and appreciated both the theoretical novelty and the supporting experiments on synthetic data and MNIST. While concerns remain regarding strong assumptions and practical feasibility of second-order methods, all reviewers ultimately leaned positive, and the consensus is that the paper makes a valuable theoretical contribution. I recommend acceptance as a poster.
Decision: Accept (poster)