Exploring the Optimal Choice for Generative Processes in Diffusion Models: Ordinary vs Stochastic Differential Equations
Authors: Yu Cao
Jingrun Chen
Yixin Luo
Xiang ZHOU
Venue: NeurIPS 2023 poster (2023)
Keywords: diffusion models; stochastic differential equations; score-based generative models; asymptotic analysis
TL;DR: We prove the guaranteed choice of diffusion coefficients in generative diffusion process, by estimating the change of KL divergence with respect to the error in learning scores.

Abstract: The diffusion model has shown remarkable success in computer vision, but it remains unclear whether the ODE-based probability flow or the SDE-based diffusion model is more superior and under what circumstances. Comparing the two is challenging due to dependencies on data distributions, score training, and other numerical issues. In this paper, we study the problem mathematically for two limiting scenarios: the zero diffusion (ODE) case and the large diffusion case. We first introduce a pulse-shape error to perturb the score function and analyze error accumulation of sampling quality, followed by a thorough analysis for generalization to arbitrary error. Our findings indicate that when the perturbation occurs at the end of the generative process, the ODE model outperforms the SDE model with a large diffusion coefficient. However, when the perturbation occurs earlier, the SDE model outperforms the ODE model, and we demonstrate that the error of sample generation due to such a pulse-shape perturbation is exponentially suppressed as the diffusion term's magnitude increases to infinity. Numerical validation of this phenomenon is provided using Gaussian, Gaussian mixture, and Swiss roll distribution, as well as realistic datasets like MNIST and CIFAR-10.

