DPM-Solver-v3: Improved Diffusion ODE Solver with Empirical Model Statistics

Suppose we have a $D$-dimensional data distribution $q_{0}(\bm{x}_{0})$. Diffusion probabilistic models (DPMs) [47, 15, 51] define a forward diffusion process $\{q_{t}\}_{t=0}^{T}$ to gradually degenerate the data $\bm{x}_{0}\sim q_{0}(\bm{x}_{0})$ with Gaussian noise, which satisfies the transition kernel $q_{0t}(\bm{x}_{t}|\bm{x}_{0})=\mathcal{N}(\alpha_{t}\bm{x}_{0},\sigma_{t}^{2}\bm{I})$, such that $q_{T}(\bm{x}_{T})$ is approximately pure Gaussian. $\alpha_{t},\sigma_{t}$ are smooth scalar functions of $t$, which are called noise schedule. The transition can be easily applied by $\bm{x}_{t}=\alpha_{t}\bm{x}_{0}+\sigma_{t}\bm{\epsilon},\bm{\epsilon}\sim\mathcal{N}(\bm{0},\bm{I})$. To train DPMs, a neural network $\bm{\epsilon}_{\theta}(\bm{x}_{t},t)$ is usually parameterized to predict the noise $\bm{\epsilon}$ by minimizing $\mathbb{E}_{\bm{x}_{0}\sim q_{0}(\bm{x}_{0}),\bm{\epsilon}\sim\mathcal{N}(\bm{0},\bm{I}),t\sim\mathcal{U}(0,T)}\left[w(t)\|\bm{\epsilon}_{\theta}(\bm{x}_{t},t)-\bm{\epsilon}\|_{2}^{2}\right]$, where $w(t)$ is a weighting function. Sampling of DPMs can be performed by solving diffusion ODE [51] from time $T$ to time $0$:

$$\frac{\mathrm{d}\bm{x}_{t}}{\mathrm{d}t}=f(t)\bm{x}_{t}+\frac{g^{2}(t)}{2\sigma_{t}}\bm{\epsilon}_{\theta}(\bm{x}_{t},t),$$

(1)

where $f(t)=\frac{\mathrm{d}\log\alpha_{t}}{\mathrm{d}t}$, $g^{2}(t)=\frac{\mathrm{d}\sigma_{t}^{2}}{\mathrm{d}t}-2\frac{\mathrm{d}\log\alpha_{t}}{\mathrm{d}t}\sigma_{t}^{2}$ [23]. In addition, the conditional sampling by DPMs can be conducted by guided sampling [10, 16] with a conditional noise predictor $\bm{\epsilon}_{\theta}(\bm{x}_{t},t,c)$, where $c$ is the condition. Specifically, classifier-free guidance [16] combines the unconditional/conditional model and obtains a new noise predictor $\bm{\epsilon}_{\theta}^{\prime}(\bm{x}_{t},t,c)\coloneqq s\bm{\epsilon}_{\theta}(\bm{x}_{t},t,c)+(1-s)\bm{\epsilon}_{\theta}(\bm{x}_{t},t,\emptyset)$, where $\emptyset$ is a special condition standing for the unconditional case, $s>0$ is the guidance scale that controls the trade-off between image-condition alignment and diversity; while classifier guidance [10] uses an extra classifier $p_{\phi}(c|\bm{x}_{t},t)$ to obtain a new noise predictor $\bm{\epsilon}_{\theta}^{\prime}(\bm{x}_{t},t,c)\coloneqq\bm{\epsilon}_{\theta}(\bm{x}_{t},t)-s\sigma_{t}\nabla_{\bm{x}}\log p_{\phi}(c|\bm{x}_{t},t)$.

In addition, except for the noise prediction, DPMs can be parameterized as score predictor $\bm{s}_{\theta}(\bm{x}_{t},t)$ to predict $\nabla_{\bm{x}}\log q_{t}(\bm{x}_{t},t)$, or data predictor $\bm{x}_{\theta}(\bm{x}_{t},t)$ to predict $\bm{x}_{0}$. Under variance-preserving (VP) noise schedule which satisfies $\alpha_{t}^{2}+\sigma_{t}^{2}=1$ [51], “v” predictor $\bm{v}_{\theta}(\bm{x}_{t},t)$ is also proposed to predict $\alpha_{t}\bm{\epsilon}-\sigma_{t}\bm{x}_{0}$ [45]. These different parameterizations are theoretically equivalent, but have an impact on the empirical performance when used in training [23, 59].

