Year: 2022

Abstract:
Denoising Diffusion Probabilistic Models (DDPMs) can generate high-quality samples such as image and audio samples. However, DDPMs require hundreds to thousands of iterations to produce final samples. Several prior works have successfully accelerated DDPMs through adjusting the variance schedule (e.g., Improved Denoising Diffusion Probabilistic Models) or the denoising equation (e.g., Denoising Diffusion Implicit Models (DDIMs)). However, these acceleration methods cannot maintain the quality of samples and even introduce new noise at a high speedup rate, which limit their practicability. To accelerate the inference process while keeping the sample quality, we provide a fresh perspective that DDPMs should be treated as solving differential equations on manifolds. Under such a perspective, we propose pseudo numerical methods for diffusion models (PNDMs). Specifically, we figure out how to solve differential equations on manifolds and show that DDIMs are simple cases of pseudo numerical methods. We change several classical numerical methods to corresponding pseudo numerical methods and find that the pseudo linear multi-step method is the best in most situations. According to our experiments, by directly using pre-trained models on Cifar10, CelebA and LSUN, PNDMs can generate higher quality synthetic images with only 50 steps compared with 1000-step DDIMs (20x speedup), significantly outperform DDIMs with 250 steps (by around 0.4 in FID) and have good generalization on different variance schedules.1

heading: 1 INTRODUCTION
text: Denoising Diffusion Probabilistic Models (DDPMs) (Sohl-Dickstein et al., 2015; Ho et al., 2020) is a class of generative models which model the data distribution through an iterative denoising process reversing a multi-step noising process. DDPMs have been applied successfully to a variety of applications, including image generation (Ho et al., 2020; Song et al., 2020b), text generation (Hoogeboom et al., 2021; Austin et al., 2021), 3D point cloud generation (Luo & Hu, 2021), textto-speech (Kong et al., 2021; Chen et al., 2020) and image super-resolution (Saharia et al., 2021).
Unlike Generative Adversarial Networks (GANs) (Goodfellow et al., 2014), which require careful hyperparameter tuning according to different model structures and datasets, DDPMs can use similar model structures and be trained by a simple denoising objective which makes the models fit the noise in the data. To generate samples, the iterative denoising process starts from white noise and progressively denoises it into the target domain according to the noise predicted by the model at every step. However, a critical drawback of DDPMs is that DDPMs require hundreds to thousands of iterations to produce high-quality samples and need to pass through a network at least once at every step, which makes the generation of a large number of samples extremely slow and infeasible. In contrast, GANs only need one pass through a network.
There have been many recent works focusing on improving the speed of the denoising process. Some works search for better variance schedules, including Nichol & Dhariwal (2021) and Watson et al. (2021). Some works focus on changing the inference equation, including Song et al. (2020a)
∗Corresponding author 1Our implementation is available at https://github.com/luping-liu/PNDM.
and Song et al. (2020b). Denoising Diffusion Implicit Models (DDIMs) (Song et al., 2020a) relying on a non-Markovian process accelerate the denoising process by taking multiple steps every iteration. Probability Flows (PFs) (Song et al., 2020b) build a connection between the denoising process and solving ordinary differential equations and use numerical methods of differential equations to accelerate the denoising process. Additionally, we introduce more related works in Appendix A.1.
However, this direct connection between DDPMs and numerical methods (e.g., forward Euler method, linear multi-step method and RungeKutta method (Timothy, 2017)) has weaknesses in both speed and effect (see Section 3.1). Some numerical methods are straightforward, like the forward Euler method, but they can only trade quality for speed. Some numerical methods can accelerate the reverse process without loss of quality, like the Runge-Kutta method, but
they need to propagate forward more times along a neural network at every step. Furthermore, we also notice that numerical methods can introduce noticeable noise at a high speedup rate, which makes high-order numerical methods (e.g., Runge-Kutta method) even less effective than DDIMs. This phenomenon is also mentioned in Salimans & Ho (2022).
To figure out the reason for the performance degradation in classical numerical methods, we conduct some analyses and find that classical numerical methods may sample data far away from the main distribution area of the data, and the inference equations of DDPMs do not satisfy a necessary condition of numerical methods at the last several steps (see Section 3.2).
To tackle these problems, we design new numerical methods called pseudo numerical methods for diffusion models (PNDMs) to generate samples along a specific manifold in Rn, which is the highdensity region of the data. We first compute the corresponding differential equations of diffusion models directly and self-consistently, which builds a theoretical connection between DDPMs and numerical methods. Considering that classical numerical methods cannot guarantee to generate samples on certain manifolds, we provide brand-new numerical methods called pseudo numerical methods based on our theoretical analyses. We also find that DDIMs are simple cases of pseudo numerical methods, which means that we also provide a new way to understand DDIMs better. Furthermore, we find that the pseudo linear multi-step method is the fastest method for diffusion models under similar generated quality.
Besides, we provide a detailed theoretical analysis of our new theory and give visualization results to support our theory intuitively. According to our experiments, our methods have several advantages:
• Our methods combine the benefits of DDIMs and high-order numerical methods successfully. We theoretically prove that our new methods PNDMs are second-order convergent while DDIMs are first-order convergent, which makes PNDMs 20x faster without loss of quality on Cifar10 and CelebA.
• Our methods can reduce the best FID of pre-trained models with even shorter sampling time. With only 250 steps, our new denoising process can reduce the best FID by around 0.4 points Cifar10 and CelebA. We achieve a new SOTA FID score of 2.71 on CelebA.
• Our methods work well with different variance schedules, which means that our methods have a good generalization and can be used together with those works introducing better variance schedules to accelerate the denoising process further.

heading: 2 BACKGROUND
text: In this section, we introduce some backgrounds. Firstly, we present the classical understanding of DDPMs. Then we provide another understanding based on Song et al. (2020b), which inspires us to use numerical methods to accelerate the denoising process of diffusion models. After that, we introduce some background on numerical methods used later in this paper.

