PAPER KEY: 57UHJ2KG
TITLE: DPM-Solver: A Fast ODE Solver for Diffusion Probabilistic Model Sampling in Around 10 Steps
AUTHORS: Lu, Cheng; Zhou, Yuhao; Bao, Fan; Chen, Jianfei; Li, Chongxuan; Zhu, Jun

DPM-Solver: A Fast ODE Solver for Diffusion Probabilistic Model Sampling in Around 10 Steps
Cheng Lu†, Yuhao Zhou†, Fan Bao†, Jianfei Chen†∗, Chongxuan Li‡, Jun Zhu†∗ †Dept. of Comp. Sci. & Tech., Institute for AI, BNRist Center, THBI Lab †Tsinghua-Bosch Joint ML Center, Tsinghua University, Beijing, 100084 China ‡Gaoling School of Artificial Intelligence, Renmin University of China, ‡Beijing Key Laboratory of Big Data Management and Analysis Methods, Beijing, China {lucheng.lc15, yuhaoz.cs}@gmail.com; bf19@mails.tsinghua.edu.cn chongxuanli@ruc.edu.cn; {jianfeic, dcszj}@tsinghua.edu.cn
Abstract
Diffusion probabilistic models (DPMs) are emerging powerful generative models. Despite their high-quality generation performance, DPMs still suffer from their slow sampling as they generally need hundreds or thousands of sequential function evaluations (steps) of large neural networks to draw a sample. Sampling from DPMs can be viewed alternatively as solving the corresponding diffusion ordinary differential equations (ODEs). In this work, we propose an exact formulation of the solution of diffusion ODEs. The formulation analytically computes the linear part of the solution, rather than leaving all terms to black-box ODE solvers as adopted in previous works. By applying change-of-variable, the solution can be equivalently simplified to an exponentially weighted integral of the neural network. Based on our formulation, we propose DPM-Solver, a fast dedicated high-order solver for diffusion ODEs with the convergence order guarantee. DPM-Solver is suitable for both discrete-time and continuous-time DPMs without any further training. Experimental results show that DPM-Solver can generate high-quality samples in only 10 to 20 function evaluations on various datasets. We achieve 4.70 FID in 10 function evaluations and 2.87 FID in 20 function evaluations on the CIFAR10 dataset, and a 4 ∼ 16× speedup compared with previous state-of-the-art training-free samplers on various datasets.2
1 Introduction
Diffusion probabilistic models (DPMs) [1–3] are emerging powerful generative models with promising performance on many tasks, such as image generation [4, 5], video generation [6], text-to-image generation [7], speech synthesis [8, 9] and lossless compression [10]. DPMs are defined by discretetime random processes [1, 2] or continuous-time stochastic differential equations (SDEs) [3], which learn to gradually remove the noise added to the data points. Compared with the widely-used generative adversarial networks (GANs) [11] and variational auto-encoders (VAEs) [12], DPMs can not only compute exact likelihood [3], but also achieve even better sample quality for image generation [4]. However, to obtain high-quality samples, DPMs usually need hundreds or thousands of sequential steps of large neural network evaluations, thereby resulting in a much slower sampling speed than the single-step GANs or VAEs. Such inefficiency is becoming a critical bottleneck for the adoption of DPMs in downstream tasks, leading to an urgent request to design fast samplers for DPMs.
∗Corresponding Author.
2Code is available at https://github.com/LuChengTHU/dpm-solver
36th Conference on Neural Information Processing Systems (NeurIPS 2022).
arXiv:2206.00927v3 [cs.LG] 13 Oct 2022