Among the methods for solving the diffusion ODE (1), recent works [56, 31, 32, 58] find that ODE solvers based on exponential integrators [18] are more efficient and robust at a small number of function evaluations (<50). Specifically, an insightful observation by Lu et al. [31] is that, by change-of-variable from $t$ to $\lambda_{t}\coloneqq\log(\alpha_{t}/\sigma_{t})$ (half of the log-SNR), the diffusion ODE is transformed to

$$\frac{\mathrm{d}\bm{x}_{\lambda}}{\mathrm{d}\lambda}=\frac{\dot{\alpha}_{\lambda}}{\alpha_{\lambda}}\bm{x}_{\lambda}-\sigma_{\lambda}\bm{\epsilon}_{\theta}(\bm{x}_{\lambda},\lambda),$$

(2)

where $\dot{\alpha}_{\lambda}\coloneqq\frac{\mathrm{d}\alpha_{\lambda}}{\mathrm{d}\lambda}$. By utilizing the semi-linear structure of the diffusion ODE and exactly computing the linear term [56, 31], we can obtain the ODE solution as Eq. (3) (left). Such exact computation of the linear part reduces the discretization errors [31]. Moreover, by leveraging the equivalence of different parameterizations, DPM-Solver++ [32] rewrites Eq. (2) by the data predictor $\bm{x}_{\theta}(\bm{x}_{\lambda},\lambda)\coloneqq(\bm{x}_{\lambda}-\sigma_{\lambda}\bm{\epsilon}_{\theta}(\bm{x}_{\lambda},\lambda))/\alpha_{\lambda}$ and obtains another ODE solution as Eq. (3) (right). Such solution does not need to change the pretrained noise prediction model $\bm{\epsilon}_{\theta}$ during the sampling process, and empirically outperforms previous samplers based on $\bm{\epsilon}_{\theta}$ [31].

$$\frac{\bm{x}_{t}}{\alpha_{t}}=\frac{\bm{x}_{s}}{\alpha_{s}}-\int_{\lambda_{s}}^{\lambda_{t}}e^{-\lambda}\bm{\epsilon}_{\theta}(\bm{x}_{\lambda},\lambda)\mathrm{d}\lambda,\quad\frac{\bm{x}_{t}}{\sigma_{t}}=\frac{\bm{x}_{s}}{\sigma_{s}}+\int_{\lambda_{s}}^{\lambda_{t}}e^{\lambda}\bm{x}_{\theta}(\bm{x}_{\lambda},\lambda)\mathrm{d}\lambda$$

(3)

However, to the best of our knowledge, the parameterizations for sampling are still manually selected and limited to noise/data prediction, which are not well-studied.

As mentioned in Sec. 2.2, it is promising to explore the semi-linear structure of diffusion ODEs for fast sampling [56, 31, 32]. Firstly, we reveal one key insight that we can choose the linear part according to Rosenbrock-type exponential integrators [19, 18]. To this end, we consider a general form of diffusion ODEs by rewriting Eq. (2) as

$$\frac{\mathrm{d}\bm{x}_{\lambda}}{\mathrm{d}\lambda}=\underbrace{\left(\frac{\dot{\alpha}_{\lambda}}{\alpha_{\lambda}}-\bm{l}_{\lambda}\right)\bm{x}_{\lambda}}_{\text{linear part}}-\underbrace{\vphantom{\int_{\lambda_{s}}^{\lambda}e^{-\int_{\lambda_{s}}^{r}\bm{s}_{\tau}\mathrm{d}\tau}\bm{b}_{r}\mathrm{d}r}(\sigma_{\lambda}\bm{\epsilon}_{\theta}(\bm{x}_{\lambda},\lambda)-\bm{l}_{\lambda}\bm{x}_{\lambda})}_{\text{non-linear part},\coloneqq\bm{N}_{\theta}(\bm{x}_{\lambda},\lambda)},$$