heading: 2.1 DENOISING DIFFUSION PROBABILISTIC MODELS
text: DDPMs model the data distribution from Gaussian distribution to image distribution through an iterative denoising process. Let x0 be an image, then the diffusion process is a Markov process and the reverse process has a similar form to the diffusion process, which satisfies:
xt+1 ∼N ( √ 1− βtxt, βtI), t = 0, 1, · · · , N − 1. xt−1 ∼N (µθ(xt, t), βθ(xt, t)I), t = N,N − 1, · · · , 1. (1)
Here, βt controls the speed of adding noise to the data, calling them the variance schedule. N is the total number of steps of the denoising process. µθ and βθ are two neural networks, and θ are their parameters.
Ho et al. (2020) get some statistics estimations of µθ and βθ. According to the properties of the conditional Gaussian distribution, we have:
q(xt|x0) = N ( √ ᾱtx0, (1− ᾱt)I), q(xt−1|xt, x0) = N (µ̄t(xt, x0), β̄tI). (2)
Here, αt = 1 − βt, ᾱt = ∏t i=1 αi, µ̄t = √ ᾱt−1βt 1−ᾱt x0 + √ αt(1−ᾱt−1) 1−ᾱt xt and β̄t = 1−ᾱt−1 1−ᾱt βt. Then this paper sets βθ = β̄t and designs a objective function to help neural networks to represent µθ.
Objective Function The objective function is defined by: Lt−1 = Eq [ ||µ̄t(xt, x0)− µθ(xt, t)||2 ] = Ex0,ϵ [ || 1√
αt
( xt(x0, ϵ)−
βt√ 1− ᾱt ϵ
) − µθ(xt(x0, ϵ), t)||2 ] = Ex0,ϵ [ β2t
αt(1− ᾱt) ||ϵ− ϵθ(
√ ᾱtx0 + √ 1− ᾱtϵ, t)||2 ] .
(3)
Here, xt(x0, ϵ) = √ ᾱtx0+ √ 1− ᾱtϵ, ϵ ∼ N (0, 1), ϵθ is an estimate of the noise ϵ. The relationship between µθ and ϵθ is µθ = 1√αt (xt − βt√ 1−ᾱt
ϵθ). Because ϵ ∼ N (0, 1), we assume that the mean and variance of ϵθ are 0 and 1.

heading: 2.2 STOCHASTIC DIFFERENTIAL EQUATION
text: According to Song et al. (2020b), there is another understanding of DDPMs. The diffusion process can be treated as solving a certain stochastic differential equation dx = ( √ 1− β(t) − 1)x(t)dt +√
β(t)dw. According to Anderson (1982), the denoising process also satisfies a similar stochastic differential equation:
dx = ( ( √ 1− β(t)− 1)x(t)− β(t)ϵθ(x(t), t) ) dt+ √ β(t)dw̄. (4)
This is Variance Preserving stochastic differential equations (VP-SDEs). Here, we change the domain of t from [1, N ] to [0, 1]. When N tends to infinity, {βi}Ni=1, {xi}Ni=1 become continuous functions β(t) and x(t) on [0, 1]. Song et al. (2020b) also show that this equation has an ordinary differential equation (ODE) version with the same marginal probability density as Equation (4):
dx = ( ( √ 1− β(t)− 1)x(t)− 1
2 β(t)ϵθ(x(t), t)
) dt. (5)
This different denoising equation with no random item and the same diffusion equation together is Probability Flows (PFs). These two denoising equations show us a new possibility that we can use numerical methods to accelerate the reverse process. As far as we know, DDIMs first try to remove this random item, so PFs can also be treated as an acceleration of DDIMs, while VP-SDEs are an acceleration of DDPMs.

heading: 2.3 NUMERICAL METHOD
text: Many classical numerical methods can be used to solve ODEs, including the forward Euler method, Runge-Kutta method and linear multi-step method (Timothy, 2017).
Forward Euler Method For a certain differential equation satisfying dxdt = f(x, t). The trivial numerical method is forward Euler method satisfying xt+δ = xt + δf(xt, t).
Runge-Kutta Method Runge-Kutta method uses more information at every step, so it can achieve higher accuracy 2. Runge-Kutta method satisfies: k1 = f(xt, t) , k2 = f(xt + δ 2k1, t+ δ 2 ) k3 = f(xt + δ 2k2, t+ δ 2 ) , k4 = f(xt + δk3, t+ δ)
xt+δ = xt + δ 6 (k1 + 2k2 + 2k3 + k4).
(6)
Linear Multi-Step Method Linear multi-step method is another numerical method and satisfies:
xt+δ = xt + δ
24 (55ft − 59ft−δ + 37ft−2δ − 9ft−3δ), ft = f(xt, t). (7)

heading: 3 PSEUDO NUMERICAL METHOD FOR DDPM
text: In this section, we first compute the corresponding differential equations of diffusion models to build a direct connection between DDPMs and numerical methods. As a byproduct, we can directly use pre-trained models from DDPMs. After establishing this connection, we provide detailed analyses on the weakness of classical numerical methods. To solve the problems in classical numerical methods, we dive into the structure of numerical methods by dividing their equations into a gradient part and a transfer part and define pseudo numerical methods by introducing nonlinear transfer parts. We find that DDIMs can be regarded as simple pseudo numerical methods. Then, We explore the pros and cons of different numerical methods and choose the linear multi-step method to make numerical methods faster. Finally, we summarize our findings and analyses and safely propose our novel pseudo numerical methods for diffusion models (PNDMs), which combine our proposed transfer part and the gradient part of the linear multi-step method. Furthermore, we analyze the convergence order of pseudo numerical methods to demonstrate the effectiveness of our methods theoretically.

heading: 3.1 FORMULA TRANSFORMATION
text: According to Song et al. (2020a), the reverse process of DDPMs and DDIMs satisfies:
xt−1 = √ ᾱt−1
( xt −
√ 1− ᾱtϵθ(xt, t)√
ᾱt
) + √ 1− ᾱt−1 − σ2t ϵθ(xt, t) + σtϵt. (8)
Here, σt controls the ratio of random noise. If σt equals one, Equation (8) represents the reverse process of DDPMs; if σt equals zero, this equation represents the reverse process of DDIMs. And only when σt equals zero, this equation removes the random item and becomes a discrete form of a certain ODE. Theoretically, the numerical methods that can be used on differential equations with random items are limited. And Song et al. (2020b) have done enough research in this case. Empirically, Song et al. (2020a) have shown that DDIMs have a better acceleration effect when the number of total steps is relatively small. Therefore, our work concentrate on the case σt equals zero.
To find the corresponding ODE of Equation (8), we replace discrete t− 1 with a continuous version t − δ according to (Song et al., 2020a) and change this equation into a differential form, namely, subtract xt from both sides of this equation:
xt−δ − xt = (ᾱt−δ − ᾱt)
( xt√
ᾱt( √ ᾱt−δ + √ ᾱt) − ϵθ(xt, t)√ ᾱt( √ (1− ᾱt−δ)ᾱt + √ (1− ᾱt)ᾱt−δ)
) .
(9) Because δ is a continuous variable from 0 to t, we can now compute the derivative of the generation
data xt and get that lim δ→0
xt−xt−δ δ = −ᾱ ′(t) ( x(t) 2ᾱ(t) − ϵθ(x(t),t) 2ᾱ(t) √ 1−ᾱ(t) ) . Here, ᾱ(t) is the continuous
version of {ᾱi}Ni=1 like the definition of x(t). Therefore, the corresponding ODE when δ tends to 2To achieve higher accuracy, more information is not enough. The reason why these methods achieve higher accuracy can be found in Appendix A.2
zero of Equation (9) is:
dx dt = −ᾱ′(t)
( x(t)
2ᾱ(t) − ϵθ(x(t), t) 2ᾱ(t) √ 1− ᾱ(t)
) . (10)

