Quasi-Taylor Samplers for Diffusion Generative Models based on Ideal Derivatives
Authors: Hideyuki Tachibana
Mocho Go
Muneyoshi Inahara
Yotaro Katayama
Yotaro Watanabe
Venue: Submitted to ICLR 2023 (2023)
Keywords: diffusion models
score based models
neural generative models

Abstract: Diffusion generative models have emerged as a new challenger to popular deep neural generative models such as GANs, but have the drawback that they often require a huge number of neural function evaluations (NFEs) during synthesis unless some sophisticated sampling strategies are employed. This paper proposes new efficient samplers based on the numerical schemes derived by the familiar Taylor expansion, which directly solves the ODE/SDE of interest. In general, it is not easy to compute the derivatives that are required in higher-order Taylor schemes, but in the case of diffusion models, this difficulty is alleviated by the trick that the authors call ``ideal derivative substitution,'' in which the higher-order derivatives are replaced by tractable ones. To derive ideal derivatives, the authors argue the ``single point approximation,'' in which the true score function is approximated by a conditional one, holds in many cases, and considered the derivatives of this approximation. Applying thus obtained new quasi-Taylor samplers to image generation tasks, the authors experimentally confirmed that the proposed samplers could synthesize plausible images in small number of NFEs, and that the performance was better or at the same level as DDIM and Runge-Kutta methods. The paper also argues the relevance of the proposed samplers to the existing ones mentioned above. 

Under review as a conference paper at ICLR 2023
QUASI-TAYLOR SAMPLERS FOR DIFFUSION GENERA- TIVE MODELS BASED ON IDEAL DERIVATIVES
Anonymous authors Paper under double-blind review
ABSTRACT
Diffusion generative models have emerged new challenger to popular deep neu- ral generative models such as GANs, but have the drawback that they often require a huge number of neural function evaluations (NFEs) during synthesis unless some sophisticated sampling strategies are employed. This paper proposes new efficient samplers based on the numerical schemes derived by the familiar Taylor expansion, which directly solves the ODF/SDE of interest. In general, is not easy to compute the derivatives that are required in higher-order Taylor schemes, but in the case of diffusion models, this difficulty is alleviated by the trick that the authors call "ideal derivative substitution," " in which the higher-order derivatives are replaced by tractable ones. To derive ideal derivatives, the authors argue the "single point approximation." in which the true score function is approximated by a conditional one, holds in many cases, and considered the derivatives of this approximation. Applying thus obtained new quasi-Taylor samplers to image generation tasks, the authors experimentally confirmed that the proposed samplers could synthesize plausible images in small number of NFEs, and that the performance was better or at the same level as DDIM and Runge- Kutta methods. The paper also argues the relevance of the proposed samplers to the existing ones mentioned above.
INTRODUCTION
Generative modeling based on deep neural networks is an important research subject for both fundamental and applied purposes, and has been a major trend in machine learning studies for several years, To date, various types of neural generative models have been studied including GANs (Goodfellow al.) 2014), VAEs (Kingma et al. 2021 Kingma & Welling 2019), normalizing flows (Rezende & Mohamed, 2015). and autoregressive models (van den Oord et al. 2016b a). In addition to these popular models, a class of novel generative models based on the idea of iteratively refinement using the diffusion process has been rapidly gaining attention recently as a challenger that rivals the classics above (Sohl-Dickstein et al. 2015 Song & Ermon 2019 Song et al., 20206]
Song & Ermon 2020; Ho et |l.]| |020 Dhariwal & Nichol, (021). The diffusion-
generative
models have recently been showing impressive results in many fields including image (Ho ct al.
2020
Vahdat al. 2021 Saharia et al.
2021 Ho et al.
2021 Sasaki et al. 2021), video (Ho et al.
2022
text-to-image NNichol et al. 2021
Ramesh et al. 2022). speech (Chen et al. 2020 [021 Kong
al. 2021 Popov al. [202] K Kameoka et al. 2020). symbolic music (Mittal al. 2021). natural
language (Hoogeboom et al. 2021] Austin et al. 2022]) chemoinformaties (Xu et al. 2022). etc.
However, while the diffusion models have good synthesis quality, it has been said that they have fatal drawback that they often require a very large number of iterations (refinement steps) during synthesis, thousand, In particular, the increase in refinement steps critically reduces
ranging from hundreds the synthesis speed, as each step involves at least one neural function evaluation (NFE). Therefore, it has been a common research question how to establish a systematic method to stably generate good data from diffusion models in a relatively small number of refinement steps, or NFEs in particular. From this motivation, there have already been some studies aiming at reducing the NFEs (See $ 2) Among these, Probability Flow ODE (PF-ODE) (Song et al. 2020b) enable efficient and deterministic sampling, and is gaining attention. This framework has the merit of deriving a simple ODE by a straightforward conceptual manipulation of diffusion process. However, the ODE is eventually solved by using a black- box Runge- Kutta solver in the original paper, which requires several NFEs per step and is clearly costly. Another PF-ODE solver includes DDIM (Song al. 2020a), and is also

---