NFE = 10 NFE = 15 NFE = 20 NFE = 100 NFE = 10
(a) DDIM [19] (b) DPM-Solver (ours)
Figure 1: Samples by DDIM [19] with 10, 15, 20, 100 number of function evaluations (NFE), and DPM-Solver (ours) with only 10 NFE, using the pre-trained DPMs on ImageNet 256×256 with classifier guidance [4].
Existing fast samplers for DPMs can be divided into two categories. The first category includes knowledge distillation [13, 14] and noise level or sample trajectory learning [15–18]. Such methods require a possibly expensive training stage before they can be used for efficient sampling. Furthermore, their applicability and flexibility might be limited. It might require nontrivial effort to adapt the method to different models, datasets, and number of sampling steps. The second category consists of training-free [19–21] samplers, which are suitable for all pre-trained DPMs in a simple plug-andplay manner. Training-free samplers include adopting implicit [19] or analytical [21] generation process, advanced differential equation (DE) solvers [3, 20, 22–24] and dynamic programming [18]. However, these methods still require ∼ 50 function evaluations [21] to generate high-quality samples (comparable to those generated by plain samplers in about 1000 function evaluations), thereby are still time-consuming.
In this work, we bring the efficiency of training-free samplers to a new level to produce high-quality samples in the “few-step sampling” regime, where the sampling can be done within around 10 steps of sequential function evaluations. We tackle the alternative problem of sampling from DPMs as solving the corresponding diffusion ordinary differential equations (ODEs) of DPMs, and carefully examine the structure of diffusion ODEs. Diffusion ODEs have a semi-linear structure — they consist of a linear function of the data variable and a nonlinear function parameterized by neural networks. Such structure is omitted in previous training-free samplers [3, 20], which directly use black-box DE solvers. To utilize the semi-linear structure, we derive an exact formulation of the solutions of diffusion ODEs by analytically computing the linear part of the solutions, avoiding the corresponding discretization error. Furthermore, by applying change-of-variable, the solutions can be equivalently simplified to an exponentially weighted integral of the neural network. Such integral is very special and can be efficiently approximated by the numerical methods for exponential integrators [25].
Based on our formulation of solutions, we propose DPM-Solver, a fast dedicated solver for diffusion ODEs by approximating the above integral. Specifically, we propose first-order, second-order and third-order versions of DPM-Solver with convergence order guarantees. We further propose an adaptive step size schedule for DPM-Solver. In general, DPM-Solver is applicable to both continuoustime and discrete-time DPMs, and also conditional sampling with classifier guidance [4]. Fig. 1 demonstrates the speedup performance of a Denoising Diffusion Implicit Models (DDIM) [19] baseline and DPM-Solver, which shows that DPM-Solver can generate high-quality samples with as few as 10 function evaluations and is much faster than DDIM on the ImageNet 256x256 dataset [26]. Our additional experimental results show that DPM-Solver can greatly improve the sampling speed of both discrete-time and continuous-time DPMs, and it can achieve excellent sample quality in around 10 function evaluations, which is much faster than all previous training-free samplers of DPMs.
2 Diffusion Probabilistic Models
We review diffusion probabilistic models and their associated differential equations in this section.
2


2.1 Forward Process and Diffusion SDEs
Assume that we have a D-dimensional random variable x0 ∈ RD with an unknown distribution q0(x0). Diffusion Probabilistic Models (DPMs) [1–3, 10] define a forward process {xt}t∈[0,T ] with T > 0 starting with x0, such that for any t ∈ [0, T ], the distribution of xt conditioned on x0 satisfies
q0t(xt|x0) = N (xt|α(t)x0, σ2(t)I), (2.1)
where α(t), σ(t) ∈ R+ are differentiable functions of t with bounded derivatives, and we denote them as αt, σt for simplicity. The choice for αt and σt is referred to as the noise schedule of a DPM. Let qt(xt) denote the marginal distribution of xt, DPMs choose noise schedules to ensure
that qT (xT ) ≈ N (xT |0, σ ̃2I) for some σ ̃ > 0, and the signal-to-noise-ratio (SNR) αt2/σt2 is strictly decreasing w.r.t. t [10]. Moreover, Kingma et al. [10] prove that the following stochastic differential equation (SDE) has the same transition distribution q0t(xt|x0) as in Eq. (2.1) for any t ∈ [0, T ]:
dxt = f (t)xtdt + g(t)dwt, x0 ∼ q0(x0), (2.2)
where wt ∈ RD is the standard Wiener process, and
f (t) = d log αt
dt , g2(t) = dσt2
dt − 2 d log αt
dt σ2
t . (2.3)
Under some regularity conditions, Song et al. [3] show that the forward process in Eq. (2.2) has an equivalent reverse process from time T to 0, starting with the marginal distribution qT (xT ):
dxt = [f (t)xt − g2(t)∇x log qt(xt)]dt + g(t)dw ̄t, xT ∼ qT (xT ), (2.4)
where w ̄t is a standard Wiener process in the reverse time. The only unknown term in Eq. (2.4) is the score function ∇x log qt(xt) at each time t. In practice, DPMs use a neural network θ(xt, t) parameterized by θ to estimate the scaled score function: −σt∇x log qt(xt). The parameter θ is optimized by minimizing the following objective [2, 3]:
L(θ; ω(t)) := 1
2
∫T
0
ω(t)Eqt(xt)
[
‖θ(xt, t) + σt∇x log qt(xt)‖2
2
]
dt
=1
2
∫T
0
ω(t)Eq0(x0)Eq()
[
‖θ(xt, t) − ‖2
2
]
dt + C,
where ω(t) is a weighting function,  ∼ q() = N (|0, I), xt = αtx0 + σt, and C is a constant independent of θ. As θ(xt, t) can also be regarded as predicting the Gaussian noise added to xt, it is usually called the noise prediction model. Since the ground truth of θ(xt, t) is −σt∇x log qt(xt), DPMs replace the score function in Eq. (2.4) by −θ(xt, t)/σt and define a parameterized reverse
process (diffusion SDE) from time T to 0, starting with xT ∼ N (0, σ ̃2I):
dxt =
[
f (t)xt + g2(t)
σt
θ(xt, t)
]
dt + g(t)dw ̄t, xT ∼ N (0, σ ̃2I). (2.5)
Samples can be generated from DPMs by solving the diffusion SDE in Eq. (2.5) with numerical solvers, which discretize the SDE from T to 0. Song et al. [3] proved that the traditional ancestral sampling method for DPMs [2] can be viewed as a first-order SDE solver for Eq. (2.5). However, these first-order methods usually need hundreds of or thousands of function evaluations to converge [3], leading to extremely slow sampling speed.
2.2 Diffusion (Probability Flow) ODEs
When discretizing SDEs, the step size is limited by the randomness of the Wiener process [27, Chap. 11]. A large step size (small number of steps) often causes non-convergence, especially in high dimensional spaces. For faster sampling, one can consider the associated probability flow ODE [3], which has the same marginal distribution at each time t as that of the SDE. Specifically, for DPMs, Song et al. [3] proved that the probability flow ODE of Eq. (2.4) is
dxt
dt = f (t)xt − 1
2 g2(t)∇x log qt(xt), xT ∼ qT (xT ), (2.6)
3