heading: 3.2 CLASSICAL NUMERICAL METHOD
text: After getting the target ODE, the easiest way to solve it is through classical numerical methods. However, We notice that classical numerical methods can introduce noticeable noise at a high speedup rate, making high-order numerical methods (e.g., Runge-Kutta method) even less effective than DDIMs. This phenomenon is also mentioned in Salimans & Ho (2022). To make better use of numerical methods, we analyze the differences between Equation (10) and usual differential equations and find two main problems when we directly use numerical methods with diffusion models.
The first problem is that the neural network ϵθ and Equation (10) are well-defined only in a limited area. Equation (2) shows that the data xt is generated along a curve close to an arc. According to Figure 2, most of xt is concentrated in a band with a width of around 0.1, namely the red area in Figure 2. This means that the neural network ϵθ cannot get enough examples to fit the noise successfully away from this area. Therefore, ϵθ and Equation (10), which contains ϵθ, are only well-defined in this limited area. However, all classical numerical methods generate results along a straight line instead of an arc. The generation process may generate samples away from the well-defined area and then introduce new errors. In Section 4.3 we will give more visualization results to support this.
The second problem is that Equation (10) is unbounded at most cases. We find that for most linear variance schedules βt, Equation (10) tends to infinity when t tends to zero (see Appendix A.4), which does not satisfy the condition of numerical methods mentioned in Section 2.3. This is an apparent theoretical weakness that previous works have not explored. On the contrary, in the original DDPMs and DDIMs, the prediction of the sample xt and the noise ϵθ in the data are more and more precise as the index t tends to zero (see Appendix A.5). This means original diffusion models do not make a significant error in the last several steps, whereas using numerical methods on Equation (10) does. This explains why DDIMs are better than higher-order numerical methods.

heading: 3.3 PSEUDO NUMERICAL METHOD ON MANIFOLD
text: The first problem above shows that we should try to solve our problems on certain manifolds. Here, target manifolds are the high-density region of the data xt of DDPMs, which is defined by xt(x0, ϵ) = √ ᾱtx0 + √ 1− ᾱtϵ, ϵ ∼ N (0, 1). Ernst & Gerhard (1996) show several numerical methods to solve differential equations on manifolds that have analytic expressions. Unfortunately, it’s challenging to use the above expression of manifolds. Because we do not know the target x0 in the reverse process and random items ϵ are hard to handle, too.
In this paper, we design a different way that we make our new equation of denoising process more fits with the equation of original DDIMs to make their results share similar data distribution. Firstly, we divide the classical numerical methods into two parts: gradient and transfer parts. The gradient part determines the gradient at each step, while the transfer part generates the result at the next step. For example, linear multi-step method can be divided into the gradient part f ′ = δ24 (55ft−59ft−δ+ 37ft−2δ − 9ft−3δ) and the transfer part xt+δ = xt + δf ′. All classical numerical methods have the same linear transfer part, while gradient parts are different.
We define those numerical methods which use a nonlinear transfer part as pseudo numerical methods. And an expected transfer part should have the property that when the result from the gradient part is precise, then the result of the transfer part is as close to the manifold as possible and the error of this result is as small as possible. We find that Equation (9) satisfies such property.
Property 3.1 If ϵ is the precise noise in xt, then the result of xt−δ from Equation (9) is also precise.
And we put the proof of this property in Appendix A.5. Therefore, we use:
ϕ(xt, ϵt, t, t− δ) = √ ᾱt−δ√ ᾱt xt − (ᾱt−δ − ᾱt)√ ᾱt( √ (1− ᾱt−δ)ᾱt + √ (1− ᾱt)ᾱt−δ) ϵt (11)
as the transfer part and ϵθ as the gradient part. That if ϵθ is precise, the result of xt−δ is also precise, which means that ϵθ can determine the direction of the denoising process to generate the final results. Therefore, such a choice also satisfies the definition of a gradient part. Now, we have our gradient part ϵθ and transfer part ϕ.
This combination solves the two problems mentioned above successfully. Firstly, our new transfer parts do not introduce new errors. This property also means that it keeps the results at the next step on the target manifold because generating samples away is a kind of error. This shows that we solve the first problem. Secondly, we know that the prediction of ϵθ is more and more precise in the reverse process in the above subsection. And our new transfer part can generate precise results according to the precise prediction of ϵθ. Therefore, our generation results are more and more precise using pseudo numerical methods, while classical numerical methods can introduce obvious error at the last several steps. This shows that we solve the second problem, too. We also find that their combination ϕ(xt, ϵθ(xt, t), t, t − 1) is just the inference equation used by DDIMs, so DDIMs is a simple case of pseudo numerical methods. Here, we define DDIMs as DDIMs*, emphasizing that it is a pseudo numerical method.