(4)

where $\bm{l}_{\lambda}$ is a $D$-dimensional undetermined coefficient depending on $\lambda$. We choose $\bm{l}_{\lambda}$ to restrict the Frobenius norm of the gradient of the non-linear part w.r.t. $\bm{x}$:

$$\bm{l}^{*}_{\lambda}=\operatornamewithlimits{argmin}_{\bm{l}_{\lambda}}\mathbb{E}_{p^{\theta}_{\lambda}(\bm{x}_{\lambda})}\|\nabla_{\bm{x}}\bm{N}_{\theta}(\bm{x}_{\lambda},\lambda)\|_{F}^{2},$$

(5)

where $p^{\theta}_{\lambda}$ is the distribution of samples on the ground-truth ODE trajectories at $\lambda$ (i.e., model distribution). Intuitively, it makes $\bm{N}_{\theta}$ insensitive to the errors of $\bm{x}$ and cancels all the “linearty” of $\bm{N}_{\theta}$. With $\bm{l}_{\lambda}=\bm{l}_{\lambda}^{*}$, by the “variation-of-constants” formula [1], starting from $\bm{x}_{\lambda_{s}}$ at time $s$, the exact solution of Eq. (4) at time $t$ is

$$\bm{x}_{\lambda_{t}}=\alpha_{\lambda_{t}}e^{-\int_{\lambda_{s}}^{\lambda_{t}}\bm{l}_{\lambda}\mathrm{d}\lambda}\bigg{(}\frac{\bm{x}_{\lambda_{s}}}{\alpha_{\lambda_{s}}}-\int_{\lambda_{s}}^{\lambda_{t}}e^{\int_{\lambda_{s}}^{\lambda}\bm{l}_{\tau}\mathrm{d}\tau}\underbrace{\bm{f}_{\theta}(\bm{x}_{\lambda},\lambda)}_{\text{approximated}}\mathrm{d}\lambda\bigg{)},$$

(6)

where $\bm{f}_{\theta}(\bm{x}_{\lambda},\lambda)\coloneqq\frac{\bm{N}_{\theta}(\bm{x}_{\lambda},\lambda)}{\alpha_{\lambda}}$. To calculate the solution in Eq. (6), we need to approximate $\bm{f}_{\theta}$ for each $\lambda\in[\lambda_{s},\lambda_{t}]$ by certain polynomials [31, 32].

Secondly, we reveal another key insight that we can choose different functions to be approximated instead of $\bm{f}_{\theta}$ and further reduce the discretization error, which is related to the total derivatives of the approximated function. To this end, we consider a scaled version of $\bm{f}_{\theta}$ i.e., $\bm{h}_{\theta}(\bm{x}_{\lambda},\lambda)\coloneqq e^{-\int_{\lambda_{s}}^{\lambda}\bm{s}_{\tau}\mathrm{d}\tau}\bm{f}_{\theta}(\bm{x}_{\lambda},\lambda)$ where $\bm{s}_{\lambda}$ is a $D$-dimensional coefficient dependent on $\lambda$, and then Eq. (6) becomes

$$\bm{x}_{\lambda_{t}}=\alpha_{\lambda_{t}}e^{-\int_{\lambda_{s}}^{\lambda_{t}}\bm{l}_{\lambda}\mathrm{d}\lambda}\bigg{(}\frac{\bm{x}_{\lambda_{s}}}{\alpha_{\lambda_{s}}}-\int_{\lambda_{s}}^{\lambda_{t}}e^{\int_{\lambda_{s}}^{\lambda}(\bm{l}_{\tau}+\bm{s}_{\tau})\mathrm{d}\tau}\underbrace{\bm{h}_{\theta}(\bm{x}_{\lambda},\lambda)}_{\text{approximated}}\mathrm{d}\lambda\bigg{)}.$$

(7)

