PFDiff: Training-free Acceleration of Diffusion Models through the Gradient Guidance of Past and Future

Venue: NeurIPS_cc_2024_Conference
Published: 2024-05-11 18:12:28

Abstract: Diffusion Probabilistic Models (DPMs) have shown remarkable potential in image generation, but their sampling efficiency is hindered by the need for numerous denoising steps. Most existing solutions accelerate the sampling process by proposing fast ODE solvers. However, the inevitable discretization errors of the ODE solvers are significantly magnified when the number of function evaluations (NFE) is fewer. In this work, we propose PFDiff, a novel training-free and orthogonal timestep-skipping strategy, which enables existing fast ODE solvers to operate with fewer NFE. Based on two key observations: a significant similarity in the model's outputs at time step size that is not excessively large during the denoising process of existing ODE solvers, and a high resemblance between the denoising process and SGD. PFDiff, by employing gradient replacement from past time steps and foresight updates inspired by Nesterov momentum, rapidly updates intermediate states, thereby reducing unnecessary NFE while correcting for discretization errors inherent in first-order ODE solvers. Experimental results demonstrate that PFDiff exhibits flexible applicability across various pre-trained DPMs, particularly excelling in conditional DPMs and surpassing previous state-of-the-art training-free methods. For instance, using DDIM as a baseline, we achieved 16.46 FID (4 NFE) compared to 138.81 FID with DDIM on ImageNet 64x64 with classifier guidance, and 13.06 FID (10 NFE) on Stable Diffusion with 7.5 guidance scale.

Main paper: 
 introduction: 
 In recent years, Diffusion Probabilistic Models (DPMs) [1] [2] [3] [4] have demonstrated exceptional modeling capabilities across various domains including image generation [5] [6] [7] , video generation [8] , text-to-image generation [9, 10] , speech synthesis [11] , and text-to-3D generation [12, 13] . They have become a key driving force advancing deep generative models. DPMs initiate with a forward process that introduces noise onto images, followed by utilizing a neural network to learn a backward process that incrementally removes noise, thereby generating images [2, 4] . Compared to other generative methods such as Generative Adversarial Networks (GANs) [14] and Variational Autoencoders (VAEs) [15] , DPMs not only possess a simpler optimization target but also are capable of producing higher quality samples [5] . However, the generation of high-quality samples via DPMs requires hundreds or thousands of denoising steps, significantly lowering their sampling efficiency and becoming a major barrier to their widespread application.Existing techniques for rapid sampling in DPMs primarily fall into two categories. First, trainingbased methods [16] [17] [18] [19] , which can significantly compress sampling steps, even achieving single-step sampling [19] . However, this compression often comes with a considerable additional training cost, and these methods are challenging to apply to large pre-trained models. Second, training-free samplers [20] [21] [22] [23] [24] [25] [26] [27] [28] [29] [30] , which typically employ implicit or analytical solutions to Stochastic Differential Equations Submitted to 38th Conference on Neural Information Processing Systems (NeurIPS 2024). Do not distribute. (b) Results from Guided-Diffusion [5] on ImageNet 64x64 [32] (Classifier Guidance, s = 1.0)OursFigure 1 : Sampling by conditional pre-trained DPMs [5, 9] using DDIM [20] and our method PFDiff (dashed box) with DDIM as a baseline, varying the number of function evaluations (NFE).(SDE)/Ordinary Differential Equations (ODE) for lower-error sampling processes. For instance, Lu et al. [21, 22] , by analyzing the semi-linear structure of the ODE solvers for DPMs, have sought to analytically derive optimally the solutions for DPMs' ODE solvers. These training-free sampling strategies can often be used in a plug-and-play fashion, compatible with existing pre-trained DPMs.However, when the NFE is below 10, the discretization error of these training-free methods will be significantly amplified, leading to convergence issues [21, 22] , which can still be time-consuming.To further enhance the sampling speed of DPMs, we have analyzed the potential for improvement in existing training-free accelerated methods. Initially, we observed a notably high similarity in the model's outputs for the existing ODE solvers of DPMs when time step size ∆t is not extremely large, as illustrated in Fig. 2a . This observation led us to utilize the gradients that have been computed from past time steps to approximate current gradients, thereby reducing unnecessary estimation of noise network. Furthermore, due to the similarities between the sampling process of DPMs and Stochastic Gradient Descent (SGD) [33] as noted in Remark 1, we incorporated a foresight update mechanism using Nesterov momentum [34] , known for accelerating SGD training. Specifically, we ingeniously employ prior observation to predict future gradients, then utilize the future gradients as a "springboard" to facilitate larger update step size ∆t, as shown in Fig. 2b . where η in DDPM [2] comes from σt in Eq. (6) . Solid lines: ODE solvers, dashed lines: SDE solvers. (b) Comparison of partial sampling trajectories between PFDiff-1 and a first-order ODE solver, where the update directions are guided by the tangent direction of the sampling trajectories.sampling quality for speed [28] , we prove that PFDiff corrects for errors in the sampling trajectories of first-order ODE solvers. This improves sampling quality while reducing unnecessary NFE in existing ODE solvers, as illustrated in Fig. 2b . To validate the orthogonality and effectiveness of PFDiff, extensive experiments were conducted on both unconditional [2, 4, 20] and conditional [5, 9] pre-trained DPMs, with the visualization experiment of conditional DPMs depicted in Fig. 1 . The results indicate that PFDiff significantly enhances the sampling performance of existing ODE solvers.Particularly in conditional DPMs, PFDiff, using only DDIM as the baseline, surpasses the previous state-of-the-art training-free sampling algorithms. 
diffusion sdes: 
 Diffusion Probabilistic Models (DPMs) [1] [2] [3] [4] aim to generate D-dimensional random variables x 0 ∈ R D that follow a data distribution q(x 0 ). Taking Denoising Diffusion Probabilistic Models (DDPM) [2] as an example, these models introduce noise to the data distribution through a forward process defined over discrete time steps, gradually transforming it into a standard Gaussian distribution x T ∼ N (0, I). The forward process's latent variables {x t } t∈[0,T ] are defined as follows:EQUATIONwhere α t is a scalar function related to the time step t, with α 2 t + σ 2 t = 1. In the model's reverse process, DDPM utilizes a neural network model p θ (x t-1 | x t ) to approximate the transition probability q(x t-1 | x t , x 0 ),EQUATIONwhere σ 2 θ (t) is defined as a scalar function related to the time step t. By sampling from a standard Gaussian distribution and utilizing the trained neural network, samples following the data distributionp θ (x 0 ) = T t=1 p θ (x t-1 | x t ) can be generated.Furthermore, Song et al. [4] introduced SDE to model DPMs over continuous time steps, where the forward process is defined as:EQUATIONwhere w t represents a standard Wiener process, and f and g are scalar functions of the time step t.It's noteworthy that the forward process in Eq. ( 1) is a discrete form of Eq. ( 3), where f (t) = d log αt dt and g 2 (t) = dσ 2 t dt -2 d log αt dt σ 2 t . Song et al. [4] further demonstrated that there exists an equivalent reverse process from time step T to 0 for the forward process in Eq. ( 3):EQUATIONwhere w denotes a standard Wiener process. In this reverse process, the only unknown is the score function ∇ x log q t (x t ), which can be approximated through neural networks.In DPMs based on SDE, the discretization of the sampling process often requires a significant number of time steps to converge, such as the T = 1000 time steps used in DDPM [2] . This requirement primarily stems from the randomness introduced at each time step by the SDE. To achieve a more efficient sampling process, Song et al. [4] utilized the Fokker-Planck equation [35] to derive a probability flow ODE related to the SDE, which possesses the same marginal distribution at any given time t as the SDE. Specifically, the reverse process ODE derived from Eq. ( 3) can be expressed as:EQUATIONUnlike SDE, ODE avoids the introduction of randomness, thereby allowing convergence to the data distribution in fewer time steps. Song et al. [4] employed a high-order RK45 ODE solver [36] ,achieving sample quality comparable to SDE at 1000 NFE with only 60 NFE. Furthermore, research such as DDIM [20] and DPM-Solver [21] explored discrete ODE forms capable of converging in fewer NFE. For DDIM, it breaks the Markov chain constraint on the basis of DDPM, deriving a new sampling formula expressed as follows:EQUATIONwhereσt = η (1 -α t-1 ) / (1 -α t ) 1 -α t /α t-1, and α t corresponds to α 2 t in Eq. ( 1). When η = 1, Eq. ( 6) becomes a form of DDPM [2] ; when η = 0, it degenerates into an ODE, the form adopted by DDIM [20] , which can obtain high-quality samples in fewer time steps.Remark 1. In this paper, we regard the gradient dx t , the noise network output ϵ θ (x t , t), and the score function ∇ x log q t (x t ) as expressing equivalent concepts. This is because Song et al. [4] demonstrated that ϵ θ (x t , t) = -σ t ∇ x log q t (x t ). Moreover, we have discovered that any first-order solver of DPMs can be parameterized as x t-1 = xt -γ t dx t + ξϵ t . Taking DDIM [20] as an example, wherext = αt-1 αt x t , γ t = αt-1 αt -α t-1 - √ 1 -α t-1 , dx t = ϵ θ (x t , t), and ξ = 0.This indicates the similarity between SGD and the sampling process of DPMs, a discovery also implicitly suggested in the research of Xue et al. [30] and Wang et al. [37] .3 Method 
