DPM-Solver-v3: Improved Diffusion ODE Solver with Empirical Model Statistics
Authors: Kaiwen Zheng
Cheng Lu
Jianfei Chen
Jun Zhu
Venue: NeurIPS 2023 poster (2023)
Keywords: diffusion models
fast sampling
ODE solver
TL;DR: We propose a fast ODE solver for diffusion models, achieving SOTA acceleration performance especially in 5 to 10 NFEs.

Abstract: Diffusion probabilistic models (DPMs) have exhibited excellent performance for high-fidelity image generation while suffering from inefficient sampling. Recent works accelerate the sampling procedure by proposing fast ODE solvers that leverage the specific ODE form of DPMs. However, they highly rely on specific parameterization during inference (such as noise/data prediction), which might not be the optimal choice. In this work, we propose a novel formulation towards the optimal parameterization during sampling that minimizes the first-order discretization error of the ODE solution. Based on such formulation, we propose \textit{DPM-Solver-v3}, a new fast ODE solver for DPMs by introducing several coefficients efficiently computed on the pretrained model, which we call \textit{empirical model statistics}. We further incorporate multistep methods and a predictor-corrector framework, and propose some techniques for improving sample quality at small numbers of function evaluations (NFE) or large guidance scales. Experiments show that DPM-Solver-v3 achieves consistently better or comparable performance in both unconditional and conditional sampling with both pixel-space and latent-space DPMs, especially in 5$\sim$10 NFEs. We achieve FIDs of 12.21 (5 NFE), 2.51 (10 NFE) on unconditional CIFAR10, and MSE of 0.55 (5 NFE, 7.5 guidance scale) on Stable Diffusion, bringing a speed-up of 15\%$\sim$30\% compared to previous state-of-the-art training-free methods. Code is available at \url{https://github.com/thu-ml/DPM-Solver-v3}.

DPM-Solver-v3: Improved Diffusion ODE Solver with Empirical Model Statistics
Kaiwen Zheng*t1, Cheng Lu"!, Jianfei Chen', Jun Zhu 123
Dept. of Comp. Sci. & Tech., Institute for AI, BNRist Center, THBI Lab Tinghuua Joint ML Center, Tsinghua University, Beijing, China "Shengshu Technology, Beijing Pazzou Lab (Huangpu), Guangzhou, China {zkwthu, lucheng. 1c15)@gmail.cmmmm com; {jianfeic, dcszj)@tsinghua. edu,
arXiv:submit/520042 [cs.CV] 27 Oct 2023
Abstract
Diffusion probabilistic models (DPMs) have exhibited excellent performance for high-fidelity image generation while suffering from inefficient sampling. Re- cent works accelerate the sampling procedure by proposing fast ODE solvers that leverage the specific ODE form of DPMs. However, they highly rely on specific parameterization during inference (such as noise/data prediction), which might not be the optimal choice. In this work, we propose novel formulation towards the optimal parameterization during sampling that minimizes the first- order discretization error the ODE solution. Based on such formulation, we propose DPM-Solver-yy coefficients efficiently computed on the pretrained model, which we call empir- ical model statistics. We further incorporate multistep methods and predictor- corrector framework, and propose some techniques for improving sample quality small numbers of function evaluations (NFR) large guidance scales. Ex- periments show that DPM-Solver-v3 achieves consistently better or comparable performance in both unconditional and conditional sampling with both pixel- space and latent-space DPMs, especially in 5~ 10 NFEs. We achieve FIDs of 12.21 (5 NFE), 2.51 (10 NFE) unconditional CIFARIÓ, and MSE of 0.55 NFE, 7.5 guidance scale) on Stable Diffusion, bringing a speed-up of 15%~~30% compared to previous taate-of-the-asd training- free methods. Code is available https: /etp//ggttttbbt com/thu ttt//thummmmm mwtpmmmm-PPPPM---- Solver-
fast ODE solver for DPMs by introducing several
Introduction
Diffusion probabilistic models (DPMs) [47, 15, are a class of state-of-the-art image generators. By training with a strong encoder, large model, and massive data well as leveraging techniques such as guided sampling, DPMs are capable of generating high-resolution photorealistic and artistic images text-to-image tasks. However, generate high-quality visual content, DPMs usually require hundreds of model evaluations to gradually remove noise using pretrained model [15]. which is much more time-consuming compared to other deep generative models such as generative adversarial networks (GANs) [13]. The sampling overhead of DPMs emerges as crucial obstacle hindering their integration into downstream tasks.
To accelerate the sampling process of DPMs, one can employ training- based methods [37, 53, 45] or training- free samplers [48, 51, 28, 3, 52, 56]. We focus the latter approach since it requires
Work done during an internship at Shengshu Technology "Equal contribution
*Corresponding author
37th Conference on Neural Information Processing Systems (NcurIPS 2023).

---

DMM-Solvrr++ [32] (MSE 0.60)
(b) UniPC [58] (MSF (MSF..66))
(c) DPM-Solver-v3 (Ours, 0.55)
Random samples of Stable- Diffusion [43] with a classifier-free guidance scale 7.5, using only number of function evaluations (NFE) and text prompt "A beautifl castle beside a waterfall in the woods, by Josef Thoma, matte painting, trending on artstation HQ".
Figure 1:
extra training and is more flexible. Recent advanced training- mhining-ffree samplers [56, 31, 32, 58] mainly rely on the ODE form of DPMs, since its absence of stochasticity is essential for high-quality samples in around 20 steps. Besides, ODE solvers built on exponential integrators [18] converge faster. To change the original diffusion ODE into the form of exponential integrators, they need to cancel its linear term and obtain an ODE solution, where only the noise predictor needs to be approximated, and other terms can be exactly computed. Besides, Lu et al. [32] find that is better use another ODE solution where instead the data predictor needs to be approximated. How to choose such model parameterization (e.g., noise/data prediction) in the sampling of DPMs remains unrevealed.
In this work, we systematically study the problem of model parameterization and ODE formulation for fast sampling of DPMs. We first show that by introducing three types of coefficients, the original ODE solution can be reformulated to an equivalent one that contains a new model parameterization. Besides, inspired by exponential Rosenbrock-
methods [19] and first-order discretization error analysis, we also show that there exists an optimal set of coefficients efficiently computed on the pretrained model, which we call empirical model statistics (EMS). Building upon our novel ODE formulation, we further propose a new high-order solver for diffusion ODEs named DPM-Solver- v3, which includes a multistep predictor-correcto framework any order, as well as some novel techniques such as pseudo high-order method to boost the performance at extremely few steps or large guidance scale.
We conduct extensive experiments with both pixel-space and latent-space DPMs to verify the effec- tiveness of DPM-Solver-v3. Compared previous fast samplers, DPM-Solver-v3 can consistently improve the sample quality in 5~ steps, and make significant advancement within 10 steps.
2 Background
2.1 Diffusion Probabilistic Models
Suppose we have a D-dimensional data distribution go(z)). Diffusion probabilistic models (DPMs) [47, 15, 51] define a forward diffusion process {qr)< to gradually degenerate the data 2n ~ 9o (2)) with Gaussian noise, which satisfies the transition kernel gor (z,(0n) = N(ayz I).). such that (r(2)) is approximately pure Gaussian. Og2 oz smooth scalar functions of t,, which are called noise schedule. The transition can be easily applied by 2; ag + a,c, e ~ (N(0, I). To train DPMs, neural network (2 is usually parameterized to predict the noise by min- imizing Enzorgg (aoorgg(mm)))eeee (onummm)e-Wii ee (=. (), [w[t)]lco t) <l13) where w(t) is weighting function. Sampling of DPMs can be performed by solving diffusion ODE [51] from time T to time 0:
dx;
f(t) + (t)
co(
(I)
dt
20,
where S(1) o. (2(1) de
dof 2 a logga of [23]. addition, the conditional sampling by
di
DPMs can be conducted by guided sampling [10, 16] with a conditional noise predictor ep (z,, l, c).
where e the condition. Specifically, classifier guidance [16] combines the unconditional/conn- ditional model and obtains a new noise predictor e (xt, t, c) = sea( t, c) + (1 s)Ea (x, 0).

---

where / is special condition standing for the unconditional case, > 0 is the guidance scale that controls the trade-off between image-condition alignment and diversity; while classifier guidance [10] uses an extra classifier (olcaa 1) to obtain a new noise predictor (,,,,, (, = c0(2;; 1) - sa,V log Po (c)
In addition, except for the noise prediction, DPMs can be parameterized as score predictor 8, (2te, t) predict V log (x,,,,) t).
or data predictor 20(#;, 1) to predict 20. Under prraeccepreserving (VP) noise schedule which satisfies of + of = 1 [51]. "v" predictor Va (2r t) is also proposed to predict aye o [45]. These different parameterizations are theoretically equivalent, but have an impact on the empirical performance when used in training [23, 59].
2.2 Fast Sampling of DPMs with Exponential Integrators
Among the methods for solving the diffusion ODE (1), recent works [56, 31, 58] find that ODE solvers based on exponential integrators [18] are more efficient and robust at a small number of function evaluations (<50). Specifically, insightful observation by Lu al. [31] is that, by change-
of-variable from & to X = log(ar/o, ox) (half of the log-SNR), the diffusion ODE is transformed
dæ
â> -ax xxE
dx
ax
where da := dos.
dx 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 [3]]. Moreover, by leveraging the equivalence of
different parameterizations, DPM-Solver++ [32] rewrites Eq. (2) by the data predictor 2& (zx A) = (azx Eg x aax and obtains another ODE solution Eq. (3) (right). Such solution does not need to change the pretrained noise prediction model Eg during the sampling process, and empirically outperforms previous samplers based co [31].
Et at
-==- Ox
Ea A)dA,
-
Ee o.
S/a ex AJdA
at
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.
Method
We now present our method. We start with new formulation of the ODE solution with extra coefficients, followed by our high-order solver and some practical considerations. In the following discussions, we assume all the products between vectors are element-wise, and (A) (xx, A)
d f(mx.A) is the k-th order total derivative of any function f w.r.t. X.
dxd
3.1 Improved Formulation of Exact Solutions of Diffusion ODEs
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-tyer exponential integrators [19, 18]. To this end, we consider a general form of diffusion ODEs by rewriting Eq. (2)
dæ dx
á' xx (oxEo A) 1xxx);
αx
lincar part
non linear part.:= (o(
where 1x is D-dimensional undetermined coefficient depending on A. We choose Ly to restrict the Frobenius norm of the gradient of the non-linear part w.r.t. æ::
1 argmin EE, Egg (mx) V NoNNO 1)
(5)
l
where o is the distribution samples on the ground-truth ODE trajectories X (i.e., model distribution). Intuitively, it makes N₂ insensitive to the errors of and cancels all the "lincarty" of

---

No. With LA = 13 by the "variation-of-constaat formula [1], starting from 2x, at time the exact solution of Eq. (4) at time t
xx = axxe 12g ladA
Ih.. - ax
-Sa eli, bedr folax.))as)
(6)
approximated
where fo( No(w..)) To calculate the solution in Eq. (6), we need to approximate fo for 0*
each A E [N,, Ao] by certain polynomials [31, 32].
Secondly, we reveal another key insight that we can choose different functions to be approximated instead of fo and further reduce the discretization error, which related to the total derivatives of the approximated function. To this end, we consider scaled version of fo i.e., ho (2x, =
e [2, a-dr fa(a;,,) A) where 88 is a D-dimensional coefficient dependent on X, and then Eq. (6)
becomes
xy = ox,e 12 Iadd
xx. e. (1, ts, )dr ho X)
a22
approximated
Comparing with Eq. (6), we change the approximated function from fo he by using an additional scaling term related to sx. As we shall see, the first-order discretization error is positively related to
the norm the first-order derivative hgg (1) = e 12, e, (fe") sxfe). Thus, we aim to minimize
|| fololl,, in order reduce | and the discretization error. As fo is fixed function
depending on the pretrained model, this minimization problem essentially finds linear function of further introduce bias term
fe to approximate f/". To achieve better linear approximation,
bx RR and construct function go satisfying 99 (1) e- J2 sn (d(((())) sxfo bx), which gives
ge( == e- -Sh dr fo(xx, /// e Sx, dr dr.
e-
(8)
With go. Eq. (7) becomes
xA,- = 022 e 121 LadA
xx
(am de
(a, go(2x,, A) 2)+ S.a e -- x, dr )ndr))ay (9)
an
linear cueflicient
scaling coeflicient approximated
coefficient
Such formulation is equivalent Eq. (3) but introduces three types of coefficients and new parameterization go- We show in Appendix I.1 that the generalized parameterization go in Eq. (8)
can cover wide range of parameterization families in the form of (bo( A) = axA)ez (æxx, (a)))((x)))) (())=x + (A). We aim to reduce the discretization error by finding better coefficients than previous works [2].
Now we derive the concrete formula for analyzing the first-order discretization error. By replacing g0 A) with g 2s) in Eq. (9), we obtain the first-order approximation J elxt (l++r))rr (gx( do) + S². - Jr a. b,dr) dx). As
Φx. ax,e 1 ladd
Td ax2
(1) (x), + O(A - xx)22 by Taylor expansion, it follows that the first-order discretization error can be expressed
x A.) = (x A) + (A, - A)ge
J Add (J,) ex l,dr (A- (A-2a) 51 (xx,
snfelax, A) - -)) bs) bs)da) =0(h3)
c xR. = axx.
(10)
where h = 21 - Ag. Thus, given the optimal Ix = 1% in Eq. (5), the discretization error ex. xx.
mainly depends s(1)) SrSO bx. Based on this insight, we choose the coefficients 8x, by by solving
x bj argmin g (as. s fo A) by
(11)
sx.bs
For any A, 1. s.. bū all have analytic solutions involving the Jacobian Jacobian-wector wwwp---eeeooooorrooee product of the pretrained
model ep, and they can unbiasedly evaluated on a few datapoints (a²(") x ~ p (mxx) via Monte- Carlo estimation (detailed in Section 3.4 and Appendix C.1.1). Therefore, we call I,, by empirical model statistics (EMS). In the following sections, we'll show that by approximating go with Taylor expansion, we can develop our high-order solver for diffusion ODEs.

---

Predictor
(1...))
Ap
Fi-
estimation
9.
:
9--2
(1..9h.))
Concetor
(a) Local Approximation
(b) Multistep Predictor-Comeet
Figure Illustration of our high-order solver. (a) (n + 1)-th order local approximation from time to time f, provided n extra function values of go- (b) Multistep predictor-corrcen procedure as our global solver. A combination of second-order predictor and sccond-order corrector is shown. a,,i-1, b,,i-1, C,.-1 abbreviations of coefficients' in Eq. (8).
3.2 Developing High-Order Solver
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] 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 A().,Ae) = - d. (A) = (L. Bx,
K e- a-dr bydr. Though high-order ODE solvers essentially share the same mathematical principles, since we utilize a more complicated parameterization go and ODE formulation 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 understand, but also neat and general for any order.
3.2.1 Local Approximation
Firstly, we derive formulas and properties of the local approximation, which describes how we transit locally from time to time t. It can be divided into two parts: discretization of the integral in Eq. (9) and estimating the high-order derivatives in the Taylor expansion.
Discretization. Denote ga (k) == 9 (xx,, A²). For 2 O. to obtain the (n++))-th order discretization of Eq. (9), we take the n-th order Taylor expansion of go (=x,, w.r.t À at A,: (0 X)
(k) En (A-2,)) gs + O((A A.) Substituting it into Eq. (9) yields
80- Ae / Ex, (A) as J
d E (2)
(A Aa)a
= A(),;2))
Σ
g
kl dd +0(h+)2
at
(12)
Here only need to estimate the k-thh order total derivatives go other terms are determined once given As, A and Ix, 8x, bx, which we'll discuss next.
(k) (xx A₂) for 0 n, since the
High-order derivative estimation. For (n+1) 1)-th order approximation,
use the finite difference of
(k) go (x A) at previous n++1 steps >e,,...., ., Aig >, to estimate each go (mx,; Ag). Such derivation is to match the coefficients of Taylor expansions. Concretely, denote dg := >is Aas Ii, == go Aix):
and the estimated high-order derivatives g! can be solved by the following linear system:
81
8
8
g
9h.
9s
=
(13)
8n
52
g:..
9
5

---