Comparing with Eq. (6), we change the approximated function from $\bm{f}_{\theta}$ to $\bm{h}_{\theta}$ by using an additional scaling term related to $\bm{s}_{\lambda}$. As we shall see, the first-order discretization error is positively related to the norm of the first-order derivative $\bm{h}_{\theta}^{(1)}=e^{-\int_{\lambda_{s}}^{\lambda}\bm{s}_{\tau}\mathrm{d}\tau}(\bm{f}^{(1)}_{\theta}-\bm{s}_{\lambda}\bm{f}_{\theta})$. Thus, we aim to minimize $\|\bm{f}^{(1)}_{\theta}-\bm{s}_{\lambda}\bm{f}_{\theta}\|_{2}$, in order to reduce $\|\bm{h}_{\theta}^{(1)}\|_{2}$ and the discretization error. As $\bm{f}_{\theta}$ is a fixed function depending on the pretrained model, this minimization problem essentially finds a linear function of $\bm{f}_{\theta}$ to approximate $\bm{f}^{(1)}_{\theta}$. To achieve better linear approximation, we further introduce a bias term $\bm{b}_{\lambda}\in\mathbb{R}^{D}$ and construct a function $\bm{g}_{\theta}$ satisfying $\bm{g}_{\theta}^{(1)}=e^{-\int_{\lambda_{s}}^{\lambda}\bm{s}_{\tau}\mathrm{d}\tau}(\bm{f}^{(1)}_{\theta}-\bm{s}_{\lambda}\bm{f}_{\theta}-\bm{b}_{\lambda})$, which gives

$$\bm{g}_{\theta}(\bm{x}_{\lambda},\lambda)\coloneqq e^{-\int_{\lambda_{s}}^{\lambda}\bm{s}_{\tau}\mathrm{d}\tau}\bm{f}_{\theta}(\bm{x}_{\lambda},\lambda)-\int_{\lambda_{s}}^{\lambda}e^{-\int_{\lambda_{s}}^{r}\bm{s}_{\tau}\mathrm{d}\tau}\bm{b}_{r}\mathrm{d}r.$$

(8)

With $\bm{g}_{\theta}$, Eq. (7) becomes

$$\bm{x}_{\lambda_{t}}=\alpha_{\lambda_{t}}\underbrace{\vphantom{\int_{\lambda_{s}}^{\lambda}e^{-\int_{\lambda_{s}}^{r}\bm{s}_{\tau}\mathrm{d}\tau}\bm{b}_{r}\mathrm{d}r}e^{-\int_{\lambda_{s}}^{\lambda_{t}}\bm{l}_{\lambda}\mathrm{d}\lambda}}_{\text{linear coefficient}}\bigg{(}\frac{\bm{x}_{\lambda_{s}}}{\alpha_{\lambda_{s}}}-\int_{\lambda_{s}}^{\lambda_{t}}\underbrace{\vphantom{\int_{\lambda_{s}}^{\lambda}e^{-\int_{\lambda_{s}}^{r}\bm{s}_{\tau}\mathrm{d}\tau}\bm{b}_{r}\mathrm{d}r}e^{\int_{\lambda_{s}}^{\lambda}(\bm{l}_{\tau}+\bm{s}_{\tau})\mathrm{d}\tau}}_{\text{scaling coefficient}}\Big{(}\underbrace{\vphantom{\int_{\lambda_{s}}^{\lambda}e^{-\int_{\lambda_{s}}^{r}\bm{s}_{\tau}\mathrm{d}\tau}\bm{b}_{r}\mathrm{d}r}\bm{g}_{\theta}(\bm{x}_{\lambda},\lambda)}_{\text{approximated}}+\underbrace{\int_{\lambda_{s}}^{\lambda}e^{-\int_{\lambda_{s}}^{r}\bm{s}_{\tau}\mathrm{d}\tau}\bm{b}_{r}\mathrm{d}r}_{\text{bias coefficient}}\Big{)}\mathrm{d}\lambda\bigg{)}.$$

(9)

Such formulation is equivalent to Eq. (3) but introduces three types of coefficients and a new parameterization $\bm{g}_{\theta}$. We show in Appendix I.1 that the generalized parameterization $\bm{g}_{\theta}$ in Eq. (8) can cover a wide range of parameterization families in the form of $\bm{\psi}_{\theta}(\bm{x}_{\lambda},\lambda)=\bm{\alpha}(\lambda)\bm{\epsilon}_{\theta}(\bm{x}_{\lambda},\lambda)+\bm{\beta}(\lambda)\bm{x}_{\lambda}+\bm{\gamma}(\lambda)$. We aim to reduce the discretization error by finding better coefficients than previous works [31, 32].