Exploring the Optimal Choice for Generative Processes in Diffusion Models: Ordinary vs Stochastic Differential Equations
Yu Cao
Jingrun Chen
Yixin Luo
Xiang Zhou $
Abstract
The diffusion model has shown remarkable success in computer vision, but it re- mains unclear whether the ODE-based probability flow or the SDE-based diffusion model is more superior and under what circumstances. Comparing the two is challenging due dependencies on data distributions, score training, and other numerical issues. In this paper, we study the problem mathematically for two limiting scenarios: the zero diffusion (ODE) case and the large diffusion case. We first introduce a pulse-shape error to perturb the score function and analyze error ac- cumulation sampling quality, followed by a thorough analysis for generalization to arbitrary error. Our findings indicate that when the perturbation occurs at the end of the generative process, the ODE model outperforms SDE model with a large diffusion coefficient. However, when the perturbation occurs earlier, the SDE model outperforms the ODE model, and we demonstrate that the error sample generation due to such pulse-shape perturbation is exponentially suppressed as the diffusion term's 's magnitude increases to infinity. Numerical validation this phenomenon is provided using Gaussian, Gaussian mixture, and Swiss roll distribution, as well as realistic datasets like MNIST and CIFAR-10.
Introduction
Diffusion models have achieved remarkable success in various artificial intelligence context generation tasks, particularly in computer vision [13]. This technique is rapidly evolving with industrial-level products like DALL-E series. The diffusion model first proposed and studied by Sohl-Dickse -ic-seictt et al. [26] in 2015. Later, Song and Ermon [29] proposed score matching with Langevin dynamics (SMLD) and Ho al. [16] further explored the Denoising Diffusion Probabilistic Models (DDPM). Both formalisms can be interpreted as time-discretization of stochastic differential equations (SDFs) [30]. Since the publication of these seminal works, many techniques have been proposed to improve the efficiency and accuracy of diffusion models, such as DDIM [28]. Analyti-DPPM [4], gDDIM [36], EDM [19], and consistency model [31], among others.
The score- score-based diffusion model involves two steps [30, 17J. Firstly, one estimates the score function, which is the gradient of the logarithm of the probability density function, in the form of neural network. This step uses trajectories of an Seastiin-Uhlebeckk (OU) process starting with given data
Institute of Natural Sciences and School of Mathematical Sciences, Shanghai Jiao Tong University, Shanghai 200240, China, pocao@ejtu.ert edu. en Uniiersstty Science and Technology of China, Hefei 230026. China: Suzhou Institute Advanced Research, University of Science and Technology China, Suzhou 215123, China, jww.grnnnnnnnnusste.... edu. cn
China; Suzhou Institute of Advanced Research, University of Science and Technology of China, Suzhou 215123, China, seeinggmail. ustc. edu.
University of Science and Technology of China, Hefei 230026,
SScoool of Data Science and Department of Mathematics, City University of Hong Kong, Kowloon, Hong Kong SAR, xizhou@cityu. edu.
37th Conference on Neural Information Processing Systems (NcurIPS 2023).

---

samples. This process of injecting noise into structured data is usually referred to as the inference process. Secondly. new samples are generated by simulating a time-reversed SDE, with drift term depending on the learned score function from the first step. This is known as the generative process.
In general, there are two diffusion coefficients g in the inference process and h in the generative process; see $ 2. Regardless of the choice of h and it is always possible to design an SDE in the generative process that matches the forward inference process in the weak sense, i.c., the probability density functions match for both processes. We highlight that this function h (unlike g) does not only appear in the diffusion term, but also enters the drift term in the generative SDE. The choice of g is equivalent to time re-scaling (see Appx. B. B.2), while the choice of h is an important topic in practice. Two common choices of are Probability Flow h 0 [9. 32]. which refers as an ODE, and an SDE-basee diffusion model with h = g [16, 29, 30J. When the score training is accurate, the choice of this function h does not affect the sample generation quality in the continuous- time setting.
In practice, numerical error is inevitable during training the score function. Recent theoretical works [10, 7] have shown that the sample generation quality are affected by three aspects: (1) the truncation of simulation time to a finite T;; (2) the inexact score function; (3) the time-discretization error. The first error is not significantly since the forward OU process converges the equilibrium Gaussian measure exponentially fast in T. The third error can be reduced systematically by more efficient numerical schemes [21], such as exponential integrator proposed in [35, 36]. The inexact training of score function has a few important but subtle consequences. Recent works [10. 7] analyzed the convergence rate of diffusion models, provided that the score training error is sufficient small. However, once the score training is not accurate, the nice equivalence of the generated distribution free of the generative diffusion coefficient h no longer holds as in the idealized situation of exact score function. This raises a key question of our interest about how the choice of h can affect the sample quality in the face of the inexact score training error. Qualitatively, there are two distinctive cases: h- = O or h is large. An important question to ask is: in the presence of non-negligible score training errors, which will produce better sampling quality? Is it the probability flow (h- = 0) or the SDE? More quantitatively, what magnitude of h is optimal?
Related works The impact of h on the generative process seems not yet fully investigated in recent literature, as most experiments used the default choice of this parameter. However, some authors have reported related empirical observations. For example, Song et al. [30] empirically observed that the choice of h = 9 produces better sample generation quality than the ODE case (h = 0) with real datasets. On the other hand, Denoising Diffusion Implicit Models (DDIM) in [28] includes both deterministic and stochastic samplers and points out that the probability flow (h = 0) can produce better samples with improved numerical schemes for the generative process. [36] generalized the DDIM and tried to explain the advantages of a deterministic sampling scheme over the stochastic one for fast sampling. Moreover, Karras et. al., [19] had empirically searched for optimal coefficients which had shown to bring practical advantages. None of these empirical results delivered comprehensive investigations on the influence of the diffusion coefficient, and a consistent and affirmative answer to question still awaits. Recently, there has been rapid progress in theoretical works on error analysis for diffusion models, as seen in [10, 7] and references therein. However, these analyses usually assume specific settings of h, such as h 9. Furthermore, it seems that directly analyzing upper bounds based on these error estimations cannot provide adequate information about choosing the optimal h; see Appx. B.4. Albergo al. proposed a unified framework known as stochastic interpolants and slightly discussed the optimal choice between the probability flow and diffusion models [1, Sec. 2.4]. It is interesting to see how our theoretical analysis below can generalize to their promising unified settings [1].
Our approach To investigate the effect of the diffusion coefficient h on sampling quality, we adopt the continuous-time framework, which precludes time discretization errors. We measure sample quality by the KL divergence between the data distribution po and the distribution of the generative SDE at the terminal time T. Given the assumption that the score function carries numerical errors, consider h as controller and aim to minimize the KL divergence with respect to h. While the optimization problem is straightforward to set up, it is challenging to draw valuable theoretical insights in general setting of approximate score functions. Therefore, choose the asymptotic approach, assuming that the error from the training score is reasonably small with a magnitude of c. Under this assumption, the leading-order term of the KL divergence takes the form
error of sample generation in KL divergence L(h) 2² + O(2)..

---

This e² order is known in [7. 10), but the dependence of this Gateaux differential L(h) on h and other factors has yet to be understood at all. Our contribution amalyze how L(h) behaves as varies; in particular, by considering the constant / in two limiting situations: h 0 and h > 1.
Main Contributions
We summarize main contributions below:
We prove that when the error in score function approximation is a time-localized function only at the beginning of the inference step (i.e., at the end of the generative process), the ODE case (h = 0) outperforms the SDE case (h + oc); see Prop. 3.5, If this (time-localized) error occurs in the middle, then the SDE case has an exponentially smaller error than the ODE case (h 0), as h + 20 (see Prop. 3.4). See Appx. E.4 for reasons behind the time-localized choice.
For general score training error, we prove that h oc, the leading-order term L(h) above converges exponentially fast to a constant, which only depends on the distribution po and the score training error the end of the generative process; see Prop. 3.6. The conclusion about the optimal h depends on how the score training error is distributed over the time horizon [0, T7.
Numerically, we validate the above phenomenon for ID Gaussian, 2D Gaussian mixture, and Swiss roll distribution, as well as realistic datasets like MNIST and CIFAR-10. Due to the tight connection between the distribution of score training error and h, our results may suggest backwardly modifying loss functions during training to adapt to a particular diffusion coefficient of interest. This is topic of independent interest and we report some preliminary experiments in Appx. to validate potential applications of our theoretical analysis. A comprehensive investigation will be left as future works.
Notation convention The time duration T' > 0 is fixed parameter. For any time-dependent function (t, 2) + fe(x), where x E Rd and t E [o. T), we denote f (z) = fr-((z). Sometimes we directly use the "arrowed" variables ft (x) to highlight the direction of time is from reference noise to the data distribution (i.e., the generative direction) even without referring to f first. The notation / 9 means that ʃ S eg where e + in a certain limit, i.e lim sup 5/9 S 1; f f~ g means lim()//) ===- 1. The asymptotic parameter will be explained below explicitly. When two matrices A 3 B, it means B A is positive semidefinite. Ia is the d-dimensional identity matrix; T means an indicator function of a set A; Id is the identity operator. For a random variable X, law(X) means the distribution of X. Some important quantities are summarized in Appx. A.
2 Background
Score-based generative models Suppose we have a collection of data from an unknown distribution with density po- we can inject noise into data via the following SDEs:
dX, = (X,) dt + g: dW,, law(Xo) = po:
(I)
where the drift coefficient ())()): Rd R RRd is dime-dependent vector field, and the diffusion
coefficient 9(-) I + R. is a scalar-valued function (for simplicity). A widely used example is variance-
SDE (VP-SDE) with fe(z) 1/2 gi and g > 0 is typically chosen non-decreasing function in literature [30]. Without loss of generality, we can assume g: = 1 since for any non-zero g. its effect is simply to re-scale the time; see Appx. B.2.
Denote the probability density of X, as Pt and the score function is defined as V log Pt. One main innovation in diffusion models is to find "backward" SDE Y, such that Y, drives the state with distribution pr back to Po- We adopt the arrow of time in this backward direction now and write Y, as
dY. = A (Yr) dt + hit dWe,
law(Y)) = PT,
(2)
where hy is an arbitrary real-valued function of time. The distribution of Y, is denoted as qt. Provided that the score function is available, we can select A such that q: is the same as pr- t, in particular, qr = Po. It is not hard to derive that we can ensure q Pr- if we choose
(gt ja (1 log Pt (1)).
(x)= = (a)
(3)
2
A self-contained proof is provided in Appx. B. When h- = n², it refers to the backward SDE used [30]; when ht 0, it refers to the probability flow therein. More general interpolation of diffusion and flow can be found in, e.g., [1].
3

---

Training of score functions The above score function (t, H V Vlogp (x) needs to be trained from data. Denoising score matching [33] refers to the following score-matching loss (SML) function to train the score whose parameterized architecture is denoted as G:
min /' Wt Ex Ex.- (-IXo) 6
(|,,(()) - v log (ooo(
(4)
where (1 (rp) is the transition probability of the state 20 at time O towards the state Ty at time t for the forward process (1). The function wi 2 0 is weight function. The default choice in many
literature is that gt v Bo + (31 Bo),, 3o < B1 are parameters, and one chooses the weight function as follows:
Default weight: W/ = 5
3 w.= V e - (
(5)
The quantity wa has the meaning as the standard deviation of X; conditioned on a fixed Xo in the forward process. See $ 3.7 and Appx. I for more weight functions.
Source of errors
There is usually intrinsic error due to an inexact score function. It is not negligible in many scenarios, e.g., there is only finite amount of samples of po available only small neural network architecture is achievable. However, is reasonable to assume that this non-negligible error is reasonably small, and we decompose the trained score function G into
6 (x) - log P (x) + EE, (z),
(6)
where is small, 65 is assumed to be ((1) and the total error is 66, The generative process used in practice has to use law(Y)) = (0, Ia) instead since is intractable in (2). By choosing T also
large enough so that pr (0, La). we can neglect this error due to the finite T. Besides, need some numerical schemes to simulate this generative SDE, which also leads into discretization errors. In summary, there are three sources of errors (1) Pr * (O, la): this is the source of errors in the initial distribution of the generative process; (2) 6 f O: error from imperfect score function from training: (3) numerical discretization of the generative process. The third error can be systemically eliminated by choosing a high-order scheme [21] or an extremely small time step. It has been observed that by choosing more accurate numerical scheme, e.g., exponential integrator, one can reduce the computational costs [35, 36). As for the first error, if one chooses the OU process for (1), pr converges to NO, Ta) exponentially fast and thus = (ln())) where 8 the error between Po and the distribution of generated samples. Therefore, the choice of T is, in practice, not hard to manage. More details about these three error sources can be found in, c.g., [7, 10] or Appx. B.3.
Score Error w.r.t. Time
Sampling Error w.r.t. Diffusion Coefficient h
Given a fixed score error (black). study the effect from blue and green regions separately
by Prop. 3.5 Prop. 3.6
by Prop.
is "additive" from different regions
Pu
Figure 1: Schematic illustration of the main message: the distribution of the score error 6,- w.r.t. time also matters, in addition to the score- motee-mathinng// loss (4). Asymptotically, the score error can be viewed as "additive" and the error from two time regions (blue and green) might decay or magnify as the magnitude of diffusion coefficient ht increases (see the right picture). The yellow region is transitional region whether the effect of h is not easy to decide.
3 Asymptotic analysis of terminal errors
We use the KI. divergence between the data distribution Po and the distribution of generated samples to quantify the performance of generative model, which depends on the error of score function c

---

the diffusion coefficient h and data distribution Po. To extract the main feature, we first let c + 0 and estimate
sample generation error in KL divergence = L(h , , 0()))
Whenever Po and 6t are obvious from the context, we simply write L(h-) L(h", , .PO). Next, we formulate the main problem setup and assumptions in $ 3.1. The expression of T, is shown in Prop. 3.2. Then we let h = h be independent of time, and study how the leading order function depends on h in various settings of the error function _,- Firstly, we consider 6; (x) = F where = (ODE case) and h - (SDE case with large diffusion). When 6,- (x) م (a) is time-localized function at the end of the generative process (i.e., 8 close to T), the ODE case will outperform the SDE case (see Prop. 3.5); otherwise, the SDE case has an exponentially smaller error than the ODE case as + oo (see Prop. 3.4). Secondly, by combing Prop. 3.4 and Prop. 3.5, the tail behavior of + L(h) for general 6 is described in Prop. 3.6. The reasons behind considering this pulse-shape 6,- will be discussed in Appx. E.4.
as a time-localized function and two limiting scenarios: h
3.1 Set-up and assumptions
We fix the time ' and consider the following SDE for the generative process in t E [0,7],
(- (- (Y:) (g - + (h) )2 G (7i) di hy dW, law(Y)) pr, dÝ. n
(7)
2
which can be regarded as perturbed equation (2) of Y, Denote the distribution of Y; as qt Note that Qr depends on both 1ttT and E6, (hidden inside G V po + e however, when by (3), Ør qr po for any h² We quantify the sample generation quality via
(L = log(po/Tr).
Due to the presence of c6 with non-zero in general KL (L > 0. Assumption 3.1. Throughout this section, we assume that:
(1) For the forward process. we assume f((x) = - 9r = 1 without loss of generality.
(2) The data distribution has the density Po-
(3) There exists cy R such that Uo(x) Ief/2 cv. for any x E Rdd where Uo == - -lg po-
Recall that generic 9t is equivalent to the time re-scaling (Appx. B.2). So this choice of 9t = refers to the generic choice in VP-SDE [30]. The second assumption is also mild; in practice, po is a delta distribution, a common practice is that one tries to learn the mollified version (x) := JRd (2xo2ya7z exp (z (z-v)²/20²) poo)))y instead, as in the GAN [27] and early-stop
techniques [7, 10]. The third assumption is not restrict, for example, Uo(2) = |o| and cv As many realistic datasets are almost compactly supported, we expect that Po = é decays faster than a Gaussian, namely, Pu(x) < Cal1=//2 for some Č > 0, which reduces to the third one.
3.2 Asymptotic expansion of the KI, divergence with respect to e
We introduce a time-dependens operator
c (w)(x) (=))))=== (G 1+h,2
hr² Voggi (()=((()))())))))))))))))))))))) + Au(x),
(8)
2
2
which is the generator in the Fokker-Planne equation of gr for (2). Define an operator & as follows: given any function μr, define (h' ((n) to be the solution at time / of the following Fokker- Planck equation with r E [s, 4]: 0,0,- = Lh '(0,), and with initial condition 0, = We define & (h*) (w) = O,, Properties of this operator are collected Appx. C.I.
8,8
We use the notation 7't to refer to ÿ since we need to calculate the its derivative for the small parameter 6. The dependence of ạ on h is suppressed for short notations. We have the following asymptotic result with the proof given in Appx. C.3.
5

---

Proposition 3.2. Define = d as the first-order perturbation of We have
KL(pul ( L(h~) -(+))=++((
(9)
where
1 (x) -(h-)=
2 Ja. Po(
dæ,
Ur = 2 2 (1+
1 (h h
2
1,T
(V (pr (x-)))
dt.
(10)
3.3 The role of h in the Fokker-Planck operator 4/
c(h²)
Let the potential Ug log Pt, and by the notation of time-reversal, U UT When 1+ >
we can rewrite
|a
(w)= Ap V. (VV.
(x) (1 (1+ (x)
(11)
2h,
We now introduce a probability distribution induced by the potential V
P5 x exp(-V;;))
(12)
By convection, we also have V; = V and Pt PT- Note that V depends h and when ht - DC, have v, U, The role of h- in E. (h now be viewed as the time re-scaling and the effect of c(hh at local time t can be viewed as evolving the Fokker-Planck equation associated with the overdamped Langevin dynamics in the time-dependent potential V for (272 amount of time. When ht - oo, _(h²)) can be roughly viewed as constructing an "almost quasi-static" thermodynamies [6] bridging the initial pr and the qquasi-jequilibrimm Po = -- for any distribution 14t (probably far away from pt"), within short time period At slightly larger than (1/hh ?).
we have Mita: ~ Pita provided that s + μe evolves according to C. ;; see Appx. C.2. This key finding will guide our analysis of the solution operator $ (h*) t.T
3.4 Score function is perturbed by a pulse
From (10), we know that vr combines the averaged effect of क (h" a first result, we consider 6,-(2) = E(x)o,_, for fixed time instance [), where 8,-5 Dirac function. In this case, uy no longer involves time integration and have
(pF 8(-) for various
(1+ 2
(h
(pî E)).
uT
2
S.T
(13)
(1 + h 2)2
(1-) IV. (p_ E)))
L(h*)
8
Rd
po
To proceed, we need to make additional assumptions: Assumption 3.3.
(1) For any E [0, 71, Ug = - log is assumed to be strongly convex and the Hessian of U₃ is bounded by two positive numbers m, and M, below
m,la (2U((a) Mcla,
Va Rd.
(14)
Moreover, assume that
> 1, E [0, T).
and
mo > 1.
(15)
(2) For all t e [0, we choose h = h as constant.
Introduce
B = (1 11+h h- h h h- = -1+ h (1+h-2)m7-+ m -2
(16)

---

which characterizes the Hessian lower bound of V (11). Note n. 2 m when h > 1. We would like to explain the reason behind the above assumptions, in particular, their practical relevance. Part (1) Strong convexity: this is a common assumption for Langevin sampling analysis [11, 14]. As the role of c( is essentially simulating a Langevin dynamies with timedeependent potential, is reasonable to use this assumption as starting point. Moreover, the algorithmic improvement in gDDIM [10] is highly inspired by a form with assuming the data distribution as a Gaussian; Fréchet inception distance (FID) [15], widely used metric to evaluate the quality of generative model, essentially treats the data (in the feature space) as Gaussians. Therefore, we believe that this assumption can still capture some main features of realistic datasets. Part (2) my 2 1 for any t E []: The second assumption m, 2 means that pt is more localized (smaller variance) than the standard Gaussian (unit variance), which is compatible with Assumption 3.1 (3). It can also ensure that V2 is strongly convex with positive Hessian lower bounds, i.e., Kj 2 1, for any & E [0, T]
and h (0, oo).
Proposition 3.4. Under Assumptions 3.1 and 3.3, suppose that 6, (x) ((2)) = for some fixed we have the
s E [), h 2 has == max 1/2,, ho(1)), V maxx
cy. lu[2] suPrE[e.T]] C+² (01) }
upper bound of L(h) (13):
L(h) S Ch(1 h2) exp (- (-f((ht (h² CF (x)) O..(.dr).
(17)
where Cn = (v.0r. 105, Cr 20+-30M; m
see Appx. E.1 for details about C and ho(2/)); CU E 1R. comes from Assumption 3.1. Moreover, lim Ch. lom C. exist.
See Appx. E.1 for proofs. We remark that the above bound focuses on capturing the scaling with respect to h2 but may not be tight for other parameters. It remains interesting to see how we can improve the above upper bound. The main conclusion is that: if the error function 6,- is a pulse at time t == s, then for large h, L(h) will decay to zero exponentially fast with respect to h. For ID Gaussian case, we can clearly observe such an exponential decay in Fig. 2a, where we pick Et = TICO.937 V Vlogr, The intuition behind this exponential suppressed prefactor is that for large h, "(A"" a,T can be viewed as an almost quasi-static thermodynamics dragging any positive measure towards pf po, mentioned above in $ 3.3. As V (p= E) has measure zero, we can split it into positive and negative parts: // μ³ V- aassmme f Dو 1 WLOG). Each term (h'
(h (v+) PT = 1 a.T (v-), which explains that ob aTT (vv) = 0 for large h. a.T
L(h)
L(h)
L(h)
10a
10
- 0.4
0. 2 0.6
10-2
0.1
10"
Og 0.6
10-7
0,
0.2
10"2
.... Og 0.8
02 0.3
50
100
50
(a) 6:" Lego.sor log P
(b) 6, -120.9967 V logpt
(c) 6:- = V log P
Figure 2: 1D Gaussian pa = NO ..) (with different on smaller than one) and T = 2. Panel (a)
validates the exponential decay of T.(h) when the score function has no error near L R T, similar to
Prop. 3.4. Panel (b) validates Prop. 3.5 that the ODE model (h 0) outperforms the SDE model when there is a large score error at f. s T. Panel (c) validates Prop. 3.6 that lim T.(h) exists.
3.5 Score function is only perturbed near the end of the generative process
Proposition 3.5. Under Assumptions 3.I and 3.3, suppose that 8, (I) = IIe[T E(x)) where < 1. Then when < 1 and h > 1, asymptotically.
2
(
F))2
L(h) x (1 - e ND
(V F2)2
I(0)-
Rd
Po
2
Rd
Pu
7

---

with = (1 + 1/h3)m 1/h2 as in (16). The upper bound of L(h) is tight asymptotically.
We remark that we made no assumption elh=e-ah"ro/ r0/2) /a² Case (I): If ah² > 1, then L(N//100 S raa which is large as a < 1. Case (II): If ah2 < 1, then L(h)/2(0) & h², which is still large. In either case, 4(h)/((0) is expected to be large for general E and the ODE model is preferred in this case. The intuition is that there
how ah? scales. I > 1, 10h/110) S
almost no time for the operator C.Ph to suppress the error E, so the prefactor + (hg )² in ur (10) dominates (which is the key difference compared with Prop. 3.4). The proof is postponed to Appx. E.3. In Fig. 2b, we consider ID Gaussian and only perturb the score function at the end of the generative process (E, 420.9957 V log p ); clearly, the SDE-based models have comparatively larger error.
3.6 General error in score function
We can generalize Prop. 3.4 and Prop. 3.5 a general error function 6t, and observe that L(h) will converge to a constant exponentially fast as h - 00.
Proposition 3.6. Under Assumptions 3.1 and 3.3, we consider a general error function (t, x) + Et (x). Let 7 = infnztu infre/0,77 KÇ. For any E (0,1) 1) and B E (0, 2), when h >> 1,
L(h)= x (1 + T + (1 (1++*
C Cr (1+ h²) exp 9)) -
L(h) 2 ( -1
exp
h2 P7))
where C given in Appx. F. and TT (only depending on Po and 6 ) is upper bounded by
(V
()
2 1 / (Vlogro of + V VV.T Po- (18)
OSTS
2mg JR"
Po
2mg
In the limit h + 00. 00. (1 (1+ (7 where a E (0, 1) is arbitrary. Hence, the tail behavior is that (hh) converges to T exponentially fast h xc. If we assume that 65 V log Po. Po = (O, o aa) in d-dimension, then the above upper bound is simply T S d, which is independent of σ (see Appx. F.2). For ID Gaussian Fig. 2c, we can indeed observe that lim (h) exists, and is bounded above by d = 1; see Appx. G for more types of eгror functions.
The above upper bound has an interesting tight connection to the generator of (overdamped) Langevin dynamics with drift- VU, = Vpo. If we adopt constrained score models [23, 25], namely, parame- terizing log Pt instead of the score function V log Pt during training, the error 67 Vip for some scalar-valued function yp Then the above upper bound becomes
1
/
/ (C'
(Ap VUo Velee- Uo
2mg Je
(19)
2m
where (*(x) Ap- V Vip whose adjoint operator C(u) V (VUnH) + Au the Fokker- Planck generator of the following Langevin dynamics dxX, = -U dt + vV2 dW, where Wa is the standard Brownian motion. We remark that this formula (19) is general for constrained score models [23, 25]; see also F.2 for elaborations. An interesting open question is whether and how we can take the above upper bound into consideration when designing the loss function.
3.7 An application: exploring the effect of training weight
The above theoretical discussions suggest that diffusion models with large diffusion coefficients are more negatively affected by score error near data's side, whereas the ODE model is more negatively affected by the score eтror near the noise end. This leads us to conjecture that if we can control the training (e.g., by optimizing the training weight wx), so that the score error distribution near the noise end is reduced and meanwhile the score-matching loss is not significantly impacted, then it will very likely improve the ODE models. We report preliminary numerical experiments to support this idea in Appx. I, whereas comprehensive study will be left future works.
8

---

4 Numerical experiments
We present experiments on 2D Gaussian mixture model, Swiss roll, CIFARIO to support our theoreti- cal results: when numerical discretization error is not dominating, the sampling quality increases as h increases, a reminiscent of Prop. 3.4. Experimental details are postponed to Appx. H, as well as more numerical results (e.g., results about ID Gaussian mixture and MNIST), Results by adopting various weight functions, a technique arising from theoretical predictions, are postponed to Appx. Source codes available at https: thub. com wwccayyco/ppmmmaaifffi...cmmmeee
2D 4-mode Gaussian mixture. We verify the theoretical results on 2D Gaussian mixture with a specified score error. In Fig. 3, clear trend is that higher h produces generated distributions that better match the marginal densities of po. and it is numerically verified by the purple line of Fig. 4a. In Fig. 4, with multiple settings of g and e, we observe consistent phenomenon that as increases, the KL divergence of true data and generated data converges exponentially fast, thus validating Prop. 3.6. It worths noticing that in all three settings of & by simply adopting larger diffusion coefficient h in (3), we can obtain better generative models without any extra training costs.
Example 1:
h2 0.0 g2
h² 1.0 g2
h² =.0
h²= 3.0
h² 4.0
0
6
2
0
Figure 3: Visualization of marginal densities of 2D 4-mode Gaussian mixture, where & = V log pf and 0.2. The top shows the marginal distribution of the first coordinate and the bottom row for the second coordinate. The red lines are the exact marginal distributions of po and the histograms (blue) visualize the empirical densities of generated samples.
KL divergence
divergence
KL divergence
02
0.1
0.
0.0
10-2
50
10 2
0/2
10-2
10-3
10 3
2 h?1g?
3
2
1
4
h//02
h²/g²
(a) E
V log
(b) 6
14sin(2n1/T)
v logpi
(c) 61-
= Teco.asr V log
2
Figure 4: Numerical results of 2D 4-mode Gaussian mixture. The above panels show that the KL divergences between the true distribution and the generated samples overall decay as h increases for various types of error perturbation of score functions.
Example 2: Swiss roll. We consider Swiss roll, more complex distribution without exact score function available. We train the score function with the denoising score-matching objective [33] (Appx. H). The first plot in Fig. 5a shows the difference between true data and generated data measured by Wasserstein distance, which decays to zero exponentially fast, verifying Prop. 3.6. In the second plot of Fig. 5a, we tested various and time steps for the generative process. The ODE model (h 0) does not improve, as the number of time discretization steps increases, but near the continuous-time limit, all SDE cases (h > 0) are better than the ODE model. suggests that our conclusions here is limited to the continuous-tine setting. The exploration of time discretization errors will be future works. The generated data results in Fig. 5b demonstrate that with increasing h. the sample quality is improved; see Appx. H.6 for more results.

---

Wassersler disance
M-00
h''1og 2--20
2.0 x10
4
10
10 x
17
10
10
Discretbation stops
(a) Improved model performance with larger
(b) Visualization results selective
Figure 5: Numerical results of Swiss roll. Panel (a) shows the decay of Wasserstein distance between the true distribution and the generated samples with increasing h and 20,000 uime-dscretization steps, and the decay of Wasserstein distance with the increasing number of dinediscretizzaton steps and different h. Panel (b) shows generated samples with different h.
Example CIFAR-10. When using a large amount of parameters for score matching in practice, we observe that SDEs appear to perform better than ODEs as discretization error descreases, result similarly observed on Swiss roll. This implicates the practical applications on generating samples of better quality under given (pre-trained) score-mathhig model.
Table CIFAR-10:
We evaluate FIDs with different discretization steps and h2 h/88 on a pre-trained checkpoint entitled "vp/cifar10 ddpmpp_ ccrtiiuuuunuutuu in [30]. We do not use any correctors and
evaluate FIDs on 10a samples, thus the results for h² //2²" = 0, are different from the reported values.
Discretization steps
100
200
500
1000
2000
3000
4000
h²/g" = 0
22.43
8.12
7.20
6.89
6.98
7.27
7.33
h2/g"
31.72
15.72
7.23
6.70
6.71
6.95
7.08
h2 h2/g² =2
52.77
26.68
10.78
6.99
6.70
6.69
6.98
h²/g" =
92.83
46.11
20.47
10.17
7.20
7.09
7.01
Summary and outlook
In this work, we study the effect of the diffusion coefficient the quality of overall sample generation the generative process. Theoretically, we provide understandings of scenarios in which the ODE- based model and the SDE-based model is superior than the other; see Prop. 3.4 and Prop. 3.5. Numerically, these results are validated via toy examples as well as benchmark datasets.
There are many interesting directions for continuing works. (1) As we focus on the asymptotic case, time-dependent hi- with large magnitude (i.e., h >> for all () is essentially no different from constant h = h with h > 1. Whether it is possible to use time- denedeeenddntt + h to combine the advantages ODE and large diffusion cases dealing with different types of error of score functions in the non-asymptotic region is still an open question. (2) Can we design a practical criterion that directly learn the optimal magnitude the noise-level function h;- by looking at the score-training error distribution? Can we develop certain theoretical understanding of the empirical results in [19]? (3) How can we find a stable and accurate numerical scheme to deal with the fast diffusion case? (4) How can we generalize the above theoretical results by removing the convexity assumption, and including the low-dimensional feature datasets into the theory [8. 5, 18]?
Concerning the broad impact, though we don't foresee any negative social impact of this work, the potential improvement of generative model might relate creation of "deep fakes".
Acknowledgments and Disclosure of Funding
YC is sponsored by Shanghai Pujiang Program, and acknowledges the support from City University of Hong Kong during his visit. JC and YL are supported by the NSFC Major Research Plan Interpretable and General-purpose Next-generation Artificial Intelligence (92370205). XZ is supported by General Research Funds from the Research Grants Council of the Hong Kong Special Administrative Region, China (Project No. 11308121 12.1 11318522) and the NSFC/RGC Joint Research Scheme |RGC Project No. N-022YU102/20 and NSFC Project No. 12061 160462].
10

---

References
[1] Michael S. Albergo, Nicholas M. Boffi, and Eric Vanden-Fijnden. Stochastic Interpolants: A Unifying Framework for Flows and Diffusions, 2023. arXiv:2303.08797 [cond-mat].
[2] D. Bakry and M. Emery. Diffusions hypercontractives. In Lecture Notes in Mathematics, pages 177-206. Springer Berlin Heidelberg, 1985. doi: 10.1007/btb00055447
[3] Dominique Bakry, Ivan Gentil, and Michel Ledoux. Analysis and Geometry of Markov Diffission Operators. Springer International Publishing, 2014. doi: 10.1007/978-33319-00077
[4] Fan Bao, Chongxuan Li, Jun Zhu, and Bo Zhang. Analytic-DPM, an analytic estimate of the optimal reverse variance in diffusion probabilistic models. In International Conference Learning Representations, 2022.
[5] Valentin De Bortoli. Convergence of denoising diffusion models under the manifold hypothesis. Transactions on Machine Learning Research, 2022.
[6] Herbert B Callen. Thermodynamics and an Introduction to Thermostatistics. John Wiley & Sons, Nashville, TN, 2 edition, 1985.
[7] Hongrui Chen, Holden Lee, and Jianfeng Lu. Improved analysis of score-based generative modeling: User-friendly bounds under minimal smoothness assumptions. In Proceedings of the 40th International Conference on Machine Learning, volume 202, pages 4735-47665 2023.
[8] Minshuo Chen, Kaixuan Huang, Tuo Zhao, and Mengdi Wang. Score approximation, estimation and distribution recovery of diffusion models on low-dimensional data. In Proceedings of the 40th International Conference Machine Learning, volume 202, pages 6672-4712 2023.
[9] Ricky T. Q. Chen, Yulia Rubanova, Jesse Bettencourt, and David K Duvenaud. Neural Ordinary Differential Equations. In Advances in Neural Information Processing Systems, volume 31.
Curran Associates, Inc., 2018.
[10] Sitan Chen, Sinho Chewi, Jerry Li, Yuanzhi Li, Adil Salim, and Anru Zhang, Sampling is as easy as learning the score: theory for diffusion models with minimal data assumptions. In The Eleventh International Conference on Learning Representations, 2023.
[11] Xiang Cheng and Peter Bartlett. Convergence of Langevin MCMC in KL-divergeece. In Algorithmic Learning Theory, pages 186-211, 2018.
[12] Jooyoung Choi, Jungbeom Lee, Chachun Shin, Sungwon Kim. Hyunwoo Kim, and Sungroh Yoon. Perception prioritized training of diffusion models. In 2022 IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR), pages 1 11462-11471, 2022. doi: 10.1109/ CVPR52688.2022.00
[13] Florinel-Alin Croitoru, Vlad Hondru, Radu Tudor lonescu, and Mubarak Shah. Diffusion models in vision: survey. IEFE Transactions on Pattern Analysis and Machine Intelligence, pages 1-20, 2023. doi: 10. 1t..199//AAMM..023322666 1988.
[14] Arnak S. Dalalyan and Lionel Kiou-Durand, On sampling from log-concave density using kinetic Langevin diffusions. Bernoulli, 26(3):1956-1985 2020. doi: (0.3110/11-BBB 178.
[15] Martin Heusel, Hubert Ramsauer, Thomas Unterthiner, Bernhard Nessler, and Sepp Hochreiter. GANs Trained by a Two Time-Scale Update Rule Converge to a Local Nash Equilibrium. In Advances in Neural Information Processing Systems, volume 30. Curran Associates, Inc., 2017.
[16] Jonathan Ho, Ajay Jain, and Pieter Abbeel. Denoising diffusion probabilistic models. In Advances in Neural Information Processing Systems, volume 33, pages 6840-6851 Curran Associates, Inc., 2020.
[17] Chin-Wei Huang, Jae Hyun Lim, and Aaron Courville. A variational perspective on diffusion- based generative models and score matching. In Advances Neural Information Processing Systems, 2021.
11

---

[18] Ian T. Jolliffe and Jorge Cadima. Principal component analysis: review and recent devel- opments. Philosophical Transactions of the Royal Society A: Mathematical,, Physical and Engineering Sciences, 374(20655:211502222 2016. doi: 10. 0099/rsta.2015..00
[19] Tero Karras, Miia Aittala, Timo Aila, and Samuli Laine. Elucidating the design space of diffusion-basee generative models. In Advances Neural Information Processing Systems, 2022.
[20] Fima C Klebaner. Introduction to Stochastic Calculus with Applications. Imperial College Press, 1998. doi: 10.1 10. 10.1142/p110.
[21] Peter E. Kloeden and Eekhard Platen. Numerical Solution of Stochastic Differential Equations. Springer Berlin Heidelberg, 1992. doi: 10.1007/978-33662-1283 1-616-5.
[22]
Alex Krizhevsky. Learning multiple layers of features from tiny images. 2009.
[23]
Chieh-Hsin Lai, Yuhta Takida, Naoki Murata, Toshimitsu Uesaka, Yuki Mitsufuji, and Stefano Ermon. FP-Diffusion: Improving score-based Diffusion Models by Enforcing the Underlying Score Fokker-Planck Equation. In Proceedings of the 40 th International Conference on Machine Learning, 2023.
[24] Otto and Villani. Generalization of an Inequality by Talagrand and Links with the Logarithmic Sobolev Inequality. Journal of Functional Analysis, 173(2)3361 - 400, 2000. doi: http:///di..org1000 1ttg//001006/ffan 1006/jfan. 1999.3557.
[25] Tim Salimans and Jonathan Ho. Should EBMs model the energy or the score? In Energy Based Models Workshop-ICI.R 2021, 2021.
[26] Jascha Gohl-Dickstein, Eric Weiss, Niru Maheswaranathan, and Surya Ganguli. Deep unsuper- vised learning using nonequilibrium thermodynamies. In Proceedings of the 32nd International Conference Machine Learning, volume 37 of Proceedings of Machine Leurning Research, pages 2265, Lille, France, 2015. PMLR.
[27] Casper Kaae Sønderby, Jose Caballero, Lucas Theis, Wenzhe Shi, and Ferenc Huszár. Amor- tised MAP inference for image super-resolution. In International Conference on Learning Representations, 2017.
[28] Jiaming Song, Chenlin Meng, and Stefano Ermon. Denoising diffusion implicit models. In International Conference Learning Representations, 2021.
[29] Yang Song and Stefano Ermon. Generative Modeling by Estimating Gradients of the Data Distribution. In Advances in Neural Information Processing Systems, volume 32. Curran Associates, Inc., 2019.
[30] Yang Song, Jascha Sohl-Dickstein, Diederik P Kingma, Abhishek Kumar, Stefano Ermon, and Ben Poole. Score-based generative modeling through stochastic differential equations. In International Conference Learning Representations, 2021.
[31] Yang Song, Prafulla Dhariwal, Mark Chen, and Ilya Sutskever. Consistency models. In Proceedings of the 40th International Conference on Machine Learning, volume 202, pages 22211-32252. PMLR. 2023.
[32] Esteban G. Tabak and Eric Vanden-Eijnden. Density estimation by dual ascent of the log- likelihood. Communications in Matheratical Sciences, 8(1):217-233, 2010.
[33] Pascal Vincent. A connection between score matching and denoising autoencoders. Neural Comput., 13(7):1661-1674 2011. doi: 10.1162/NECO 00142.
[34] LeCun Yann. MNIST, 1998. URL //ytt: www.euwwoo//eeeeeeeeeeee//////ee hemm///dbmmmis/..//////mmmm
[35] Qinsheng Zhang and Yongxin Chen. Fast sampling of diffuuson models with exponential integrator. In The Eleventh International Conference on Learning Representations, 2023.
[36] Qinsheng Zhang, Molei Tao, and Yongxin Chen. gDDIM: Generalized denoising diffusion implicit models. In The Eleventh International Conference on Learning Representations, 2023.
12

---

Supplementary Material for "Exploring the Optimal Choice for Generative Processes in Diffusion Models: Ordinary vs Stochastic Differential Equations"
A
Notation Convention.
Table 2: Summary of important quantities in this paper
Quantity
Notation
Notes
Forward process
X
0: data distribution
Backward/generative process
t 0: noise distribution
Backward process with inexact
Y2
Distribution of forward process
== law (X.)
Distribution of backward process
Qt
law(Y,) = St pr- pt
Distribution of approximated backward process
qt
Qt = (aw(Y)), qt = 4: when
eπrо c 0
Error function of the score
E,,
0 S 6 < and & = (1))
The exact potential
U2
Pt -
Modified potential
Vi=
see (11)
Normalizing constants
Zy = ʃ -v
is arbitrary
B Discussion and proof for $ 2
B.1 Proof of (3)
We re-state the conclusion in (3) in the following lemma:
Lemma B.1. For any function h if one chooses A as in (3), then we have qe(x)) = (r))(x) for any t E [0 and x E Rd.
Proof. We can easily write down the Fokker-Plande equation of (1)
(.)(()) V (s((z),(z)) of. Ap((x).
Since need to ensure 9: Pr require
(x) V. (Ji (x)) -
(g((a)
V (1r (z)
-;- h
v(((z) + Aq((z)
hr2
2
2
V (UFF (z)
gi
Vlog (logq(())) (loga(z)))((1)) (xx)
ht²
Aq((x)
2
V (UT
Di
+2
+
+h,
Vlogri (x) (z))
h
Ag((x)
2
2
h
(47(4))
Aqi(z),
2
by noting that A (x) specified in (3). This equation is exactly the Pokker-Plande equation of (2).
B.2
The role of gt
The function ++ 9r as the diffusion coefficient in the forward Fokker-Planck equation (1), in fact, plays role as time re-scaling as long 9i 0 for any t. Indeed, if we have SDE in the following
13

---

form (the VP-SDE in [30])
dX, = 96 X, dt + 9r dW,, 2
E [0, T).
then by introducing т RR + RR via the ODE r'()) (9-(1) -2 =(0) = 0, we know that X X-(t) satisfies the following SDE [20]:
dx, = 2 X, dt + dW,, E [0,
where 0 :== (T) means the inverse function of r at time T. The Brownian motion W in X, is not the same Brownian motion X,, i.e. the driven-noise in the last two equations are not the same technically; we slightly abuse the notation for the simplicity of notations. By Lem. B.1, this SDE has the backward process as follows:
dY. Y
1 h
V log (Y.)) de + h dW, te [0,0]
(20)
1 h log Po (Y,) dt h dW,,
2 Y,
where P, is the density function of X, by notation conventions in 2. By Lem. B.1, far know that
law (ow(Y--)) Lom (X,) lla((X (1)))
Let fo = - (T- s) and Y. == ,, where E [0, T). By straightforward computation, know
dY, = f/ (/12x,, Mr +
1 h+ fo Loggo- ( (I))) at + v 1/ h2- aW,
2
S/ -
1 1 h2 -1 Vlogpr- di S3 hう dW..
To get the second line, we used the fact that
Da-f_ law law(Xx- (.) law (X,- law (Xw(X = pr-s-
By chain rules, it is easy to compute that
1
1' (TT (9T-s)
-6
Therefore, the above SDE for Y, has the form
2 (T (9i- Y. -7-
+
dY, This matches the form in Lem. B.1 by choosing
v log p5 (Ys) dt - dW..
h = 9r -r = 9r hT-1(T-s)-
Hence, if we simply pick = c as a constant function (20), it has the same effect as choosing h cg." where E R+ The former (h is used in 3 for simplicity in theoretical analysis, and the latter (hg- cg.-) is used in numerical experiments in $ 4 to match previous literature in practice (namely, general g). In conclusion, the above discussion justifies the consistency of notation and set-up between our theoretical analysis and numerical experiments.
B.3 Existing analysis sample generation quality
In practice, one simulates the following SDE:
dz, = (- ft (Z,) + (9) 2 + (hr 6, (Zi) dt ht dWe
law(Zo) = N(O0 (21)
2
14

---

From the analysis by Chen et al. [7, Theorem 2.1] for the case htT =g- one has
(pulll ( = (L(prlIW(O, + o(T)) + O(T)
where T is the distribution of Zr after applying the exponential integrator to (21) and N is the number of time-discretization steps. The first term (L(prllW(O, - (M2 + d)e - where
M2 Epo |x| is the second moment of data distribution.
Therefore, to ensure that (poll)) S 6, it is sufficient to choose
T O (log ((Ma+d)/5).
The fine-discretization error be eliminated by choosing N cc. What so far less clear in literature is the term (
B.4 Existing error bounds appear to fail to characterize the optimal /t
Note that the target dynamics qe in (2) and the approximated dynamics qe in (7) only differ in the drift term. Recall that how the quantity KL.(gel (+)) changes for / E (0, T). By [7, Lemma C.I]. for any E R,
and to estimate KL ( = (L( (gg]lTTT) we can simply estimate
h I (gel +
2 + h
(geL(
E6i (Xi) V log
(22)
2
where the Fisher information
I (pllg) = ff dp
v log //
By Sauchy-Schward inequality,
0.KL(gellee)
hr²
I (qellat) + Dt
- h B (x)))
1
log de (x()) (Y,
12] |
2
20
(23)
(h+ + 9)
2)2
8h1 E E161, (Y)))?².
where we chose 2² h 240
the last line. This bound captures the scaling extremely well when = which helps to establish the effectiveness of score-based diffusion method in [7]. However, this bound not directly applicable for general h², as it is clear that this bound fails to provide useful information when ht z O (namely, the probability flow), as well as the large
2h;
diffusion case (h* + co). From directly optimizing the upper bound, i.e., minimizing (h--++, ends up with the choice ht g², which has been used in many literature. We acknowledge that ht = g is an effective choice; however, as we show in $ 4, this is not really the optimal case in general, and the above argument cannot justify choosing the diffusion model over the probability flow.
sh:-
C Discussion for $ 3 and proof of Prop. 3.2
The operator @ (h*) 3,8
C.1
Lemma C.1.
( satisfies the following property, i.., any s, r e IR
(h
क (h* s,t
(h²)
arr
(24)
Moreover, (h Id is an identity operator for any l. For any 8, RR.
0, (4
(h") (4)) - (c. (b))
(1*))
4)
a. th" (W))
-1th
100 (4)). (25)
(h-)
15

---

Proof. The structure in (24) is easy to imagine and is thus omitted herein. Next, we shall only prove the first result in (25) for illustration:
B (h' )
s468, s468, (1)
th'
(())
lim 0, ( n,t 8a-401
(h (a))
ôs
th x )
sit otos, (11)
0.t (w))
lim
68-101
ōs
(
5,t s455,5 ldd (1)
( (M
lim
68-101
8s
5,t (Id 8oL. (h-) (d)
(
lim
68-201
ōs
Φ (h- (E. (mi).
Lemma C.2. Suppose (t, + pc(x). (t, (x(x)) are time-dependent functions. Suppose
0.p4 c (mx) + 8₂ with pa = 0, then
PT 1"d
(1)))))) 1T (0c) dl.
0
Proof. Note that
(h*) 29 (1**) ) 0.(4 m1) 1.0 (-0001(2 (-c
(1²))
(6s) + (0000))))) (c(h))) Me Ne+00 0,) 10()(()))) (0.).
Therefore,
(h' ) (UT) = So" D (0.)) dl
$
T,O
& (h)) 0,7 to both sides, we have
By applying the operator
1.0 ((0,) dt (24 (" 2 (0,) dt.
th'
PT
Por
C.2 Time re-scaling of L(h and connections to almost quasi-static Langevin process
By time re-scaling, it immediate to obtain the following result:
Lemma C.3. Suppose hū- foo any t E [0, for simplicity. For any probability distribution the probability distribution & (hh (with < t) the final state of the following PDE on the time interval F E [0, (t-o)b²/2))
O = ATF V (VV
Too =
(26)
It corresponds to the following Langevin with time deneddepnndent potential:
dX; -VV- Taapa (X;) dF + /2 dWr,
law(Xo) p..
Proof. To obtain $ (h (n), we simply solve the following PDE:
h2 2
0,
Aur V (VV,
H1. μ..
By change of variables F F= (r >)b²/// /2 and jop == Pr M4425/22 we immediately have (26), and the
corresponding Langevin dynamics easily follows.
When + oo, the potential in the Langevin dynamics F + V
aayaa evolves extremely slowly. From (11), we also know that V */2r//a s U when F 0 and is approximately U.- when (t-ajh²
Therefore, the Langevin dynamics can be viewed as an almost quasi-static process [6] approximately transforming the state p to the state p over extremely long time period though.
16

---

C.3 Proof of Prop. 3.2
Denote
24 (x(x) == d=
G( == 0 (=)
Namely, we expand of via the following
(x) 2 (xv)(() G(x)) +
Eqe
Recall that when = 0. & gci in particular, 97 = qr Po- The cost function can be easily expanded via Taylor's formula:
Lem.
2²
+ [x) po(
(LolllT = -E Jea (T( dx
(GT(x)
da (c²)).
2 Jgd
Since f = 1, it is easy to know that J fOr = Cr = 0 and thus
e² 1
(pollo)) 2
dr + O
Rd (o(
Next we need to study er. Recall from (7) that the Fokker- Planck equation of ge is
1 1 h
O.q (h
2
By taking derivatives with respect to c on both sides of this equation and then passing c - 0,
0,00 c. (0r) V
1 + h
6, 9t).
By Lem. C.2,
1 T VT
2 / (1 hi 2) (1*)
/
1
1.T
V
qt
61-)
dt.
(27)
EPT
D Preliminary results
We collect two lemmas which will be useful in proving Prop. 3.5 and Prop. 3.6 later. Lemma D.1.
1. The operator K((f) := Af + V. (VV. f) is a Hermittinnselff-ddioiin operator in the space
(:. (1-0)1)66-1) where & e V. Therefore, K. has the eigen-decompositiinn in this weighted
(p= space.
2. Assume that (2V(x) a any Rdd for some Then the operator -. is positive operator on the space {f f f 0} with spectrum gap at least m.
Proof. For any f. note that
(f. (f.K)) = JRc v(f/p) Vgg/p)p
Therefore, IC is Hermitian operator. the space {f f 0}, we know
-(IKKf11a-)) = [ |IV//D)0
Poincan incg. m
/ -
The validity of Poincaré inequality under the strong convexity assumption well-known in literature; see, c.g., [2, 24] and [3, Corollary 4.8.2].
17

---

Lemma D.2. Suppose the operator K. is defined as K(f) := + V (VV . (f)) with v2V(2)2
ml for any E R. Then for any C, RR, and any function yp with /p 0,
J - ( exp (c)
achp
ak
cd
k=1
aco
4²
cm
P
where p x e {(Ax ox) Hm are eigen pairs for the operator Remark. Note that as K. is an operator, exp (ctK) : =
Kc and
= Ek= aLok.
[etK]* is the operator exponential.
M!
Proof. Let us first denote the eigenvalue decompasition of K. as K(Øk) = d where A 2 m for k E N and (Ø); @4) (x)))))))) = 65,6 for any j. E N by Lem. D.1. Then we can decompose yp by
Σx agoi is not hard verify that
06
exp (ct
Σ
ad agok.
Then
L. exp (o) dt
00
1 -
accd
ag agok.
cAk
Hence,
1 1 ! p (L cxp (ct de) at)²
0
eaack
=
Rd
P
CAA
20
-ockk
OKk
Cλk
E
Σ
em
aacm
acm
cm
E Proof of Prop. 3.4, Prop. 3.5, and a discussion on the pulse-shape error
E.1 Proof of Prop. 3.4
For convenience, we summarize some notations below; see also Appx. A.
Denote the global minimum of U. as X and denote the global minimum of V as V When -> 20, we know
lim Vv + Xr h
Denote the normalizing constant Zy := == -V for an arbitrary potential V. Recall that the distribution of the exact dynamics Y, is de = pT-1 = pt and the distribution of the approximated dynamics Y, is qr-
Recall that p: := -U: and Pt x e in (12).
Lemma F.1. Under Assumption 3.3, we have:
(i) Assumption 3.1 (3) is valid, i.e., the existence of cu.
18

---

(ii) For any t E [0, T the probability distribution Ps defined in (12) satisfies the Poincaré inequality with constant Kt (16):
J Vpl dp: 2 Ktt f x² dpes
Vp with [ up
(28)
(iii) Zvo :- -g e Vo is both upper and lower bounded: for any d E (0. 1). there exists ho(d) > 0 such that whenever 2 ho(o),
1- S Zvo Se Se- -0
and
P0 e-h
_ -
Po
(iv) If we pick 8 = for (hx V ma bounded:
or in(2]
}. the function po/po is uniformly
Pu
S
(29)
Po
Proposition E.2. For any 2 1/2, we have
Upee r² d Pi
/
6² pt 0t
/
VpppE
pF.
Vip CR
Rd
Rd
where
20d 5t
c- = şt + +2 - Xx,1
85
20 + 30M,2
CT c-
mi
mi
(30)
5d(1 + Mt)
5
5t 1r log Zv,-
11,, 1,,1-11
S 5( 3M
Zv.- Se-t, Mr and m. are shown in Assumption 3.3, and X,- and yi are the global minimum points of functions x + U₁ (x) and x + V (x), respectively.
Remark. A similar result holds for any h > 0; choose h 1/2 in order to simplify the above constants (30).
Note that the constant (2) is independent of h. We remark that when h >>>> 1, CF mately behaves follows: Lemma E.3. When h - we have
(1) approxi-
5d(1 + M,)
+ 1x (50 + 75M, 7
lim C C1
h-000
2
m
We now proceed to finish the proof of Prop. 3.4. The detailed proofs of Lem. E. Prop. E.2, and Lem. E.3 are postponed to Appx. E.2.
Proof of Prop. 3.4. By Lem. E.
( (h)) T.T
V (p_ E"))
(1 + h222
I(h)
(3)
Po
Rd
Po
Pu
(31)
$ (h aT
V (p_ E¹))
(29 1
(1 + h222
2
Rd
Po
19

---

To simplify notations, let
A=2 ga Pt
/
A = o 5.
(ps E)),
(h,
I LE
18,
After taking the time-derivative and using Lemma C.I (in the first line below), and expressions (11) in the second line below, and integration by parts in the third line, the expression of pi (12) in the fourth line, we have
/
At 4
(AF)
AF² (p=)² 2² Pt
J(h t) =2
dt
JBd
Pf
Rd
@ h2
Ar
(DAF + V (VV;-A1)) Sfe
Ar2 d
(pt) dt Pt
pt
(32)
AF
+ +V, A5) - / (pt de
A,-
=- h2
V
(VAF
Rd
Pi
Ar
2
AF2
h2
V
Pt
Pt
Rd (p=)) dt
P:
The major challenge is to estimate
fiZd
A72 (ot)²
dePi
d
By Prop. E.2,
di [h, -- - -h2
V Pt
A5 2
dop + tF (1)
/
2
A
dpt
Rd
Rd
Pi
/
A
-C )tt + c.,
)²
dpt.
Pi
To verify the condition in Prop. E.2, we can readily confirm that
A
(h 2
(v. = / v (v, E) - 0,
where the third equality in the last equation comes the fact that the Lokker-Planck solution operator (h')) preserves the total mass. Note that we need h c² 2 0 in order to ensure the correct direction when applying the Poincaré inequality above, which explains the lower bound that needs
s.t
to satisfy in Prop. 3.4.
By Gröwnwall's inequality, for any 2 8,
J(h, < 8),, s) exp (- ((h2 -2--())) (0)rr)
(33)
Finally, by (31), we have
(1 + h2)²
J(h, T) C exp (-/, T ((h² - 12012) - c..())rr)
L(h) <
2
where
Cn - J(h, 8) /
(V (p E))
(34)
p
Recall that Cx=(2)) (30) is independent of h. By Lem. E.3, we already discussed that limp- C (1)(1 exists. When h -> oo, we know that p3_ - p3_ (11)(12), and thus,
lim C₂ / 00000
(V
(pt E))
gg
p
This completes the proof of Prop. 3.4.
20

---

E.2 Proof of Lem. E.1, Prop. E.2, and Lem. E.3
Proof of Lem. E.I.
(i) From mo > 1 in Assumption 3.3, x + Uo(x) - x²/2 is strongly convex and thus is surely bounded from below by some cy IR.
(ii) The second result is classical result by Bakry-Emery criterion 3], V (12) is strongly convex with Hessian lower bound Ki (16).
(ii) To prove the third one, note that
Zvo
Vo (1))
-1+h 2)00(())h dr
Rd
Rd
e-Uo(2) -h (Uo(o)))1 dæ
Rdd
Se 200
e Uoo - 2eu.
To prove the lower bound.
Zvo 2 e -1 dx e { e ( dr. (35)
[ nd
2d
Recall from the beginning of this Appendix, Xn defined as the global minimum of x + (z)).
Then, Uo(x)) 2 Uo) for any and we know h + -1+h (Vo(s)))(((x)) is monotone increasing. By monotone convergence theorem,
lim f e (UC(())) de / lim e (1+h ()( dr h-000
j h
U6())) dr vo(())
Rd
Thus, the limit the right-hand side of (35) is when h + co. Hence, the lower bound of Zv, follows immediately. Finally, when h 2 ho),
Po (x) c (1+h P) 9|x12/2 Po
(36)
Zvyevu Uo
8
(iv) When we pick 8 = 1/2 and choose as specified, we can immediately obtain (29).
Before we prove Prop. E.2, we need the following three lemmas. Lemma E.4. Under Assumption 3.3, have
(.U[(()) < d(1 M,)
x 2
3M2
|r -x,/2
Proof. Since Pt follows the Fokker-Planck equation 0,p, V. satisfies 0,U,(z)
( xpt) + 1 Ap, then U, = log Pe
(d (U.(=) (1) + (.)
)) Then by Assumption 3.3,
d
1
r vU.(2)) AUl(s) + 1
(0.)(())
(triangle inequality)
2
2
d
1
(l=2= + v). AU)(x) + 1 (V)22
2 +
(Cauchy inequality)
d
1
[Vi((=) - vUi(2))2 + d.v.
(by VU,(E) 0)
2 1
4
d(1 + M.)
|x|2
3M? Ir - X
2
Recall that X, is defined as the global minimum of U, and thus VU.(K) 0.
21

---

Lemma E.5. When h 2 1 2'
d
de Teli
(2) s Et + 1 |o
(37)
Pj
where
5d(1 + _M
5
ςt = | log Zv. +
+
1
2
2
1
3.
c = 5( 2 +
4
Proof. By direct calculation from the definitions of Pt in (12) and V in (11),
dep = -8, log Zv-))r - (1 + h hh-))ppt ) 0,U,
Hence, by Lem. E.4 and h 2 1/
Oept
d(1 + MM,)
x
2
3M,²
(x) S log Zv;
+ (1+h- (
2
X, 12)
Pt
2
4
4
<(10, log log Zx:-
5d(1 + +M,)
+
2
12-1
2
18 2,
X
2
3.M, -
+ 5(
+
x - X,
2
2
4
1
3M,
= + 562
1x:
X
Lemma E.6. Assume that P x e -V and & decays fast enough such that
Jaam V (v') +0
as the radius R - 00 (Bo(R) is the ball centered at with radius RJ, e.g., E CR(R)). Then
f 1V 4 /
/ x² AVP
(38)
Rd
Rd
Proof. The identity V (p*Vp) Vipa. Up+ 6² Ap and the trivial facts Vp -PVV, Ap ([VV - AV)p show that
/ V 0² Vp p2 Ap
0
2 /
(VV. Vp)pp + / |V +VV)3
Canchy incq.
1
/ |VV1202> / / / @V / /VViT'
= / /VV1" - / Vip'0 &2AVO
Therefore,
/ / Φ² AVP.
22

---

Proof of Prop. E.2.
The above three lemmas show that
1
d
(² + 5.- x X, Pp=
Pr
Rd
dt
nd
/
St
6² Pt
+ +t
0/11x - -V- +3² - 12.
Rd
Cautyy incy, -
(5t + +2x D/= 1²)
4² PT + +24 1
1 PP PE
Rd
(4)
(5t + 2, 5 X 19)
4²p1 24
VVIVV, (0) - -V
m.
24
= (4) + 2 1.tt xx P)
4 P
m
3 (7) + 1 x P)
/
4² PT
24
(4
/
PuPPE +2 jx
/ 4² AV, oF)
m
jan
04
2 (ET 255 1 xT-1)
24
(4
P + +10M5 d JR
(2p(t) /12
m
20M dst
85 m
= (st + 25 DJ=t X 1
Pt
m
RRd
To get the fifth line, we used the definition of Vr= 2 is the global minimum of V," so that (V+ = To get the second last line, we used the fact that
20h I 5dM
Dyh 1/21
AV = =1 (1+h 2) h- s (1+hh d.,
Proof of Lem. E.3. Note that 5;- does not depend on h. When + oc, as + (11), we
know Zv; converges to Zu: by the definition of p == e U2 we immediately know that
lim | log Zv. = O. h+000
Then
5d(1 MM
lim &
h
2
Moreover, V1 + X, as - 00, we have
20dM 55
lim C.- lim &t + lim (D. - xx-1+
h
m
5d(1 + M )
(50 + 35,M,-)dMM;;
+ 2
| 12
2
m
E.3 Proof of Prop. 3.5
Case I: h 0. For this ODE case, (8) becomes
1
(1)
(39)
and the formula of leading order term T.(0) is
L(0) = - / (z) 2 Rd po(x)
1
T
dx,
VT
( (p(=)) dt.
2 T-a
23

---

since we perturb the o,- only when t ~ T, that is, 6, (x) = [e (7) for small positive Then we know that
1 T V (poli) dt 2 a V. (pol))
UT
2 T-a
Recall that Po = p 2 pi when R T; moreover, iT ~ $ (0) term of I((0) when < 1 is
T,T = Id. Therefore, the leading order
a²
(V (D)
I(0) 8 R Po
Case II: h + cc. For the SDE case,
1 h
UT
①
(h) 1,
V (p, Ex-))
dt
2
1+h
4 1T
(h)
V (po dt
2
T-a
1+
T
(T 1 KT
V (pnE)) dt
(by a < 1)
2
T-
1 h2
c Kt
v (po de,
2
where
KF Ap V (VV, 14).
(40)
Then when we further assume h
>>>> 1,
(UT)
(UT)
(by Lem. E.I (iii)
Po
Po
(1 h222
e-a 22 No
(V. (p)
(by Lem. D.2)
8
h2 Kin
Rd
Pa
(1 + h2)2 (1
e a 2 no)
(V . E))
h4
2nd
Rd
Po
(1 e 002/
F))
2nd
J2d
Po
(1 e- 2 No
E))
2xD
Po
To apply Lem. D.2 in the second line, we remark that KF (4) = Ap + V (VV n). V has Hessian lower bound Kt = Kn (16).
E,4
Discussion on the pulse-shape error
Though the error of the score function time-dependent, we can always divide the error function into a linear combination of step functions:
(x) E
where X is the indicator function and Øj is the function value of ε on the time interval tj.l- Without loss of generality, assume that the time-discretization is uniform, i.e., -1 = A/ for any j. If we don't worry about ill-behaved functions, this decomposition can be made arbitrarily accurate by choosing a smaller time interval; it is not difficult to make this approximation mathematically rigorous. It is hard to assume that this error function at time c.g., t 0.1 has some subtle connection
24

---

with its eтrоr at e.g., t = 0.8 in general; how the error function exactly looks like depend on vast amount hyper-parametery
training. Therefore, we might as well treat each oj as "mostly independent" in order to handle the worst case situation. If we consider how each Φ contributes to the final sample generation error (quantified via KI. divergence), then we might as well study each one independently. We remark that this is simply reasonable assumption in order to treat generic error types.
More specifically, recall that in Prop. 3.2, VT linearly depends on each Øj
(h' UT /' (1
(pf gT-)) di (n Φz (pt,o)) Al,
~ E (1 + 2) ky.
(h*) V
and recall that the leading order term (10) is simply
(z) dx.
L(h
2 Jgd po(x)
Therefore, after plugging the decomposition o, into L(h*), we have
h22 (At)2
(h
Pr. Ø3
T (h-)
V
L(h*)
8
J.k
We had chosen h,- = h for all t E [0, T] for simplicity. Either by Cauchy-Schwartz inequality, or by assumptions on the independence of each o, as discussed above, the main feature is the following term
+ h2
④ (h
Pr 0; 0,
At
V
add Pu
Therefore, it makes sense to study how this quantity scales for each j. and for two asymptotic regions h 0 and h + oo. For each fixed j. we can easily observe that the above quantity arises from choosing the error function as pulse-shape error, i.c., ΦE if + j (with certain normalization rescaling).
Proof of Prop. 3.6 and the discussion of the large diffusion limit
Proof Prop. 3.6
Let us pick an a = h & 1. We split ür (10) into two parts:
1
T
(h*)
UT
Tb L.T
V
(pi Eú
dt
2
=:
1 h
T
T (b+)
(TT) di
1 h2
T T-a
T (h-
((TT))
dt.
2
i.T
T-o
25

---

Upper bound
By Cauchy-Schwarte inequality ((x + y)² x² + x² 2ry S (1+ a² a²)x² + +1+ a 2)y) for any
I, and a >> 0),
T(h) 1 2
2 Po
$ (1 + a2) (1+h?)²
(S7
1// (T )dt)²
R2
po
=
(1 + h2)2
(So (1'-a (h-)
(T )dt)
+ a1+
1,
8
Rd
Po
=
For the second term T2.
T- o (h" ) (TF )dt)
(1 h2)
T2
8
Rd
Po
Cauchy ineq
(1+ 2
a)
T T-a
$ 1.T
(hT)
(TT )2
dt
8
0
Rd
Po
(29)
h
T
$
(
LT
(TT))
dt
Rd
Po
(1+ h2)2T
T-a
r+-
(-1," ((h2 - de
exp
2
R
Pt
(1 hT
T
00 + Cr (r) frre -Pa(eT-- dt
sup
exp
2
1E[0,71
Rd
Pt
(1 h2)T
r+2
(2 + +C, c. (1 + h²)e
sup
exp
2h
(E(0.2]
JRd
Pt
which decays exponentially fast + oo as long as h² >> a1 > 1. To get the second line above, used Souchy-Schwarz incquality and Fubini's theorem. Recall the definition of y Rt in the statement of Prop. 3.6.
Overall, when h 1
L(h)= < (1+ 2) c (1 h2) h2) exp (- (-h2-87)
+ (1 (1+aa)TT
7
where
C
(1+ h2)
r tE(O,T] zd PE
sup
) exp ( =xppl
() + (00
2h
< T
exp 1-d 10 + lim CF h
sup
2
(E(T.T]
g pt
Please refer to Lem. E.3 for limp (1)
Lower bound
The lower bound can be proved in the same way: by Cauchy-Schwarz inequality again,
L(h) 2 (1 1 aT)TI + (1 a1 22
N (1-a?)TT (a - - 1)C (1 + h?) exp ( h22 hh2007)
The remaining task is to estimate T
26

---

Asymptotic limit of Ti
We only provide an asymptotic result below. As < 1,
(JT- 1 (A- (TT)de)
T h1)
8
JRd
Po
(1 + h2)
)
dt) (TF))
8
Rd po
T-a
(1 + h2)²
h2
| dt) 2
cxp
(T t)KF Ir
8
2
T-a
+ h2)²
T
h2
()k dt) di)² ?
exp
8
Rd po
c <h> h²
where we used 17 2 17, when 2 T the second line; we used
when t T in
the third line; L7 (hh) h2/2 KG (40) herein. In the last line, we used the fact that po Po when 00.
By Lem. D.2.
62 X
(1+
2 (h)
T
a
8
h2 (h) A
k=1
(by h 9 h >>>>> 1).
2
(h)
where (A(h) Lk are eigen pairs of Ki- and F E ok ok When h oo, know that
(h)
(h) (h)
the eigenvalues of KGT, which depends on h and has the form
K (z) = Ap(x)= + V (VUo(z))(z))+hr +h ((VV())) - -2)
should converge to K defined as
(m) = Au + V (VUlH).
Therefore, in the limit,
lim T
Σ
/x(
0-000
2
k=1
where (X(00), (lo0)) are eigen pairs of 1000 and I =
(oc) (00)
0660
kk=1 0
Upper bound of Ti
By the upper bound in Lem. D.2 (together with h- directly to the asymptotic of Ti, when h >> 1,
and h >>>>> 1) or by applying X
2 Ko
r
1
(pT E1 ))2
(V.
(poo)) ))²
T
2kD
JRd
Po
2ND
JRd
Po
2mg 2
Ju
Po
The term T Prop. 3.6 is simply the limit of Tí.
F22 Remark on the large diffusion limit
Score function parameterization in constrained score models
We would like to remark on the parameterization score function. is common in literature to directly parameterize G, (4) via some neural network. However, in general, this practice can
27

---

rarely guarantee that 6, (which is supposed to approximate V log p.) has the gradient form as well [23, 25J. If we instead parameterize log p. (x) N, where (::0) is some (scalar-valued) neural network with parameter 0, we can simply use the following ansatz. in training
(z)z) (2;(2;)).
Consequently, the error of score function estimate also has gradient form:
E6 (x) ,, (2) - V logg (a))
= VN, (z; 0) log PÊ (z)
(N, (z; log P. (x)
This validates the discussion in $ 3.6. Apart from the apparent benefit in preserving the gradient form, this ansatz can help us identify some interesting structure for the the upper bound of lith- L(h) in (18), discussed in 3.6.
More discussion on the upper bound
Because the score error e comes from approximating V log Po. it is reasonable to consider the case = log Po. the upper bound (18), or equivalently (19), becomes
(41)
2m Rd
This quantity heuristically characterizes how difficult the probability distribution Pu can be learn by score-based diffusion models with large diffusion coefficient.
of interest to further explore how to utilize the above upper bound improve the training of score function, provide theoretical understanding about what types of probability distributions are easy/difficult to learn via score-based SDEs. We shall leave an extensive study of these important and interesting questions for future work. Below, we shall provide simple example.
Example: d-dimensional Gaussian
When po -Uo [M, og), the quantity mo (14) my
1. oB
Then this upper bound (18)
od
d
-o =d.
This is compatible with the insight that more difficult learn a high dimensional probability distribution.
G Numerical experiment: 1D Gaussian case
We demonstrate and validate main findings via simple ID Gaussian. Suppose f(()) - and the exact data distribution is a Gaussian po = (0, oū). The error is quantified by KL ( We investigate this eтror as function of the magnitude of h (which is chosen constant function h h for any E T.).
G.1 Explicit formulas for 1D Gaussian case
From solving the SDE for the forward process, we know that
Pt NO, of),
of = ofe + 1
Vlogp( = x
(42)
of
The backward dynamics (7) on & [0, is
1 1+
Y
dř, 2
(
+ c6 (8.(()))) dt + hh dWe
o
law(Yr) = PT.
28

---

In this example, there is only one source of error, which is 6", as we know explicitly the distribution Pr and we can choose the time step small enough such that numerical discretization error is negligible. Due to the structure of the score function, it is reasonable to consider the following ansatz.
C,
(43)
This example has explicit formulas: Yr is Gaussian distribution with B[Yr] = 0 and
Var(Yr) GT Var(Po) + S' GT2 dt,
where
1
h
1
G exp
+ EaL. ds).
2
07
The sample generation error quantified by KL divergence also has an explicit formula:
Var(Yr);
o
KI- (p-(pol) = 2 log
2Vrr(Yr) 2
G.2
Experiment 1: Fix error magnitude and consider various error types
We choose T = 2, h = for all E [0,7], 7), and the error function chosen as
case
-1
case 2;
+ sin (2
0.02,
(x) V llggp(())
case 3;
(44)
2
case 4;
I29.9997
case 5.
For these choices, we only perturb the true score function by a bounded prefactor. We can observe that the eиror (L(pollğr) is a complicated function of general in Fig. 6. When we only perturb the score function during the initial period of the generative process (case 4), we can clearly observe that the sampling error decays exponentially fast with respect to h2, which numerically validates Prop. 3.4. When we only perturb the score function near the end of the generative process (case 5), increasing the diffusion coefficient will actually increase the error, which is predicted by Prop. 3.5. For general error (case 1, 2, 3), overall we can still expect that increasing diffusion coefficient will generally suppress the error when do small.
Table Approximated value of L(h) when h2 Rs 20. can observe that in each column (referring to specific error type) the value is almost independent of the do in particular when < 1.
Case
Case 2
Case
0.2
0.2567
0,3032
0.0658
0.4
0.2569
0.3028
0.0636
0.6
0.2570
0.3018
0.0597
0.8
0.2569
0.3023
0.0544
1.5
0.2570
0,3017
0.0321
2,0
0.2564
0.3022
0.0198
3.0
0.2566
0.3020
0.0121
G.3 Experiment 2: consider L(h) for various error types
We further numerically approximate L(h) = L(h, 6- by linear regression and study how do (i.e., the data distribution Pu). error type 6 and diffusion coefficient h affect the leading order term L(h, _ .Po). Recall that we use L(h) as a simplified notation when d and pa are clear from context; see 3. As remark, approximate L(h). used the leading-order approximation that KL(pollãs l L(h)e² + Oe²): we choose few values and use linear regression estimate L(h).
29

---

In Fig. 7, we can clearly observe that L(h) converges to a constant extremely fast when h increases for cases dg < The value of 6 for the case is only 1/2 of that for case and case 2. By Prop. 3.6, we know that limp- T.(h) for the case 3 should be approximated 1 1/4 of that for cases and 2. This is also (approximately) numerically observed in T. 3.
-3.5
-4.0
-3.0
-3.0
-4.4
-366
-3.5
-6
-4.0
-4.5%
-4.8
-3.6
-4.0
-3.0
-3.2
-4.0
-4
-4.8
335
-366
-5.6
-4.0
-4.0
20
20
20
20
-3.2
-3.2
-4.0
3.6
-3.6
-4.4-
4,00
20
20
20
5
20
-3.3
-3.3
-4.0
-6
-3.9
-3.9
-4.8
-7
-5.6
-3.75
-3.7
-6.4
-3.00
-3.90
-72
20
20
5
20
-3.8
-3.8
-4.4
-6.0
-3.9
-3.9
-3.9
66.
-4.0
-4.0
-5.2
-4.1
-7.0
20
20
20
5
10
20
-4.00
-4.00
-4.00
-6.8
-4,8
-4.25
4.25
-4.25
7.00
-5.1
4.50
-4,50
20
15
Figure 6: We show log10 (polpol)gr) as a function of h2 for the ID Gaussian model. T = different do and error functions in (44) are considered.
H More details about numerical experiments in $ 4
In this section, we discuss datasets, network architectures, evaluation metrics, numerical schemes (exponential integrator), and the default weight in denoising score matching.
H.1 Datasets
1D 2-mode Gaussian mixture: Po(x) = Σ 0.5N (- 0.01).
2D 4-mode Gaussian mixture: Po(x) Σi (1: ((-1000))((1111)))) (..0)))...052
Swiss roll: Swiss roll generates samples by (z, v) (tsin(t), cos((t) with t drawn from the uniform distribution U( 3x 3
9x
30

---

case
case
case
0.0
0g 0.2
0.6
0.5
-0.
0.0
0.0
-0.5
0.5
0.4
dg =04
-0.4
0.0
-0.8
-0.5
-0.4
1.2
10
0
10
10
0.4
0.3
-0.4
0g = 0.6
0.0
0.0
-0.8
-0.3
0
0.38.8
0.0
0.00
-0.8
-0.3
-025
-1.2
-0.6
-0.50
10
0
10
10
-0.8
00.30
-0.3
0.5
-1.2
-0.45
-0.4
-1.6
-055
-0.4
2.0
-0.45
-1.2
-0.5
-060
-1.6
-0.6
0.50
-1.
3.0
-0.75
-0.5
-1.5
-1.00
-1.00
-1.8
10
10
10
Figure 7: We show logno (L(h)) as a function h2 for the ID Gaussian model. and error functions in (44) are considered.
2, different do
31

---

MNIST: MNIST [34] contains 60,000 28 28 gray-scale images with hand-written digits. CIFAR-10: CIFAR-IO [22] contains 60,000 32 x 32 RGB images with ten categories.
H.2 Network architectures and other parameters
For experiments on 1D/2D Gaussian mixtures, the exact scores can be obtained analytically if we use VP-SDE. We set I' = 4 and 9. for t E [0, T]. For the time discretization when solving the reverse SDE with Euler-Maruyama method, we apply 40,000 steps and 80,000 steps for ID and 2D Gaussian mixtures, respectively.
For experiments on Swiss roll, we apply a three-layer neural network for score matching. where the width of each layer is set 50, 50, and 2 and apply ReLU the nonlinear activation for two hidden layers. We set 7 = and g: / 01 - 1) + 20/ for / E [0. 7]. The learning rate is set as 0.01 and decays by 0.5 every 8,000 steps. The batch size is set 400. We train the neural network for 20,000 steps. For the time discretization when solving the reverse SDE with Euler-Maruyand, method, we apply 20,000 steps.
For experiments on MNIST, we apply the net architecture in [17] for score matching, where we use two resolution blocks in U-net and set the multipliers of channels be one and two. We set T 1.4 and St = (1 1) + 20/ for E [0, T). The number of iteration is 20, 000, the batch size is set as 64. We solve the reverse SDE with exponential integrator; see also Appx. H.4.
For experiments CIFAR-10, use their pretrained checkpoint in for score estimation. We use the same setting as [30], i.e., T = 1 and 9: = 0.1(1 - 1) + 20/ for / E [0, TJ. For the time discretization when sampling with Euler- Maruyama method, we apply 100, 200, 500, 1000, 2000, 3000 and 4000 steps.
apply the DDPM++ cont. in [30] as the net architecture, and
H.3 Evaluation metrics
For experiments of 1D/2D Gaussian mixtures and Swiss roll, we apply approximated divergences for evaluating the performances. Specifically, we discretize the space into 100 bins in each dimension, then obtain the empirical densities of 10,000 true samples and 10,000 generated samples, and use Jensen-Shanne divergence, Bullback-Leible divergence and Wasserstein distance between both empirical densities as the metrics for evaluation.
H.4 Exponential integrator
We shall explain the exponential integrator for (2) with t E [0, T), i.e., the following equation
dY, ff- (Yi) dt (9. )2
2 V Vlogp; (V)) dt h aw.
= 96 2,, + (gt 12+ (hT))
Voggt (Y,) dt + ht dW..
2
Then for any time & E [A,, 4++1 1.. and given Yı,, we approximate the above dynamics by
a² ~ 99 28, di + (g + (h
V Voogii (Yr.) di + hy dW.-
2
This dynamics is linear SDE and solve it exactly
[1,, flat et 1.++1 9 do 91 2 h4
Y ē Sit at aey,, +
2 dl S
=1
Ne+s
e Sittt gi 2 dW,,
+
1x
:=9s
S = V V Prg (Yte).
32

---

Therefore,
Sma. (F2) dt Zk.
Yee 91 + V
Zu ~~~~ ~o,,
we pick
h = agi
a E R++,
9: v B0 + (B1 Bo)t
which implies that
gt v Bo (B1 Bo)(T t),
then after straightforward calculations,
PeNI = Vor + (1+ -1)s + vae(1) -1)22 -k=
where 0g = 16+1 4x and 7/k = exp 6n (280+(211 -27 27'-2------ 16.0 ( 21l)
H.S
Default training weight
Given the fixed Xo. explicitly solve (1) (with the choice fe(x)) a:):
X, = (G,) Xo+ [ G. GT dW..
where G, exp (Jo oh 2 ds). The standard deviation of the Brownian motion term above is
V.. G. 92 ds.
When 9: / + (B1 foot used in many literature [30, 17], explicitly solve for we;
V 1 exp
1 -- (o) - 1,30
In literatures, to balance the noise over time in training the score function, a common practice is to use the default weight W/ = wu (5) in the (denoising) score-matching loss function (4).
H.6 Additional numerical results for experiments in
We present more numerical results on ID Gaussian mixture, Swiss roll, and MNIST to further verify theoretical results.
ID Gaussian mixture
We evaluate the performances of generative models under different values for 1D Gaussian mixture, and present the visualization and numerical results in Fig. 8 and Fig. 9, respectively. In Fig. a clear trend shows that with increasing h, the empirical density of generated samples better matches the true density function. This trend is more quantitatively captured in Fig. 9, from which clearly observe that the distance between the empirical and the true density function decreases to the numerical threshold exponentially fast. Numerical threshold means the error of various distances when O; due to the space discretization when computing various distances, the numerical values of various distances are not exactly zero even when we use the exact score function. However, increasing h can help us to almost reach this limit and this phenomenon is theoretically described in Prop. 3.4.
33

---

h? g?
h² h2-1.0.
h2 2.0
h² 3.0
4.0
(a) E -..2.
h2 0.0 g²
h1=100 g²
-202.0
h? 3.0 g²
4.0 g2
3.5
130
3.0
17:
130
130
300
(h) = - 0.1.
h² 0.0
11=1.0g
3.0
h2 = 4,0
345
130
(c) 0.0. h²=2.0
= 0.0 g²
g2
=.0 g²
4.0 g²
43
24
(d) e 0.1. -20
h2 0.0
h2 3.0 g?
h2 4,0 g2
34
140
-5
45
-5
(e) 0.2.
Figure 8: Visualization of ID GMM, where (x) v log (x). We can observe that increasing h can help to generate samples with a distribution closer to the true Po-
34

---

divergence
divergence
KL divergence
50
-02
-01
003
20
10
04
1''9"
P
(b) 6
1
log
(c)
1c00.95 V log
2
Figure Numerical results of 1D 2-mode Gaussian mixture. We can observe that the distance between the empirical density and the true density function decreases as h increases, when &f is not extremely large.
Swiss roll
In Fig. 10, we provide additional figures to discuss the effect of time-discretizatinn When the numerical error is negligible, we can observe that the generative process with a larger can provide clearer picture of Swiss roll, as shown in Fig. 10c. However, when the discretization error cannot be ignored, the conclusion may be reversed. Therefore, is necessary to design and employ more accurate numerical methods for models with large h in order to fully benefit from diffusion models with large diffusion coefficient.
In Fig. 11, we display JS and KL divergences between the true density Po and the generated samples. Since the data distribution of Swiss roll highly localized (data concentrated on curve embed in 2D), accurately computing the KI. divergence poses significant numerical challenge. That is why we use Wasserstein distance instead Fig. 5a. Based on more robust symmetric metrics (i.e., the IS and Wasserstein distance herein), we can observe that a larger h can indeed diminish the error in sample generation.
h 0.0 g?
h²= 1.0 g²
h² =.0
=²03.0
4..00
(a) diseretization steps 100 hº = 2.0 g² h² g²
0000
1.=
4.0
(b) discretization steps 1,000 ht 2.0
020
10=00
h² 3.0 g²
h2 4.0
discrctization steps 10,000
Figure 10: Visualization results of Swiss roll with different number of time steps.
35

---

is divergence
diverg ence
5.0 x10-
x10-2
x10-1
x10-2
3
h9/g"
h2/g*
Figure 11: Numerical results of Jensen-Shannon divergence and Kullback-Leibler divergence of Swiss roll.
MNIST
We conduct further experiments explore the effect of ht and the time discretization steps. is important to note that if one has an extremely well trained score function, then the effect of h is indeed negligible, as shown in (3). To somewhat magnify the score training error for MNIST (but with reasonable score function), we increase the time T and use smaller architecture with fewer parameters; the reason of the occasional failure to generate clear MNIST images later figures comes from this deliberate experimental design.
For an existing pre-trained score function with non non-negligible error, we notice that increasing h can almost ubiquitously improve the quality of sample generation; see Fig. 12 (as well as Fig. 15 and Fig. 16 under different training setup). This improvement is supported by our theoretical results, particularly Prop. 3.4. Notably, even when choosing h = 0 (ODE) and ht g" (the default diffusion model many studies) occasionally fail generate an image, simply by increasing the magnitude of h- (possibly at the cost of more computational resources), we have larger chance of generating an image with reasonable quality; for instance, see the last row in Fig. 12 (particularly see the last row in Fig. 15). The ODE-based model sometimes fails to generate hand-writing digits even when the step number is 103, whereas the diffusion model (with large h²* does not encounter this issue. This conclusion might be reversed when the time step size is large, which is similarly observed in the Swiss roll.
The computational cost of SDE-based models with large diffusson (h* > gg*) is relatively high due to the necessity of larger number of time steps: larger h² has the similar effect as running Langevin for a larger time horizon as discussed in $ 3.3 and longer time simulation is expected be more expensive and also its accuracy largely relies on a well-chosen time step. However, this can be offset by the ability to use lightweight architecture that possibly speeds up the generative process. A detailed comparison of various diffusion-bart models with the same computational budget constraint is challenging and is slightly beyond the scope of this work, and we will leave this task future work.
36

---

weight type=default, trial=1
=
h²/g2= = 0.5
b²/g²=1.0 1.0
h²/g²= = 1.5
h21g² 2.0
step 30
r 7
1999
6 4
5
8
step 60
A 7 0 2 R 9
/ 6
2 5 92
2 1 /
7 2
2 1207
7 2 7
G 2 d
V
23
D 2 005094
D 0
A D 3663
C D
F 4 4 3
8 4 P3 P
19 1998440 16366111
2 D 4 2/44
2 7D944724
2 1ED00/10 17011220
A 2200 7
4 84482222 2
206011220
8969:
19622 9 F
V FE3
O 4 4 5
3 F 5555 6 3
O 9 & 6
step 100
7
A
/ 2 $
O 7 2 / 2
R S é A 6
9 7
O
A 3 9 1
F A 9 9 A 4
K 4 9 8 2
A
è A 4 f
A 2
A 8 9
denei 0
8 07763
PED44420
/6997/40
" 5 4
/ <
50000
5 4490
12
X
5 4
E 8 0 5
7
2 274
2 / 6476
2 d S F 2 4
2 9 L 2 6
step 200
8 6 2
D 32
0 4 3
O
A 4
C A
1 9 102
G 8 2
4 A
R 6 6
R9 69699999
2299 290 19901
6 1
1 2
7 F 1 13334150
2 0
2 0
4
P 222229J0
C P 122200900
C ( 0
7 3
S 6 3 6
a Y 5
step= 1000
2
A O 9 2
C 2 2
21
6 W 8
YYYO
2 7 7 0
8
4
O O
7 7
& 8 7
7
7
5
5
382
5 7
L 200957
6 9 0D0059
D022
2 8
2A 3
2 2 83884400
8
Figure 12: Visualization of sample generation in MNIST where we used default weight w (5). We use different random seed to ensure robustness and this is the result for trial 1. Results for trial is not included this work as the overall tendency remains the same.
37

---

I Numerical experiments for adopting different weight in training
We have reasons explore the effects of various weight functions in training.
The first reason that theoretically, any positive scalar-valued functions w/ on (0. T) is a valid candidate. This prompts the question of whether such default weight function (5) is optimal in designing the loss function. We fully acknowledge that the default choice adopted most literatures is very effective one. However, the mathematical reason behind it is still not satisfactory in our opinion. This motivates us to ask whether the default choice is really the optimal one, at least in certain circumstances.
The second (and actually the primary reason) comes from our theoretical predictions discussed in $ 3.7. As there is a tight connection between the optimal reverse-ting generative process and the iimediitributtin of the score error (see Prop. 3.4, Prop. 3.5 and Prop. 3.6), if we are willing to train re-train the score function and are interested in using the ODE-based model (for fast sample generation), it appears that we should focus more on the noise's s end comparatively. To achieve this goal, namely, control the distribution of score error, we adopt different weight schemes in the loss function only for training: default weight (5), data-driven case (more weight in the data side) and noise-driven case (more weight in the noise side):
noise-drivens or simply referred as "noise""
W
(45)
data-driven, simply referred "data".
0.25 + w.
There is no theoretical reasons behind the noise-driven and data-driven choices in the last equation. We merely experiment with reasonable choices and at the same time they expected to help us control the time distribution of the score error.
After the initial camera-ready version of this work, notice that such weight design (in particular the noise-driven weight) has also been studied in [12]. Their experiments seem to provide further evidence about the validity of our arguments. However, we would like emphasize that the approaches and perspectives towards proposing such a weight design are different: in [12], their conjecture is based on an insightful perception argument, whereas take mathematical approach and come up with the design based on the above error accumulation analysis.
1.1 Our guess
Based on our theoretical results, we guess that if a different weight function can really achieve our expected goal (namely, control the time distribution of the score function), the noise-driven case should be more suitable for ODE models, and the data-driven weight should give us the worst performance for ODE models.
We remark that this conjecture is surely not universal, and its validity remains to be fully validated by more benchmark experiments. Nevertheless, our numerical experiments below suggest its potential usefulness and it prompts an interesting question to explore and design the optimal score-matching loss function, which is rarely studied in literature. In what follows, report numerical experiments for MNIST and CIFAR-10.
1.2 MNIST
We used different weights to train the score function and then visualize their generated samples in Fig. 13. We can clearly observe that the score function trained by the noise-driven weight produces comparatively better samples in ODE models. We plot the denoising score-matching loss for score functions obtained from training with various weights in Fig. 14:
to make the comparison more straightforward, we visualize the dine-distribution of relative score-matching loss rather than the absolute value, namely, we demonstrate:
[[())() V log Pr (xX,
Exo Ex. Xx)
+
i E {noise, data],
E Expr E Ex (X))
(ddefault) ) (X,)
log (X.,
(46)
38

---

and () refers to the score function obtained from training using weight i and we use test datasets to approximately represent Pn;
to ensure robustness, we independently conduct two trials with different neural network initializa- tions: For each trial, the initial neural network the same and the only difference in training is to adopt different weights in the loss function.
In Fig. 14, we can clearly observe that adopting different weights indeed help us to control how the score error is distributed over time as we expect (e.g., noise-driven weight helps us reduce the error near t = T and has an opposite effect near f = 0); as shown in Fig. 13, their numerical performances terms of sample generation also match our guess above.
We further visualize how the diffusion coefficient / and time step size affect the sample generation quality for score functions obtained by data-driven weight (see Fig. 15) and noise-driven weight (see Fig. 16). The conclusion is the same as in the default weight case. This further validates our theoretical results in 3.
weight type = data
weight type = default
weight type noise
3
3
8
4
4
7
3
8
1
B
8
6
7
Q B
8
0
5
8
9
/
7
8
5 2
D
a
6
2
9
3
6
2 A
9 0
5
2
9
0
7
6
7
3
4
7
9
7
a
9
3 0
5
s
4 Q
8
S
6
0
8 8
G
i 2
8
S
8
010
9
2
9
8
5
0
8
9
A
71
9 &
5
0
7
2
2 J / 7
2
83
/.
f
6
f
7
o
O
9
3
o
9
4
6
3
4 9
0
7
2
8
5
6
9
Q
9
6
9
0
7
9 9
9
0
9
0
5
9
5
5 200
(a) Trial number
weight type data
weight type default
weight type noise
4
4
6
9
E
6 4 0 0
L
U
4 6
U
0
//
U
0
6
3
5 4
4
3
$
2
4
3
90
9
9
6
'1
4
7
1
6
75
O
4 6
3
8 0
+
90
i 6
3
9
4
1
93
d
0
7
5
1
/
6
0
5
4
?
1
9
k
5
$
1
4
4
5
7
9 2
0
0
9
6
2
9
I
5 5
a
6
2
0 3
3
5 8
4
9
6
3
4 c
2 8
6
4 8
8
2
9
7
A
L
7
1
A
2
3
3
7
6
7
3
3
Q
3
6
7
7
9
(b) Trial number
Figure 13: We show generated samples trained for MNIST using three different weight functions in the loss function: we use k²/g" 0 (ODE), time step is 1000, the trial number is the index for independent random initialization.
39

---

Relative score-matching loss
data, trial noise, trial data, trial
1.6
1.4
noise, trial 2
1.2
reference
1.0
0.8
0.6
0.0
0.2
0.4
0.6
0.8
1.0
t/T
Figure 14: Score- matching loss for data-driven and noise-driven on MNIST, compared with the default case
Score-matching
40

---

weight type= trial=1
=
h²/g2= 0.5
h²/g²= 1.0
h²/g² = 1.5
h21g² 2.0
2
step 30
8
8
Kepp 60
/
>
2
/ 5 y
6
8
b 4 48 6
7 6 & $
2 1729940
2 2 Q 11
24
L
V & F54 3
O 4 65278
9 C
$.4
step 100
3 M 9 4
/ & 9
O R I 7 7 2270
2
O 2
4 2 119199 11919A1 19911
A 9 1
A 7 3 1
A 2 4
A 8 171112
A 60149
4366
S 1000300
1770
0 0770
Q 077
CC97744VS
" 1661517 5445
J 6 1 KUS
b O 119011211
& 6
6 D 3 /2
F I 11199709 77770989 D
0 7 097779 D
E
D
7
5 2 6
2 S S 6 6 A
2 J 9 6 ML
step step = 200
699
D 4 8
A O 2
D A
E
G 7
C 9 3 93
2 62966487
2 10099430 519
2
4 9 9
E7701314
L 07703555 4
A 2229944
< AJ A.A-
C- (-LF)3/0O 0
N
Y
step=1000
A t 0
P 9 K
C FI
O H
4 > 87 2 7 9 2
2 L 7
L 7 4
6 7
6 8 4
7
7040
2
2
2
F 00050
9 / D
A DC
C 05
7 P 4 3
2 3
9 A.A.A1 e
Figure 15: Visualization sample generation in MNIST where we used data-driven weight w (45).
We use different random seed to ensure robustness and this is the result for trial 1. Results for trial is not included in this work as the overall tendency remains the same.
41

---

weight nope= trial=1
=
h4/92=1.0 1.0
h²/g² = 1.5
h3 =20
8J & 0 7
step 30
7276679
2 8 8 0
72009 9 2883
99900
6.
3 1 7 2
F K 8
1892
5 6
step 60
5 1
q 3 0
2 2 8
0 962
/ 5
/ 5
2
7 2
120077
7
7 az
3
8948
A
D 0539
O 0
004029 8/3
4 84
V O 4
7 F 4
11900 8875
18044000
2 29994
2 S 15299 R
1200
148287741
8 R 8 94457
3 6 V A
O V 5
OO7 5
9
6
9
4
step Step=100 100
S 9 2
9
2 G 7 f &
J
7 2
7 2 8
/ 4 8
9 7
7
8 0
4 3 9 2
E A D 96100
A D D 2 9
2 2
3
283
è A 4 7
ARAF882 22172
8
A
54644
5900T
300
d99e 2996
6 77703
8 69
2 DED4
4 7694774VC7
" 1 4
5 5
5 50 50997 14413
5 8 1.8
8
7
8 8 & A 4
9 7055 5
F 77689 778977950
9 97797799 5
7
2 3 1 2
/ 9 2742
O /
44
O J 4 6076
2 5 V L 607
step step=200 200
& 0 27
8
2
D 862121
4 O 000939 3
D 4
4
S 28
4 A3
4 0 V 8 8
9 DVODSTTO 17
G 02409277
4 A 9
R229 722963814 33341
2 029963898
2099 137V1
2 2 3 7
777035M4
2 07/0330 0703540 4464
2 0 4
9939
2
C G4P322 P4PP30qU0
C -L-28 20030000
A7/3
A?H2
20
3
3 1
S 6Y
8
step= Sxep=1000 1000
9 2
A 0 2
F 17044430 4
C 4
4
& & 2
6 9 17
2 4 7 R 1
7 0 7 0
6
6
0
C 4 3
7 9767
6
4 7
V
4
7 54
7 17ADD88
7
577777
0
4902
7 1 202240000
/
7
9
9
9 200957
9 20095
6 D095 6
DO25
5
2 8
2 6 3
6 4 4
L 8
Figure 16: Visualization of sample generation MNIST where we used noise-driven weight w (45).
We use different random seed to ensure robustness and this is the result for trial 1. Results for trial is not included in this work as the overall tendency remains the same.
42

---

1.3 CIFAR-10
We carry out a similar experiment for CIFAR-10 10 to test the effect of training weights wi: we used the same initialization (referred to as the trial number below) and all other hyper-pee parameters, except that we employ different weights in score- matching loss (4). The same architecture is used as [17] for CIFAR-10 and the detailed hyper-parameters can be found in source codes.
In Fig. 17, observe that overall over the whole training period, the noise-driver weight leads into score function estimate no worsen than that by the default weight: due to stochastic fluctuations and other uncertainties (in particular we adopt the mixed precision training). there no guarantee that the noise-driven one is always better, but the overall tendency is still observable and clear. In Fig. 18, the noise-driven one actually has a slightly larger score-matching loss (SML) than the default one (which could be possibly explained by how we measure the SML in Fig. 18). What is interesting is that score functions trained by the noise-driven weight and the data-driven weight have similar SML, which both decay at the similar pace: however, the FID values for the score function estimate by noise-driven weight are much smaller than that by the data-driven weight. This apparent gap clearly explains that apart from the total score-matching loss (which does matter), the time distribution of the score error plays an important role in determining the final sampling error, echoing our theoretical results in $ 3.
For instance, in Fig. 17, if we consider the first experiment (i.e., trial=0) with floatl mixed-precision training (i.e., mixed=True), we notice that relatively near 80k and 120k training iterations, the performance of noise-driven one much better than the default one, which is consistent with Fig. 19 that the relative loss near t 2 T is more minimised for iterations 80k and 120k, compared with other iteration stages. Moreover, for the same experiment in Fig. 17. the data-driven one has a much worsen FID value at iteration 200k, which is compatible with the increasing relative error near the noise's end (i.e., f 7) in the last row of Fig. 19. For the remaining three experimental setup (either different initialization or training precision), we notice similar consistency between how the time-distiibution of the SML behaves and how FID values change. This relation, so far, still cannot be used rigorous quantitative indicator to predict one based on the other quantity, but qualitatively, the above explained relationship does appear to be numerically valid and theoretically sound.
In summary, if we have two score functions 1 and Gg from training:
If the SML for G is much larger than the SML for 62, then we can probably confidently expect that 62 a more accurate estimate.
However, when the total SMI. (4) for both are close, then the time-distribution of the score- matching loss together with which generative dynamics is chosen will play a significant role in determining the final sample generation quality, which is probably largely overlooked in current literature as far as we know, A full investigation and in particular whether it is possible to adapt this observation to achieve the state-of-art models will be left to the next stage of research.
43

---

trial=0, True
trial=0,, mixed=
10a
102 1x
6 10
10"
FID
4 x
4 x 10
3 x
3 x 104
2 x
2 x
40k
80k
120k
160k
200k
40k
80k
120k
160k
200k
trial=1, mixed=True
trial=1, mixed=False
data
10
6 10
noise
default
FID
4 x
10a
x0
3 x 10
2 x 10a
40k
80K
120k
160k
200k
40k
80k
120k
160k
200k
We visualize FIDs of ccore-training estimates by various weights with respect to the number of training iteration. We consider ODE models (i.e., h = 0) when computing FIDs. Trial = 0, 1 refers to different neural network initialization; mixed=True means floatló mixed-precision for training, mixed=False means we use 32-bit precision.
Figure 17:
44

---

trial=0,
trial=0, mixed False
77
77
76
76
SML
5ML
75
75
74
74
40k
80k
120k
160k
200k
40k
80k
120k
160k
200k
trial=1, mixed=True
trial=1, maxed=False
77
77
data
noise
76
76
default
SML
5ML
75
75
74
74
40k
80k
120k
160k
200k
40k
80k
120k
160k
200k
Figure 18: We visualize score-matching loss (4) for score function trained using various weights with respect to the number of training iteration's s. Trial = 0, refers to different training initialization; mixed=True means we use 6 precision. When we compare SML for different score functions above, we consistently use the "default" weight for fairer comparison, even though some score functions are trained using different weight functions; we also use the test dataset to approximately represent Po-
for training, mixed=False means we use 32-bit
45

---

trial=0, mixed=True
Relative loss
Loss
train step 40k
noise
10a
1.5
data
104
default
1.0
10a
0.00
0.25
0.50
0.75
1.00
0.00
0.25
0.50
0.75
1.00
train step 200k train step = 160k train step 120k train step 80k
2.0
102
10a
0.5
0.00
0.25
0.50
0.75
1.00
0.00
0.25
0.75
1.00
2.0
10a
1.0
109
0.5
0.00
0.25
0.50
0.75
1.00
0.00
0.25
0.50
0.75
1.00
10a
1.5
10a
10
0.00
0.25
0.50
0.75
1.00
0.00
0.25
0.50
0.75
1.00
2.0
10a
1.5
10a.
0.5
0.50
0.00
0.75
t/T
Figure 19: The first independent experiment with the float16 mixed precision training: We visualize the relative time distribution the SML (46) (on the left) and the time distribution of SML (on the right) for various training weights and for various training stages (each row).
46

---

trial=O, False
Relative loss
Loss
train step 40k
noise
10a
1.5
data
104
default
1.0
10a
0.00
0.25
0.50
0.75
1.00
0.00
0.25
0.50
0.75
1.00
train step 200k train step = 160k train step 120k train step 80k
2.0
102
10a
10a
0.5
0.00
0.25
0.50
0.75
1.00
0.00
0.25
0.75
1.00
2.0
10a
1.0
10a
0.5
0.00
0.25
0.50
0.75
1.00
0.00
0.25
0.50
0.75
1.00
10a
1.5
10a
1.0
109
0.00
0.25
0.50
0.75
1.00
0.00
0.25
0.50
0.75
1.00
2.0
10a
1.5
10a
0.5
0.50
0.00
0.75
t/T
Figure 20: The first independent experiment with the full 32-bit precision training: We visualize the relative time distribution of the SML (46) (on the left) and the time distribution of SML (on the right) for various training weights and for various training stages (each row).
47

---

trial=1, mixed=True
Relative loss
Loss
train step 40k
noise
10a
1.5
data
10*
default
1.0
109
0.00
0.25
0.50
0.75
1.00
0.00
0.25
0.50
0.75
1.00
train step 200k train step = 160k train step 120k train step 80k
2.0
10a
101
0.5
0.00
0.25
0.50
0.75
1.00
0.00
0.25
0.75
1.00
2.0
10a
1.0
10a
0.5
0.00
0.25
0.50
0.75
1.00
0.00
0.25
0.50
0.75
1.00
10a
1.5
101
10a
0.00
0.25
0.50
0.75
1.00
0.00
0.25
0.50
0.75
1.00
2.0
10a
1.5
109
0.5
0.50
0.00
0.75
t/T
Figure 21: The second independent experiment with the float16 mixed precision training: We visualize the relative time distribution the SML (46) (on the left) and the time distribution of SML (on the right) for various training weights and for various training stages (each row).
48

---

trial=1, maxed=False
Relative loss
Loss
train step 40k
noise
10a
1.5
data
default
1.0
10a
0.00
0.25
0.50
0.75
1.00
0.00
0.25
0.50
0.75
1.00
train step 200k train step = 160k train step 120k train step 80k
2.0
10*
10a
0.5
0.00
0.25
0.50
0.75
1.00
0.00
0.25
0.50
0.75
1.00
2.0
10a
1.0
10a
0.5
0.00
0.25
0.50
0.75
1.00
0.00
0.25
0.50
0.75
1.00
10a
1.5
101
10
0.00
0.25
0.50
0.75
1.00
0.00
0.25
0.50
0.75
1.00
2.0
10a
1.5
10a
0.5
0.50
0.00
0.75
t/T
Figure 22: The second independent experiment with the full 32-bit precision training: We visualize the relative time distribution the SML (46) (on the left) and the time distribution of SML (on the right) for various training weights and for various training stages (each row).
49


Reviews:
[{"review_id": "Rh7ymm26q9", "reviewer": "Reviewer_wiU7", "summary": "The paper explores the difference between ODE-based probability flow and SDE-based diffusion models when score training errors are present. Specifically, they investigate how setting the generative diffusion coefficient h impacts sample quality.", "questions": "N/A", "limitations": "- The authors briefly address potential negative societal impacts of their work. \n- I did not see any mention of limitations of their analysis. ", "rating": 3, "confidence": 1, "soundness": 2, "presentation": 2, "contribution": 2, "strengths": "- The question the authors are trying to answer seems to be a good research question. \n- It's nice to see experiments validating some of their assumptions (Fig 1.).\n", "weaknesses": "- As is, the paper is very difficult to follow. \n- Many equations are presented, but few of them are explained. \n- With so many variables presented in the paper, it is difficult to keep track of all of them, making the paper challenging to comprehend. \n", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for taking the time to provide feedback on our paper. We value your insights regarding\nthe presentation of equations and variables. We are dedicated to enhancing the comprehensibility of the\nrevised manuscript by incorporating more detailed explanations and motivations. It’s important to note\nthat this paper serves a technical purpose, primarily centered around mathematical and asymptotic\nanalyses. As such, a substantial number of equations and variables are often unavoidable to precisely\nconveying our intended message, which may require a high level of understanding of mathematical\nanalysis. We acknowledge that we have already dedicated Appendix A to define clearly the most\nimportant variables in the last submitted manuscript, to hopefully facilitate a more streamlined\ntracking of the variables.\n\nRegrettably, we cannot concur with the outlined limitations in the review report:\n + We have indeed addressed the potential social impact in our paper, as we believe it aligns\nwith the submission requirements, particularly when considering the potential enhancements\nto train generative models. As a result, we want to clarify that the negative societal impact\nshould not be regarded as a limitation of our study.\n+  We want to emphasize that we have taken explicit steps to address certain potential limitations.\nThis includes clearly outlining our assumptions, highlighting the interests of incorporating low-\ndimensional manifold information in the summary section (a step we hadn’t previously taken\nin this work), and explicitly indicating that certain theoretical results hold under asymptotic\nconditions. We believe that these aspects could be construed as limitations in our analysis,\nand we have already presented and discussed limitations."}, {"review_id": "SjOa6LEghV", "reviewer": "Reviewer_LAed", "summary": "The paper provides a theoretical analysis of the estimation error of SDE and ODE methods along with some numerical experiments.", "questions": "see weaknesses.", "limitations": "yes", "rating": 3, "confidence": 4, "soundness": 2, "presentation": 2, "contribution": 1, "strengths": "By perturbing the score function, the authors study how the estimation error changes in ODE and SDE methods. ", "weaknesses": "It is commonly known that the sample generation error consists of three parts: discretization error, estimation error, and initialization error, see e.g. [1, 2]. It is also commonly known that usually the ODE works better with fewer NFE and SDE works better with more NFE due to the interplay between discretization error and estimation error. The estimation error might not even dominates the total error. I don't think only studying the continuous time model provides very useful insight into how to choose those two methods optimally. There is a significant lack of connection between this theory and what happens in reality. Please comment on this.\n\n[1] Chen et al., Sampling is as easy as learning the score: theory for diffusion models with minimal data assumptions\n\n\n[2] Lee et al., Convergence of score-based generative modeling for general data distributions", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for taking the time to provide feedback on our paper.\n\nIt's noteworthy that our findings align with certain factual aspects you raised. In the ensuing discussion, we will predominantly focus on two pivotal matters: \n- the significance of refining score function estimation, with the intention to surpass existing conventions; \n- the inherent naturalness and illuminative quality of the continuous-time model within our framework.\n\n### Clarification of our results.\nWe agree with the reviewer's observation that ODE models tend to perform better with fewer function evaluations (NFE), while SDE models excel with higher NFE—a crucial finding. This prompts the question: why does increased NFE potentially boost SDE models over ODE models **whose framework lies in the continuous-time model**, and under what specific conditions? This forms the core of our study.\n\nImportantly, there is no definitive explanation for this phenomenon according to current literature. Our research is fueled by this recognized but mathematically unverified insight. Our objective is to bridge this gap by delving into underlying dynamics, illuminating the conditions that empower SDE models to surpass ODE models with heightened NFE.\n\n### Score function estimation is important\nIt's crucial to emphasize that the estimation error in score training holds comparable, if not greater, significance than the design of an accurate discretization scheme for the inference process. While the reviewer correctly notes that the estimation error may not always dominate the total error, this holds true for numerous scenarios. Notably, situations with limited training data or a preference for lighter-size architectures (e.g., for memory conservation) fall within this category. Score function estimation forms the nucleus of score-based diffusion models, and we believe it's essential to refrain from assuming that training error can be easily minimized or disregarded in most circumstances.\n\n### The continuous-time model is useful.\nWe acknowledge the reviewer's valid point that a continuous-time model does not offer a finalized solution. Nonetheless, despite these acknowledged gaps, the continuous-time model retains its utility:\n-  **Many seminal works are essentially based on continuous-time models.** It's notable that prevailing theoretical studies heavily focus on the continuous-time model or are built upon it. The references highlighted by the reviewer also lean toward continuous-time analysis rather than discrete-time mappings. When we mention the continuous-time model, we encompass both a fully continuous-time representation and a discretized model using a highly accurate numerical scheme. This contrasts with a fully discrete model like DDPM with $10$-NFE or optimal transport maps. Furthermore, we spotlight the seminal work [Score-Based Generative Modeling through Stochastic Differential Equations, ICLR 2021] focusing on continuous-time models. This exemplifies the continued relevance and influence of this framework within the research community.\n- **Continuous-time model is still useful.** If the reviewer acknowledges the significance of comprehending the mentioned wisdom—namely, the interplay between time-discretization error and score estimation error—we believe that a comprehensive response necessitates an understanding of two pedagogical scenarios: (1) absence of score estimate error, addressed in numerous numerical analysis textbooks; (2) absence of discretization error, requiring an exploration of continuous-time models.\nGiven the limited analyses exploring score error's impact on the inference process, our study initiates this inquiry as far as we know, and validates it through numerical experiments. It's crucial to clarify that we don't directly extend findings from continuous-time models to real-world scenarios. Nonetheless, these findings hold significance. For instance, our numerical experiments indicate that for ODE models, the score error during initial inference stages significantly affects final sample errors— typically not a big challenge for SDE models (proved in Sec. 3.4). Prioritizing score error refinement in the noise end could potentially enhance sample generation quality. This insight is illustrated with examples like MNIST and supported by preliminary CIFAR-10 experiments (see the additional PDF).\n- **The perception of discretization error being uncontrollable perhaps has evolved.** Efficient numerical schemes for inference have been thoroughly explored through collaborative research efforts. Notable examples encompass DDIM, gDDIM, and several others. The control of numerical discretization error has evolved from a formidable challenge to a manageable one. Our raised question yet still lacks comprehensive systematic exploration, as far as we know.\n\n### Future direction. \nWe have indeed undertaken some efforts to employ analogous techniques utilized in this paper to comprehend a discrete-time model, characterized by a modest number of function evaluations (e.g., around 10-20 NFE)."}, {"review_id": "h8XW73GpGQ", "reviewer": "Reviewer_GpQ2", "summary": "The authors focus on reverse diffusion process in the presence of non-negligible error in the score function, and estimate KL divergence between the data distribution and the distribution generated by reverse process. They analyze how the KL divergence varies with the diffusion coefficient and demonstrate that a large diffusion coefficient is beneficial when the error of score function is concentrated near beginning of generation process, while a small diffusion coefficient is beneficial in the opposite case. The findings are validated through numerical experiments on Gaussian mixture model and MNIST data set.", "questions": "Regarding the scaling hyperparameter in the loss function of conventional diffusion model, which aims to balance the error at different time of the generative process, can the practical choice of the scaling be associated with the theoretical results presented in this work?", "limitations": "The authors acknowledge limitation of their work. One additional limitation is lack of large-scale experiments to further validate the findings. ", "rating": 7, "confidence": 4, "soundness": 4, "presentation": 4, "contribution": 3, "strengths": "The paper introduces a novel analysis of how error accumulation occurs in diffusion models when the error in the score function is not negligible. This is in contrast to conventional analysis that bound error between target and generated distributions when the error in the score function is small.\nThe analysis is only employed in exploring the diffusion coefficient in reverse generative process. However, the analysis has potential implications for various applications where the correction or modification of score function is important. For instance, more theoretical treatments could be possible on image editing(classifier-free-guidance), fair generation(discriminator-guidance), or so on.", "weaknesses": "In Line 178, relation between the assumption and low-dimensionality of data distribution require further explanation to enhance understanding.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for taking the time to provide valuable feedback on our paper.\n\n### Q1: Need to clarify low-dimensional data distribution.\n\nWe acknowledge that the previous description lacks informative detail, and we aim to provide further elucidation as follows: \nGiven that many realistic datasets exhibit effective compact support, it becomes possible to identify a constant denoted as $\\widetilde{C}$, such that the distribution follows the inequality:\n$$\n\\rho_0(x) \\le \\widetilde{C} \\exp(-|x|^2/2),\\qquad \\forall x.\n$$\nConsequently, we can deduce:\n$$\ne^{-U_0(x)} \\le \\widetilde{C} \\exp(-|x|^2/2)\\ \\implies U_0(x) \\ge \\frac{|x|^2}{2} - \\log{\\widetilde{C}}.\n$$\nThis assertion corresponds to the third condition in Assumption 3.1.\n\nTo delve deeper, in situations where $\\rho_0$ is predominantly supported within $B_0(R)$, it is feasible to ascertain an adequately large $\\widetilde{C}$ such that $\\rho_0(x) \\le \\widetilde{C} e^{-|x|^2/2}$ within $B_0(R)$. Outside of this domain, the aforementioned assumption suggests that the decay surpasses that of a standard Gaussian distribution. In instances where data conforms to a low-dimensional subspace, a common practice is to pump in certain noise for the degenerate coordinate, e.g., approximate the distribution in the degenerate coordinate via $N(0, \\sigma^2)$ with small $\\sigma$ but NOT infinitesimally small. The above equation can hold true by selecting an appropriate $\\widetilde{C}$. We will incorporate these clarifications in our manuscript to enhance accuracy and clarity.\n\n### Q2: The practical implication of the score-training weight hyperparameter based on our analysis\n\nOur theoretical analysis indicates that for a large diffusion SDE, the sample quality is much less unaffected by score error occurring near the noise end, while the quality  is quite sensitive to score error near the data side ($p_0$). In contrast, for the ODE model,    the conclusion is reversed. This finding immediately suggests a practical hint for enhancing ODE flow about how to train the score function by adjusting the  weight  to favor ODEs:  one would like less error in the noise end and  the  noise-driven  weight (stated in the Appendix of manuscript, which use larger weights near the noise for score training) will likely enhance ODE flows.\n\nOur previous submission clearly demonstrated the effectiveness of this practical strategy  for MNIST (in Fig 5). The supplementary PDF features loss function distributions under various training weight functions, verifying our conjecture. We've also included initial CIFAR-10 results, subject to further exploration and incorporation in the revised version.\n\n### Q3: about limitations for larger datasets\n\nWe have conducted  preliminary assessments on CIFAR-10, and the findings are documented in the supplementary PDF. And our numerical outcomes on CIFAR-10 echo those observed in the MNIST dataset. We will undertake a more comprehensive exploration toward  a broader array of results to be incorporated in  the forthcoming revised manuscript."}, {"review_id": "PAVnH2Nmwl", "reviewer": "Reviewer_pYKm", "summary": "In this work the authors try to understand the impact of noise in the reverse process of diffusion models in the presence of an approximate score network. In particular, they look at how the Kullback-Leibler divergence between the true data distribution and the denoised distribution evolves w.r.t. the diffusion term $h$ of the backward (generative) process, with an asymptotic expansion of this KL divergence and studying the first non-zero term.\nThey then show that if the score has a pulse (delta) error early on in the generative sampling, higher values of $h$ are able to recover the density trajectory better, whilst with and error in the late stage of the denoising process it's the opposite.\nThey empirically back these results on synthetic experiments (With heuristic errors) and on the swiss-roll and MNIST dataset with trained score network.", "questions": "- Eq 5: missing the score? Or this is the process of the error?\n- Section 3.3: typo in title 'Placn'\n- Equation 9: This does not assuming $h \\gg 0$ right? Worth moving it above in Section 3.2 to avoid any ambiguity perhaps?\n- line 194: 'By convection [...]' what does this mean? Doesn't this simply come from the definition of the time-reversal?\n- line 196-200: So when $h \\rightarrow \\infty $ not only $V_t \\rightarrow U_t$ but because of the prefactor $h_t^2/2$ in Eq 9 (left) the distribution $\\mu_t$ converges extremely fast to to $\\mu_0$? Not sure to know what '“almost quasi-static” thermodynamics' refers to. Likely worth expanding a bit more on this in the main paper?\n- Can this be seen trough the perspective of Langevin dynamics as a corrector, which would 'project' $\\rho_t$ back to $p_t$ at $t>s$?\n- line 233: '$L(h)$ will decay to zero exponentially fast' with $h$?\n- Section 3.5: Here the setting in that the error accumulates over time in contract with the 'pulse' of Section 3.4. The results of Section 3.5 are that with error at the end of the generative process the smaller $h$ the better, yet if the pulse is also near the end, wouldn't large $h$ help (as the larger $h$ the faster it would correct this er§ror?).\n- MNIST: This section is not really self-contained as the weighting functions are in the appendix. What is the motivation for this cubic 'noise' weighting?\n- appendix: 'The computational cost of SDE-based models with large diffusion (h > g) is relatively high due to the necessity of a larger number of time steps.' What is the reason for this? That is quite key and worthy of being further developed (in the main text).", "limitations": "- Figure 4: It shows that indeed higher values of $h^2/g^2$ can lead to smaller error yet it also requires a lot more step. At $10^3$ steps, which is quite a reasonable number of steps already, the probability flow ODE ($h=0$) perhaps best and the curves starts overlapping around ~2000-3000 steps. So although the theoretical conclusion is that larger values of $h^2/g^2$ is better, for practical implementations there is a trade-off between larger number of discretisation steps (with larger values of $h$) and larger score network architecture.\n- The authors look at different types of errors on the score network, both theoretically with pulse-shape error to simplify the analysis, but also experimentally with simple errors (weighting the true score for some time range) which makes it easier to interpret. It would be interesting to look at the empirical error of trained score network (e.g. on image dataset), depending on the weighting $w(t)$ in the loss, but also the scheduling $g(t)$, and in particular its (average) evolution through the denoising time. Is the error spread uniformly (wrt time)? Writing the score as $\\epsilon / \\sigma(t)$ (with $\\sigma(t)$ the marginal variance at time $t$) couldn't one argue that with uniform error in the estimation of noise $epsilon$, the error in the score will blow up near $t=0$ thus the setting where the error is large close to the end of the generative process is the most likely setting (and therefore $h=0$ should work best)? Perhaps this is related with the fact that practitioners stop the denoising process at time $t=\\eta$ with $\\eta$ a small hyperparameter (e.g. \\eta=e-3)?\n- Would suggest increasing the fontsize of Figures (e.g. Fig 4)", "rating": 7, "confidence": 3, "soundness": 4, "presentation": 3, "contribution": 3, "strengths": "- I believe this paper tackle a really important question, which is the role of noise in the generative process. In particular with the success of flow matching which shows that continuous normalising flows, which are (deterministic) ODE based generative models, can be trained akin to denoising score matching models.\n- The analysis seems sounds and rigorous, and the assumptions reasonable.\n- The take home message (if I understood correctly) is that apart if the score error is concentrated at the end, the higher the level of noise the better (with infinite computational budget to actually discretise the backward process). \n- They link they theoretical findings with empirical results, both on simple settings which is nice to unambiguously be able to conclude, and more realistic settings when the score network is trained.\n", "weaknesses": "- No major weakness, but some of the interpretation of the theoretical results could perhaps be enhanced (cf following questions).\n- The fact that higher values of $h$ requires more discretisation steps is something that should be further stressed I believe.\n", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for taking the time to provide detailed and valuable feedback on our paper. The summary from the reviewer regarding the take-home message is accurate. For a detailed error distribution illustration in score-matching loss over time, see the supplementary PDF. We'll enhance presentation in the revised version based on suggestions. In the following, we address technical queries and concerns.\n\n### Q1: The fact that higher values of $h$ requires more discretisation steps is something that should be further stressed.\n\nIndeed, if one uses a large $h$ in the generative model, it holds crucial significance for time discretization. We will  give more explicit  statement about this practical limitation for large $h$ in the revised manuscript. We would like to highlight that this work does not say a large $h$ is favored in practice: it discussed both small and large $h$, but under different cases of score error.  Our analysis also indicates under what circumstance, the use of ODE model is preferred, which seems more efficient with less NFE. \n\n### Q2: Missing of score in Eq. 5\n\nThanks for pointing out this typo. In the second line below eq (5), ``$\\mathfrak{S}^{\\leftarrow}_t = \\epsilon\\mathcal{E}_t^{\\leftarrow}$'' should be $\\mathfrak{S}^{\\leftarrow}_t = \\nabla\\log p_t^{\\leftarrow} + \\epsilon\\mathcal{E}_t^{\\leftarrow}$. We will correct this.\n\n### Q3: About eq. (9)\n\nWe don't assume that $h^{\\leftarrow}_t \\gg 1$ (we only need $h^{\\leftarrow}_t > 0$) in Eq. (9). We will remove $h^{\\leftarrow}_t \\gg 1$ in the title of Sec. 3.3 to avoid confusion.\n\n### Q4: In the line 194\n\nWe mean \"by notational convention\". We'll rephrase it to \"by the notation of time-reversal\".\n\n### Q5: About Lines 196-200 and quasi-static dynamics\n\nIndeed, rapid convergence is unattainable unless the prefactor $\\frac{(h^{\\leftarrow}_t)^2}{2}\\to\\infty$.\nWe use \"quasi-static\" to indicate that the system will remain in proximity to the equilibrium $\\rho_t^{\\leftarrow}$ at time $t$, regardless of pulse perturbation. We will incorporate the following explanation in the revision: For any distribution ${\\mu}_t^{\\leftarrow}$ (even quite deviating from $\\rho_t^{\\leftarrow}$),  over a brief time $\\Delta t$ slightly exceeding $\\mathcal{O}(1/(h_t^{\\leftarrow})^{2})$, we can anticipate that\n\n $\\mu_{t+\\Delta t}^{\\leftarrow}$ is close to $\\rho_{t+\\Delta t}^{\\leftarrow},$\nassuming the evolution of $\\mu_s^{\\leftarrow}$ aligns with $\\mathcal{L}_{s}^{(h^{\\leftarrow})}$.\n\n### Q6: Relation to the Langevin corrector\n\nThis resembles a Langevin corrector: both our findings and the Langevin corrector are blessed by the same convergence property. However, nuanced distinctions exist: In literature, the Langevin corrector addresses errors from discretization schemes, without referring to the score estimation errors, whereas our work deals with diminishing effect from score error.\n\n### Q7: In line 233\n\nIndeed, it decays exponentially fast with respect to $\\mathsf{h}$.\n\n### Q8: Regarding Sec 3.5\n\nWhen the pulse occurs near (but not at) the generative process's end, larger $h$ leads to smaller sampling errors (as in Sec 3.4). For instance, an error like $\\mathbb{I}_{[T-2\\delta, T-\\delta]}(t) E(x)$, \nwhere $\\delta>0$ is fixed and $E:\\mathbb{R}^d\\to\\mathbb{R}^d$ still belongs to the case in Sec. 3.4. However, the findings of Sec 3.4 become inapplicable if $\\mathcal{E}_T^{\\leftarrow} \\neq 0$, addressing the subtleties in Sec 3.5. \n\nWe isolate the key aspect of $\\mathcal{E}^{\\leftarrow}$ as the following $\\mathbb{I}_{[T-a, T]}(t)E(x)$\nwith a small $a \\ll 1$. Larger $h$ may not be advantageous here and the tricky part is the complex scaling relationship between $a$ and $\\mathsf{h}$. For simplicity, let us say $a = \\frac{1}{\\mathsf{h}^2}$.\n\nExamining the score function error $\\mathcal{E}^\\leftarrow_t$ within $t\\in [T-a, T]$, the effective time of Langevin dynamics is only $\\mathcal{O}\\big(a \\frac{\\mathsf{h}^2}{2}\\big) = \\mathcal{O}(1)$, insufficient to return to equilibrium and reduce score estimate error. The $\\frac{\\mathsf{h}^2}{2}$ scaling comes from the prefactor in eq (9). Heuristically, the key distinct is that the SDE, with constant diffusion $\\mathsf{h}$, can address error $\\mathcal{E}_{t}^\\leftarrow$ in the interval $[0,T-\\frac{c}{\\mathsf{h}^2}]$ (section 3.4), but not in $[T-\\frac{c}{\\mathsf{h}^2}, T]$ (section 3.5), where $c$ is a constant irrelevant herein.\n\n### Q9: The cubic training weight\n\nPlease refer to our global discussion regarding the relationship between training weight, score error, dynamics selection, and numerical schemes. Our analysis suggests that one might enhance ODE via improving its score training near the noise end, and cubic weight function, increasing monotonically from 0 to 1 as $t\\in [0,T]$, reduces data-side loss contribution compared to the quadratic default in literature (i.e., this prioritizes noise-side effects in training loss). For more details, see the supplementary PDF.\n\n### Q10: Computational cost of large $h$\nThe rationale for introducing additional discretization steps stems from the amplified magnitude of both drift and diffusion terms (see Eq. (6)) when using larger $h$.  Enhancing the differential equation's magnitude is akin to prolonging dynamics, inevitably resulting in larger errors in most cases. This necessitates the choice of an improved scheme or a reduction in step size. \n\n### Q11: Distribution of loss over time and connection to the training weight\n\nSupplementary figures depict score-matching loss over time. We explored loss distribution with varying training weights $\\omega_t$. Investigating the impact of $g$ practically remains an interesting future endeavor. Errors are unevenly distributed across time, especially pronounced at the data side under conventional training loss. It appears that the current training weight scheme favors ODE over SDE. Indeed, such a score error distribution over time also connects to the adoption of the early-stopping techniques."}, {"review_id": "se5KSlg3W8", "reviewer": "Reviewer_Luhj", "summary": "The paper studies diffusion models, which comprise a class of generative models based on stochastic differential equations. In diffusion models, data is first transformed into Gaussian noise via a stochastic differential equation (usually an Ornstein-Uhlenbeck process), and then a backward process, which turns noise into synthetic data, is learned via the score function. In most works, the diffusion coefficient of the backward process is either chosen to be zero, or to be a constant matching the diffusion term of the forward. The first case is an ODE implementation, while the second case is an SDE corresponding to the exact time reversal of the forward process. However, any function of time can be chosen for the diffusion coefficient of the backward process, and in this work the authors explore the consequences of different choices of diffusion coefficient. In particular, they ask how the optimal choice of diffusion coefficient depends on the error induced by the estimated score function.\n\nThe main results of the paper are as follows. Let $h$ denote the diffusion coefficient of the backward process. Assume that the score function estimator has the for $s_t(x)= \\nabla \\log p_t(x) + \\epsilon E_t(x)$, where the first term is the exact score function and the function $E_t(x)$ is a bounded error function. Let $\\hat{q}_T$ denote the law of the backward process at time $T$ with $s_t(x)$ in the place of the true score function. Then under certain assumptions on $p_0$, it holds that $$ KL(p_0 \\parallel \\hat{q}_T = L(h) \\epsilon^2 + O(\\epsilon^3),$$ where the leading order term $L(h)$ depends on $h$, $p_0$, $E_t(x)$ and $T$. The leading order term $L(h)$ is computed explicitly, which allows the authors to analyze several specific cases when $h$ is constant. In particular, they show that if the error $E_t(x)$ is time-localized (i.e., only nonzero at $t=s$ for some $s > 0$) in the middle of the time interval, the leading order term $L(h)$ decays exponentially fast to zero as $h \\rightarrow \\infty$. This suggests that the SDE implementation can significantly outperform the ODE implementation in some cases. They also show that for certain forms of the error function, the opposite is true: the ODE implementation has better KL divergence error than the SDE implementation.\n\n\n\n--------------\nAfter rebuttal: Thank the authors for the clarifications and comments, which addressed my earlier concerns. I will keep my rating unchanged.", "questions": "1.\tCan you explain in what setting, if any, one would ever encounter a score estimator where the error term was localized in time? Or was the purpose of this choice of error term simply to illustrate that different choices of diffusion coefficient can lead to very different performance? It seems like a very strong assumption that would never be realized in practice.\n2.\tIt would be good to state in the main result any dependence on the sub-Gaussian constant $c_U$. As it is stated in the paper, one often works with a mollified version of the empirical distribution, which leads to a tradeoff between bias from smoothing and faster convergence for learning smoother distributions. Thus, rather than treat the constant $c_U$ as hidden, it could be beneficial to discuss how it balances with other terms.\n3.\tIn this paper, only constant choices of the diffusion term $h$ are discussed. Is this just for simplicity? Could it be of interest to consider certain parameterizations of nonconstant functions for $h$? If so, this could be mentioned for future work.", "limitations": "None", "rating": 7, "confidence": 4, "soundness": 4, "presentation": 2, "contribution": 2, "strengths": "Overall, the paper presents a novel and interesting question with some practical implications. The question of choosing the diffusion coefficient in the backward SDE seems largely unstudied because, as the authors mention, the ODE and reverse-time SDE approaches are standard conventions. The fact that they show that different choices of the diffusion coefficient can lead to vastly different performances suggests that the problem is worth studying further. \n\nThe authors provide a nice explicit characterization of the leading order term in the KL divergence error, the derivation of which requires a lot of work. Namely, they show that as $h \\rightarrow \\infty$, the term $L(h)$ converges exponentially fast to a constant $\\Tau$ which is upper bounded by $$ \\Tau \\lesssim \\frac{1}{2m_0^2} \\int_{\\mathbb{R}^d} \\frac{(\\nabla \\cdot (p_0 E_t))^2}{p_0} dx.$$ Here, $m_0$ is the second moment of $p_0$. This suggests the quantity on the right-hand side above plays a significant role in the error analysis of diffusion models, and thus it could be insightful to understand it for different classes of probability distributions.\n\nThe authors also leave a lot of room for interesting future work. Most interestingly, the authors pose the question of how to infer the optimal choice of $h$ from data, which could lead to a lot of subsequent research.", "weaknesses": "One weakness of the paper is that the upper bound on the constant $\\Tau$ is not further explained through discussion or examples. It seems like the upper bound plays a big role in bounding the KL divergence error for diffusion models. It would be very interesting to discuss explicit examples of distributions $p_0$ where the upper bound on $\\Tau$ can be explicitly computed or controlled, other than Gaussian distributions. This could be mentioned as a direction for future work. It would also be helpful to describe intuitively what the bound represents about a distribution (what does a large/small $\\Tau$ say about a distribution?).\n\nIn addition, the numerical experiment on the 2D 4-mode Gaussian mixture could benefit from additional explanation. It doesn’t seem like a machine learning problem, since you already know the score of the distribution from which you are trying to sample. Perhaps the point is simply to illustrate the convergence of $L(h)$ to a constant as $h \\rightarrow \\infty$, in which case the purpose of the experiment should be stated more clearly.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your appreciation and valuable feedback of our paper.  \n\n### Q1: More understanding about the upper bound.\n\nThis suggestion is insightful. We are also interested in exploring the structure of this upper bound. But we are afraid that a simple and heuristic  interpretation of this formula for upper bound is not available. Quantifying the upper bound's structure poses challenges due to the dependence of the score function error, which relies on multiple hyper-parameters and is problem dependent. We've observed some incomplete but interesting aspects regarding this question, and we'll discuss these findings in the revised manuscript.\n\nConsider employing constrained score models, where we parameterize the log probability $\\log p_t(x)$ during training, rather than the score function $\\nabla\\log p_t(x)$. The training process remains unchanged, differing only in how the score function is parameterized via a neural network. In this scenario, the error $\\mathcal{E}^\\leftarrow_{T} = \\nabla\\varphi$ for a scalar-valued function $\\varphi$ (a result derived from the nature of energy-based models), leading to the following modified upper bound expression:\n$$\n\\frac{1}{2m_0^2} \\int_{\\mathbb{R}^d} (\\Delta \\varphi - \\nabla U_0 \\cdot \\nabla\\varphi)^2 e^{-U_0} = \\frac{1}{2m_0^2} \\int_{\\mathbb{R}^d} (\\mathcal{L}^{*}\\varphi)^2 e^{-U_0},\\qquad p_0 := e^{-U_0}\n$$\n\nwhere $\\mathcal{L}^{*}(\\varphi) := \\Delta \\varphi - \\nabla U_0 \\cdot \\nabla\\varphi$. We can readily verify that its adjoint operator $\\mathcal{L}(\\mu) = \\nabla\\cdot(\\nabla U_0 \\mu) + \\Delta \\mu$ corresponds to the Fokker-Planck generator of the following Langevin dynamics $d X_t = -\\nabla U_0(X_t)\\ d t + \\sqrt{2}\\ d W_t$ where $W_t$ is the Brownian motion.\n\nGiven that the error $\\mathcal{E}^\\leftarrow_{T} = \\nabla \\varphi$ originates from approximating $\\nabla\\log p_0$, it is logical to consider the scenario where $\\varphi = \\log p_0$. In this context, the aforementioned upper bound can be expressed as follows:\n$$\\mathcal{T} \\lesssim \\frac{1}{2 m_0^2} \\int_{\\mathbb{R}^d} \\big(|\\nabla U_0|^2 - \\Delta U_0\\big)^2 e^{-U_0}.$$\n\nThe above two expressions appear to be slightly more informative than the formula in manuscript.\n\n### Q2: Gaussian mixture experiment\n\nThank you for the suggestion. Indeed, the primary use of GMM is for verifying the theory. We will provide clearer clarification of its purpose in the revised version.\n\n### Q3: About the practicality of pulse shape error.\n\nThank you for your question. A true pulse shape error might not manifest in real-world scenarios; it is, in essence, a theoretical construct. Nevertheless, we find value in studying it, which also hold practical significance:\n- The time-localized ansatz is a foundational element of our broader findings. It serves as a crucial component in establishing the general result presented in Prop. 3.6, which addresses a generic score error function.\n- The pulse-shape error emerges from analyzing the impact of score estimation errors at each time on the final sample generation error, rather than examining an averaged score-matching loss over time. This concept contributes to a deeper comprehension of the significance of both the overall score-matching loss and the distribution of score function errors across time.\n- To briefly elucidate the origin of the pulse-shape error, we would like to recall Prop. 3.2, where the leading-order term $L$ has a quadratic dependence on $v_T$, which in turn has a linear dependence on $\\mathcal{E}^{\\leftarrow}_t$. Since $\\mathcal{E}^{\\leftarrow}_t$ can be decomposed as a linear combination of step functions without loss of generality, it is effectively equivalent to investigating the influence of individual values of $\\mathcal{E}^{\\leftarrow}_t$ for a fixed $t$ on the final leading-order term $L(h)$—essentially examining a pulse-shape error. We will provide full technical details of this argument in revision.\n\n### Q4: Scaling of error with respect to $c_U$.\n\nBelow are our understandings. We shall discuss three instances where the role of $c_U$ becomes apparent:\n- In Prop. 3.4, $c_U$ explicitly appears in the main result. It influences the lower bound of $\\mathsf{h}$ required to ensure exponential decay. \n- In Propositions 3.5 and 3.6, where we address the asymptotic case, the scaling involving $c_U$ becomes harder to track. A lot more effort is required to develop non-asymptotic results and we will defer this to the next stage of research.\n- Another potential context in which $c_U$ could play a role is in determining the upper bound $\\mathcal{T}$ as mentioned by the reviewer. Characterizing how $\\mathcal{T}$ relates to $p_0$ (and consequently $c_U$) could offer valuable insights. However, one challenge lies in understanding the specific form of the score function error (which is problem-dependent). We acknowledge this intriguing question and plan to explore it further in future research endeavors.\n\n### Q5: discussion on the constant $h$\n\nThe primary reason for the time-independent choice is its simplicity. Furthermore,  the complexity increases significantly for non-asymptotic analysis; when we resort to  asymptotic cases ($h_t = 0$ and $h_t=\\infty$), time-dependence becomes less significant. In the case of large diffusion ($h_t \\gg 1$), little distinction arises from using time-independent $h_t = \\mathsf{h}$. Practicality also guides our choice, as a time-independent $h_t$ is easier to tune: only a single scalar parameter needs adjustment for potential sample generation error reduction. While acknowledging the intrigue of time-dependent cases, we plan to address this in future work, which will be noted and discussed in the revision."}]

Author rebuttal:
We would like to thank all reviewers for comments and suggestions. We are glad by the positive reception of the proposed question, and we hold respect and appreciation for the critical and diverse perspectives provided.

### About practicality and motivation:
A concern or question is about the discretization error (or practicality). It's well-acknowledged that the interplay between score function estimation error and discretization error collectively shapes the ultimate quality of sample generation.
Our overarching perspective encompasses four sequential facets:
1. We choose hyper-parameters for training.
2. Training scheme introduces score estimation error.
3. Score error guides dynamics selection.
4. Subsequently, a numerical scheme needs to be determined.

Rapid advances in inference solvers and diverse algorithmic forms pose many challenges in directly integrating discretization error into optimal dynamics study. Our current focus centers on the interplay between score training error and noise-level within SDE (ODE or large diffusion), i.e., the relation between (2) and (3) above. This is necessary in order to comprehend the whole four facets. Our objective is not to favor specific dynamics but to reveal their relationship and underlying mechanism. This approach may offer a practical avenue to explore optimal dynamics selection under finite NFE and a fixed numerical scheme, which awaits future research pursuits to integrate more factors.

### Explanation of supplementary PDF:
Regarding the supplementary PDF, we highlight the following: our theoretical analysis indicates that large diffusion SDE is unaffected by score error near the noise end (start of inference process) and is sensitive to score error near the data side ($p_0$), with the conclusion reversed for ODEs. This suggests the potential for enhancing ODE flow through training weight adjustment: less error in the noise end (namely, $N(0,I_d)$) is  beneficial for ODEs. Our previous submission demonstrated this for MNIST. Additionally, the supplementary PDF features how score-matching loss function distributions with respect to time under various training weight functions, in line with the above theoretical conjecture. We've also included initial CIFAR-10 results, subject to further exploration and incorporation in the revised version.

Meta-review: This paper studies the effect of diffusion coefficient on the accumulation of score approximation error, and in turn the error of generated samples as measured in KL divergence. There has been many excellent theory papers in the past few years on diffusion generative models, but most of these focus on discretization error, while simply assuming that the score approximation error is bounded by a constant. On the other hand, the error of score-approximation is equally important, but much less understood. 

The authors establish a connection between the scale of Gaussian noise in the reverse process, and the final error due to score-approximation. Their analysis shows that SDE can outperform ODE in certain cases, such as when the score approximation is localized to the middle of the generative process. The theory is supported by experiments on some simple synthetic data, as well as on MNIST and CIFAR-10. 

There are limitations, for instance, much of the analysis is asymptotic, and  assumptions such as time-localized error are not realistic, and the experiments focuses on relatively simple datasets. Nonetheless, I think the analysis in this paper is an important contribution for understanding score error accumulation, both for theorists, as well as for practitioners in tuning parameters of the reverse process.

Finally, I note that I largely discounted the reviews by wiU7 and LAed; the former because the criticisms raised are largely orthogonal to the quality of the paper; the latter because many excellent papers have indeed focused on analyzing only the discretization error, or the approximation error. For instance, the best existing bounds on discretization error are completely unchanged whether we assume the score error is 0 or bounded by some constant.
Decision: Accept (poster)