solving for reverse process diffusion odes: 
 By substituting ϵ θ (x t , t) = -σ t ∇ x log q t (x t ) [4] , Eq. ( 5) can be rewritten as:EQUATIONGiven an initial value x T , we define the time steps {t i } T i=0 to progressively decrease from t 0 = T to t T = 0. Let xt0 = x T be the initial value. Using T steps of iteration, we compute the sequence{x ti } T i=0to obtain the solution of this ODE. By integrating both sides of Eq. ( 7), we can obtain the exact solution of this sampling ODE.EQUATIONFor any p-order ODE solver, Eq. ( 8) can be discretely represented as:EQUATIONwhere t0 = t i-1 , tp = t i , and ∆ t = tn+1 -tn denote the time step size. The function h represents the different solution methodologies applied by various p-order ODE solvers to the function s. For the Euler-Maruyama solver [38] , h is the identity mapping of s. Further, we defineϕ(Q, xti-1 , t i-1 , t i ) := xti-1 + p-1 n=0 h(ϵ θ (x tn , tn ), xt n , tn ) • ∆ t.Here, ϕ is any p-order ODE solver, and buffer Q = ϵ θ (x tn , tn )p-1 n=0 , t i-1 , t i , where t0 = t i-1 and tp = t i .When using the ODE solver defined in Eq. ( 9) for sampling, the choice of T = 1000 leads to significant inefficiencies in DPMs. The study on DDIM [20] first revealed that by constructing a new forward sub-state sequence of length M + 1 (M ≤ T ), {x ti } M i=0 , from a subsequence of time steps [0, . . . , T ] and reversing this sub-state sequence, it is possible to converge to the data distribution in fewer time steps. However, as illustrated in Fig. 2a , for ODE solvers, as the time step ∆t = t i -t i-1increases, the gradient direction changes slowly initially, but undergoes abrupt changes as ∆t → T .This phenomenon indicates that under minimal NFE (i.e., maximal time step size ∆t) conditions, the discretization error in Eq. ( 9) is significantly amplified. Consequently, existing ODE solvers, when sampling under minimal NFE, must sacrifice sampling quality to gain speed, making it an extremely challenging task to reduce NFE to below 10 [21, 22] . Given this, we aim to develop an efficient timestep-skipping sampling algorithm, which reduces NFE while correcting discretization errors, thereby ensuring that sampling quality is not compromised, and may even be improved.For any p-order timestep-skipping sampling algorithm for DPMs, the sampling process can be reformulated according to Eq. ( 9) as follows:EQUATIONwhere bufferQ = ϵ θ (x tn , tn ) p-1 n=0 , t i-1 , t i and [1, . . . , M ] is an increasing subsequence of [1, . . . , T ].As illustrated in Fig. 2a , when the time step size ∆t (i.e., t i -t i-1 ) is not excessively large, the MSE of the noise network, defined as 1T -∆t T -∆t-1 t=0 ∥ϵ θ (x t , t) -ϵ θ (x t+∆t , t + ∆t)∥ 2 , isremarkably similar. This phenomenon is especially pronounced in ODE-based sampling algorithms, such as DDIM [20] and DPM-Solver [21] . This observation suggests that there are many unnecessary time steps in ODE-based sampling methods during the complete sampling process (e.g., when T = 1000), which is one of the reasons these methods can generate samples in fewer steps. Based on this, we propose replacing the noise network of the current timestep with the output from a previous timestep to reduce unnecessary NFE without compromising the quality of the final generated samples.Initially, we store the output of the previous timestep's noise network in a buffer as follows:EQUATIONThen, in the current timestep, we directly use the noise network output saved in the buffer from the previous timestep to replace the current timestep's noise network output, thereby updating the intermediate states to the next timestep, as detailed below:EQUATIONBy using this approach, we can effectively accelerate the sampling process, reduce unnecessary NFE, and ensure the quality of the samples is not affected. The convergence proof is in Appendix B.1.As stated in Remark 1, considering the similarities between the sampling process of DPMs and SGD [33] , we introduce a foresight update mechanism of Nesterov momentum, utilizing future gradient information as a "springboard" to assist the current intermediate state in achieving more efficient leapfrog updates. Specifically, for the intermediate state xti+1 predicted using past gradients as discussed in Sec. 3.2, we first estimate the future gradient and update the current buffer as follows:EQUATIONSubsequently, leveraging the concept of foresight updates, we predict a further future intermediate state xti+2 using the current intermediate state xti along with the future gradient information corresponding to xti+1 , as shown below:EQUATIONFurthermore, Zhou et al. [39] (13) . Furthermore, we propose that in a single iteration process, xti+2 in Eq. ( 14) can be modified to xt i+(k+1) , achieving a k-step skip to sample more distant future intermediate states. Additionally, when k ̸ = 1, the buffer Q, which acts as an intermediate "springboard" from Eq. ( 13), has various computational origins. This can be accomplished by modifying xti+1 in Eq. ( 12) to xt i+l . We collectively refer to this multi-step skipping and different "springboard" selection strategy as PFDiff-k_l (l ≤ k). Further algorithmic details can be found in Appendix C. Finally, through the comparison of sampling trajectories between PFDiff-1 and a first-order ODE sampler, as shown in Fig. 2b , PFDiff-1 showcases its capability to correct the sampling trajectory of the first-order ODE sampler while reducing the NFE.Proposition 3.1. For any given DPM first-order ODE solver ϕ, the PFDiff-k_l algorithm can describe the sampling process within an iteration cycle through the following formula:EQUATIONAlgorithm 1 PFDiff-1Require: initial value x T , NFE N , model ϵ θ , any p-order solver ϕ 1: Define time steps{t i } M i=0 with M = 2N -1p 2: xt0 ← x T 3: Q buffer ← ---ϵ θ (x tn , tn ) p-1 n=0 , t 0 , t 1 , where t0 = t 0 , tp = t 1 ▷ Initialize buffer 4: xt1 = ϕ(Q, xt0 , t 0 , t 1 ) 5: for i ← 1 to M p -2 do 6:if (i -1) mod 2 = 0 then where the value of ϵ θ (x t i-(k-l+1) , t i-(k-l+1) ) can be directly obtained from the buffer Q, without the need for additional computations. The iterative process defined by Eq. (15) ensures that the sampling outcomes converge to the data distribution consistent with the solver ϕ, while effectively correcting errors in the sampling process (Proof in Appendix B).7: xti+1 = ϕ(Q, xtiIt is noteworthy that, although the PFDiff is conceptually orthogonal to the SDE/ODE solvers of DPMs, even when the time size ∆t is relatively small, the MSE of the noise network in the SDE solver exhibits significant differences, as shown in Fig. 2a . Consequently, PFDiff shows marked improvements on the ODE solver, and our experiments are almost exclusively based on ODE solvers, with exploratory experiments on SDE solvers referred to Sec. 4.1.Relationship with p-order solver [21, 22, 27] . According to Eq. ( 10), a single iteration of the p-order solver can be represented as:EQUATIONA single iteration of the p-order solver uses p NFE to predict the next intermediate state. The intermediate step gradients obtained during this process can be considered as an approximation of future gradients. This approximation is implicitly contained within the sampling guided by future gradients that we propose. Furthermore, as shown in Eq. ( 15), a single iteration update of PFDiff based on a first-order solver can be seen as using a 2-order solver with only one NFE. 
experiments: 
 In this section, we validate the effectiveness of PFDiff as an orthogonal and training-free sampler through a series of extensive experiments. This sampler can be integrated with any order of ODE solvers, thereby significantly enhancing the sampling efficiency of various types of pre-trained DPMs.To systematically showcase the performance of PFDiff, we categorize the pre-trained DPMs into two main types: conditional and unconditional. Unconditional DPMs are further subdivided into discrete and continuous, while conditional DPMs are subdivided into classifier guidance and classifier-free guidance. In choosing ODE solvers, we utilized the widely recognized first-order DDIM [20] ,Analytic-DDIM [23] , and the higher-order DPM-Solver [21] as baselines. For each experiment, we use the Fréchet Inception Distance (FID↓) [40] as the primary evaluation metric, and provide the experimental results of the Inception Score (IS↑) [41] in the Appendix D.7 for reference. Lastly, apart from the ablation studies on parameters k and l discussed in Sec. 4.3, we showcase the optimal results of PFDiff-k_l (where k = 1, 2, 3 and l ≤ k) across six configurations as a performance demonstration of PFDiff. As described in Appendix C, this does not increase the computational burden in practical applications. All experiments were conducted on an NVIDIA RTX 3090 GPU.For unconditional DPMs, we selected discrete DDPM [2] and DDIM [20] , as well as pre-trained models from continuous ScoreSDE [4] , to assess the effectiveness of PFDiff. For these pre-trained models, all experiments sampled 50k instances to compute evaluation metrics.For unconditional discrete DPMs, we first select first-order ODE solvers DDIM [20] and Analytic-DDIM [23] as baselines, while implementing SDE-based DDPM [2] and Analytic-DDPM [23] methods for comparison, where η = 1.0 is from σt in Eq. ( 6). We conduct experiments on the CIFAR10 [42] and CelebA 64x64 [43] cannot achieve. Furthermore, in Appendix D.2, by varying η from 1.0 to 0.0 in Eq. ( 6) to control the scale of noise introduced by SDE, we observe that as η decreases (reducing noise introduction), the performance of PFDiff gradually improves. This once again validates our assumption proposed in Sec. 3.2, based on Fig. 2a , that there is a significant similarity in the model's outputs at the time step size that is not excessively large for the existing ODE solvers.For unconditional continuous DPMs, we choose the DPM-Solver-1, -2 and -3 [21] as the baseline to verify the effectiveness of PFDiff as an orthogonal timestep-skipping algorithm on the first and higher-order ODE solvers. We conducted experiments on the CIFAR10 [42] using quadratic time steps, varying the NFE. The experimental results using FID↓ as the evaluation metric are shown in Fig. 3c . More experimental details can be found in Appendix D.3. We observe that PFDiff consistently improves the sampling performance over the baseline with fewer NFE settings, particularly in cases where higher-order ODE solvers fail to converge with a small NFE (below 10) [21] .For conditional DPMs, we selected the pre-trained models of the widely recognized classifier guidance paradigm, ADM-G [5] , and the classifier-free guidance paradigm, Stable-Diffusion [9] , to validate the effectiveness of PFDiff. We employed uniform time steps setting and the DDIM [20] ODE solver as a baseline across all datasets. Evaluation metrics were computed by sampling 50k samples on the ImageNet 64x64 [32] dataset for ADM-G and 10k samples on other datasets, including ImageNet 256x256 [32] in ADM-G and MS-COCO2014 [31] in Stable-Diffusion.For conditional DPMs employing the classifier guidance paradigm, we conducted experiments on the ImageNet 64x64 dataset [32] with a guidance scale (s) set to 1.0. For comparison, we implemented DPM-Solver-2 and -3 [21] , and DPM-Solver++(2M) [22] , which exhibit the best performance on conditional DPMs. Additionally, we introduced the AutoDiffusion method [26] using DDIM as a baseline for comparison, noting that this method incurs additional search costs. We compared FID↓ scores by varying the NFE as depicted in Fig. 4a , with corresponding visual comparisons shown in Fig. 1b . We observed that PFDiff reduced the FID from 138.81 with 4 NFE in DDIM to 16.46, achieving an 88.14% improvement in quality. The visual results in Fig. 1b further demonstrate that, at the same NFE setting, PFDiff achieves higher-quality sampling. Furthermore, we evaluated PFDiff's sampling performance based on DDIM on the large-scale ImageNet 256x256 dataset [32] . Detailed results are provided in Appendix D.4.For conditional, classifier-free guidance paradigms of DPMs, we employed the sd-v1-4 checkpoint and computed the FID↓ scores on the validation set of MS-COCO2014 [31] . We conducted experiments with a guidance scale (s) set to 7.5 and 1.5. For comparison, we implemented DPM-Solver-2 and -3 [21] , and DPM-Solver++(2M) [22] methods. At s = 7.5, we introduced the state-of-the-art method reported in DPM-Solver-v3 [27] for comparison, along with DPM-Solver++(2M) [22] , UniPC [29] , and DPM-Solver-v3(2M) [27] . The FID↓ metrics by varying the NFE are presented in Figs. 4b and 4c , with additional visual results illustrated in Fig. 1a . We observed that PFDiff, solely based on DDIM, achieved state-of-the-art results during the sampling process of Stable-Diffusion, thus demonstrating the efficacy of PFDiff. Further experimental details can be found in Appendix D.5.We conducted ablation experiments on the six different algorithm configurations of PFDiff mentioned in Appendix C, with k = 1, 2, 3 (l ≤ k). Specifically, we evaluated the FID↓ scores on the unconditional and conditional pre-trained DPMs [2, 4, 5, 9] by varying the NFE. Detailed experimental setups and results can be found in Appendix D.6.1. The experimental results indicate that for various pre-trained DPMs, the choice of parameters k and l is not critical, as most combinations of k and l within PFDiff can enhance the sampling efficiency over the baseline. Moreover, with k = 1 fixed, PFDiff-1 can significantly improve the baseline's sampling quality within the range of 8∼20 NFE.For even better sampling quality, one can sample a small subset of examples (e.g., 5k) to compute evaluation metrics or directly conduct visual analysis, easily identifying the most effective k and l combinations.To validate the PFDiff algorithm as mentioned in Sec. 3.3, which necessitates the joint guidance of past and future gradients for first-order ODE solvers, and only past gradients for higher-order ODE solvers, offering a more effective means of accelerating baseline sampling. This study employs the first-order ODE solver DDIM [20] as the baseline, isolating the effects of both past and future gradients, and uses the higher-order ODE solver DPM-Solver [21] as the baseline, removing the influence of future gradients for ablation experiments. Specific experimental configurations and results are shown in Appendix D.6.2. The results indicate that, as described by the PFDiff algorithm in Sec. 3.3, it is possible to further enhance the sampling efficiency of ODE solvers of any order. 
conclusion: 
 In this paper, based on the recognition that the ODE solvers of DPMs exhibit significant similarity in model outputs when the time step size is not excessively large, and with the aid of a foresight update mechanism, we propose PFDiff, a novel method that leverages the gradient guidance from both past and future to rapidly update the current intermediate state. This approach effectively reduces the unnecessary number of function evaluations (NFE) in the ODE solvers and significantly corrects the errors of first-order ODE solvers during the sampling process. Extensive experiments demonstrate the orthogonality and efficacy of PFDiff on both unconditional and conditional pre-trained DPMs, especially in conditional pre-trained DPMs where PFDiff outperforms previous state-of-the-art training-free sampling methods.Limitations and broader impact Although PFDiff can effectively accelerate the sampling speed of existing ODE solvers, it still lags behind the sampling speed of training-based acceleration methods and one-step generation paradigms such as GANs. Moreover, there is no universal setting for the optimal combination of parameters k and l in PFDiff; adjustments are required according to different pre-trained DPMs and NFE. It is noteworthy that PFDiff may be utilized to accelerate the generation of malicious content, thereby having a detrimental impact on society. 
experimental setting/details: 
 Question: Does the paper specify all the training and test details (e.g., data splits, hyperparameters, how they were chosen, type of optimizer, etc.) necessary to understand the results?Answer: [Yes] Justification: Our method is training-free, but we also report the hyperparameters used when evaluating our proposed method. Details can be found in section 4 and Appendix D.Guidelines:• The answer NA means that the paper does not include experiments.• The experimental setting should be presented in the core of the paper to a level of detail that is necessary to appreciate the results and make sense of them.• The full details can be provided either with the code, in appendix, or as supplemental material. 
experiment statistical significance: 
 Question: Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?Answer: [No] Justification: We generate 50k images for 32x32, 64x64 datasets and 10k for 256x256 to evaluate the FID metric. According to previous works [20, 21, 26, 27] , when evaluating with the generated samples mentioned above, the standard deviation of the FID evaluations is rather small (mainly less than 0.01). These small standard deviations do not change the conclusions.Guidelines:• The answer NA means that the paper does not include experiments.• The authors should answer "Yes" if the results are accompanied by error bars, confidence intervals, or statistical significance tests, at least for the experiments that support the main claims of the paper.• The factors of variability that the error bars are capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, or overall run with given experimental conditions).• The method for calculating the error bars should be explained (closed form formula, call to a library function, bootstrap, etc.)• The assumptions made should be given (e.g., Normally distributed errors).• It should be clear whether the error bar is the standard deviation or the standard error of the mean.• It is OK to report 1-sigma error bars, but one should state it. The authors should preferably report a 2-sigma error bar than state that they have a 96% CI, if the hypothesis of Normality of errors is not verified.• For asymmetric distributions, the authors should be careful not to show in tables or figures symmetric error bars that would yield results that are out of range (e.g. negative error rates).• If error bars are reported in tables or plots, The authors should explain in the text how they were calculated and reference the corresponding figures or tables in the text. 