Now we derive the concrete formula for analyzing the first-order discretization error. By replacing $\bm{g}_{\theta}(\bm{x}_{\lambda},\lambda)$ with $\bm{g}_{\theta}(\bm{x}_{\lambda_{s}},\lambda_{s})$ in Eq. (9), we obtain the first-order approximation $\hat{\bm{x}}_{\lambda_{t}}=\alpha_{\lambda_{t}}e^{-\int_{\lambda_{s}}^{\lambda_{t}}\bm{l}_{\lambda}\mathrm{d}\lambda}\left(\frac{\bm{x}_{\lambda_{s}}}{\alpha_{\lambda_{s}}}-\int_{\lambda_{s}}^{\lambda_{t}}e^{\int_{\lambda_{s}}^{\lambda}(\bm{l}_{\tau}+\bm{s}_{\tau})\mathrm{d}\tau}\left(\bm{g}_{\theta}(\bm{x}_{\lambda_{s}},\lambda_{s})+\int_{\lambda_{s}}^{\lambda}e^{-\int_{\lambda_{s}}^{r}\bm{s}_{\tau}\mathrm{d}\tau}\bm{b}_{r}\mathrm{d}r\right)\mathrm{d}\lambda\right)$. As $\bm{g}_{\theta}(\bm{x}_{\lambda_{s}},\lambda_{s})=\bm{g}_{\theta}(\bm{x}_{\lambda},\lambda)+(\lambda_{s}-\lambda)\bm{g}_{\theta}^{(1)}(\bm{x}_{\lambda},\lambda)+\mathcal{O}((\lambda-\lambda_{s})^{2})$ by Taylor expansion, it follows that the first-order discretization error can be expressed as

$$\hat{\bm{x}}_{\lambda_{t}}-\bm{x}_{\lambda_{t}}=\alpha_{\lambda_{t}}e^{-\int_{\lambda_{s}}^{\lambda_{t}}\bm{l}_{\lambda}\mathrm{d}\lambda}\bigg{(}\int_{\lambda_{s}}^{\lambda_{t}}e^{\int_{\lambda_{s}}^{\lambda}\bm{l}_{\tau}\mathrm{d}\tau}(\lambda-\lambda_{s})\Big{(}\bm{f}^{(1)}_{\theta}(\bm{x}_{\lambda},\lambda)-\bm{s}_{\lambda}\bm{f}_{\theta}(\bm{x}_{\lambda},\lambda)-\bm{b}_{\lambda}\Big{)}\mathrm{d}\lambda\bigg{)}+\mathcal{O}(h^{3}),$$

(10)

where $h=\lambda_{t}-\lambda_{s}$. Thus, given the optimal $\bm{l}_{\lambda}=\bm{l}_{\lambda}^{*}$ in Eq. (5), the discretization error $\hat{\bm{x}}_{\lambda_{t}}-\bm{x}_{\lambda_{t}}$ mainly depends on $\bm{f}^{(1)}_{\theta}-\bm{s}_{\lambda}\bm{f}_{\theta}-\bm{b}_{\lambda}$. Based on this insight, we choose the coefficients $\bm{s}_{\lambda},\bm{b}_{\lambda}$ by solving

$$\bm{s}_{\lambda}^{*},\bm{b}_{\lambda}^{*}=\operatornamewithlimits{argmin}_{\bm{s}_{\lambda},\bm{b}_{\lambda}}\mathbb{E}_{p^{\theta}_{\lambda}(\bm{x}_{\lambda})}\left[\left\|\bm{f}^{(1)}_{\theta}(\bm{x}_{\lambda},\lambda)-\bm{s}_{\lambda}\bm{f}_{\theta}(\bm{x}_{\lambda},\lambda)-\bm{b}_{\lambda}\right\|^{2}_{2}\right].$$

(11)