heading: 3.4 GRADIENT PART
text: Because we split numerical methods into two parts, we can use the same gradient part from different classical numerical methods freely (e.g., linear multi-step method), although we change the transfer part of our inference equation. Our theoretical analyses and experiments show that the gradient part from different classical methods can work well with our new transfer part (see Section 3.6, 4.2). By using the same gradient part of the linear multi-step method, we have: et = ϵθ(xt, t) e′t = 1 24 (55et − 59et−δ + 37et−2δ − 9et−3δ)
xt+δ = ϕ(xt, e ′ t, t, t+ δ).
(12) By using the same gradient part of Runge-Kutta method, we have: e1t = ϵθ(xt, t) x1t = ϕ(xt, e 1 t , t, t+ δ 2 ) e2t = ϵθ(x 1 t , t+ δ 2 ) x2t = ϕ(xt, e 2 t , t, t+ δ 2 ) e3t = ϵθ(x 2 t , t+ δ 2 ) x3t = ϕ(xt, e 3 t , t, t+ δ) e4t = ϵθ(x 4 t , t+ δ) e′t = 1 6 (e 1 t + 2e 2 t + 2e 3 t + e 4 t )
xt−δ = ϕ(xt, e ′ t, t, t+ δ).
(13)
Algorithm 1 DDIMs 1: xT ∼ N (0, I) 2: for t = T − 1, · · · , 1, 0 do 3: xt = ϕ(xt+1, ϵθ(xt+1, t+ 1), t+ 1, t) 4: end for 5: return x0
Algorithm 2 PNDMs 1: xT ∼ N (0, I) 2: for t = T − 1, T − 2, T − 3 do 3: xt, et = PRK(xt+1, t+ 1, t) 4: end for 5: for t = T − 4, · · · , 1, 0 do 6: xt, et = PLMS(xt+1, {ep}p>t, t+1, t) 7: end for 8: return x0
Abbreviate Equation (12) and (13) as xt+δ, et = PLMS(xt, {ep}p<t, t, t + δ), xt+δ, e1t = PRK(xt, t, t+ δ).
Here, we have provided three kinds of pseudo numerical methods. Although advanced numerical methods can accelerate the denoising process, some may have to compute the gradient part ϵθ more times at every step, like the Runge-Kutta method. Propagating forward four times along a neural network makes the denoising process slower. However, we find that the linear multi-step method can reuse the result of ϵθ four times and only compute ϵθ once at every step. And theoretical analyses tell us that the Runge-Kutta and linear multi-step method have the same convergence order and similar results.
Therefore, we use the gradient part of the linear multi-step method and our new transfer part as our main pseudo numerical methods for diffusion models (PNDMs). In Table 1, we show the relationship between different numerical methods. Here, we can see PNDMs combine the benefits of higher-order classical numerical methods (in the gradient part) and DDIMs (in the transfer part).

heading: 3.5 ALGORITHM
text: We can provide our whole algorithm of the denoising process of DDIMs now. According to Song et al. (2020a), the algorithm of the original method satisfies Algorithm 1. And our new algorithm of DNPMs uses the pseudo linear multi-step and pseudo Runge-Kutta method, which satisfies Algorithm 2. Here, we cannot use linear multi-step initially because the linear multi-step method cannot start automatically, which needs at least three previous steps’ information to generate results. So we use the Runge-Kutta method to compute the first three steps’ results and then use the linear multi-step method to calculate the remaining.
We also use the gradient parts of two second-order numerical methods to get another pseudo numerical method. We introduce the details of this method in Appendix A.3. We call it S-PNDMs, because its gradient part uses information from two steps at every step. Similarly, we also call our first PNDMs F-PNDMs, which use data from four steps, when we need to distinguish them.

heading: 3.6 CONVERGENCE ORDER
text: Change the transfer part of numerical methods may introduce unknown error. To determine the influence of our new transfer part theoretically, we compute the local and global error between the theoretical result of Equation (10) x(t+ δ) and our new methods, we find that x(t+ δ)−xDDIM(x+ δ) = O(δ2) and x(t+ δ)− xS/F-PNDM(x+ δ) = O(δ3). (14) If the target ODE satisfies Lipschitz condition and local error elocal = O(δk), then there are C and h such that the global error eglobal satisfies eglobal ≤ Cδk(1 + eh + e2h + · · · ) ≤ C ′δk−1. And we have that the convergence order is equal to the order of the global error. The detailed proof can be found in Appendix A.6. Therefore, we get the following property:
Property 3.2 S/F-PNDMs have third-order local error and are second-order convergent.

heading: 4 EXPERIMENT

heading: 4.1 SETUP
text: We conduct unconditional image generation experiments on four datasets: Cifar10 (32 × 32) (Krizhevsky et al., 2009), CelebA (64 × 64) (Liu et al., 2015), LSUN-church (256 × 256) and LSUN-bedroom (256× 256) (Yu et al., 2016). According to the analysis in Section 3.1, we can use pre-trained models from prior works in our experiments. The pre-trained models for Cifar10, LSUNchurch and LSUN-bedroom are taken from Ho et al. (2020) and the pre-trained model for CelebA is taken from Song et al. (2020a). In these models, the number of total steps N is 1000 and the variance schedule is linear variance schedule. And we also use a pre-trained model for Cifar10, which uses a cosine variance schedule from improved denoising diffusion probabilistic models (iDDPMs (Nichol & Dhariwal, 2021)).

heading: 4.2 SAMPLE EFFICIENCY AND QUALITY
text: To analyze the acceleration effect, we test Fenchel Inception Distance (FID (Heusel et al., 2018)) on different datasets under different steps and different numerical methods, including DDIMs, S-PNDMs, F-PNDMs and classical fourth-order numerical methods (FONs) (e.g., Runge-Kutta
dataset model
FID step
10 20 50 100 250 1000 time
Cifar10 DDIM 13.4 6.84 4.67 4.16 4.04PF 13.8 3.89 3.69 3.71 3.72
Cifar10 (linear) DDIM* 18.5 10.9 6.99 5.52 4.52 4.00 0.337 FON 13.1 7.41 5.26 4.65 4.12 3.71 0.390 S-PNDM 11.6 7.56 5.18 4.34 3.91 3.80 0.344 F-PNDM 7.03 5.00 3.95 3.72 3.60 3.70 0.391
Cifar10 (cosine) DDIM 14.5 8.79 5.86 4.92 4.30 3.69 0.505 S-PNDM 8.64 5.77 4.46 3.94 3.71 3.38 0.517 F-PNDM 7.05 4.61 3.68 3.53 3.49 3.26 0.595
CelebA DDIM 17.3 13.7 9.17 6.53 3.51
CelebA (linear) DDIM* 16.9 13.4 8.95 6.36 4.44 3.41 1.237 FON 16.0 11.6 8.13 6.70 5.14 4.17 1.431 S-PNDM 12.2 9.45 5.69 4.03 3.19 2.99 1.258 F-PNDM 7.71 5.51 3.34 2.81 2.71 2.86 1.433
Table 2: Image generation measured in FID on Cifar10 and CelebA. PFs use black box ODE solvers and we use the number of score function evaluations as the step of PFs. DDIM* is a retest of DDIM. The bold results mean the best ones using the same pretrained model. We use the 50-step, 512 batch size experiment on an RTX3090 to test the computational cost and the column time is the average computational cost per step in seconds. And we put the results of standard deviation in Appendix A.12
method and linear multi-step method). On Cifar10 and CelebA, we first provide the results of previous works DDIMs. Then, we use the same pre-trained models to test numerical methods mentioned in this paper and put the results in Cifar10 / CelebA (linear). We also use models from iDDPMs to test nonlinear variance schedules and put the results in Cifar10 (cosine). Song et al. (2020b) do not provide detailed FID results of probability flows (PFs) under different steps, so we retest the results using its pretrained models by ourselves.
Efficiency Our two baselines are DDIM and PF. DDIM is a simple case of pseudo numerical methods, and PF is a case of classical numerical methods. However, PF uses a much bigger model than DDIM and uses some tricks to improve the sample quality. To ensure the experiment’s fairness, we use fourth-order numerical methods on Equation (10) and the model from DDIM. In Table 2, we find that the performance of FON is limited when the number of steps is small. By contrast, our new methods, including S-PNDM and F-PNDM, can improve the generated results regardless of whether the number of steps used is large or small. According to Cifar10 / CelebA (linear), F-PNDM can achieve lower FID than 1000 steps DDIM using only 50 steps, making diffusion models 20x faster without losing quality.
We draw a line chart of computation cost with FID according to the results of Cifar10 (linear) above in Figure 3. Because F-PNDM uses the pseudo Runge-Kutta method to generate the first three steps, it is slower than other methods at the first several steps. Therefore, S-PNDM can achieve the best FID initially, then F-PNDM becomes the best and the acceleration is significant.
Quality When the number of steps is relatively big, the results of FON become more and more similar to that of pseudo numerical methods. This is because all the methods are solving Equation (10), and their convergent results should be the same. However, pseudo numerical methods still work better using a large number of steps empirically. F-PNDM can improve the best FID around 0.4 using pretrained models and achieves a new SOTA FID score of 2.71
on CelebA, which shows that our work can not only accelerate diffusion models but also improve the sample quality topline. We also notice that the FID results of F-PNDM converge after more than 250 steps. The FID results will fluctuate around a value then. This phenomenon is more pronounced when we test our methods on LSUN (see Table 5, 6).
According to Cifar10 (cosine), the cosine variance schedule can lower FID using a relatively large number of steps. More analyses about variance schedule can be found in Appendix A.7. What’s more, we test our methods on other datasets and provide our FID results in Appendix A.9 and image results in Appendix A.10. We can draw similar conclusions on our methods’ acceleration and sampling quality, regardless of the datasets and the size of the images.