where the marginal distribution of xt is also qt(xt). By replacing the score function with the noise prediction model, Song et al. [3] defined the following parameterized ODE (diffusion ODE):
dxt
dt = hθ(xt, t) := f (t)xt + g2(t)
2σt
θ(xt, t), xT ∼ N (0, σ ̃2I). (2.7)
Samples can be drawn by solving the ODE from T to 0. Comparing with SDEs, ODEs can be solved with larger step sizes as they have no randomness. Furthermore, we can take advantage of efficient numerical ODE solvers to accelerate the sampling. Song et al. [3] used the RK45 ODE solver [28] for the diffusion ODEs, which generates samples in ∼ 60 function evaluations to reach comparable quality with a 1000-step SDE solver for Eq. (2.5) on the CIFAR-10 dataset [29]. However, existing general-purpose ODE solvers still cannot generate satisfactory samples in the few-step (∼ 10 steps) sampling regime. To the best of our knowledge, there is still a lack of training-free samplers for DPMs in the few-step sampling regime, and the sampling speed of DPMs is still a critical issue.
3 Customized Fast Solvers for Diffusion ODEs
As highlighted in Sec. 2.2, discretizing SDEs is generally difficult in high dimensions [27, Chap. 11] and it is hard to converge within few steps. In contrast, ODEs are easier to solve, yielding a potential for fast samplers. However, as mentioned in Sec. 2.2, the general black-box ODE solver used in previous work [3] empirically fails to converge in few steps. This motivates us to design a dedicated solver for diffusion ODEs to enable fast and high-quality few-step sampling. We start with a detailed investigation of the specific structure of diffusion ODEs.
3.1 Simplified Formulation of Exact Solutions of Diffusion ODEs
The key insight of this work is that given an initial value xs at time s > 0, the solution xt at each time t < s of diffusion ODEs in Eq. (2.7) can be simplified into a very special exact formulation which can be efficiently approximated.
Our first key observation is that a part of the solution xt can be exactly computed by considering the particular structure of diffusion ODEs. The r.h.s. of diffusion ODEs in Eq. (2.7) consists of two parts:
the part f (t)xt is a linear function of xt, and the other part g2(t)
2σt θ(xt, t) is generally a nonlinear function of xt because of the neural network θ(xt, t). This type of ODE is referred to as semi-linear ODE. The black-box ODE solvers adopted by previous work [3] are ignorant of this semi-linear structure as they take the whole hθ(xt, t) in Eq. (2.7) as the input, which causes discretization errors of both the linear and nonlinear term. We note that for semi-linear ODEs, the solution at time t can be exactly formulated by the “variation of constants” formula [30]:
xt = e
∫t
s f (τ )dτ xs +
∫t
s
(
e
∫t
τ f (r)dr g2(τ )
2στ
θ(xτ , τ )
)
dτ. (3.1)
This formulation decouples the linear part and the nonlinear part. In contrast to black-box ODE solvers, the linear part is now exactly computed, which eliminates the approximation error of the linear term. However, the integral of the nonlinear part is still complicated because it couples the coefficients about the noise schedule (i.e., f (τ ), g(τ ), στ ) and the complex neural network θ, which is still hard to approximate.
Our second key observation is that the integral of the nonlinear part can be greatly simplified by introducing a special variable. Let λt := log(αt/σt) (one half of the log-SNR), then λt is a strictly decreasing function of t (due to the definition of DPMs as discussed in Sec. 2.1). We can rewrite g(t) in Eq. (2.3) as
g2(t) = dσt2
dt − 2 d log αt
dt σ2
t = 2σ2
t
( d log σt
dt − d log αt
dt
)
= −2σ2
t
dλt
dt . (3.2)
Combining with f (t) = d log αt/dt in Eq. (2.3), we can rewrite Eq. (3.1) as
xt = αt
αs
xs − αt
∫t
s
( dλτ dτ
) στ
ατ
θ(xτ , τ )dτ. (3.3)
As λ(t) = λt is a strictly decreasing function of t, it has an inverse function tλ(·) satisfying t = tλ(λ(t)). We further change the subscripts of x and θ from t to λ and denote xˆλ := xtλ(λ), ˆθ(xˆλ, λ) := θ(xtλ(λ), tλ(λ)). Rewrite Eq. (3.3) by “change-of-variable” for λ, then we have:
4