For any $\lambda$, $\bm{l}_{\lambda}^{*},\bm{s}_{\lambda}^{*},\bm{b}_{\lambda}^{*}$ all have analytic solutions involving the Jacobian-vector-product of the pretrained model $\bm{\epsilon}_{\theta}$, and they can be unbiasedly evaluated on a few datapoints $\{\bm{x}_{\lambda}^{(n)}\}_{K}\sim p_{\lambda}^{\theta}(\bm{x}_{\lambda})$ via Monte-Carlo estimation (detailed in Section 3.4 and Appendix C.1.1). Therefore, we call $\bm{l}_{\lambda},\bm{s}_{\lambda},\bm{b}_{\lambda}$ empirical model statistics (EMS). In the following sections, we’ll show that by approximating $\bm{g}_{\theta}$ with Taylor expansion, we can develop our high-order solver for diffusion ODEs.

In this section, we propose our high-order solver for diffusion ODEs with local accuracy and global convergence order guarantee by leveraging our proposed solution formulation in Eq. (9). The proposed solver and analyses are highly motivated by the methods of exponential integrators [17, 18] in the ODE literature and their previous applications in the field of diffusion models [56, 31, 32, 58]. Though the EMS are designed to minimize the first-order error, they can also help our high-order solver (see Appendix I.2).

For simplicity, denote $\bm{A}(\lambda_{s},\lambda_{t})\coloneqq e^{-\int_{\lambda_{s}}^{\lambda_{t}}\bm{l}_{\tau}\mathrm{d}\tau}$, $\bm{E}_{\lambda_{s}}(\lambda)\coloneqq e^{\int_{\lambda_{s}}^{\lambda}(\bm{l}_{\tau}+\bm{s}_{\tau})\mathrm{d}\tau}$, $\bm{B}_{\lambda_{s}}(\lambda)\coloneqq\int_{\lambda_{s}}^{\lambda}e^{-\int_{\lambda_{s}}^{r}\bm{s}_{\tau}\mathrm{d}\tau}\bm{b}_{r}\mathrm{d}r$. Though high-order ODE solvers essentially share the same mathematical principles, since we utilize a more complicated parameterization $\bm{g}_{\theta}$ and ODE formulation in Eq. (9) than previous works [56, 31, 32, 58], we divide the development of high-order solver into simplified local and global parts, which are not only easier to understand, but also neat and general for any order.

In this section, we introduce some practical techniques that further improve the sample quality in the case of small NFE or large guidance scales.

Pseudo-order solver for small NFE. When NFE is extremely small (e.g., 5$\sim$10), the error at each timestep becomes rather large, and incorporating too many previous values by high-order solver at each step will cause instabilities. To alleviate this problem, we propose a technique called pseudo-order solver: when estimating the $k$-th order derivative, we only utilize the previous $k+1$ function values of $\bm{g}_{\theta}$, instead of all the $n$ previous values as in Eq. (13). For each $k$, we can obtain $\hat{\bm{g}}^{(k)}_{s}$ by solving a part of Eq. (13) and taking the last element:

$$\begin{pmatrix}\delta_{1}&\delta_{1}^{2}&\cdots&\delta_{1}^{k}\\ \vdots&\vdots&\ddots&\vdots\\ \delta_{k}&\delta_{k}^{2}&\cdots&\delta_{k}^{k}\end{pmatrix}\begin{pmatrix}\cdot\\ \vdots\\ \frac{\hat{\bm{g}}^{(k)}_{s}}{k!}\end{pmatrix}=\begin{pmatrix}\bm{g}_{i_{1}}-\bm{g}_{s}\\ \vdots\\ \bm{g}_{i_{k}}-\bm{g}_{s}\end{pmatrix},\quad k=1,2,\dots,n$$

(15)

In practice, we do not need to solve $n$ linear systems. Instead, the solutions for $\hat{\bm{g}}^{(k)}_{s},k=1,\dots,n$ have a simpler recurrence relation similar to Neville’s method [36] in Lagrange polynomial interpolation. Denote $i_{0}\coloneqq s$ so that $\delta_{0}=\lambda_{i_{0}}-\lambda_{s}=0$, we have