heading: 4.3 SAMPLE ON MANIFOLDS
text: Here, we design visualization experiments to show the effort of our new methods and support our analyses. Because it is hard to visualize high-dimensional data, we use the change of a global characteristic norm and a local characteristic pixel to show the change of the data under different steps. For pixel, we randomly choose two positions p1, p2. Then for a series of images xT , xT−k, · · · , x0 derived from the reverse process, we denote ykt as the value of xt at position p
k. Then we draw a polyline (y1t , y 2 t )t=T,··· in R2. For norm, we first count the distribution of the norm of the training datasets under different steps and use this to make a heat map as the background. After that, we draw the norm of our generated results using different methods and steps above this heat map.
In Figure 4, we can see that the FON may run far away from the high-density area of the data, which explains why FON may introduce noticeable noise. However, PNDM can avoid this problem and appropriately fit the target result. More visualization results supporting our analysis can be found in Appendix A.11. What’s more, we design a toy example to test our new methods without the influence of neural networks and get similar conclusions as to the real cases above. We put the detailed results in Appendix A.8.

heading: 5 DISCUSSION
text: In this paper, we provided DNPMs, a new numerical method suitable for solving the corresponding ODEs of DDPMs. DNPMs can generate high-quality images using fewer steps without loss of quality successfully. Based on the idea of this work, further improvement can be explored in our future works: 1) find a better variance schedule for DNPMs: although we tested DNPMs on linear variance schedule and cosine variance schedule in this work, there might be another variance schedule more suitable for our proposed numerical methods. 2) Find higher-order convergent pseudo numerical methods: we analyzed the convergence order of S/F-DNPMs, which are both second-order convergent. However, F-DNPMs achieve better FID than S-DNPMs in most cases. We think this is because the result between our transfer part and target ODE has a higher-order error, which limits the convergence order of F-DNPMs. This error from the change of the transfer part is theoretical but does not influence the quality of images according to the property of Equation (11). However, making the transfer part higher-order convergent and finding the influence of such change is still exciting and needs more research. 3) Extend PNDMs to more general applications: when we proved the convergence order of DNPMs, we found that other kinds of transfer parts can keep the convergence order unchanged too, which means that pseudo numerical methods can be used on more general applications, like certain neural ODEs (Chen et al., 2019; Dupont et al., 2019).

heading: A APPENDIX

heading: A.1 RELATED WORK
text: DDPMs have been well developed in the last few years. Some works concentrate on improving the quality and the speed of DDPMs and making DDPMs more practical. Song et al. (2020a) introduce new inference equations to accelerate DDPMs. Nichol & Dhariwal (2021), Watson et al. (2021) and Kong & Ping (2021) choose to find better variance schedules to improve the images quality. Vahdat et al. (2021) combine the advantages of DDPMs and Variational Autoencoders and get better results. Furthermore, Kim et al. (2021) try to solve an existing bottleneck that the inference equations of DDPMs are unbounded in some situations.
Song et al. (2020b) find the similarity between DDPMs and noise conditional score networks (NCSNs (Song & Ermon, 2020b)), which is that they both use a process similar to Langevin dynamics to produce samples. Therefore, some works (Song & Ermon, 2020a; Kim et al., 2021) that can improve the results of NCSNs also can be used in DDPMs. Additional, Song et al. (2020b) combine DDPMs and NCSNs under the framework of neural differential equations (Chen et al., 2019; Dupont et al., 2019). Therefore, numerical methods widely used in neural differential equations can also be applied to accelerate DDPMs. Our work successfully combines the advantages of Song et al. (2020a) and Song et al. (2020b). We use a transfer part from DDIMs and use different gradient parts from different numerical methods. Although we and Song et al. (2020b) solve certain differential
equations derived from DDPMs, we use different target different equations and different numerical methods, which get better results.
The application of DDPMs is not limited to unconditional image generation. Some works apply DDPMs to various types of data successfully, including text-to-speech (Chen et al., 2020; Lam et al., 2021), singing voice (Liu et al., 2021), 3D Point Cloud (Luo & Hu, 2021), text generation (Austin et al., 2021). Additionally, DDPMs can also be used to generate conditional samples, too (Jeong et al., 2021; Choi et al., 2021).