Then by substituting (k) into Eq. (12) and dropping the 2) error terms, we obtain the (x+1)--tt order local approximation:
i; a2
A()o.
( Ex, (A) (A)dA
Σ
R
( 1x) _ k!
d Ex, (x, /
di (14)
as
where 9o (I>, ==)== 9s. Eq. (13) and Eq. (14) provide an update rule to transit from time s to time t and get an approximated solution é;, when we already have the solution 2,. For (n + 1 1)-th order approximation, we need extra solutions 2 BAs, and their corresponding function values gix We illustrate the procedure in Fig. 2(a) and summarize in Appendix C.2. In the following theorem, we show that under some assumptions, such local approximation has a guarantee of order of accuracy.
Theorem 3.1 (Local order of accuracy, proof in Appendix B.2. Suppose DA are exact (i.e., the ground-truth ODE trajectory passing xx)) for k = I,..., then under some regularity conditions detailed in Appendix B.I, 1, the local truncation error &t - De = (li+t2), which means the local approximation has (n 1)-th order of accuracy.
Besides, we have the following theorem showing that, whatever the order is, the local approximation unbiased given our choice of x by in Eq. (1I). In practice, the phenomenon of reduced bias can be empirically observed (Section 4.3).
Theorem 3.2 (Local unbiasedness, proof in Appendix B.4), Given the optimal 8 b3, in Eq. (11), For the (n + 1) 1)-th order approximation, suppose Dx;;"" ..... FA,, are on the ground-trust ODE trajectory passing ax,, then Eyg (a.) [ig =.. 0.
3.2.2 Global Solver
Given M + 1 time steps [4 starting from some initial value, we can repeat the local approxi- mation M times to make consecutive transitions from each 4-11 to i, until we reach an acceptable solution. At each step, we apply multistep methods [1] by caching and reusing the previous n values timesteps fi-1-n..... f4-2, which proven to be more stable when NFE small [32, 56]. Moreover, we also apply the predictor- predictor-corrector method [58] to refine the approximation at each step without introducing extra NFE. Specifically, the (n 1)-th order predictor is the case of the local
approximation when we choose (L.,....,, lig 8. (4i-1-m:- 44-2,44 44 44), and the (n + 1)-th order corrector is the case when we choose (t.......,,44 8,,,,, 8,1) = (4,-,,--. .......,,. 14,-,, 4,, 12-1, 1.). We present our (n + 1)-th order multistep predictor-corrector procedure in Appendix C.2. We also illustrate a second-order case in Fig. 2(b). Note that different from previous works, in the local transition from 4 4, the previous function values 9i, (1 < k S n)) used for derivative estimation are dependent on and are different during the sampling process because go is dependent on the current 4-1 (see Eq. (8)). Thus, we directly cache i; é; and reuse them to compute 9; in the subsequent steps. Notably, our proposed solver also has a global convergence guarantee, as shown in the following theorem. For simplicity, we only consider the predictor case and the case with corrector can also be proved by similar derivations in [58].
Theorem 3.3 (Global order of convergence, proof in Appendix B.2.2). For (n+ 1)-th order predictor, if we iteratively compute a sequence {æx M to approximate the true solutions {ax).M at (t, ( M
then under both local and global assumptions detailed in Appendix B. the final error @MM EM Ohh=+1), where
denotes the element-wise absolute value, and h max MaxxssiMMl (A, A
3.3 Practical Techniques
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~ 10), the error 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 technique called pseudo- order solver: when estimating the k-th order derivative, we only utilize the previous + 1 function
values of ge. instead of all the previous values as in Eq. (13). For each k. we can obtain g by (k)

---

solving part of Eq. (13) and taking the last element:
87
9h 9s
1,2,.....
(15)
9i.
9a
In practice, we do not need to solve n linear systems. Instead, the solutions for g: A = 1,.......... .............................. n have simpler recurrence relation similar to Neville's method [36] in Lagrange polynomial interpolation.
Denote ig == 8 so that 6n- = Xio A3 = O, we have
Theorem 3.4 (Pseudo-order solver). For each ki the solution Eq. (15) is g k! where
D 9h1: = 0.1,...,n
(k-1) D((
(16)
D( 141
61+2 61
Proof in Appendix B.3. Note that the pseudo-order solver of order > 2 no longer has the guarantee of n-th order of accuracy, which not so important when NFE small. In our experiments, we n-th order predictor, then combine it with n-th
mainly rely two combinations: when order corrector or (n + 1)-th pseudo-order corrector.
Half-corree calf-corrector for large guidance scales. find that corrector may have negative effects on the sample quality. We propose useful technique called half-corrector by using the corrector only in the time region 0.5. Correspondingly, the strategy that we use corrector at each step is called fill-corrector.
When the guidance scale is large in guided sampling, we
3.4 Implementation Details
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 X, the EMS 1 S by can be estimated by firstly drawing K (1024-4096) datapoints æy ~ p (xx) with 200-step DPM-Solver++ [32] and then analytically computing some terms related eø (detailed in Appendix C.1.1). In practice, we find both convenient and effective to choose the distribution of the dataset qu to approximate ph. Thus, without further specifications, we directly use samples from qo-
Estimating integrals of EMS. We estimate EMS on N (120 ~ 1200) timesteps > Aj.....,>...... Ajn and trapezoidal rule estimate the integrals in Eq. (12) (see Appendix 1.3 for the estimation error analysis). We also apply some pre-computation for the integrals to avoid extra computation costs during sampling, detailed Appendix C.1.2.
Adaptation guided sampling. Empirically, we find that within common range of guidance scales, can simply compute the EMS on the model without guidance, and it work for both unconditional sampling and guided sampling cases. See Appendix for more discussions.
3.5 Comparison with Existing Methods
By comparing with existing diffusion ODE solvers that are based 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 x, sx, by specific values. (2) Our first-order discretization can be seen as improved DDIM. See more details in Appendix A.
Experiments
In this section, we show that DPM-Solver-v3 can achieve consistent and notable speed-up for both unconditional and conditional sampling with both pixel-space and latent-space DPMs. We conduct extensive experiments on diverse image datasets, where the resolution ranges from 256. First, present the main results of sample quality comparison with previous state-of-the-an training- free
7

---

HouR's V
DMM
selvers2
(a) CIFARIO
(b) CIFARIO (EDM, Pixel DPM)
(c) -SU--eddroom
(ScoreSDE, Pixel DPM)
(Latent-Diffusion, Latent DPM)
Figure 3: Unconditional sampling results. We report the FID of the methods with different numbers of function evaluations (NFE), evaluated on 50k samples. We borrow the results reported in their original paper directly.
199i bulers
UFFCC
Soleen)
611
ImageNet-256
(b) CSSCC002014
(c) MS-COC02014
(Guided-Diffusee Pixel DPM)
(Stable-Diffusion, Latent DPM) (Stable-Diffusion Latent DPM) (Classifier-Free Guidance, 8 = 1.5) (Classifier-Free Guidance, & = 7.5)
(Classifier Guidance, = 2.0)
Conditional sampling results. We report the FIDI or MSEJ of the methods with different numbers of function evaluations (NFE), evaluated on 10k samples.
Figure 4:
methods. Then we illustrate the effectiveness of our method by visualizing the EMS and samples. Additional ablation studies are provided in Appendix G. On each dataset, we choose sufficient number of timesteps N and datapoints K for computing the EMS to reduce the estimation error, while the EMS can still be computed within hours. After we obtain the EMS and precompute the integrals involving them, there is negligible extra overhead in the sampling process. We provide the runtime comparison in Appendix E. We refer to Appendix for more detailed experiment settings.
4.1 Main Results
We present the results in 5 ~ 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-an 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 't show the results for other methods such DDIM [48]. PNDM [28]. since they have already been compared in previous works and have inferior performance. The quantitative results on CIFARIO [24] are listed Table 1, and more detailed quantitative results presented in Appendix F.
Unconditional sampling We first evaluate the unconditional sampling performance of different methods on CIFARIO [24] and LSUN-Bedroom [55]. For CIFARIO we use two pixel-space DPMs,
(a) CIFARIO (ScoreSDE)
(b) CIFARIO (EDM)
(c) MS-COC02014
(Stable-Diffusioon = 7.5)
Figure 5: Visualization of the EMS Lx, 8x, by w.r.t. À estimated on different models.
8

---

(a) DPM Solver++ [32] (FID 18.59)
(b) UniPC [58] (FID 12.24)
(c) DPM Solver v3 (Ours) (FID 7.54)
Figure 6: Random samples of Latent- Diffusion [43] on
[55] with only NFE = 5.
Table I: Quantitative results on CIFARIO [24]. We report the FID, of the methods with different numbers of function evaluations (NFE), evaluated on 50k samples. We borrow the results reported their original paper directly.
NFF
Method
Model
6
8
10
12
15
20
25
!DIS [56]
15.37
4.17
3.37
2.86
DPM-Solver++ [32]
28.53
13.48
5.34
4.01
4.04
3.32
2.90
2.76
ScoreSDE [51]
UniPC [58]
23.71
10.41
5.16
3,93
3.88
3.05
2.73
2.65
DPM-Solvvr--3
12.76
7.40
3.94
3.40
3.24
2.91
2.71
2.64
Heun's 2nd [21]
320.80
103.86
39.66
16.57
7.59
4.76
2.51
2.12
DPM-Solver++ [32]
24.54
11,85
4.36
2.91
2.45
2.17
2.05
2.02
EDM [21]
UniPC [58]
23.52
11.10
3.86
2.85
2.38
2.08
2.01
2.00
DPM-Solver-v3
12.21
8.56
3.50
2.51
2.24
2.10
2.02
2.00
one is based on ScoreSDE [51] which is widely adopted model by previous samplers, and another is the latent-
based EDM [21] which achieves the best sample quality. For L.SUN-eddroo, model [43]. We apply the multistep 3rd-order version for DPM-Solver++: space DPM-Solver-v3 by default, which performs best in the unconditional setting. For UniPC, UniPC and we report the better result of their two choices B (h) = h and (h) = e at each NFE. For DPM-Solver-v3,
tune the strategies of whether use the pseudo-order predictor/corrector at each NFE on CIFARIO, 10, and use the pseudo-order corrector on LSUN- -UUN-BBrroom.. As shown in Fig. find that DPM-SolPe PPM--Soler-v3 can achieve consistently better FID, which especially notable when NFE is 5~10.. For example, we improve the FID on CIFARIO 10 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. with around 12 minutes computing of the EMS, DPM-Solver-v3 converges to the FID of 3.06 with NFE, which approximately 60% % sampling cost of the previous best training- free method (20 NFE by UniPC).
On LSUN-Bedrood
Conditional sampling. We then evaluate the conditional sampling performance, which more widely used since it allows for controllable generation with user-customized conditions. We choose conditional settings, is classifier guidance pixel-space Suided-Diffussion [10] model trained on ImageNet-2566 dataset [9] with 1000 class labels as conditions; the other is classifier- free guidance latent-space Stable-Diffusic 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 [10]; and the latter the guidance scale of 1.5 (following the original paper) or 7.5 (following the official code) with prompts random selected from IS-COC02014 validation [26]. Note that the FID Stable-Diffusion samples saturates 15.0~16 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 generated latent code Ộ and the ground-truth solution xf (i.e., - xx" / D) to evaluate convergence, starting from the same Gaussian noise. We obtain x² by 200-step OPM-Solvert+, which is enough to ensure the convergence.

---

We apply the multistep 2nd-order version for DPM-Solvert+. UniPC and DPM-Solver-v3, which performs best in conditional setting. For UniPC, we only apply the choice (2 = ex which performs better than B₁ (h). For our 3 we use the pseudo-order corrector by default, and report the best results between using cefffcoorrcctofull-orreeee 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 extra corrector, performs even worse than DPM-Solver++ when NFE<IO on Stable-Diffusion (s = 7.5), With the combined effect the EMS and the half- corrector technique, we successfully outperform DPM-Solver++ in such case. Detailed comparisons can be found in the ablations in Appendix G.
on Stable-Diffusion (s = 7.5).
4.2 Visualizations of Estimated EMS
We further visualize the estimated EMS in Fig. 5. Since LA, 82, by are D-dimensional vectors, we average them over all dimensions to obtain a scalar. From Fig. 5, we find that 4x gradually changes from to near 0 as the sampling proceeds, which suggests we should gradually slide from data prediction to noise prediction. As for 8x, by they are more model-specific and display many fluctuations, especially for ScoreSDE model [51] on CIFARIO. Apart from the estimation error of the EMS, we suspect that it comes from the fluctuations of ( which caused by the periodicity of trigonometric functions in the positional embeddings of the network input. It's worth noting that the fluctuation 8x, by will not cause instability in our sampler (see Appendix J).
4.3
Visual Quality
We present some qualitative comparisons in Fig. 6 and Fig. 1. We can find that previous methods tend to have 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 large guidance scales (in the left-top image, a part of the castle becomes hill; in the bef-bootom image, the hill is translucent and the castle is hanging in the air), while ours won't t. Additional samples are provided in Appendix K.
Conclusion
We study the ODE parameterization problem for fast sampling of DPMs. Through theoretical analysis, we find a novel ODE formulation with empirical model statistics, which is towards the optimal one to minimize the first-order discretization error. Based on such improved ODE formulation, we propose a new fast solver named DPM-Solver-33, which involves a multistep predictor-correctt framework of any order and several techniques for improved sampling with few steps large guidance scale. Experiments demonstrate the effectiveness of DPM-Solver-v? in both unconditional conditional sampling with both pixel-space latent-space pre-trained DPMs, and the significant advancement of sample quality in 5~100 steps.
Limitations and broader impact Despite the great speedup in small numbers of steps. DPM-Solver- v3 still lags behind training-bast based methods and is not fast enough for real-time applications. Besides, we conducted theoretical analyses of the local error, but didn't 't explore the global design spaces, such as the design of timestep schedules during sampling. And commonly, there are potential undesirable effects that DPM-Solver-v3 may be abused to accelerate the generation of fake and malicious content.
Acknowledgements
This work was supported by the National Key Research and Development Program of China (No. 2021ZD0110502) NSFC Projects (Nos. 62061 136001, 62106123, 62076147, 119A2081, 61972224, 62 (2111020 120), Tsinghua Institute for Guo Qiang, and the High Performance Computing Center, Tsinghua University. J.Z. is also supported by the XPlorer Prize.
10

---

References
[1] Kendall Atkinson, Weimin Han, and David Stewart. Numerical solution of ordinary differential equations, volume 108. John Wiley & Sons, 2011.
[2] Fan Bao, Chongxuan Li, Jiacheng Sun, Jun Zhu, and Bo Zhang. Estimating the optimal covariance with imperfect mean in diffusion probabilistic models. In International Conference on Machine Learning, pages 1555-1584 PMLR, 2022.
[3] Fan Bao, Chongxuan Li, Jun Zhu, and Bo Zhang. Analytic-DPM: An analytic estimate of the optimal reverse variance in diffusion probabilistic models. In International Conference Learning Representations, 2022.
[4] Costas Bekas, Effrosyni Kokiopoulou, and Yousef Saad. An estimator for the diagonal of 2007.
matrix. Applied numerical mathematics,
[5] Mari Paz Calvo and César Palencia. A class of explicit multistep exponential integrators for semilinear problems. Numerische Mathematik, 102:367-381, 2006.
[6] Nanxin Chen, Yu Zhang, Heiga Zen, Ron Weiss, Mohammad Norouzi, and William Chan. Wavegrad: Estimating gradients for waveform generation. In International Conference on Learning Representations, 2021.
[7] Hyungjin Chung, Jeongsol Kim, Michael Thompson Mccann, Mare Louis Klasky, and Jong Chul Ye. Diffusion posterior sampling for general noisy inverse problems. In The Eleventh Interna- tional Conference on Learning Representations, 2022.
[8] Guillaume Couairon, Jakob Verbeek, Holger Schwenk, and Matthieu Cord. Diffedit: Diffusion- based semantic image editing with mask guidance. In The Eleventh International Conference on Learning Representations, 2022.
[9] Jia Deng, Wei Dong, Richard Socher, Li-Jia Li, Kai Li, and Li Fei-Fei. ImageNet: A large- scale hierarchical image database. In 2009 IEEE Conference on Computer Vision and Pattern Recognition, pages 248-255. IEEE, 2009,
[10] Prafulla Dhariwal and Alexander Quinn Nichol. Diffusion models beat GANs image synthesis. In Advances in Neural Information Processing Systems, volume 34, pages 8780- 8794, 2021.
[11] Tim Dockhorn, Arash Vahdat, and Karsten Kreis. Genie: Higher-order denoising diffusion solvers. In Advances in Neural Information Processing Systems, 2022.
[12] Alfredo Eisinberg and Giuseppe Fedele. On the inversion of the vandermonde matrix. Applied mathematics and computation, 174(2):1384-139) 2006.
[13] Ian J. Goodfellow, Jean Pouget-Abadie Mehdi Mirza, Bing Xu, David Warde-Farley. Sherjil Ozair, Aaron C. Courville, and Yoshua Bengio. Generative adversarial nets. In Advances in Neural Information Processing Systems, volume 27, pages 2672-2680, 2014.
[14] Jonathan Ho, William Chan, Chitwan Saharia, Jay Whang, Ruiqi Gao. Alexey Gritsenko, Diederik Kingma, Ben Poole. Mohammad Norouzi, David J Fleet, et al. Imagen video: High definition video generation with diffusion models. arXiv preprint arXhv:2210.023033 2022.
[15] Jonathan Ho, Ajay Jain, and Pieter Abbeel. Denoising diffusion probabilistic models. In Advances in Neural Information Processing Systems, volume 33, pages 6840-6851, 2020.
[16] Jonathan Ho and Tim Salimans, Classifier-free diffusion guidance. In NeurIPS 2021 Workshop on Deep Generative Models and Downstream Applications, 2021.
[17] Marlis Hochbruck and Alexander Ostermann. Explicit exponential Kunge-Kuttar methods for semilinear parabolic problems. SIAM Journal on Numerical Analysis, 1369):0069-1090 2005.
[18] Marlis Hochbruck and Alexander Ostermann, Exponential integrators. Acta Numerica, 19:209- 286, 2010.

---

[19] Marlis Hochbruck, Alexander Ostermann, and Julia Schweitzer. Exponential rosenbrock typp methods. SIAM Journal on Numerical Analysis, 47(1):786-880) 803, 2009.
[20] Michael Hutchinson. A stochastic estimator of the trace of the influence matrix for laplacian smoothing splines. Communications in -itiitics--imlatiion and Computation, 18(3): 1059- 1076, 1989.
[21] Tero Karras, Miika Aittala, Timo Aila, and Samuli Laine. Elucidating the design space of diffusion-based generative models. In Advances in Neural Information Processing Systems, 2022.
[22] Bahjat Kawar, Michael Elad, Stefano Ermon, and Jiaming Song. Denoising diffusion restoration models. In Advances Neural Information Processing Systems, 2022.
[23] Diederik Kingma, Tim Salimans, Ben Poole, and Jonathan Ho. Variational diffusion models. In Advances in Neural Information Processing Systems, 2021.
[24] Alex Krizhevsky, Geoffrey Hinton, et al. Learning multiple layers of features from tiny images. 2009.
[25] Shengmeng Li, Luping Liu, Zenghao Chai, Runnan Li, and Xu Tan. Era-solver: Error- robust adams solver for fast sampling of diffusion probabilistic models. arXiv preprint Stti::2301.19335 2023.
[26] Tsung-Y Yi Lin, Michael Maire, Serge Belongie, James Hays, Pietro Perona, Deva Ramanan, Piotr Dollár, and C Lawrence Zitnick. Microsoft coco: Common objects in context. In Computer Vision-ECCV 2014: 13th Europeun Conference, Zurich, Switzerland, September 6-12, 2014, Proceedings, Purt 13, pages 740-755. Springer, 2014.
[27] Jinglin Liu, Chengxi Li, Yi Ren, Feiyang Chen, and Zhou Zhao. Diffsinger: Singing voice synthesis via shallow diffusion mechanism. In Proceedings of the AAAI Conference on Artificial Intelligence, volume 36, pages 11020- 11028, 2022.
[28] Luping Liu, Yi Ren, Zhijie Lin, and Zhou Zhao. Pseudo numerical methods for diffusion models on manifolds. In International Conference on Learning Representations, 2021.
[29] Xingehao Liu, Chengyue Gong, et al. Flow straight and fast: Learning to generate and transfer data with rectified flow. In The Eleventh International Conference on Learning Representations, 2022.
[30] Cheng Lu, Kaiwen Zheng, Fan Bao, Jianfei Chen, Chongxuan Li, and Jun Zhu. Maximum likelihood training for score-based diffusion odes by high order denoising score matching. In International Conference on Machine Learning, pages 14429-14460. PMLR, 2022.
[31] Cheng Lu, Yuhao Zhou, Fan Bao, Jianfei Chen, Chongxuan Li, and Jun Zhu. Dpm-solver: A fast ode solver for diffusion probabilistic model sampling in around 10 steps. In Advances in Neural Information Processing Systems, 2022.
[32] Cheng Lu, Yuhao Zhou, Fan Bao, Jianfei Chen, Chongxuan Li, and Jun Zhu. Dpm- solver++: Fast solver for guided sampling of diffusion probabilistic models. arXiv preprint 1:Xiv:221..01095 2022.
[33] Eric Luhman and Troy Luhman. Knowledge distillation in iterative generative models for improved sampling speed. arXiv preprint arXiv:2101.02388. 2021.
[34] Chenlin Meng, Ruiqi Gao, Diederik P Kingma, Stefano Ermon, Jonathan Ho, and Tim Salimans. On distillation of guided diffusion models. In NeurlPS 2022 Workshop Score-Based Methods, 2022.
[35] Chenlin Meng, Yang Song, Jiaming Song, Jiajun Wu, Jun-Yan Zhu, and Stefano Ermon. SDEdit: Image synthesis and editing with stochastic differential equations. In International Conference on Learning Representations, 2022.
[36] Eric Harold Neville. Iterative interpolation. St. Joseph's ' IS Press, 1934.
12

---

[37] Alexander Quinn Nichol and Prafulla Dhariwal. Improved denoising diffusion probabilistic models. In International Conference on Machine Learning, pages 8162-8171. PMLR, 2021.
[38] Alexander Quinn Nichol, Prafulla Dhariwal, Aditya Ramesh, Pranav Shyam, Pamela Mishkin, Bob Megrew, Ilya Sutskever, and Mark Chen. Glide: Towards photorealistic image genera- tion and editing with text text-guided diffusion models. In International Conference on Machine Learning, pages 16884-16804. PMIR, 2022.
[39]
Adam Paszke, Sam Gross, Francisco Massa, Adam Lerer, James Bradbury, Gregory Chanan, Trevor Killeen, Zeming Lin, Natalia Gimelshein, Luca Antiga, Alban Desmaison, Andreas Kopf, Edward Yang, Zachary Viti,o, Martin Raison, Alykhan Tejani, Sasank Chilamkurthy, Benoit Steiner, Lu Pang, Junjie Bai, and Soumith Chintala. Pytorch: An imperative style, high- performance deep learning library. In H. Wallach, H. Larochelle, A. Beygelzimer, F. d'Alché- Bue, E. Fox, and R. Garnett, editors, Advances in Neural Information Processing Systems, volume 32, Curran Associates, Inc., 2019.
[40] Vadim Popov, Ivan Vovk, Vladimir Gogoryan, Tasnima Sadekova, Mikhail Kudinov, and Jiansheng Wei. Diffusion-based voice conversion with fast maximum likelihood sampling scheme. In International Conference on Learning Representations, 2022.
[41] Alec Radford, Jong Wook Kim, Chris Hallacy, Aditya Ramesh, Gabriel Goh, Sandhini Agarwal, Girish Sastry, Amanda Askell, Pamela Mishkin, Jack Clark, et al. Learning transferable visual models from natural language supervision. In International conference on machine learning. pages 8748 8763. PMLR, 2021.
[42] Aditya Ramesh, Prafulla Dhariwal, Alex Nichol, Casey Chu, and Mark Chen. Hierarchical text-conditional image generation with CLIP latents. arXiv preprint arXiv:2204.0611 06122......, 2022.
[43] Robin Rombach, Andreas Blattmann, Dominik Lorenz, Patrick Esser, and Björ Ommer. High- resolution image synthesis with latent diffusion models. In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 10684-10699 2022.
[44] Chitwan Saharia, William Chan, Saurabh Saxena, Lala Li, Jay Whang, Emily Denton, Seyed Kamyar Seyed Ghasemipour, Raphael Gontijo-Lopes Burcu Karagol Ayan, Tim Salimans, et al. Photorealistic text-to-image diffusion models with deep language understanding. In Advances in Neural Information Processing Systems, 2022.
[45] Tim Salimans and Jonathan Ho. Progressive distillation for fast sampling of diffusion models. In International Conference on Learning Representations, 2022.
[46] Christoph Schuhmann, Romain Beaumont, Richard Vencu, Cade W Gordon, Ross Wightman, Mehdi Cherti, Theo Coombes, Aarush Katta, Clayton Mullis, Mitchell Wortsman, et al. Laion- 5b: An open large-scale dataset for training next generation image-text models. In Thirty-sixth Conference on Neural Information Processing Systems Datasets and Benchmarks Track, 2022.
[47] Jascha Sohl-Dicke Dickstein, Eric Weiss, Niru Maheswaranathan, and Surya Ganguli. Deep unsuper- vised learning using nonequilibrium thermodynamics. In International Conference on Machine Learning, pages 2256 2265. PMLR, 2015.
[48] Jiaming Song. Chenlin Meng, and Stefano Ermon. Denoising diffusion implicit models. In International Conference on Learning Representations, 2021.
[49] Yang Song, Prafulla Dhariwal, Mark Chen, and Ilya Sutskever. Consistency models. arXiv preprint arXiv:2303.01446, 2023.
[50] Yang Song, Conor Durkan, Iain Murray, and Stefano Ermon. Maximum likelihood training of score-based diffusion models. In Advances in Neural Information Processing Systems, volume 34, pages 1415-1428, 2021.
[51] Yang Song, Jascha Sohl-Dickstein, Diederik P. Kingma, Abhishek Kumar, Stefano Ermon, and Ben Poole. Score-based generative modeling through stochastic differential equations. In International Conference on Learning Representations, 1021.
13

---

[52] Hideyuki Tachibana, Mocho Go. Muneyoshi Inahara, Yotaro Katayama, and Yotaro Watanabe. Itô-taylor sampling scheme for denoising diffusion probabilistic models using ideal derivatives. arXiv e-prints, pages arXiv-2112, 2021.
[53] Daniel Watson, William Chan, Jonathan Ho, and Mohammad Norouzi. Learning fast samplers for diffusion models by differentiating through sample quality. In International Conference on Learning Representations, 2022.
[54] Suttisak Wizadwongsa and Supasorn Suwajanakom. Accelerating guided diffusion sampling with splitting numerical methods. In The Eleventh International Conference on Learning Representations, 2022.
[55] Fisher Yu, Ari Seff, Yinda Zhang, Shuran Song, Thomas Funkhouser, and Jianxiong Xiao. LSUN: Construction of scale image dataset using deep learning with humans in the loop. arXiv preprint arXiv: 5006 03360, 2015.
[56] Qinsheng Zhang and Yongxin Chen. Fast sampling of diffusion models with exponential integrator. In The Eleventh International Conference on Learning Representations, 2022.
[57] Min Zhao, Fan Bao, Chongxuan Li. and Jun Zhu. Egsde: Unpaired image-to-image transla- tion via energy-guided stochastic differential equations. Advances Neural Information Processing Systems, 2022.
[58] Wenliang Zhao, Lujia Bai, Yongming Rao, Jie Zhou, and Jiwen Lu. UniPC: A unified predictor- corrector framework for fast sampling of diffusion models. arXiv preprint arXiv2230.004867, 2023.
[59] Kaiwen Zheng, Cheng Lu, Jianfei Chen, and Jun Zhu. Improved techniques for maximum likelihood estimation for diffusion odes. arXiv preprint 2305:22055.0933 2023.
14

---

A Related Work
Diffusion probabilisticc models (DPMs) [47, 15, S1], also known score-based generative models (SGMs), have achieved remarkable generation ability on image domain [10, 21J, yielding extensive applications such as speech, singing and video synthesis [6, 27, 14], controllable image genera- tion, translation and editing [38, 42, 43, 35, 57, 8], likelihood estimation [50, 23, 30, 59], data
compression [23] and inverse problem solving [7. 22].
A.1
Fast Sampling Methods for DPMs
Fast sampling methods based on extra training or optimization include learning variances in the reverse process [37, 2], learning sampling schedule [53], learning high-order model derivatives [11], model refinement [29] and model distillation [45, 33, 49]. Though distillation-ban
methods can generate high-quality samples in less than 5 steps, they additionally bring onerous training costs. Moreover, the distillation process will inevitably lose part of the information of the original model, and is hard to be adapted to pre-trained large DPMs [43, 44, 42] and conditional sampling [34]. Some of iiiillation-bbsed methods also lack the ability to make flexible trade-off's between sample speed and sample quality.
In contrast, training- free samplers are more lightweight and flexible. Among them, samplers based on diffusion ODEs generally require fewer steps than those based on diffusion SDEs [51, 40, 3, 48, 32], since SDEs introduce more randomness and make the denoising harder in the sampling process. Previous samplers handle the diffusion ODEs with different methods, such as Heun's methods [21], splitting numerical methods [54], pseudo numerical methods [28], Adams methods [25] or exponential integrators [56, 31, 32, 58].
A.2 Comparison with Existing Solvers Based on Exponential Integrators
Table 2: Comparison between DPM-Solver-v3 and other high-order diffusion ODE solvers based on exponential integrators.
DEIS [56]
DPM-Solver
DPM-Solver++
UniPC [58]
DPM-Solver-v3 (Ours)
[31]
[32]
First-Order
DDIM
DDIM
DDIM
DDIM
Improved DDIM ge for À Multistep
Taylor Expanded Predictor
69 for Multistep
eg for Singlestep
Ty for X Multistep
3g for À Multistep
Solver Type (High-Order) Applicable for Guided Sampling Corrector Supported Model-Specific
V
X
/
X
X
X
V
X
X
X
X
In this section, we make some theoretical comparisons between DPM-Solver-v3 and existing diffusion ODE solvers that are based on exponential integrators [56, 31, 32, 58]. We summarize the results in Table and provide some analysis below.
Previous ODE formulation as special cases of ours. First, we compare the formulation of the ODE solution. Our reformulated ODE solution in Eq. (9) involves extra coefficients I,, 8x, bs. and it corresponds to a new predictor ge to be approximated by Taylor expansion. By comparing ours with previous ODE formulations Eq. (3) and their corresponding noise/data prediction, we can easily figure out that they are special cases of ours by setting l,, 8x, by to specific values:
Noise prediction: 1 0, 8 -
0
Data prediction: L 1, 8) 0, b 0
It's worth noting that, though our ODE formulation can degenerate to previous ones, it may still result in different solvers, since previous works conduct some equivalent substitution of the same order in the local approximation (for example, the choice of B, (h) or B2(h) = el 1 in UniPC [58]).
We never conduct such substitution, thus saving the efforts to tune it.
15

---

Moreover, under our framework, we find that DPM-Solver++ is a model-agnostic approximation of DPM-Solver-v3, under the Gaussian assumption. Specifically. according Eq. (5), we have
X argmin E ((x)) (ax)) Vaco(2x, - 11 12
(17)
If we assume (xx)= = Maxslaxa I for some fixed Fo. then the optimal noise predictor is
Es(2), x - Vv (xx))((x)) 2AA 0r0200
(18)
σx
follows that a (= ~ thus 1 R 1 by Eq. (5), which corresponds the data prediction model used in + Moreover, for small enough À (i.e., near to T), the Gaussian assumption is almost true (see Section 4.2), thus the data-prediction DPM-Solver++ approximately computes all the linear terms at the initial stage. To the best of knowledge, this is the first explanation for the reason why the data-prediction DPM-Solver++ outperforms the noise-prediction DPM-Solver.
First-order discretization as improved DDIM Previous methods merely use noise/data parameteri- zation, whether or not they change the time domain from to X. While they differ in high-order cases, they are proven to coincide in the first-order case, which is DDIM [48] (deterministic case, '/ = 0):
at C. a2
0=
0z ar
Eg(
(19)
Oks
Os
However, the first-order case of our method is
(1+th, Ln, /"'
d
L/m 0s
Ar Ex, (J)d)) coller)he)
Ot 2. A(dor as
Ex. (N)
E.
(20)
- A(
Ex, (A)
which is DDIM since we choose better parameterization by the estimated EMS. Empirically, our first-order solver performs better than DDIM, as detailed in Appendix G.
B Proofs
B.1
Assumptions
In this section, we will give some mild conditions under which the local order of accuracy of Algorithm and the global order of convergence of Algorithm (predictor) are guaranteed.
B.1.1 Local
First, we will give the assumptions to bound the local truncation error.
Assumption B.1. The total derivatives of the noise prediction model d' d' k = 1,..., exist and are continuous.
d
Assumption B.2. The coefficients Lx, 8x, by are continuous and bounded. 1...., exist and are continuous.
d dxF:
d²
d'bs
axF k
dxF:
Assumption B.3. & O X.), 1,.
Assumption B.1 is required for the Taylor expansion which is regular in high-order numerical methods. Assumption B.2 requires the boundness the coefficients as well as regularizes the coefficients' smoothness to enable the Taylor expansion for (xx, 2), which holds in practice given the smoothness of ea ( A) and p (æx). Assumption B.3 makes sure dn and X X, is of the same order, i.e., there exists some constant fik = O(1) so that òn = TA (Xe As), which is satisfied in regular multistep methods.
16

---

B.1.2 Global
Then we will give the assumptions to bound the global error.
Assumption B.4. The noise prediction model E, is Lipschitz w.r.t. to
Assumption B.5. h maxigisM(A; -A--1 (L/MM).
Assumption B.6. The starting values æ,, 1 S i - satisfies ã. -
Assumption B.4 is common in the analysis of ODEs, which assures () ( Assumption B.5 implies that the step sizes are rather uniform. Assumption B.6 common in the convergence analysis of multistep methods [5].
B.2
Order of Accuracy and Convergence
In this section, we prove the local and global order guarantees detailed Theorem 3.1 and Theo- rem 3.3.
B.2.1
Local
Proof. (Proof of Theorem 3.1) Denote h == d >a. Subtracting the Taylor-expanded exact solution in Eq. (12) from the local approximation in Eq. (14), we have
(k) 9p (x),, 2a) 9o
(k) (ax,
No)
A: -A
8 - A A
k!
1m
(21)
First we examine the order of A(),, Ae) and 1.: Ex, (A)(A dd Under Assumption B.2, there exists some constant C₁, C2 such that -lx < C1,,4 + 82 < C2. So
12 1.dr
A()
< Cih
(22)
= 011)
1/m -2)(2
frm de eB (1 (A 2o)
(23)
(2 - -2o)
O(h+t))
Next we examine the order (1x since Do is elementary function of Eg and lx SA, by we know / are continuous. Adopting the notations in Eq. (13), by Taylor expansion, we have
(k) (w)) Under Assumption B.1 and Assumption B.2, ( (zx, A.), k I,..., n exist and
9:: 9. + 01g (1) + g. + +... (m) + 0(8++))
9:, g. + 82g (1) + 8gg (2) + +... +.... 1
(24)
91., 9, + 5. (1) + gg
(2) + 6 + 0(6++1)
Comparing it with Eq. (13), we have
8
8th
(
82
8
g g
o ( +1
2
(25)
8 n
5
(8
n
17

---

From Assumption B.3, we know there exists some constants so that 6% = rkh, h = =1, I..... n. Thus
09
r1
r r
1
(ha+l
82
5%
r2
12
2
h²
( : (
O
:
6
r
h
(
(26)
And finally, have
g
h
r1
r
Ohhat2)
h-
2
O(h+t2)
T2
12
r2
9 (n)
9>
(1)
n
Tn
((hat2)
(27)
O(h))
O(h-11)
(
Substitute Eq. (22), Eq. (23) and Eq. (27) into Eq. (21), can conclude that 2. O(h+t2).
B.2.2 Global
First, we provide a lemma that gives the local truncation error given inexact previous values when estimating the high-order derivatives.
Lemma B.7. (Local truncation error with inexact previous values) Suppose inexact values T k 1...., n and C, are used in Eq. (13) to estimate the high-order derivatives, then the local truncation error of the local approximation Eq. (14) satisfies
A a,A(),, Ad) A + 0(L) (O(A.) Σ +((x))+(()t))
a
(28)
Og
where A
C x. h = d;- >..
Proof. By replacing æ. with in Eq. (13) and subtracting Eq. (12) from Eq. (14), the expression for the local truncation error becomes
A A(s,, de) A. a.A(A,, Ar) (go ( As) -00 Ao)) / dr
EA,AA
Oss
(29)
n
9o
(k)
(2),, 2o)
90 (k)
(xX. Ao) / pd
(J((A X,)=x -0h
k!
1=1
And the linear system for solving 90
(k) (xx, Aa), k
n2 becomes
61
8
87
9
g.
82
9i,
9,
2
(30)
8
9i.. 9
where 9. _ X.). Under Assumption B.4, we know 9. 9 O(Ax). Thus, under Assumption B.I, Assumption B.2 and Assumption B.3, similar to the deduction the last section, we have
9.
h
T1
r
O( + A + Ax.,) 0(ha+1 + A + Axiz)
g (2)
h2
r2
T2
n
2
(m)
r2
r/
0 + A. + AM...
ml
(31)
18

---

Besides, under Assumption B.2, the orders of the other coefficients are the same we obtain in the last section:
=
Jx.
(,a 1(h
(32)
Thus
xx,, Aa) 99
(k) (oXx Ao)
A:
(A)(AAA-o)
k!
Ex. (A)AA dA
(1)
99
9
[xx
Ex, (A)( ddA
(
2
IN EA, (A)(A Xd)'dA
9
O(P3)
h--
O( + A, + +Anr
O(13)
h-
12
O + A + Arin
)(
h-n
Tn
In
0 + D, + A
O(h)
F1 x2
0Oth+1 A. Axi,
O(h)
12
O(hat A A
O(h))
F'
O( + A, + AA...)
O(h)O(h** + + Ax,,)
(33)
Combining Eq. (29), Eq. (32) and Eq. (33), we can obtain the conclusion in Eq. (28).
Then prove Theorem 3.3 below.
Proof. (Proof of Theorem 3.3)
As we have discussed, the predictor step from 1m to fm is a special case of the local approximation (fm-n-1:- nl 1...., f, fmm-2. fm- (m). By Lemma B.7 we have
Eq. (14) with ....:...,,,,....
Ot A A(Aem dem)
Am- + O(h)
Σ 20 +
(34)
O't,.
follows that there exists constants C, Co irrelevant to h, so that
Am S
attm
A(Xtm-11 Aem
Ch) Am-1l Ch Σ
Am Coh
(35)
atim11
k=0
Denote fm = maxosism Adl we then have
O't A(dim-1s dem) + Cih) fm-1 + Coh
A
(36)
Otmm-
Since atn A(Aem dem)
when h + 0 and it has bounded first-order derivative due to As-
atm-1
sumption B.2, there exists a constant C2, so that for any C. 2 C2. atm A(he m- prem) Ch > 1 for
am
sufficiently small h. Thus, by taking C3- max(C1, C2), we have
O't A(Xem 11 dem) +Cgh)) fm-1 Coh
fm <
(37)
Oltm-1
19

---

orn Aldm
Denote Am-1 =
Czh, by repeating Eq. (37), we have
alm-1
M
M
M
fM
II
A fn
Σ 1=n+1
I
A
Cohn+2
(38)
By Assumption B.5, h (1/M), have
M-1
Oesy A(drn dese
M-1
On, a
II 1
A
1+
Cte A A(N A1,))
at
==
Ct (een+ dese)
M-1
0tt-1 Ci
1+ at A( A
(39)
ata
im
Ctr
(e desel
(1
M-n
ala
M
I Csea
ac,-1 Co
at A(de,iiza))
where maxnsigM-1 on A(Ar;-1 (
Then denote B = maxp+1SiSM
have
Σ
M-1
II A
Σ
Otay
A(de, dewe (1+
M-i
at2
M
1=
ix+1
M-n-1
3
Σ
+
(40)
M
1=0
3M
(1
M-n
(1+
M
S C (E - 1) M
Then substitute Eq. (39) and Eq. (40) into Eq. (38). Note that M = O(1/h) by Assumption B.5, and f. = Il by Assumption B.6, finally we conclude that A S fM 1(
B.3 Pscudo-Order Solver
First, we provide lemma that gives the explicit solution to Eq. (15),
Lemma B.S. The solution to Eq. (15)
(k)
g
gip gio
(41)
k!
(6p 8g)
0
Proof. Denote
51
82
82
82
85
R²
(42)
6
8%
Then the solution to Eq. (15) can be expressed as
g
Σ (R pm1
(R (gn - -
(43)
hd
where (R kp is the element of R at the k-th row and the p-th column. From previous studies of the inversion of the Vandermonde matrix [12], we know
RR
kpp
(44)
6 11== (6p 5,)) (6, 50) 11z=1
(6p 5,)
Substituting Eq. (44) into Eq. (43), we finish the proof.
20

---

Then we prove Theorem 3.4 below:
Proof. (Proof of Theorem 3.4) First, mathematical induction prove that
k
D( = -1
9
g
1 k:
(45)
0=000p 0utp 6149)
For = 1. Eq. (45) holds by the definition of ( Suppose the equation holds for ki we then prove it holds for k +
Define the Lagrange polynomial which passes (6i+p: Putp
9i:) for p k:
Plu = Σ giitp g4.
II q=0,p#p
1 < n, snsnssnn-
(46)
ôitp
81+9
Then 1/44) (x) is the coefficients before the highest-order term x² in P(A)(() We then prove that (())(a) satisfies the following recurrence relation:
(zx*)
(x)= (x) (x - P1) (x) - (x- 61th ( (p(
(47)
81+b - 81
By definition, P p(k-1)) (x) is the (k 1) order polynomial which passes (Ōi+p: 9u+r 9.)) for
1 S ps & k, and P. (x) is the (k - 1h order polynomial which passes (ōi-ps Putt - gin for
k 1.
Thus, for S p k have
(81=p 5,) p)) (b(+p))
(ōttp P( (81tp)
p(()) (8(+p)
-ii+
gi:
(48)
61tk 51
For p 0, have
p(k) (8) (ōr 61) (8t) ()6)) (ō) ((k-))
= In -
(49)
8144 - 81
for p k. we have
p(O) (Br4x) (6L+x (
(6(+k) (5L+k P(t girik
8141 81
g:.
(50)
Therefore, p(()) (a)) is the k-th order polynomial which passes + distince points (51+p- Pi-p 9
for 0 I Sp S k. Due to the uniqueness of the Lagrange polynomial, we can conclude that p (kk
(x:)
((k) (x). By taking the coefficients of the highest-order term, we obtain
((--1) D((-1)
D(M) D41
(51)
814b 81
where by the induction hypothesis we have (k-1) (kk-1) ",(="",((*")) = D(k-1) Comparing
Eq. (51) with the recurrence relation of ((A) in Eq. (16), follows that (()) = -(()) which
completes the mathematical induction.
Finally, by comparing the expression for (((0) in Eq, (45) and the expression for (() in Lemma B.8, (k) (k)
conclude that g:
21

---

B.4 Local Unbiasedness
Proof. (Proof of Theorem 3.2) Subtracting the local exact solution in Eq. (9) from the (n + 1)-th order local approximation in Eq. (14), we have the local truncation ertor
A Ex, (2)g A)dA
go (I Xo) L Ex, (A) (A-12)) 2a)
Σ
dd
h!
= -As((A))A
(x, (20 - go -3) dd
A: 9i (ax, (o) /.A
Ex, EA, (A- Aa)A (1.)
- -a Nx)
k!
1=1
=o,((A,,Ar)
(g) go Aa)) dd
(alg And go go (ol)) L./ (2)
(2 -2)+ 2a)a
- -,,((H,
Σ (
dA
k!
(52)
where œ is on the ground-truth ODE trajectory passing 2x,, and (R,, )a is the element of the inverse matrix R, the k-th row and the I-th column, as discussed in the proof of Lemma B.8. By Newton-Leibniz theorem, we have
/ A
go ( A) - go , Aa) ) = 9o (x,,, T)dT
(53)
, Also, since ,
I,.... . are on the ground-truth ODE trajectory passing æx,, have
Au
go ( Az,) - 220
I9 (x,,T)dr
(1)
(54)
where
(1) (I., T) = Sx, endr (2)) T) - 8, fzfo T) - b,) 90
(55)
Note that sx.l the solution to the least square problem in Eq. (11), which - fo( b,
makes sure (xp((r,) (1))) (x,, -
0. It follows that
(ex,)
(exx,) ( (IT,, T) 8 s-fo(a;,; by
O, since is on the ground-truth ODE
trajectory passing Øx,- Therefore, we have Bpg, (mx,) A) (o( = 0 and
Egg, (mx,)
go(x An)) - go(2x, Ae) = 0. Substitute them into Eq. (52), we conclude that
Epg, (mxx)
L. - xc) =0.
C
Implementation Details
C.1 Computing the EMS and Related Integrals in the ODE Formulation
The ODE formulation and local approximation require computing some complex integrals involving lx, x,, bx In this section, we'll give details about how estimate 1 8 bš, on few datapoints, and how to use them to compute the integrals efficiently.
C.1.1 Computing the EMS
First for the computing of 1 in Eq. (5), note that VaNo(Ex. = αλ Vec(Ex diag(lA)
(56)
Since diag(Ix) is diagonal matrix, minimizing E (mx)
No A) equivalent to mini-
mizing (xg(sr) (V (VaNNg(xx, x) = (x)(
||diag
- (ax V A)) 1x [1]] where
22

---

diag" - denotes the operator that takes the diagonal of a matrix as a vector. Thus we have K = = ((x)))) (ax (a V A))]-
However, this formula for 1 requires computing the diagonal of the full Jacobian of the noise prediction model, which typically has (d²) time complexity for d-dimensional data and is unaccept- able when d is large. Fortunately, the cost can be reduced to (dd) by utilizing stochastic diagonal estimators and employing the efficient Jacobian-vecte vect--o-ooor--dueett operator provided by forward- -ode automatic differentiation in deep learning frameworks.
For d-by-d matrix D, its diagonal can be unbiasedly estimated by [4]
diag -1 -(D)
Σ (Db) O Uk
Σ
2k 2k
(57)
where p(p) are d-dimensional i.i.d. samples with zero mean, @ is the element-wise multiplica- tion i.e., Hadamard product, and @ is the element-wise division. The stochastic diagonal estimator is analogous to the famous Hutchinson's 's trace estimator [20]. By taking p(v) as the Rademacher distribution, we have dy C 2 = 1, and the denominator can be omitted. For simplicity, we use regular multiplication and division symbols, assuming they are element-wise between vectors. Then Ix can be expressed as:
K = = 1(xx Ea Ca A)u)al
(58)
which an unbiased estimation when we replace the expectation with mean on finite samples TA p (=x)) ~ p(v). The process for estimating I can easily be paralleled on multiple devices by computing [oxV (x A)v) u on separate datapoints and gather them the end.
Next, for the computing of x by in Eq. (11), note that it's a simple least square problem. By taking partial derivatives w.r.t. 82, by and set them O, we have
(Ox) (Fo (xx, sifo(xx, A) -b;) bb; fo(zx. A) E ((x) f6 (2x, - (z)o 1) -bi =0
(59)
And we obtain the explicit formula for sx by
Bog(oN)) fo(2x,, A) (2) (zx,, 2) E (mx) [fo A)] (ox) (ex)
sx Ex( So (x A) (x (2A, A)] Bpg((r)) (2xx (a B (REa) (2(x
(60)
6y- Bp((=))) xxx A)] B (o) A)]
(61)
which are unbiased least square estimators when we replace the expectation with mean on finite samples 2y ~ pxx (xx). Also. the process for estimating 8x 6y can be paralleled on multiple devices
by computing fa, Ese") Soso, E.f on separate datapoints and gather them the end. Thus, the estimation of s3 by involving evaluating fo and f(") æ fo is direct transformation of Eg and requires a single forward pass. For ((L))
we have
f(2)) (xx,,) = @fo(ox,1) A)
dx dA
+ V A)
OX
= (1) (2x, co(2x),))
isas - -irs
la
( Ox xs - -ox0((0x:11
af
as
as
= ((x(((1))00(xs + t (
aA
(62)
After we obtain 4, i can be estimated by finite difference. To compute (1) (I), A), we have
dæy dA
6 (2x, A) = Oxco(xx, A) co ))
(63)
( xx - oxxe(((())
as
= = A) + VeEa(EX,)1)
αx
23

---

which can also be computed with the Jacobian-vecton secptn-eeettt---roeoe pector-product operator.
In conclusion, for any X, 5x 8 by can be efficiently and unbiasedly estimated by sampling few datapoints ω p (x)) and using the Jacobian-
C.1.2 Integral Precomputing
In local approximation in Eq. (14), there are three integrals involving the EMS, which are A(),, Ae). de 1xx Ex, (A) Σx, (A) (A-A.)* (A)(~>>)) dA. Define the following terms, which also evaluated Asa: Ap..... Aix and can be precomputed in ((V) time:
L
1,
Ar
{ 8,d7
SA
f b,
d
Bx
e
-Jn e-dr
b
(64)
Ar
Cx (Fiy((++a.)))
A
e -Jiy Big a
a b,dr du fod el Budu
AA
dr Ir cl
I (f
Then for any Asy de, we can verify that the first two integrals can be expressed as
A()o,de) cL,,, Lx.
(65)
EA, (A)d - (Cx, -C2- (Cx,
(Ix - Ix.))
which can be computed in ((1) time. For the third and last integral, denote it (K) ie.,
((k)
S
A Ex,
(A) )a)* Yd
b, Ae =
(66)
k!
Au
We need to compute it for < S n and for every local transition time pair (A Ac) in the sampling process. For & = O, have
E (0) Anide
L., Sx, (Ix, Ixa)
(67)
which can also be computed in O(1) time. But for k >> 0. we no longer have such a simplification technique. Still, for any fixed timestep schedule { during the sampling process, we can use
lazy precomputing strategy: compute (k)
1,A,1 1 < i M when generating the first sample, store it with unique key (k, and retrieve in ((1) in the following sampling process.
C.2 Algorithm
We provide the pseudocode the local approximation and global solver in Algorithm and Algo- rithm 2, which concisely describes how we implement DPM-Solver-v3.

---

Algorithm (n + 1)-th order local approximation: LUpdaten+11
Require: noise schedule ae, dt, coefficients L2, Sx, by
Input: transition time pair (s, 1), Vs. n extra timesteps {Lig Ik=1 go values (gi.. ................................ 9.) at
(axx, by J== and (æa, 6)
Input Format: {ts, I..
{ts,: 9t, {8, Te. got.t Ex, (A) (A)dA, Sr. Ex, ( dd (Appendix C.1.2)
1: Compute A(,,, Az), S
A.
2: 6k
Xxx
k
9
81
82
81
-1
9h1
9.
g (2)
82
0
912
9:
3:
2!
(Eq. (13))
8n
9i. 9
X1
Σ 9 Jo A
12
--
4: Zt a,e()e,
Ex, (A) (A)dA
Ex,
di
a3
h!
(Eq. (14))
Output: 2t
Algorithm 2 (n + 1)-th order multistep predictor-corrector algorithm
Require: noise prediction model 6g. noise schedule 02, ot coefficients L, 8s, bu. cache Q1, Q2 Input: timesteps (44 J.M initial value 22
1: Q1 caphe
2: Q2 < EB (00, 40)
cache
3: for
to M do
4: 7m < min (n + 1, m)
fm fetch f
5: 4m nm
Q1
ém fetch
6: ém nm
Q2
91 Shin dr ON - 11 tt
dr
e -im- o.dr bydr, 1= m Mmms 1
7:
Ot);
Am-1
(Eq. (8))
8:
om + LUpdaten (m1m-mm Im- b.... (4m- 9m (1m-1 2m- 9m-1) lm)
9:
if m + M then
10:
Em 5- co c0(2m.tm) tm)
Im +- Bam d, Are Cm
Ix.m Lm
Am
- -hm-i 5rdT
11:
dr (Eq. (8))
a
Am-1
12:
te + LUpdate um (1tm nom +1 Im - am am +1
ftm Im 22n-2 1o,,ttte frmm 9m {m 1 B 1 9 fm- Lm)
13:
E +- Em + 1ym (#e - km)/ (to ensure g. = 9m)
Q1 coche 2i
Q2 cache ec
16: end if
end for
Output: &M
D Experiment Details
In this section, we provide more experiment details each setting, including the codebases and the
configurations for evaluation, EMS computing and sampling, Unless otherwise stated, we utilize the forward-mode mode automatic differentiation (torch, autograd autograd.ffrrard forward _d) provided by Py PyTorch [39] to compute the
(JVPs). Also, stated in Section 3.4, we draw datapoints
25

---

2 from the marginal distribution yA defined by the forward diffusion process starting from some data distribution 9o, instead of the model distribution p
D.1 ScoreSDE on CIFAR10
Codebase and evaluation For unconditional sampling on CIFARIO [24], one experiment set- ting based on the pretrained pixel-space diffusion model provided by ScoreSDE [51]. We use their official codebase of PyTorch implementation, and their checkpoint checkpoint_ 8. 9th under vp/cifarl __iffarr_dddpppppeeee________ cetppppddee__ootinoottee continuous config. We adopt their own statistie file and code for computing FID.
EMS computing We estimate the EMS N 1200 uniform timesteps Ajo+ A Ajss..., dix by drawing K = 4096 datapoints Xxo ~ 9o. where is the distribution of the training set. We compute two sets of EMS, corresponding to start time c 10-3 (NFE< 10) and 10 (NFE>10 10) in the sampling process respectively. The total time for EMS computing is ~ on 8 GPU cards of NVIDIA A40.
Sampling Following previous works [31, 32, 58]. we use start time = 10-3 (NFE< 10) and 10-1 (NFE>10); end time T and adopt the uniform logSNR timestep schedule. For DPM-Solver-v3, we use the 3rd-order predictor with the 3rd-order corrector by default. In particular, we change to the pseudo 3rd-order predictor 5 NFE to further boost the performance.
D.2 EDM on CIFAR1O
Codebase and evaluation For unconditional sampling on CIFARIO [24], another experiment setting based on the pretrained pixel-space diffusion model provided by EDM [21]. We use their official codebase of Py PoTcc implementation, and their checkpoint edm- cdm-cifaalln 1elffr10-32233 12x33-uncood--o- pk1. For consistency, we borrow the statistic file and code from ScoreSDE [51] for computing FID.
EMS computing Since the pretrained models of EDM are stored within the pickles, fail to torch. autograd. forward_ad ad for computing JVPs. Instead, torch. autograd. www.ciiona...ii......o ward trick. We estimate jvp. which is much slower since it employs the double back- of EMS. One corresponds to N 1200 uniform timesteps Ajo A and K = 1024 datapoints A ~ qo. where Yo is the distribution of the training set. The other corresponds to N 120, K 4096. They used when NFE<10 and NFE> 10 respectively. The total time for EMS computing is 3.55h on 8 GPU cards of NVIDIA A40.
Sampling Following EDM.
start time tini 0.002 and end time max 80.0, but adopt the uniform logSNR timestep schedule which performs better in practice. For --M--olver--3, we use the 3rd-order predictor and additionally employ the 3rd-order corrector when NFES 6. In particular, we change to the pseudo 3rd-order predictor at 5 NFE to further boost the performance.
D.3 Latent-Diffusion on
Codebase and evaluation The unconditional sampling on LSUN-Bedroom [55] is based the pretrained latent-space diffusion model provided by Latent-Diffusion [43]. We use their official codebase of PyTorch implementation and their default checkpoint. We borrow the statistic file and code from Guided-Difflusion [10] for computing FID.
EMS computing We estimate the EMS at N = 120 uniform timesteps A Ajs ........ Aix by drawing K = 1024 datapoints I>, ~ qu. where go is the distribution of the latents the training set. The total time for EMS computing is ~12min on GPU cards of NVIDIA A40.
Sampling Following previous works [58], we use start time e 10-3, end time 1 and adopt the uniform / timestep schedule. For DPM-Solver-v3,, we use the 3rd-order predictor with the pseudo 4th-order corrector.
D.4 Guided-Diffune on ImageNet-256
Codebase and evaluation The conditional sampling ImageNet-2566 [9] is based on the pre- trained pixel-space diffusion model provided by Guidee-Diffusson [10]. We use their official
codebase of Porch implementation and their two checkpoints; the conditional diffusion model
26

---

256x256_diffusion www.ffusiioonoooo and the classifier 256x256 pt. We adopt their own statistic file and code for computing FID.
EMS computing We estimate the EMS at N = 500 uniform timesteps jjo- >...... Aix by drawing K 1024 datapoints AA qo- where go is the distribution of the training set. Also, we find that the FID metric on the ImageNet-256 dataset behaves specially, and degenerated L (I), = 1) performs better. The total time for EMS computing is ~9.5h on GPU cards of NVIDIA A40.
Sampling Following previous works [31, 32, 58], we use start time = 10 33, end time T' = 1 and the 2nd-order predictor with the
adopt the uniform t timestep schedule. For JPM-Solver-v3, pseudo 3rd-order corrector.
D.5 Stable-Diffusion on MS-COC02014 prompts
Codebase and evaluation The text-to-image sampling on MS-COC02014 [26] prompts based on the pretrained latent-space diffusion model provided by Stable-Diffusion [43]. We use their official codebase of PyTorch implementation and their checkpoint sd-v1-4. We compute MSE on
randomly selected captions from the MS-COC02014 validation dataset, detailed in Section 4.1.
EMS computing We estimate the EMS at N 250 uniform timesteps Aju: A Ay....., AiN by drawing is trained on the LAION-5B dataset [46], there a gap between the images in the MS-COC02014 validation dataset and the images generated by Stable- Sifle-Diffuuson with certain guidance scale. Thus, we choose Qo to be the distribution of the latents generated by Diffussonn with corresponding guidance scale, using 200-step DPM- Solver++ [32]. We generate these latents with random captions and Gaussian noise different from those we use to compute MSE. The total time for EMS computing is ~1 Ih on 8 GPU cards of NVIDIA A40 for each guidance scale.
K = 1024 datapoints X ~ 9o- Since
Sampling Following previous works [32, 58], we start time 10-3, end time T 1 and adopt the uniform t timestep schedule. For DPM-Solver-yy we use the 2nd-order predictor with the pseudo 3rd-order corrector.
D.6
License
Table 3: The used datasets, codes and their licenses.
Name
URL
Citation
License
CIFARIO
bttpa://wxw.. www.:ca.ttoooooooooooooo toronto..o.oo...... edu/ "etu/"rri/// cifar .tml
[24]
L.UNN-Bedmur
[55]
ImageNet-256
httpa: w/w...inggg..........m
[9]
MS-COC02014
https: //toco/otatte.ee...cccccccc
[26]
CC BC 4.0
ScoreSDE
https: ggt.huu.com/yaaae.
[S]]
Apache-2.0
EDM
https: huop::///itthiee /tt.thgt.ccom//t////
[2]]
BY-NC-SA 4.0
Cuuded--Diffusion
https:// http//ggbbccooo//eeeeeeee
diffusion
[10]
MIT
Latent-Diffusion
https: https://github... ggithubb.com/Com/iiii/ii....
fiuui
[43]
MIT
Stuble-Diflusio
sttpa://gither
[43]
CreativeML Open RAIL-M
DPM-Solver++
hteps://github http::///gtabbboc////eeeegg///aa
[32]
MIT
UniPC
Attpe://gitae
com/ aolu/vl-aae zhazz Un1PC
[58]
We list the used datasets, codes and their licenses in Table 3.
E Runtime Comparison
As we have mentioned in Section 4, the runtime of DPM--olver-v3 is almost the same as other solvers (DDIM [48], DPM-Solver [31], DPM-Solver++ [32], UniPC [58], etc.) long they use the same NFE. This is because the main computation costs are the serial evaluations of the large neural network eø. and the other coefficients either analytically computed [48, 31, 32, 58], or precomputed (DPM-Solver-v3), thus having neglectable costs.
Table 4 shows the runtime of DPM-Solver-v3 and some other solvers on single NVIDIA A40 under different settings. We use torch. cuda. Event and torch. cuda. synchronize accurately compute the runtime. We evaluate the runtime 8 batches (dropping the first batch since it contains
27

---

Table 4: Runtime of different methods to generate single batch (second batch, =std) on a single NVIDIA A40, varying the number of function evaluations (NFE). We don't include the runtime of the decoding stage for latent-space DPMs.
NFE
Method
5
10
15
20
CIFARIO [24], ScoreSDE [51] (batch size 128)
DPM-Solver++ [32]
1.253(++0.001))
2.503(+0.0011)
3.754(#0.0022)
5.010(==0.0048) 5.080(=0.0049)
UniPC [58]
2.532(=000018) (.0000))))0
3.803(=0.0037)
DPM-Solver-v3
1.273(110.0005)
2.540(#0.0023)
3.826(=0.0039)
5.108(20.0055)
CIFARIO 124], EDM [21] (batch size 128)
DoM-Solver++ [32]
1.137(=0.0011) 2.278(=000015))
3,426(#0.0024)
4.569(=0.0031)
UniPC [58]
2.289(+0.0019) 1.142(10.0016)
3,441(#0.0035)
4.590(=0.0021)
DPM-Solver-v3
1.146(+0.0010)
2.293(#0.0015)
3,448(20.0018)
4.600(+0.002))
ASUN-Bedromd [55], Latent-Diffusion [43] (batch size = 32)
DPM-Solver++ [32]
1.302(+0:0009)
2.608(+000010)
3.921(=0.0023)
5.236(10.0045)
UniPC [58]
2.616(+:0.011) 1.302(=:0.0100)
1.305(+:0.0005)
3,934(=0.0033)
5.244(::00043))
DPM-Solver-v3
2.620(#0.0027)
3.982(=0.0028)
5.290(=0.0030)
ImageNer256 [9], Guided-Diffistee [10] (batch size 4)
DPM-Solver++ [32]
1.594(+000011)
3.194(+0.0018)
4.792(+0.0031)
6.391(+000045)
UniPC [58]
1.606(+0.0026)
3.205(+:00025))
4.814(20.0049
6.427(=000060)
DPM-Solver-v3
1.601(10.0059)
3.229(+0.0031)
4.807(#0.0068)
6.458((1.0255)
MS-COC02014 [26], Stable-Diffusion [43] (batch size = 4)
DPM-Solver++ [32]
1.732(+:.00122)
3.464(+0.0020)
6.974(20.0013) 5.229(=0.0027)
UniPC [58]
1.735(::0.0122)
3.484(+0.0064)
5.212(=0.0015)
6.988(1000035)
DPM-Solver-v3
1.731(+000008) 3.471(+0.0011)
5.211(=0.0030) 6.945(=0.0022))
extra initializations) and report the mean and std. We can see that the runtime is proportional to NFE and has a difference of about +1% for different solvers, which confirms our statement. Therefore, the speedup for the NFE is almost the actual speedup of the runtime.
Quantitative Results
Table 5:
Quantitative results on LSUN-Bedroom [55]. We report the FID, of the methods with different numbers of function evaluations (NFE), evaluated on 50k samples.
NFE
Method
Model
5
6
8
10
12
15
20
DPM-Solver++ [32]
18.59
8.50
4.19
3.63
3.43
3.29
3.16
UniPC [58]
Satent-Diffusion [43]
12.24
6.19
4.00
3.56
3.34
3.18
3.07
DPM-Solver-v3
7.54
4.79
3.53
3.16
3.06
3.05
3.05
We present the detailed quantitative results of Section 4.1 for different datasets in Table 1, Table 5, Table 6 and Table 7 respectively. They clearly verify that DPM-Solver-v3 achieves consistently better comparable performance under various settings, especially in 5~10 NFEs.
G
Ablations
In this section, we conduct some ablations to further evaluate and analyze the effectiveness of DPM-Solver-v3.
28

---

Table 6: Quantitative results on ImageNet-256 [9]. We report the FID., of the methods with different numbers of function evaluations (NFE), evaluated on 10k samples.
NFB
Method
Model
10
12
15
20
+ [32]
16.87
13.09
9.95
8.72
8.13
7.73
7.48
Siiled-Diffusin [10]
UniPC [58]
15.62
11.91
9.29
8.35
7.95
7.64
7.44
2.0)
DPM-Solver-v3
15.10
11.39
8.96
8.27
7.94
7.62
7.39
Table
Quantitative results on MS-COC02014 [26] prompts. We report the MSEJ of the methods with different numbers of function evaluations (NFE), evaluated on 10k samples.
NFE
Method
Model
5
DPM-Solver++ [32]
0.076
0.056
0.028
0.016
0.012
0.009
0.006
[43]
UniPC [58]
0.055
0.039
0.024
0.012
0.007
0.005
0.002
== 1.5)
DPM-Solver-v3
0.037
0.027
0.024
0.007
0.005
0.001
0.002
DPM-Solver++ [32]
0.60
0.65
0.50
0.46
0.42
0.38
0.30
Stable-Diffusion [43]
UniPC [58]
0.65
0.71
0.56
0.46
0.43
0.35
0.31
7.5)
DPM-Solver-v3
0.55
0.64
0.49
0.40
0.45
0.34
0.29
G.1 Varying the Number of Timesteps and Datapoints for the EMS
Table 8: Ablation of the number of timesteps N and datapoints K for the EMS, experimented with ScoreSDE [51] on CIFARIO [24]. We report the FIDJ with different numbers of function evaluations (NFE), evaluated 50k samples.
NFE
N
K
5
8
10
12
20
1200
512
18.84
7.90
4.49
3.74
3.88
3.52
3.12
1200
1024
15.52
7.55
4.17
3.56
3.37
3.03
2.78
120
4096
13.67
7.60
4.09
3.49
3.24
2.90
2.70
250
4096
13.28
7.56
4.00
3.45
3.22
2.92
2.70
1200
4096
12.76
7.40
3.94
3.40
3.24
2.91
2.71
First, we'd like to investigate how the number of timesteps and the number datapoints for computing the EMS affects the performance. We conduct experiments with the DPM ScoreSDE [51] on CIFAR10 [24], by decreasing N and K from our default choice N 1200, K 4096.
We list the FID results using the EMS of different N and K in Table 8. We can observe that the number of datapoints K is crucial the performance, while the number of timesteps N is less significant and affects mainly the performance in 5~100 NFEs. When NFE>10, we can decrease N as little as 50, which gives even better FIDs. Note that the time cost for computing the EMS is proportional to NK, so how to choose appropriate N and K for both efficiency and accuracy is worth studying.
G.2 First-Order Comparison
As stated in Appendix A, the first-order case of DPM-SSlver-v3 (DPM-Solver--3-I) is different from DDIM [48], which is the previous best first-order solver for DPMs. Note that 1oM--olver--31 applies no corrector, since any corrector has an order of at least 2.
In Table 9 and Figure 7, we compare DPM-Solver-w3-1 with DDIM both quantitatively and qual- itatively, using the DPM ScoreSDE [51] on CIFARTO [24]. The results verify our statement that DPM-Solver-v3-1 performs better than DDIM.
29

---

Table 9: Quantitative comparison of first-order solvers (DPM-Solver-v3-1 and DDIM [48]), ex- perimented with ScoreSDE [51] on CIFAR10 [24]. We report the FID with different numbers of function evaluations (NFE), evaluated on 50k samples.
NFE
Method
10
12
15
20
DDIM [48]
54.56
41.92
27.51
20.11
15.64
12.05
9.00
IPM--olver-v3-1
39.18
29.82
20.03
14.98
11.86
9.34
7.19
=
NFE
NFE
20.11
DPM-Solver-v3-1 (Ours)
HID 39.18
FID 14.98
FID 7.19
Figure 7: Random samples by first-order solvers (DPM-Solver-v3.) and DDIM [48]) of ScoreSDE [51] on CIFARIÓ dataset [24], using 5, 10 and 20 NFE.
G.3 Effects of Pseudo-Order Solver
We now demonstrate the effectiveness of the pseudo-order solvee, including the pseudo-order predictor and the pseudo-order corrector.
Pseudo-order predictor The pseudo-order predictor is only applied in one case (at 5 NFE on CIFARIO [24]) to achieve maximum performance improvement. In such cases, without the pseudo- order predictor, the FID results will degenerate from 12.76 to 15.91 for ScoreSDE [51]. and from 12.21 12.72 for EDM [21]. While they are still better than previous methods, the pseudo-order predictor is proven to further boost the performance at NPEs as small as 5.
Pseudo-order corrector We show the comparison between true and pseudo-order corrector in Table 10. We can observe consistent improvement when switching to the pseudo-order corrector. Thus, it suggests that if we n-th order predictor, we'd better combine with pseudo (n + 1)-th order corrector rather than (n + 1)-th order corrector.
G.4 Effects of Half-Corrector
We demonstrate the effects of half-corrector in Table using the popular Stable-Diffusiin
model [43]. We can observe that under the relatively large guidance scale of 7.5 which is necessary
for producing samples high quality, the corrector adopted by UniPC [58] has a negative effect on the convergence to the ground-truth samples, making UniPC even worse than DPM-Solver++ [32]. When we employ the half-corrector technique, the problem is partially alleviated. Still, it lags behind DPM- APM--oller-v3 v3, since we further incorporate the EMS.
30

---

Table 10: Effects of pseudo-order corrector under different settings. We report the FID4 with different numbers function evaluations (NFE).
NFE
Method
5
6
8
10
12
15
20
ISUN--edroom [55], Latent-Diffusion
[43]
4th-order corrector
8.83
5.28
3.65
3.27
3.17
3.14
3.13
>pseudo (default)
7.54
4.79
3.53
3.16
3.06
3.05
3.05
ImageNet-256 [9], Quided-Difffusien
[10] (s
2.0)
3rd-order corrector
15.87
11.91
9.27
8.37
7.97
7.62
7.47
>pseudo (default)
15.10
11.39
8.96
8.27
7.94
7.62
7.39
MS-COC02014 [26]. Stable-Diffusion
[43] (s 1.5)
3rd-order corrector
0.037
0.028
0.028
0.014
0.0078
0.0024
0.0011
pseudo (default)
0.037
0.027
0.024
0.0065
0.0048
0.0014
0.0022
Table 11:
Ablation of half--oorectoofuullccorretoo on MS-COC02014 [26] prompts with Stable- Diffusion model [43] and guidance scale 7.5. We report the MSEJ of the methods with different numbers of function evaluations (NFE), evaluated on 10k samples.
NFE
Method
Corrector Usage
6
10
12
15
20
DPM-Solver++ [32]
no corrector
0.60
0.65
0.50
0.46
0.42
0.38
0.30
full-corrector
0.65
0.71
0.56
0.46
0.43
0.35
0.31
UniPC [58]
+half-corrector
0.59
0.66
0.50
0.46
0.41
0.38
0.30
full-corrector
0.65
0.67
0.49
0.40
0.47
0.34
0.30
DPM-Solver-v3
hhaff-corrector
0.55
0.64
0.51
0.44
0,45
0.36
0.29
G.5
Singlestep vs. Multistep
Table 12: Quantitative comparison of single-step methods (S) multi-step methods (M), perimented with ScoreSDE [51] on CIFAR10 [24]. We report the FID,, with different numbers of function evaluations (NFE), evaluated on 50k samples.
NFE
Method
10
15
20
25
DPM-Solver (S) [31]
290.51
23.78
23.51
4.67
4.97
3.34
2.85
2.70
DPM-Solver (M) [32]
27.40
17.85
9.04
6.41
5.31
4.10
3.30
2.98
DPM-Solver++ (S) [32]
51.80
38.54
12.13
6.52
6.36
4.56
3.52
3.09
MPM--Slver++ (M) [32]
28.53
13.48
5.34
4.01
4.04
3.32
2.90
2.76
DPM-Solver-y3 (S)
21.83
16.81
7.93
5.76
5.17
3.99
3.22
2.96
DPM-Solver-v3 (M)
12.76
7.40
3.94
3.40
3.24
2.91
2.71
2.64
As we stated in Section 3.2.2, multistep methods perform better than singlestep methods. To study the relationship between parameterization and these solver types, we develop singlestep version of DPM-Solver-v3 in Algorithm Note that since (n + 1)-th order singlestep solver divides each step into n + substeps, the total number of timesteps is multiple of + 1. For flexibility, we follow the adaptive third-order strategy of DPM-Solver [31], which first takes third-order steps and then takes first-order second-orde steps in the end.
31

---

Algorithm 3 (n + 1)-th order singlestep solver
Require: noise prediction model co, noise schedule at, ot, coefficients 4x, SA, b2, cache Q1, Q2
Input: timesteps {t4 {tt)(++)))) (M+1)M iinitial value 20
Q1 cache
30
2: Q2 cache
EB (Zo, to) 3: for m = 1
do 4: to (n + (m 1)
5: for = 0 to do
fetch Q1
6:
, Dio-t t
selch
7:
Eig Eipts
Q2
12 o.dr 0 - -x,
- Ar - x, b, 1 = io...... 20....,20 + (Eq. (8))
8:
aA
9:
Digti+1 f (1
-10-1h..... ftiots Poo ftig Dio gio t1otit1)
if (i+ 1)m # (n 1)M then
11:
bitit1 t Eo( ligtit1)))
12:
Q1 t Ebptit1
noche
Q2 coche Eigtitl
13:
14:
end if
15:
end for end for
16:
Output: (n+1)M
Since UniPC [58] based on a multistep predictor-corredtee framework and has no singlestep version, we only compare with DPM-Solver [31] and DPM-Solver++ [32], which uses noise prediction and data prediction respectively. The results ScoreSDE model [51] on CIFARIO [24] are shown in Table 12, which demonstrate that different parameterizzaions have different relative performance under singlestep and multistep methods.
Specifically, for singlestep methods, DPM-Solver-v3 (S) outperforms DPM-Solverss (S) across NFEs, but DPM-Solver (S) even better than DPM-Solver-v's (S) when NFE2 10 (though when NFE<10, DPM-Solver (S) has worst performance), which suggests that noise prediction is best strategy in such scenarios; for multistep methods, we have 3PM--Slver-v33 (M) DPM-Solver++ (M) > DPM-Solver (M) across NFEs. Moreover, DPM-Solver (S) even outperforms DPM-Solver++ (M) when NFE>20,, and the properties of different parameterizations in singlestep methods are left for future study. Overall, DPM-Solver-v3 (M) achieves the best results among all these methods, and performs the most stably under different NFEs.
H FID/CLIP Score on Stable-Diffusion
Sample quality and text-image alignment performance on MS-COC02014 [26] prompts with model [43] and guidance scale 7.5. We report the FID, and CLIP score the methods with different numbers of function evaluations (NFE), evaluated on 10k samples.
Table 13:
NFE
Method
Metric
10
12
15
20
FID
18.87
17.44
16.40
15.93
15.78
15.84
15.72
DPM-Solver++ [32]
CLIP score
0.263
0.265
0.265
0.265
0.266
0.265
0.265
FID
18.77
17.32
16.20
16.15
16.09
16.06
15.94
UniPC [58]
CLIP score
0.262
0.263
0.265
0.265
0.265
0.265
0.265
FID
18.83
16.41
15.41
15.32
15.13
15.30
15.23
DPM-Solver-v3
CLIP score
0.260
0.262
0.264
0.265
0.265
0.265
0.265
32

---

In text-to- toiaage generation, since the sample quality is affected not only by discretization error of the sampling process, but also by estimation error of neural networks during training, low MSE (faster convergence) does not necessarily imply better sample quality. Therefore, we choose the MSCOC02014 [26] validation set as the reference, and additionally evaluate DPM-Solver-v3 on Stable-Diffusion model [43] by the standard metrics FID and CLIP score [41] which measure the sample quality and text-image alignment respectively. For DPM-Solver-v3, we use the full-corrector strategy when NFE<10, and no corrector when NO.
The results in Table 13 show that DPM-Solver-v3 achieves consistently better FID and similar CLIP scores. Notably, we achieve an FID of 15,4 in 8 NFE, close to the reported FID of Stable-Diffussi v1.4.
Still, we claim that FID is not a proper metric for evaluating the convergence of latent-space diffusion models. As stated in DPM-Solver++ and Section 4.1,
see that the FIDs quickly achieve 15.0~16.0 within 10 steps, even the latent code does not converge, because of the strong image decoder. Instead, MSE in the latent space is direct way measure the convergence. By comparing the MSE, our sampler does converge faster to the ground-truth samples of Stable Diffusion itself.
I More Theoretical Analyses
L.1 Expressive Power of Our Generalized Parameterization
Though the introduced coefficients L,, 8x, by seem limited guarantee the optimality of the pa- rameterization formulation itself, we claim that the generalized parameterization ge in Eq. (8) can actually cover wide range of parameterization families in the form of the (I A) = (AA)ca (2xx, A)+ ((A)a> + 7(A). Considering the paramerization on [>.. Xr). by rearranging the terms, Eq. (8) can be written
A))
lx A) = S2 s, σx co( - e Sh. *, dr
xx ax
e JN, s, b,dr A
(68)
We compare the coefficients before Eg and a in ge and the to figure how l, 8A, by corre- sponds to ((A), 7A)). In fact, we can not directly let go equal the since when X = A. we have Sr (-) 0, and the coefficient before En in go is fixed. Still, the can be equalized go by linear transformation, which only depends on A and does not affect our analyses of the discretization егror and solver.
Specifically, assuming the = wr.90 EA., by corresponding the coefficients we have
Wx, (As)
a(A)= = wx,e" 12 dr 0s as 3(1) = -wse -hh, o, la
Ex, (A,))
l = - cx
B(A)
(69)
7(A) = -w, SA, - Sr o,dr d + Ex.
= - (2)
by - c A (A)
((
Therefore, long (xA) and (A) have first-order derivatives, our proposed parameteri- zation ge holds the same expressive power as æhe, while at the same time enabling the neat optimality criteria of lx, sx, by in Eq. (5) and Eq. (11).
1.2 Justification of Why Minimizing First-order Discretization Error Can Help Higher-order Solver
The EMS 87 bj in Eq. (11) are designed to minimize the first-order discretization error Eq. (10). However, the high-order solver is actually more frequently adopted in practice (specifically, third- order in unconditional sampling, second-order in conditional sampling), since it incurs lower sampling errors by taking higher-order Taylor expansions to approximate the predictor go.
33

---

In the following, we show that the EMS can also help high-order solver. By Eq. (52) in Appendix B.4, the (n + 1)-th order local error can be expressed as
2t 2, 10-2(=0,A1 a,A()s))))) A
Ah Ex, Ex, (90(22)1)) - go Xo)) dd
Σ (unaoaul Au) go ( Anl) S/m
n
Ex (2) (A
Aa)*
- -a,A(ds, )x)
da
k!
(70)
By Newton-Leibniz theorem, it equivalent to
Ar ④t - x2 = a0)((22;))) (/m (N) (1, 10000 (x t)ar) (trr) ad
(71)
n
( 1 (dlne l.m Ex, (2) (A) Jy, I0
12)x dd
-aiA(Ao1
Σ (R )nl
k!
We assume that the estimated EMS are bounded (in the order of (1)), Assumption B.2 in Ap- pendix B.I), which is empirically confirmed as in Section 4.2. By the definition of go in Eq. (8), we have gø (I,, = e Jz, d. (x,, r) - + fo +) - by). Therefore, Eq. (11) controls
Igg 2 lle and further controls || x since other terms are only dependent on the EMS and are bounded.
1.3 The Extra Error EMS Estimation and Integral Estimation
Analysis of EMS estimation error practice, the EMS in Eq. (5) and Eq. (11) are estimated on finite datapoints by the explicit expressions in Eq. (58) and Eq. (63), which may differ from the true 13, x, bz. Theoretically, on one hand, the order and convergence theorems in Section 3.2 are irreevant to the EMS estimation error: The ODE solution in Eq. (9) is correct whatever Lx, 8x, bx are, and we only need the assumption that these coefficients are bounded (Assumption B.2 in Appendix B) to prove the local and global order; on the other hand, the first-order discretization error in Eq. (10) is vulnerable to the EMS estimation error, which relates to the performance at few steps. Empirically, to enable fast sampling, EMS (see ablations in Table 8), and we find that our method is robust to the EMS estimation error given only 1024 datapoints in most cases (see EMS computing configs in Appendix D).
need to ensure the number of datapoints for estimating
Analysis of integral estimation error Another source of error is the process of integral estima- tion (f²,e Ex, (A)dX and 1a Ex, (A) (1) (- dA in Eq. (14)) by trapezoidal rule. We can analyze the estimation error by the error bound formula of trapezoidal rule: suppose we use uniform discretization (a, with interval h to estimate f² f(z)dx then the error E satisfies
(b
a)h²
E
max (x))
(72)
12
34

---

Under Assumption B.2, the EMS and their first-order derivative are bounded. Denote fi(A) Ex, Bx) f2(2) - Ex, (2) (a-A.)"
then
kd
fi()) eR (l++add) // T- --J bndr b,dr
he
e Ii o b,dr + bxeli ldr S( e Si s b,dr
f((A) - (4 + sx)ees, (r+a+)dr
e-Ii,
A
e Sr ondr b,dr + (b2 + 8x) Rh fr E e-Sin or by
(((A)) - (1'+ + 85)esB (o+
by b,dr b`y b'e lodr ( e- or bydr
+ (l2 sxubretee Lidr
ardr e Sh e-
balxedi do fx Sh bydr + bze (1, -a)) fr e eBi a. a,
(73)
f2(A) = Bi (dr (A Xa)a
k!
15(A) = (L2 + 8x)e (dr (A Aa)
+ e- (l,ta)) ()- 2a)a-1
k!
(I: - 1)!
(2 = (1's + 25) (1- ta-)d- (A - >x,
+ (1x + 0^) eels (++++) )dr (A Aa)
(74)
k!
k!
(1a sx)e (r+er)dr (A xx)
+ (l2 + sx)e (l++a)ddr (A - \₂) 2a)
k 1!
(k 1)!
ed2, (1, + )dr (A- (A-A))*-2 (A: 2)!
Since L2, sx, bx, 12, S b2 are all O(1), we conclude that (A) and the errors of Sh Ex, (x) (A)dA, 1N) (A) (2- dd are 0(1312), respec- tively, where ho is the stepsize of EMS discretization (ho = d= n corresponds to 120~1200 timesteps for our EMS computing), and h A As. Therefore, the extra error of integral estimation is under high order and ignorable.
J. More Discussions
J.1 Extra Computational and Memory Costs
The extra memory cost of DPM-Solver-v3 is rather small. The extra coefficients lx 8x, by are discretized and computed at N timesteps, each with dimension D same the diffused data. The extra memory cost is O(ND), including the precomputed terms in Appendix 1.1.2, and is rather small compared to the pretrained model (e.g. only 125M in total Stable-Diffusion, compared to 44G of the model itself).
The pre-computation time for estimating EMS is rather short. The EMS introduced by our method can be effectively estimated on around 1k datapoints within hours (Appendix D), which is rather short compared to the long training/distillation time of other methods. Moreover, the integrals of these extra coefficients are just some vector constants that can be pre-computed within seconds, as shown in Appendix 1.2. The precomputing is done only once before sampling.
The extra computational overhead of P during sampling is negligible. Once we obtain the estimated EMS and their integrals at discrete timesteps, they can be regarded as constants. Thus, during the subsequent sampling process, the computational overhead the same as previous training-fie free methods (such as DPM-Solver++) with negligible differences (Appendix E).
35

---

By a.tr
10
Rydt-
-15
Figure 8: Visualization of the EMS 8x, by and their integrals w.r.t X, estimated ScoreSDE [51] on CIFARIO [24]. 8x, by are rather fluctuating, but their integrals are smooth enough to ensure the stability of DPM-Solver-v3.
J.2 Flexibility
The pre-computed EMS can be applied for any time schedule during sampling without re-computing EMS. Besides, we compute EMS on unconditional models and it can be used for wide range of guidance scales (such as cfg=7.5 in Stable Diffusion). In short, EMS is flexible and easy to adopt in downstream applications.
Time schedule: The choice of N, K in EMS is disentangled with the timestep scheduler in sampling. Once we have estimated the EMS at N (e.g., 1200) timesteps, they can be flexibly adapted any schedule (uniform A/uniform t...) during sampling, by corresponding the actual timesteps during sampling to the N bins. For different time schedule, we only need re-precompute E(k) Appendix C.1.2, and the time cost within seconds.
Guided sampling: We compute the EMS on the unconditional model for all guided cases. Empirically, the EMS computed on the model without guidance (unconditional part) per- forms more stably than those computed on the model with guidance, and can accelerate the sampling procedure in a wide range of guidance scales (including the common guidance scales used in pretrained models). We think that the unconditional model contains some common information (image priors) for all the conditions, such as color, sketch, and other image patterns. Extracting them helps correct low saturation level and lack of details. In contrast, the conditional model is dependent on the condition and has a large variance.
common biases such shallow color,
Remark J.1. EMS computed models without guidance cannot work for extremely large guidance scales (c.g.. cfg scale 15 for Stable Diffusion), since in this case, the condition has large impact on the denoising process. Note that at these extremely large scales, the sampling quality very low (compared to cfg scale 7.5) and they are rarely used in practice. Therefore, our proposed EMS without guidance suitable enough for the common applications with the common guidance.
J.3 Stability
As shown Section 4.2, the estimated EMS 8 bx appear much fluctuating, especially for ScoreSDE on CIFARIO. We would like to clarify that the unstable x by is not issue, and our sampler is stable:
The fluctuation of 8x, by on ScoreSDE is intrinsic and not due to the estimation error. As we increase the number of samples to decrease the estimation error, the fluctuation is not reduced. We attribute it to the periodicity of trigonometric functions the positional timestep embedding stated in Section 4.2.
36

---

Moreover, we only need to consider the integrals of 8x, by in the ODE solution Eq. (9). As shown in Figure 8, the integrals of x. b2 are rather smooth, which ensures the stability of our method. Therefore, there is no need for extra smoothing.
J.4 Practical Value
When NFE is around 20, our improvement of sample quality is small because all different fast samplers based on diffusion ODEs almost converge. Therefore, what matters is that our method has faster convergence speed to good sample quality. As we observed. the less diverse the domain is, the more evidently the speed-up takes effect. For example, on up 40% faster convergence.
our method can achieve
To sum up, the practical value of DPM-Solver-v3 embodies the following aspects:
15~300% speed-up can save lots of costs for online text-to-image applications. Our speed- ups applicable large text-to-image diffusion models, which today's AIGC community. As the recent models become larger, a single NEE requires more computational resources, and 15~30% speed-up save lot of the companies' expenses for commercial usage.
important part of
2. Improvement 5~~0 NFES benefits image previewing. Since the samples are controlled by the random seed, coarse samples with 5--10 NFEs can be used to preview thousands of samples with low costs and give guidance on choosing the random seed, which can then used to generate fine samples with best-quality sampling strategies. This is especially useful for text-to-ine tottgeiggeeeeeeeee generation. Since our method achieves better quality and converges faster to the ground-truth sample, it can provide better guidance when used for preview.
K Additional Samples
We provide more visual samples Figure 9, Figure 10, Figure 11 and Table 14 to demonstrate the qualitative effectiveness of DPM-Solver-v3. outperforms previous state-of-the-ard solvers. Our method can generate images that have reduced bias (less "shallow"), higher saturation level and more visual details, mentioned Section 4.3.
be seen that the visual quality of DPM-Solver-v3
37

---

NFE 5
NFF 10
DPM-Solver++ [32]
FID 28.53
4.01
UniPC [58]
FID 23.71
FID 3,93
DPM-Solver-v3 (Ours)
FID 12.76
FID 3.40
Figure 9: Random samples of ScoreSDE [51] on CIFARIO dataset [24] with only 5 and 10 NFE.
38

---

Hcun's 2nd [21]
FID 320.80
FID 16.57
D [32]
FID 24.54
FID 2.91
UniPC [58]
FID 23.52
2.85
DPM-Solver-v3 (Ours)
FID 12.21
FID 2.51
Figure 10: Random samples of EDM [21] on CIFAR10 dataset [24] with only 5 and 10 NFE.
39

---

DPM-Salver++ [32]
(FID 11.02)
THM
UniPC [58] (FID 10.19)
DPM-Solver-v3 (Ours) (FID 9.70)
Figure 11: Random samples of Guided-Diffusion [10] ImageNet-256 dataset [9] with classifier guidance scale 2.0, using only 7 NFE. We manually remove the potentially disturbing images such as those containing snakes or insects.
40

---

Table 14: Additional samples of Stable-Diffusion [43] with a classifier-free guidance scale 7.5, using only 5 NFE and selected text prompts. Some displayed prompts truncated.
HPM-Sollert+
UniPC [58] (MSE 0.65)
DPM-Solver-v3
Text Prompis
[32] (MSR
(Ours) (MSE 0.55)
"pixar movie still portrait photo of modison beer, jessicu albu, woman, as hero catgirl cyborg woman pixar, greg nutkowski. wlop, rossdraws, artgerm, weta, marvel, rove girl, leeloo, urreal engine, glossy skin, pearlescent, wet, bright morning, anime, sci-fi, maxim magazine cover"
"oll painting with heavy impasto of pirate ship and its captain, cosmic horror painting, elegant intricate artstation concept art by craig mullins detailed"
"environment living room interior, century modern, indoor garden fountain, retm, vintage, designer furniture made wood and plastic, concrete table, wood walls, indoor potted tree, large window, outdoor forest landscape, bematiful sunset, cinematic, concept art, sunstainable architecture, octone render, stopiu, ethereal, cinematic light"
"the living rom of coty wooden house with fireplace, at night, interior design, concept art, wallpoper, warm, digital art. art by james gurney and larry elmore."
"Full page concept design how eraft life Poison, intricate details, in- fographic of alchemical, diogram of how make potions, captions, directions, ingredients, drawing, magic, waxia"
Fantasy octone render, 16h, Sk, envirmment, Wood pavilion architecture, warm led lighting, dusk, Land- scape, snow, arctic, with aqua water, silver Guggenheim museum spire, with rays of sunshine, white fabric landscape, tall building, zaha hadid and Santiago calatrava, smooth landscape, cracked ice. igloo, warm lighting, autora berialis, 3d egi, high definition, natural lighting, realis- tic, hyper realism"
back causties,
"tree house in the forest, atmospheric. hyper realistic, epic composition. cinematic, landscape vista photography by Carr Clifion & Galen Rowell, 16K resolution, Landscape veduta photo by Dustin Lefevre & idraw, de- sailed landscape painting by Ivan Shishkin, DeviantAr, Flickr, rendered
in Enscape, Miyazaki, Nausicaa Ghibli, Breath of The Wild, 4k detailed post processing, artstation, unreal engine"
"A mail through the unknown, atmospheric, hyper realistic, 8k, epic com- position, cinematic, octime render, aristation landscape vista photogra- phy by Clifton & Galen Rowell, 16K resolution, Landscape veduta photo by Dustin Lefevre tdraw, 8k resolation, detailed landscape paint- by Ivan Shishdin, DeviontArt, Flickr, rendered in Enscape, Miyazaki, Nausicaa Glubli, Breath of The Wild, 4k detailed post processing, artsta- tion, rendering octane,
engine"
postopocalsptic city turned to fracial glass, clone render, 8 exploration, cinematic, trending on artstation, by beeple, mallisic,, 3 mm camera, un- real engine, hyper detailed. photo-read Mhoto-realistic maximum detal, volumetric light, moody cinematic epic concept realistic matte painting, hyper photorealistic, concept art, cinematic epic. octane render, Sk. corona render, novie concept art, octane render, k, corona render, trending cinematic composition, ultra-detailed, reper-realistic, volumetric lighting"
"OORDDS": zoological fantasy ecosystem infographics, magazine lay- oud with typography, annotations, in the style of Elena Masci, Studio Ghibli, Caspar David Friedrich, Daniel Merriam, Doug Chiang, Ivan Aivazovsky, Herbert Bauer, Edward Tufte, David McCandless"
41


Reviews:
[{"review_id": "JbWSELnJQu", "reviewer": "Reviewer_dqUr", "summary": "This paper proposed an improved/generalized version of DPM-solver. It include three more coefficients based on the semi-linear ODE solution proposed by DPM-solver, accounting for minimizing the 'linearity' for the nonlinear component of the ODE and minimizing the first-order discretization error respectively. Such optimal coefficients can be approximately calculated by Monte Carlo samples from the learned model. Together with multistep formulation, predictor-corrector framework by UniPC, and some practical techniques, the paper empirically demonstrated the proposed sample can outperform existing ones especially when NFE is very small. ", "questions": "- Why EMS without guidance can be directly applied to model with guidance? Does it still hold if the guidance weight is really high?\n- Why for conditional setting, the best setting is 2nd-solver instead of higher-order counterparts? ", "limitations": "I'd expect to see more discussion on limitations due to the intractable integrals. ", "rating": 6, "confidence": 5, "soundness": 2, "presentation": 3, "contribution": 3, "strengths": "- The paper is well-written and easy to follow. \n- The paper built on DPM-solve formulation and generalized it with additional coefficients that are optimized based on trained model. The added coefficients are intuitive. Theoretical analysis has been carried out to ensure local/global convergence. \n- Two practical techniques has been proposed and seems to be useful in general. \n- Extensive experiments have been performed to justify the effectiveness.", "weaknesses": "- I'm not convinced by optimizing $s_\\tau$ with Eq. (11):\n  - In equation 10, the first coefficient inside the big integration should be $e^{\\int_{\\lambda_s}^\\lambda (l_\\tau + s_\\tau) d\\tau}$ instead of  $e^{\\int_{\\lambda_s}^\\lambda l_\\tau d\\tau}$. Therefore, the first-order discretization error not only depends on $f_\\theta^{(1)} - s_\\lambda f_\\theta - b_\\lambda$. There could be a chance that Eq. (11) is minimized but this term $e^{\\int_{\\lambda_s}^\\lambda (l_\\tau + s_\\tau) d\\tau}$ is amplified. \n  - Main results shown in this paper are not really by first-order solver but by higher-order solvers. A justification of why minimizing first-order discretization error can help higher-order solvers is needed. \n\n- A downside of the introduced coefficients is that the integrals in Eq. (14) become intractable, which has to be approximated by trapezoidal rule. It could introduce extra approximation error and add computational cost. I'd like to see a detailed analysis on the extra error and computational cost. \n\n-  It is understandable and acceptable that this work is heavily built on DPM-solver, but certain sentences in this paper is duplicated from DPM solver. E.g., line 94 our first key insight is..., line 105 our second key insight is ..., line 135-136: the proposed solvers and analysis ....The paper should at least rephrase these sentences. \n\n- The paper can do a better job on clarifying which is the unique contribution of this paper while which part has been proposed by previous works. For example, most techniques and theoretical conclusions in section 3.2 have already been developed in DPM-solver / DPM-solver++ / UniPC, which should be further admitted and clarified. ", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your positive comments. Below we provide detailed responses for the questions and clarify misunderstandings. We kindly request that you consider raising the score accordingly if you are satisfied.\n\n*W1: Equation 10 is wrong, the first-order discretization error may be amplified.*\n\n**A**: **We respectfully disagree. This is a misunderstanding**. Specifically, by the definition of $g_\\theta$, we have\n\n$$\ng_\\theta^{(1)}(x_\\lambda,\\lambda) = e^{-\\int_{\\lambda_s}^{\\lambda}s_\\tau d\\tau}\\left(\n        f_\\theta^{(1)}(x_\\lambda,\\lambda)-s_\\lambda f_\\theta(x_\\lambda,\\lambda)\n        - b_\\lambda\n    \\right)\n$$\n\nBy using the Taylor expansion $g_\\theta(x_{\\lambda_s},\\lambda_s) = g_\\theta(x_{\\lambda},\\lambda ) + (\\lambda_s-\\lambda)g_\\theta^{(1)}(x_{\\lambda},\\lambda) + O((\\lambda-\\lambda_s)^2)$ stated in the paper (**expand it at each $\\lambda$**), we can see that the term $e^{\\int_{\\lambda_s}^{\\lambda}s_\\tau d\\tau}$ in Eq. (9) is cancelled by the term $e^{-\\int_{\\lambda_s}^{\\lambda}s_\\tau d\\tau}$ in $g_\\theta^{(1)}(x_\\lambda,\\lambda)$, and we obtain exactly Eq. (10).\n\n*W2: A justification of why minimizing first-order discretization error can help higher-order solvers is needed*.\n\n**A**: Sure. Below we provide the corresponding theoretical justification about how minimizing first-order discretization error helps high-order cases.\n\nWe assume that the EMS are bounded (Assumption B.2 in Appendix B.1.1), which is empirically confirmed as in Section 4.2. Here we take the 2nd-order case as an example. By Eq. (52) in Appendix B.4, the 2nd-order local error is\n\n$$\n\\hat x_t-x_t=\\alpha_tA(\\lambda_s,\\lambda_t)\\int_{\\lambda_s}^{\\lambda_t}E_{\\lambda_s}(\\lambda)(g_\\theta(x_\\lambda,\\lambda)-g_\\theta(x_{\\lambda_s},\\lambda_s))d\\lambda-\\frac{\\alpha_tA(\\lambda_s,\\lambda_t)}{\\lambda_{i_1}-\\lambda_s}(g_\\theta(x_{\\lambda_{i_1}},\\lambda_{i_1})-g_\\theta(x_{\\lambda_s},\\lambda_s))\\int_{\\lambda_s}^{\\lambda_t}E_{\\lambda_s}(\\lambda)(\\lambda-\\lambda_s)d\\lambda\n$$\n\nwhich is equivalent to\n\n$$\n\\alpha_tA(\\lambda_s,\\lambda_t)\\int_{\\lambda_s}^{\\lambda_t}e^{\\int_{\\lambda_s}^{\\lambda}(l_{\\tau}+s_{\\tau})d\\tau}\\left(\\int_{\\lambda_s}^{\\lambda}g^{(1)}\\_\\theta(x_{\\tau},\\tau)d\\tau\\right)d\\lambda-\\frac{\\alpha_tA(\\lambda_s,\\lambda_t)}{\\lambda_{i_1}-\\lambda_s}\\int_{\\lambda_s}^{\\lambda_t}E_{\\lambda_s}(\\lambda)(\\lambda-\\lambda_s)d\\lambda\\int_{\\lambda_s}^{\\lambda_{i_1}}g^{(1)}\\_\\theta(x_{\\lambda},\\lambda)d\\lambda\n$$\n\nWhich is controlled by $\\|g_\\theta^{(1)}\\|$. Other terms are only dependent on the EMS and are bounded. For higher orders, the derivation is similar, since they all involve the form $g_{\\theta,t_2}-g_{\\theta,t_1}$, which is the integral of $g_\\theta^{(1)}$ . We will put the detailed analysis in the revised paper.\n\n*W3: A detailed analysis on the extra error and computational cost of the integrals w.r.t. the introduced coefficients.*\n\n**A**: Below we provide a detailed analysis correspondingly.\n- **The computational cost is ignorable**: See *Common response, Q1*.\n- **The extra error of trapezoidal rule is under higher order and ignorable**: We assume that the EMS and their derivatives are bounded, then by the error bound formula of trapezoidal rule,\n$$\n|E|\\leq \\frac{(b-a)^3}{12n^2}\\max |f''(x)|\n$$\nWe can easily conclude that the errors of $\\int_{\\lambda_s}^{\\lambda_t}E_{\\lambda_s}(\\lambda)B_{\\lambda_s}(\\lambda)d\\lambda$ and $\\int_{\\lambda_s}^{\\lambda_t}E_{\\lambda_s}(\\lambda)\\frac{(\\lambda-\\lambda_s)^k}{k!}d\\lambda$ in Eq. (14) are $O(h_0^2h)$ and $O(h_0^2h^{k-1})$ respectively, where $h_0$ is the stepsize of EMS discretization, and $h=\\lambda_t-\\lambda_s$. We'll put detailed analysis in the revised paper.\n\n*W4: Certain sentences is duplicated from DPM solver.*\n\n**A**: Thank you for pointing it out. We take these expressions since they concisely emphasize our motivation, novelty and how we get inspiration from previous works, and did not notice the duplicated expressions. We will rephrase them in the revised paper.\n\n*W5: Clarifying the unique contribution in section 3.2.*\n\n**A**: Thank you for the valuable suggestion. We coarsely admit that by saying *\"...highly motivated by ... in the field of diffusion models\"* at the beginning of Section 3.2, and the specific technical details are different. We will further clarify detailedly at each part of Section 3.2 to seperate the unique contributions of ours from the previous works in the revised version.\n\n*Q1: Why EMS without guidance can be directly applied to model with guidance? Does it still hold if the guidance weight is really high?*\n\n**A**: It is because the unconditional model contains some common information, please refer to *common response, Q2* for details. \n\nHowever, it no longer holds for extremely large guidance scale (e.g., cfg scale 15 for Stable Diffusion), since in this case the condition has a large impact on the denoising process. Note that at these extreme large scales, the sampling quality is very low (compared to cfg scale 7.5) and they are rarely used in practice. Therefore, our proposed EMS without guidance is suitable enough for the common applications with the common guidance.\n\n*Q2: Why for conditional setting, the best setting is 2nd-solver instead of higher-order counterparts?*\n\n**A**: The high-order solvers are more unstable compared to low-order, especially for large guidance scales[1]. However, a slightly large guidance scales are often preferred in practice to improve the condition-sample alignment, as shown in Imagen, Stable Diffusion. Therefore, we and previous works (e.g., DPM-Solver++, UniPC) all find that 3rd-solver will produce artifacts in conditional setting, and 2nd-solver is better.\n\n[1] Wizadwongsa, Suttisak, et al. \"Diffusion Sampling with Momentum for Mitigating Divergence Artifacts.\".\n\nOnce again, thank you for your constructive feedback and for considering our paper for acceptance. We'll revise our paper according to your suggestions.\n"}, {"review_id": "L37f9HVqFY", "reviewer": "Reviewer_yHtS", "summary": "This work proposes DPM-Solver-v3, an ODE solver for diffusion model inference. The method builds upon existing formulations involving exponential integrators, where the linear term of the ODE is canceled and only the noise predictor needs to be approximated. A new model parameterization is derived which involves an optimal set of coefficients, empirical model statistics (EMS). The authors describe methods for computing the EMS in practice, as well as a pseudo high-order method in the few NFE regime and a half-corrector technique when the guidance scale in guided sampling is large.", "questions": "- FID degenerating from 12.76 to 15.91 for ScoreSDE on CIFAR10 (line 822, G.3) in whether or not to use pseudo-order predictor/corrector seems quite substantial. In which cases is the pseudo-order predictor necessary? It might be compelling to produce a table similar to Table 10 (pseudo-order corrector experiments) examining the influence of the pseudo-order predictor.\n- How does choice of N, K for the EMS change as the timestep scheduler changes?", "limitations": "The authors adequately address the limitations of the work.", "rating": 6, "confidence": 2, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "- Although the method builds upon existing methods involving exponential integrators in the diffusion literature, the model parameterization involving EMS is novel.\n- Experiments and ablations are thorough and clearly done. Results are compelling and the proposed method consistently outperforms existing ODE solvers across datasets/models, especially in few NFE regime.\n- Theoretical results for local/global error are clearly presented and assumptions seem reasonable.", "weaknesses": "- In the experiments, the authors mention that \"we tune the strategies of whether to use pseudo-order predictor/corrector at each NFE on CIFAR10\" (lines 254-255). This seems to be an important caveat and should be reiterated in the appendix where the detailed empirical results are presented (Table 4).\n- Similarly, for Stable Diffusion, s=7.5, the authors mention that the best results among no/half/full corrector are reported. This should be reiterated in the appendix in Table 7.\n- In general, I wonder if the presentation of the paper would be more compelling if some of the quantitative results could be moved to the main body from the appendix.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "\nThank you for your detailed review, especially for reading the proofs and ablations in the Appendix. Below we provide detailed responses for the questions and clarify misunderstandings. We kindly request that you consider raising the score accordingly if you are satisfied.\n\n*W1 & W2: Some tuning strategies should be reiterated in the Appendix (Table 4/Table 7).*\n\n**A**: Thank you for carefully reading the Appendix! We fully agree that although they have been mentioned in the main text, they should be reiterated where the quantitative table results are presented. We will revise relative paragraphs in the final version.\n\n*W3: Some of the quantitative results could be moved to the main body from the appendix.*\n\n**A**: Thank you for the valuable suggestion. Due to the page limitation of 9, we will move some of them to the main body in the final version.\n\n*Q1: In which cases is the pseudo-order predictor necessary? It might be compelling to produce a table similar to Table 10 (pseudo-order corrector experiments) examining the influence of the pseudo-order predictor.*\n\n**A**: The pseudo-order predictor only provided improvement in our mentioned case (5 NFE on CIFAR10). It can be seen as a specialized trick at 5 NFE on CIFAR10.\n\n*Q2: How does choice of N, K for the EMS change as the timestep scheduler changes?*\n\n**A**: The choice of $N,K$ is **disentangled** with the timestep scheduler in sampling: once we have estimated the EMS at $N$ (e.g., 1200) timesteps, they can be **flexibly adapted** to any schedule (uniform logSNR/uniform t...) in sampling, by corresponding the actual timesteps during sampling to the $N$ bins. The choice of $N,K$ and timestep scheduler are dependent on the pretrained model. See more discussions in *common response, Q2*.\n\nOnce again, thank you for your constructive feedback and for considering our paper for acceptance. We'll revise our paper according to your suggestions.\n"}, {"review_id": "6vFRplo9Lz", "reviewer": "Reviewer_oHv2", "summary": "This work proposes a method to speed up the DPM solver by using empirical model statistics (pre-trained model) as well as employing multistep methods with a predictor-corrector framework for improving sample quality. Promising results were reported.", "questions": "Please address the comments in the Weaknesses section.\n\n- Will the speed-up of 15-30% have a practical value? Please discuss.", "limitations": "Yes", "rating": 4, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 2, "strengths": "- This work addresses an important problem of fast sampling in diffusion models without retraining.\n- The proposed parametric model is simple and seems neat with optimality criteria.", "weaknesses": "- The improvement of the proposed method over prior works such as DPM-solver++ seems incremental. It is unclear if the proposed empirical model statistics is good enough. While this work discusses an optimality under this parametric model, it is unclear if it is truly optimal in terms of modeling itself. This work may also discuss the prior work such as [H Zheng et al., Fast Sampling of Diffusion Models via Operator Learning, ICML 2023] and its arXiv version, which was able to achieve incredible results with only a single model evaluation.\n- The experiment results seem to need more improvements and should ensure fairness. For example, Figures 3, 4 indicate that the proposed method is faster than other methods in early iterations, but in order to ensure reasonable FID, one may need a certain number of NFE and in those cases, the differences of FIDs look negligible. Moreover, the comparisons in Fig 6 seems unfair - will other methods achieve similar image quality and FID scores with 15-30% more evaluations? Incorporating the factor of computation should be incorporated for fair comparisons.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for recognizing our work's importance and theoretical neatness. However, we feel that there are some fundamental misunderstandings. We will address your questions and concerns below. We kindly request that you consider raising the score accordingly if you are satisfied.\n\n*W1: The improvement seems incremental. Discuss with [H Zheng et al. 2023].*\n\n**A**: First of all, **we respectfully disagree with the comments on incremental results**. We appreciate you for informing us of the nice work of [H Zheng et al. 2023] and will add proper citations in related works, but we have to point out that [H Zheng et al. 2023] is a distillation-based method, different from what our work focuses on. We would like to make the following clarifications:\n1. We achieve the **most SOTA sample quality among all the training-free sampling methods** on a fair benchmark, under various settings, against the most recent and strong baselines.\n2. **Our method is training-free samplers, while [H Zheng et al. 2023] needs heavy training of an additional generative model.** In general, although distillation-based methods can achieve good performance even in one step, it has obvious flaws (as discussed in Appendix A.1) such as onerous extra training and lose of information. Moreover, **their application scope is usually limited to unconditional cases, and it is much costly to distill a text-to-image model**. Since our training-free approach **can be easily applied to large text-to-image models** to boost the sampling with negligible extra cost, while the distillation-based methods such as [H Zheng et al. 2023] needs much more costs to do so, we argue that our work has high practical value.\n\n*W2: It is unclear if the proposed empirical model statistics is good enough. It is unclear if the parameterization is truly optimal in terms of modeling itself.*\n\n**A**: The proposed empirical model statistics can greatly speed up the convergence and improve the sample quality, especially within 5~10 NFEs. **It has already achieves the SOTA performance among all training-free fast samplers**.\n\nMoreover, we did not claim the optimality of the modeling parameterization. We **first** systematically study the parameterization of exponential integrators for diffusion models and **first** explain the superiority of previous data-pred over noise-pred (Appendix A.2). Also, this parameterization is elegantly designed by Rosenbrock-type methods and elucidating the design space of previous parameterizations. We understand that a better form of parameterization may be derived, but our introduced parameterization is already novel, insightful and efficient.\n\n*W3: Experiment results does not have much improvement and unfair.*\n\n**A**: **We respectfully disagree with the comments on unfairness and poor results**. We would like to make the following clarifications:\n\n1. **The extra computation is negligible**. Please refer to *common response, Q1*.\n2. **Our evaluation is fair**. We feel that the expression \"early iterations\" implies a critical misunderstanding: when users generate images with a pretrained diffusion model, they predetermine a certain NFE. Therefore, comparing FID at the same NFE, or comparing NFE for the same FID is completely fair, and **are adopted by all the previous fast diffusion sampling methods**.\n3. **Our method converges much faster**. Our method converges 15-30% faster than previous samplers, thus saving 15-30% costs. **On a fair benchmark, under various settings, we achieve the best convergence speed among all the training-free samplers (Appendix F/H)**. A notable example is that, on LSUN-Bedroom, we reach the FID of 3.06 with 12 NFE, while the previous best method requires 20 NFE, which means our computation cost is approximately 60%.\n\n*Q1: Will the speed-up of 15-30% have a practical value?*\n\n**A**: Yes, the practical value embodies the following aspects:\n\n1. **15~30% speed-up can save lots of costs for online text-to-image applications**. As we have clarified, our speed-ups are applicable to large text-to-image diffusion models, which are an important part of today's AIGC community. As the recent models become larger, a single NFE requires more computational resources, and 15-30% speed-up can save a lot of the companies' expenses for commercial usage.\n2. **Improvement in 5~10 NFEs benefits image previewing**. Since the samples are controlled by the random seed, coarse samples with 5~10 NFEs can be used to ***preview*** thousands of samples with low costs and give guidance on choosing the random seed, which can be then used to generate fine samples with best-quality sampling strategies. This is especially useful for text-to-image generation. Since our method achieves better quality and converges better to the ground-truth sample, it can provide better guidance when used for preview.\n\n\nFinally, **we respectfully disagree with the comments on limited contribution**. We would like to emphasize our main contributions as follows:\n\n- Theoretically, we are the **first** to systematically study the parameterization problem of exponential integrators for diffusion sampling, the **first** to explain the superiority of previous parameterizations, and the **first** to introduce Rosenbrock-type methods into diffusion sampling. We also propose the pseudo-order solver and the half-corrector technique, which are **novel, insightful and efficient**. These ideas may inspire future works for training-free sampling of diffusion models, such as SDE-based samplers.\n- Empirically, we achieve the **most SOTA sample quality among all the training-free sampling methods** on a fair benchmark, under various settings, against the most recent and strong baselines. Our methods have **high practical value** in the AIGC community regarding the performance at small NFEs and the applicability to the text-to-image Stable Diffusion, for the usage of **generating previews** and high-quality samples of text-to-image diffusion models with low costs.\n"}, {"review_id": "pMvowp5yYO", "reviewer": "Reviewer_wxE7", "summary": "The paper proposes a new sampling algorithm, DPM-Solver-v3, to solve the diffusion ODE. Applying the Rosenbrock-type exponential integrators, DPM-Solver-v3 pre-computes several coefficients (empirical model statistics) to minimize the norm of the gradient of non-linear term, which leads to reduced discretization error. Empirical results of DPM-Solver-v3 demonstrate its effectiveness in unconditional/conditional sampling and classifier-free guidance on Stable Diffusion.", "questions": "- The author mentioned that EMS computed on the model without guidance can work for guided sampling cases within a common range of guidance scales. I'm wondering how large is the performance gap between algorithms using EMS computed on the model with guidance and the model without guidance. How is this affected by different magnitudes of guidance scales?\n\n- Theoretically and empirically, is the model robust to estimation error of EMS?\n\n- How is DPM-Solver-v3's performance affected by single-step versus multistep sampling, in comparison with the baselines?\n\n- The estimated $s_{\\lambda}$ seems to be unstable w.r.t $\\lambda$ in some cases. Would using moving average or other smoothing methods help with the sampling algorithm?\n\n- I think it would be very helpful to compare the FID score / CLIP score on Stable-Diffusion since low MSE does not necessarily imply better sample quality (the quality is affected not only by discretization error but also by estimation error of neural networks).\n\n- How well is DPM-Solver-v3 compatible with current SDE samplers?", "limitations": "yes", "rating": 5, "confidence": 4, "soundness": 3, "presentation": 2, "contribution": 2, "strengths": "- By introducing empirical model statistics arising from Rosenbrock-type exponential integrators, DPM-Solver-v3 has better model parametrization than data/noise prediction ones, which leads to lower discretization error.\n\n- Experiments show that DPM-Solver-v3 outperforms DPM-Solver++ and UniPC for both conditional and unconditional sampling. DPM-Solver-v3 also incurs a smaller MSE on Stable Diffusion, across NFEs.", "weaknesses": "- The sampling algorithm is based on the pre-computed empirical model statistics $l_{\\lambda}, s_{\\lambda}, b_{\\lambda}$ at each time step $\\lambda$, and their corresponding estimated integrals. This seems to lead to extra memory cost and makes the model not able to be flexibly adapted to new use cases (e.g., a different time schedule, guided sampling, etc.)\n\n- When developing the error order and convergence of the algorithm in section 3.2, the estimation error of EMS is not taken into account. However, the unstable $s_{\\lambda}$ w.r.t $\\lambda$ in Figure 5 (especially for the one estimated on ScoreSDE) indicates that the EMS estimation error is not ignorable. Is it possible to provide such analysis or guarantees given the error of EMS?\n\n- In the low FID region (which is closer to practical usage), the benefit of DPM-Solver-v3 seems not very significant compared with UniPC in Figure 3 and Figure 4\n\n- Applying Rosenbrock-type exponential integrators to diffusion ODE is a valid contribution, but other components (e.g. higher-order derivative estimation) of the method seem to be quite relevant to those proposed by existing works like the RK45 Butcher tableau or [1,2]. The real improvement of the Rosenbrock-type technique is unclear when compounded with many other components. Could the authors maybe ablate on its individual effect?\n\n- I think it would be very helpful to compare the FID score / CLIP score on Stable-Diffusion since low MSE does not necessarily imply better sample quality (the quality is affected not only by discretization error but also by estimation error of neural networks).\n\n- Could the authors also compares DPM-Solver-V3 with some typical stochastic samplers, like Alg 2 in EDM, Gonna Go Fast [3], or more recent works on stochastic samplers? Is it possible to integrate the solver in these stochastic samplers, as these samplers seem to provide significantly better sample quality compared to their deterministic counterparts on complex datasets?\n\n\n[1] Fast Sampling of Diffusion Models with Exponential Integrator, Zhang et al., ICML 2023\n\n[2] DPM-Solver: A Fast ODE Solver for Diffusion Probabilistic Model Sampling in Around 10 Steps, Lu et al., NeurIPS 2022\n\n[3] Gotta go fast when generating data with score-based models.  A. Jolicoeur-Martineau, K. Li, R. Piché-Taillefer, T. Kachman, and I. Mitliagkas. CoRR, abs/2105.14080, 2021.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for the interest and acknowledgment of our theoretical and empirical contributions. Below are our explanations to the questions, which we hope may clarify some misunderstandings. We kindly request that you consider raising the score accordingly if you are satisfied.\n\n*W1: The method needs extra memory cost and makes the model not able to be flexibly adapted to new use cases (e.g., a different time schedule, guided sampling, etc.)*\n\n**A**: **The extra memory cost is rather small**. Please refer to *common response, Q1*; And **the EMS are flexible for your mentioned cases**. Please refer to *common response, Q2*.\n\n*W2: In the low FID region, the benefit of DPM-Solver-v3 seems not very significant compared with UniPC in Figure 3 and Figure 4.*\n\n**A**: When NFE is around 20, our improvement of sample quality is small because all different fast samplers almost converge. However, what matters is that our method has a **faster convergence speed** to good sample quality. A notable example is that, on LSUN-Bedroom, we reach the FID of 3.06 with 12 NFE, while the previous best method requires 20 NFE, which means our computation cost is approximately 60%.\n\n*W3: Could the authors ablate on the individual effect of Rosenbrock-type exponential integrators?*\n\n**A**: Other components (e.g. higher-order derivative estimation) are just what we formulate to adapt high-order solvers to our new parameterization. Therefore, the EMS estimation itself is what we are different from existing works. If we set $l_\\lambda,s_\\lambda,b_\\lambda$ to special values, our solvers degenerate to previous ones such as DPM-Solver and DPM-Solver++ (Appendix A.2). When we compare our method to DPM-Solver++/UniPC, we are already ablating on the EMS' effect.\n\n*Q1: The performance gap between EMS computed with guidance and without guidance? How is this affected by different magnitudes of guidance scales?*\n\n**A**: Empirically, the EMS computed on the model without guidance (unconditional part) performs more stably than those computed on the model with guidance. See *common response, Q2* for details.\n\n*Q2: Theoretically and empirically, is the model robust to estimation error of EMS?*\n\n**A**: Theoretically, **the order and convergence theorems are irrelevant to the EMS estimation error**. The ODE solution Eq. (9) is correct whatever $l_\\lambda,s_\\lambda,b_\\lambda$ are, and we only need the assumption that these coefficients are bounded (Assumption B.2 in Appendix B.1.1) to prove the local and global order; Empirically, we need to ensure the amount of datapoints for estimating EMS (see examples in Appendix G.1), and we find that our method is robust to the estimation error of EMS **given only 1024 datapoints**.\n\n*Q3: Single-step vs multistep sampling?*\n\n**A**: We implement the singlestep DPM-Solver-v3 and compare it with the singlestep version of DPM-Solver++, since UniPC only has multistep version. We report the FID on CIFAR10 (ScoreSDE). See *response pdf, Table 1* for details.\n\nThe results show that DPMv3 (singlestep) is better than DPM++ (singlestep), but multistep version is better than singlestep and can achieve the best results.\n\n*Q4: The estimated $s_\\lambda$ seems to be unstable. Would moving average help with the sampling algorithm?*\n\n**A**: We want to clarify that the unstable $s_\\lambda$ is not an issue, and **our sampler is stable**:\n- The instability of $s_\\lambda$ on ScoreSDE is intrinsic and not due to the estimation error. As we increase the number of samples to decrease the estimation error, the fluctuation is not reduced. We attribute it to the periodicity of trigonometric functions as stated in Sec.4.2.\n- Moreover, as we only consider the integrals of $s_\\lambda$ and the sign of $s_\\lambda$ does not change much, the integrals of $s_\\lambda$ is rather smooth, which ensures the stability of our method.\n\nTherefore, it is no need for smoothing since the involved form is not $s_\\lambda$ but the integral of $s_\\lambda$, and it's actually already smooth.\n\n*Q5: Compare the FID score / CLIP score on Stable-Diffusion since low MSE does not necessarily imply better sample quality.*\n\n**A**: We choose the MSCOCO2014 validation set as reference, and compute the FID/CLIP score with 10k samples of cfg scale 7.5. See detailed results in *response pdf, Table 2*. \n\nThe results show that our methods achieve consistently better FID and similar CLIP scores. **Notably, we achieve an FID of 15.4 in 8 NFE, close to the reported FID of Stable-Diffusion v1.4**. \n\nStill, we claim that FID is not a proper metric for evaluating the convergence of latent-space diffusion models. As stated in DPM-Solver++, The FIDs quickly achieves 15.0~16.0 within 10 steps, even if the latent code does not converge, because of the strong image decoder. Instead, MSE in the latent space is a direct way to measure the convergence. By comparing the MSE, our sampler does converge faster to the ground-truth samples of Stable Diffusion itself.\n\n*Q6: Compare DPM-Solver-V3 with stochastic samplers? Is it possible to integrate the solver in these stochastic samplers?*\n\n**A**: \n- **Stochastic samplers have bad performance for NFE <= 20**. SDE solvers are good at best-quality sampling, not fast sampling. In the papers of SDE solvers, the reported NFE is usually greater than 100, and the sample quality is quite low at around 20 NFE, even using the recent advanced SDE solver (Appendix E.7 in [1]).\n- Yes, our EMS can also be used in SDE solvers. Note that previous SDE solvers also limit their parameterization to noise/data prediction, and we can find better parameterization by the proposed EMS. Therefore, It is promising to improve SDE solvers with our EMS for both fast and best-quality sampling.\n\n[1] Gonzalez, Martin, et al. \"SEEDS: Exponential SDE Solvers for Fast High-Quality Sampling from Diffusion Models.\"\n\nOnce again, thank you for your constructive feedback and for considering our paper for acceptance. We'll revise our paper according to your suggestions."}]

Author rebuttal:
We sincerely thank all reviewers' effort for the detailed and insightful suggestions. First we'd like to present all the paper revisions.

### Summary of paper revisions

**Theories**

- Theoretical justification of why minimizing first-order discretization error can help higher-order solvers (Reviewer dqUr, W2)
- Theoretical analysis on the extra error of the trapezoidal rule (Reviewer dqUr, W3)

**Experiments**

- DPM-Solver-v3’s performance affected by single-step versus multistep sampling (Reviewer wxE7, Q3)
- FID score / CLIP score on Stable-Diffusion (Reviewer wxE7, Q5)

**Writing**

- Rephrase the duplicate expressions from DPM-Solver (Reviewer dqUr, W4)
- Clarify which is the unique contribution in Section 3.2 (Reviewer dqUr, W5)
- Reiterate the tuning strategies in Appendix Table 4/Table 7 (Reviewer yHtS, W1/W2)
- Moving the quantitative results to the main body from the appendix (Reviewer yHtS, W3)
- Clarify the fluctuation of $s_\lambda$ is not due to estimation error, and there is no need for moving average or other smoothing methods (Reviewer wxE7, W2/Q4)
- Discuss why the EMS of unconditional model can be applied to guided sampling (Reviewer dqUr/wxE7, Q1)
- Discuss why 2nd-order is better for conditional setting (Reviewer dqUr, Q2)
- Discuss the practical value of 15\~30% speed-up and improvement in 5\~10 NFE (Reviewer oHv2, Q1)

We find there are common misunderstandings to our paper, and we'd like to clarify here again.

### Response for common questions

*Q1: Extra computation/memory costs of the proposed method? (from dqUr,oHv2,wxE7,wxE7)*

**A**: **The extra memory cost is rather small**. The extra coefficients are discretized and computed at $N$ timesteps, each with a dimension $D$ same as the diffused data. The extra memory cost is $O(ND)$, including the precomputed terms in Appendix C.1.2, and is rather small compared to the pretrained model (e.g. only ~125M in total on Stable-Diffusion, compared to ~4G of the model itself).

**The pre-computation time for estimating EMS is rather short**. The EMS introduced by our method can be effectively estimated on around 1k datapoints within hours, which is rather short comparing to the long training / distillation time of other methods. Moreover, the integrals of these extra coefficients are just some vector constants which can be pre-computed **within seconds**, as shown in Appendix C.1.2. The precomputing is done only once before sampling.

**The extra computational overhead during sampling is negligible**. Once we obtain the estimated EMS and their integrals at discrete timesteps, they can be regarded as constants. Thus, during the subsequent sampling process, **the computational overhead is the same as previous training-free methods** (such as DPM-Solver++) with negligible difference (Appendix E).

*Q2: How does EMS apply for different time schedules and guided sampling? (from dqUr,yHtS,wxE7)*

**A**: The pre-computed EMS can be applied for any time schedules during sampling without re-computing EMS. Besides, we compute EMS on unconditional models and it can be used for a wide range of guidance scales (such as cfg=7.5 in Stable Diffusion). In short, EMS is flexible and easy to adopt in downstream applications.

- **Time schedule**: The choice of $N,K$ in EMS is **disentangled** with the timestep scheduler in sampling. Once we have estimated the EMS at $N$ (e.g., 1200) timesteps, they can be flexibly adapted to any schedule (uniform logSNR/uniform t...) during sampling, by corresponding the actual timesteps during sampling to the $N$ bins. For different time schedule, we only need to re-precompute $E_{\lambda_s,\lambda_t}^{(k)}$ in Appendix C.1.2, and the time cost is within seconds.
- **Guided sampling**. We compute the EMS on unconditional model for all guided cases. Empirically, the EMS computed on the model without guidance (unconditional part) performs more stably than those computed on the model with guidance, and can accelerate the sampling procedure in a wide range of guidance scales (including the common guidance scales used in pretrained models). We think that the unconditional model contains some common information (image priors) for all the conditions, such as color, sketch, and other image patterns. Extracting them helps correct some common bias such as shallow color, low saturation level and lacking of details. In contrast, the conditional model is dependent on the condition and has large variance.
    - In addition, EMS computed on models without guidance cannot work for extremely large guidance scale (e.g., cfg scale 15 for Stable Diffusion), since in this case the condition has a large impact on the denoising process. Note that at these extreme large scales, the sampling quality is very low (compared to cfg scale 7.5) and they are rarely used in practice. Therefore, our proposed EMS without guidance is suitable enough for the common applications with the common guidance.

We think the the quality of our work has been improved a lot. We are welcome for further questions.

Meta-review: This paper proposes an improved/generalized version of DPM-solver with additional coefficients that are optimized based on trained model. Theoretical analysis has been carried out to ensure local/global convergence, and extensive experiments have been performed to justify the effectiveness. Reviewers are overall positive on the soundness, presentation, and contribution of this work.
Decision: Accept (poster)