Proof in Appendix B.3. Note that the pseudo-order solver of order $n>2$ no longer has the guarantee of $n$-th order of accuracy, which is not so important when NFE is small. In our experiments, we mainly rely on two combinations: when we use $n$-th order predictor, we then combine it with $n$-th order corrector or $(n+1)$-th pseudo-order corrector.

Half-corrector for large guidance scales. When the guidance scale is large in guided sampling, we find that corrector may have negative effects on the sample quality. We propose a useful technique called half-corrector by using the corrector only in the time region $t\leq 0.5$. Correspondingly, the strategy that we use corrector at each step is called full-corrector.

In this section, we give some implementation details about how to compute and integrate the EMS in our solver and how to adapt them to guided sampling.

Estimating EMS. For a specific $\lambda$, the EMS $\bm{l}_{\lambda}^{*},\bm{s}_{\lambda}^{*},\bm{b}_{\lambda}^{*}$ can be estimated by firstly drawing $K$ (1024$\sim$4096) datapoints $\bm{x}_{\lambda}\sim p^{\theta}_{\lambda}(\bm{x}_{\lambda})$ with 200-step DPM-Solver++ [32] and then analytically computing some terms related to $\bm{\epsilon}_{\theta}$ (detailed in Appendix C.1.1). In practice, we find it both convenient and effective to choose the distribution of the dataset $q_{0}$ to approximate $p_{0}^{\theta}$. Thus, without further specifications, we directly use samples from $q_{0}$.

Estimating integrals of EMS. We estimate EMS on $N$ ($120\sim 1200$) timesteps $\lambda_{j_{0}},\lambda_{j_{1}},\dots,\lambda_{j_{N}}$ and use trapezoidal rule to estimate the integrals in Eq. (12) (see Appendix I.3 for the estimation error analysis). We also apply some pre-computation for the integrals to avoid extra computation costs during sampling, detailed in Appendix C.1.2.

Adaptation to guided sampling. Empirically, we find that within a common range of guidance scales, we can simply compute the EMS on the model without guidance, and it can work for both unconditional sampling and guided sampling cases. See Appendix J for more discussions.

By comparing with existing diffusion ODE solvers that are based on exponential integrators [56, 31, 32, 58], we can conclude that (1) Previous ODE formulations with noise/data prediction are special cases of ours by setting $\bm{l}_{\lambda},\bm{s}_{\lambda},\bm{b}_{\lambda}$ to specific values. (2) Our first-order discretization can be seen as improved DDIM. See more details in Appendix A.

We present the results in $5\sim 20$ number of function evaluations (NFE), covering both few-step cases and the almost converged cases, as shown in Fig. 3 and Fig. 4. For the sake of clarity, we mainly compare DPM-Solver-v3 to DPM-Solver++ [32] and UniPC [58], which are the most state-of-the-art diffusion ODE solvers. We also include the results for DEIS [56] and Heun’s 2nd order method in EDM [21], but only for the datasets on which they originally reported. We don’t show the results for other methods such as DDIM [48], PNDM [28], since they have already been compared in previous works and have inferior performance. The quantitative results on CIFAR10 [24] are listed in Table 1, and more detailed quantitative results are presented in Appendix F.

Unconditional sampling We first evaluate the unconditional sampling performance of different methods on CIFAR10 [24] and LSUN-Bedroom [55]. For CIFAR10 we use two pixel-space DPMs, one is based on ScoreSDE [51] which is a widely adopted model by previous samplers, and another is based on EDM [21] which achieves the best sample quality. For LSUN-Bedroom, we use the latent-space Latent-Diffusion model [43]. We apply the multistep 3rd-order version for DPM-Solver++, UniPC and DPM-Solver-v3 by default, which performs best in the unconditional setting. For UniPC, we report the better result of their two choices $B_{1}(h)=h$ and $B_{2}(h)=e^{h}-1$ at each NFE. For our DPM-Solver-v3, we tune the strategies of whether to use the pseudo-order predictor/corrector at each NFE on CIFAR10, and use the pseudo-order corrector on LSUN-Bedroom. As shown in Fig. 3, we find that DPM-Solver-v3 can achieve consistently better FID, which is especially notable when NFE is 5$\sim$10. For example, we improve the FID on CIFAR10 with 5 NFE from $23$ to $12$ with ScoreSDE, and achieve an FID of $2.51$ with only $10$ NFE with the advanced DPM provided by EDM. On LSUN-Bedroom, with around 12 minutes computing of the EMS, DPM-Solver-v3 converges to the FID of 3.06 with 12 NFE, which is approximately 60% sampling cost of the previous best training-free method (20 NFE by UniPC).