heading: A.2 CONVERGENT ORDER OF METHOD
text: We use the forward Euler method and linear multi-step method to show what is the order of a method. Assume that xt is precise and compute the error at xt+δ . For forward Euler method, we have:
et,δ = x(t+ δ)− xt,δ
= ( x(t) + δf(x(t), t) + δ2
2 f ′′(c)
) − (x(t) + δf(xt, t))
= δ2 2 f ′′(c) ≤ δ 2 2 M
(15)
Here, we assume that f ′′ is continuous, so f ′′ is bounded in a close area. x(t + δ) is the precise result and xt,δ is the numerical result from t to t+ δ.
For linear multi-step method, we have: et,δ =x(t+ δ)− xt,δ
= ( x(t) + δ
1! f(x(t), t) +
δ2
2! f ′(x(t), t) + · · ·+ δ
4
4! f (3)(x(t), t) +O
( δ5 ))
− ( x(t) + b1δf(x(t), t) +
4∑ s=2 bsδ
( 4∑
k=0
(−(s− 1)δ)k
(k)! f (k−1)(x(t), t) +O
( δ5 )))
=O(δ5)
(16)
Here, {bs} satisfy ∑4 s=1 bs = 1 and the following equations for j ∈ {1, · · · , 3}:
(−1)jb2 + (−2)jb3 + (−3)jb4 = 1
j + 1 . (17)
We call the error at xt,δ local error, and the error at xt+Mδ (M is big enough but finite) global error. Assume the local error of our method has order k+1 and the target ODE satisfies Lipschitz condition, then:
x(t+Mδ)− xt,Mδ ≤|x(t+Mδ)− xt+(M−1)δ,δ|+ |xt+(M−1)δ,δ − xt,Mδ| ≤et+(M−1)δ,δ + eLδ|x(t+ (M − 1)δ)− xt,(M−1)δ| ≤et+(M−1)δ,δ + eLδet+(M−2)δ,δ + · · ·
≤Cδk+1 ( 1 + eLh + · · ·+ e(i−1)Lh ) =Chk+1
eiLδ − 1 eLδ − 1 ≤ Cδk+1 e iLδ − 1 Lδ
=O(δk).
(18)
Therefore, the global error will be one order lower than the local error. From Equation (15), we can see the forward Euler method has local error O(δ2) and global error O(δ), so we call it the first-order numerical method. And the linear multi-step method has local error O(δ5) and global error O(δ4), and we call it the fourth-order numerical method.
In addition, assuming that a numerical method has kth-order global error, we can compute the convergence speed of this numerical method:
lim δ→0 x2δt+T − x(t+ T ) xδt+T − x(t+ T ) = (2δ)k δk = 2k. (19)
Here, xδt+T is the result at t+T and move δ every step. This shows that the fourth-order method can converge to the exact solution faster than the first-order method when δ → 0, which means that we can use a bigger iteration interval δ to achieve similar global error and a bigger iteration interval means that we can iterate fewer times to get results with high quality.

heading: A.3 PSEUDO SECOND-ORDER METHOD
text: We introduce two second-order numerical methods. First is improved Euler method satisfying: k1 = f(xt, t) k2 = f(xt + δk1, t+ δ)
xt+δ = xt + δ 2 (k1 + k2)
(20)
Second is another linear multi-step method called second-order linear multi-step method satisfying:
xt+δ = xt + δ
2 (3ft − ft−δ) (21)
And the corresponding pseudo improved Euler methods satisfying: e1t = ϵθ(xt, t) x1t = ϕ(xt, e 1 t , t, t+ δ) e2t = ϵθ(x 1 t , t+ δ) e′t = 1 2 (e 1 t + e 2 t )
xt+δ = ϕ(xt, e ′ t, t, t+ δ)
(22)
Pseudo second-order linear multi-step method satisfying: et = ϵθ(xt, t) e′t = 1 2 (3et − et−δ)
xt+δ = ϕ(xt, e ′ t, t, t+ δ)
(23)
Similar to what we do to get F-PNDMs, We combine them to get S-PNDMs. Abbreviate Equation (22) and (23) as
xt+δ, et = PIE(xt, {ep}p<t, t, t+ δ), xt+δ, e 1 t = PLMS ′(xt, t, t+ δ).
Algorithm 3 S-PNDMs 1: xT ∼ N (0, I) 2: for t = T − 1 do 3: xt, et = PIE(xt+1, t+ 1, t) 4: end for 5: for t = T − 2, · · · , 1, 0 do 6: xt, et = PLMS′(xt+1, {ep}p>t, t+1, t) 7: end for 8: return x0

heading: A.4 THE EXISTENCE OF A DERIVATIVE
text: Because ᾱt is usually obtained by multiplying a linear variance schedule βt. So we have
ᾱt = e at2+bt+c, (24)
and ᾱ0 = 1, so c = 0. Now we have
lim δ→0 xt−δ − xt δ
= lim δ→0 ᾱt−δ − ᾱt δ
( xt√
ᾱt( √ ᾱt−δ + √ ᾱt) − ϵθ(xt, t)√ ᾱt( √ (1− ᾱt−δ)ᾱt + √ (1− ᾱt)ᾱt−δ)
)
= lim δ→0 ᾱt−δ − ᾱt δ ( xt 2ᾱt − ϵθ(xt, t) 2 √ 1− ᾱtᾱt ) = (eat 2+bt)′ ( xt 2ᾱt − ϵθ(xt, t) 2 √ 1− ᾱtᾱt ) =(2at+ b)ᾱt ( xt 2ᾱt − ϵθ(xt, t) 2 √ 1− ᾱtᾱt ) = 1 2 (2at+ b)(xt − ϵθ(xt, t)√ 1− ᾱt ).
(25)
To make limδ→0 xt−δ−xt δ |t=0 is well-defined, b must equal to zero, or (2at+ b) ϵθ(xt,t)√ 1−ᾱt will tend to infinity. This is a strong condition that most variance schedules do not satisfy. In practice, DDPMs can choose the variance schedule very freely. This means that treating DDPMs as ODEs directly is not proper and has theoretical weakness.
A.5 RELATIONSHIP BETWEEN t, ϵθ AND xt
Relationship between t and ϵθ. In Figure 5, we can see that the denoising process tends to converge, whether in the ϵθ domain or the sample/image domain when the step-index tends to zero. Therefore, we can say that the noise becomes more and more precise when step, namely t, tends to zero.
Relationship between ϵθ and xt To prove Property 3.1, assume that xt = √ ᾱtx0 + √ 1− ᾱtϵ, NN is the neural network and ϵθ = NN(xt, t). Because we assume that the gradient part is precise, then we have ϵθ = ϵ. Then for all t′ ≤ t, we have:
xt′ = √ ᾱt′
( xt −
√ 1− ᾱtϵθ√ ᾱt
) + √ 1− ᾱt′ϵθ
= √ ᾱt′
(√ ᾱtx0 + √ 1− ᾱtϵ−
√ 1− ᾱtϵ√
ᾱt
) + √ 1− ᾱt′ϵθ(xt, t)
= √ ᾱt′x0 + √ 1− ᾱt′ϵ.
(26)
Here, we can find that xt′ = √ ᾱt′x0 + √ 1− ᾱt′ϵ is also precise, so Property 3.1 is true.