Proposition 3.1 (Exact solution of diffusion ODEs). Given an initial value xs at time s > 0, the solution xt at time t ∈ [0, s] of diffusion ODEs in Eq. (2.7) is:
xt = αt
αs
xs − αt
∫ λt
λs
e−λˆθ(xˆλ, λ)dλ. (3.4)
We call the integral ∫ e−λˆθ(xˆλ, λ)dλ the exponentially weighted integral of ˆθ, which is very special and highly related to the exponential integrators in the literature of ODE solvers [25]. To the best of our knowledge, such formulation has not been revealed in prior work of diffusion models.
Eq. (3.4) provides a new perspective for approximating the solutions of diffusion ODEs. Specifically, given xs at time s, According to Eq. (3.4), approximating the solution at time t is equivalent to directly approximating the exponentially weighted integral of ˆθ from λs to λt, which avoids the error of the linear terms and is well-studied in the literature of exponential integrators [25, 31]. Based on this insight, we propose fast solvers for diffusion ODEs, as detailed in the following sections.
3.2 High-Order Solvers for Diffusion ODEs
In this section, we propose high-order solvers for diffusion ODEs with convergence order guarantee by leveraging our proposed solution formulation Eq. (3.4). The proposed solvers and analysis are highly motivated by the methods of exponential integrators [25, 31] in the ODE literature.
Specifically, given an initial value xT at time T and M + 1 time steps {ti}M
i=0 decreasing from t0 = T to tM = 0. Let x ̃t0 = xT be the initial value. The proposed solvers use M steps to iteratively
compute a sequence {x ̃ti }M
i=0 to approximate the true solutions at time steps {ti}M
i=0. In particular, the last iterate x ̃tM approximates the true solution at time 0.
In order to reduce the approximation error between x ̃tM and the true solution at time 0, we need to reduce the approximation error for each x ̃ti at every step [30]. Starting with the previous value x ̃ti−1 at time ti−1, according to Eq. (3.4), the exact solution xti−1→ti at time ti is given by
xti−1→ti = αti
αti−1
x ̃ti−1 − αti
∫ λti
λti−1
e−λˆθ(xˆλ, λ)dλ. (3.5)
Therefore, to compute the value x ̃ti for approximating xti−1→ti , we need to approximate the
exponentially weighted integral of ˆθ from λti−1 to λti . Denote hi := λti −λti−1 , and ˆ(n)
θ (xˆλ, λ) :=
dnˆθ (xˆλ,λ)
dλn as the n-th order total derivative of ˆθ(xˆλ, λ) w.r.t. λ. For k ≥ 1, the (k − 1)-th order
Taylor expansion of ˆθ(xˆλ, λ) w.r.t. λ at λti−1 is
ˆθ(xˆλ, λ) =
k−1
∑
n=0
(λ − λti−1 )n
n! ˆ(n)
θ (xˆλti−1 , λti−1 ) + O((λ − λti−1 )k),
Substituting the above Taylor expansion into Eq. (3.5) yields
xti−1→ti = αti
αti−1
x ̃ti−1 − αti
k−1
∑
n=0
ˆ(n)
θ (xˆλti−1 , λti−1 )
∫ λti
λti−1
e−λ (λ − λti−1 )n
n! dλ + O(hk+1
i ), (3.6)
where the integral ∫ e−λ (λ−λti−1 )n
n! dλ can be analytically computed by repeatedly applying n times
of integration-by-parts (see Appendix B.2). Therefore, to approximate xti−1→ti , we only need to
approximate the n-th order total derivatives ˆ(n)
θ (xˆλ, λ) for n ≤ k − 1, which is a well-studied
problem in the ODE literature [31, 32]. By dropping the O(hk+1
i ) error term and approximating the
first (k − 1)-th total derivatives with the “stiff order conditions” [31, 32], we can derive k-th-order ODE solvers for diffusion ODEs. We name such solvers as DPM-Solver overall, and DPM-Solver-k for a specific order k. Here we take k = 1 for demonstration. In this case, Eq. (3.6) becomes
xti−1→ti = αti
αti−1
x ̃ti−1 − αti θ(x ̃ti−1 , ti−1)
∫ λti
λti−1
e−λdλ + O(h2
i)
= αti
αti−1
x ̃ti−1 − σti (ehi − 1)θ(x ̃ti−1 , ti−1) + O(h2
i ).
5