Conditional sampling. We then evaluate the conditional sampling performance, which is more widely used since it allows for controllable generation with user-customized conditions. We choose two conditional settings, one is classifier guidance on pixel-space Guided-Diffusion [10] model trained on ImageNet-256 dataset [9] with 1000 class labels as conditions; the other is classifier-free guidance on latent-space Stable-Diffusion model [43] trained on LAION-5B dataset [46] with CLIP [41] embedded text prompts as conditions. We evaluate the former at the guidance scale of $2.0$, following the best choice in [10]; and the latter at the guidance scale of $1.5$ (following the original paper) or $7.5$ (following the official code) with prompts random selected from MS-COCO2014 validation set [26]. Note that the FID of Stable-Diffusion samples saturates to 15.0$\sim$16.0 even within 10 steps when the latent codes are far from convergence, possibly due to the powerful image decoder (see Appendix H). Thus, following [32], we measure the mean square error (MSE) between the generated latent code $\hat{\bm{x}}$ and the ground-truth solution $\bm{x}^{*}$ (i.e., $\|\hat{\bm{x}}-\bm{x}^{*}\|_{2}^{2}/D$) to evaluate convergence, starting from the same Gaussian noise. We obtain $\bm{x}^{*}$ by 200-step DPM-Solver++, which is enough to ensure the convergence.

We apply the multistep 2nd-order version for DPM-Solver++, UniPC and DPM-Solver-v3, which performs best in conditional setting. For UniPC, we only apply the choice $B_{2}(h)=e^{h}-1$, which performs better than $B_{1}(h)$. For our DPM-Solver-v3, we use the pseudo-order corrector by default, and report the best results between using half-corrector/full-corrector on Stable-Diffusion ($s=7.5$). As shown in Fig. 4, DPM-Solver-v3 can achieve better sample quality or convergence at most NFEs, which indicates the effectiveness of our method and techniques under the conditional setting. It’s worth noting that UniPC, which adopts an extra corrector, performs even worse than DPM-Solver++ when NFE<10 on Stable-Diffusion ($s=7.5$). With the combined effect of the EMS and the half-corrector technique, we successfully outperform DPM-Solver++ in such a case. Detailed comparisons can be found in the ablations in Appendix G.

We further visualize the estimated EMS in Fig. 5. Since $\bm{l}_{\lambda},\bm{s}_{\lambda},\bm{b}_{\lambda}$ are $D$-dimensional vectors, we average them over all dimensions to obtain a scalar. From Fig. 5, we find that $\bm{l}_{\lambda}$ gradually changes from $1$ to near $0$ as the sampling proceeds, which suggests we should gradually slide from data prediction to noise prediction. As for $\bm{s}_{\lambda},\bm{b}_{\lambda}$, they are more model-specific and display many fluctuations, especially for ScoreSDE model [51] on CIFAR10. Apart from the estimation error of the EMS, we suspect that it comes from the fluctuations of $\bm{\epsilon}_{\theta}^{(1)}$, which is caused by the periodicity of trigonometric functions in the positional embeddings of the network input. It’s worth noting that the fluctuation of $\bm{s}_{\lambda},\bm{b}_{\lambda}$ will not cause instability in our sampler (see Appendix J).

We present some qualitative comparisons in Fig. 6 and Fig. 1. We can find that previous methods tend to have a small degree of color contrast at small NFE, while our method is less biased and produces more visual details. In Fig. 1(b), we can observe that previous methods with corrector may cause distortion at large guidance scales (in the left-top image, a part of the castle becomes a hill; in the left-bottom image, the hill is translucent and the castle is hanging in the air), while ours won’t. Additional samples are provided in Appendix K.