heading: A.6 ORDER ANALYSIS OF PSEUDO METHOD
text: For the convenience of theoretical analysis, we generalize the problem. Let ϕ(x(t), ϵ, t, δ) = f(x(t), t, δ) + g(t, δ)ϵ(x(t), t) and we have the property f(x(t), t, 0) = g(t, 0) = 0. Then we have:
x(1) = x(0) + ∑ δ→0 (y(t+ δ)− y(t))
= x(0) + ∑ δ→0 (f(x(t), t, δ) + g(t, δ)ϵ(x(t), t))
= x(0) + ∫ 1 0 ( ∂f ∂δ (x(t), t, 0) + ∂g ∂δ (t, 0)ϵ(x(t), t) ) .
(27)
Now, Equation (10) becomes a special case of this more general version and, in this special case, we have:
f(x(t), t, δ) = (√ α(t+ δ)√
α − 1
) x(t)
g(t, δ) = √ 1− α(t+ δ)− √ (1− α(t))α(t+ δ)√
α(t)
(28)
Now, we compute the local error of S-PNDMs. We first compute the theoretical and numerical results of different numerical methods. We have:
x(t+ δ)
=x(t) + δ
( ∂f
∂δ (x(t), t, 0) +
∂g ∂δ (t, 0)ϵ(x(t), t)
) +
δ2
2
( ∂f
∂δ (x(t), t, 0) +
∂g ∂δ (t, 0)ϵ(x(t), t)
)′ +O(δ3)
=x(t) + δ
( ∂f
∂δ (x(t), t, 0) +
∂g ∂δ (t, 0)ϵ(x(t), t)
) +O(δ3)+
δ2
2
( ∂2f
∂δ∂t (x(t), t, 0) +
∂2f
∂δ∂x (x(t), t, 0)
( ∂f
∂δ (x(t), t, 0) +
∂g ∂δ (t, 0)ϵ(x(t), t)
)) +
δ2
2
( ∂2g
∂δ∂t (t, 0)ϵ(x(t), t) +
∂g ∂δ (t, 0)ϵ′(x(t), t)
)
(29)
and
xS-PNDM(t+ δ)
=x(t) + f(x(t), t, δ) + g(t, δ) 1
2 (ϵ(x(t), t) + ϵ(x(t) + ϕ(..., t, δ), t+ δ))
=x(t) + δ ∂f
∂δ (x(t), t, 0) +
δ2
2
∂2f
∂δ2 (x(t), t, 0) +O(δ3)+(
δ ∂g
∂δ (t, 0) +
δ2
2
∂2g ∂δ2 (t, 0)
) 1
2 (ϵ(x(t), t) + ϵ(x(t+ δ) +O(δ2), t+ δ))
=x(t) + δ ∂f
∂δ (x(t), t, 0) +
δ2
2
∂2f
∂δ2 (x(t), t, 0) +O(δ3)+(
δ ∂g
∂δ (t, 0) +
δ2
2
∂2g ∂δ2 (t, 0)
) 1
2 (ϵ(x(t), t) + ϵ(x(t+ δ), t+ δ))
=x(t) + δ ∂f
∂δ (x(t), t, 0) +
δ2
2
∂2f
∂δ2 (x(t), t, 0) +O(δ3)+(
δ ∂g
∂δ (t, 0) +
δ2
2
∂2g ∂δ2 (t, 0)
) 1
2 (ϵ(x(t), t) + ϵ(x(t), t) + δϵ(x(t), t)′)
=x(t) + δ ∂f
∂δ (x(t), t, 0) +
δ2
2
∂2f ∂δ2 (x(t), t, 0) +O(δ3)+
δ ∂g
∂δ (t, 0)
( ϵ(x(t), t) + 1
2 δϵ(x(t), t)′
) + δ2
2
∂2g ∂δ2 (t, 0)ϵ(x(t), t)
=x(t) + δ
( ∂f
∂δ (x(t), t, 0) +
∂g ∂δ (t, 0)ϵ(x(t), t)
) +O(δ3)+
δ2
2
( ∂2f
∂δ2 (x(t), t, 0) +
∂g ∂δ (t, 0)ϵ(x(t), t)′ +
∂2g ∂δ2 (t, 0)ϵ(x(t), t)
) .
(30)
Then we compute the difference between the theoretical and numerical results. We have:
x(t+ δ)− xS-PNDM(x+ δ)
= δ2
2
( ( ∂2f
∂δ∂t − ∂
2f
∂δ2 )(x(t), t, 0) +
∂2f
∂δ∂x (x(t), t, 0)
∂f ∂δ (x(t), t, 0)
) +
δ2
2
( ( ∂2g
∂δ∂t − ∂
2g
∂δ2 )(t, 0) +
∂2f
∂δ∂x (x(t), t, 0)
∂g ∂δ (t, 0)
) ϵ(x(t), t) +O(δ3)
(31)
In this special case, we compute the derivatives of some items needed in Equation (31). We have:
∂g ∂δ (t, δ) = ∂ ∂δ
(√ 1− α(t+ δ)− √ (1− α(t))α(t+ δ)√
α(t)
)
= −α′(t+ δ) 2 √ 1− α((t+ δ)) − √ 1− α(t)α′(t+ δ) 2 √ α(t)α(t+ δ)
(32)
∂f ∂δ (x(t), t, δ) = ∂ ∂δ
((√ α(t+ δ)√
α − 1
) x(t) )
= α′(t+ δ) 2 √ α(t)α(t+ δ) x(t)
(33)
∂2f ∂δ2 (x(t), t, δ)|δ=0 =
α′′(t+ δ) 2 √ α(t)α(t+ δ) x(t) + −α′(t+ δ)2 4 √ α(t)α(t+ δ)3 x(t)|δ=0 (34)
∂2f ∂δ∂t (x(t), t, 0) = α′′(t) 2α(t) x(t)− α ′(t)2 2α(t)2 x(t) (35)
∂2f
∂δ∂x (x(t), t, 0) =
α′(t) 2α(t) (36)
∂2g ∂δ2 (t, δ)|δ=0 = −α′′(t+ δ) 2 √ 1− α((t+ δ)) − √ 1− α(t)α′′(t+ δ) 2 √ α(t)α(t+ δ) +
−α′(t+ δ)2 4 √ 1− α(t+ δ) 3 +
√ 1− α(t)α′(t+ δ)2
4 √ α(t)α(t+ δ)3 |δ=0
= −α′′(t) 2 √ 1− α((t)) − √ 1− α(t)α′′(t) 2 √ α(t)α(t) +
−α′(t)2 4 √ 1− α(t) 3 +
√ 1− α(t)α′(t)2
4 √ α(t)α(t)3
(37)
∂2g
∂δ∂t (t, δ)|δ=0 =
∂
∂t
( −α′(t)
2 √ 1− α((t)) − √ 1− α(t)α′(t) 2 √ α(t)α(t) )
= −α′′(t) 2 √ 1− α((t)) − √ 1− α(t)α′′(t) 2 √ α(t)α(t) +
−α′(t)2 4 √ 1− α(t) 3 +
√ 1− α(t)α′(t)2
2α(t)2 +
α′(t)2 4 √ α2(t)(1− α(t))
(38)
Now we can compute the final result of Equation (31). We split it into three parts and the values of the first two terms. We have:
( ∂2f ∂δ∂t − ∂ 2f ∂δ2 )(x(t), t, 0) + ∂2f ∂δ∂x (x(t), t, 0) ∂f ∂δ (x(t), t, 0)
=
( α′′(t)
2α(t) x(t)− α
′(t)2
2α(t)2 x(t)
) − ( α′′(t)
2α(t) x(t) +
−α′(t)2
4α(t)2 x(t)
) + 1
2α(t)
α′(t) 2α(t) x(t)
=0
(39)
and
( ∂2g ∂δ∂t − ∂ 2g ∂δ2 )(t, 0) + ∂2f ∂δ∂x (x(t), t, 0) ∂g ∂δ (t, 0)
=
√ 1− α(t)α′(t)2
4α(t)2 +
α′(t)2 4α(t) √ 1− α(t) + α′(t) 2α(t)
( −α′(t)
2 √ 1− α((t))
− √ 1− α(t)α′(t)
2α(t)
)
= α′(t)2 4α(t)2 √ 1− α(t) + α′(t) 2α(t)
( −α′(t)
2 √ 1− α(t)α(t) ) =0
(40)
Finally, we get the final result of Equation (31):
x(t+ δ)− xS-PNDM(x+ δ) = O(δ3) (41) And the computation of the convergence order of F-PNDMs is similar, and we ignore it here. Therefore, Property 3.2 is true.
A.7 VARIANCE SCHEDULE
According to Cifar10 (cosine) in Table 2, PNDMs can be used on both linear variance schedule and cosine variance schedule. However, we also notice cosine variance schedule can make FID lower when we use relatively big generation steps, but the effort is limited when the number of steps is small. F-PNDM uses information from four consecutive steps, so the smoothness of the schedule is more important for F-PNDM than DDIM. According to this experiment, our work can be used with works that pay attention to variance schedules to improve the acceleration effect further. However, a variance schedule that fits pseudo numerical methods better remains to be found in further work.