Under review as a conference paper at ICLR 2023
commonly used. It is certainly efficient and can generate plausible images. However, it was not originally formulated as a PF-ODE solver, and the relationship between DDIM and PF-ODE is not straightforward.
From these motivations, we provide another sampler to solve the same ODE, which performs better than on par with DDIM. The derivation outline is simple and intuitive: (1) consider the Taylor expansion of the given system, and (2) replace the derivatives in the Taylor series with appropriate functions; that's all.
The contribution of this paper would be as follows: (1) We propose novel samplers for diffusion models based on Taylor expansion of PF-ODE. They outperformed, or were par with Runge- Kutta methods. (2) To derive our algorithms, we show that the derivatives of score function can be approximated by simple functions. We call this technique the ideal derivative substitution. (3) It has been known that the 1st order term of DDIM is same the Euler method for PF-ODE. This paper gives further explanation for higher order terms of DDIM: we show that the proposed Quasi-Taylor method and DDIM are identical at least up to 3rd order terms. (4) The same idea can be naturally extended to derive a stochastic solver for reverse-time SDE, which we call R-SDE this paper.
BACKGROUND AND RELATED WORK
Diffusion Process to draw new data from target density; Let us first briefly summarize the framework of the diffusion-based generative models. Following Song et al. (2020b), we describe the mechanisms using the language of continuous-time diffusion process for later convenience. Let us consider "particles" {xr} moving in d-dim space obeying the following Itô diffusion,
dx, = fxx,, t)dt + g (xxt, l)dB,, SDE:
(1)
where B, is the d-dim Brownian motion whose temporal increments obeys the standard Gaussian. The drift f(:--) -) is d-dim vector, and the diffusion coefficient g(::)) is scalar. The SDE describes the microscopic dynamics of each particle. On the other hand, the "population" of the particles obeying the above SDE, i.e. density function p(x,, x),8), (t > s). follows the following PDEs, which known as Kolmogorov's 's forward and backward equations (KFE and KBE); the former is also known as the Fokker-Planck equation (FPE), see
0,p(xe,1 / x) = Vx, f(xi, Op(xx, / xx,8) + Ax, g(x,,t) t)²
FPE:
p(xi, xs; (2)
KBE: - - 0xx((x 8) f(x,,,) Vx,P(x,, xx)8) + (xx, 8)²
Dx, p(x,; Xxy
where Ax = Vx Vx is Laplacian. (FPE also holds for p(x;, t); consider the expectation Ep(x,.s) (:l.) These PDEs enables us understand the macroscopic behavior of the particle ensemble. For example, if f(x,
V2D, where Uxx) certain potential and D constant, then may verify that the stationary solution of FPE is p(x) e -U((x)/ It means that we may draw sample x that follows the stationary density by evolving the SDE over time. This technique is often referred to the Langevin Monte Carlo method (Rossky et al. 1978 Roberts & Tweedie 1996).
(U()), g((,
Some of the diffusion generative models are based on this framework, e.g. (Song & Hrmon 2019 2020), in which the potential gradient (U(x) is approximated by neural network.
Another systematic approach is considering the reverse-time dynamics (Song et al. 2020b). An approach is based on KBE eq. (3) Roughly speaking, FPE gives information about the future from the initial density, while KBE gives information about what the past states were likely to be from the terminal density. Here, instead of using KBE directly, is useful to consider variant of it which is transformed into the form of FPE, because it has an associated SDE that enables the particle-wise backward sampling (Stratonovich) 1965
Anderson
1982. ; see also F.3.2
[(x.; 8)2
R-FPE:
-0,p(xe, 8 t) = Vx, f(x,, ()p(,,, $ Xx, t) + Ax.
p(x,, 8 X4, f)
(4)
2
R-SDE:
dx, = - f(x,, s) ds) + g(x,, s)dB..
(5)
Hereafter, let g(x,, 1) g(() for simplicity. Then the specific forms of drift and diffusion coefficients are written as follows,
R-SDE coeffs: f(x,, () = (xx: 1) = f(x1,t) =(1)2 Vx, log logp(x;, (), g(() -()) (6)
2

---

Under review as a conference paper at ICLR 2023
Starting from certain random variable Xy, then by evolving the R-SDE reverse in time, we may obtain a Xo which follows p(xo, xr, 7") (i.e. the solution of R-FPE eq. (4). Therefore, the initial density p(x,. 0) of the forward dynamies eq. (2) is the true density, then may utilize this mechanism as a generative model to draw a new sample Xo from it.
Another approach is based on FPE eq. ( By formally eliminating the diffusion term of the FPE for the forward process, we can derive another backward FPE (see also $ E..31). Being diffusion- free, the backward FPE yields a deterministic ODE, which is called the Probability Flow ODE
(PF-ODE) (Song et al. 2020b), and is an example of neural ODEs (Chen et al. 2018). The population
density obtained by evolving this system exactly the same as the above R-SDE.
PF-ODE coeffs: f(x;, () (x,,, () ((x,,)) 1 20(00 log log,(x;, (). 5(() 0.
(7)
Some extensions of this framework include as follows. Dockhorn et al. 2021)) introduced the
velocity variable considering the Hamiltonian dynamics. Another extension is the introduction of
conditioning parameter, and guidance techniques using it (Dhariwal & Nichol 2021| Ho & Salimans
2021 Choi al. I02I) to promote the dynamics to go to a specific class of images, which has achieved remarkable results in text-to-image tasks (Nichol et al. 2021
Ramesh et al. 2022).
Variance-Preserving Model (VP-SDE Model): The solution of unconditioned FPE is written as the convolution with the initial density p(xo. 0) and the fundamental solution, or the heal kernel, p(x,, f Xg. 0), which the solution of the conditional FPE under the assumption that the initial density was delta function, p(xo. 0) = 5(x0 -x,). Although it is still intractable to solve this problem in general, well-known exception is the (time-dependent)
-haetineUhleebeckk (OU) process where 2 3, and g(,, t) = / /,.. Bg 3(t) is non-negative continuous function. The
f f(x,,) t) specific form of diffusion coefficient 3r has some options: a simplest one would be the linear function, and another would be the cosine schedule proposed in (Nichol & Dhariwal (02I); see also $ D In any cases, it is the OU process, the heat kernel is simply written as follows,
1
p(x,,t t Xo, 0) (xx VI- x,00 ofI), where of = 1 - exp (-f0 (-(
(8)
Hereafter, we denote the noise variance by Vs == o7. (In some literature, the signal level at = /1 o of is used as a basic parameter instead of the variance.) This model is referred to as the variance-preserving (VP) model by [Song et al. (2020b). It has good properties such as the scale of data ||x |2 is almost homogeneous, which is advantageous in neural models. However, the variance exploding (VE) model (Song et al. 2020bb in which the norm increases is also practicable, and the theory can be developed in similar manner.
Training Objective: In diffusion-based generative models, one estimates the score function [p(x,, [ptx,,o)[P(x,, learning has been referred to as the score matching (Hyvärinen & Dayan 2005 Vincent 201 I). However, the exact evaluation of this training target clearly intractable because of the expectation
Vx, log t(),,t) = Vx, log Ep(x.,0)
0j] by neural network So (x, t). This sort of
0p []]. so it has been common to consider a Variational Bayesian surrogate loss; Ho & Salimans (2021) showed that the following loss function approximates the negative ELBO,
L [1| Vvx V log (xx,, Xo, 0) So 112) [E[
Dixo
So(x,, ()
= =[[| w So ( 4,00 + / /Pw,, ()
(10)
where the expectation in eq. 10 is taken w.r.t Xo ~ D. ~V I). and Uniform (0. T
Some variants of the score matching objectives are also studied. For example, CCen et al. (2020)
reported that the L, loss gave better results than the T-2 loss in speech synthesis. Also, Kingma et al.
(2021) argued that the weighted loss with SNR-based weights improves the performance.
It should be noted that the above loss function will actually be very close the ideal score matching loss function in practice, where the probability is not conditioned on Xo. i.e..
Liteal Ejll V xx, x; log p(xe, () (x( 1)
(11)
This is because there almost always exists a point xy on the data manifold such that Vx, log p( p(x;, f) se Vx, log p(x,,t Xp. 0) holds with very high accuracy in very high-dim cases, because of the well- known "log-sum-exp max" law. For more details, see $ 3.3 and
3

---

Under review as a conference paper at ICLR 2023
Sampling Schemes for R-SDE and PF-ODE: Thus obtained (x is expected to finely ap- proximate - v xx, x, log p(x,, t), and we may use it in eq. (5). One of the simplest numerical schemes for solving SDEs is the and many diffusion generative models are actually using it.
method (Maruyama 1955 Theorem. 1) follows,
Euler- Mader-Maruyama: - X (x,, t) vhg(t)w, where I0,,I
(12)
where / > 0 is the step size. The error of the Euler-Maruyama method is the order of O(vh) in general, though it actually (h) in our case; this is because 9t,
As a better solver for
the R-SDE, the Predictor- Corrector (PC)-based sampler was proposed in (Song et al. 20206). The PC sampler outperformed the Predictor-only strategy, but requires many NPEs in the correction process, so we will exclude it in our discussion. Another R-SDE solver is the one proposed by Matrrouuu-aatiinuu et al. (202I), whose NFE per refinement step is 2.
On the other hand, there are also deterministic samplers for PF-ODE eqs.
(5) as follows,
Euler: X(-h f x, - hf,, (x,, t)
(13)
Runge-Kutte Kutge-Kutta: Xt-h - Σm byk,, where k; = fx, Σラ= a.yky, hcg)
(14)
where {ay} J. : } are coefficients of the Runge-Kutta (RK) method (see E.5). The error of the Euler method is (h), and that of the RK method O(h)), ) < m in general (Press et al. 2007
S 16). Another deterministic sampler is DDIM (Song et al. 2020a Eq. (13)), and is also understood
as a PF-ODE solver (Salimans & Ho 2022). Its NFE per step only 1, and is capable of efficiently generate samples.
DDIM: X₃-h Gr
+ Oo " 01) So(xe, t).
(15)
In addition, as concurrent work as ours, Lu et al. (2022) proposed the DPM-solver, which is based on the Taylor expansion of PF-ODE. However, as the gradient is evaluated using several different points, the NFE per step is greater than 1 in general. Liu al. (2022) proposed sampler based on the linear multi-step method, in which the NFE/step is reduced to except initial 3 steps. Another PF-ODE solver is the DEIS (Zhang & Chen 2022) which is based on the exponential integrator with some non-trivial approximations such as the polynomial interpolation of score function.
Other techniques that aimed to make sampling faster include as follows. Song & Ermon 2020 proposed variety of techniques to accelerate the sampling. Watson et (022)) proposed a Da-based optimization method to tune noise schedules for faster sampling. [uumman & Luhman (202) and
Salimans & Ho (002)) proposed distilling the pretrained teacher model to student model that can predict teacher's several steps in single step, which is efficient during the sampling but extra training
for distillation is required. [aoo et al. (2022a b) derived some analytic expressions of reverse dynamics to enable faster sampling.
PROPOSED METHOD: QUASI-TAYLOR SAMPLERS
3.1 MOTIVATION: HIGHER-ORDEL STRAIGHTFORWARE SOLVERS FOR R-SDE AND PF-ODE
As mentioned above, DDIM already exists as an efficient solver for PF-ODE, but it can only be considered PF-ODE solver up first-order terms (Song et al.) 2020a Salimans & Ho 2022) and it would not be clear enough whether it can be considered a higher-order solver for PF-ODE
Some other techniques eLu al. 2022 Liu et al. 2022 Zhang & Chen 2022 were designed as higher-order PF-ODE solvers, though their derivations are rather sophisticated and less simple. Since
PF-ODE and R-SDE provide the basis for the diffusion generative models, it would be beneficial to
develop samplers that directly solve them through intuitive and straightforward arguments.
From these motivations, we propose simple but efficient sampler based on the Taylor expansion, a very basic technique that is familiar to many researchers and practitioners. In general, Taylor methods are not very popular as numerical schemes because they require higher-order derivatives, which are not always tractable. However, in diffusion models, the derivatives are easily and effectively evaluated, albeit approximately. The validity of this approximation requires some consideration (see $ A B). but once accepted, an efficient sampler can be derived simply by substituting this approximation formula into the Taylor series. This section describes the details of the idea, and derives solvers for both OF-ODE and R-SDE, Entire sampling procedures are summarized in F
4

---

Under review as a conference paper at ICLR 2023
3.2 TAYLOR SCHEME FOR ODE AND
SCHEME FOR SDE
Taylor Scheme for Deterministic Systems For simplicity, we consider the 1-dim case here, but can easily generalized multidimensionall cases. (See $ E.I.1) Given a ODE d a (), where the function a is sufficiently smooth, then we can consider the Taylor expansion of it, using a differential operator L,- (ā, + a(t, xe )0,)), We can write the Taylor expansion of the path 2, as follows. Ignoring o(h)) terms of the series, we obtain numerical scheme of order p.
h² 2!
h 3! L, +
Leth 21 ha( + L,)
(16)
Itō- To-Taylor Scheme for Stochastic Systems In stochastic systems, the Taylor expansion requires modifications because of the relation [[dB = dt. If obeys a stochastic system dz, a(æ,, t) b(æ,, t)d then the path is written in stochastic version of Taylor-like series, which is often called the Itô-TTyy Taylor expansion, a.k.a. Wagner-Plater expansion (Platen & Wagner [ et al.. 1994] $ 2.3. ..3.B)), (Säākkkk & Solin 2019 $ 8.2). The Itô- Taylor expansion is based on the following differential operators Lg, Gg. which are based on Itô's formula (It 1944).
1
Lg = 0 + ((xx,)),
+ tr,,) G₁ = b(x,
(17)
In (Kloeden & Platen 1992), a number of higher order numerical schemes for SDEs based on the Itô oo-Taylor expansion are presented. One of the simplest of them is as follows. See also |E.112]
Theorem (Kloeden & Platen 19922, $ 14.2): An Itô- Taylor scheme of weak order 3 = 2). Ler In obeys the above SDE, and let the differential operators Lg, Gy be given by eq. (77). Then, the following numerical scheme weakly converges with the order of 8 = 2 (see $ E.44 Furthermore, in special case where Gūb = O, the strong 1.5 convergence is also guaranteed Kloeden & Platen 1992 $ 10.4).
wi 2
h² Gyb + Iya + (weh - 2.) + &
Iith < It ha + w.b
(18)
2
where Vhwe, hy lyTT, are correlated Gaussian random variables, and wt, 2t are given by Wi = and 2₁ = 2 + 2/3 where u1, U2 ~ (0, 1) (i.id.).. The notations o, Lga, etc. are the abbreviations for a( (), (Lga) 1), etc.
3.3 SINGLE POINT APPROXIMATION OF SCORE FUNCTION
Before proceeding, let us introduce the single point approximation of score function that Xx, log,(x,,t) t) almost certainly has some point Xo on the data manifold such that the following approximation holds,
Xx, logglxx,siic Vx, log p(xx,t xo, 0)p(xo,, 0)dx, Xx, logp(xx. Xo, 0).
(19)
To date, this approximation has often been understood as a tractable variational surrogate. However, the error between the integral and the single point approximation is actually very small in practical scenarios. More specifically, the following facts can be shown under some assumptions.
1. The relative Le distance between Xx, log p(x,( t) and Vx, log,(x,, t Xg. 0) is bounded above by v(1 - vz)/uz vn for any point Xo on the "data manifold" practical scenarios.
2. When the noise level is low v 2s 0, and the data space sufficiently high-dimensional, the distant points far from do not contribute to the integral. If the data manifold is locally k-dim subspace of the entire d-dim data space, where < < d, then the relative Le distance bounded above by around 2 k/dd
Of course, the single point approximation is not always valid. In fact, the approximation tends to break down when the noise level v4 is around 0.9 (SNR = (1 - v/)/v/ is around 0.1). In this region, the single point approximation Conversely, however, it would be also said that the error is as small as this level even in the worst empirical cases. For more details this approximation, see
deviates from the true gradient by about 20% in some cases.
5

---

Under review as a conference paper at ICLR 2023
3.4 IDEAL DERIVATIVE SUBSTITUTION
In order to adopt the above Taylor schemes to our problem setting where the base SDE is eq. ( and fr, f, are given by eqs. (6). 7 we need consider the following differential operators. Note that the time evolves backward in time in our case, the temporal derivative should be -0,-
I - (f,((xx))) 1)- Vxx)+
-== 1==-0,-(A(x;)) -(f:(x,,)) Vxx) +
3e
G. C.-VBL (1 Vxi),
2
3r
3t (x,::) t), f,((x,,) 3
3e (x t).
where (x,:
X, +
x
(20)
2 2/42
v
It is not easy in general to evaluate expressions involving such many derivatives. Indeed, for example, [-(- has the derivatives of the learned score function, viz. 0,S₂ (x.: 1) and ( Vx (x, (). which are costly to evaluate exactly, whether in approaches based on finite differences (as in (Lu ct al. 2022), pack-propagation, or the JAX paradigm (Bradbury et al. 2018), because they even- tually require extra evaluation of a deeply nested function other than So So(x,t), t), and extra memory consumption. Fortunately, however, by using the trick which the authors call the "ideal derivative substituuiion", we may write all of the derivatives as simple combination of known values, only consisting of x,, So(X,, t), Dt. fe and derivatives of (3, and no extra computation is needed. Since the score function has single point approximation eq. (19) we may assume that the derivatives should ideally hold following equalities. For derivation, see $ B.1
Conjecture (Ideal Derivatives). Under assumptions in $ A i.e. the data space R2d is sufficiently high dimensional d > 1, the data manifold .M C. Rd is also sufficiently high dimensional but much smaller than the entire space (1 & dim M < d), M is bounded, M sufficiently smooth locally, and the variance parameter Vi is close to 0 or 1; - then it is likely that the following approximations hold, where a E R&d is an arbitrary vector. We call them the "ideal derivatives".
Be
So(Xe. t)
(a Vx,) So(x,t))
a
-0,So(x,,t) 2/P.
X;
(21)
VV
Vii
To confirm the accuracy of this approximation, we compared empirical and ideal derivatives using MNIST (LeCun al. 2010) and CIFARIO (Krizhevsky) 2009). As result, it was confirmed that the approximation of spatial derivative, i.e. (a ()))) is usually very accurate; the cosine similarity between the empirical and ideal derivatives is nearly always 0.99 (Figure To). On the other hand, for the time derivative 0g it was confirmed that it is quite accurate when the time parameter t (and the variance 14,) are small, but the error increases when the time parameter t (and the variance 12,) become larger (Figure 9) See B.2 for more details.
3.5 QUASI-TAYLOR AND QUASI-ITÓ-TANE TAYLOR SCHEMES WITH IDEAL DERIVATIVES
As we can see in $ B.2 the ideal derivative approximation is sometimes very accurate while sometimes not. In any case, however, the error in the ideal derivative only affects the second or higher order terms of Tayloor series, and it will not be the dominant error in the whole. As there is an overall correlation between the true and ideal derivatives, the advantages will outweigh the disadvantages on average, and we can regularly use this approximation on a speculative basis, even though there cxist some cases where the approximation is not accurate.
If we accept the ideal derivative approximation, we can formally compute the symbolic expressions for the derivatives f,)- Lg f.). Lg(g). G3 (-f,) and Gç(g) that appear in the Taylor and Itô- Taylor series by routine calculations, which can be easily automated by computer algebra systems such as SymPy (Meurer al., (017) shown in $ B.3 By substituting thus obtained symbolic expressions into the above Taylor series, we can derive Taylor schemes for both PF-ODE and R-SDE as follows.
Algorithm 1 (Quasi-Taylor Sampler with Ideal Derivatives for PF-ODE), Starting from a Gaussian noise XT ~ (N I), iterate the following refinement steps until xo is obtained.
X Ph,nx, X, + 1²,,, So (x,, t)/ /Vi, where
(22)
Pi.
=
1 +
B.h
2
+
22
4
(0?
2
B.)
+
13
4
03
12
3.
2
+
3,
3
+....
(23)
u =
2
T
2
82
+
T
22
1207
+3v.
9.
2vz
3
+
(24)

---

Under review as a conference paper at ICLR 2023
Using terms up to O(h²), the sampler will have 2nd-order convergence (henceforth referred as Taylor 2nd), and using terms up to O(h"), the sampler will 3rd-order convergent (similarly, Taylor 3rd). If we up the (h) terms, the algorithm is same the Euler method.
Algorithm 2 (Quasi-Ite I-----t---yya Taylor Sampler with Ideal Derivatives for R- R-SDE), Starting from a Gaussian noise XT (O, iterate the following refinement steps until Xo is obtained.
X pi,1 + 14,n (x t)/ Dur + mi where
(25)
Pi =1
-/2 +
122 1
B1)
p -Beh
+
(26)
4
2
A VBvVhw h3/2
(-2/3,
2//3
9, (w: zt) 8/ 1
201
Z,
(27)
The Gaussian variables wt and z have dimension-wise correlations, and each dimension is sampled similarly to Theorem
At first glance, these algorithms may appear to be very complex. However, the computational complexity hardly increases compared to the Euler or Euler- Maruyama methods, because almost all of the computational cost accounted for by the neural network So(x,, (), and the costs for scalar values Pi,,: Pī,, and noise generation ni are almost negligible. It should also be noted that these scalar values can be pre-computed and stored in the memory before synthesis. Thus the computational complexity of these methods are practically equal to Euler, Muler-Maruyama, and DDIM methods.
Computation Cost:
Error from the Exact Solution of PF-ODE: The numerical error of the Quasi-Taylor method
from the exact solution increases depending on the following factors: (1) The truncation error of the Taylor series in each step, i.e. 0(h++1,, (2) The number of the steps i.e. O(1/h), (3) The training and generalization eтror of the score function, i.e. = C, and (4) The average eтror between the true and ideal derivatives of the score function ||8]]. the factors 3 and 4 could be zero, then the numerical error is the order of O(h)). Otherwise, the expected numerical error is roughly evaluated as follows,
error = 0 (h (hC + [[]1) + 12(C (|=l) +.... l)t
=0 ([C+h(C+)) [[]]) + 1 (C + [8]]) + +... h2)
(28)
That is, the error of Euler method is O(C-1) ft), the Heun method (2nd order Kunge-Kutta will be OLC+hC+ +2²), and the Taylor-2nd method is O(C+/((C+ predominant (h) term will not disappear. Therefore, the overall order of the error will not decrease even if we increase the order of Taylor series greater than p 2 3. Nevertheless, beyond such an order evaluation, specific coefficients in higher order terms can still affect the performance, which should be validated empirically.
As long as L, ||ō]| > the
IMAGE SYNTHESIS EXPERIMENT
Experimental Configuration: In this section, conduct experiments verify the effectiveness of the methods developed in this paper. Specifically, we compare the performance of the Fuler scheme eq. (13). Taylor 2nd & Taylor 3rd (Alg. [ DDIM (Song et al. (2000a), and the Runge Kutta methods (Heun and RK4 $ E.5 these are less efficient than others because of NFEs per step)
for PF-ODE, as well as the Euler-Maruyant scheme eq. (12) and Itô-Taylor (Alg. 2 for R-SDE. The datasets we used were CIFAR- 10 (32 x 32) Krizhevsky 2009) and CelebA (64 x 64) (Liu et al. 2015).
The network structure was not novel but was based on an existing open source implementation; we used the "NCSN++" implemented the official PyTorch code by Song et al. (2020b). The network consisted of 4 levels of resolution, with the feature dimension of each level being 128 - 128 - 256 + 256 + 256. Each level consisted of BigGAN-type ResBlocks, and the number of ResBlocks in each level was 8 (CIFAR-10) 10) and 4 (CelebA). The loss function we used was the unweighted L2 loss similarly to (Ho et al. 2020). The optimizer was Adam (Kingma & Ba 2014). The machine used for training was in house Linux server dedicated to medium-scale machine learning training with four GPUs (NVIDIA Tesla V100). The batch size was 256. The number of training steps was 0.1 M steps, and the training took about day for each dataset.
7

---

Under review as a conference paper at ICLR 2023
Euler
(b) DDIM
(c) Heun (NFE=2)
(d) RK4 (NFE=4)
(e) Taylor (2nd)
(i) Taylor (3rd)
(g) Euler-Maruyama
(h) Itô- Itô-Taylor
Figure I: Comparison of the synthesis results of CelebA (64 x 64) data. The number refinement
steps is N = 12.
The noising schedule was also the same as the existing one, the default configuration of VP-SDE (Song) et al. 2020b): 3, = 0.1 + 19.9/ V - expl -0.16 9.95(²) eq. (76). The integration duration was 1, and the step size h was constant, i.e. h T/N where N is the number of refinement steps.
As quality assessment metric, we used the Fréchet Inception Distance (FID) (Heusel et al. (017), To evaluate FIDs, we used the pretrained Inception v3 checkpoint (Szegedy et al. 2016), and resized all images 299 x 299 x 3 by bilinear interpolation before feeding them to the Inception network. For each condition, 10,000 images were randomly generated to compute the FID score. Note that in this experiment, the computational resources for training were limited, and training was stopped before it fully converged (only 0.1 M steps, while in some other papers the number of training steps was e.g. 1.3 M steps in (Song al. 2020b)). Therefore, it would be necessary observe relative comparisons between samplers rather than directly comparing these FID value to those presented in other papers.
Results: Figure and Figure 2 show random samples for each sampler. More examples are available in $ G The deterministic samplers considered in this paper generated plausible images much faster than the vanilla Euler-Maruyama sampler.
Figure 3a and Figure 3b reports the FID scores. From these figures, the following observations can be made. First, the proposed Quasi-Taylor methods have about the same or slightly better than DDIM. The reason for this is discussed in the next section $ } We also found that the Runge-Kute Kutta methods reduces FID in fewer steps overall. However, they also hit bottom faster. This may be due to the effect of the singularity at the time origin (see (D) in the final step. (This can be seen in Figure [6] In the second right column, the Runge-Kurt Kutta methods produce images similar to the other deterministic samplers, but the rightmost ones seem to be slightly noisier than the others).
Even though the ideal derivatives are only approximations and contain some errors, the convergence destinations of Quasi-Tayla methods were almost the same as the Kunge-Kutte methods. This suggests that the error the ideal derivatives is actually hardly problem, because in regions where the approximation error is large, the state x; is noisy to begin with (e.g. left 2/3 figures in Figure 16 and the approximation error is negligible compared to the noise that was originally there.
The proposed stochastic sampler (Itô- (to-Taylor) also showed sufficiently competitive results, in terms of both FID scores and visual impression. Comparison the figures in $ G Figure confirms that the Itylor method empirically reaches almost the same target Fuler-Maruyama method
8

---

Under review as a conference paper at ICLR 2023
(a) DDIM
(b) RK4 (NFE=4)
(c) Taylor (3rd)
(d) Itö-Taylor
Figure 2: Comparison of the synthesis results of CIFAR-10 (32 x 32) data. The number of refinement
steps is N 30.
Refinement steps
FID
Refinement steps FID
50
Filler
Taplor 2nd
Taplor 3rd
301
DDM
Fille Martana
no Tagle
Hin (ETT.)
Rid (Nibid
(NTTsc)
1211z2
90
160
1211s20
Relinement steps
Relisement
(a) CelebA (64 x 64)
(b) CIFAR-10 (32 x 32)
Figure 3: Comparison of FID scores (L).
much more accurately, and it could be expected to be safe alternative to Euler-Maruyama method when stochastic sampling is important.
DISCUSSION: RELATIONSHIP WITH DDIM
In the above experiment, the performance of the proposed Quasi-Taylor methods are found to be almost equivalent to that of DDIM. In fact, despite having distinctly different derivation logics, the proposed method and DDIM actually agree, at least up the 3rd order terms of h. Therefore, it is not surprising the results are similar; and the smaller h is, the closer the results are. This can be quickly verified by doing Taylor expansion of the coefficients of eq. (15) i.e., or
Or-h and (Øt-n or oz).
w.r.t. h. Although it is tedious to perform this calculation by hand, the computer algebra systems e.g. SymPy immediately calculate it. For this computation, see $C
This finding that truncating DDIM at the 2nd 3rd order of h yields exactly the same algorithms as the proposed Quasi-Taylor methods may be a useful insight for DDIM users, even if it does not lead them switch the regular sampler from DDIM Quasi-Taylor. That is, it offers an option of truncating the higher-order terms of DDIM.
CONCLUDING REMARKS
This paper proposed Taylor-expansion approach for diffusion generative models, particularly the Probability Flow ODE (PF-ODE) and the reverse-time SDE (R-SDE) solvers. The assumptions to derive our sampler were minimalistic, and the derivation process was straightforward. We just substituted the derivatives in the Taylor series by ideal ones.
The obtained Quasi-Taylor and Quasi-Itó-Taylon samplers performed better than par with DDIM and Runge- Kutge-Kutta methods. This fact implicitly supports the validity our approximations. Conversely, if we could find some examples where the
methods, DDIM and RK methods gave decisively different results, we might be able to gain a deeper understanding of the structure of data manifold and the fundamentals of diffusion models by investigating the causes of diserepancy.

---

Under review as a conference paper at ICLR 2023
Reproducibility Statement Pseudocodes of the proposed methods available in s and the derivation of the proposed method is described [.I BB33 The experiment based on open source code with minimal modifications to match the proposed method, and all the data used in this paper are publicly available. Experimental conditions are elaborated in $
Ethics Statement As a final note, negative aspects of generative models are generally pointed out, such the risk of reproducing bias and discrimination in training data and the risk of being misused for deep fakes. Since this method only provides a solution to existing generative models, it does not take special measures against these problems. Maximum ethical should be taken in the practical application of this method.
REFERENCES
Brian DO Anderson. Reverse-time diffusion equation models. Stochastic Processes and their Applications, 12(3):313-326, 1982.
Jacob Austin, Daniel Johnson, Jonathan Ho, Daniel Tarlow, and Rianne van den Berg. Structured denoising diffusion models in discrete state-spaces. Advances in Neural Information Processing Systems, 34, 2021.
Fan Bao, Chongxuan Li, Jiacheng Sun, Jun Zhu, and Bo Zhang. Estimating the optimal covariance with imperfect mean in diffusion probabilistic models. arXiv preprint arXiv:2206.07338 07309, 2022a.
Fan Bao, Chongxuan Li, Jun Zhu, and Bo Zhang. Analytic-DPM: an analytic estimate of the optimal reverse variance in diffusion probabilistic models. arXiv preprint vrXiv:2201.06500 2022b.
James Bradbury, Roy Frostig, Peter Hawkins, Matthew James Johnson, Chris Leary, Dougal Maclaurin, George Necula, Adam Paszke, Jake VanderPlas, Skye Wanderman-Milne and Qiao Zhang. JAX: composable transformations of Python+NumPy programs, 2018. URL
http::/ /ttp:///gtthbb... com. ccmm/qoogle/aaa.
Nanxin Chen, Yu Zhang. Heiga Zen, Ron Weiss, Mohammad Norouzi, and William Chan. Wave- Grad: Estimating gradients for waveform generation. In International Conference on Learning Representations (ICLR), 2020.
Nanxin Chen, Yu Zhang, Heiga Zen, Ron J Weiss, Mohammad Norouzi, Najim Dehak, and William Chan. WaveGrad 2: Iterative refinement for text-to-speech synthesis. In Proceedings of INTER- SPEECH, pp. 3765-3769. ISCA, 2021.
Ricky TQ Chen, Yulia Rubanova, Jesse Bettencourt, and David Duvenaud, Neural ordinary differential equations. arXiv preprint 1rXiv:1806.073666 2018.
Jooyoung Choi, Sungwon Kim, Yonghyun Jeong, Youngjune Gwon, and Sungroh Yoon. ILVR: Con- ditioning method for denoising diffusion probabilistic models. arXiv preprint 200.v:2008.022 2021.
David Roxbee Cox and Hilton David Miller. The theory of stochastic processes. Routledge, 1965.
Prafulla Dhariwal and Alexander Quinn Nichol. Diffusion models beat GANs on image synthesis. In Advances in Neural Information Processing Systems, 2021.
Tim Dockhorn, Arash Vahdat, and Karsten Kreis. Score-based generative modeling with critically- damped Langevin diffusion. arXiv preprint arXiv:2112.07068 2021.
Ian Goodfellow, Jean Pouget-Abadie, Mehdi Mirza, Bing Xu, David Warde-Farley, Sherjil Ozair, Aaron Courville, and Yoshua Bengio. Generative adversarial nets. Advances in Neural Information Processing Systems, 27, 2014.
Martin Heusel, Hubert Ramsauer, Thomas Unterthiner, Bernhard Nessler, and Sepp Hochreiter. GANs trained by a two simesccle update rule converge to local Nash equilibrium. Advances in Neural Information Processing Systems, 30, 2017. https://ggitii com wttp:/iwco..o.o.tttt. jut.tUUUk
10

---

Under review as a conference paper at ICLR 2023
Jonathan Ho and Tim Salimans. Classifier-free diffusion guidance. In NeurlPS 2021 Workshop Deep Generative Models and Downstream Applications, 2021.
Jonathan Ho, Ajay Jain, and Pieter Abbeel. Denoising diffusion probabilistic models. Advances in
Neural Information Processing Systems, 2020.
Jonathan Ho, Chitwan Saharia, William Chan, David Fleet, Mohammad Norouzi, and Tim Salimans. Cascaded diffusion models for high fidelity image generation. arXiv preprint arXiv:2 /5282, 2021.
Jonathan Ho, Tim Salimans, Alexey Gritsenko, William Chan, Mohammad Norouzi, and David Fleet. Video diffusion models. arXiv preprint arXiv:2204.03458, 2022.
Emiel Hoogeboom, Didrik Nielsen, Priyank Jaini, Patrick Forré, and Max Welling, Argmax flows and multinomial diffusion: Towards -uooreeoreeee language models. Advances in Neural Information Processing Systems, 2021.
Aapo Hyvärinen and Peter Dayan, Estimation of non normalzzed statistical models by score matching. Journal of Machine Learning Research, 6(4), 2005.
Kiyosi Itô. Stochastic integral. Proceedings of the Imperial Academy, 20(8):519-524, 1944. doi: 10.3792/pia/1 195572786
Alexia Dolicoeur-Martinee Ke Li, Rémi Gotta go fast when generating data with score-based models. arXiv preprint arXiv:2105. 140.10000 2021.
Tal Kachman, and Joannis Mitliagkas.
Hirokazu Kameoka, Takuhiro Kaneko, Kou Tanaka, Nobukatsu Hojo, and Shogo Seki. VoiceGrad: Non-parallel any-to-many voice conversion with annealed Langevin dynamics. arXiv preprint arXiv:2010.02977 2020.
Samuel Karlin and Howard E Taylor. A second course in stochastic processes. Academic Press, 1981.
Diederik P Kingma and Jimmy Ba. Adam: method for stochastic optimization. arXiv preprint arXiv: 1412.6980, 2014.
Diederik P Kingma and Max Welling. An introduction to variational autoencoders. arXiv preprint 1eXivv:906.000 12091, 2019.
Diederik P Kingma, Tim Salimans, Ben Poole, and Jonathan Ho. Variational diffusion models. arXiv preprint urXiv:2107.006333 2021.
Peter E Kloeden and Eckhard Platen. Numerical solution of stochastic differential equations. Stochastic Modelling and Applied Probability, Applications of Mathematics. Springer, 1992. doi:
12616-5. Corrected 3rd Printing.
Peter E Kloeden, Eckhard Platen, and Henri Schurz. Numerical solutiio of SDE through computer experiments. Universitext. Springer, 1994. doi: 10.1007/978.33-442-577933. Corrected 3rd Printing.
Zhifeng Kong, Wei Ping, Jiaji Huang, Kexin Zhao, and Bryan Catanzaro. Diffwave: A versatile diffusion model for audio synthesis. In International Conference on Learning Representations (ICLR), 2021.
Alex Krizhevsky. Learning multiple layers of features from tiny images. Technical report, 2009. Dataset license: MIT.
Beatrice Laurent and Pascal Massart, Adaptive estimation of quadratic functional by model selection. 2000.
Annals of Statistics, pp.
Yann LeCun, Corinna Cortes, and CJ Burges. MNIST handwritten digit database, 2010. http: //yann. lecun. com/ exdo / mn ist Dataset license: CC-SA 3.0.
Luping Liu, Yi Ren, Zhijie Lin, and Zhou Zhao. Pseudo numerical methods for diffusion models on manifolds. arXiv preprint arXiv:2202.007778 2022.

---

Under review as a conference paper at ICLR 2023
Ziwei Liu, Ping Luo, Xiaogang Wang, and Xiaoou Tang. Deep learning face attributes in the wild. In Proceedings of the IEEE international conference on computer vision (ICCV), pp. 3738,, 2015. Dataset license: available for non-commercial research purpose only.
Cheng Lu, Yuhao Zhou, Fan Bao, Jianfei Chen, Chongxuan Li, and Jun Zhu. DPM-solver: A fast ODE solver for diffusion probabilistic model sampling in around 10 steps. arXiv preprini arXiv:2206.00022,, 2022.
Eric Luhman and Troy Luhman. Knowledge distillation in iterative generative models for improved sampling speed. arXiv preprint (0t::200110238800 2021.
Gisirö Maruyama. Continuous Markov processes and stochastic equations. Rendiconti del Circolo Matematico di Palermo, 4(1):48-90, 1955.
Aaron Meurer, Christopher Smith, Mateusz Paprocki, Ondfej Čertík, Sergey B. Kirpichev, Matthew Rocklin, AMIT Kumar. Sergiu Ivanov, Jason K. Moore. Sartaj Singh, Thilina Rathnayake, Sean Vig, Brian E. Granger, Richard P. Muller, Francesco Bonazzi, Harsh Gupta, Shivam Vats, Fredrik Johansson, Fabian Pedregosa, Matthew J. Curry, Andy R. Terrel, Stěpán Roučka, Ashutosh Saboo, Isuru Fernando, Sumith Kulal, Robert Cimrman, and Anthony Scopatz. SymPy: symbolic computing in Python. PeerJ Computer Science, 3:3 January 2017. ISSN 2376-5992. doi: 10.7717/peerj-cs. 103. URL https:: org 7ttp://1107777711ee----- 103
Gautam Mittal, Jesse Engel, Curtis Hawthorne, and Ian Simon. Symbolic music generation with diffusion models. In Proceedings of the 22nd International Society for Music Information Retrieval Conference (ISMIR), 2021.
Alex Nichol and Prafulla Dhariwal. Improved denoising diffusion probabilistic models. In Interna- lional Conference on Machine Learning (ICMI). PMLR, 2021.
Alex Nichol, Prafulla Dhariwal, Aditya Ramesh, Pranav Shyam, Pamela Mishkin, Bob McGrew, Ilya Sutskever, and Mark Chen. Glide: Towards photorealistic image generation and editing with text-guided diffusion models. arXiv preprint 12 2021.
Bernt Øksendal. Stochastic differential equations: an introduction with applications. Springer, 2013.
Eckhard Platen and Wolfgang Wagner. On Taylor formula for class of Ito processes. Probability and Mathematical Statistics, 3:37-51, 1982.
Vadim Popov, Ivan Vovk, Vladimir Gogoryan, Tasnima Sadekova, and Mikhail Kudinov. Grad- TTS: A diffusion probabilistic model for text-to-speech. In International Conference on Machine Learning (ICML), pp. $599-86088 PMLR, 2021.
William H Press, Saul A Teukolsky, William Vetterling, and Brian P Flannery. Numerical recipes 3rd edition: The art of scientific computing. Cambridge university press, 2007.
Aditya Ramesh, Prafulla Dhariwal, Alex Nichol, Casey Chu, and Mark Chen. Hierarchical text- conditional image generation with CLIP latents. arXiv preprint arXiv:2204.06125, 2022.
Danilo Rezende and Shakir Mohamed. Variational inference with normalizing flows In International conference on machine learning, pp. 1530-1538. PMLR, 2015.
Gareth O Roberts and Richard L Tweedie. Exponential convergence of Langevin distributions and their discrete approximations. Bernoulli, pp. 341-36313 1996.
Peter Rossky, Jimmie D Doll, and Harold L Friedman. Brownian dynamics as smart Monte Carlo simulation. The Journal of Chemical Physics, 69(10):4628 4633, 1978.
Chitwan Saharia, Jonathan Ho, William Chan, Tim Salimans, David J Fleet, and Mohammad Norouzi. Image super-resolution via iterative refinement. arXiv preprint 2rXiv:2104.0700 2021.
Tim Salimans and Jonathan Ho. Progressive distillation for fast sampling of diffusion models. arXiv preprint w:Xiv:2202.005/// 2022.
12

---

Under review as a conference paper at ICLR 2023
Simo Särkkä and Amo Solin. Applied stochastic differential equations, volume 10 of Institute of Mathematical Statistics Textbooks. Cambridge University Press, 2019.
Hiroshi Sasaki, Chris G Willcocks, and Toby P Breckon. UNIT-DDPM: Unpaired image translation
with denoising diffusion probabilistic models. arXiv preprini arXiv:2104.05358, 2021.
Steven E. Shreve. Stochastic calculus for finance Continuous-tine models, volume 11. Springer, 2004.
Jascha Sohl-Dickstein, Eric Weiss, Niru Maheswaranathan, and Surya Ganguli. Deep unsupervised learning using nonequilibrium thermodynamics. In International Conference on Machine Learning (ICML), pp.
PMLR, 2015.
Jiaming Song, Chenlin Meng, and Stefano Ermon. Denoising diffusion implicit models. In Interna- tional Conference on Learning Representations (ICLR), 2020a.
Yang Song and Stefano Ermon. Generative modeling by estimating gradients of the data distribution. Advances in Neural Information Processing Systems, 2019.
Yang Song and Stefano Ermon. Improved techniques for training score-based generative models.
Advances in Neural Information Processing Systems, 2020.
Yang Song, Jascha Sohl-Dickstein, Diederik P Kingma, Abhishek Kumar, Stefano Ermon, and
Ben Poole. Score-based generative modeling through stochastic differential equations. In Inter- national Conference on Learning Representations (ICLR), 2020b. https://github.ccm com/ yan.--sona song/sccoresdd__ptceeee pytorch (Apache-2.0 license).
Ruslan Leont'evich Stratonovich. Conditional Markov processes. In Non-linear transformations of stochastic processes, pp. 427-453, Elsevier, 1965.
Christian Szegedy, Vincent Vanhoucke, Sergey loffe, Jon Shlens, and Zbigniew Wojna. Rethinking the inception architecture for computer vision. In Proceedings of the IEEE conference on computer vision and pattern recognition, pp. 2818-2826, 2016.
Arash Vahdat, Karsten Kreis, and Jan Kautz. Score-based generative modeling in latent space. Advances in Neural Information Processing Systems, 2021.
Aaron van den Oord, Sander Dieleman, Heiga Zen, Karen Simonyan, Oriol Vinyals, Alex Graves, Nal Kalchbrenner, Andrew Senior, and Koray Kavukcuoglu. WaveNet: A generative model for raw audio. arXiv preprint arXiv: 06096666999999 201
Aaron van den Oord, Nal Kalchbrenner, and Koray Kavukcuoglu. Pixel recurrent neural networks. In International conference machine learning, pp. 1747-1756. PMLR, 2016b.
Pascal Vincent. A connection between score matching and denoising autoencoders. Neural computa- fion, 23(7):1661-1674 2011.
Daniel Watson, Jonathan Ho, Mohammad Norouzi, and William Chan. Learning to efficiently sample from diffusion probabilistic models. arXiv preprint arXiv:2106.0388 038000 2021.
Minkai Xu, Lantao Yu, Yang Song. Chence Shi, Stefano Ermon, and Jian Tang. Geodiff: A geometric diffusion model for molecular conformation generation. arXiv preprint (XXix:203300222. 2022.
Qinsheng Zhang and Yongxin Chen. Fast sampling of diffusion models with exponential integrator. arXiv preprint arXiv:2204.13902 2022.
13

---

Under review as a conference paper at ICLR 2023
A ON THE APPROXIMATION Xx, (x,, t) ~ V x; log ( xo, 0)
Although Vx, log t) and Xx, log p(x,, t Xp. 0) are clearly distinct concepts, they are nearly equivalent in practical situations that the data space is high-dimensional and the data are distributed in a small subset Now-dimensional manifold) of the space. That is, such a case, the integrated gradient Xx, log p(x,, t), given by
Vx, log p(x,, t) - Xx, log / p(xe, & xo, 0)p(xo: 0)dxo
xx, / p(x,, No. ()p(Xu, 0)dxa
J /p(x,, 1 No. (xp 0)dxo
Ep(xo) Vx. Pxx, Xo. 0)]
Ep(xx) (p( (xx;1 xo, 0)]
Ep(xo)
X x--V--ix p(x,, xo. 0)]
E p((x)) [p(x,,t Nu, 0)]
X vI 1/:00
|x - VI- mixol /1? yu
Ep(xa)
(29)
14
Ep(xx) (e- 11x [xx--II vixll 12
=Vx, log f
almost always has single point approximation which is written as Vx, log p(x,, t x 0), where
(6)
(2)
X is certain point the data manifold. There are phases depending the noise scale
ag = Vi, and both phases have different reasons for the validity the approximation.
A.1 PHASE (1): ANYONE CAN BE A REPRESENTATIVE (ā, 1/4 ())~0) 0)
If x; is far from any of { (VI- ( xx p(xo)): ). the gradients from each scaled data points almost the same. Therefore, the integrated one viz. V log p(x,, t) can approximated by any of intuitively illustrates the reason.
V log p(x,, (i) 0) for any xx") p(xu). Figure 4
More quantitatively, noting that the weight function q(-) satisfies Ep(xo) (x(x xx)] = 1, the L2
distance between log pxx,, t) and log, xx 0) bounded above as follows,
V log pxx+, t) V log p(xe,
xx0
(i) 0)
(4) 2
Ex(xo)
X; VI 14X0
q(xo xx)
X1
2
Xt VI V1X0
X: /I 1220
(i)
2
Ep(xo)
q(xo
xx)
V4
1/2
Bp(xoo)
(x50 - xx0 (x)) xxi)
u
V/
Ex(xx)
||x00 -x0 1129(x0 xx)
(t. Jensen's inequality)
V4
max Xo
||xu xx (1)
(p la(xo xx)]
12
14
|XX0- x6" 11
max Xo
4(1
V/ )
R2
(30)
where RM is the radius of the smallest ball that covers the data manifold M. The radius RM will be finite constant in most practical scenarios we are interested in. For example, if the data space is d-dim box, i.e., M C [0. 1]d, then R. is bounded above by RM vd/2.
14

---

Under review as a conference paper at ICLR 2023
gradient log p(x,, t) V Ploop(xe xxf
scaled datapoints (0) {vI 1/X9 [i=1
Figure 4: Intuitive reasons why a single point approximation is valid. The case where - RJ 0
Noting that | log p(x,, (0) 0)ll2 vd/u (this is actually random variable that follows the x-distribution D.
the relative L2 error is evaluated as follows,
V VVogp(x,, V log p(x, t X (
,0)l2
1 u
(relative L2 error)
(31)
Vlogp(x,,, x 0)ll2
Vt
Similarly, the cosine similarity is evaluated as
cossim V log p(x,, t), V log 1x 00) 2 1-
1 (1
1+
1
1 1 U
(32)
Vt
1/2
These bounds suggest that the relative error between the true gradient and the single point approxi- mation approaches to 0, and the cosine similarity goes to 1, as + 1 (āz - 0). That is, the single point approximation largely valid when the noise level 24 is high.
Proofs for some facts used in this section
Jensen's 's inequality: A sketchy proof of Jensen's 's inequality for Lz-norm is as follows:
2
2
Σ wn xm)
Σ
Σ
Wn (m)
Σ
Wrnt (11)
(n)
Σ Wn. (=))))== = Σ wn |x (33)
Wn (x)
Cosine similarity: Given the relative distance between x and as ||x - y |||l/|||||2 = r (for
simplicity, let < 1) then the cosine similarity between the vectors evaluated follows,
(x,y)
||x||2
(x,y)
||x|||
||x|
+ |yl1 - ||x y1
cossim(x.y)
||||||||||
|yy|a
|x|
2
l|yll2
|
||x|| [1
|
|Xx||
|ylla
|
= 2 | +
(34)
2 ||y
|
|
la
Noting that the triangle inequality implies that the Le norm of y bounded as follows,
[--r|||/|2 llyll2 (1 x|||2
(35)
we can put |y 12yl2/110 l2/l///////22///// = 1 (0 < |0| < < 1), and can evaluate the cosine similarity as follows,
1
cossim(x,y)
1 + +(1+ 8)
1+8"
- (1 ((1 8) + (1+ - (1 (1- > - 2 ~
(36)
'Given Gaussian random variables X, (O0 1), then the squared sum of them ΣL X² follows the x" distriibutio with d degrees of freedom. is well known that x² distribution converges to (d,, 2d) as d increases. Therefore, If is sufficiently large, ||x||2 = (I x/2 v Vd /2d R vd + v1/2. Thus the Le norm of d-dim Gaussian variable is approximated by vd, and the enror is as small as the scale of / 1/2. implies Gaussian distribution looks like sphere in high-dimensional space, contrary to our low-dimensional intuition. Indeed, the following incquality holds (Laurent Massam 2000 Lem. 1):
p((1 /2 ||x|/// < (1 2///d d 2y/d) /22) 2 1 2e
15

---

Under review as a conference paper at ICLR 2023
(a) Discrete case
(h) Continuous case
Figure 5: Intuitive reasons why single point approximation is valid. The case where - 14, R3
A.2 PIASE (2): WINNER TAKES ALL (ā, = 1 > 0)
The above bounds suggest that the single point approximation is valid only when the noise level v, is high (vt s 1). However, we can also show that the approximation is also valid because of another reason when the noise level V, is low (Vz 0).
If p(xo) is discrete distribution, and if i s O, the weight factor g(x) certainly "one-hot vector" such that
Xx) in eq. (9) is almost
q(x)) (6) xx) 1, and q(x)
X
= 0 i).
(37)
That is, the entropy of (xxo Xe) is nearly zero (see Figure
8 It is almost certain (prob R 1), since
(i) the converse is quite rare; i.e., there exists another data point x near X, p(x,, ()) 0) in such situation Figure 5a intuitively shows the reason.
If p(xo) is a continuous distribution, very small area (almost a single point) of the data manifold contributes the integration eq. (29) Figure Seguur/Sbjintuuue intiituity shows the situation. The weight function will become nearly the delta function, and the integration will almost the same as the single point approximation,
(xp xx) x 6(x, - x
(38)
where xij is a certain point close to the perpendicular projection of x; on the data manifold. These are related to the well-known formulae e log Σ exp(-/s) and log f exp(-/e) goes to max(-) when e + O. Naturally, this more likely to be true for high-dimensional space and is expected to break down for liw-dimensional toy data.
Now lct us elaborate the above discussion. As the scale we are interested now is very small, the data manifold is approximated as flat k-dim space. In addition, as q(: xx) is Gaussian function, decays rapidly with distance from the center (perpendicular projection of x2 E Rd to the k-dim
Let us evaluate the probability more qualitatively using toy model. Let denote k-dim ball of radius r by B² (r), and let D; be the L2 distance between x, and /1 vxxu
(3) i.e. D; := -=|| - /1 xxx, (1) ||a P
v /p.d. The question is whether there exists another data point ()) [j # v) in the discrete data distribution D = {xb }s such that Dj < D.. In other words, whether there exists another data point (Vvvv) centered at 111111/ [1-vixje vxxa (6) If it is the case, q(x))) (1) x.) has significantly large value compared to q(x.)" x.). and it will have strong effects on the result of integration eq. (29). Now let us show that it is unlikely in high-dimensional data space. As it is difficult evaluating the probability for general cases, let us consider toy model that the all the points are uniformly distributed in k-dim ball (k & d). That is, we assume that the data manifold
M = (* which satisfies D C M C Rd. Then, the probability we are interested in is roughly evaluated
D2) x1 (1 (1-vol//- u. -
(vlkk) =1
-1-(1---
vollB* vollBB (vnk) vol[B* vn)k)]]
( D.)
vollV1
Remembering that volume of k-dim ball given by vol[B* (r)] mated follows.
h the probability is approxi- F(K/241)"
k/2
r.hs= (1-
Vek
k/2
x/2
if Ds
mi)l
Thus, when in is sufficiently small the data manifold sufficiently high-dimensional. it almost unlikely that there exists such point x(5)) (j of i) unless very large.
16

---

Under review as a conference paper at ICLR 2023
subspace). For this reason, only small ball around the perpendicular foot of x; contributes to the integration. As the majority of k-dim Gaussian variables of variance o'I are distributed on sphere of radius /ka, would be sufficient to consider a ball of radius (vl: + 2)
Now the integration eq.
29 is approximated follows,
eq.
Vx B
lxg,(xx, |1 xo. 0) xx)p(xo)] dxo.
where B = {xo |l II - 1:X0 x1112 (vk+2) (39)
As B is small, the gradient vector Xx, p(x,, Xu, 0) is almost constant in this region. Quanti- tatively, let x
(1) xx² B, then the distance between these points is evidently bounded above by ||x0 (2) l2 s 2(vT / which implies
V log (( x
200000 ,0) - V Ploop(xx ((x, x
2
4xx0 (1)
(2)
X
X,
1/
1/
2
V v(x1) (1)
X
-I-v
22V2 + 2)
Vt
11 1 1-V
2(V + 2)
(40)
11
Noting that // log Xo, 5)l22 dd/vr, the relative L2 error between two gradients is bounded above by
V log p(x,,
X (1)
0) - Voogp(x,, x .00l 2 + 2)V1 v4
2 d
(41)
V logp(xx Xt,
xx= 0)112
Vd
By integrating both sides w.r.t. XD i.e. 5()(2))) xx)x(xx))) (2) ddxe (1) inequality for any point x E B.
we obtain the following
V log,((;,1) 1) V 1xx, / x5. 0)ll2
& 2 d'
(42)
V log,(x,, t x 0)
Now, let recall that a noised data ~ p(x; t) is generated in the following procedure.
1. Draw a data x from p(xo, 0).
(3) 2. Draw Gaussian x2 ~ p(x,, x 0) WVV - v.La).
Here, the data point perpendicular foot of x, is about vurki which is clearly smaller than (Vh- 2T. Thus we can write the following inequality.
(
E E is in the ball B, because the distance between VT 4xx0 and the
(8)
logp(xx,() 1) V log 1(x,, x .))l2 2
2 d'
(43)
|V log & x and we may conclude that the following approximation largely valid if v Rs 0, k d.
0)ll2
v Vlogp(x,,t) t) V Vlogp(x,,t 1x( 0).
(44)
We used the following relation: [[x]f(xx) yll2 |[x(f(xx yl1a {xll((() y112)-
17

---

Under review as a conference paper at ICLR 2023
A.3 COMPARISON OF THE EMPIRICAL SCORE FUNCTION AND THE SINGLE POINT APPROXIMATION
Let us empirically validate the accuracy of single point approximation using real data follows,
. D (MNIST (LeCun et al.
2010) 60,000 samples),
. D [CIFAR-IL 10 (Krizhevsky 2009) 50,000 samples).
Since the true score function cannot be determined without knowing the true density (which will be possible with synthetic data, but discussing such data will not be very interesting here), the empirical score function was calculated using the real data 'D above as follows,
True Score
log p(g, pxx(x) (x x x V log,(xx;1 xo, 0)]
xx) V log Xo, 0)] =: Empirical Score. [x(x
(45)
X
The evaluation of empirical score function using the entire dataset is unrealistic if the dataset D is large, it is feasible is small dataset like MNIST and CIFAR-10.
In order to evaluate the accuracy of single point approximation, we evaluated following three metrics.
Relative T.2 error between the empirical score function and V log & Xp, 0),
Cosine similarity between the empirical score function and V Vlogp(x,,1 Xo, 0),
Entropy of q(xo xx).
F shows the relative L2 distance, for both datasets. Figure 7 similarly show the distribution (random 10,000 trials) of the cosine similarity, and Figure g shows the entropy. Dashed curves indicate the bounds evaluated in eq. (1 and eq. 32
These figures show that the range of intermediate region between Phase (1) and Phase (2) will not have impact in practical situations since do not evaluate the neural network Sa(;-)) in this range so many times (i.e., a₂ ~ 10-3 to 10-1 <> Vt 0.999 to 0.9). Moreover, the approximation accuracy still very high even in this region. Furthermore, although MNIST and CIFAR-10 quite "low din-eimensional" for real-world images, approximations are established with such high accuracy. Therefore, is expected be established with higher accuracy for more realistic images.
18

---

Under review as a conference paper at ICLR 2023
Relative L2 Ervor (MNIST)
Relative Error (CIFAR-10)
0.25
0.25
Phase
Phase
Poase
Phase
0.30
Relaive La Ervor
0.90
0.10
0.10
0000
0000
0.00
0.00
(a) MNIST training data. 784-dim (= 28 28)
(b) CIFARIO data. 3072-dim (= 32 32 x 3)
Figure 6: Relative L2 distance (4) between the empirical score function eq. (45 and the single point approximation, depending the noise level 14. Dashed curves indicate the upper bounds evaluated in eq. 31
cosine similarity (MNIST)
cosine similarity (CIFAR-10)
1.00
1.00
sm
0.98
0.98
0.97
Phose
0.97
Phose (1i
Plase 0
10
10
10
10
10
(a) MNIST training data. 784-dim (= 28 28)
(b) CIFARIO data. 3072-dim (= 32 32 x 3)
Figure 7: Cosine similarity (1) between the empirical score function and the single point approx- imation depending on the noise level 2/. Dashed curves indicate the lower bounds evaluated in eq. (32
Entrgy of g(xg N,) (MNIST)
Entropy of g(xp x,) (CIFAR-10)
10.0
Phase (1)
Phase (2)
Plase
7.0
Earoogy
Dnnrogy
2-
0.0
0-
(a) MNIST training data. 784-dim (= 28 28)
(b) CIFARIO data. 3072-dim (= 32 x 32 x 3)
Figure 8: Entropy g(xo) Xx) depending on the noise level 2/-
19

---

Under review as a conference paper at ICLR 2023
B ON THE IDEAL DERIVATIVE APPROXIMATION
Thus, we can assume that the single point approximation almost always holds practically.
So(Xi,t) t) model
almost equal
Xx logp(xe, / x ,) X, v1 1/x00
(1)
Xx, logp ((x,,)) 1)
~
Vut
Vt
Therefore, we may also expect that the similar approximation will be valid for their derivatives.
Of course, strictly speaking, such an expectation is mathematically incorrect. For example, let g(x) = f()) + sin wa, then the difference (zx) - f(z) = E sin wT goes to zero as + 0, but the difference of derivatives (x) - (x(x) = ew cos war: does not if w - oc faster than 1/e. the error between them in the Fourier domain is written
[(w) G(w) - F(w), then the T2 error ([w)()) [2 x const (Parseval's theorem). other words, the single point approximation does not necessarily imply the ideal derivative approximation. to be mathematically rigorous, it must be supported by other nontrivial knowledge on the data manifold. This nontrivial leap is the most important "conjecture" made in this paper and its theoretical background should be more closely evaluated in the future.
between the derivatives is [[g (x))] =
B.1 DERIVATION OF THE "IDEAL DERIVATIVES"
Because of the discussion in $ A the true score function Vx, log p(x,, t) is finely approximated by a single point approximation Vx, log p(x,, xo, 0). Now may also assume that the derivatives of both will also be close.
In this paper, we interested in the Taylor expansion of the following form (see also $ E.1.1
h*
(xt, h) (xo, 0) +
k! (0, + a(x,,t)
Vx.) (xxx, t)
(46)
1=
t=
If the function ub(x,, () is separable in each dimension (i.e., 0z thy 0 for i * j). the following relation holds,
(a(x,,)) Vx, ()...)) = a(x;, t) o Vx, O (x,, t).
(47)
where @ the element-wise product or operation. If a(x,, t) is also separable each dimensior the Taylor series is formally rewritten as follows,
(x(x,, () = (xo. 0) + Σ (10, + a(xx, o Ox.)) ((xx,))
(48)
t=
where 0x, == Vx, O is the element-wise derivative operator. This is formally the same as the 1-dim Taylor series. Therefore, it is sufficient to consider the I-dim Taylor series first, and parallelize each dimension later. Thus the derivatives actually need are the following two.
0x, So I) Vx, ④ So(Xe, (), 0,So(x::) () = (10,) + So(x,,
(49)
B.1.1 SPATIAL DERIVATIVE 0z (x t) = Vx. @ (x(,
Let us first compute the spatial derivative of the conditional score function.
(a Vxa Vx)((-/1) logp(x;, Xo. 0))
Σ Carday
x,
14
Viz
E, (d,a, + a,),)), If is separable in each dimension, the (i j) terms vanish, and (a v)² -, Σ, ard 0,0,), If the function (xx,, t) is separable in each dimension, then (a V) =,(a,,,a + 0,,00, )on aa v.. Thus
"In general, (a 2 a,,) = E, (,) a,
formally write (a V) (a Vevv@ataee a
V)
② (V a + O V & V) O a v (a o V) @ &b (a @ vol' (a 0 2d (Note that operator (a V) is scalar while (a ax) is d-dim vector.) We can similarly show (a V) 5b (a @a for k > 3.
20

---

Under review as a conference paper at ICLR 2023
Ex a.0.)))) (xi - v1 1420)
1
v
(E,,de.) (x, v
a.0 [.0
(x0 v1 M230 )
1
:
v
[I a. (xte" v1 vI-1 11000 )
(a10 (2 1 -V1- - -V1-vzo'
1
V4
agd 22d
(Ire" -111 v1 - v200
1
Vv2 a=
1.
(50)
V/z
Vuz
ad
Here, we used the notation x, to denotes the i-th component of a vector x Note that up to this point in the discussion, there have been no approximations, but strict ones.
Now 1et us consider the approximation. Because of the single point approximation, we may assume that the derivative the integrated score function will also be approximated by the derivative of the conditional score function, i.c.,
(a
logp(x., ()) (a vx.)(-vi, Vx, log 1og,( Xo, 0)).
(51)
As the neural network So(x,, is trained that it approximates the integrated score function, can also assume the following relation,
1
(a V (x,, (a Vx,( V t xo. 0)) =
(52)
Thus we have obtained the ideal spatial derivative of the neural network.
We can also formally write the spatial derivative follows using the above notation,
1
C (0x, So t)) = 2
1.
(53)
Vit
We can also write it as
1
0x, Sx(
(54)
Vtr
1.1.2 TIME DERIVATIVE -x,,
Next, let us compute the relation
O.(-vv.V Vx, logp(xe: / Xo: 0)). During the computation, Xo is replaced by
X
VI- (xt + us Vx, log px,,( Xx, 0))
(55)
We also use the following relations between Vt, Sr, which is immediately obtained from the definition of Dt,
i (1 3r.
(56)
Using the above information, may compute the temporal derivative of the conditional score function as follows.
0 Vr Vx, 1xx,( Xo, 0))
= - Xx
ViX0
21

---

Under review as a conference paper at ICLR 2023
(-2mer 3/2
/ ve) -1/x1/20)) - -xx: - -v- mxxo)
Vu
vI X0 - -xx - V1 -1- v1-1xx0)
11
20 3/2
1
-x +
-I- Vt X0
201
3/2
Vt
1
X1 +
(Xt Vx, p(xe: / xo. 0))
21 3/2
VI - Vt VI -1-W Vt
it
-1+ 1- V4
+ 1- (ve logp(xe,( Xo. 0))
2 3/2
ip
(vix² + 14 Vx, log p(x,, Xo: 0))
20/ 3/2
1- 1/2
Be (vixe + u. Vx logp(xx, / xo. 0))
2v 3/22"
Be
(x, + Xx, logp(x,, Xp, 0)).
(57)
2
(Note that this caleulation is exact, and no approximation injected.) Because of the single point approximation, we may assume
- vv xx, logp(x;,1)) - - 0,(-von x, logp(xe; / xo, 0))
Be (x, Vx. log 1ogg((x;:) xo: 0))
2/4
Br (x+ + Vx, log (x),(xx;
(58)
2/vV
and therefore, also assume that the temporal derivative of the neural network is approximated as
a, a,So(xx, ~
Be
xt
So(xI.1)
(59)
2/12
Vit
The "derivatives" " have some good points. For example, the partial derivatives commute,
0x, O,S 1) = -,,0x. So(x,,
(60)
22

---

Under review as a conference paper at ICLR 2023
B.2 COMPARISON OF THE EMPIRICAL SCORE DERIVATIVES AND IDEAL DERIVATIVES
Let us empirically validate that idela approximation using real data similarly as above. However, since the equations will become very complicated if evaluate the exact empirical score derivatives, we instead used finite differences as the ground truths. That is, let S(x, () be the routine that computes the empirical score function follows,
S(x, 1)
[a(xa
xx) V (xx, Xo, 0)];
(61)
and we evaluated the empirical score derivatives by the finite differences as follows
0,5 R2 S(x,, + E) S(x,, t)
Empirical Deriv:
(62)
(a Vx )5x S(x, + sa, -(xx, () where ~ (O,
Empirical x; Deriv:
(63)
where e should be sufficiently small value, and used E empirical derivatives with the ideal derivatives using MNIST and CIFAR-10.
10-3 here We compared these
Be
1 So(().1)) V
Be
X,
1/X00
Ideal É Deriv:
0.S0
Xt
2 2/4
2 2/
44
Ideal X Deriv:
(a. Vx )s
Vii
As the ideal derivatives require the specific function forms diffusion and variance schedules, we tested following noise schedules.
Linear schedule We first tested on the linear schedule eq. 76 where Bo = 0.1 and B = 9.95.
This is the same schedule as the one used in the main text.
Figure D shows the relative Lg error and the cosine similarity between the ideal derivative eq.
21 and
the empirical derivative eq. (62), in which is observed that they are very close when 0 - 0.5, while the approximation accuracy decreases as / increases. However, even in that case, there tends to be an overall positive correlation. can also be observed that there is an error that seems to originate from the singularity of time origin when sz 0. (See also [.2)
For the x derivative (Figure D on the other hand, we can confirm that the errors between the ideal x derivative eq. (1 and empirical x derivative eq. 62 are generally very highly correlated, except around t 0.5.
Modified tanh schedule We also tested another noise schedule, the modified tanh schedule eq. (79) which does not have the singularity at the time origin. The parameters A, k were determined so that Po 0.001 and /1 0.999. Figure II and Figure 12 [ the results. In this case, the overall trend is similar to the linear schedule, but we can observe that the singularity of the time origin of the t derivative is eliminated.
'To verify the empirical x, derivative, let us consider simple case of variabe function fa, y, 2). As its total derivative is df f.fdæ + fyy + fzz. we have f(x + a, y + b, + c) f(x v. (od: + 60, + =c0:)(zz, for small a, b. Let a ea' = et' and = ed' then f(r + ea' . + ([x,v. E(a'a, b'a d' 2). Therefore,
ab', + ec') (a'a, + B'a + cd..) lim (f(x ea', y eb', , 2 + sc() f(z, y, 2)).
write the spatial derivative as
23

---

Under review as a conference paper at ICLR 2023
Relative Error of deriv. (MNIST)
Relative L2 Error of deriv. (CIFAR-10)
Relasve Io Ervor
cosine similarity of deriv. (MNIST)
cosine similarity of f deriv. (CIFAR- 10)
0.5
0.5
sim
0.0
0.0
0.0
0.2
0.4
0.5
0.5
1.0
0.0
0.2
0.4
0.5
0.5
1.0
Figure 9: Relative L2 Error ()) and Cosine similarity (1) between the empirical & derivative and the ideal derivative. (Linear schedule.)
Relative 1.2 Froor of x deriv. (MNIST)
Relative 1.2 Frror deriv. (CIFAR-I0)
0.4
Relative J. Eioor
0..
0.2
0.2
0.0
0.0
0.6
0.5
1.0
0.0
0.2
0.4
cosine similarity of deriv, (MNIST)
cosine similarity deriv, (CIFAR 10)
1.00
0.98
0.58
0.06
0.96
006
0.94
0.94
0.0
0.2
0.4
0.5
0.8
1.0
0.0
0.2
0.4
0.5
0.8
1.0
Figure 10: Relative L2 Error (.) and Cosine similarity (1) between the empirical x derivative and the ideal x derivative. (Linear schedule.)
24

---

Under review as a conference paper at ICLR 2023
Relative Error of deriv. (MNIST)
Relative L2 Error of deriv. (CIFAR-10)
Relasve Io Ervor
Roluuse
cosine similarity of deriv. (MNIST)
cosine similarity of f deriv. (CIFAR 10)
0.5
0.5
0.0
0.0
0.0
0.2
0.4
0.5
0.5
1.0
0.0
0.2
0.4
0.5
0.8
1.0
Figure 11: Relative Le Error (L) and Cosine similarity (1) between the empirical derivative and the ideal derivative. (Modified tanh schedule.)
Relative 1.2 Froor of deriv. (MNIST)
Relative 1.2 Error deriv. (CIFAR-IO)
0.4
Relative J. Eiror
01.4
0.2
0.2
0.0
0.0
0.4
0.5
1.0
0.0
0.2
0.4
cosine similarity of deriv, (MNIST)
cosine similarity deriv, (CIFAR 10)
1.00
0.98
0.58
0.06
0.96
006
0.94
0.94
0.0
0.2
0.4
0.5
0.8
1.0
0.0
0.2
0.4
0.5
0.8
1.0
Figure 12: Relative L2 Error (.) and Cosine similarity (1) between the empirical x derivative and the ideal x derivative. (Modified tanh schedule.)
25

---

Under review as a conference paper at ICLR 2023
B.3 THE DERIVATIVES (s (-f,), L3( f.), L:( Gx(- Gz(9)
The computation of the derivative (- f,,,, Le( (g)g, Ga(-Ea), Gg(9) does not require any particular nontrivial process. All we have to do is rewrite a term every time we encounter derivative of So(#,, 1) or Vt, and the rest is the level of elementary exercises in introductory calculus. To execute this symbolic computation, the use of computer algebra systems will be a good option. should be noted, however, that some implementation tricks process such custom derivatives are required (in other words, the term-rewriting system should be customized).
The results shown below. Although these expressions appear complex first glance, the code generation system can automatically generate code for such expressions.
3
Be
Bi
T (=(-)) =
St
+
So(
t)
(64)
4
2
2/
Aue
3/2
12
3
B₂
42( (x),))) =
14
+
So (Ic.t)
(65)
4
2
V
G₁( (z) =
3/22
(66)
2
3r
Leg(t) 2/B3
(67)
Geg(t) = O.
(68)
We may also compute higher order derivatives, though we do not use them in this paper except T (-f,,,
83
3,
Be
(=)
+
It
4
2
B8( 19(-12 1/ + +34
3)
33,B,
Be
So(1e;
(69)
Su
5/2
Av
3/2
2
183
3,,,9
B
33 + 4,,
L,4((E4((00,))
+
It
So (o(22,t)
8
4
2
4
VBL
Vi
(2
+ 3/,)
3v B
4,Gx(-44
4
30
(,4((-6((2), = VB.
V4
GGG((-E)((rn))) =
Lebbg(()= 233 33
4/33/2
L,G:g(t) = 0
G.Lıg(t) = 0
G2G.g(t) O.
As we can see, no factors other than integers, Tt, So f). V4, B, and derivatives of 3 appear. This also true for higher order derivatives, which can be easily shown.
26

---

Under review as a conference paper at ICLR 2023
SymPy Code Snippet for Automatic Symbolic Computation of Derivatives The following code snippet is minimalistic example of SymPy code to compute the above derivatives using the customized derivative method. We used SymPy 1.11 to test the following code snippet.
fiom sympy Import Functicn, symbols, squt, simplify
symools ("x bT) 2,
B Function("betee
define customized decivatives of Pi
mufennction):
def fdiff (self, argindex=1);
t, self.args
seturn (2 (u(L)) B(E) & & (1 vz)0.
define customized derivatives of (a(r,t)
class 9_theta (Function):
def Iselé, argindex=1):
t salf.arge argindex 1: 0/8r
1 / sqrt (nult)) elif ergindex + 0/01
return 8_thataix, (nugg:((uutt))) B(tl (2 sqgt
define
class f_flat (Puncticn): aclasthedddd
def (cls, X, thf return Bit)
8_theta(x, B(L) (nu
define differential operator I class f_latt (Functien!:
Bclassmathed def eva (cls, fxt1:
zeturn (xL.diCI(L) (x U. (xt.diCI(x)
show each derivative print (f_ fltt (s, ti)
print (simplify fl.ftat (s_flat
(,6,(z,, t): eq.
print (simplify f1at (La Ff (x,)t)))) +,, (, (); see eg.
can similarly define fg 4, G2 and compute other decivatives.
The result will look like
xB)
So (),t))(t
[Out
2 vvv(t)
(t)
d 3(1)
So(x,))
1o(1,
[Out
2
4t²(1)
2/v(1)
(t)
3() (((1)330) 3(1)
IH B(1) x
Sa(x,t)33(1
[Out 3]
+
2
8vv(1)
38. (z,
(., (1) 3S, t)ß(t) 4 3(t) St²(t)
So(z, t) l dea 3(t)
82i(t) (t)
4v
vvt(t)
and so on. Some additional coding techniques can further improve the readability of these expressions, but there will be no need to go any deeper into such subsidiary issues here.
Thus obtained symbolic expressions can be automatically converted into executable code in practical programming languages including Python and C++ using code generator, though the authors hand-coded the obtained expressions in Python for the experiments in this paper.
27

---

Under review as a conference paper at ICLR 2023
C TRUNCATED DDIM IS EQUIVALENT TO THE QUASI-TAYLOR SAMPLER
Using SymPy, we can easily compute the Taylor expansion of given function. For example, the following code
sympy. (eryy.eerte((((t-h)) h, 0, 4)
yields the result like
h² de
3(1))
d h3
((1))
((t) + h
3(E1)
+
dep
E1=t
SiEt
+ 0 (1i")
dei
2
f1=t
Similarly, using the relation i4 (1 24 Sr, we can easily compute the Taylor expansion of 1/- h as follows.
sympy. series (t-bh), h,
v(t)
4 a5.
((1))
d B(E))
((()((())
(((()
2E,
Pt- = -00)))(0)00000)0000000)0)000)))1
Ei
6
+0 +0(19)
2
2
2
Using this functionality of SymPy, we can easily compute the Taylor expansion the DDIM (Song et al. 2020a). Let recall that the DDIM algorithm is given by (15). and using our notation V/ and o vv, it can be written as follows,
DDIM: X
Ve-h X
(Vu-h
- Vi-h Ms
So(,
1 14
DODM
Then using SymPy, the Taylor expansion of Pt.h
DDIM DDIM and Mrgh
computed as follows,
Br
1
B3
3
3,3
3
Pt.
DDIM
+
h
31
+
+
o(h3),
(70)
2
4
2
4
12
2
3
12
82
3
53
83
93
B
Vi
DDIM
h
31
+
+
o(h²).
2
4
2v:
4
12
4
4
2Pr
3
(71)
Although it has been known that DDIM corresponds to the Euler method up to 1st order terms (Song et al., 22200a Salimans & Ho. 2022). this expansion gives better understanding of higher order terms.
That these exactly equivalent our deterministic Quasi-Taylo sampler eq. (23) and eq. (24)
up to 3rd-order terms. This fact may suggest that the assumptions behind the DDIM derivation will be logically equivalent our assumptions of ideal derivatives.
The advantage of the proposed Quasi-Taylor method is that we can decide the hyperparameter at
which order the Taylor expansion is truncated. the other hand, DDIM automatically incorporates terms of much higher order, leaving no room for order tuning.
28

---

Under review as a conference paper at ICLR 2023
D ON THE NOISE SCHEDULE
D.1 BACKGROUND: PICARD-LINDELON THEOREM
Let us consider I-dim deterministic system *(t) a(x(tt, f). It is well known that this ODE has a unique solution if a(x, t) is Lipschitz continuous w.r.t x and continuous w.r.t. (Picard-Lindelöf Theorem). Otherwise, ODEs often behave less favorably. (Similar Lipschitz conditions are also required for SDEs.)
Example 1. For example, the ODE à: = x00) 1 has the solution 1/(1 - 1) when and it blows up at t = 1. It is usually impossible to consider what happens after > 1 in ordinary contexts.
Example 2. Another well-known example is i = - = 0. has a solution x = -12/4, but 0 is also a solution. It actually has infinitely many solutions (ift S ln). (t to) /4 (ift > to), where to 2 0 is an arbitrary constant.
Example 3. Let us consider the following ODE
z00) 1,
(72)
1 e 1t-
which is simplified model of the Linear schedule eq.
(76
The exact solution is as follows,
Ve-I
773)
Ve(1-1)²
which diverges at = 1. In this case, a(x,t) = - x. (t-1)((1-e-t) is not Lipschitz continuous,
as the Taylor expansion of the denominator is ((111) (1-1) 1)2 + ((l 1)²), and alx, () is approximately -x((1-1) near 1, 1.
In these cases, the coefficient -) is not Lipschitz continuous. Even these seemingly simplest ODEs behave very complexly unless the coefficients are carefully designed.
In PF-ODE, the Lipschitz condition is written as follows,
Lip(f,) = 0x. (2
B: B.
So (Ze.
< co.
(74)
2/4r
Using the ideal derivative of (1 t), this condition translates
Lip(f,) 13, 1/21) = in
<0.
(75)
D.2 SPECIFIC SCHEDULES
Including this point, the necessary conditions for a variance schedule Ve will be summarized as follows.
un so that the initial density p(x,, 0) close to the true data density.
2. Mr == 1 so that the terminal density p(xr, T') is close to the Gaussian.
3. Sufficiently smooth that 8 de log(1 1)) is well defined.
In addition, Br should also be smooth so that the Taylor schemes can be used.
4. Monotonic (s < = v, S vt) to make 3, non-negative.
5. Preferably, make the drift coefficient f, Lipschitz continuous so that PF-ODE has a unique solution, i.c., Li,(E,) R //./vll < oc.
The following two scheduling functions which are common in diffusion generative models satisfy the conditions 1, 2, 4 above (the linear schedule also satisfies the 3rd condition),
Linear:
V/ 1 - E -
BB Bu + 2,34
(76)
[t/[++5)
#
tan
1/++s T4c
if O
Cosine:
V: = 1 C cos²
B
(77)
2 1155
θ
I' <TST
29

---

Under review as a conference paper at ICLR 2023
where 0 is small constant, C 1/ cos? ( +)) is constant to make My 0, and the threshold constant is θ = Br
However, these common schedules do not satisfy the 5th condition that the drift coefficient f, is Lipschitz continuous. Indeed, it is casily verified that lim4-so 1//1/ oo in both cases, since 1/-- = but v/ 0. Nevertheless, t 0 is the only singular point, and since function value derivative at = 0 is evaluated by numerical methods (except by the Runge-Kuut Kutta method), this point can practically be ignored.
Note that, we can also consider some other schedule functions such as the sigmoid function and the hyperbolic tangent, which satisfy the condition 2, 3, 4, 5 but do not satisfy the 1st condition rigorously (but if Po is less than or equal to the level of the quantization error in the data, we may consider the first condition to essentially satisfied).
Sigmoid:
V1 =
1 1+ E A(t-R))
Br Avi,
(78)
Modified Tanh:
V1 = Lanh 2 (A(1)/2),
B. = X(1) (anh(X(1)22,
(79)
where the parameter function X() has some options, such as X(t) log(1 + Aett), and A > 0,
>> are hyperparameters.
D.3 How TO AVOID THE TIME ORIGIN SINGULARITY IN THE RUNGE- RUNGE-KUTTA METHODS
When using the Heun and Classical RK4 methods, the function f, (x,, t) is evaluated time However, since the function f, (x,,t) t) contains the term proportional to 1/v it will diverge at time t 0 if the linear eq. 7 cosine schedule eq. (77) is used. The simplest way to avoid this is to replace the function f, (x. 0) with (xe, where e >> 0 is a sufficiently small constant, only when the need to evaluate the function at time = 0 arises.
The same thing could happen at f = T if the cosine schedule and DDIM were used simultaneously, but this can be handled in the same way.
If we use the sigmoid eq. (78) or modified tanh schedules, eq. 779 these problems do not occur unless the hyperparameters and k are chosen be very extreme values.
30

---

Under review as a conference paper at ICLR 2023
E SUPPLEMENT ON FUNDAMENTALS
For convenience, let us summarize some basics behind the ideas this paper. The contents of this section not particularly novel, but the authors expect that this section will give better understanding of the ideas of this paper and the continuous-time approach to diffusion generative models.
E.1 TAYLOR EXPANSION AND ITÔ-TAYLOR EXPANSION
E.1.1 TAYLOR EXPANSION OF DETERMINISTIC SYSTEMS
1-dimensional case Let first consider I 1-dim deterministic system a(() = ((x( l), where a(s)) sufficiently smooth, and let us derive the Taylor series expression the solution of this ODE. Let (px((t),t) t) be differentiable function. Its total derivative is written
04p dt + de dx 0t dr òt
1 04 dt + Oe dx dt Ot
04
dip
Or a(x, 1) dt
dt
0 (xx,t)
0
pdt.
(80)
dx
=,,
By integrating both sides from t, we have
((x(1),0)= -(((),,)) I, (L))) s)ds.
(81)
We use this formula recursively to obtain the Taylor series of the above system. Let Let((x((),)) x((). then we have
x(t) z00) + [ (L,x))(((s) s)ds
= xx + f a(x((s)
(82)
Let ((x((),)) =(x(1),)) 1), then we have
a(x(t).t) a()(0),)) I'' (L)a)(
(83)
Using the above two equations, have
xtt)= =(00 + /'' (dlx(0).)) + 1/" (2t)) dtl
=x(0)= / a(x(0), 0)dt; + S 1'' (Lta) dtyddd
(0) ta(x(0).0) 0) J/
L/d (Lsa)(a((t), dtt
(84)
We expand again the term inside of the integral by using the following relation,
(L,a)((())) = ((00((t())))0(::((((0000)) + J' (IFa))((s)) s)ds,
(85)
and obtain the expansion as follows,
<(t) (0) ta() (Z))) J/ [. 1. (o)t
(86)
Applying this argument recursively, can obtain the Taylor series of the deterministic system, if a(;;) is sufficiently smooth.
(tf) = x0)) + ta(x(0),)) +
(I,
1 ())))))
(87)
k
31

---

Under review as a conference paper at ICLR 2023
Since the above discussion is valid when the integration interval is (t, h) instead of (0, t), it can be written as follows,
h Ā.. (4
xlt + h) - x(1) ha(a()))))
(88)
Multi-dimensional case Let us consider the multi-dimensional ODE & a(x, 1). The total derivative of a smooth scalar function cx(x, 1) is written as
de di dt
Σ de dz,
04 dt + Σ
Oe dæ dr, dt
dt dt "ai(x t)
0e
0
gdt
dzi
ez.
= -02g 0t
04 + (a(x.t) (V))p dt
0 + a(x,t) v) pdt.
(89)
Let (x(,() () be vector-valued smooth function, then we immediately have
dp (0, + a(x, V) pdtt.
(90)
Using the scalar operator L, = (0: + a(x,t) t) V), we can obtain the following Taylor expansion similarly to the 1-dim case,
22
x(t) = x(0) + t(L,x)
+
(Lxx)
+...
2
1=0
t=
12
x0) + (a, a(x,t) V)x
(0, a(x,t) v)
+...
1=0
1=0
1
x(0) + la (0) 0)
(0, aa(x, 1) v)a
+..
(91)
2
1=0
The second order term (a V l is written follows,
(x,xxt) t)
Ta
0g0 24 a1
(a V)a
Σ az.z.
(92)
as(x, 1)
a ad
050 ad
In special where each dimension is separable, i.e. 0g aj diagonal, and we have
0(i # j), the above d x d matrix is
aaiO, ai
ai
mm
@1
(a v)a=
V
a
(93)
add 21 ad
ad
ir
ad
where is the element- wise product or operation. In this case, it is sufficient to consider each dimension separately, and it is formally equivalent to the 1-dim case.
1.1 ITÓ-TAYLOR EXPANSION OF STOCHASTIC SYSTEMS
For the above reasons, it sufficient to consider the 1-dim case even in the stochastic case. As is well known, Taylor expansion is not valid for stochastic systems x, a(x,, t)dt b(x,, f)dB, This because the relation "d., ~ di". This effect taken into account the celebrated Ité's lemma, ie.. the stochastic version of eq. (80).
② + a azz t) dr 0t
b(x, 1)²
8²
θ
dip
sdt +
b(xr, t)
Bd
(94)
2
dx
dr
=;G,
By using Itô's formula recursively, we can obtain the following higher-order expansion of stochastic system, which is called the Itô- -6-Taylor expansion.
zo + ('' alxe, t) + 1'' b(x,, t)dB,
(95)
32

---

Under review as a conference paper at ICLR 2023
/ (a(za:0) L (Iaa) + S (1s))((ss))d1))) sydB, dl
hh
20 +
+ //² (b(z0.0) + f (Ls)) o)ds + (Gab) s)dB. dB.
20 + a(ro, 0) I dt Jo
ph [ (I.a) + I0a th
(Gaa)
! dB + fo !aa I (..))
f"" f (G.b)(E,,
+ bb(ro.00 0)
/ dt + b(100,0) 0) f" dB.
20 + a(z0, 0)
h 1. (Z,a) 0)+ dsdt [- J. ((C,o)(0000) 0) dd
+ I.' [:' (L,b) 0) dsdBi+ {/'' 1
f'' (IG.D)(00.))+ 0) dB,dB, d.,....BddBB
20 + al 0)h + b(5g, B)
+ (zo, 0) 1 [ dadt + (Go,) 0) 1/¹ I dB,dt
(Lab)((ro: 0) f-¹ h f dsdB, + (Trb) 0) /-' S dB,dBB
+ Remainder,
(96)
where the remainder consists of triple integrals as follows,
Remainder = S.¹ te ['' [ (Taa)) u)dudsdt +....
(97)
We ignore these terms now. If we also ignore the double integrals, obtain the Euler-Maruyama scheme.
Evaluation of Each Integral in eq.
966) Let us evaluate the double integrals.
J dodt
(deterministic)
f' dB,dl
(stochastic
dsd
(stochastic 2)
dB.dB,
(stochastic 3)
Deterministic: The deterministic one, 1thh for dsdt, is easy to evaluate.
[ dsdtt 1.' h tat h2
[[dstt
(98)
Stochastic 1: Other integrals contain stochastic integrations. Let denote the first one by ž.
$ = (() I dB,dt 1"' B,dt.
(99)
33

---

Under review as a conference paper at ICLR 2023
As is the limit of a sum of Gaussian variables, i.e.,
h
lim
Bui/n
lim
Bni/nn Bn(i-1)/aa)
i=0
== ==1
n-1 lim ΣΣ W W ~ NN0, m 11
(100)
==
so 2 is also Gaussian, whose mean is O. The variance, however, requires some discussions, which shall be seen later. In addition, we shall also see that is correlated with B₁.
E Bx] +
(101)
Stochastic The second one has the correlation with the first one as follows. Here, we use the integral-by-part formula. See e.g. (Øksendal 2013 Theorem 4.1.5).
Jo" dodBi SaaB,,-/f -1," tdB, = tB: -1 Bidt hB₁ - 2.
(102)
Stochastic 3: The third one is computed as follows, using a famous formula of Itō integral.
J'" L dB.dB. = 1-²² BidB, (B) h). This is derived by substituting p(x,t) = x² and xr. B, (i.e., dæ; 0- + 1 dB,) into Itô's formula eq. (94..
(103)
Substituting the Integrals to the Itô Taylor Expansion: Let us denote w = B,, ~ NO, h), and let us shift the integral interval from (0, h) to (t, t + h). Then, may rewrite the above second order expansion as follows, which we have already seen in the main text,
1ith = xx + ha(x,, 1) wb(x,, t)
+ (L4a) 1) 12/2 + (G;a)(zr, () + (L,b)(Z1,t) t) (wh 2) + (z,,) (w² h2/2.
(104)
Covariance the Random Variables w, 8 Next, let us evaluate the the variance of ž, and the correlation between the Gaussian variables v and ž. Let first calculate the variance of (wh 2). By Itô's isometry, we have
[ sdB, |
B[(wh =) =B
= E
d
8
- 11
(105)
3
The correlation between ip and is similarly evaluated by Itô's isometry (see Oksendal 2013 of Lemma 3.1.5)) follows,
Proof
[1, dB.) (["s117)) sdB, ] / sds
E
EE
Bu
sdB, EE
h
(106)
Jo
Using the above variance and covariance, calculate the variance of as follows,
1 1/ 133-122 2 2h
h²2
E E[
[wih 2)² h2 2hr02]
h3
(107)
3
We need to find random variables w, 2 that satisfy the requirements for (co)variances that E( = h, R[wia] h22 h2/2 and R[ 1/3 ///3, and can easily verify that the following do,
i =
Vh
(108)
hNTo/2 hTT//233
122
34

---

Under review as a conference paper at ICLR 2023
Lid. 1). (O,
where u1, U2
Let us compute the covariance matrix just to sure.
w
i
=]
=
lNvi/2
Vh 0
hvh/2v3 3
U2
U1
[w]
W2]
Vh
hvh/2v3
hvh/2
Vh hv//2v/3 hv/h/2
Vh
hyh/2
E
U
U2
U1
W2
0
hvh/2v3
VI 11/h/2
vh
hyT/22
h/ hv/22/33
0
hyh/2v/3
h 12/21
(109)
12/2
12 13
35

---

Under review as a conference paper at ICLR 2023
E.2 DIFFUSION PROCESS AND FOKKER-PLANCK EQUATION
Derivation of Fokker-Planck Equation In this section, let us derive the Fokker-Plann equation
Oep(x,, () = -x. (f(x,, ()p((x))))) + Vx
(110)
from the Itô SDE,
dx, f(x,, t)dt (t))dB.
(111)
For simplicity, we consider a 1-dim case here. The following approach using the Fourier methods (characteristic function) will be easy and intuitive. See also e.g. (Cox & Miller 1965 $ 5).(Karlin &
Taylor 1981 $ 15). (Shreve 2004] $ 6) and (Särkkä & Solin 2019 $ 5).
Let us consider infinitesimal time step h. Then Fi+th is written as follows,
Te-h It + Wt We
(112)
Let us consider the characteristic function (z(w) = Eptxr (exp iwa] (where i VTI, which the Fourier transform of the density function p(ar). Because of the convolution theorem Φx+z (w) = (w)) (w), the characteristic functions of p( + h), p(x,, t), and p(w)) have the following relation,
Saitn (w) Own, (w) (w).
(113)
is easily shown that the characteristic function of the above Gaussian is given by
1 Ow, (w) = exp iiw t)h 110-55000 ( hwo?
(114)
Expanding the r...s. up to the first order terms of h, we have
1 Own (w) 1 + iw f thh gt()))h + 0(h2).
(115)
Thus we obtain
Pasta (w) Dxe (w)
(-(-iw) + (- (e(w)) + 0
(116)
h
When h + O,
0, (w) -(-iww ()) + (-iw)² 9(1)
)= (w)
(117)
2
Since (-iw) in the Fourier domain corresponds to the spatial derivative 02, in the real domain translates as follows,
ap((I. t) = -0, (f(z,, t)p((z,)) t)) + 8., g())
2 p(x,, t),
(118)
and thus we have obtained the Fokker-Plann equation. In particular, if f 0, this equation is called the heat equation, which was also first developed by Fourier.
Example: Overdamped Langevin Dynamics When the drift term is the gradient of a potential function U(-), the SDE is often called the overdampled Langevin equatior"
dx,- V x, U (x:)dt + v2DdB,
(119)
p(x)) = O. Formally writing, Integral by parts and p(-co)
e jux d
# p(n) [e ww
dx p(z)dd= = -(x)dx= ie twe p(x)dr = (-iw) le Jw
7 'The Langevin equation should actually be the following equation system, which is called the underdamped Langevin equation,
X = MV. = -vv, - (x.) 22D V2DdB
B: dt
where v. is the velocity (momentum) variahle, A is the mass of particle, and is constant called friction viscosity coefficient. In this case, the energy function should also be modified E = !M ||| + (x). Assuming that the mass is very small compared the friction, the derivative momentum V,,, can ignored (i.c., if the force F is constant, the ODE MD -yu + F has the solution -t-Mr/M + F/y, and the velocity immediately converges to f/yitt if M & 7 and the overdamped equation is obtained.
36

---

Under review as a conference paper at ICLR 2023
Gradiint feeld. -(0())
Potenial function, U(x)
Probabiliyy póx) (npl-U(x):
Figure 13: Force field (potential gradient, score function) - (x)(),, potential function U (x). and the
Boltzmann distribution p(x) ze (xx(x)) The scalar potential ((xx) given by eq. 122
Figure 14: Numerical simulation of Langevin and Fokker-Plance equations. (Top: microscopic
particle-wise dynamies) Time evolution of 1000 random samples {x; (9) 10000 obeying the Langevin equation. Paths of four particular samples are shown in color. The time evolves from left to right. (Bottom: macroscopic population dynamies) Time evolution of the density function p(x,t) obeying the Fokker-Plande equation.
where a scalar constant. Its associated Fokker-Planck equation is
0,p(x,, Vx, (Vx, p(x,, + Vx, Dp(x,, t)) = -xx Vvx, (U(x,) + D log pox,( t)) p(xx, f)]
(120)
If U (x;) + Diogologepe,. 1) is constant, the r.h.s. will be zero, and therefore, 0,p(x.,1) = 0. That is,
1 -U()// where = / -U(x)/DDds
p(x) Z
(121)
the stationary solution of the FPE, and it does longer evolve over time. Therefore, it is expected that the particles obeying the Langevin equation will eventually follow this Boltzmann distribution.
Let compare the microscopic dynamics of each particle obeying the Langevin equation with the macroscopic dynamics of the population obeying the FPE using a 2-dim toy model. The energy function we here is the following Gaussian mixture model,
(a
cos
2km 2 + (v sin
2kr 5
)a
Ux, v) log
exp
(122)
20
where a = 0.1. The diffusion parameter was D 5. Figure 13 shows the force field (U() where X (x. y), potential function U(x), and the Boltzmann distribution p(x) = -
Figure 14 shows the time evolution of Langevin and Fokker-Planck equations, where the initial density was p(Xo. 0) (0, I). The time step size was h 5 x 10-5, and the figures are plotted every 50 steps. These figures will be helpful to intuitively understand that the SDE (Langevin equation) and the FPE are describing the same phenomenon from different perspectives: individual description population description.
37

---

Under review as a conference paper at ICLR 2023
E.3 DERIVATION THE PROBABILITY FLOW ODE AND REVERSE- SDE
E.3.1 DERIVATION OF THE PROBABILITY FLOW ODE
For convenience, the derivation of the PF-ODE is briefly described here. For more general case details where g(t) is a matrix- matrix-valued coefficient which is dependent on xx, i.e. G t), please refer to the original paper (Song et al. 2020b).
Firstly, let us consider an SDE,
dx. fxx,, t)dt + g(t)dB,.
(123)
The associated FPE is,
0,p(x,,t) = - (f((x, t)p((x,))) t)) + A 2
(124)
Here, let us forcibly incorporate the diffusion term (i.e., the second-order derivative) into the drift term (first-order derivative). Noting that the Laplacian is the inner product of two V-s, i.e. A = V.
r.h.s. V f(x,, ()(x,,) 1)
2 g2) p(x)))(xx b)
- [(((x))1) p(x,,t) V "2 1)
((xxx)))((ss)))) plarss))
-- [ [(((x)))) - bog t) p(xx,t) t)
(125)
(x,()x))))))
Thus we have obtained a diffusion-fres version of the FPE of the form,
,, V (x,, t)p(x,, t)) A (0 (0x,pxx,1)..
(126)
Since it is an 'FPE', it has an associated 'SDE' as follows,
dx, = f, (x,, t)dt I.dB..
(127)
Although this 'SDE' does not give the same particle-wist wise dynamics as the original SDE, its macro- scopic dynamics of population exactly the same as the original one, since the associated FPE has not changed from the original one at all. That is, the density evolves from p(Xo, 0) to p(xr, exactly the same way.
Since the obtained 'SDE' is actually a deterministic ODE, its time-reversal is simply obtained by flipping the sign, i.c.,
dx, (-f, (x,, (-dt)
(128)
Thus the Probability Flow ODE is obtained. If the terminal random variables are drawn from p(xr, I'), then the particles obeying this ODE will reconstruct the initial density p(xo. 0) as a population.
Note that the vector field of the form j := pv is often referred to as the flux, particularly in physical contexts e.g. fluid dynamies, where p is density and is vector field (such as the velocity field), and the PDE of the form
0, V --
(129)
often called the continuity equation, which is closely related to the conservation laws; in the fluid case, the mass conservation law. Another well known example is the charge conservation law in electromagnetism, where P is the charge density and j is the current density. In the present eq. 126). j := pf, is understood as the flux and therefore, f, could he understood as a sort of velocity field, and the conservation law corresponds to the fact that the sum of probability mass is constant. This will give an intuitive understanding why the method is called the Probability "Flow" ODE.
E.3.2 DERIVATION OF R-FPE FROM KBE
In this section, let us show that R-FPE eq. (4 is the same as the KBE eq. (3) Let us consider the case where the diffusion coefficient is independent of the spatial variable, i.... g(x,,) g(t). The overall discussion is based on (Anderson, 1982 $ 5), though the following calculation is not explicitly written in the paper.
38

---

Under review as a conference paper at ICLR 2023
Firstly, the KBE is
-,p(xe, Xs; f(xs; 8) Vx.P(xe, Xx, 8) g(s)²
Ax,p(xe, Xs,
and what we need to show is that this is equivalent as the R-FPE,
s))p 1) t)xxx. g(s)2 p(x,, s xx,
-0,p(xe, xx, Vx, () Xx, p(g.
For simplicity, we derive the former from the latter. By the Bayes theorem, we may substitute each p(x,, s xt,t) with p(xe, t Xs, s)/(xx, t), and obtain
0, p(x,, Xs, s)p(x. n(x(:t)
Vx, (f(xs,t) -((s) logp(xs, 8)) p(Xx,t xx, s)p(s, s)
n(xx;;]]
+ g(5) p(xe, t Ns s)p(x,,s)
(130)
p(xitt)
To simplify the notation, let us denote p. p(x,; s). Pr|s = p(x,, Xg. and drop the variables x, and (that is, @ always means 0., and V always means x.,). Let us also recall the Leibniz rule of the divergence of vector field. Let a be vector-valued function and let u be scalar function, then the following relation hold
Leibniz rule:
(au) Vu + (V aju
(131)
Now let us evaluate each term. The l.h.s. is
(p(,,,, x,, s)p(x,, 8)) -P(p))PPs) 1000=-0000 (OPp)s)p. Pp.
(Opr)e)Ps - Pije -V. (fpa) + A 9² _ Ps (P.FPE)
AA
g²
= =0(0p:))pp - Pt|s
f. (Vps) - (V f)pp. +
(:' Leibniz rule).
(132)
The drift term is
V [ff 9² V logps)p((sPs
[ f (f- 9² Vps
V.
Pt|sPs
Ps
V [EPtiPP= (Vpa)g" Peps]
V fpttoPa] [V. (Vp.)g²²ptt P+t=]
[ (PPe.P.) + (V f)PttaPa] [(Vp.) (Vg)p(j)) + (pp))9Ppaa]
= [f (p+ Vprs + Pr|s Vp.) + (V (f) [9" (Vps) (Vpraa) + g" (A(Pp))p
(133)
Finally, the diffision term is
4
P(joPs) V. v P(joPs)
V. P(V(Pee))) (Vp.)Pua)
2
(p,Apri + P: App + 2Vp, Vpp.).
(134)
2
TThis is also written as div(au) = agaddu + u(dva)). This is just expressing the relation, Σ, (((a,))
Σ, (a,()) (a...)u).
39

---

Under review as a conference paper at ICLR 2023
Now let us simplify the R-FPE using the above relations.
Opalc + V. [[f Pogg + A g²
2 Pale
x(P((PPa) V. [[f g'V A ge
Pr(oPs)
92
=(0pq))p + PPls -P-VP.. fips
2 Ps
+ [r. (ps Vptj= + VTa) + (V PijpPs
-[92(Vpa)- (app)~(VPpyy) -VVp()) + -2(App)
g² (p. Apt= + Ap. + 2Vp. Vptee)
2
g²
(8pt)a)ps + (p. Vptss) + (Apeis)Pe
2
= Ps (pp f Upps + g²
2 (AP))))
(135)
implies that if the following relation holds, then the R-FPE also holds,
O.Pule + Vpq== (Apps) 0.
This equation is nothing other than the KBE. Hence, KBE implies R-FPE.
We may also do the same thing for the case where the diffusion coefficient g(x,, () is dependent on the spatial variable x,. this case, the target R-FPE is modified as
Opelt v. v. [(r
1 v(g2P))) Ps
g² 2 Pole
Pole
gg
= V [(f - Vg² loggog) p.) Paje] + A
Pals
(136)
2
but we can also show that the KBE implies this relation. The use of computer algebra systems such as SymPy can reduce the calculation effort for checking this fact.
40

---

Under review as a conference paper at ICLR 2023
E.4 ON THE CONVERGENCE OF NUMERICAL SDE SCHEMES.
Let us introduce two convergence concepts which are commonly used in numerical SDE studies. See also (Kloeden al. 1994 $ 3.3, $ 3.4)
Definition (Strong Convergence). Let X2 be the path at the continuous limit h 0, and X2 be the discretized numerical path, computed by numerical scheme with the step size h > O. Then, is said the numerical scheme has the strong order convergence " the following inequality holds for a certain constant K₃ O,
E|[x x11l S K2h
(137)
Definition 2 (Weak Convergence). Similarty, it is said that the scheme has the weak order of convergence B, if the following inequality holds for any test functions o(.) in a certain class of functions, and certain constant Kw > O,
[[)((x)] -l(o()) K.h
(138)
It is known that the Euler-Maruyama scheme has the strong convergence of order 0.5, and weak order of | in general. However, for more specific cases including the diffusion generative models that the diffusion coefficient g(t) is not dependent on x,, the Buler-Maruyama scheme has little better strong convergence of order y = 1.
The strong convergence is concerned with the precision of the path, while the weak convergence with the precision of the moments. In our case, we are not much interested whether a data No generated using finite h 0 approximates the continuous limit & (h - 0) driven by the same Brownian motion. Instead, we are more interested in whether the density p(x.) of the samples generated with finite step size h 0 approximates the ideal density p(xn), (h - 0) which is supposed to approximate the true density p(Xo). In this sense, the concept of strong convergence is much important for us, but the weak convergence would be sufficient.
41

---

Under review as a conference paper at ICLR 2023
E.5 RUNGE-KUTIA METHODS
Let us briefly introduce the derivation of Runge- Kutta methods for 1-dim ODE å = f(x, t). We are interested in deriving the following derivative-free formula
Σ hbik; + o(hP), where k [
Σ haigki, hoa, p. (139)
5ith Tt
The array of coefficients are often shown in the following form, which called the Butcher tableau.
C2
a21
(140)
Gn
an1
an2
On
b1
b2
bn-1
bn
Now the problem is how to design each value in this tableau. First, for simplicity, let us consider the following 2nd order case,
Tith = It + hbyk1 + hb2k2 + o(h²)
b f(z,,t)
K2 = f(x + haz1ki, ho2),
or,
Tith = xa + hhbff t) + hb2ff + haz1 far, t), ho2) + o(h²).
(141)
Noting that the Taylor expansions of the l.h.s. and the third term in r.h.s. are written as follows,
h2
(l.h.s.) xz hf(xx,
2 (f ff') + o(h?)
(142)
(third term of r.h.s.) hbb2 (s(r)s)) 1) + haz1 h (x(,,))))))))))))))))) + he2 he2l =((1).
(143)
we obtain the following relation.
2 (f (fr)
2e + 1f( t) +
2
It + hhihhhhfaaaa f t) + hdz (f(z,,t) + ffr, f) + hezf(zz, t)) + o(h?). (144)
Comparing both sides, can find that the derivatives in both sides can be eliminated if the following equation are satisfidd
hf(ar, 1) = hb, S(xe; 1) + hbrs(xr;
(145)
1 (f + ()) = hb2 (han f f t) + (xe)((
(146)
By simplifying the above equations, we obtain the following relations of coefficients.
b + b2 = b2021 2'
8202 =
(147)
and thus we obtain the following Butcher tableau.
0
0
C2
021
C2
C2
(148)
6
b2
1(c2)
1/(202)
(When C2 the method is particularly called Heun's method.) Thus we have confirmed that the second order Taylor expansion of 244h expressed as eq. (14I), which is independent of any derivatives of f.
42

---

Under review as a conference paper at ICLR 2023
We can similarly consider higher-order methods. In the 3rd-order case, following relations automatically obtained after a little effort symbolic math programming using e.g. SymPy,
1
11 b2 by 1, b202 bycs
b2c2 byc
032021
032
2
3'
6²
6'
1
021b2 + (a31 + 032)bs = 2
a by + (032 + a31) Pb
021b2c2 + (a31 + 032 33)3
3
3
By solving this equation system using SymPy, the following Butcher tableau is obtained.
0
0
c2(7 2/3)
C2
021
cs(# 02)
eg (3e Seg ca)
Og (C2 Ca)
(149)
C3
031
a32
(3rg
2)
b
b2
b3
20g
2ng
30ges
6o2 (02 C3)
6og (0g 03)
It will naturally become much more complicated when we consider higher-ordes methods, though can be simpler if substituting specific values. In the 4th-order case, when the Butcher tableau is written as follows, the method is particularly called the Classical 4th-order Runge- Kutta (RK4) method.
0
1/2
1/2
C3
031
032
1/2
1/2
(150)
042
1
0
0
1
b1
b2
by
b1
1/6
1/3 1/3
1/6
Numerical Example Let us consider toy example 5 xiint, vo 1. The exact solution is It e--cost Figure 15 compares the numerical solutions of the Euler, Heun, Classical RK4 and Taylor 2nd methods, where the step size is h 0.5. It is clearly observed that the Euler method immediately deviates significantly from the exact solution but higher-order methods follow it for longer periods; while 2nd-order methods (Heun and Taylor 2nd) gradually deviate from the exact solution, the Classical 4th-order Runge-Kutta method more finely approximates the exact solution.
Exact
Hean
Toylor 2nd
10
15
Figure 15: Comparison of numerical solutions of the ODE = x sin 4, 40 = 1.
Note that the Taylor 2nd method is given follows,
h2 (1x cost + sin"
F xz Iu(2a sin
(151)
derivative
In this case the computation of the derivative (0, + sin 10 t0,) (ar t) tractable. However, will be infeasible if this term is little more complicated. The Runge- Kutta method is advantagcous if the derivatives are not easy to compute.
43

---

Under review as a conference paper at ICLR 2023
F PSEUDOCODE
F.1 QUASI-TAYLOOR SAMPLER
The following pseudocode shows the proposed Quasi-Taylos sampler.
Algorithm 1 (Quai-TTalor Sampling Scheme with Ideal Derivatives
Require:
Trained neural network model Sa (x,,,,, c)
The conditioning information is optional.
Data size > 0 Int
Total time T
Float
Number steps > Int
Variance schedule function in ... Float - Float
Diffusion schedule function B = ai log(1 - v1) (Optional) Conditioning information
... Float - Float
Step size schedule /1, AN Begin
List of Float
x NV(0,, La)
Draw d-dimensional Gaussian noise with unit variance
for to N do
[Step size and time parameter]
h + hn
h
[Compute the Noise Level]
Compute, or get from lookup table, the values of vi, Br, 3r, B.
[Compute Coefficients
Compute p² using 23
Compute 1 using eq. 24
[Update Data]
x p'x + H Sa(x L./VVv
eq.
22
end for Clip outliers of that c.g. -1 11
End
Output: x
Note that, Ps and μb, as well as Vt, At, 3r and their derivatives, are only dependent on T' and h. Therefore, they can be pre-computed before the actual synthesis.
44

---

Under review as a conference paper at ICLR 2023
F.2
QUASI-ITÓ- QUAS--IÓÓ-TAYLOR TAYLOR SAMPLER
The following pseudocode shows the proposed Quasi-Itó-Taylor sampler.
Algorithm 2 QQusi--tô-Taalln Sampling Scheme with Ideal Derivatives
Require:
Trained neural network model Sa (.))
The conditioning information is optional.
Data size >> 0
Int
Total time T > 0 ... Flout
Number of steps N > .-- Int
Variance schedule function
Float - Float log(1 Vs).
Diffusion schedule function
Float Float
(Optional) Conditioning information
Step size schedule h1, ,hN Begin
... List of Float
X ((2,,, Ia)
Draw d-dimensional Gaussian noise with unit variance
for n to N do
[Step size and time parameter]
1t f- hn
L 1 T
hu
[Compute the Noise Level]
Compute, get from lookup table, the values of in Br, 3
[Compute Coefficients]
Compute of using eq.
Compute 4d using eq.
26
[Draw Correlated Driving Noise]
n N then n' = 02
I No noise is injected at the final step
else
w~ N,O,, La),
u~ NO La), z=
2 w+ 2/3
See eq. 108 and Theorem
Compute n² using eq, (27) and above.
end if
[Update Data]
p'x + n t,
eq. 25
end for
Clip outliers of that e.g. 5
End
Output: x
45

---

Under review as a conference paper at ICLR 2023
F.3 TRAINING
The following pseudocode shows an example for the training of diffusion-hasee generative models.
Algorithm 3 Training of Diffusion models
Require:
Data size > 0 ...
Training Data D {xx] , where x E Rd
The conditioning information is optional.
Neural network model t,) parameterized by 8.
Neural network optimizer (Adam) and parameters (c.g. learning rate)
Total time T >> ...Fa
t necessarily same as the one for the synthesis.
Variance schedule function ... Float + Float
# Not necessarily the same as the one for the synthesis.
Batch size b > 0 ... Int
Begin
for sufficiently many times until convergence do [Draw Batch]
batch + []]
# Empty list
for b times (batch size) do (Xu: c)
Draw data from the of training data; the conditioning is optional.
f t~ Uniform(O, T)
Draw a time parameter & from the uniform distribution.
Compute vs
W ~ (O2,, Ia)
Draw d-dimensional Gaussian noise with unit variance.
Append the tuple (Xo, t, M4 w) to batch end for
Foowwrdd Computation] ([xo), [c). [t). [w]
<- batch
[s] 4 So([v v) @ [xo] [vv] @ [w], [t], [c]
See also eq. 110
[Back-propagatton]]
L Enanch [[][ [[] [s][[]]
Loss. See also eq. 10
Compute VoL by propaproop..tion.
0 +- Adam(O, VeL)
I Update the parameters Sp()))))
end for
End
Output: The network parameter
46

---

Under review as a conference paper at ICLR 2023
G ADDITIONAL RESULTS OF IMAGE SYNTHESIS
G.1 SAMPLING EXAMPLES
Fuler
DDIM
Heun
RK4
Taylor (2nd)
Taylor (3rd)
Buler-Maruyama
Itō-Tage l6yTaylor
Buler
DDIM
Heun
RK4
Taylor (2nd)
Taylor (3rd)
Euler-Maruyama
Itô- No-Taylor
Euler
DDIM
Heun
RK4
Taylor (2nd)
Taylor (3rd)
Euler-Maruyama
Itô-Taylor
Figure 16: Comparison of the image synthesis process. The dataset CelebA (64 x 64). All the conditions but the sampling algorithm are the same including the random seed. The number of refinement steps is 16.
47

---

Under review as a conference paper at ICLR 2023
G.2
RANDOM SAMPLES
(a) Euler
(b) DDIM
(c) Heun (NFE=2)
(d) RK4 (NFE=4)
(e) Taylor (2nd)
(f) Taylor (3rd)
(g) Mallurraa
(h) Itö-Taylor
Figure 17: CelebA (64 x 64) synthesis samples. N= 4.
Euler
DDIM
(c) Heun (NFE=2)
(d) RK4 (NFE=44
(e) Taylor (2nd)
(I) Taylor (3rd)
(g) Buler-Maruyee
(h) Itô-Taylor
Figure 18: CelebA (64 64) synthesis samples. N
48

---

Under review as a conference paper at ICLR 2023
Euler
DDIM
(c) Heun (NFE=2)
(d) RK4 (NFE=4)
(e) Taylor (2nd)
(i) Taylor (3rd)
(g) Euler-Maruyama
(h) Itô- Itô-Taylor
Figure 19: CelebA (64 64) synthesis samples. 16.
(a) Euler
(b) DDIM
(c) Heun (NFE=2)
(d) RK4 (NFE=4)
(e) Taylor (2nd)
(f) Taylor (3rd)
(g) Fuler-Ma
(h) Itô-Taylor
Figure 20: CelebA (64 64) synthesis samples. N 20.
49

---

Under review as a conference paper at ICLR 2023
Euler
DDIM
(c) Heun (NFE=2)
(d) RK4 (NFE=4)
(e) Taylor (2nd)
(i) Taylor (3rd)
(g) Euler-Maruyama
(h) Itô- Itô-Taylor
Figure 21: CelebA (64 64) synthesis samples. 30..
(a) Euler
(b) DDIM
(c) Heun (NFE=2)
(d) RK4 (NFE=4)
(e) Taylor (2nd)
(f) Taylor (3rd)
(g) Fuler Maruyama
(h) Itô-Taylor
Figure 22: CelebA (64 x 64) synthesis samples. = 100.
50

---

Under review as a conference paper at ICLR 2023
Euler
DDIM
(c) Heun (NFE=2)
(d) RK4 (NFE=4)
(e) Taylor (2nd)
(i) Taylor (3rd)
(g) Euler-Maruyama
(h) Itô-Taylor
Figure 23: CIFAR-10(32 x 32) synthesis samples. N =44 4.
(a) Euler
(b) DDIM
(c) Heun (NFE=2)
(d) RK4 (NFE=4)
(e) Taylor (2nd)
(f) Taylor (3rd)
(g) Fuler-Maruyamy
(h) Itô-Taylor
Figure 24: CIFAR-10 10 (32 x 32) synthesis samples. N = 8.
51

---

Under review as a conference paper at ICLR 2023
Euler
DDIM
(c) Heun (NFE=2)
(d) RK4 (NFE=4)
(e) Taylor (2nd)
(i) Taylor (3rd)
(g) Euler-Maruyama
(h) Itô-Taylor
Figure 25: CIFAR-10 (32 x 32) synthesis samples. N 20.
(a) Euler
(b) DDIM
(c) Heun (NFE=2)
(d) RK4 (NFE=4)
(e) Taylor (2nd)
(f) Taylor (3rd)
(g) Fuler- Maruyama
(h) Itô-Taylor
Figure 26: CIFAR-10 (32 32) synthesis samples. 100.
52


Reviews:
[{"review_id": "xSts9RBEriG", "reviewer": "Reviewer_rnag", "summary": "", "questions": "", "limitations": "", "rating": 3, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "In this paper, the authors deal with the acceleration of denoising diffusion models. In particular, they propose improved samplers with higher order in order to reduce the number of steps requires at sampling times. The acceleration proposed in that paper is described for both Ordinary Differential Equation (ODE) and Stochastic Differential Equation (SDE) flows. It is based on the higher-order integrators [1]. In order to efficiently compute the derivatives the authors identify the score $\\nabla \\log p(x_t)$ with the conditional score $\\nabla \\log p(x_t|x_0)$ which has a tractable expression. Doing so they are able to compute efficient approximations of the derivatives. The theoretical and methodological study is complemented with experiments on CelebA 64x64 and CIFAR10.\n\n[1] Kloeden, Platen, Schurz - Numerical Solution of SDE through computer experiments", "summary_of_the_review": "While I think that the issue of accelerating diffusion models is interesting, I think that the paper has several important flaws from a methodological and experimental point of view. ", "strength_and_weaknesses": "STRENGTHS:\n* One of the strenght of this paper is that it introduces a novel way to sample from diffusion models using the theoretical of higher order integrator. To the best of my knowledge this approach is novel.\n* I also found the use of the approximation of the score by the conditional score to be quite interesting.\n* The paper is well-written and the ODE, SDE and reverse ODE, SDE are clearly introduced.\n\nWEAKNESSES:\n* I think the authors make misleading claim regarding the Backward Kolmogorov Equation (BKE) (by the way, maybe I missed it in the main text but FPE and BKE are never defined, I just assumed that FPE was Fokker-Planck Equation and BKE was Backward Kolmogorov Equation). I could not find in [1] any reference to the backward Kolmgorov Equation but maybe I missed something here. However I disagree that the reverse-time dynamics is somehow associated with BKE. In fact FPE and BKE are dual of each other but there is no connection here with the time-reversal. The time-reversal dynamics also satisfies FPE and BKE evolutions but these are not related to the BKE of the forward process (even though one can use the BKE to establish the time-reversal SDE, see [1]).\n* More important maybe, I don't really understand why the authors say that the derivatives of the learned score functions cannot be computed? It is easy to perform autodifferenciation w.r.t. the (one-dimensional) time variable and one can use vector-Jacobian-product to compute the  term $a . \\nabla s_\\theta$. I might be missing something but it seems that these issues can be leveraged by the careful use of automatic differenciation.\n* Proposition 1 is very poorly worded. What are these \"many cases\"? I understand that here the authors are trying to justify their method but this is too hand-wavey. As of now the statement is too imprecise. What are the edge cases? How could we extend this proposition to hold in a more general setting?\n* The deterioriation of the method noted by the author \"the proposed Quasi-Taylor methods tend to give good results around $N=16,20$ and the FID deteriorates from there\" is quite worrisome. This deterioriation which to me is a key limitation of the work should have been investigated in greater depth.\n* The FID results provided by the authors are quite surprising. For CelebA a FID of 20 something is very large. In the original DDIM paper [2], the authors reported way better results with FID around 17 even for 10 steps and FID around 6 for 100 steps. Can the authors explain this stricking discrepancy between the announced numbers and the reported numbers in [2]? \n* I think that the authors do not discuss important part of the literature dealing with acceleration of diffusion models. I understand that the authors are not going to compare themselves with knowledge distillation approaches like the one of [3] because these approaches are quite different from the one considered in that paper which is focused in better samplers but improved samplers have already been proposed in the literature, like [4,5] for instance. Comparisons with these methods are important and omitted here.\n\n[1] Song, Sohl-Dickstein, Kingma, Kumar, Ermon, Poole - Score-Based Generative Modeling through Stochastic Differential Equations\n\n[2] Song, Meng, Ermon - Denoising Diffusion Implicit Models \n\n[3] Luhman, Luhman - Knowledge distillation in iterative generative models for improved sampling speed\n\n[4] Liu, Ren, Lin, Zhao - Pseudo Numerical Methods for Diffusion Models on Manifolds\n\n[5] Zhang, Chen - Fast Sampling of Diffusion Models with Exponential Integrator", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "3: Some of the paper’s claims have minor issues. A few statements are not well-supported, or require small changes to be made correct.", "technical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work.", "empirical_novelty_and_significance": "2: The contributions are only marginally significant or novel."}, "extra_text": {}}, {"review_id": "YkdcIu5rLiR", "reviewer": "Reviewer_hEQV", "summary": "", "questions": "", "limitations": "", "rating": 5, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "This paper proposes an approach to improving the sample efficiency of diffusion models by incorporating high-order derivatives. As high-order derivatives are often expensive to compute, the authors propose an approximation. Empirically, the proposed approach is able to generate images using small number of sampling steps. ", "summary_of_the_review": "Given that the clarity of the paper is limited, the experiments are not very impressive or convincing, the entire paper needs to be improved for acceptance.", "strength_and_weaknesses": "**Strength:**\n1. This work considers an important direction of diffusion models: as diffusion models are known for their inefficiency in sampling, improving the sampling speed is of great interest to the community.\n2. Using high-order derivates to improve the sampling speed is an interesting and promising direction.\n\n\n**Weakness:**\n1. The writing needs to be improved. There are many equations with parameters/notations that are not properly defined in Section 3.4. The main algorithms (Algorithms 1, 2) are not properly justified or explained with rigorous proof. The method formulation is hard to follow.\n2. Proposition 1 is not rigorously stated. For instance, what does \"in many cases\" mean? What does an \"arbitrary vector\" mean? What would be the dimension of the vector? What would be the domain?\n3. Although in the left figure in figure 4, the proposed approach (Taylor 2nd, 3rd) is able to achieve better performance than the baselines, the performance degrades after 20 steps. Why is that the reason? At the same time, the DDIM performance reported in figure 2, 4 is much worse than the one reported in the original DDIM paper: in the original paper, DDIM is able to achieve an FID of 13.36 on CIFAR-10, and 17.33 on CelebA 64 using 10 steps---a performance better than all of the proposed methods using 10 steps. However, in both figure 2 and figure 4, the reported performance for DDIM is much worse. My understanding is that this can be caused by using different noise scheduling. If that is the case, the comparison in figure 2 and 4 might not be fair for the baselines.\n5. It seems that noise scheduling will affect the performance. How do you select the noise scheduling parameters?\n6. It is hard to tell the difference between the samples from DDIM and the ones from the proposed method in Figure 1. At the same time, the samples seem to have shifted color. If it is because of not having enough sampling steps, then it would be better to visualize samples with higher quality but using more steps.", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "2: Several of the paper’s claims are incorrect or not well-supported.", "technical_novelty_and_significance": "2: The contributions are only marginally significant or novel.", "empirical_novelty_and_significance": "2: The contributions are only marginally significant or novel."}, "extra_text": {}}, {"review_id": "WVfTbg-LomA", "reviewer": "Reviewer_Qqb3", "summary": "", "questions": "", "limitations": "", "rating": 5, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "This work proposes high order methods for solving both SDE and ODE, by approximating high order gradients in Taylar expansion. It seems that this method can be extended to arbitrary order by computing some constants. Experiments are conducted on CIFAR10 and CelebA 64.\n", "summary_of_the_review": " ", "strength_and_weaknesses": "Strength:\n\nThis work proposes high order methods for solving both SDE and ODE. It seems that this method can be extended to arbitrary order without much extra cost.\n\nWeakness:\n\n1. The approximation in Eq.(19) and Eq.(21) looks rough. Is there any bound of the approximation error?\n\n2. It seems this method needs to use a specific noise schedule function in Eq.(28), which makes it less flexible. Would this method work on other noise schedule functions?\n\n3. There are only experiments on CIFAR10 and CelebA 64. Besides, the used diffusion model is also weak, which has a FID around 10 on CIFAR10.\n\n4. Missing related works on speeding up diffusion models, such as [1, 2, 3].\n\n[1] Analytic-DPM: an Analytic Estimate of the Optimal Reverse Variance in Diffusion Probabilistic Models\n\n[2] Estimating the Optimal Covariance with Imperfect Mean in Diffusion Probabilistic Models\n\n[3] DPM-Solver: A Fast ODE Solver for Diffusion Probabilistic Model Sampling in Around 10 Steps\n", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "3: Some of the paper’s claims have minor issues. A few statements are not well-supported, or require small changes to be made correct.", "technical_novelty_and_significance": "2: The contributions are only marginally significant or novel.", "empirical_novelty_and_significance": "2: The contributions are only marginally significant or novel."}, "extra_text": {}}, {"review_id": "5cvim8z9QGD", "reviewer": "Reviewer_jtaQ", "summary": "", "questions": "", "limitations": "", "rating": 3, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "The authors propose a novel solver of the Probability Flow ODE for diffusion models introduced in [1]. By taking a Taylor expansion of the ODE and including higher-order terms, the solver can take larger steps. This speeds up sampling, which is a well-known computational bottleneck in diffusion models. However, the higher order terms in the Taylor expansion are themselves expensive to compute. Thus the authors propose substituting these terms with \"ideal derivatives\" which involve \n\n[1] Y. Song, J. Sohl-Dickstein, D. P. Kingma, A. Kumar, S. Ermon, and B. Poole. Score-based generative modeling through stochastic differential equations. 2021\n", "summary_of_the_review": "Overall, I find the motivation of the paper relevant. However, I find key steps of the derivation of the method problematic (i.e., the \"idealized derivative\" assumption, the \"independent dimensions\" assumption). Moreover, empirical evaluations are lacking. From Figures 2 and 4, it appears that image quality in terms of FID is not competitive compared to other faster diffusion models.", "strength_and_weaknesses": "Strengths:\n* The work tackles a relevant problem in current score-based diffusion models: slow data generation.\nWeaknesses:\n* The empirical results in the main Figure 2 are not very compelling, compared to, e.g., [2, 3, 4].\n* There is limited discussion on why the \"ideal derivatives\" are a suitable replacement for the true derivatives. Arguments (in the appendix) are mainly intuitive, with no proofs.\n* The extra theory is dense, and it is unclear whether it is worth the empirical improvements.\n\n[2] Salimans, T. and Ho, J., 2022. Progressive distillation for fast sampling of diffusion models. arXiv preprint arXiv:2202.00512.\n[3] Song, J., Meng, C. and Ermon, S., 2020. Denoising diffusion implicit models. arXiv preprint arXiv:2010.02502.\n[4] Kong, Z. and Ping, W., 2021. On fast sampling of diffusion probabilistic models. arXiv preprint arXiv:2106.00132.", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "2: Several of the paper’s claims are incorrect or not well-supported.", "technical_novelty_and_significance": "2: The contributions are only marginally significant or novel.", "empirical_novelty_and_significance": "2: The contributions are only marginally significant or novel."}, "extra_text": {}}]

Author rebuttal:

Decision: Reject