experiments compute resources: 
 Question: For each experiment, does the paper provide sufficient information on the computer resources (type of compute workers, memory, time of execution) needed to reproduce the experiments?Answer: [Yes]Justification: In section 4, we mentioned that all experiments were conducted on an NVIDIA RTX 3090 GPU.Guidelines:• The answer NA means that the paper does not include experiments.• The paper should indicate the type of compute workers CPU or GPU, internal cluster, or cloud provider, including relevant memory and storage.• The paper should provide the amount of compute required for each of the individual experimental runs as well as estimate the total compute.• The paper should disclose whether the full research project required more compute than the experiments reported in the paper (e.g., preliminary or failed experiments that didn't make it into the paper). 
code of ethics: 
 Question: Does the research conducted in the paper conform, in every respect, with the NeurIPS Code of Ethics https://neurips.cc/public/EthicsGuidelines? 
broader impacts: 
 Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?Answer: [Yes] Justification: See section 5.Guidelines:• The answer NA means that there is no societal impact of the work performed.• If the authors answer NA or No, they should explain why their work has no societal impact or why the paper does not address societal impact.• Examples of negative societal impacts include potential malicious or unintended uses (e.g., disinformation, generating fake profiles, surveillance), fairness considerations (e.g., deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations.• The conference expects that many papers will be foundational research and not tied to particular applications, let alone deployments. However, if there is a direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that an improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that a generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster.• The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly, harms that could arise when the technology is being used as intended but gives incorrect results, and harms following from (intentional or unintentional) misuse of the technology.• If there are negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition to attacks, mechanisms for monitoring misuse, mechanisms to monitor how a system learns from feedback over time, improving the efficiency and accessibility of ML). 
safeguards: 
 Question: Does the paper describe safeguards that have been put in place for responsible release of data or models that have a high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?Answer: [NA] Justification: Our method is a training-free accelerated sampling approach. It relies on existing pre-trained DPMs and does not involve the release of data or models, thus there is no need for safeguards.Guidelines:• The answer NA means that the paper poses no such risks.• Released models that have a high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere to usage guidelines or restrictions to access the model or implementing safety filters.• Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images.• We recognize that providing effective safeguards is challenging, and many papers do not require this, but we encourage authors to take this into account and make a best faith effort. Guidelines:• The answer NA means that the paper does not use existing assets.• The authors should cite the original paper that produced the code package or dataset.• The authors should state which version of the asset is used and, if possible, include a URL.• The name of the license (e.g., CC-BY 4.0) should be included for each asset.• For scraped data from a particular source (e.g., website), the copyright and terms of service of that source should be provided.• We recognize that the procedures for this may vary significantly between institutions and locations, and we expect authors to adhere to the NeurIPS Code of Ethics and the guidelines for their institution.• For initial submissions, do not include any information that would break anonymity (if applicable), such as the institution conducting the review. 
 
 Appendix: While the solvers for Diffusion Probabilistic Models (DPMs) are categorized into two types, SDE and ODE, most current accelerated sampling techniques are based on ODE solvers due to the observation that the stochastic noise introduced by SDE solvers hampers rapid convergence. ODE-based solvers are further divided into training-based methods [16] [17] [18] [19] and training-free samplers [20] [21] [22] [23] [24] [25] [26] [27] [28] [29] [30] . Trainingbased methods can notably reduce the number of sampling steps required for DPMs. An example of such a method is the knowledge distillation algorithm proposed by Song et al. [19] , which achieves one-step sampling for DPMs. This sampling speed is comparable to that of GANs [14] and VAEs [15] . However, these methods often entail significant additional costs for distillation training. This requirement poses a challenge when applying them to large pre-trained DPM models. Therefore, our work primarily focuses on training-free, ODE-based accelerated sampling strategies.Training-free accelerated sampling techniques based on ODE can generally be applied in a plugand-play manner, adapting to existing pre-trained DPMs. These methods can be categorized based on the order of the ODE solver-that is, the NFE required per sampling iteration-into first-order [20, [23] [24] [25] and higher-order [21, 22, 27, 29, 36] . Typically, higher-order ODE solvers tend to sample at a faster rate, but may fail to converge when the NFE is low (below 10), sometimes performing even worse than first-order ODE solvers. In addition, there are orthogonal techniques for accelerated sampling. For instance, Li et al. [26] build upon existing ODE solvers and use search algorithms to find optimal sampling sub-sequences and model structures to further speed up the sampling process; Ma et al. [28] observe that the low-level features of noise networks at adjacent time steps exhibit similarities, and they use caching techniques to substitute some of the network's low-level features, thereby further reducing the number of required time steps.The algorithm we propose belongs to the class of training-free and orthogonal accelerated sampling techniques, capable of further accelerating the sampling process on the basis of existing first-order and higher-order ODE solvers. Compared to the aforementioned orthogonal sampling techniques, even though the skipping strategy proposed by Ma et al. [28] effectively accelerates the sampling process, it may do so at the cost of reduced sampling quality, making it challenging to reduce the NFE below 50. Although Li et al. [26] can identify more optimal subsampling sequences and model structures, this implies higher search costs. In contrast, our proposed orthogonal sampling algorithm is more efficient in skipping time steps. First, our skipping strategy does not require extensive search costs. Second, we can correct the sampling errors of first-order ODE solvers while reducing the number of sampling steps required by existing ODE solvers, achieving more efficient accelerated sampling.In this section, we prove the convergence of PFDiff and elaborate on how it theoretically corrects first-order ODE solver errors. To delve deeper into PFDiff, we propose the following assumptions: Assumption B.1. When the time step size ∆t = t i -t i-(k-l+1) is not excessively large, the output estimates of the noise network based on the p-order ODE solver at different time steps are approximately the same, that is, ϵ θ (x tn , tn )p-1 n=0 , t i , t i+l ≈ ϵ θ (x tn , tn ) p-1 n=0 , t i-(k-l+1) , t i .Assumption B.2. For the integral from time step t i to t i+(k+1) , originally applies to real-valued functions and does not directly apply to vector-valued functions [45] . However, the study by Zhou et al. [39] using Principal Component Analysis (PCA) on the trajectories of the ODE solvers for DPMs demonstrates that these trajectories are primarily distributed on a two-dimensional plane, which allows the Mean Value Theorem for Integrals to approximately hold under Assumption B.2.t i+(k+1) ti s(ϵ θ (x t , t), x t , t)dt, there exist intermediate time steps t s, t s ∈ (t i , t i+(k+1) ) such that t i+(k+1) ti s(ϵ θ (x t , t), x t , t)dt = s(ϵ θ (x ts , t s), x ts , t s) • (t i+(k+1) -t i ) = h(ϵ θ (x ts , t s ), x ts , t s ) • (t i+(k+1) -t i ) holds,Starting from Eq. ( 8), we consider an iteration process of a p-order ODE solver from xti to xt i+l , where l is the "springboard" choice determined by PFDiff-k_l. This iterative process can be expressed as: xt i+l = xti + t i+l ti s(ϵ θ (x t ,EQUATIONwhere ϕ is any p-order ODE solver. Eq. (B.3) demonstrates that under Assumption B.1, for any porder ODE solver ϕ, PFDiff-k_l utilizes past gradient information as a substitute for current gradient information to update the current intermediate state. This method not only reduces the NFE but also approximates the solution of xt i+l , ensuring convergence to the data distribution corresponding to the solver ϕ. It is important to note that the sampling process described in Eq. (B.3) relies solely on past gradient information and does not estimate the output of the noise network based on the current intermediate state.In particular, within Proposition 3.1 for any first-order (p = 1) ODE solver ϕ, according to Eq.(B.3), we can approximate xt i+l ≈ ϕ(ϵ θ (x t i-(k-l+1) , t i-(k-l+1) ), xti , t i , t i+l ). Thus, Eq. ( 15) can be rewritten as: Based on Eq. ( 8), we consider an iteration process of a first-order (p = 1) ODE solver from xti to xt i+(k+1) , which can be expressed as: xt i+(k+1) ≈ ϕ(ϵ θ (x t i+l , t i+l ), xti , t i , t i+(k+1) ). (B.EQUATIONxt i+(k+1) = xti + t i+(k+1) ti s(ϵ θ (x t , t), x t , t)dt = xti + s(ϵ θ (x ts , t s), xts , t s) • (t i+(k+1) -t i ) = xti + h(ϵ θ (x ts , t s ), xts , t s ) • (t i+(k+1) -t i ) = ϕ(ϵ θ (x ts , t s ), xti , t i , t i+(k+1) ) ≈ ϕ(ϵ θ (x t i+l , t i+l ), xti , t i , t i+(k+1) ), (B.6)where k and l are determined by the selected PFDiff-k_l and the second and third lines are obtained According to the definition of PFDiff-k_l, t i+l also lies within the interval (t i , t i+(k+1) ), and for different combinations of k and l, this means searching and approximating the ground truth t s within the interval (t i , t i+(k+1) ). Among the six different versions of PFDiff-k_l defined in Appendix C, we believe that the optimal t s has been approximated. Compared to the direct discretization of ϵ θ (x t , t) in Eq. (B.5), we corrected the sampling error of the first-order ODE solver and proved its convergence by guiding sampling based on the future gradient information ϵ θ (x t i+l , t i+l ) under Assumption B. As described in Sec. 3.4, during a single iteration, we can leverage the foresight update mechanism to skip to a more distant future. Specifically, we modify Eq. ( 14) to xt i+(k+1) ≈ ϕ(Q, xti , t i , t i+(k+1) )to achieve a k-step skip. We refer to this method as PFDiff-k. Additionally, when k ̸ = 1, the computation of the buffer Q, originating from Eq. ( 13), presents different selection choices. We modify Eq. ( 12) to xt i+l ≈ ϕ(Q, xti , t i , t i+l ), l ≤ k to denote different "springboard" choices with the parameter l. This strategy of multi-step skips and varying "springboard" choices is collectively termed as PFDiff-k_l (l ≤ k). Consequently, based on modifications to parameters k and l in Eq. ( 12) and Eq. ( 14), Eq. ( 13) is updated toQ buffer ← ---ϵ θ (x tn , tn )p-1 n=0 , t i+l , t i+(k+1) , and Eq. ( 11)is updated to Q buffer ← --- ϵ θ (x tn , tn ) p-1 n=0 , t i-(k-l+1) , t i . When k = 1, since l ≤ k, then l = 1,and PFDiff-k_l is the same as PFDiff-1, as shown in Algorithm 1 in Sec. Furthermore, when k = 3, this forms three different versions of PFDiff-3, as shown in Algorithm 3. In this study, we utilize the optimal results from the six configurations of PFDiff-k_l (k = 1, 2, 3(l ≤ k))to demonstrate the performance of PFDiff. As described in Appendix B.2, this is essentially an approximation of the ground truth t s . Through these six different algorithm configurations, we approximately search for the optimal t s . It is important to note that despite using six different algorithm configurations, this does not increase the computational burden in practical applications. This is because, by visual analysis of a small number of generated images or computing specific evaluation metrics, one can effectively select the algorithm configuration with the best performance.Moreover, even without any selection, directly using the PFDiff-1 configuration can achieve significant performance improvements on top of existing ODE solvers, as shown in the ablation study results in Sec. 4.3.Require: initial value x T , NFE N , model ϵ θ , any p-order solver ϕ, skip type l 1: Define time steps{t i } M i=0 with M = 3N -2p 2: xt0 ← x T 3: Q buffer ← ---ϵ θ (x tn , tn ) p-1 n=0 , t 0 , t 1 , where t0 = t 0 , tp = t 1 ▷ Initialize buffer 4: xt1 = ϕ(Q, xt0 , t 0 , t 1 ) 5: for i ← 1 to M p -3 do 6:if (i -1) mod 3 = 0 then 7: if l = 1 then ▷ PFDiff-2_1 8: xti+1 = ϕ(Q, xtiIn this section, we provide further supplements to the experiments on both unconditional and conditional pre-trained Diffusion Probabilistic Models (DPMs) as mentioned in Sec. 4. Through these additional supplementary experiments, we more fully validate the effectiveness of PFDiff as an orthogonal and training-free sampler. Unless otherwise stated, the selection of pre-trained DPMs, choice of baselines, algorithm configurations, GPU utilization, and other related aspects in this section are consistent with those described in Sec. 4.In this section, we list the used datasets, codes, and their licenses in Table 1 .Table 1 : The used datasets, codes, and their licenses.Table 2 : Sample quality measured by FID↓ on the CIFAR10 [42] , CelebA 64x64 [43] , LSUNbedroom 256x256 [44] , and LSUN-church 256x256 [44] using unconditional discrete-time DPMs, varying the number of function evaluations (NFE). Evaluated on 50k samples. PFDiff uses DDIM [20] and Analytic-DDIM [23] as baselines and introduces DDPM [2] and Analytic-DDPM [23] with η = 1.0 from Eq. ( 6 3 , illustrate that with the reduction of η from 1.0 (SDE) to 0.0 (ODE), PFDiff's sampling performance significantly improves at fewer time steps (NFE≤20). The experiment results regarding FID variations with NFE as presented in Table 3 , align with the trends of MSE of noise network outputs with changes in time step size ∆t as depicted in Fig. 2a . This reaffirms the motivation we proposed in Sec. 3.2.In this section, we supplement the specific FID↓ scores for the unconditional, continuous pretrained DPMs models with first-order and higher-order ODE solvers, DPM-Solver-1, -2 and -3, [21] as baselines, as shown in Table 4 . For all experiments in this section, we conducted tests on the CIFAR10 dataset [42] , using the checkpoint checkpoint_8.pth under the vp/cifar10_ddpmpp_deep_continuous configuration provided by ScoreSDE [4] . For the hyperparameter method of DPM-Solver [21] , we adopted singlestep_fixed; to maintain consistency with the discrete-time model in Appendix D.2, the parameter skip was set to time_quadratic (i.e., quadratic time steps). Unless otherwise specified, we used the parameter settings recommended by DPM-Solver. The results in Table 4 show that by using the PFDiff method described in Sec. 3.4 and taking DPM-Solver as the baseline, we were able to further enhance sampling performance on the basis of first-order and higher-order ODE solvers. Particularly, in the 6∼12 NFE range, PFDiff significantly improved the convergence issues of higher-order ODE solvers under fewer NFEs. For instance, at 9 NFE, PFDiff reduced the FID of DPM-Solver-3 from 233.56 to 5.67, improving the sampling quality by 97.57%. These results validate the effectiveness of using PFDiff with first-order or higher-order ODE solvers as the baseline.Table 5 : Sample quality measured by FID↓ on the ImageNet 64x64 [32] and ImageNet 256x256 [32] , using ADM-G [5] model with guidance scales of 1.0 and 2.0, varying the number of function evaluations (NFE). Evaluated: ImageNet 64x64 with 50k, ImageNet 256x256 with 10k samples. * We directly borrowed the results reported by AutoDiffusion [26] , and AutoDiffusion requires additional search costs. * We directly borrowed the results reported by AutoDiffusion [26] , and AutoDiffusion requires additional search costs. "\" represents missing data in the original paper.Step In this section, we provide the specific FID scores for pre-trained DPMs in the conditional, classifier guidance paradigm on the ImageNet 64x64 [32] and ImageNet 256x256 datasets [32] , as shown in Table 5 . We now describe the experimental setup in detail. For the pre-trained models, we used the ADM-G [5] provided 64x64_diffusion.pt and 64x64_classifier.pt for the ImageNet 64x64 dataset, and 256x256_diffusion.pt and 256x256_classifier.pt for the ImageNet 256x256 dataset. All experiments were conducted with uniform time steps and used DDIM as the baseline [20] .We implemented the second-order and third-order methods from DPM-Solver [21] for comparison and explored the method hyperparameter provided by DPM-Solver for both singlestep (corresponding to "Single" in Table 5 ) and multistep (corresponding to "Multi" in Table 5 ). Additionally, we implemented the best-performing method from DPM-Solver++ [22] , multi-step DPM-Solver++(2M), as a comparative measure. Furthermore, we also introduced the superior-performing AutoDiffusion In this section, we supplemented the specific FID↓ scores for the Stable-Diffusion [9] (conditional, classifier-free guidance paradigm) setting with a guidance scale (s) of 7.5 and 1.5. Specifically, for the pre-trained model, we conducted experiments using the sd-v1-4.ckpt checkpoint provided by Stable-Diffusion. All experiments used the MS-COCO2014 [31] validation set to calculate FID↓ scores, with uniform time steps. PFDiff employs the DDIM [20] method as the baseline. Initially, under the recommended s = 7.5 configuration by Stable-Diffusion, we implemented DPM-Solver-2 and -3 as comparative methods, and conducted searches for the method hyperparameters provided by DPM-Solver as singlestep (corresponding to "Single" in Table 6 ) and multistep (corresponding to "Multi" in Table 6 ). Additionally, we introduced previous state-of-the-art training-free methods, including DPM-Solver++(2M) [22] , UniPC [29] , and DPM-Solver-v3(2M) [27] for comparison.The experimental results are shown in Table 6 . † We borrow the results reported in DPM-Solver-v3 setting, demonstrated consistent performance improvements as seen in the s = 7.5 setting, as shown in Table 6 .Table 6 : Sample quality measured by FID↓ on the validation set of MS-COCO2014 [31] using Stable-Diffusion model [9] with guidance scales of 7.5 and 1.5, varying the number of function evaluations (NFE). Evaluated on 10k samples. † We borrow the results reported in DPM-Solver-v3 [27] directly.Step In this section, we extensively investigate the impact of the hyperparameters k and l on the performance of the PFDiff algorithm, supplemented by a series of ablation experiments regarding their configurations and outcomes. Specifically, we first conducted experiments on the CIFAR10 dataset Table 7 : Ablation of the impact of k and l on PFDiff in CIFAR10 [42] , ImageNet 64x64 and MS-COCO2014 using DDPM [2] , ScoreSDE [4] , ADM-G [5] and Stable-Diffusion [9] Particularly when k = 1 is fixed, PFDiff-1 significantly improves the sampling performance of the DDIM baseline within the range of 8∼20 NFE. For practical applications requiring higher sampling quality at fewer NFE, optimal combinations of k and l can be identified by fixing NFE and sampling a small number of samples for visual analysis or computing specific metrics, without significantly increasing the computational burden. However, as discussed in Sec. 5, although the experimental results presented in Table 7 demonstrate the excellent performance of the combinations of k and l under various pre-trained DPMs and NFE settings, no universally optimal configuration exists. This finding somewhat limits the generality of the proposed PFDiff algorithm and sets objectives for our future research.To further investigate the impact of gradient guidance from the past or future on the rapid updating of current intermediate states, this section supplements experimental results and analysis using first-order and higher-order ODE solvers as baselines. Specifically, as described in Sec. 3.3, PFDiff uses a first-order ODE solver as a baseline, where future gradient guidance corrects sampling errors, with detailed proofs provided in Appendix B.2. Hence, using the first-order ODE solver DDIM [20] as a baseline, we removed past and future gradients separately and employed quadratic time steps. This was based on the pre-trained model from DDPM [2] on the CIFAR10 [42] dataset, evaluating the FID↓ metric by changing the number of function evaluations (NFE). For higher-order ODE solvers, the solving process implicitly utilizes future gradients as mentioned in Sec. 3.5, and the additional explicit introduction of future gradients increases sampling error. Therefore, when using higher-order ODE solvers as a baseline, PFDiff accelerates the sampling process using only past gradients. Specifically, for higher-order ODE solvers, we selected DPM-Solver-2 and -3 [21] as the baseline, also employing quadratic time steps, and based on the ScoreSDE [4] pre-trained model on CIFAR10 [42] . Only the future gradients were removed to validate the effectiveness of the PFDiff algorithm by changing the NFE and evaluating the FID↓ metric. As shown in Table 8 , the experimental results indicate that using the first-order ODE solver DDIM as a baseline, employing only past gradients (similar to DeepCache [28] ), or only future gradients, only slightly improves the baseline's sampling performance; however, combining both significantly enhances the baseline's sampling performance. Meanwhile, using higher-order ODE solvers DPM-Solver-2 and -3 as the baseline, because the algorithm inherently contains future gradients, continuing to explicitly introduce future gradients increases the overall error. Therefore, using only past gradients (PFDiff) significantly improves the baseline's sampling efficiency, especially under fewer steps (NFE<10), where PFDiff markedly ameliorates the non-convergence issues of the higher-order ODE solvers.Table 8 : Ablation of the impact of the past and future gradients on PFDiff, using different orders of ODE Solver as the baseline, in CIFAR10 [42] using DDPM [2] and ScoreSDE [4] To evaluate the effectiveness of the PFDiff algorithm and the widely used Fréchet Inception Distance (FID↓) metric [40] in the sampling process of Diffusion Probabilistic Models (DPMs), we have also incorporated the Inception Score (IS↑) metric [41] for both unconditional and conditional pre-trained DPMs. Specifically, for the unconditional discrete-time pre-trained DPMs DDPM [2] , we maintained the experimental configurations described in Table 2 of Appendix D.2, and added IS scores for the CIFAR10 dataset [42] . For the unconditional continuous-time pre-trained DPMs ScoreSDE [4] ,the experimental configurations are consistent with Table 4 in Appendix D.3, and IS scores for the CIFAR10 dataset were also added. For the conditional classifier guidance paradigm of pre-trained DPMs ADM-G [5] , the experimental setup aligned with Table 5 in Appendix D.4, including IS scores for the ImageNet 64x64 and ImageNet 256x256 datasets [32] . Considering that the computation of IS scores relies on features extracted using InceptionV3 pre-trained on the ImageNet dataset, calculating IS scores for non-ImageNet datasets was not feasible, hence no IS scores were provided for the classifier-free guidance paradigm of Stable-Diffusion [9] . The experimental results are presented in Table 9 . A comparison between the FID↓ metrics in Tables 2, 4 , and 5 and the IS↑ metrics in Table 9 shows that both IS and FID metrics exhibit similar trends under the same experimental settings, i.e., as the number of function evaluations (NFE) changes, lower FID scores correspond to higher IS scores. Further, Figs. 1a and 1b , along with the visualization experiments in Appendix D.8, demonstrate that lower FID scores and higher IS scores correlate with higher image quality and richer To demonstrate the effectiveness of PFDiff, we present the visual sampling results on the CIFAR10 [42] , CelebA 64x64 [43] , LSUN-bedroom 256x256 [44] , LSUN-church 256x256 [44] , ImageNet 64x64 [32] , ImageNet 256x256 [32] , and MS-COCO2014 [31] datasets in Figs. 5-10. These results illustrate that PFDiff, using different orders of ODE solvers as a baseline, is capable of generating samples of higher quality and richer detail on both unconditional and conditional pre-trained Diffusion Probabilistic Models (DPMs). Figure 9 : Random samples by DDIM [20] and PFDiff (baseline: DDIM) with 4 and 8 number of function evaluations (NFE), using the same random seed, uniform time steps, and pre-trained Guided-Diffusion [5] on ImageNet 256x256 [32] with a guidance scale of 2.0.Justification: We provide code in the supplemental materials, which will be made available on the GitHub platform.Guidelines:• The answer NA means that paper does not include experiments requiring code.• Please see the NeurIPS code and data submission guidelines (https://nips.cc/ public/guides/CodeSubmissionPolicy) for more details.• While we encourage the release of code and data, we understand that this might not be possible, so "No" is an acceptable answer. Papers cannot be rejected simply for not including code, unless this is central to the contribution (e.g., for a new open-source benchmark).• The instructions should contain the exact command and environment needed to run to reproduce the results. See the NeurIPS code and data submission guidelines (https://nips.cc/public/guides/CodeSubmissionPolicy) for more details.• The authors should provide instructions on data access and preparation, including how to access the raw data, preprocessed data, intermediate data, and generated data, etc.• The authors should provide scripts to reproduce all experimental results for the new proposed method and baselines. If only a subset of experiments are reproducible, they should state which ones are omitted from the script and why.• At submission time, to preserve anonymity, the authors should release anonymized versions (if applicable).• Providing as much information as possible in supplemental material (appended to the paper) is recommended, but including URLs to data and code is permitted.Justification: Our work is conducted with the NeurIPS code of ethics.Guidelines:• The answer NA means that the authors have not reviewed the NeurIPS Code of Ethics.• If the authors answer No, they should explain the special circumstances that require a deviation from the Code of Ethics.• The authors should make sure to preserve anonymity (e.g., if there is a special consideration due to laws or regulations in their jurisdiction). 


Decision: This work addressed the common discretization errors of fast ODE solvers in accelerating the sampling of Diffusion Probabilistic Models. These errors typically  grow as the number of function evaluations decreases, and the  proposed method, namely PFDiff, is presented as a training-free with an  orthogonal timestep-skipping strategy, which effectively reduces the number of function evaluations relative to existing fast ODE solvers. The strategy is specifically rooted in exploiting two key observations: a strong closeness of  the model's outputs at relatively close time step size during the denoising process of existing ODE solvers, and a resemblance between the denoising process and SGD. The proposed PFDiff uses a  gradient replacement from past time steps and foresight updates to  rapidly update intermediate states, and  reduce unnecessary NFE which affect the discretization errors in first order ODEs. A solid experimental evaluation validates the claimed speed improvement , demonstrating applicability across a number pre-trained Diffusion Probabilistic Models, to surpass methods for   conditional DPMs and  previous state-of-the-art training-free methods. 
Five reviews were  submitted with 2 W. Accpt, 1 BL Accpt. and 2 BL. Rejects. The authors  have been very meticulous in their rebuttal with additional experiments and comparisons, and point to point rebuttal: In particular, Reviewer 7w4h claimed no novelty and provided  references, the authors have addressed each and every reference showing that the reviewer was either confused or not totally familiar with the references provided. They have also provided answers to all the questions to the Reviwer's satisfaction  and who had promised to revisit the grade. Well he did not.  The other BL Reject by Reviwer FcEN has raised a number of points most of which turned out to be the result of his/her  confusion. Specifically, FcEN asking for further explanations why using past gradients was helpful in guiding the new faster sampling approach, which of course was not the case and the authors politely pointed that out. Having requested an evaluation and comparison based on ImageNet 256, with an openness to review his grading. All that was provided by the authors with no follow up with several reminders sent by the AC.  
The AC is reluctantly ranking this paper as reject in light of the two "weak" BRs, The AC quickly read the paper and read the additional edits proposed by the authors, and estimates that the authors did their part, and while the paper is perhaps not award winning, it is acceptable as a Poster.

Reviews: review: paper_summary: paper_summary: This work proposes a training-free time step-skipping method that can be used with existing ODE solvers for reduced NFE. The method was motivated by two observations: 1) a significant similarity in the model's outputs at time step size during the denoising process and 2) a high resemblance between the denoising process and SGD. The proposed method employed gradient replacement from past time steps and rapidly updated intermediate states inspired by Nesterov momentum. The proposed method yielded promising results.
strength_weakness: strength_weakness: - There have been a lot of prior works on accelerating diffusion models for sampling. While this manuscript cited many, it still missed important prior works - some of them look quite similar to the proposed method. Thus, the novelty of the proposed method is unclear in the current form of this manuscript. For example, using Nesterov acceleration for fast diffusion models is not really new (e.g., R Li et al., Hessian-Free High-Resolution Nesterov Acceleration For Sampling, ICML 2022). Eq (15) of this work can be seen as a special case of the following prior works such as [R1], DeepCache [28] (using past), [R3] (using three moments or future) or [R2] (using all). Some recent work like [R4] even used partial caching instead of using the whole results. A more theoretically grounded work on using Nesterov momentum for sampling can be found in [R5]. 
[R1] M Xia et al., Towards More Accurate Diffusion Model Acceleration with A Timestep Tuner, CVPR 2024.
[R2] A Pokle et al., Deep Equilibrium Approaches to Diffusion Models, NeurIPS 2022.
[R3] H Guo et al., Gaussian Mixture Solvers for Diffusion Models, NeurIPS 2023.
[R4] F Wimbauer et al., Cache Me if You Can: Accelerating Diffusion Models through Block Caching, CVPR 2023.
[R5] R Li et al., Hessian-Free High-Resolution Nesterov Acceleration For Sampling, ICML 2022.
- A number of acceleration works for diffusion models also investigated the feasibility of the parallel computation. Will the proposed method be parallelized for computation? 
- It is unclear if the proposed method was compared with other methods in terms of computation. Will 1 NFE of the proposed method take the same computation time as 1 NFE of other methods since the proposed method contains multiple evaluations of the neural network as in Eq. (15).
- The notation and explanation are quite confusing, so it is not easy to understand the whole idea as well as the algorithm itself.
questions: questions: See the weaknesses.
limitations: limitations: N/A
paper_summary: paper_summary: The paper proposes a new training-free acceleration method for the inference of diffusion probabilistic models. The key components of the presented time-skipping strategy are the use of past and future gradients to eliminate redundant neural function evaluations (NFE). The proposed method is shown effective compared to other training-free acceleration methods, leading to solid performance improvements especially for ODE solvers with less than 10 NFEs.
strength_weakness: strength_weakness: *    The performance gap compared to training-based methods is still apparent, especially considering the latest distillation techniques resulting in one-step models.
*    The mathematical notations are a bit hard to follow up. I would advise the authors to add a schematic clarifying for a given setting of hyperparameters, which timepoints are being evaluated and which are being skipped.
*    The optimal setting of hyperparameters *k,l* is model/dataset dependent and it is not clear apriori how to set these. Therefore, this requires empirical experimentation which makes it time-consuming to get optimal performance when using the method out-of-the-box.
questions: questions: *    How does your samples fare in terms of diversity compared to the original model with enough NFE? In other words, does the faster sampling somehow come at the cost of reduced diversity?
limitations: limitations: The authors were upfront about the limitations of their method.
paper_summary: paper_summary: The paper proposes PFDiff, a training-free approach for accelerating diffusion models. Motivated by the high similarity of the diffusion network outputs at adjacent timesteps on the sampling trajectory, PFDiff utilizes past and future information for sampling with time-skipping, 
and decreases the number of function evaluations (NFEs) significantly. Experiments on various settings show significant acceleration, especially in the low NFE regime.
strength_weakness: strength_weakness: - The highly concise writing and complex notations might be a bit confusing. Additional illustrations for certain local algorithm procedures can be helpful for understanding the overall idea.
- There are fundamental mistakes in the writing. Eqn. (8) (9) are represented as Euler discretizations of the original PF-ODE. However, both DDIM and the series of DPM-Solvers rely on exponential integrators to transform the PF-ODE into other forms, so that the linear term $x_t$ is cancelled. Though this does not mean the method is wrong, such simplified writing can be misleading. The authors are obligated to correct this, or I will be forced to reject this paper.
- It will be more convincing to include experiments on EDM, the SOTA diffusion model on CIFAR-10 and ImageNet 64x64.
questions: questions: - Are there any insights why PFDiff is more effective on first-order solvers (DDIM+PFDiff even outperforms high-order solvers)?
limitations: limitations: Yes.
paper_summary: paper_summary: This paper proposes PFDiff, a fast training-free sampler for diffusion models. PFDiff updates the current state with both the past score network evaluation and the future score network evaluation. It can achieve good sample quality with less than 10 NFE. The authors showcase the effectiveness of PFDiff on various pre-trained diffusion models.
strength_weakness: strength_weakness: 1. Flawed justification for future gradient: The authors' claim that using future gradient information is better than using current gradient information is based on the mean value theorem (lines 164-167 and Appendix B.2). However, this theorem only guarantees the existence of an optimal point within an interval, not its specific location. Therefore, the mean value theorem itself doesn't justify the preference for future gradients. 
2. Missing justification for the approximation: While the authors claim that their approximation is better, there is no theoretical justification for it. The proof in Appendix B.2 assumes that the optimal point is already known, which is not informative. The manuscript will benefit from a further approximation error analysis. 
3. Expensive and case-specific tuning: the proposed method essentially defines a set of candidate points and searches for the optimal point by tuning parameters $k$ and $l$. This tuning process can be computationally expensive and needs to be done for each specific case, limiting its practicality.
questions: questions: Line 163-164: the mean value theorem always holds for any continuous differentiable function. What does it even mean for the mean value theorem to hold "approximately"? Also, it's unclear how the observations mentioned in lines 161-163 indicate the mean value theorem.
limitations: limitations: The algorithm performance depends heavily on parameters $k$ and $l$ as shown in Table 7.  Optimal values for $k$ and $l$ vary based on the pre-trained model and the number of function evaluations. This necessitates extensive parameter tuning when applying the proposed method in practice.
paper_summary: paper_summary: To accelerate the sampling speed in diffusion models, this paper proposes a training-free denoising method, dubbed PFDiff. 
Concretely, PFDiff employs the gradient from past time steps to update intermediate states, aiming to reduce unnecessary NFEs while correcting for discretization errors.
In this manner, PFDiff enables to improve classic samplers without any training computation.
Importantly, experimental results demonstrate the effectiveness of the proposed PFDiff.
strength_weakness: strength_weakness: 1. In my humble opinion, the theoretical analysis part is naive. Can you provide more explanation about why previous gradients is helpful to guide current sampling direction? Since different noise levels correspond to different gradients, is there any harm in denoising images with the proposed method?

2. Many works investigate using previous gradients to improve sampling speed, so the contribution is limited.
questions: questions: 1. Is there any memory overhead for saving the gradients?

2. Can you present more experiments on ImageNet 256? Including conditional and unconditional.

3. How about the comparison with DEIS sampler?

4. Can the proposed method connect with distillation models?

5. If more metrics are tested, such as recall and precision, it would be better.
limitations: limitations: Please see in Weaknesses and Questions. If all of my concerns are addressed, I will improve my score.

Comments: comment: Thank you again for your response and for raising the rating!

Our proof shows that future gradients can reduce discretization errors caused by neglecting higher-order terms, but the errors from remainder terms still exist. We would like to emphasize that the current gradients **neglect, rather than effectively address,** the errors from remainder terms and higher-order derivative terms. In a few cases, The interaction of the two discretization errors cancels each other out, leading to the exact result. In practical applications of DPMs sampling, such situations generally do not occur due to the complexity of neural network functions. 

Thank you again for the counterexamples you provided. We will further elaborate such special cases in the proof process of our manuscript. We are very pleased to discuss if you have any further questions!
Thanks for the response. I've raised my score to 5. 

Regarding your latest response, the approximation based on future time might be more accurate for neural networks in practice. But I think the proof is not rigorous. One can always construct a higher-dimensional and higher-order function example like the one constructed above. We can add as many higher-order terms as we want, such as the current gradient is optimal. For example, $s_{\theta}(t)=1+ \sum_{k=1}^{K} (2kx^{2k-1}-(2k+1)x^{2k})$. The approximation based on the current time step is always exact, but PFDiff with finite order is not always exact depending on the choice of $r$.
Thank you very much for your response! We are glad that the rebuttal has addressed almost all of your concerns. Let us know if further discussion is needed.
We sincerely appreciate your consideration of raising the rating. In terms of the new issue in the analysis, first, we want to emphasize that utilizing current gradients introduces two discretization errors: **neglecting remainder terms and neglecting higher-order derivative terms**. The future gradients aim to reduce the impact of neglecting higher-order derivative terms, rather than remainder term $\int_{t_{i-1}}^{t_i} R_n(t)dt$.

We provide a more detailed response here:

We carefully checked the example you provided, and found it is just a "**coincidence**" that occurs in low-dimensional, low-order functions. In our proof process, there is a step involving approximation:

$\int_{t_{i-1}}^{t_i}\left[\sum_{n=0}^{\infty} \frac{s_\theta^{(n)}\left(x_{r}, r \right)}{n!}(t-r)^n+R_n(t)\right] dt \approx
\frac{1}{(n+1)!} \sum_{n=0}^{\infty} s_\theta^{(n)}\left(x_{r}, r\right)\left[\left(t_i-r\right)^{n+1}-\left(t_{i-1}-r\right)^{n+1}\right].$

Whether using current or future gradients, the impact of the remainder term $\int_{t_{i-1}}^{t_i} R_n(t)dt$ is neglected. Further, in the discretization of integrals for a first-order ODE solver, higher-order derivative terms (containing $s_\theta^{(n)}\left(x_{r}, r \right), n \ge 1$) are directly neglected, thus our focus is on analyzing the impact of higher-order derivative terms. In our previous response, we demonstrated using future gradients can reduce the discretization errors caused by neglecting higher-order terms (since the coefficients' absolute values are smaller).

Regarding the further analysis of the example $s_{\theta }(x_{t}, t)=1+4t-6t^{2}$. **Firstly**, based on the future time point $t=\frac{2}{3}$ (i.e., PFDiff-2_2), $\int_{0}^{1} s_{\theta }(x_{t}, t) dt=s_{\theta }(\frac{2}{3} )\int_{0}^{1} dt=1$ is also exact. **Secondly**, when $n\ge 3$, the higher-order derivative terms in the example are all zero. However, in the sampling process of DPMs, this corresponds to the derivatives of high-dimensional neural networks $s_\theta^{(n)}\left(x_{r}, r\right)$, which generally are not zero. Therefore, the future gradients are dedicated to complex neural network functions in practical applications rather than low-dimensional, low-order functions. **Lastly**, the "correct" solution is merely a "coincidence" in the example, as the current gradients neglect rather than effectively address the errors from remainder terms and higher-order derivative terms. The interaction of two errors leads to the "correct" result.

We have updated the proof process in the paper based on the above points. We are very pleased to discuss If you have further questions. Thank you again for your response!
I acknowledge that I have read the rebuttal. The rebuttal addressed almost all of my concerns. I will finalize my score after the discussion with other reviewers.
Thanks for the detailed rebuttal. It addressed most of my concerns. I'm leaning toward borderline acceptance now. I appreciate the effort in approximation analysis. However, I think there is an issue in the analysis.
- The new approximation error analysis is based on the absolute values of the coefficients of higher-order derivatives. It does not necessarily imply that the approximation error $\int R_n(t)dt$ is smaller. In fact, using the current gradient could be optimal. For example, consider the case where $s_{\theta}(x_t, t)=1+4t-6t^2$ and $t_i=1, t_{i-1}=0$. The approximation error based on the current point $\int_0^1s_{\theta}(x_t, t)dt=s_{\theta}(0)\int_0^1dt=1$ is exact.
We are glad to know that your questions have been addressed! We greatly appreciate your valuable suggestions, which help to improve our paper. Thanks again for your recognition and for maintaining a positive score!
Thanks. The authors have addressed my questions, and I'm happy to retain my positive score.
We are so glad to hear that your concerns were well addressed! Thanks again for recognizing our work!
Thank you for the detailed responses. My concerns are well addressed, and I keep the score of leaning towards acceptance in the current reviewing stage.