heading: A.8 TOY EXAMPLE
text: Here, we design a toy example to test our new methods without the influence of neural networks. We randomly generate the initial input x1 = (m1,m2), mi ∼ U(0, 1) and use a simple analytic equations ϵθ(x) = (sin x[0], cosx[1]) to replace the neural networks in real cases. Let ϕ in Equation (11) is unchanged and ᾱt = α(t) = 1− t, then we get:
ϕ(xt, ϵθ(xt), t, t− δ)
= √ ᾱt−δ√ ᾱt xt − (ᾱt−δ − ᾱt)√ ᾱt( √ (1− ᾱt−δ)ᾱt + √ (1− ᾱt)ᾱt−δ) ϵt
= √ 1− (t− δ)√
1− t xt − δ √ 1− t (√ (t− δ)(1− t) + √ t(1− (t− δ))
)ϵθ(xt) (42)
Here, we use three different numerical methods to generate x0.
For DDIM, we have: xt−δ = xt + ϕ(xt, ϵθ(xt), t, t− δ) (43)
For FON, according to Equation (10), we have:
e′t = ᾱ ′(t)
( xt
2ᾱ(t) − ϵθ(xt) 2ᾱ(t) √ 1− ᾱ(t)
)
= − (
xt 2(1− t) − ϵθ(xt) 2(1− t) √ t ) xt−δ = xt + δ
24 (55e′t − 59e′t+δ + 37e′t+2δ − 9e′t+3δ)
(44)
For F-PNDM, we have:
e′ = 1
24 (55ϵθ(xt)− 59ϵθ(xt+δ) + 37ϵθ(xt+2δ)− 9ϵθ(xt+3δ))
xt−δ = xt + ϕ(xt, e ′, t, t− δ)
(45)
Then we draw the corresponding generation curves in Figure 6. We find that the result is similar to the real cases. The main difference here is that FON can correct its results while the real case cannot. The reason is that the gradient is well-defined everywhere, while in real cases, the gradient is meaningful on the high-density region of the data xt of DDPMs.

heading: A.9 MORE FID RESULTS
text: Here, We provide our more detailed FID results on Cifar10, CelebA, LSUN-church and LSUNbedroom in Table 3, 4, 5, 6.

heading: A.10 MORE IMAGE RESULTS
text: Here, we show more generated images on Cifar10, CelebA, LSUN-church and LSUN-bedroom in Figure 7, 9, 8, 10, 11, 12.

heading: A.11 MORE VISUALIZATION RESULTS
text: Here, we put more visualization results similar to Figure 4 in Figure 13.

heading: A.12 FID RESULT WITH STANDARD DEVIATION
text: Here, we report the mean and standard deviation of FID results, tested over four sampling runs.

Review summaries:
The paper proposes a new efficient method for denoising diffusion probabilistic models (DDPM) (generative models that optimize for the closest solution on a manifold) based on the observation that this can be seen a solving a set of differential equations on a manifold. This allows efficient pseudo numerical methods to be applied here which have many advantageous properties over classical optimization methods, including less optimization steps and guaranteed manifold solutions due to separating the gradient part from the transfer part in the optimization. Results are shown on four datasets comparing to two reasonable but not sota baselines. 

Highlighting the high computational complexity for sampling from Denoising Diffusion Probabilistic Models (DDPMs) (e.g. wrt GANs), authors build on the connection between diffusion processes and ODEs to propose efficient (pseudo-)numerical methods so as to sample data from the data manifold. The main idea is to combine the discrete update proposed in DDIMs, with a fourth-order gradient estimation given by the Runge-Kutta or linear multi-step methods. The motivation being that such gradient estimator should yield trajectories that stay closer to the data manifold. They empirically assess their method(s) on CIFAR10 and CelebA in terms of sample quality (measured by FID), and show that they get a ~x20 speed-up wrt DDIMs or a significant improvement in FID with the same number of steps.