By dropping the high-order error term O(h2
i ), we can obtain an approximation for xti−1→ti . As k = 1 here, we call this solver DPM-Solver-1, and the detailed algorithm is as following.
DPM-Solver-1. Given an initial value xT and M + 1 time steps {ti}M
i=0 decreasing from t0 = T
to tM = 0. Starting with x ̃t0 = xT , the sequence {x ̃ti }M
i=1 is computed iteratively as follows:
x ̃ti = αti
αti−1
x ̃ti−1 − σti (ehi − 1)θ(x ̃ti−1 , ti−1), where hi = λti − λti−1 . (3.7)
For k ≥ 2, approximating the first k terms of the Taylor expansion needs additional intermediate points between t and s [31]. The derivation is more technical so we defer it to Appendix B. Below we propose algorithms for k = 2, 3 and name them as DPM-Solver-2 and DPM-Solver-3, respectively.
Algorithm 1 DPM-Solver-2.
Require: initial value xT , time steps {ti}M
i=0, model θ 1: x ̃t0 ← xT 2: for i ← 1 to M do 3: si ← tλ
( λti−1 +λti 2
)
4: ui ← αsi
αti−1
x ̃ti−1 − σsi
(
e hi
2 −1
)
θ(x ̃ti−1 , ti−1)
5: x ̃ti ← αti
αti−1
x ̃ti−1 − σti
(ehi − 1) θ(ui, si)
6: end for 7: return x ̃tM
Algorithm 2 DPM-Solver-3.
Require: initial value xT , time steps {ti}M
i=0, model θ
1: x ̃t0 ← xT , r1 ← 1
3 , r2 ← 2
3
2: for i ← 1 to M do 3: s2i−1 ← tλ
(λti−1 + r1hi
) , s2i ← tλ
(λti−1 + r2hi
)
4: u2i−1 ← αs2i−1
αti−1
x ̃ti−1 − σs2i−1
(er1hi − 1) θ(x ̃ti−1 , ti−1)
5: D2i−1 ← θ(u2i−1, s2i−1) − θ(x ̃ti−1 , ti−1)
6: u2i ← αs2i
αti−1
x ̃ti−1 − σs2i
(er2hi − 1) θ(x ̃ti−1 , ti−1) − σs2i r2
r1
( er2hi −1
r2hi − 1
)
D2i−1
7: D2i ← θ(u2i, s2i) − θ(x ̃ti−1 , ti−1)
8: x ̃ti ← αti
αti−1
x ̃ti−1 − σti
(ehi − 1) θ(x ̃ti−1 , ti−1) − σti
r2
( ehi −1
h −1
)
D2i
9: end for 10: return x ̃tM
Here, tλ(·) is the inverse function of λ(t), which has an analytical formulation for the practical noise schedule used in [2, 16], as shown in Appendix D. The chosen intermediate points are (si, ui) for DPM-Solver-2 and