DPM-Solver: A Fast ODE Solver for Diffusion Probabilistic Model Sampling in Around 10 Steps

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\sim 16\times$ speedup compared with previous state-of-the-art training-free samplers on various datasets.

1 Introduction

Diffusion probabilistic models (DPMs) [1, 2, 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 discrete-time 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.

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, 16, 17, 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, 20, 21] samplers, which are suitable for all pre-trained DPMs in a simple plug-and-play manner. Training-free samplers include adopting implicit [19] or analytical [21] generation process, advanced differential equation (DE) solvers [3, 20, 22, 23, 24] and dynamic programming [18]. However, these methods still require $\sim$ 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 continuous-time 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.

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$\times$256 with classifier guidance [4].

2 Diffusion Probabilistic Models

We review diffusion probabilistic models and their associated differential equations in this section.

2.1 Forward Process and Diffusion SDEs

Assume that we have a $D$-dimensional random variable ${𝒙}_0\in {ℝ}^D$ with an unknown distribution $q_0({𝒙}_0)$. Diffusion Probabilistic Models (DPMs) [1, 2, 3, 10] define a forward process ${\{{𝒙}_t\}}_{t\in [0,T]}$ with $T>0$ starting with ${𝒙}_0$, such that for any $t\in [0,T]$, the distribution of ${𝒙}_t$ conditioned on ${𝒙}_0$ satisfies

$$q_{0t}({𝒙}_t|{𝒙}_0)=𝒩({𝒙}_t|\alpha (t){𝒙}_0,{\sigma }^2(t)𝑰),$$ (2.1)

where $\alpha (t),\sigma (t)\in {ℝ}^{+}$ are differentiable functions of $t$ with bounded derivatives, and we denote them as ${\alpha }_t,{\sigma }_t$ for simplicity. The choice for ${\alpha }_t$ and ${\sigma }_t$ is referred to as the noise schedule of a DPM. Let $q_t({𝒙}_t)$ denote the marginal distribution of ${𝒙}_t$, DPMs choose noise schedules to ensure that $q_T({𝒙}_T)\approx 𝒩({𝒙}_T|\mathrm{𝟎},{\tilde{\sigma }}^2𝑰)$ for some $\tilde{\sigma }>0$, and the signal-to-noise-ratio (SNR) ${\alpha }_t^2/{\sigma }_t^2$ 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 $q_{0t}({𝒙}_t|{𝒙}_0)$ as in Eq. (2.1) for any $t\in [0,T]$:

$$\mathrm{d}{𝒙}_t=f(t){𝒙}_t\mathrm{d}t+g(t)\mathrm{d}{𝒘}_t,{𝒙}_0\sim q_0({𝒙}_0),$$ (2.2)

where ${𝒘}_t\in {ℝ}^D$ is the standard Wiener process, and

$$f(t)=\frac{d\log {\alpha }_t}{\mathrm{d}t},g^2(t)=\frac{d{\sigma }_t^2}{\mathrm{d}t}-2\frac{d\log {\alpha }_t}{\mathrm{d}t}{\sigma }_t^2.$$ (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 $q_T({𝒙}_T)$:

$$\mathrm{d}{𝒙}_t=[f(t){𝒙}_t-g^2(t){\nabla }_{𝒙}\log q_t({𝒙}_t)]\mathrm{d}t+g(t)d{\overline{𝒘}}_t,{𝒙}_T\sim q_T({𝒙}_T),$$ (2.4)

where ${\overline{𝒘}}_t$ is a standard Wiener process in the reverse time. The only unknown term in Eq. (2.4) is the score function ${\nabla }_{𝒙}\log q_t({𝒙}_t)$ at each time $t$. In practice, DPMs use a neural network ${\mathit{ϵ}}_{\theta }({𝒙}_t,t)$ parameterized by $\theta$ to estimate the scaled score function: $-{\sigma }_t{\nabla }_{𝒙}\log q_t({𝒙}_t)$. The parameter $\theta$ is optimized by minimizing the following objective [2, 3]:

$ℒ(\theta ;\omega (t))$ $≔{\displaystyle \frac{1}{2}}{\displaystyle {\int }_0^T}\omega (t){𝔼}_{q_t({𝒙}_t)}\left[{\Vert {\mathit{ϵ}}_{\theta }({𝒙}_t,t)+{\sigma }_t{\nabla }_{𝒙}\log q_t({𝒙}_t)\Vert }_2^2\right]dt$

$={\displaystyle \frac{1}{2}}{\displaystyle {\int }_0^T}\omega (t){𝔼}_{q_0({𝒙}_0)}{𝔼}_{q(\mathit{ϵ})}\left[{\Vert {\mathit{ϵ}}_{\theta }({𝒙}_t,t)-\mathit{ϵ}\Vert }_2^2\right]dt+C,$

where $\omega (t)$ is a weighting function, $\mathit{ϵ}\sim q(\mathit{ϵ})=𝒩(\mathit{ϵ}|\mathrm{𝟎},𝑰)$, ${𝒙}_t={\alpha }_t{𝒙}_0+{\sigma }_t\mathit{ϵ}$, and $C$ is a constant independent of $\theta$. As ${\mathit{ϵ}}_{\theta }({𝒙}_t,t)$ can also be regarded as predicting the Gaussian noise added to ${𝒙}_t$, it is usually called the noise prediction model. Since the ground truth of ${\mathit{ϵ}}_{\theta }({𝒙}_t,t)$ is $-{\sigma }_t{\nabla }_{𝒙}\log q_t({𝒙}_t)$, DPMs replace the score function in Eq. (2.4) by $-{\mathit{ϵ}}_{\theta }({𝒙}_t,t)/{\sigma }_t$ and define a parameterized reverse process (diffusion SDE) from time $T$ to $0$, starting with ${𝒙}_T\sim 𝒩(\mathrm{𝟎},{\tilde{\sigma }}^2𝑰)$:

$$\mathrm{d}{𝒙}_t=\left[f(t){𝒙}_t+\frac{g^2(t)}{{\sigma }_t}{\mathit{ϵ}}_{\theta }({𝒙}_t,t)\right]\mathrm{d}t+g(t)d{\overline{𝒘}}_t,{𝒙}_T\sim 𝒩(\mathrm{𝟎},{\tilde{\sigma }}^2𝑰).$$ (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

$$\frac{\mathrm{d}{𝒙}_t}{\mathrm{d}t}=f(t){𝒙}_t-\frac{1}{2}{g}^2(t){\nabla }_{𝒙}\log q_t({𝒙}_t),{𝒙}_T\sim q_T({𝒙}_T),$$ (2.6)

where the marginal distribution of ${𝒙}_t$ is also $q_t({𝒙}_t)$. By replacing the score function with the noise prediction model, Song et al. [3] defined the following parameterized ODE (diffusion ODE):

$$\frac{\mathrm{d}{𝒙}_t}{\mathrm{d}t}={𝒉}_{\theta }({𝒙}_t,t)≔f(t){𝒙}_t+\frac{g^2(t)}{2{\sigma }_t}{\mathit{ϵ}}_{\theta }({𝒙}_t,t),{𝒙}_T\sim 𝒩(\mathrm{𝟎},{\tilde{\sigma }}^2𝑰).$$ (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 $\sim$ 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 ($\sim$ 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 ${𝒙}_s$ at time $s>0$, the solution ${𝒙}_t$ at each time 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 ${𝒙}_t$ 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){𝒙}_t$ is a linear function of ${𝒙}_t$, and the other part $\frac{g^2(t)}{2{\sigma }_t}{\mathit{ϵ}}_{\theta }({𝒙}_t,t)$ is generally a nonlinear function of ${𝒙}_t$ because of the neural network ${\mathit{ϵ}}_{\theta }({𝒙}_t,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 ${𝒉}_{\theta }({𝒙}_t,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]:

$${𝒙}_t=e^{{\int }_s^{t}f(\tau )d\tau }{𝒙}_s+{\int }_s^t\left(e^{{\int }_{\tau }^{t}f(r)dr}\frac{g^2(\tau )}{2{\sigma }_{\tau }}{\mathit{ϵ}}_{\theta }({𝒙}_{\tau },\tau )\right)d\tau .$$ (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(\tau ),g(\tau ),{\sigma }_{\tau }$) and the complex neural network ${\mathit{ϵ}}_{\theta }$, 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 ${\lambda }_t≔\log \left({\alpha }_t/{\sigma }_t\right)$ (one half of the log-SNR), then ${\lambda }_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

$g^2(t)={\displaystyle \frac{d{\sigma }_t^2}{\mathrm{d}t}}-2{\displaystyle \frac{d\log {\alpha }_t}{\mathrm{d}t}}{\sigma }_t^2=2{\sigma }_t^2\left({\displaystyle \frac{d\log {\sigma }_t}{\mathrm{d}t}}-{\displaystyle \frac{d\log {\alpha }_t}{\mathrm{d}t}}\right)=-2{\sigma }_t^2{\displaystyle \frac{d{\lambda }_t}{\mathrm{d}t}}.$ (3.2)

Combining with $f(t)=d\log {\alpha }_t/\mathrm{d}t$ in Eq. (2.3), we can rewrite Eq. (3.1) as

${𝒙}_t={\displaystyle \frac{{\alpha }_t}{{\alpha }_s}}{𝒙}_s-{\alpha }_t{\displaystyle {\int }_s^t}\left({\displaystyle \frac{d{\lambda }_{\tau }}{d\tau }}\right){\displaystyle \frac{{\sigma }_{\tau }}{{\alpha }_{\tau }}}{\mathit{ϵ}}_{\theta }({𝒙}_{\tau },\tau )d\tau .$ (3.3)

As $\lambda (t)={\lambda }_t$ is a strictly decreasing function of $t$, it has an inverse function $t_{\lambda }(\cdot )$ satisfying $t=t_{\lambda }(\lambda (t))$. We further change the subscripts of $𝒙$ and ${\mathit{ϵ}}_{\theta }$ from $t$ to $\lambda$ and denote ${\widehat{𝒙}}_{\lambda }≔{𝒙}_{t_{\lambda }(\lambda )}$, ${\widehat{\mathit{ϵ}}}_{\theta }({\widehat{𝒙}}_{\lambda },\lambda )≔{\mathit{ϵ}}_{\theta }({𝒙}_{t_{\lambda }(\lambda )},t_{\lambda }(\lambda ))$. Rewrite Eq. (3.3) by “change-of-variable” for $\lambda$, then we have:

Proposition 3.1 (Exact solution of diffusion ODEs).

Given an initial value ${𝐱}_s$ at time $s>0$, the solution ${𝐱}_t$ at time $t\in [0,s]$ of diffusion ODEs in Eq. (2.7) is:

$${𝒙}_t=\frac{{\alpha }_t}{{\alpha }_s}{𝒙}_s-{\alpha }_t{\int }_{{\lambda }_s}^{{\lambda }_t}{e}^{-\lambda }{\widehat{\mathit{ϵ}}}_{\theta }({\widehat{𝒙}}_{\lambda },\lambda )d\lambda .$$ (3.4)

We call the integral $\int e^{-\lambda }{\widehat{\mathit{ϵ}}}_{\theta }({\widehat{𝒙}}_{\lambda },\lambda )d\lambda$ the exponentially weighted integral of ${\widehat{\mathit{ϵ}}}_{\theta }$, 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 ${𝒙}_s$ at time $s$, According to Eq. (3.4), approximating the solution at time $t$ is equivalent to directly approximating the exponentially weighted integral of ${\widehat{\mathit{ϵ}}}_{\theta }$ from ${\lambda }_s$ to ${\lambda }_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 ${𝒙}_T$ at time $T$ and $M+1$ time steps ${\{t_i\}}_{i=0}^M$ decreasing from $t_0=T$ to $t_M=0$. Let ${\widetilde{𝒙}}_{t_0}={𝒙}_T$ be the initial value. The proposed solvers use $M$ steps to iteratively compute a sequence ${\{{\widetilde{𝒙}}_{t_i}\}}_{i=0}^M$ to approximate the true solutions at time steps ${\{t_i\}}_{i=0}^M$. In particular, the last iterate ${\widetilde{𝒙}}_{t_M}$ approximates the true solution at time $0$.

In order to reduce the approximation error between ${\widetilde{𝒙}}_{t_M}$ and the true solution at time $0$, we need to reduce the approximation error for each ${\widetilde{𝒙}}_{t_i}$ at every step [30]. Starting with the previous value ${\widetilde{𝒙}}_{t_{i-1}}$ at time $t_{i-1}$, according to Eq. (3.4), the exact solution ${𝒙}_{t_{i-1}\to t_i}$ at time $t_i$ is given by

$${𝒙}_{t_{i-1}\to t_i}=\frac{{\alpha }_{t_i}}{{\alpha }_{t_{i-1}}}{\widetilde{𝒙}}_{t_{i-1}}-{\alpha }_{t_i}{\int }_{{\lambda }_{t_{i-1}}}^{{\lambda }_{t_i}}{e}^{-\lambda }{\widehat{\mathit{ϵ}}}_{\theta }({\widehat{𝒙}}_{\lambda },\lambda )d\lambda .$$ (3.5)

Therefore, to compute the value ${\widetilde{𝒙}}_{t_i}$ for approximating ${𝒙}_{t_{i-1}\to t_i}$, we need to approximate the exponentially weighted integral of ${\widehat{\mathit{ϵ}}}_{\theta }$ from ${\lambda }_{t_{i-1}}$ to ${\lambda }_{t_i}$. Denote $h_i≔{\lambda }_{t_i}-{\lambda }_{t_{i-1}}$, and ${\widehat{\mathit{ϵ}}}_{\theta }^{(n)}({\widehat{𝒙}}_{\lambda },\lambda )≔\frac{d^n{\widehat{\mathit{ϵ}}}_{\theta }({\widehat{𝒙}}_{\lambda },\lambda )}{d{\lambda }^n}$ as the $n$-th order total derivative of ${\widehat{\mathit{ϵ}}}_{\theta }({\widehat{𝒙}}_{\lambda },\lambda )$ w.r.t. $\lambda$. For $k\ge 1$, the $(k-1)$-th order Taylor expansion of ${\widehat{\mathit{ϵ}}}_{\theta }({\widehat{𝒙}}_{\lambda },\lambda )$ w.r.t. $\lambda$ at ${\lambda }_{t_{i-1}}$ is

$${\widehat{\mathit{ϵ}}}_{\theta }({\widehat{𝒙}}_{\lambda },\lambda )=\sum _{n=0}^{k-1}\frac{{(\lambda -{\lambda }_{t_{i-1}})}^n}{n!}{\widehat{\mathit{ϵ}}}_{\theta }^{(n)}({\widehat{𝒙}}_{{\lambda }_{t_{i-1}}},{\lambda }_{t_{i-1}})+𝒪({(\lambda -{\lambda }_{t_{i-1}})}^k),$$

Substituting the above Taylor expansion into Eq. (3.5) yields

$${𝒙}_{t_{i-1}\to t_i}=\frac{{\alpha }_{t_i}}{{\alpha }_{t_{i-1}}}{\widetilde{𝒙}}_{t_{i-1}}-{\alpha }_{t_i}\sum _{n=0}^{k-1}{\widehat{\mathit{ϵ}}}_{\theta }^{(n)}({\widehat{𝒙}}_{{\lambda }_{t_{i-1}}},{\lambda }_{t_{i-1}}){\int }_{{\lambda }_{t_{i-1}}}^{{\lambda }_{t_i}}{e}^{-\lambda }\frac{{(\lambda -{\lambda }_{t_{i-1}})}^n}{n!}d\lambda +𝒪(h_i^{k+1}),$$ (3.6)

where the integral $\int e^{-\lambda }\frac{{(\lambda -{\lambda }_{t_{i-1}})}^n}{n!}d\lambda$ can be analytically computed by repeatedly applying $n$ times of integration-by-parts (see Appendix B.2). Therefore, to approximate ${𝒙}_{t_{i-1}\to t_i}$, we only need to approximate the $n$-th order total derivatives ${\widehat{\mathit{ϵ}}}_{\theta }^{(n)}({\widehat{𝒙}}_{\lambda },\lambda )$ for $n\le k-1$, which is a well-studied problem in the ODE literature [31, 32]. By dropping the $𝒪(h_i^{k+1})$ 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

${𝒙}_{t_{i-1}\to t_i}$ $={\displaystyle \frac{{\alpha }_{t_i}}{{\alpha }_{t_{i-1}}}}{\widetilde{𝒙}}_{t_{i-1}}-{\alpha }_{t_i}{\mathit{ϵ}}_{\theta }({\widetilde{𝒙}}_{t_{i-1}},t_{i-1}){\displaystyle {\int }_{{\lambda }_{t_{i-1}}}^{{\lambda }_{t_i}}}{e}^{-\lambda }d\lambda +𝒪(h_i^2)$

$={\displaystyle \frac{{\alpha }_{t_i}}{{\alpha }_{t_{i-1}}}}{\widetilde{𝒙}}_{t_{i-1}}-{\sigma }_{t_i}(e^{h_i}-1){\mathit{ϵ}}_{\theta }({\widetilde{𝒙}}_{t_{i-1}},t_{i-1})+𝒪(h_i^2).$

By dropping the high-order error term $𝒪(h_i^2)$, we can obtain an approximation for ${𝒙}_{t_{i-1}\to t_i}$. 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 ${𝒙}_T$ and $M+1$ time steps ${\{t_i\}}_{i=0}^M$ decreasing from $t_0=T$ to $t_M=0$. Starting with ${\widetilde{𝒙}}_{t_0}={𝒙}_T$, the sequence ${\{{\widetilde{𝒙}}_{t_i}\}}_{i=1}^M$ is computed iteratively as follows:

$${\widetilde{𝒙}}_{t_i}=\frac{{\alpha }_{t_i}}{{\alpha }_{t_{i-1}}}{\widetilde{𝒙}}_{t_{i-1}}-{\sigma }_{t_i}(e^{h_i}-1){\mathit{ϵ}}_{\theta }({\widetilde{𝒙}}_{t_{i-1}},t_{i-1}),\text{where}{h}_i={\lambda }_{t_i}-{\lambda }_{t_{i-1}}.$$ (3.7)

For $k\ge 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.

1:initial value ${𝒙}_T$, time steps ${\{t_i\}}_{i=0}^M$, model ${\mathit{ϵ}}_{\theta }$

2:${\widetilde{𝒙}}_{t_0}\leftarrow {𝒙}_T$

3:for $i\leftarrow 1$ to $M$ do

4: $s_i\leftarrow t_{\lambda }\left(\frac{{\lambda }_{t_{i-1}}+{\lambda }_{t_i}}{2}\right)$

5: ${𝒖}_i\leftarrow \frac{{\alpha }_{s_i}}{{\alpha }_{t_{i-1}}}{\widetilde{𝒙}}_{t_{i-1}}-{\sigma }_{s_i}\left(e^{\frac{h_i}{2}}-1\right){\mathit{ϵ}}_{\theta }({\widetilde{𝒙}}_{t_{i-1}},t_{i-1})$

6: ${\widetilde{𝒙}}_{t_i}\leftarrow \frac{{\alpha }_{t_i}}{{\alpha }_{t_{i-1}}}{\widetilde{𝒙}}_{t_{i-1}}-{\sigma }_{t_i}\left(e^{h_i}-1\right){\mathit{ϵ}}_{\theta }({𝒖}_i,s_i)$

7:end for

8:return ${\widetilde{𝒙}}_{t_M}$

Algorithm 2 DPM-Solver-3.

1:initial value ${𝒙}_T$, time steps ${\{t_i\}}_{i=0}^M$, model ${\mathit{ϵ}}_{\theta }$

2:${\widetilde{𝒙}}_{t_0}\leftarrow {𝒙}_T$, $r_1\leftarrow \frac{1}{3}$, $r_2\leftarrow \frac{2}{3}$

3:for $i\leftarrow 1$ to $M$ do

4: $s_{2i-1}\leftarrow t_{\lambda }\left({\lambda }_{t_{i-1}}+r_1{h}_i\right),s_{2i}\leftarrow t_{\lambda }\left({\lambda }_{t_{i-1}}+r_2{h}_i\right)$

5: ${𝒖}_{2i-1}\leftarrow \frac{{\alpha }_{s_{2i-1}}}{{\alpha }_{t_{i-1}}}{\widetilde{𝒙}}_{t_{i-1}}-{\sigma }_{s_{2i-1}}\left(e^{r_1{h}_i}-1\right){\mathit{ϵ}}_{\theta }({\widetilde{𝒙}}_{t_{i-1}},t_{i-1})$

6: ${𝑫}_{2i-1}\leftarrow {\mathit{ϵ}}_{\theta }({𝒖}_{2i-1},s_{2i-1})-{\mathit{ϵ}}_{\theta }({\widetilde{𝒙}}_{t_{i-1}},t_{i-1})$

7: ${𝒖}_{2i}\leftarrow \frac{{\alpha }_{s_{2i}}}{{\alpha }_{t_{i-1}}}{\widetilde{𝒙}}_{t_{i-1}}-{\sigma }_{s_{2i}}\left(e^{r_2{h}_i}-1\right){\mathit{ϵ}}_{\theta }({\widetilde{𝒙}}_{t_{i-1}},t_{i-1})-\frac{{\sigma }_{s_{2i}}{r}_2}{r_1}\left(\frac{e^{r_2{h}_i}-1}{r_2{h}_i}-1\right){𝑫}_{2i-1}$

8: ${𝑫}_{2i}\leftarrow {\mathit{ϵ}}_{\theta }({𝒖}_{2i},s_{2i})-{\mathit{ϵ}}_{\theta }({\widetilde{𝒙}}_{t_{i-1}},t_{i-1})$

9: ${\widetilde{𝒙}}_{t_i}\leftarrow \frac{{\alpha }_{t_i}}{{\alpha }_{t_{i-1}}}{\widetilde{𝒙}}_{t_{i-1}}-{\sigma }_{t_i}\left(e^{h_i}-1\right){\mathit{ϵ}}_{\theta }({\widetilde{𝒙}}_{t_{i-1}},t_{i-1})-\frac{{\sigma }_{t_i}}{r_2}\left(\frac{e^{h_i}-1}{h}-1\right){𝑫}_{2i}$

10:end for

11:return ${\widetilde{𝒙}}_{t_M}$

Here, $t_{\lambda }(\cdot )$ is the inverse function of $\lambda (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 ($s_i$, ${𝒖}_i$) for DPM-Solver-2 and $(s_{2i-1},{𝒖}_{2i-1})$ and $(s_{2i},{𝒖}_{2i})$ for DPM-Solver-3. As shown in the algorithm, DPM-Solver-$k$ requires $k$ function evaluations per step for $k=1,2,3$. Despite the more expensive steps, higher-order solvers ($k=2,3$) are usually more efficient since they require much fewer steps to converge, due to their higher convergence order. We show that DPM-Solver-$k$ is $k$-th-order solver, as stated in the following theorem. The proof is in Appendix B.

Theorem 3.2 (DPM-Solver-$k$ as a $k$-th-order solver).

Assume ${\mathbf{ϵ}}_{\theta }({𝐱}_t,t)$ follows the regularity conditions detailed in Appendix B.1, then for $k=1,2,3$, DPM-Solver-$k$ is a $k$-th order solver for diffusion ODEs, i.e., for the sequence ${\{{\widetilde{𝐱}}_{t_i}\}}_{i=1}^M$ computed by DPM-Solver-$k$, the approximation error at time $0$ satisfies ${\widetilde{𝐱}}_{t_M}-{𝐱}_0=𝒪(h_{\max }^k)$, where $h_{\text{max}}={\max }_{1\le i\le M}({\lambda }_{t_i}-{\lambda }_{t_{i-1}})$.

Finally, solvers with $k\ge 4$ need much more intermediate points as shown by previous work [31, 32] for exponential integrators. Therefore, we only consider $k$ from $1$ to $3$ in this work, while leaving the solvers with higher $k$ for future study.

3.3 Step Size Schedule

The proposed solvers in Sec. 3.2 need to specify the time steps ${\{t_i\}}_{i=0}^M$ in advance. We propose two choices of the time step schedule. One choice is handcrafted, which is to uniformly split the interval $[{\lambda }_T$, ${\lambda }_0$], i.e. ${\lambda }_{t_i}={\lambda }_T+\frac{i}{M}({\lambda }_0-{\lambda }_T)$, $i=0,\mathrm{\dots },M$. Note that this is different from previous work [2, 3] which chooses uniform steps for $t_i$. Empirically, DPM-Solver with uniform time steps ${\lambda }_{t_i}$ can already generate quite good samples in few steps, where results are listed in Appendix E. As the other choice, we propose an adaptive step size algorithm, which dynamically adjusts the step size by combining different orders of DPM-Solver. The adaptive algorithm is inspired by [20] and we defer its implementation details to Appendix C.

For few-step sampling, we need to use up all the number of function evaluations (NFE). When the NFE is not divisible by $3$, we firstly apply DPM-Solver-3 as much as possible, and then add a single step of DPM-Solver-1 or DPM-Solver-2 (dependent on the reminder of $K$ divided by $3$), as detailed in Appendix D. In the subsequent experiments, we use such combination of solvers with the uniform step size schedule for NFE $\le 20$, and otherwise the adaptive step size schedule.

3.4 Sampling from Discrete-Time DPMs

Discrete-time DPMs [2] train the noise prediction model at $N$ fixed time steps ${\{t_n\}}_{n=1}^N$, and the noise prediction model is parameterized by ${\widetilde{\mathit{ϵ}}}_{\theta }({𝒙}_n,n)$ for $n=0,\mathrm{\dots },N-1$, where each ${𝒙}_n$ is corresponding to the value at time $t_{n+1}$. We can transform the discrete-time noise prediction model to the continuous version by letting ${\mathit{ϵ}}_{\theta }(𝒙,t)≔{\widetilde{\mathit{ϵ}}}_{\theta }(𝒙,\frac{(N-1)t}{T})$, for all $𝒙\in {ℝ}^d,t\in [0,T]$. Note that the input time of ${\widetilde{\mathit{ϵ}}}_{\theta }$ may not be integers, but we find that the noise prediction model can still work well, and we hypothesize that it is because of the smooth time embeddings (e.g., position embeddings [2]). By such reparameterization, the noise prediction model can adopt the continuous-time steps as input, and thus we can also use DPM-Solver for fast sampling.

4 Comparison with Existing Fast Sampling Methods

Here, we discuss the relationship and highlight the difference between DPM-Solver and existing ODE-based fast sampling methods for DPMs. We further briefly discuss the advantage of training-free samplers over those training-based ones.

4.1 DDIM as DPM-Solver-1

Denoising Diffusion Implicit Models (DDIM) [19] design a deterministic method for fast sampling from DPMs. For two adjacent time steps $t_{i-1}$ and $t_i$, assume that we have a solution ${\widetilde{𝒙}}_{t_{i-1}}$ at time $t_{i-1}$, then a single step of DDIM from time $t_{i-1}$ to time $t_i$ is

$${\widetilde{𝒙}}_{t_i}=\frac{{\alpha }_{t_i}}{{\alpha }_{t_{i-1}}}{\widetilde{𝒙}}_{t_{i-1}}-{\alpha }_{t_i}\left(\frac{{\sigma }_{t_{i-1}}}{{\alpha }_{t_{i-1}}}-\frac{{\sigma }_{t_i}}{{\alpha }_{t_i}}\right){\mathit{ϵ}}_{\theta }({\widetilde{𝒙}}_{t_{i-1}},t_{i-1}).$$ (4.1)

Although motivated by entirely different perspectives, we show that the updates of DPM-Solver-1 and Denoising Diffusion Implicit Models (DDIM) [19] are identical. By the definition of $\lambda$, we have $\frac{{\sigma }_{t_{i-1}}}{{\alpha }_{t_{i-1}}}=e^{-{\lambda }_{t_{i-1}}}$ and $\frac{{\sigma }_{t_i}}{{\alpha }_{t_i}}=e^{-{\lambda }_{t_i}}$. Plugging these and $h_i={\lambda }_{t_i}-{\lambda }_{t_{i-1}}$ to Eq. (4.1) results in exactly a step of DPM-Solver-1 in Eq. (3.7). However, the semi-linear ODE formulation of DPM-Solver allows for principled generalization to higher-order solvers and convergence order analysis.

Recent work [13] also show that DDIM is a first-order discretization of diffusion ODEs by differentiating both sides of Eq. (4.1). However, they cannot explain the difference between DDIM and the first-order Euler discretization of diffusion ODEs. In contrast, by showing that DDIM is a special case of DPM-Solver, we reveal that DDIM makes full use of the semi-linearity of diffusion ODEs, which explains its superiority over traditional Euler methods.

4.2 Comparison with Traditional Runge-Kutta Methods

One can obtain a high-order solver by directly applying traditional explicit Runge-Kutta (RK) methods to the diffusion ODE in Eq. (2.7). Specifically, RK methods write the solution of Eq. (2.7) in the following integral form:

$${𝒙}_t={𝒙}_s+{\int }_s^t{𝒉}_{\theta }({𝒙}_{\tau },\tau )d\tau ={𝒙}_s+{\int }_s^t\left(f(\tau ){𝒙}_{\tau }+\frac{g^2(\tau )}{2{\sigma }_{\tau }}{\mathit{ϵ}}_{\theta }({𝒙}_{\tau },\tau )\right)d\tau ,$$ (4.2)

and use some intermediate time steps between $[t,s]$ and combine the evaluations of ${𝒉}_{\theta }$ at these time steps to approximate the whole integral. The approximation error of explicit RK methods depends on ${𝒉}_{\theta }$, which consists of the error corresponding to both the linear term $f(\tau ){𝒙}_{\tau }$ and the nonlinear noise prediction model ${\mathit{ϵ}}_{\theta }$. However, the error of the linear term may increase exponentially because the exact solution of the linear term has an exponential coefficient (as shown in Eq. (3.1)). There are many empirical evidence [25, 31] showing that directly using explicit RK methods for semi-linear ODEs may suffer from unstable numerical issues for large step size. We also demonstrate the empirical difference of the proposed DPM-Solver and the traditional explicit RK methods in Sec. 5.1, which shows that DPM-Solver have smaller discretization errors than the RK methods with the same order.

4.3 Training-based Fast Sampling Methods for DPMs

Samplers that need extra training or optimization include knowledge distillation [13, 14], learning the noise level or variance [15, 16, 33], and learning the noise schedule or sample trajectory [17, 18]. Although the progressive distillation method [13] can obtain a fast sampler within 4 steps, it needs further training costs and loses part of the information in the original DPM (e.g., after distillation, the noise prediction model cannot predict the noise (score function) at every time step between $[0,T]$). In contrast, training-free samplers can keep all the information of the original model, and thereby can be directly extended to the conditional sampling by combining the original model and an external classifier [4] (e.g. see Appendix D for the conditional sampling with classifier guidance).

Beyond directly designing fast samplers for DPMs, several works also propose novel types of DPMs which supports faster sampling. For instance, defining a low-dimensional latent variable for DPMs [34]; designing special diffusion processes with bounded score functions [35]; combining GANs with the reverse process of DPMs [36]. The proposed DPM-Solver may also be suitable for accelerating the sampling of these DPMs, and we leave them for future work.

Table 1: FID $\downarrow$ on CIFAR-10 for different orders of Runge-Kutta (RK) methods and DPM-Solvers, varying the number of function evaluations (NFE). For RK methods, we evaluate diffusion ODEs w.r.t. both $t$ (Eq. (2.7)) and $\lambda$ (Eq. (E.1)). We use uniform step size in $t$ for RK ($t$), and uniform step size in $\lambda$ for RK ($\lambda$) and DPM-Solvers.

Sampling method $\backslash$ NFE	12	18	24	30	36	42	48
RK2 ($t$)	16.40	7.25	3.90	3.63	3.58	3.59	3.54
RK2 ($\lambda$)	107.81	42.04	17.71	7.65	4.62	3.58	3.17
DPM-Solver-2	5.28	3.43	3.02	2.85	2.78	2.72	2.69
RK3 ($t$)	48.75	21.86	10.90	6.96	5.22	4.56	4.12
RK3 ($\lambda$)	34.29	4.90	3.50	3.03	2.85	2.74	2.69
DPM-Solver-3	6.03	2.90	2.75	2.70	2.67	2.65	2.65

(a) CIFAR-10 (continuous)

(b) CIFAR-10 (discrete)

(c) CelebA 64x64 (discrete)

(d) ImageNet 64x64 (discrete)

(e) ImageNet 128x128 (discrete)

(f) LSUN bedroom 256x256 (discrete)

Figure 2: Sample quality measured by FID $\downarrow$ of different sampling methods for DPMs on CIFAR-10 with both continuous-time and discrete-time models, CelebA 64x64, ImageNet 64x64, ImageNet 128x128 and LSUN bedroom 256x256 with discrete-time models, varying the number of function evaluations (NFE). The method ^†GGDM [18] needs extra training to optimize the sample trajectory, while other methods are training-free. To get the strongest baseline, we use the quadratic step size for DDIM on CelebA, which has a better FID than that of the uniform step size in the original paper [19].

5 Experiments

In this section, we show that as a training-free sampler, DPM-Solver can greatly speedup the sampling of existing pre-trained DPMs, including both continuous-time and discrete-time ones, with both linear noise schedule [2, 19] and cosine noise schedule [16]. We vary different number of function evaluations (NFE) which is the number of calls to the noise prediction model ${\mathit{ϵ}}_{\theta }({𝒙}_t,t)$, and compare the sample quality between DPM-Solver and other methods. For each experiment, We draw 50K samples and use the widely adopted FID score [37] to evaluate the sample quality, where lower FID usually implies better sample quality.

Unless explicitly mentioned, we always use the solver combination with the uniform step size schedule in Sec. 3.3 if the NFE budget is less than 20, and otherwise the DPM-Solver-3 with the adaptive step size schedule in Sec. 3.3. We refer to Appendix D for other implementation details of DPM-Solver and Appendix E for detailed settings.

5.1 Comparison with Continuous-Time Sampling Methods

We firstly compare DPM-Solver with other continuous-time sampling methods for DPMs. The compared methods include the Euler-Maruyama discretization for diffusion SDEs [3], the adaptive step size solver for diffusion SDEs [20] and the RK methods for diffusion ODEs [3, 28] in Eq. (2.7). We compare these methods for sampling from a pre-trained continuous-time “VP deep” model [3] on the CIFAR-10 dataset [29] with the linear noise schedule.

Fig. 2(a) shows the efficiency of compared solvers. We use uniform time steps with 50, 200, 1000 NFE for the diffusion SDE with Euler discretization, and vary the tolerance hyperparameter [3, 20] for the adaptive step size SDE solver [20] and RK45 ODE solver [28] to control the NFE. DPM-Solver can generate good sample quality within around 10 NFE, while other solvers have large discretization error even in 50 NFE, which shows that DPM-Solver can achieve $\sim$5 speedup of the previous best solver. In particular, we achieve 4.70 FID with 10 NFE, 3.75 FID with 12 NFE, 3.24 FID with 15 NFE, and 2.87 FID with 20 NFE, which is the fastest sampler on CIFAR-10.

As an ablation study, we also compare the second-order and third-order DPM-Solver and RK methods, as shown in Table 1. We compare RK methods for diffusion ODEs w.r.t. both time $t$ in Eq. (2.7) and half-log-SNR $\lambda$ by applying change-of-variable (see detailed formulations in Appendix E.1). The results show that given the same NFE, the sample quality of DPM-Solver is consistently better than RK methods with the same order. The superior efficiency of DPM-Solver is particularly evident in the few-step regime under 15 NFE, where RK methods have rather large discretization errors. This is mainly because DPM-Solver analytically computes the linear term, avoiding the corresponding discretization error. Besides, the higher order DPM-Solver-3 converges faster than DPM-Solver-2, which matches the order analysis in Theorem 3.2.

5.2 Comparison with Discrete-Time Sampling Methods

We use the method in Sec. 3.4 for using DPM-Solver in discrete-time DPMs, and then compare DPM-Solver with other discrete-time training-free samplers, including DDPM [2], DDIM [19], Analytic-DDPM [21], Analytic-DDIM [21], PNDM [22], FastDPM [38] and Itô-Taylor [24]. We also compare with GGDM [18], which uses the same pre-trained model but needs further training for the sampling trajectory. We compare the sample quality by varying NFE from 10 to 1000.

Specifically, we use the discrete-time model trained by $L_{\text{simple}}$ in [2] on the CIFAR-10 dataset with linear noise schedule; the discrete-time model in [19] on CelebA 64x64 [39] with linear noise schedule; the discrete-time model trained by $L_{\text{hybrid}}$ in [16] on ImageNet 64x64 [26] with cosine noise schedule; the discrete-time model with classifier guidance in [4] on ImageNet 128x128 [26] with linear noise schedule; the discrete-time model in [4] on LSUN bedroom 256x256 [40] with linear noise schedule. For the models trained on ImageNet, we only use their “mean” model and omit the “variance” model. As shown in Fig. 2, on all datasets, DPM-Solver can obtain reasonable samples within 12 steps (FID 4.65 on CIFAR-10, FID 3.71 on CelebA 64x64 and FID 19.97 on ImageNet 64x64, FID 4.08 on ImageNet 128x128), which is $4\sim 16\times$ faster than the previous fastest training-free sampler. DPM-Solver even outperforms GGDM, which requires additional training.

6 Conclusions

We tackle the problem of fast and training-free sampling from DPMs. We propose DPM-Solver, a fast dedicated training-free solver of diffusion ODEs for fast sampling of DPMs in around 10 steps of function evaluations. DPM-Solver leverages the semi-linearity of diffusion ODEs and it directly approximates a simplified formulation of exact solutions of diffusion ODEs, which consists of an exponentially weighted integral of the noise prediction model. Inspired by numerical methods for exponential integrators, we propose first-order, second-order and third-order DPM-Solver to approximate the exponentially weighted integral of noise prediction models with theoretical convergence guarantee. We propose both handcrafted and adaptive step size schedule, and apply DPM-Solver for both continuous-time and discrete-time DPMs. Our experimental results show that DPM-Solver can generate high-quality samples in around 10 function evaluations on various datasets, and it can achieve $4\sim 16\times$ speedup compared with previous state-of-the-art training-free samplers.

Limitations and broader impact Despite the promising speedup performance, DPM-Solver is designed for fast sampling, which may be not suitable for accelerating the likelihood evaluations of DPMs. Besides, compared to the commonly-used GANs, diffusion models with DPM-Solver are still not fast enough for real-time applications. In addition, like other deep generative models, DPMs may be used to generate adverse fake contents, and the proposed solver may further amplify the potential undesirable influence of deep generative models for malicious applications.

Checklist

1. For all authors…

   1. Do the main claims made in the abstract and introduction accurately reflect the paper’s contributions and scope? [Yes]

   2. Did you describe the limitations of your work? [Yes] See section 6.

   3. Did you discuss any potential negative societal impacts of your work? [Yes] See section 6.

   4. Have you read the ethics review guidelines and ensured that your paper conforms to them? [Yes]

2. If you are including theoretical results…

   1. Did you state the full set of assumptions of all theoretical results? [Yes] See Appendix B.

   2. Did you include complete proofs of all theoretical results? [Yes] See Appendix B.

3. If you ran experiments…

   1. Did you include the code, data, and instructions needed to reproduce the main experimental results (either in the supplemental material or as a URL)? [Yes] Code is attached in the supplemental materials, with the appendix.

   2. Did you specify all the training details (e.g., data splits, hyperparameters, how they were chosen)? [Yes] Our method is training-free. But we also report the hyperparameters for evaluations used in our proposed solver.

   3. Did you report error bars (e.g., with respect to the random seed after running experiments multiple times)? [No] We observe that the standard deviation of the FID evaluations of DPM-Solver are rather small (mainly less than 0.01) because the FID is already averaged over 50K samples, following existing work [18, 20, 21]. The small standard deviation does not change the conclusion.

   4. Did you include the total amount of compute and the type of resources used (e.g., type of GPUs, internal cluster, or cloud provider)? [Yes] The GPU type and amount is detailed in Appendix E.

4. If you are using existing assets (e.g., code, data, models) or curating/releasing new assets…

   1. If your work uses existing assets, did you cite the creators? [Yes]

   2. Did you mention the license of the assets? [Yes] See Appendix E

   3. Did you include any new assets either in the supplemental material or as a URL? [Yes] We include our code in the supplemental materials.

   4. Did you discuss whether and how consent was obtained from people whose data you’re using/curating? [No] All of the datasets used in the experiments are publicly available.

   5. Did you discuss whether the data you are using/curating contains personally identifiable information or offensive content? [Yes] We mentioned the human privacy issues of the ImageNet dataset in Appendix E.

5. If you used crowdsourcing or conducted research with human subjects…

   1. Did you include the full text of instructions given to participants and screenshots, if applicable? [N/A]

   2. Did you describe any potential participant risks, with links to Institutional Review Board (IRB) approvals, if applicable? [N/A]

   3. Did you include the estimated hourly wage paid to participants and the total amount spent on participant compensation? [N/A]

Appendix A Sampling with Invariance to the Noise Schedule

Table 2: Formulations that are invariant to the choice of the noise schedules. The maximum likelihood training loss w.r.t. $\lambda$ is equivalent to the objectives in [10, 41], and the exact solution of the diffusion ODEs are proposed in Proposition 3.1.

Method	Invariance Formulation
Maximum likelihood training	${\displaystyle {\int }_{{\lambda }_T}^{{\lambda }_0}}{𝔼}_{q_0({𝒙}_0)}{𝔼}_{\mathit{ϵ}\sim 𝒩(\mathrm{𝟎},𝑰)}\left[{\Vert {\mathit{ϵ}}_{\theta }({\widehat{𝒙}}_{\lambda },\lambda )-\mathit{ϵ}\Vert }_2^2\right]d\lambda$
Sampling by diffusion ODEs	${\widehat{𝒙}}_{{\lambda }_t}={\displaystyle \frac{{\widehat{\alpha }}_{{\lambda }_t}}{{\widehat{\alpha }}_{{\lambda }_s}}}{\widehat{𝒙}}_{{\lambda }_s}-{\widehat{\alpha }}_{{\lambda }_t}{\displaystyle {\int }_{{\lambda }_s}^{{\lambda }_t}}{e}^{-\lambda }{\widehat{\mathit{ϵ}}}_{\theta }({\widehat{𝒙}}_{\lambda },\lambda )d\lambda$

In this section, we discuss more about the exact solution in Proposition 3.1 and give some insights about the formulation. Below we firstly restate the proposition w.r.t. $\lambda$ (i.e. the half-logSNR).

Proposition 3.1 (Exact solution of diffusion ODEs).

Given an initial value ${\widehat{𝐱}}_{{\lambda }_s}$ at time $s$ with the corresponding half-logSNR ${\lambda }_s$, the solution ${\widehat{𝐱}}_{{\lambda }_t}$ at time $t$ of diffusion ODEs in Eq. (2.7) with the corresponding half-logSNR ${\lambda }_t$ is:

$${\widehat{𝒙}}_{{\lambda }_t}=\frac{{\alpha }_t}{{\alpha }_s}{\widehat{𝒙}}_{{\lambda }_s}-{\alpha }_t{\int }_{{\lambda }_s}^{{\lambda }_t}{e}^{-\lambda }{\widehat{\mathit{ϵ}}}_{\theta }({\widehat{𝒙}}_{\lambda },\lambda )d\lambda .$$ (A.1)

In the following subsections, we will show that such formulation decouples the model ${\mathit{ϵ}}_{\theta }$ from the specific noise schedule, and thus is invariant to the noise schedule. Moreover, such change-of-variable for $\lambda$ in Proposition 3.1 is highly related to the maximum likelihood training of diffusion models [10, 41]. We show that both the maximum likelihood training and the sampling of diffusion models have invariance formulations that are independent of the noise schedule.

A.1 Decoupling the Sampling Solution from the Noise Schedule

In this section, we show that Proposition 3.1 can decouples the exact solutions of the diffusion ODEs from the specific noise schedules (i.e. choice of the functions ${\alpha }_t=\alpha (t)$ and ${\sigma }_t=\sigma (t)$). Namely, given a starting point ${\lambda }_s$, a ending point ${\lambda }_t$, an initial value ${\widehat{𝒙}}_{{\lambda }_s}$ at ${\lambda }_s$ and a noise prediction model ${\widehat{\mathit{ϵ}}}_{\theta }$, the solution of ${\widehat{x}}_{{\lambda }_t}$ is invariant of the noise schedule between ${\lambda }_s$ and ${\lambda }_t$.

We firstly consider the VP type diffusion models, which is equivalent to the original DDPM [2, 3]. For VP type diffusion models, we always have ${\alpha }_t^2+{\sigma }_t^2=1$, so defining the noise schedule is equivalent to defining the function ${\alpha }_t=\alpha (t)$ (For example, DDPM [2] uses a noise schedule such that $\beta (t)=\frac{d\log {\alpha }_t}{dt}$ is a linear function of $t$, and i-DDPM [16] uses a noise schedule such that $\beta (t)=\frac{d\log {\alpha }_t}{dt}$ is a cosine function of $t$). As ${\lambda }_t=\log {\alpha }_t-\log {\sigma }_t$, we have ${\alpha }_t=\sqrt{\frac{1}{1+e^{-2{\lambda }_t}}}$ and ${\sigma }_t=\sqrt{\frac{1}{1+e^{2{\lambda }_t}}}$. Thus, we can directly compute the ${\alpha }_t$ and ${\sigma }_t$ for a given ${\lambda }_t$. Denote ${\widehat{\alpha }}_{\lambda }≔\sqrt{\frac{1}{1+e^{-2\lambda }}}$, we have

$${\widehat{𝒙}}_{{\lambda }_t}=\frac{{\widehat{\alpha }}_{{\lambda }_t}}{{\widehat{\alpha }}_{{\lambda }_s}}{\widehat{𝒙}}_{{\lambda }_s}-{\widehat{\alpha }}_{{\lambda }_t}{\int }_{{\lambda }_s}^{{\lambda }_t}{e}^{-\lambda }{\widehat{\mathit{ϵ}}}_{\theta }({\widehat{𝒙}}_{\lambda },\lambda )d\lambda .$$ (A.2)

We should notice that the integrand $e^{-\lambda }{\widehat{\mathit{ϵ}}}_{\theta }({\widehat{𝒙}}_{\lambda },\lambda )$ is a function of $\lambda$, so its integral from ${\lambda }_s$ to ${\lambda }_t$ is only dependent on the starting point ${\lambda }_s$, the ending point ${\lambda }_t$ and the function ${\widehat{\mathit{ϵ}}}_{\theta }$, which is independent of the intermediate values. As other coefficients (${\widehat{\alpha }}_{{\lambda }_s}$ and ${\widehat{\alpha }}_{{\lambda }_t}$) are also only dependent on the starting point ${\lambda }_s$ and the ending point ${\lambda }_t$, we can conclude that ${\widehat{𝒙}}_{{\lambda }_t}$ is invariant of the specific choice of the noise schedules. Intuitively, this is because we converts the original integral of time $t$ in Eq. (3.1) to the integral of $\lambda$, and the functions $f(t)$ and $g(t)$ are converted to an analytical formulation $e^{-\lambda }$, which is invariant to the specific choices of $f(t)$ and $g(t)$. Finally, for other types of diffusion models (such as the VE type and the subVP type), they are all equivalent to the VP type by equivalently rescaling the noise prediction models, as proved in [10]. Therefore, the solutions of these types also have such property.

In summary, Proposition 3.1 decouples the solution of diffusion ODEs from the noise schedules, which gives us an opportunity to design tailor-made samplers for DPMs. In fact, as shown in Sec. 3.2, the only approximation of the proposed DPM-Solver is about the Taylor expansion of the neural network ${\widehat{\mathit{ϵ}}}_{\theta }$ w.r.t. $\lambda$, and DPM-Solver analytically computes other coefficients (which are corresponding to the specific noise schedules). Intuitively, DPM-Solver keeps the known information as much as possible, and only approximates the intractable integral of the neural network, so it can generate comparable samples within much fewer steps.

A.2 Choosing Time Steps for $\lambda$ is Invariant to the Noise Schedule

As mentioned in Appendix A.1, the formulation of Proposition 3.1 decouples the sampling solution from the noise schedule. The solution depends on the starting point ${\lambda }_s$ and the ending point ${\lambda }_t$, and is invariant to the intermediate noise schedule. Similarly, the updating equations of the algorithm of DPM-Solver are also invariant to the intermediate noise schedule. Therefore, if we have chosen the time steps ${\{{\lambda }_i\}}_{i=0}^M$, then the solution of DPM-Solver is also determined and is invariant to the intermediate noise schedule.

A simple way for choosing time steps for $\lambda$ is uniformly splitting $[{\lambda }_T,{\lambda }_{ϵ}]$, which is the setting in our experiments. However, we believe that there exists more precise ways for choosing the time steps, and we leave it for future work.

A.3 Relationship with the Maximum Likelihood Training of Diffusion Models

Interestingly, the maximum likelihood training of diffusion SDEs in continuous time also has such invariance property [10]. Below we briefly review the maximum likelihood training loss of diffusion SDEs, and then propose a new insight for understanding diffusion models.

Denote the data distribution as $q_0({𝒙}_0)$, the distribution of the forward process at each time $t$ as $q_t({𝒙}_t)$, the distribution of the reverse process at each time $t$ as $p_t({𝒙}_t)$ with $p_T=𝒩(\mathrm{𝟎},𝑰)$. In [3], it is proved that the KL-divergence between $q_0$ and $p_0$ can be bounded by a weighted score matching loss:

$$D_{\mathrm{KL}}(q_0\parallel p_0)\le D_{\mathrm{KL}}(q_T\parallel p_T)+\frac{1}{2}{\int }_0^T\frac{g^2(t)}{{\sigma }_t^2}{𝔼}_{q_0({𝒙}_0)}{𝔼}_{\mathit{ϵ}\sim 𝒩(\mathrm{𝟎},𝑰)}\left[{\Vert {\mathit{ϵ}}_{\theta }({𝒙}_t,t)-\mathit{ϵ}\Vert }_2^2\right]dt+C,$$ (A.3)

where ${𝒙}_t={\alpha }_t{𝒙}_0+{\sigma }_t\mathit{ϵ}$ and $C$ is a constant independent of $\theta$. As shown in Sec. 3.1, we have

$g^2(t)={\displaystyle \frac{d{\sigma }_t^2}{\mathrm{d}t}}-2{\displaystyle \frac{d\log {\alpha }_t}{\mathrm{d}t}}{\sigma }_t^2=2{\sigma }_t^2\left({\displaystyle \frac{d\log {\sigma }_t}{\mathrm{d}t}}-{\displaystyle \frac{d\log {\alpha }_t}{\mathrm{d}t}}\right)=-2{\sigma }_t^2{\displaystyle \frac{d{\lambda }_t}{\mathrm{d}t}},$ (A.4)

so by applying change-of-variable w.r.t. $\lambda$, we have

$$D_{\mathrm{KL}}(q_0\parallel p_0)\le D_{\mathrm{KL}}(q_T\parallel p_T)+{\int }_{{\lambda }_T}^{{\lambda }_0}{𝔼}_{q_0({𝒙}_0)}{𝔼}_{\mathit{ϵ}\sim 𝒩(\mathrm{𝟎},𝑰)}\left[{\Vert {\mathit{ϵ}}_{\theta }({\widehat{𝒙}}_{\lambda },\lambda )-\mathit{ϵ}\Vert }_2^2\right]d\lambda +C,$$ (A.5)

which is equivalent to the importance sampling trick in [41, Sec. 5.1] and the continuous-time diffusion loss in [10, Eq. (22)]. Compared to Proposition 3.1, we can find that the sampling and the maximum likelihood training of diffusion models can both be converted to an integral w.r.t. $\lambda$, such that the formulation is invariant to the specific noise schedules, and we summarize it in Table 2. Such invariance property for both training and sampling brings a new insight for understanding diffusion models. For instance, we can directly define the noise prediction model ${\mathit{ϵ}}_{\theta }$ w.r.t. the (half-)logSNR $\lambda$ instead of the time $t$, then the training and sampling for diffusion models can be done without further choosing any ad-hoc noise schedules. Such finding may unify the different ways of the training and the inference of diffusion models, and we leave it for future study.

Appendix B Proof of Theorem 3.2

B.1 Assumptions

Throughout this section, we denote ${𝒙}_s$ as the solution of the diffusion ODE Eq. (2.7) starting from ${𝒙}_T$. For DPM-Solver-$k$ we make the following assumptions:

Assumption B.1.

The total derivatives $\frac{d^j{\widehat{\mathit{ϵ}}}_{\theta }({\widehat{𝒙}}_{\lambda },\lambda )}{d{\lambda }^j}$ (as a function of $\lambda$) exist and are continuous for $0\le j\le k+1$.

Assumption B.2.

The function ${\mathit{ϵ}}_{\theta }(𝒙,s)$ is Lipschitz w.r.t. to its first parameter $𝒙$.

Assumption B.3.

$h_{max}=𝒪(1/M)$.

We note that the first assumption is required by Taylor’s theorem Eq. (3.6), and the second assumption is used to replace ${ϵ}_{\theta }({\widetilde{𝒙}}_s,s)$ with ${ϵ}_{\theta }({𝒙}_s,s)+𝒪({𝒙}_s-{\widetilde{𝒙}}_s)$ so that the Taylor expansion w.r.t. ${\lambda }_s$ is applicable. The last one is a technical assumption to exclude a significantly large step-size.

B.2 General Expansion of the Exponentially Weighted Integral

Firstly, we derive the Taylor expansion of the exponentially weighted integral. Let and then ${\lambda }_t>{\lambda }_s$. Denote $h≔{\lambda }_t-{\lambda }_s$, and the $k$-th order total derivative ${\widehat{\mathit{ϵ}}}_{\theta }^{(k)}({\widehat{𝒙}}_{\lambda },\lambda )≔\frac{d^k{\widehat{\mathit{ϵ}}}_{\theta }({\widehat{𝒙}}_{\lambda },\lambda )}{d{\lambda }^k}$. For $n\ge 0$, the $n$-th order Taylor expansion of ${\widehat{\mathit{ϵ}}}_{\theta }({\widehat{𝒙}}_{\lambda },\lambda )$ w.r.t. $\lambda$ is

$${\widehat{\mathit{ϵ}}}_{\theta }({\widehat{𝒙}}_{\lambda },\lambda )=\sum _{k=0}^n\frac{{(\lambda -{\lambda }_s)}^k}{k!}{\widehat{\mathit{ϵ}}}_{\theta }^{(k)}({\widehat{𝒙}}_{{\lambda }_s},{\lambda }_s)+𝒪(h^{n+1}).$$ (B.1)

To expand the exponential integrator, we further define [31]:

$${\phi }_k(z)≔{\int }_0^1{e}^{(1-\delta )z}\frac{{\delta }^{k-1}}{(k-1)!}d\delta ,{\phi }_0(z)=e^z$$ (B.2)

and it satisfies ${\phi }_k(0)=\frac{1}{k!}$ and a recurrence relation ${\phi }_{k+1}(z)=\frac{{\phi }_k(z)-{\phi }_k(0)}{z}$. By taking the Taylor expansion of ${\widehat{\mathit{ϵ}}}_{\theta }({\widehat{𝒙}}_{\lambda },\lambda )$, the exponential integrator can be rewritten as

$${\int }_{{\lambda }_s}^{{\lambda }_t}{e}^{-\lambda }{\widehat{\mathit{ϵ}}}_{\theta }({\widehat{𝒙}}_{\lambda },\lambda )d\lambda =\frac{{\sigma }_t}{{\alpha }_t}\sum _{k=0}^n{h}^{k+1}{\phi }_{k+1}(h){\widehat{\mathit{ϵ}}}_{\theta }^{(k)}({\widehat{𝒙}}_{{\lambda }_s},{\lambda }_s)+𝒪(h^{n+2}).$$ (B.3)

So the solution of ${𝒙}_t$ in Eq. (3.4) can be expanded as

$${𝒙}_t=\frac{{\alpha }_t}{{\alpha }_s}{𝒙}_s-{\sigma }_t\sum _{k=0}^n{h}^{k+1}{\phi }_{k+1}(h){\widehat{\mathit{ϵ}}}_{\theta }^{(k)}({\widehat{𝒙}}_{{\lambda }_s},{\lambda }_s)+𝒪(h^{n+2}).$$ (B.4)

Finally, we list the closed-forms of ${\phi }_k$ for $k=1,2,3$:

${\phi }_1(h)$ $={\displaystyle \frac{e^h-1}{h}},$ (B.5)

${\phi }_2(h)$ $={\displaystyle \frac{e^h-h-1}{h^2}},$ (B.6)

${\phi }_3(h)$ $={\displaystyle \frac{e^h-h^2/2-h-1}{h^3}}.$ (B.7)

B.3 Proof of Theorem 3.2 when $k=1$

Proof.

Taking $n=0,t=t_i,s=t_{i-1}$ in Eq. (B.4), we obtain

$${𝒙}_{t_i}=\frac{{\alpha }_{t_i}}{{\alpha }_{t_{i-1}}}{𝒙}_{t_{i-1}}-{\sigma }_t(e^{h_i}-1){\mathit{ϵ}}_{\theta }({𝒙}_{t_{i-1}},t_{i-1})+𝒪(h_i^2).$$ (B.8)

By Assumption B.2 and Eq. (3.7), it holds that

${\widetilde{𝒙}}_{t_i}$ $={\displaystyle \frac{{\alpha }_{t_i}}{{\alpha }_{t_{i-1}}}}{\widetilde{𝒙}}_{t_{i-1}}-{\sigma }_{t_i}(e^{h_i}-1){\mathit{ϵ}}_{\theta }({\widetilde{𝒙}}_{t_{i-1}},t_{i-1})$

$={\displaystyle \frac{{\alpha }_{t_i}}{{\alpha }_{t_{i-1}}}}{\widetilde{𝒙}}_{t_{i-1}}-{\sigma }_{t_i}(e^{h_i}-1)\left({\mathit{ϵ}}_{\theta }({𝒙}_{t_{i-1}},t_{i-1})+𝒪({\widetilde{𝒙}}_{t_{i-1}}-{𝒙}_{t_{i-1}})\right)$

$={\displaystyle \frac{{\alpha }_{t_i}}{{\alpha }_{t_{i-1}}}}{𝒙}_{t_{i-1}}-{\sigma }_{t_i}(e^{h_i}-1){\mathit{ϵ}}_{\theta }({𝒙}_{t_{i-1}},t_{i-1})+𝒪({\widetilde{𝒙}}_{t_{i-1}}-{𝒙}_{t_{i-1}})$

$={𝒙}_{t_i}+𝒪(h_{max}^2)+𝒪({\widetilde{𝒙}}_{t_{i-1}}-{𝒙}_{t_{i-1}}).$

Repeat this argument, we find that

$${\widetilde{𝒙}}_{t_M}={𝒙}_{t_0}+𝒪(M{h}_{max}^2)={𝒙}_{t_0}+𝒪(h_{max}),$$

and thus completes the proof. ∎

B.4 Proof of Theorem 3.2 when $k=2$

We prove the discretization error of the general form of DPM-Solver-2 in Algorithm 4.

Proof.

First, we consider the following update for .

$s_1$ $=t_{\lambda }\left({\lambda }_s+r_{1}h\right),$ (B.9a)

$\overline{𝒖}$ $={\displaystyle \frac{{\alpha }_{s_1}}{{\alpha }_s}}{𝒙}_s-{\sigma }_{s_1}\left(e^{r_{1}h}-1\right){\mathit{ϵ}}_{\theta }({𝒙}_s,s),$ (B.9b)

${\overline{𝒙}}_t$ $={\displaystyle \frac{{\alpha }_t}{{\alpha }_s}}{𝒙}_s-{\sigma }_t\left(e^h-1\right){\mathit{ϵ}}_{\theta }({𝒙}_s,s)-{\displaystyle \frac{{\sigma }_t}{2r_1}}(e^h-1)({\mathit{ϵ}}_{\theta }(\overline{𝒖},s_1)-{\mathit{ϵ}}_{\theta }({𝒙}_s,s)).$ (B.9c)

Note that the above update is the same as a single step of DPM-Solver-2 with $s=t_{i-1}$ and $t=t_i$, except that ${\tilde{x}}_{t_{i-1}}$ is replaced with the exact solution ${𝒙}_{t_{i-1}}$. Once we have proven that ${\overline{𝒙}}_t={𝒙}_t+𝒪(h^3)$, we can show that ${\widetilde{𝒙}}_{t_i}={𝒙}_{t_i}+𝒪(h_{max}^3)+𝒪({\widetilde{𝒙}}_{t_{i-1}}-{𝒙}_{t_{i-1}})$ by a similar argument as in Appendix B.3, and therefore completes the proof.

In this remaining part we prove that ${\overline{𝒙}}_t={𝒙}_t+𝒪(h^3)$.

Taking $n=1$ in Eq. (B.4), we obtain

$${𝒙}_t=\frac{{\alpha }_t}{{\alpha }_s}{𝒙}_s-{\sigma }_{t}h{\phi }_1(h){\mathit{ϵ}}_{\theta }({𝒙}_s,s)-{\sigma }_t{h}^2{\phi }_2(h){\widehat{\mathit{ϵ}}}_{\theta }^{(1)}({\widehat{𝒙}}_{{\lambda }_s},{\lambda }_s)+𝒪(h^3).$$ (B.10)

From Eq. (B.1), we have

${\overline{𝒙}}_t$ $={\displaystyle \frac{{\alpha }_t}{{\alpha }_s}}{𝒙}_s-{\sigma }_t\left(e^h-1\right){\mathit{ϵ}}_{\theta }({𝒙}_s,s)-{\displaystyle \frac{{\sigma }_t}{2r_1}}(e^h-1)({\mathit{ϵ}}_{\theta }(\overline{𝒖},s_1)-{\mathit{ϵ}}_{\theta }({𝒙}_s,s))$

$={\displaystyle \frac{{\alpha }_t}{{\alpha }_s}}{𝒙}_s-{\sigma }_t\left(e^h-1\right){\mathit{ϵ}}_{\theta }({𝒙}_s,s)-{\displaystyle \frac{{\sigma }_t}{2r_1}}\left(e^h-1\right)\left[{\mathit{ϵ}}_{\theta }(\overline{𝒖},s_1)-{\mathit{ϵ}}_{\theta }({𝒙}_{s_1},s_1)\right]$

$-{\displaystyle \frac{{\sigma }_t}{2r_1}}\left(e^h-1\right)\left[({\lambda }_{s_1}-{\lambda }_s){\widehat{\mathit{ϵ}}}_{\theta }^{(1)}({\widehat{𝒙}}_{{\lambda }_s},{\lambda }_s)+𝒪(h^2)\right].$

Note that by the Lipschitzness of ${\mathit{ϵ}}_{\theta }$ w.r.t. $𝒙$ (Assumption B.2),

$$\Vert {\mathit{ϵ}}_{\theta }(\overline{𝒖},s_1)-{\mathit{ϵ}}_{\theta }({𝒙}_{s_1},s_1)\Vert =𝒪(\Vert \overline{𝒖}-{𝒙}_{s_1}\Vert )=𝒪(h^2),$$

where the last equation follows from a similar argument in the proof of $k=1$. Since $e^h-1=𝒪(h)$, the second term of the above display is $𝒪(h^3)$.

As ${\lambda }_{s_1}-{\lambda }_s=r_{1}h$, ${\phi }_i(h)=(e^h-1)/h$ and ${\phi }_2(h)=(e^h-h-1)/h^2$, we find

${𝒙}_t-{\overline{𝒙}}_t$ $={\sigma }_t\left[h^2{\phi }_2(h)-(e^h-1){\displaystyle \frac{{\lambda }_{s_1}-{\lambda }_s}{2r_1}}\right]{\widehat{\mathit{ϵ}}}_{\theta }^{(1)}({\widehat{𝒙}}_{{\lambda }_s},{\lambda }_s)+𝒪(h^3).$

Then, the proof is completed by noticing that

$$h^2{\phi }_2(h)-(e^h-1)\frac{{\lambda }_{s_1}-{\lambda }_s}{2r_1}=(2e^h-h-2-h{e}^h)/2=𝒪(h^3).$$

∎

B.5 Proof of Theorem 3.2 when $k=3$

Proof.

As in Appendix B.4, it suffices to show that the following update has error ${\overline{𝒙}}_t={𝒙}_t+𝒪(h^4)$ for and $h={\lambda }_s-{\lambda }_t$.

$s_1$ $=t_{\lambda }\left({\lambda }_s+r_{1}h\right),s_2=t_{\lambda }\left({\lambda }_s+r_{2}h\right),$ (B.11a)

${\overline{𝒖}}_1$ $={\displaystyle \frac{{\alpha }_{s_1}}{{\alpha }_s}}{𝒙}_s-{\sigma }_{s_1}\left(e^{r_{1}h}-1\right){\mathit{ϵ}}_{\theta }({𝒙}_s,s),$ (B.11b)

${𝑫}_1$ $={\mathit{ϵ}}_{\theta }({\overline{𝒖}}_1,s_1)-{\mathit{ϵ}}_{\theta }({𝒙}_s,s),$ (B.11c)

${\overline{𝒖}}_2$ $={\displaystyle \frac{{\alpha }_{s_2}}{{\alpha }_s}}{𝒙}_s-{\sigma }_{s_2}\left(e^{r_{2}h}-1\right){\mathit{ϵ}}_{\theta }({𝒙}_s,s)-{\displaystyle \frac{{\sigma }_{s_2}{r}_2}{r_1}}\left({\displaystyle \frac{e^{r_{2}h}-1}{r_{2}h}}-1\right){𝑫}_1,$ (B.11d)

${𝑫}_2$ $={\mathit{ϵ}}_{\theta }({\overline{𝒖}}_2,s_2)-{\mathit{ϵ}}_{\theta }({𝒙}_s,s),$ (B.11e)

${\overline{𝒙}}_t$ $={\displaystyle \frac{{\alpha }_t}{{\alpha }_s}}{𝒙}_s-{\sigma }_t\left(e^h-1\right){\mathit{ϵ}}_{\theta }({𝒙}_s,s)-{\displaystyle \frac{{\sigma }_t}{r_2}}\left({\displaystyle \frac{e^h-1}{h}}-1\right){𝑫}_2.$ (B.11f)

First, we prove that

$${\overline{𝒖}}_2={𝒙}_{s_2}+𝒪(h^3).$$ (B.12)

Similar to the proof in Appendix B.4, since $\frac{e^{r_{2}h-1}}{r_{2}h}-1=𝒪(h)$ and ${\overline{𝒖}}_1={𝒙}_{s_1}+𝒪(h^2)$, then

${\overline{𝒖}}_2$ $={\displaystyle \frac{{\alpha }_{s_2}}{{\alpha }_s}}{𝒙}_s-{\sigma }_{s_2}\left(e^{r_{2}h}-1\right){\mathit{ϵ}}_{\theta }({𝒙}_s,s)$

$-{\sigma }_{s_2}{\displaystyle \frac{r_2}{r_1}}\left({\displaystyle \frac{e^{r_{2}h}-1}{r_{2}h}}-1\right)\left({\mathit{ϵ}}_{\theta }({𝒙}_{s_1},s_1)-{\mathit{ϵ}}_{\theta }({𝒙}_s,s)\right)+𝒪(h^3)$

$={\displaystyle \frac{{\alpha }_{s_2}}{{\alpha }_s}}{𝒙}_s-{\sigma }_{s_2}\left(e^{r_{2}h}-1\right){\mathit{ϵ}}_{\theta }({𝒙}_s,s)$

$-{\sigma }_{s_2}{\displaystyle \frac{r_2}{r_1}}\left({\displaystyle \frac{e^{r_{2}h}-1}{r_{2}h}}-1\right){\mathit{ϵ}}_{\theta }^{(1)}({𝒙}_s,s)({\lambda }_{s_1}-{\lambda }_s)+𝒪(h^3).$

Let $h_2=r_{2}h$, then following the same line of arguments in the proof of Appendix B.4, it suffices to check that

${\phi }_1(h_2)h_2$ $=e^{h_2}-1,$

${\phi }_2(h_2)h_2^2$ $={\displaystyle \frac{r_2}{r_1}}\left({\displaystyle \frac{e^{h_2}-1}{h_2}}-1\right)({\lambda }_{s_1}-{\lambda }_s)+𝒪(h^3),$

which holds by applying Taylor expansion.

Using ${\overline{𝒖}}_2={𝒙}_{s_2}+𝒪(h^3)$ and ${\lambda }_{s_2}-{\lambda }_s=r_{2}h=\frac{2}{3}h$, we find that

${\overline{𝒙}}_t$ $={\displaystyle \frac{{\alpha }_t}{{\alpha }_s}}{𝒙}_s-{\sigma }_t\left(e^h-1\right){\mathit{ϵ}}_{\theta }({𝒙}_s,s)-{\sigma }_t{\displaystyle \frac{1}{r_2}}\left({\displaystyle \frac{e^h-1}{h}}-1\right)\left({\mathit{ϵ}}_{\theta }({\overline{𝒖}}_2,s_2)-{\mathit{ϵ}}_{\theta }({𝒙}_s,s)\right)$

$={\displaystyle \frac{{\alpha }_t}{{\alpha }_s}}{𝒙}_s-{\sigma }_t\left(e^h-1\right){\mathit{ϵ}}_{\theta }({𝒙}_s,s)-{\sigma }_t{\displaystyle \frac{1}{r_2}}\left({\displaystyle \frac{e^h-1}{h}}-1\right)\left({\mathit{ϵ}}_{\theta }({𝒙}_{s_2},s_2)-{\mathit{ϵ}}_{\theta }({𝒙}_s,s)\right)+𝒪(h^4)$

$={\displaystyle \frac{{\alpha }_t}{{\alpha }_s}}{𝒙}_s-{\sigma }_t\left(e^h-1\right){\mathit{ϵ}}_{\theta }({𝒙}_s,s)$

$-{\sigma }_t{\displaystyle \frac{1}{r_2}}\left({\displaystyle \frac{e^h-1}{h}}-1\right)\left({\mathit{ϵ}}_{\theta }^{(1)}({𝒙}_s,s)r_{2}h+{\displaystyle \frac{1}{2}}{\mathit{ϵ}}_{\theta }^{(2)}({𝒙}_s,s)r_2^2{h}^2\right)+𝒪(h^4).$

Comparing with the Taylor expansion in Eq. (B.4) with $n=2$:

$${𝒙}_t=\frac{{\alpha }_t}{{\alpha }_s}{𝒙}_s-{\sigma }_{t}h{\phi }_1(h){\mathit{ϵ}}_{\theta }({𝒙}_s,s)-{\sigma }_t{h}^2{\phi }_2(h){\mathit{ϵ}}_{\theta }^{(1)}({𝒙}_s,s)-{\sigma }_t{h}^3{\phi }_3(h){\mathit{ϵ}}_{\theta }^{(2)}({𝒙}_s,s)+𝒪(h^4),$$

we need to check the following conditions:

$h{\phi }_1(h)$ $=e^h-1,$

$h^2{\phi }_2(h)$ $=\left({\displaystyle \frac{e^h-1}{h}}-1\right)h,$

$h^3{\phi }_3(h)$ $=\left({\displaystyle \frac{e^h-1}{h}}-1\right){\displaystyle \frac{r_2{h}^2}{2}}+𝒪(h^4).$

The first two conditions are clear. The last condition follows from

$h^3{\phi }_3(h)$ $=e^h-1-h-{\displaystyle \frac{h^2}{2}}={\displaystyle \frac{h^3}{6}}+𝒪(h^4)=\left({\displaystyle \frac{e^h-1}{h}}-1\right){\displaystyle \frac{r_2{h}^2}{2}}.$

Therefore, ${\overline{𝒙}}_t={𝒙}_t+𝒪(h^4)$. ∎

B.6 Connections to Explicit Exponential Runge-Kutta (expRK) Methods

Assume we have an ODE with the following form:

$$\frac{d{𝒙}_t}{dt}=\alpha {𝒙}_t+𝑵({𝒙}_t,t),$$

where $\alpha \in ℝ$ and $𝑵({𝒙}_t,t)\in {ℝ}^D$ is a non-linear function of ${𝒙}_t$. Given an initial value ${𝒙}_t$ at time $t$, for $h>0$, the true solution at time $t+h$ is

$${𝒙}_{t+h}=e^{\alpha h}{𝒙}_t+e^{\alpha h}{\int }_0^h{e}^{-\alpha \tau }𝑵({𝒙}_{t+\tau },t+\tau )d\tau .$$

The exponential Runge-Kutta methods [25, 31] use some intermediate points to approximate the integral $\int e^{-\alpha \tau }𝑵({𝒙}_{t+\tau },t+\tau )d\tau$. Our proposed DPM-Solver is inspired by the same technique for approximating the same integral with $\alpha =1$ and $𝑵={\widetilde{\mathit{ϵ}}}_{\theta }$. However, DPM-Solver is different from the expRK methods, because their linear term $e^{\alpha h}{𝒙}_t$ is different from our linear term $\frac{{\alpha }_{t+h}}{{\alpha }_t}{𝒙}_t$. In summary, DPM-Solver is inspired by the same technique of expRK for deriving high-order approximations of the exponentially weighted integral, but the formulation of DPM-Solver is different from expRK, and DPM-Solver is customized for the specific formulation of diffusion ODEs.

Appendix C Algorithms of DPM-Solvers

We firstly list the detailed DPM-Solver-1, 2, 3 in Algorithms 3, 4, 5. Note that DPM-Solver-2 is the general case with $r_1\in (0,1)$, and we usually set $r_1=0.5$ for DPM-Solver-2, as in Sec. 3.

Algorithm 3 DPM-Solver-1.

1:initial value ${𝒙}_T$, time steps ${\{t_i\}}_{i=0}^M$, model ${\mathit{ϵ}}_{\theta }$

2:def dpm-solver-1$({\widetilde{𝒙}}_{t_{i-1}},t_{i-1},t_i)$:

3: $h_i\leftarrow {\lambda }_{t_i}-{\lambda }_{t_{i-1}}$

4: ${\widetilde{𝒙}}_{t_i}\leftarrow \frac{{\alpha }_{t_i}}{{\alpha }_{t_{i-1}}}{\widetilde{𝒙}}_{t_{i-1}}-{\sigma }_{t_i}\left(e^{h_i}-1\right){\mathit{ϵ}}_{\theta }({\widetilde{𝒙}}_{t_{i-1}},t_{i-1})$

5: return ${\widetilde{𝒙}}_{t_i}$

6:${\widetilde{𝒙}}_{t_0}\leftarrow {𝒙}_T$

7:for $i\leftarrow 1$ to $M$ do

8: ${\widetilde{𝒙}}_{t_i}\leftarrow \text{dpm-solver-1}({\widetilde{𝒙}}_{t_{i-1}},t_{i-1},t_i)$

9:end for

10:return ${\widetilde{𝒙}}_{t_M}$

Algorithm 4 DPM-Solver-2 (general version).

1:initial value ${𝒙}_T$, time steps ${\{t_i\}}_{i=0}^M$, model ${\mathit{ϵ}}_{\theta }$, $r_1=0.5$

2:def dpm-solver-2$({\widetilde{𝒙}}_{t_{i-1}},t_{i-1},t_i,r_1)$:

3: $h_i\leftarrow {\lambda }_{t_i}-{\lambda }_{t_{i-1}}$

4: $s_i\leftarrow t_{\lambda }\left({\lambda }_{t_{i-1}}+r_1{h}_i\right)$

5: ${𝒖}_i\leftarrow \frac{{\alpha }_{s_i}}{{\alpha }_{t_{i-1}}}{\widetilde{𝒙}}_{t_{i-1}}-{\sigma }_{s_i}\left(e^{r_1{h}_i}-1\right){\mathit{ϵ}}_{\theta }({\widetilde{𝒙}}_{t_{i-1}},t_{i-1})$

6: ${\widetilde{𝒙}}_{t_i}\leftarrow \frac{{\alpha }_{t_i}}{{\alpha }_{t_{i-1}}}{\widetilde{𝒙}}_{t_{i-1}}-{\sigma }_{t_i}(e^{h_i}-1){\mathit{ϵ}}_{\theta }({\widetilde{𝒙}}_{t_{i-1}},t_{i-1})-\frac{{\sigma }_{t_i}}{2r_1}(e^{h_i}-1)({\mathit{ϵ}}_{\theta }({𝒖}_i,s_i)-{\mathit{ϵ}}_{\theta }({\widetilde{𝒙}}_{t_{i-1}},t_{i-1}))$

7: return ${\widetilde{𝒙}}_{t_i}$

8:${\widetilde{𝒙}}_{t_0}\leftarrow {𝒙}_T$

9:for $i\leftarrow 1$ to $M$ do

10: ${\widetilde{𝒙}}_{t_i}\leftarrow \text{dpm-solver-2}({\widetilde{𝒙}}_{t_{i-1}},t_{i-1},t_i,r_1)$

11:end for

12:return ${\widetilde{𝒙}}_{t_M}$

Algorithm 5 DPM-Solver-3.

1:initial value ${𝒙}_T$, time steps ${\{t_i\}}_{i=0}^M$, model ${\mathit{ϵ}}_{\theta }$, $r_1=\frac{1}{3}$, $r_2=\frac{2}{3}$

2:def dpm-solver-3$({\widetilde{𝒙}}_{t_{i-1}},t_{i-1},t_i,r_1,r_2)$:

3: $h_i\leftarrow {\lambda }_{t_i}-{\lambda }_{t_{i-1}}$

4: $s_{2i-1}\leftarrow t_{\lambda }\left({\lambda }_{t_{i-1}}+r_1{h}_i\right),s_{2i}\leftarrow t_{\lambda }\left({\lambda }_{t_{i-1}}+r_2{h}_i\right)$

5: ${𝒖}_{2i-1}\leftarrow \frac{{\alpha }_{s_{2i-1}}}{{\alpha }_{t_{i-1}}}{\widetilde{𝒙}}_{t_{i-1}}-{\sigma }_{s_{2i-1}}\left(e^{r_1{h}_i}-1\right){\mathit{ϵ}}_{\theta }({\widetilde{𝒙}}_{t_{i-1}},t_{i-1})$

6: ${𝑫}_{2i-1}\leftarrow {\mathit{ϵ}}_{\theta }({𝒖}_{2i-1},s_{2i-1})-{\mathit{ϵ}}_{\theta }({\widetilde{𝒙}}_{t_{i-1}},t_{i-1})$

7: ${𝒖}_{2i}\leftarrow \frac{{\alpha }_{s_{2i}}}{{\alpha }_{t_{i-1}}}{\widetilde{𝒙}}_{t_{i-1}}-{\sigma }_{s_{2i}}\left(e^{r_2{h}_i}-1\right){\mathit{ϵ}}_{\theta }({\widetilde{𝒙}}_{t_{i-1}},t_{i-1})-\frac{{\sigma }_{s_{2i}}{r}_2}{r_1}\left(\frac{e^{r_2{h}_i}-1}{r_2{h}_i}-1\right){𝑫}_{2i-1}$

8: ${𝑫}_{2i}\leftarrow {\mathit{ϵ}}_{\theta }({𝒖}_{2i},s_{2i})-{\mathit{ϵ}}_{\theta }({\widetilde{𝒙}}_{t_{i-1}},t_{i-1})$

9: ${\widetilde{𝒙}}_{t_i}\leftarrow \frac{{\alpha }_{t_i}}{{\alpha }_{t_{i-1}}}{\widetilde{𝒙}}_{t_{i-1}}-{\sigma }_{t_i}\left(e^{h_i}-1\right){\mathit{ϵ}}_{\theta }({\widetilde{𝒙}}_{t_{i-1}},t_{i-1})-\frac{{\sigma }_{t_i}}{r_2}\left(\frac{e^{h_i}-1}{h}-1\right){𝑫}_{2i}$

10: return ${\widetilde{𝒙}}_{t_i}$

11:${\widetilde{𝒙}}_{t_0}\leftarrow {𝒙}_T$

12:for $i\leftarrow 1$ to $M$ do

13: ${\widetilde{𝒙}}_{t_i}\leftarrow \text{dpm-solver-3}({\widetilde{𝒙}}_{t_{i-1}},t_{i-1},t_i,r_1,r_2)$

14:end for

15:return ${\widetilde{𝒙}}_{t_M}$

Then we list the adaptive step size algorithms, named as DPM-Solver-12 (combining 1 and 2; Algorithm 6) and DPM-Solver-23 (combining 2 and 3; Algorithm 7). We follow [20] to let the absolute tolerance ${ϵ}_{\text{atol}}=\frac{{𝒙}_{\text{max}}-{𝒙}_{\text{min}}}{256}$ for image data, which is $0.0078$ for VP type DPMs. We can tune the relative tolerance ${ϵ}_{\text{rtol}}$ to balance the accuracy and NFE, and we find that ${ϵ}_{\text{rtol}}=0.05$ is good enough and can converge quickly.

In practice, the inputs of the adaptive step size solvers are batch data. We simply choose $E_2$ and $E_3$ as the maximum value of all the batch data. Besides, we implement the comparison $s>ϵ$ by $|s-ϵ|>{10}^{-5}$ to avoid numerical issues.

Algorithm 6 (DPM-Solver-12) Adaptive step size algorithm by combining DPM-Solver-1 and 2.

1:start time $T$, end time $ϵ$, initial value ${𝒙}_T$, model ${\mathit{ϵ}}_{\theta }$, data dimension $D$, hyperparameters ${ϵ}_{\text{rtol}}=0.05$, ${ϵ}_{\text{atol}}=0.0078$, $h_{\text{init}}=0.05$, $\theta =0.9$

2:the approximated solution ${𝒙}_{ϵ}$ at time $ϵ$

3:$s\leftarrow T$, $h\leftarrow h_{\text{init}}$, $𝒙\leftarrow {𝒙}_T$, ${𝒙}_{\text{prev}}\leftarrow {𝒙}_T$, $r_1\leftarrow \frac{1}{2}$, $\text{NFE}\leftarrow 0$

4:while $s>ϵ$ do

5: $t\leftarrow t_{\lambda }({\lambda }_s+h)$

6: ${𝒙}_1\leftarrow \text{dpm-solver-1}(𝒙,s,t)$

7: ${𝒙}_2\leftarrow \text{dpm-solver-2}(𝒙,s,t,r_1)$ (Share the same function value ${\mathit{ϵ}}_{\theta }(𝒙,s)$ with dpm-solver-1.)

8: $𝜹\leftarrow \max ({ϵ}_{\text{atol}},{ϵ}_{\text{rtol}}\max (|{𝒙}_1|,|{𝒙}_{\text{prev}}|))$

9: $E_2\leftarrow \frac{1}{\sqrt{D}}{\Vert \frac{{𝒙}_1-{𝒙}_2}{𝜹}\Vert }_2$

10: if $E_2\le 1$ then

11: ${𝒙}_{\text{prev}}\leftarrow {𝒙}_1$, $𝒙\leftarrow {𝒙}_2$, $s\leftarrow t$

12: end if

13: $h\leftarrow \min (\theta h{E}_2^{-\frac{1}{2}},{\lambda }_{ϵ}-{\lambda }_s)$

14: $\text{NFE}\leftarrow \text{NFE}+2$

15:end while

16:return $𝒙$, NFE

Algorithm 7 (DPM-Solver-23) Adaptive step size algorithm by combining DPM-Solver-2 and 3.

1:start time $T$, end time $ϵ$, initial value ${𝒙}_T$, model ${\mathit{ϵ}}_{\theta }$, data dimension $D$, hyperparameters ${ϵ}_{\text{rtol}}=0.05$, ${ϵ}_{\text{atol}}=0.0078$, $h_{\text{init}}=0.05$, $\theta =0.9$

2:the approximated solution ${𝒙}_{ϵ}$ at time $ϵ$

3:$s\leftarrow T$, $h\leftarrow h_{\text{init}}$, $𝒙\leftarrow {𝒙}_T$, ${𝒙}_{\text{prev}}\leftarrow {𝒙}_T$, $r_1\leftarrow \frac{1}{3}$, $r_2\leftarrow \frac{2}{3}$, $\text{NFE}\leftarrow 0$

4:while $s>ϵ$ do

5: $t\leftarrow t_{\lambda }({\lambda }_s+h)$

6: ${𝒙}_2\leftarrow \text{dpm-solver-2}(𝒙,s,t,r_1)$

7: ${𝒙}_3\leftarrow \text{dpm-solver-3}(𝒙,s,t,r_1,r_2)$ (Share the same function values with dpm-solver-2.)

8: $𝜹\leftarrow \max ({ϵ}_{\text{atol}},{ϵ}_{\text{rtol}}\max (|{𝒙}_2|,|{𝒙}_{\text{prev}}|))$

9: $E_3\leftarrow \frac{1}{\sqrt{D}}{\Vert \frac{{𝒙}_2-{𝒙}_3}{𝜹}\Vert }_2$

10: if $E_3\le 1$ then

11: ${𝒙}_{\text{prev}}\leftarrow {𝒙}_2$, $𝒙\leftarrow {𝒙}_3$, $s\leftarrow t$

12: end if

13: $h\leftarrow \min (\theta h{E}_3^{-\frac{1}{3}},{\lambda }_{ϵ}-{\lambda }_s)$

14: $\text{NFE}\leftarrow \text{NFE}+3$

15:end while

16:return $𝒙$, NFE

Appendix D Implementation Details of DPM-Solver

D.1 End Time of Sampling

Theoretically, we need to solve diffusion ODEs from time $T$ to time $0$ to generate samples. Practically, the training and evaluation for the noise prediction model ${\mathit{ϵ}}_{\theta }({𝒙}_t,t)$ usually start from time $T$ to time $ϵ$ to avoid numerical issues for $t$ near to $0$, where $ϵ>0$ is a hyperparameter [3].

In contrast to the sampling methods based on diffusion SDEs [2, 3], We don’t add the “denoising” trick at the final step at time $ϵ$ (which is to set the noise variance to zero), and we just solve diffusion ODEs from $T$ to $ϵ$ by DPM-Solver, since we find it performs well enough.

For discrete-time DPMs, we firstly convert the model to continuous time (see Appendix D.2), and then solver it from time $T$ to time $t$.

D.2 Sampling from Discrete-Time DPMs

In this section, we discuss the more general case for discrete-time DPMs, in which we consider the 1000-step DPMs [2] and the 4000-step DPMs [16], and we also consider the end time $ϵ$ for sampling.

Discrete-time DPMs [2] train the noise prediction model at $N$ fixed time steps ${\{t_n\}}_{n=1}^N$. In practice, $N=1000$ or $N=4000$, and the implementation of the 4000-step DPMs [16] converts the time steps of 4000-step DPMs to the range of 1000-step DPMs. Specifically, the noise prediction model is parameterized by ${\widetilde{\mathit{ϵ}}}_{\theta }({𝒙}_n,\frac{1000n}{N})$ for $n=0,\mathrm{\dots },N-1$, where each ${𝒙}_n$ is corresponding to the value at time $t_{n+1}$. In practice, these discrete-time DPMs usually choose uniform time steps between $[0,T]$, thus $t_n=\frac{nT}{N}$, for $n=1,\mathrm{\dots },N$.

However, the discrete-time noise prediction model cannot predict the noise at time less than the smallest time $t_1$. As the smallest time step $t_1=\frac{T}{N}$ and the corresponding discrete-time noise prediction model at time $t_1$ is ${\widetilde{\mathit{ϵ}}}_{\theta }({𝒙}_0,0)$, we need to “scale” the discrete time steps $[t_1,t_N]=[\frac{T}{N},T]$ to the continuous time range $[ϵ,T]$. We propose two types of scaling as following.

Type-1. Scale the discrete time steps $[t_1,t_N]=[\frac{T}{N},T]$ to the continuous time range $[\frac{T}{N},T]$, and let ${\mathit{ϵ}}_{\theta }(\cdot ,t)={\mathit{ϵ}}_{\theta }(\cdot ,\frac{T}{N})$ for $t\in [ϵ,\frac{T}{N}]$. In this case, we can define the continuous-time noise prediction model by

$${\mathit{ϵ}}_{\theta }(𝒙,t)={\widetilde{\mathit{ϵ}}}_{\theta }(𝒙,1000\cdot \max (t-\frac{T}{N},0)),$$ (D.1)

where the continuous time $t\in [ϵ,\frac{T}{N}]$ maps to the discrete input $0$, and the continuous time $T$ maps to the discrete input $\frac{1000(N-1)}{N}$.

Type-2. Scale the discrete time steps $[t_1,t_N]=[\frac{T}{N},T]$ to the continuous time range $[0,T]$. In this case, we can define the continuous-time noise prediction model by

$${\mathit{ϵ}}_{\theta }(𝒙,t)={\widetilde{\mathit{ϵ}}}_{\theta }(𝒙,1000\cdot \frac{(N-1)t}{NT}),$$ (D.2)

where the continuous time $0$ maps to the discrete input $0$, and the continuous time $T$ maps to the discrete input $\frac{1000(N-1)}{N}$.

Note that the input time of ${\widetilde{\mathit{ϵ}}}_{\theta }$ may not be integers, but we find that the noise prediction model can still work well, and we hypothesize that it is because of the smooth time embeddings (e.g., position embeddings [2]). By such reparameterization, the noise prediction model can adopt the continuous-time steps as input, and thus we can also use DPM-Solver for fast sampling.

In practice, we have $T=1$, and the smallest discrete time $t_1={10}^{-3}$. For fixed $K$ number of function evaluations, we empirically find that for small $K$, the Type-1 with $ϵ={10}^{-3}$ may have better sample quality, and for large $K$, the Type-2 with $ϵ={10}^{-4}$ may have better sample quality. We refer to Appendix E for detailed results.

D.3 DPM-Solver in 20 Function Evaluations

Given a fixed budget $K\le 20$ of the number of function evaluations, we uniformly divide the interval $[{\lambda }_T,{\lambda }_{ϵ}]$ into $M=(\lfloor K/3\rfloor +1)$ segments, and take $M$ steps to generate samples. The $M$ steps are dependent on the remainder $R$ of $K$ mod $3$ to make sure the total number of function evaluations is exactly $K$.

• If $R=0$, we firstly take $M-2$ steps of DPM-Solver-3, and then take $1$ step of DPM-Solver-2 and 1 step of DPM-Solver-1. The total number of function evaluations is $3\cdot (\frac{K}{3}-1)+2+1=K$.

• If $R=1$, we firstly take $M-1$ steps of DPM-Solver-3 and then take $1$ step of DPM-Solver-1. The total number of function evaluations is $3\cdot (\frac{K-1}{3})+1=K$.

• If $R=2$, we firstly take $M-1$ steps of DPM-Solver-3 and then take $1$ step of DPM-Solver-2. The total number of function evaluations is $3\cdot (\frac{K-2}{3})+2=K$.

We empirically find that this design of time steps can greatly improve the generation quality, and DPM-Solver can generate comparable samples in 10 steps and high-quality samples in 20 steps.

D.4 Analytical Formulation of the function $t_{\lambda }(\cdot )$ (the inverse function of $\lambda (t)$)

The costs of computing $t_{\lambda }(\cdot )$ is negligible, because for the noise schedules of ${\alpha }_t$ and ${\sigma }_t$ used in previous DPMs (“linear” and “cosine”) [2, 16], both $\lambda (t)$ and its inverse function $t_{\lambda }(\cdot )$ have analytic formulations. We mainly consider the variance preserving type here, since it is the most widely-used type. The functions of other types (variance exploding and sub-variance preserving type) can be similarly derived.

Linear Noise Schedule [2]. We have

$$\log {\alpha }_t=-\frac{({\beta }_1-{\beta }_0)}{4}{t}^2-\frac{{\beta }_0}{2}t,$$

where ${\beta }_0=0.1$ and ${\beta }_1=20$, following [3]. As ${\sigma }_t=\sqrt{1-{\alpha }_t^2}$, we can compute ${\lambda }_t$ analytically. Moreover, the inverse function is

$$t_{\lambda }(\lambda )=\frac{1}{{\beta }_1-{\beta }_0}\left(\sqrt{{\beta }_0^2+2({\beta }_1-{\beta }_0)\log \left(e^{-2\lambda }+1\right)}-{\beta }_0\right).$$

To reduce the influence of numerical issues, we can compute $t_{\lambda }$ by the following equivalent formulation:

$$t_{\lambda }(\lambda )=\frac{2\log \left(e^{-2\lambda }+1\right)}{\sqrt{{\beta }_0^2+2({\beta }_1-{\beta }_0)\log \left(e^{-2\lambda }+1\right)}+{\beta }_0}.$$

And we solve diffusion ODEs between $[ϵ,T]$, where $T=1$.

Cosine Noise Schedule [16]. Denote

$$\log {\alpha }_t=\log \left(\cos \left(\frac{\pi }{2}\cdot \frac{t+s}{1+s}\right)\right)-\log \left(\cos \left(\frac{\pi }{2}\cdot \frac{s}{1+s}\right)\right),$$

where $s=0.008$, following [16]. As [16] clipped the derivatives to ensure the numerical stability, we also clip the maximum time $T=0.9946$. As ${\sigma }_t=\sqrt{1-{\alpha }_t^2}$, we can compute ${\lambda }_t$ analytically. Moreover, given a fixed $\lambda$, let

$$f(\lambda )=-\frac{1}{2}\log \left(e^{-2\lambda }+1\right),$$

which computes the corresponding $\log \alpha$ for $\lambda$. Then the inverse function is

$$t_{\lambda }(\lambda )=\frac{2(1+s)}{\pi }\arccos \left(e^{f(\lambda )+\log \cos \left(\frac{\pi s}{2(1+s)}\right)}\right)-s.$$

And we solve diffusion ODEs between $[ϵ,T]$, where $T=0.9946$.

D.5 Conditional Sampling by DPM-Solver

DPM-Solver can also be used for conditional sampling, with a simple modification. The conditional generation needs to sample from the conditional diffusion ODE [3, 4] which includes the conditional noise prediction model. We follow the classifier guidance method [4] to define the conditional noise prediction model as ${\mathit{ϵ}}_{\theta }({𝒙}_t,t,y)≔{\mathit{ϵ}}_{\theta }({𝒙}_t,t)-s\cdot {\sigma }_t{\nabla }_{𝒙}\log p_t(y|{𝒙}_t;\theta )$, where $p_t(y|{𝒙}_t;\theta )$ is a pre-trained classifier and $s$ is the classifier guidance scale (default is 1.0). Thus, we can use DPM-Solver to solve this diffusion ODE for fast conditional sampling, as shown in Fig. 1.

D.6 Numerical Stability

As we need to compute $e^{h_i}-1$ in the algorithm of DPM-Solver, we follow [10] to use expm1($h_i$) instead of exp($h_i$)-1 to improve numerical stability.

Appendix E Experiment Details

We test our method for sampling the most widely-used variance-preserving (VP) type DPMs [1, 2]. In this case, we have ${\alpha }_t^2+{\sigma }_t^2=1$ for all $t\in [0,T]$ and $\tilde{\sigma }=1$. In spite of this, our method and theoretical results are general and independent of the choice of the noise schedule ${\alpha }_t$ and ${\sigma }_t$.

For all experiments, we evaluate DPM-Solver on NVIDIA A40 GPUs. However, the computation resource can be other types of GPU, such as NVIDIA GeForce RTX 2080Ti, because we can tune the batch size for sampling.

E.1 Diffusion ODEs w.r.t. $\lambda$

Alternatively, the diffusion ODE can be reparameterized to the $\lambda$ domain. In this section, we propose the formulation of diffusion ODEs w.r.t. $\lambda$ for VP type, and other types can be similarly derived.

For a given $\lambda$, denote ${\widehat{\alpha }}_{\lambda }≔{\alpha }_{t(\lambda )}$, ${\widehat{\sigma }}_{\lambda }≔{\sigma }_{t(\lambda )}$. As ${\widehat{\alpha }}_{\lambda }^2+{\widehat{\sigma }}_{\lambda }^2=1$, we can prove that $\frac{d\lambda }{d{\widehat{\alpha }}_{\lambda }}=\frac{1}{{\widehat{\alpha }}_{\lambda }{\widehat{\sigma }}_{\lambda }^2}$, so $\frac{d\log {\widehat{\alpha }}_{\lambda }}{d\lambda }={\widehat{\sigma }}_{\lambda }^2$. Applying change-of-variable to Eq. (2.7), we have

$$\frac{d{\widehat{𝒙}}_{\lambda }}{d\lambda }={\widehat{𝒉}}_{\theta }({\widehat{𝒙}}_{\lambda },\lambda )≔{\widehat{\sigma }}_{\lambda }^2{\widehat{𝒙}}_{\lambda }-{\widehat{\sigma }}_{\lambda }{\widehat{\mathit{ϵ}}}_{\theta }({\widehat{𝒙}}_{\lambda },\lambda ).$$ (E.1)

The ODE Eq. (E.1) can be also solved directly by RK methods, and we use such formulation for the experiments of RK2 ($\lambda$) and RK3 ($\lambda$) in Table 1.

E.2 Code Implementation

We implement our code with both JAX (for continuous-time DPMs) and PyTorch (for discrete-time DPMs), and our code is released at https://github.com/LuChengTHU/dpm-solver.

E.3 Sample Quality Comparison with Continuous-Time Sampling Methods

Table 3: Sample quality measured by FID $\downarrow$ on CIFAR-10 dataset with continuous-time methods, varying the number of function evaluations (NFE).

Sampling method $\backslash$ NFE			10	12	15	20	50	200	1000
CIFAR-10 (continuous-time model (VP deep) [3], linear noise schedule)
SDE	Euler (denoise) [3]	$ϵ={10}^{-3}$	304.73	278.87	248.13	193.94	66.32	12.27	2.44
		$ϵ={10}^{-4}$	444.63	427.54	395.95	300.41	101.66	22.98	5.01
	Improved Euler [20]	$ϵ={10}^{-3}$	82.42(NFE=48), 2.73(NFE=151), 2.44(NFE=180)
ODE	RK45 Solver [28, 3]	$ϵ={10}^{-3}$	19.55(NFE=26), 17.81(NFE=38), 3.55(NFE=62)
		$ϵ={10}^{-4}$	51.66(NFE=26), 21.54(NFE=38), 12.72(NFE=50), 2.61(NFE=62)
	DPM-Solver (ours)	$ϵ={10}^{-3}$	4.70	3.75	3.24	3.99	3.84 (NFE = 42)
		$ϵ={10}^{-4}$	6.96	4.93	3.35	2.87	2.59 (NFE = 51)

Table 3 shows the detailed FID results, which is corresponding to Fig. 2(a). We use the official code and checkpoint in [3], the code license is Apache License 2.0. We use their released “checkpoint_8” of the “VP deep” type. We compare methods for $ϵ={10}^{-3}$ and $ϵ={10}^{-4}$. We find that the sampling methods based on diffusion SDEs can achieve better sample quality with $ϵ={10}^{-3}$; and that the sampling methods based on diffusion ODEs can achieve better sample quality with $ϵ={10}^{-4}$. For DPM-Solver, we find that DPM-Solver with less than 15 NFE can achieve better FID with $ϵ={10}^{-3}$ than $ϵ={10}^{-4}$, while DPM-Solver with more than 15 NFE can achieve better FID with $ϵ={10}^{-4}$ than $ϵ={10}^{-3}$.

For the diffusion SDEs with Euler discretization, we use the PC sampler in [3] with “euler_maruyama” predictor and no corrector, which uses uniform time steps between $T$ and $ϵ$. We add the “denoise” trick at the final step, which can greatly improve the FID score for $ϵ={10}^{-3}$.

For the diffusion SDEs with Improved Euler discretization [20], we follow the results in their original paper, which only includes the results with $ϵ={10}^{-3}$. The corresponding relative tolerance ${ϵ}_{rel}$ are $0.50$, $0.10$ and $0.05$, respectively.

For the diffusion ODEs with RK45 Solver, we use the code in [3], and tune the atol and rtol of the solver. For the NFE from small to large, we use the same atol = rtol = $0.1$, $0.01$, $0.001$ for the results of $ϵ={10}^{-3}$, and the same atol = rtol = $0.1$, $0.05$, $0.02$, $0.01$, $0.001$ for the results of $ϵ={10}^{-4}$, respectively.

For the diffusion ODEs with DPM-Solver, we use the method in Appendix D.3 for NFE $\le 20$, and the adaptive step size solver in Appendix C. For $ϵ={10}^{-3}$, we use DPM-Solver-12 with relative tolerance ${ϵ}_{\text{rtol}}=0.05$. For $ϵ={10}^{-4}$, we use DPM-Solver-23 with relative tolerance ${ϵ}_{\text{rtol}}=0.05$.

E.4 Sample Quality Comparison with RK Methods

Table 1 shows the different performance of RK methods and DPM-Solver-2 and 3. We list the detailed settings in this section.

Assume we have an ODE with

$$\frac{\mathrm{d}{𝒙}_t}{\mathrm{d}t}=𝑭({𝒙}_t,t),$$

Starting with ${\widetilde{𝒙}}_{t_{i-1}}$ at time $t_{i-1}$, we use RK2 to approximate the solution ${\widetilde{𝒙}}_{t_i}$ at time $t_i$ in the following formulation (which is known as the explicit midpoint method):

$h_i$ $=t_i-t_{i-1},$

$s_i$ $=t_{i-1}+{\displaystyle \frac{1}{2}}{h}_i,$

${𝒖}_i$ $={\widetilde{𝒙}}_{t_{i-1}}+{\displaystyle \frac{h_i}{2}}𝑭({\widetilde{𝒙}}_{t_{i-1}},t_{i-1}),$

${\widetilde{𝒙}}_{t_i}$ $={\widetilde{𝒙}}_{t_{i-1}}+h_i𝑭({𝒖}_i,s_i).$

And we use the following RK3 to approximate the solution ${\widetilde{𝒙}}_{t_i}$ at time $t_i$ (which is known as “Heun’s third-order method”), because it is very similar to our proposed DPM-Solver-3:

$h_i$ $=t_i-t_{i-1},r_1={\displaystyle \frac{1}{3}},r_2={\displaystyle \frac{2}{3}},$

$s_{2i-1}$ $=t_{i-1}+r_1{h}_i,s_{2i}=t_{i-1}+r_2{h}_i,$

${𝒖}_{2i-1}$ $={\widetilde{𝒙}}_{t_{i-1}}+r_1{h}_i𝑭({\widetilde{𝒙}}_{t_{i-1}},t_{i-1}),$

${𝒖}_{2i}$ $={\widetilde{𝒙}}_{t_{i-1}}+r_2{h}_i𝑭({𝒖}_{2i-1},s_{2i-1}),$

${\widetilde{𝒙}}_{t_i}$ $={\widetilde{𝒙}}_{t_{i-1}}+{\displaystyle \frac{h_i}{4}}𝑭({\widetilde{𝒙}}_{t_{i-1}},t_{i-1})+{\displaystyle \frac{3h_i}{4}}𝑭({𝒖}_{2i},s_{2i}).$

We use $𝑭({𝒙}_t,t)={𝒉}_{\theta }({𝒙}_t,t)$ in Eq. (2.7) for the results with RK2 ($t$) and RK3 ($t$), and $𝑭({\widehat{𝒙}}_{\lambda },\lambda )={\widehat{𝒉}}_{\theta }({\widehat{𝒙}}_{\lambda },\lambda )$ in Eq. (E.1) for the results with RK2 ($\lambda$) and RK3 ($\lambda$). For all experiments, we use the uniform step size w.r.t. $t$ or $\lambda$.

E.5 Sample Quality Comparison with Discrete-Time Sampling Methods

Table 4: Sample quality measured by FID $\downarrow$ on CIFAR-10, CelebA 64$\times$64 and ImageNet 64$\times$64 with discrete-time DPMs, varying the number of function evaluations (NFE). The method ^†GGDM needs extra training, and some results are missing in their original papers, which are replaced by “$\backslash$”.

Sampling method $\backslash$ NFE		10	12	15	20	50	200	1000
CIFAR-10 (discrete-time model [2], linear noise schedule)
DDPM [2]	Discrete	278.67	246.29	197.63	137.34	32.63	4.03	3.16
Analytic-DDPM [21]	Discrete	35.03	27.69	20.82	15.35	7.34	4.11	3.84
Analytic-DDIM [21]	Discrete	14.74	11.68	9.16	7.20	4.28	3.60	3.86
†GGDM [18]	Discrete	8.23	$\backslash$	6.12	4.72	$\backslash$	$\backslash$	$\backslash$
DDIM [19]	Discrete	13.58	11.02	8.92	6.94	4.73	4.07	3.95
DPM-Solver (Type-1 discrete)	$ϵ={10}^{-3}$	6.37	4.65	3.78	4.28	3.90 (NFE = 44)
	$ϵ={10}^{-4}$	11.32	7.31	4.75	3.80	3.57 (NFE = 46)
DPM-Solver (Type-2 discrete)	$ϵ={10}^{-3}$	6.42	4.86	4.39	5.52	5.22 (NFE = 42)
	$ϵ={10}^{-4}$	10.16	6.26	4.17	3.72	3.48 (NFE = 44)
CelebA 64$\times$64 (discrete-time model [19], linear noise schedule)
DDPM [2]	Discrete	310.22	277.16	207.97	120.44	29.25	3.90	3.50
Analytic-DDPM [21]	Discrete	28.99	25.27	21.80	18.14	11.23	6.51	5.21
Analytic-DDIM [21]	Discrete	15.62	13.90	12.29	10.45	6.13	3.46	3.13
DDIM [19]	Discrete	10.85	9.99	7.78	6.64	5.23	4.78	4.88
DPM-Solver (Type-1 discrete)	$ϵ={10}^{-3}$	7.15	5.51	4.28	4.40	4.23 (NFE = 36)
	$ϵ={10}^{-4}$	6.92	4.20	3.05	2.82	2.71 (NFE = 36)
DPM-Solver (Type-2 discrete)	$ϵ={10}^{-3}$	7.33	6.23	5.85	6.87	6.68 (NFE = 36)
	$ϵ={10}^{-4}$	5.83	3.71	3.11	3.13	3.10 (NFE = 36)
ImageNet 64$\times$64 (discrete-time model [16], cosine noise schedule)
DDPM [2]	Discrete	305.43	287.66	256.69	209.73	83.86	28.39	17.58
Analytic-DDPM [21]	Discrete	60.65	53.66	45.98	37.67	22.45	17.16	16.14
Analytic-DDIM [21]	Discrete	70.62	54.88	41.56	30.88	19.23	17.49	17.57
†GGDM [18]	Discrete	37.32	$\backslash$	24.69	20.69	$\backslash$	$\backslash$	$\backslash$
DDIM [19]	Discrete	67.07	52.69	40.49	30.67	20.10	17.84	17.73
DPM-Solver (Type-1 discrete)	$ϵ={10}^{-3}$	24.44	20.03	19.31	18.59	17.50 (NFE = 48)
	$ϵ={10}^{-4}$	27.74	23.66	20.09	19.06	17.56 (NFE = 51)
DPM-Solver (Type-2 discrete)	$ϵ={10}^{-3}$	24.40	19.97	19.23	18.53	17.47 (NFE = 57)
	$ϵ={10}^{-4}$	27.72	23.75	20.02	19.08	17.62 (NFE = 48)

Table 5: Sample quality measured by FID $\downarrow$ on ImageNet 128$\times$128 with classifier guidance and on LSUN bedroom 256$\times$256, varying the number of function evaluations (NFE). For DDIM and DDPM, we use uniform time steps for all the experiments, except that the experiment^† uses the fine-tuned time steps by [4]. For DPM-Solver, we use the uniform logSNR steps as described in Appendix D.3.

Sampling method $\backslash$ NFE		10	12	15	20	50	100	250
ImageNet 128$\times$128 (discrete-time model [4], linear noise schedule, classifier guidance scale: 1.25)
DDPM [2]	Discrete	199.56	172.09	146.42	119.13	49.38	23.27	2.97
DDIM [19]	Discrete	11.12	9.38	8.22	7.15	5.05	4.18	3.54
DPM-Solver (Type-1 discrete)	$ϵ={10}^{-3}$	7.32	4.08	3.60	3.89	3.63	3.62	3.63
	$ϵ={10}^{-4}$	13.91	5.84	4.00	3.52	3.13	3.10	3.09
LSUN bedroom 256$\times$256 (discrete-time model [4], linear noise schedule)
DDPM [2]	Discrete	274.67	251.26	224.88	190.14	82.70	34.89	†2.02
DDIM [19]	Discrete	10.05	7.51	5.90	4.98	2.92	2.30	2.02
DPM-Solver (Type-1 discrete)	$ϵ={10}^{-3}$	6.10	4.29	3.30	3.09	2.53	2.46	2.46
	$ϵ={10}^{-4}$	8.04	4.21	2.94	2.60	2.01	1.95	1.94

We compare DPM-Solver with other discrete-time sampling methods for DPMs, as shown in Table 4 and Table 5. We use the code in [19] for sampling with DDPM and DDIM, and the code license is MIT License. We use the code in [21] for sampling with Analytic-DDPM and Analytic-DDIM, whose license is unknown. We directly follow the best results in the original paper of GGDM [18].

For the CIFAR-10 experiments, we use the pretrained checkpoint by [2], which is also provided in the released code in [19]. We use quadratic time steps for DDPM and DDIM, which empirically has better FID performance than the uniform time steps [19]. We use the uniform time steps for Analytic-DDPM and Analytic-DDIM. For DPM-Solver, we use both Type-1 discrete and Type-2 discrete methods to convert the discrete-time model to the continuous-time model. We use the method in Appendix D.3 for NFE $\le 20$, and the adaptive step size solver in Appendix C for NFE $>20$. For all the experiments, we use DPM-Solver-12 with relative tolerance ${ϵ}_{\text{rtol}}=0.05$.

For the CelebA 64x64 experiments, we use the pretrained checkpoint by [19]. We use quadratic time steps for DDPM and DDIM, which empirically has better FID performance than the uniform time steps [19]. We use the uniform time steps for Analytic-DDPM and Analytic-DDIM. For DPM-Solver, we use both Type-1 discrete and Type-2 discrete methods to convert the discrete-time model to the continuous-time model. We use the method in Appendix D.3 for NFE $\le 20$, and the adaptive step size solver in Appendix C for NFE $>20$. For all the experiments, we use DPM-Solver-12 with relative tolerance ${ϵ}_{\text{rtol}}=0.05$. Note that our best FID results on CelebA 64x64 is even better than the 1000-step DDPM (and all the other methods).

For the ImageNet 64x64 experiments, we use the pretrained checkpoint by [16], and the code license is MIT License. We use the uniform time steps for DDPM and DDIM, following [19]. We use the uniform time steps for Analytic-DDPM and Analytic-DDIM. For DPM-Solver, we use both Type-1 discrete and Type-2 discrete methods to convert the discrete-time model to the continuous-time model. We use the method in Appendix D.3 for NFE $\le 20$, and the adaptive step size solver in Appendix C for NFE $>20$. For all the experiments, we use DPM-Solver-23 with relative tolerance ${ϵ}_{\text{rtol}}=0.05$. Note that the ImageNet dataset includes real human photos and it may have privacy issues, as discussed in [42].

For the ImageNet 128x128 experiments, we use classifier guidance for sampling with the pretrained checkpoints (for both the diffusion model and the classifier model) by [4], and the code license is MIT License. We use the uniform time steps for DDPM and DDIM, following [19]. For DPM-Solver, we only use Type-1 discrete method to convert the discrete-time model to the continuous-time model. We use the method in Appendix D.3 for NFE $\le 20$, and the adaptive step size solver DPM-Solver-12 with relative tolerance ${ϵ}_{\text{rtol}}=0.05$ (detailed in Appendix C) for NFE $>20$. For all the experiments, we set the classifier guidance scale $s=1.25$, which is the best setting for DDIM in [4] (we refer to their Table 14 for details).

For the LSUN bedroom 256x256 experiments, we use the unconditional pretrained checkpoint by [4], and the code license is MIT License. We use the uniform time steps for DDPM and DDIM, following [19]. For DPM-Solver, we only use Type-1 discrete method to convert the discrete-time model to the continuous-time model. We use the method in Appendix D.3 for DPM-Solver.

E.6 Comparing Different Orders of DPM-Solver

We also compare the sample quality of the different orders of DPM-Solver, as shown in Table 6. We use DPM-Solver-1,2,3 with uniform time steps w.r.t. $\lambda$, and the fast version in Appendix D.3 for NFE less than 20, and we name it as DPM-Solver-fast. For the discrete-time models, we only compare the Type-2 discrete method, and the results of Type-1 are similar.

As the actual NFE of DPM-Solver-2 is $2\times \lfloor \text{NFE}/2\rfloor$ and the actual NFE of DPM-Solver-3 is $3\times \lfloor \text{NFE}/3\rfloor$, which may be smaller than NFE, we use the notation ^† to note that the actual NFE is less than the given NFE. We find that for NFE less than 20, the proposed fast version (DPM-Solver-fast) is usually better than the single order method, and for larger NFE, DPM-Solver-3 is better than DPM-Solver-2, and DPM-Solver-2 is better than DPM-Solver-1, which matches our proposed convergence rate analysis.

Table 6: Sample quality measured by FID $\downarrow$ of different orders of DPM-Solver, varying the number of function evaluations (NFE). The results with ^† means the actual NFE is smaller than the given NFE because the given NFE cannot be divided by $2$ or $3$. For DPM-Solver-fast, we only evaluate it for NFE less than 20, because it is almost the same as DPM-Solver-3 for large NFE.

Sampling method $\backslash$ NFE		10	12	15	20	50	200	1000
CIFAR-10 (VP deep continuous-time model [3])
$ϵ={10}^{-3}$	DPM-Solver-1	11.83	9.69	7.78	6.17	4.28	3.85	3.83
	DPM-Solver-2	5.94	4.88	†4.30	3.94	3.78	3.74	3.74
	DPM-Solver-3	†18.37	5.53	4.08	†4.04	†3.81	†3.78	†3.78
	DPM-Solver-fast	4.70	3.75	3.24	3.99	$\backslash$	$\backslash$	$\backslash$
$ϵ={10}^{-4}$	DPM-Solver-1	11.29	9.07	7.15	5.50	3.32	2.72	2.64
	DPM-Solver-2	7.30	5.28	†4.23	3.26	2.69	2.60	2.59
	DPM-Solver-3	†54.56	6.03	3.55	†2.90	†2.65	†2.62	†2.62
	DPM-Solver-fast	6.96	4.93	3.35	2.87	$\backslash$	$\backslash$	$\backslash$
CIFAR-10 (DDPM discrete-time model [2]), DPM-Solver with Type-2 discrete
$ϵ={10}^{-3}$	DPM-Solver-1	16.69	13.63	11.08	8.90	6.24	5.44	5.29
	DPM-Solver-2	7.90	6.15	†5.53	5.24	5.23	5.25	5.25
	DPM-Solver-3	†24.37	8.20	5.73	†5.43	†5.29	†5.25	†5.25
	DPM-Solver-fast	6.42	4.86	4.39	5.52	$\backslash$	$\backslash$	$\backslash$
$ϵ={10}^{-4}$	DPM-Solver-1	13.61	10.98	8.71	6.79	4.36	3.63	3.49
	DPM-Solver-2	11.80	6.31	†5.23	3.95	3.50	3.46	3.46
	DPM-Solver-3	†67.02	9.45	5.21	†3.81	†3.49	†3.45	†3.45
	DPM-Solver-fast	10.16	6.26	4.17	3.72	$\backslash$	$\backslash$	$\backslash$
CelebA 64$\times$64 (discrete-time model [19], linear noise schedule), DPM-Solver with Type-2 discrete
$ϵ={10}^{-3}$	DPM-Solver-1	18.66	16.30	13.92	11.84	8.85	7.24	6.93
	DPM-Solver-2	5.89	5.83	†6.08	6.38	6.78	6.84	6.85
	DPM-Solver-3	†11.45	5.46	6.18	†6.51	†6.87	†6.84	†6.85
	DPM-Solver-fast	7.33	6.23	5.85	6.87	$\backslash$	$\backslash$	$\backslash$
$ϵ={10}^{-4}$	DPM-Solver-1	13.24	11.13	9.08	7.24	4.50	3.48	3.25
	DPM-Solver-2	4.28	3.40	†3.30	3.17	3.19	3.20	3.20
	DPM-Solver-3	†49.48	3.84	3.09	†3.15	†3.20	†3.20	†3.20
	DPM-Solver-fast	5.83	3.71	3.11	3.13	$\backslash$	$\backslash$	$\backslash$
ImageNet 64$\times$64 (discrete-time model [16], cosine noise schedule), DPM-Solver with Type-2 discrete
$ϵ={10}^{-3}$	DPM-Solver-1	32.84	28.54	24.79	21.71	18.30	17.45	17.18
	DPM-Solver-2	29.20	24.97	†22.26	19.94	17.79	17.29	17.27
	DPM-Solver-3	†57.48	24.62	19.76	†18.95	†17.52	17.26	17.27
	DPM-Solver-fast	24.40	19.97	19.23	18.53	$\backslash$	$\backslash$	$\backslash$
$ϵ={10}^{-4}$	DPM-Solver-1	32.31	28.44	25.15	22.38	19.14	17.95	17.44
	DPM-Solver-2	33.16	27.28	†24.26	20.58	18.04	17.46	17.41
	DPM-Solver-3	†162.27	27.28	22.38	†19.39	†17.71	†17.43	†17.41
	DPM-Solver-fast	27.72	23.75	20.02	19.08	$\backslash$	$\backslash$	$\backslash$

E.7 Runtime Comparison between DPM-Solver and DDIM

Theoretically, for the same NFE, the runtime of DPM-Solver and DDIM are almost the same (linear to NFE) because the main computation costs are the serial evaluations of the large neural network ${\mathit{ϵ}}_{\theta }$ and the other coefficients are analytically computed with ignorable costs.

Table 7 shows the runtime of DPM-Solver and DDIM on a single NVIDIA A40, varying different datasets and NFE. We use torch.cuda.Event and torch.cuda.synchronize for accurately computing the runtime. We use the discrete-time pretrained diffusion models for each dataset. We evaluate the runtime for 8 batches and computes the mean and std of the runtime. We use 64 batch size for LSUN bedroom 256x256 due to the GPU memory limitation, and 128 batch size for other datasets.

For DDIM, we use the official implementation. We find that our implementation of DPM-Solver reduces some repetitive computation of the coefficients, so under the same NFE, DPM-Solver is slightly faster than DDIM of their implementation. Nevertheless, the runtime evaluation results show that the runtime of DPM-Solver and DDIM are almost the same for the same NFE, and the runtime is approximately linear to the NFE. Therefore, the speedup for the NFE is almost the actual speedup of the runtime, so the proposed DPM-Solver can greatly speedup the sampling of DPMs.

Table 7: Runtime of a single batch (second / batch, $\pm$std) on a single NVIDIA A40 of DDIM and DPM-Solver for sampling by discrete-time diffusion models, varying the number of function evaluations (NFE).

Sampling method $\backslash$ NFE	10	20	50	100
CIFAR-10 32$\times$32 (batch size = 128)
DDIM	0.956($\pm$0.011)	1.924($\pm$0.016)	4.838($\pm$0.024)	9.668($\pm$0.013)
DPM-Solver	0.923($\pm$0.006)	1.833($\pm$0.004)	4.580($\pm$0.005)	9.204($\pm$0.011)
CelebA 64$\times$64 (batch size = 128)
DDIM	3.253($\pm$0.015)	6.438($\pm$0.029)	16.132($\pm$0.050)	32.255($\pm$0.044)
DPM-Solver	3.126($\pm$0.003)	6.272($\pm$0.006)	15.676($\pm$0.008)	31.269($\pm$0.012)
ImageNet 64$\times$64 (batch size = 128)
DDIM	5.084($\pm$0.018)	10.194($\pm$0.022)	25.440($\pm$0.044)	50.926($\pm$0.042)
DPM-Solver	4.992($\pm$0.004)	9.991($\pm$0.003)	24.948($\pm$0.007)	49.835($\pm$0.028)
ImageNet 128$\times$128 (batch size = 128, with classifier guidance)
DDIM	29.082($\pm$0.015)	58.159($\pm$0.012)	145.427($\pm$0.011)	290.874($\pm$0.134)
DPM-Solver	28.865($\pm$0.011)	57.645($\pm$0.008)	144.124($\pm$0.035)	288.157($\pm$0.022)
LSUN bedroom 256$\times$256 (batch size = 64)
DDIM	37.700($\pm$0.005)	75.316($\pm$0.013)	188.275($\pm$0.172)	378.790($\pm$0.105)
DPM-Solver	36.996($\pm$0.039)	73.873($\pm$0.023)	184.590($\pm$0.010)	369.090($\pm$0.076)

E.8 Conditional Sampling on ImageNet 256x256

For the conditional sampling in Fig. 1, we use the pretrained checkpoint in [4] with classifier guidance (ADM-G), and the classifier scale is $1.0$. The code license is MIT License. We use uniform time step for DDIM, and the fast version for DPM-Solver in Appendix D.3 (DPM-Solver-fast) with 10, 15, 20 and 100 steps.

Fig. 3 shows the conditional sample results by DDIM and DPM-Solver. We find that DPM-Solver with 15 NFE can generate comparable samples with DDIM with 100 NFE.

NFE = 10

NFE = 15

NFE = 20

NFE = 100

DDIM

[19]

DPM-Solver

(ours)

DDIM

[19]

DPM-Solver

(ours)

Figure 3: Samples by DDIM [19] and DPM-Solver (ours) with 10, 15, 20, 100 number of function evaluations (NFE) with the same random seed, using the pre-trained DPMs on ImageNet 256$\times$256 with classifier guidance [4].

E.9 Additional Samples

Additional sampling results on CIFAR-10, CelebA 64x64, ImageNet 64x64, LSUN bedroom 256x256 [40], ImageNet 256x256 are reported in Figs. 4-8.

NFE = 10

NFE = 12

NFE = 15

NFE = 20

DDIM
[19]
DPM-Solver
(ours)

Figure 4: Random samples by DDIM [19] (quadratic time steps) and DPM-Solver (ours) with 10, 12, 15, 20 number of function evaluations (NFE) with the same random seed, using the pre-trained discrete-time DPMs [2] on CIFAR-10.

NFE = 10

NFE = 12

NFE = 15

NFE = 20

DDIM
[19]
DPM-Solver
(ours)

Figure 5: Random samples by DDIM [19] (quadratic time steps) and DPM-Solver (ours) with 10, 12, 15, 20 number of function evaluations (NFE) with the same random seed, using the pre-trained discrete-time DPMs [19] on CelebA 64x64.

NFE = 10

NFE = 12

NFE = 15

NFE = 20

DDIM
[19]
DPM-Solver
(ours)

Figure 6: Random samples by DDIM [19] (uniform time steps) and DPM-Solver (ours) with 10, 12, 15, 20 number of function evaluations (NFE) with the same random seed, using the pre-trained discrete-time DPMs [16] on ImageNet 64x64.

NFE = 10

NFE = 12

NFE = 15

NFE = 20

DDIM
[19]
DPM-Solver
(ours)

Figure 7: Random samples by DDIM [19] (uniform time steps) and DPM-Solver (ours) with 10, 12, 15, 20 number of function evaluations (NFE) with the same random seed, using the pre-trained discrete-time DPMs [4] on LSUN bedroom 256x256.

NFE = 10

NFE = 12

NFE = 15

NFE = 20

DDIM
[19]
DPM-Solver
(ours)

Figure 8: Random class-conditional samples (class: 90, lorikeet) by DDIM [19] (uniform time steps) and DPM-Solver (ours) with 10, 12, 15, 20 number of function evaluations (NFE) with the same random seed, using the pre-trained discrete-time DPMs [4] on ImageNet 256x256 with classifier guidance (classifier scale: 1.0).