Response to Weakness:

• [assumption 2]: We acknowledge that Assumption 2 may appear stringent compared to previous works. The reason we introduce this assumption is tied to the specific algorithm we are analyzing: SDE-DPM-2. At each update step, this algorithm requires the first derivative of the score function, as described above in equation (6). This naturally leads to a requirement for a close estimation of the score’s first derivative. Furthermore, as discussed under Assumption 2, Meng et al. (2021) empirically demonstrated that the first derivative of the score can be learned effectively under Gaussian mixture models. Therefore, despite the stringency of Assumption 2, our approach remains empirically reasonable, supported by experimental results.

• [assumptions 3 and 4]: We would like to clarify that in Appendix B, we specifically verify that Assumptions 3 and 4 hold in the case of a Gaussian mixture distribution. For other types of distributions, such as those supported on a low-dimensional manifold, it is currently unclear whether these assumptions hold as $t$ approaches 0. This is something we hope to explore further in future work.

• [inconsistent]: We would like to clarify that in equation (6), the first-order derivative represents the total derivative of the score function with respect to $t$. As explained below equation (12), we use the chain rule to express total derivatives in terms of partial derivatives. Regarding the statements in Lines 201 and 401, both refer to the same concept: the key difference between EI and SDE-DPM-2 lies in how the score function is approximated within the update scheme at each time step. We have revised the manuscript to ensure that this distinction is more clearly explained and standardize the notation to avoid any confusion.

• [rk-2]: We apologize for any misunderstanding. Corollary 3.3 indeed only indicates that, under our current analytical method, Runge-Kutta does not demonstrate higher efficiency than SDE-DPM-2.

  Regarding the experimental results: we have redrawn the comparison of KL divergence against the number of discretization steps on a log scale. While RK-2 achieves a slightly lower KL divergence than SDE-DPM, the rate of decrease for both methods is similar and slower than that of SDE-DPM-2. This observation aligns with our theoretical findings. Moreover, we acknowledge that there is a gap between existing theoretical analyses and empirical observations, as experiments show that the KL divergence for each method decreases more rapidly than the theoretical bounds. This suggests potential directions for future research to further improve the convergence rate of these methods.

• [VE]: In VP-type diffusion models, the initial error term decays exponentially fast, thus contributing minimally to overall convergence. On the other hand, as shown in equation (17), the VE-type diffusion model has an initial error term that decays at a polynomial rate, specifically $1/T$, which becomes the dominant factor relative to the discretization error term of $1/T^2$. Therefore, in the case of the VE-type diffusion model, when the initial error is overlooked (i.e., assuming the reverse process distribution $q_0$ at time $0$ matches the forward process distribution $p_T$ at time $T$), we can directly compare the algorithm’s discretization error. This key point has been emphasized in the introduction of the revised manuscript.

• [Proposition 4.2]: We have provided a more detailed explanation of the proof of Proposition 4.2 in Appendix C.1 of the revised manuscript to ensure that the argument is clear and self-contained.

Response to Question:

• [ better convergence rate for rk-2]: Our analysis focuses on the stochastic RK-2 algorithm described in Equations (9-10). The claim regarding the lower efficiency of RK-2 is primarily based on its direct discretization of the linear component of the drift term in the SDE. If the error bound for this discretization could be further reduced, it might be possible to achieve a better convergence rate for Runge-Kutta methods.

Response to Weakness:

• [assumptions]: As discussed under Assumption 2, experimental work by Meng et al. (2021) designed a loss function that approximates both the score and its derivative, demonstrating that solvers using the estimated score and score derivative can achieve fast convergence. This suggests that the score's derivative can be well-approximated. Additionally, experiments by Meng et al. (2021) on Gaussian mixture models have shown that this method can effectively learn more accurate the score derivatives.

  Furthermore, while some studies have theoretically proven that $L^2$ accuracy in score estimation is achievable when the target distribution is a simple Gaussian mixture (e.g., symmetric bimodal Gaussian mixtures), the theoretical understanding of the learnability of scores and second-order scores for more complex distributions remains limited. This is an area that merits further exploration.

• [ figure 1]: We have redrawn Figure 1 with a logarithmic scale to better illustrate the theoretical convergence rates of the different methods. The revised figure clearly shows that SDE-DPM-2 achieves a faster convergence rate than SDE-DPM and RK-2. Moreover, we acknowledge a gap between existing theoretical analyses and empirical observations, as experiments indicate that the KL divergence for each method decreases more rapidly than predicted by theoretical bounds. This discrepancy suggests potential directions for future research aimed at improving the convergence rates of these methods.

• [improve dimension dependence]: Thank you for your insightful comment. We acknowledge that [2] achieves an improvement to sample complexity scaling from $O(d^{3/2})$ to $O(d)$ by leveraging stochastic localization, which reduces the dependency of $\mathbb{E}\|\nabla^2 \log q_t\|^2_F$ on the data dimension $d$ . However, the solver we analyze in this work, SDE-DPM-2, is fundamentally different from [1] and [2]: SDE-DPM-2 relies on the first derivative of the score function during sampling.

  In our analysis (as detailed in Appendix C.1), we employ a Taylor expansion to bound discretization error. Even when assuming a constant second derivative of the score, this approach introduces a dependency of $O(d^{3/2})$. Therefore, at present, it is challenging to reduce the sample complexity's dimension dependency for SDE-DPM-2 to $O(d)$.

• [literature review]: We sincerely appreciate the reviewer’s suggestions and for highlighting several important studies related to diffusion models, which will enhance our theoretical understanding. Our primary focus is on the sampling complexity of diffusion model solvers, and we find that [4, 5, 9] are more closely related to our work. It is worth mentioning that [4] was published after our submission and introduces a novel stochastic RK-2 method. The assumptions in [5] may be a little stringent, requiring the target distribution to have a strongly convex potential, be four-times differentiable, and have Lipschitz continuous first three derivatives. Study [9] discusses parallel sampling methods, but each sampling step involves multiple evaluations of the score function, which may be computationally demanding. We have included discussions on these studies in our revised manuscript. The other studies suggested by the reviewer are mainly empirical or focus on ODE-based samplers, which we believe may not require extensive discussion in the related work section.

Response to Questions:

• [assumptions 3 and 4]: Unlike in [2], where the discretization error's derivative is directly bounded using expectations of $|\nabla \log p_t(x_t)|^2$ and $|\nabla^2 \log p_t(x_t)|^2$, our analysis focuses on the convergence of the SDE-DPM-2 method. This method introduces an additional term—the first derivative of the score function—into the discretization error. The presence of this extra term complicates the direct bounding of the error’s derivative as done in [2].

  To address this, we decompose the discretization error into two components: temporal and spatial. We bound these components independently using Assumptions 3 (for the partial derivative with respect to $t$) and Assumption 4 (for the partial derivative with respect to $x_t$). We acknowledge that directly bounding the derivative of the discretization error in the same way as [2] may not be feasible due to the additional term in the SDE-DPM-2 discretization error, and we plan to explore this further in future work.

• [if assumption 2 is replaced]: Due to the nature of our proof technique, Assumption 2 is necessary for the theoretical validation. However, it is worth noting that in our experimental comparisons between SDE-DPM-2 and SDE-DPM, we used commonly available pre-trained models (e.g., those trained for DDPM). Our results indicated that SDE-DPM-2 outperformed SDE-DPM even without specifically training to satisfy Assumption 2.

  Similar to the common training objective of Denoising Score Matching (DSM), Assumption 2 can be validated in cases where the ground truth score and its derivative are explicitly computable, such as for Gaussian mixture distributions. This validation was demonstrated in the experiments by Meng et al. (2021). However, for more general distributions, neither Assumption 2 nor the assumption from [1] can be empirically validated during training at present.

Chenlin Meng et al. "Estimating High Order Gradients of the Data Distribution by Denoising." (2021).

Response to Weakness

• Wu et al. (2024) introduced a novel stochastic RK-2 method that differs from the stochastic RK-2 algorithm analyzed in our paper (Equations (9-10)). Our claim regarding the lower efficiency of RK-2 is based on its direct discretization of the linear component of the drift term in the SDE, which results in a higher error compared to the SDE-DPM-2 method. In contrast, Wu et al. (2024) developed a new stochastic RK-2 method that effectively incorporates the linear term of the drift, making it more efficient than the RK-2 method discussed in our work. We have also demonstrated that the Rk-2 analyzed in our paper is less efficient than SDE-DPM-2, as shown in Figure 1. We have clarified this distinction in the revised manuscript to avoid any potential confusion.

• As discussed in Appendix A, our focus is on the existing SDE-DPM-2 method [Lu et al., 2022] and its theoretical significance. While Li et al. (2024) designed a different sampler and provided sampling complexities of $O(d^3/\sqrt{\epsilon})$ to achieve $\tilde{O}(\epsilon)$ in KL divergence and $O(d^3/\epsilon)$ to achieve $\tilde{O}(\epsilon)$ in TV, our contribution lies in analyzing the SDE-DPM-2 method. SDE-DPM-2 has been experimentally validated for strong performance but lacked comprehensive theoretical analysis—a gap that our work fills. Our analysis demonstrates that SDE-DPM-2 achieves comparable complexities, specifically $O(d^{1.5}/\sqrt{\epsilon})$ in KL divergence and $O(d^{1.5}/\epsilon)$ in TV, matching Li et al. (2024) in terms of $\epsilon$ dependence but with improved dependence on the data dimension $d$.

  Furthermore, the SDE-DPM-2 method is empirically more efficient: while both algorithms in Wu et al. (2024) and Li et al. (2024) achieve FID scores around 100 (Figure 2 in Wu et al., 2024), SDE-DPM-2 achieves an FID score of around 18 when sampling from the CIFAR-10 dataset with 20 steps. We have discussed this comparison in more detail in the revised manuscript to provide a clearer perspective on the theoretical contributions of our work.

• Thank you for pointing out the comparison with these works. As we noted in our introduction, the first referenced work primarily provides convergence analysis for an ODE-based solver, whereas our focus is on SDE-based solvers. The second referenced work introduces a novel algorithm, RTK-ULD, and provides a sampling complexity of $\tilde{O}(L^2 d^{1/2} / \epsilon)$ to achieve $\tilde{O}(\epsilon)$ accuracy in total variation distance. It is important to emphasize that our analysis centers on existing algorithms, specifically SDE-DPM-2. While the second work reports a better dependence on the dimension $d$, our results maintain consistency regarding the $\epsilon$ dependence. Notably, the second paper includes dependencies on the Lipschitz constant, which is often associated with the data dimension, suggesting that their dimensional dependence might still be incomplete.

Response to Weakness:

• [computational cost]: The primary factor influencing the computational cost of SGMs is the evaluation of the estimated score function, $s_theta$, typically represented by a large neural network. SDE-DPM-Solver++(2M) efficiently updates the derivative of the score function by reusing previously stored evaluations, thus matching the computational cost of a first-order method. In contrast, as shown in Equations (9-10), the RK-2 method requires evaluating the score function twice per time step, resulting in significantly higher computational cost. Consequently, both SDE-DPM-Solver++(2M) and EI have similar computational costs, which are notably lower than that of RK-2. We have clarified this comparison in the revised manuscript.
• [ comparision with Chen 2023]: Chen et al. (2023) conducted a convergence analysis of the Exponential Integrator/SDE-DPM method. In contrast, our work focuses on the convergence analysis of the SDE-DPM-2 method, a second-order extension of the SDE-DPM method. We provide the first convergence analysis of SDE-DPM-2, demonstrating that it achieves a faster convergence rate compared to the original SDE-DPM method. We have compared the empirical performance of SDE-DPM-2 with SDE-DPM as demonstrated in Table 1 and Figure 1. We revised the manuscript to emphasize this distinction and clearly delineate our contributions from those of Chen et al. (2023).
• [assumptions]: We have verified that under Gaussian mixture models, the score function does not exhibit highly singular behavior near $t=0$, confirming that the conditions outlined in Assumptions 3 and 4 hold. Since Gaussian mixtures are universal approximators for smooth probability density functions, these assumptions can be extended to a wide range of distributions that are effectively approximated by Gaussian mixtures. To clarify this, we have included a more detailed discussion in the revised manuscript. However, we acknowledge that the potential singular behavior of the score function for other types of distributions near $t=0$ remains less explored and presents an interesting direction for future research.

Response to Questions:

• [$\nabla^3$]: We use the operator $\nabla^3$ in Assumption 4 to bound the discretization error term in Appendix C.1, the proof of Lemma 4.3. We have clarified this point in the revised manuscript to ensure that the use of the operator $\nabla^3$ is explicitly stated and justified.
• [error bar]: We run each method multiple times and give the result in Table 1 of the revised manuscript, which further substantiates the improved experimental performance of SDE-DPM-2 over SDE-DPM.

Response to Weakness:

• [Structure]: We agree that the structure you suggested is indeed more logical and easier to follow. We have revised the manuscript accordingly to reflect this improved structure.

• [Acronyms for Discretization Methods]: We appreciate your feedback and will make sure to define all acronyms for the discretization methods used in the paper to avoid any confusion for the readers.

• [Notation for $x_k$]: We will clearly point out that $x_k$ refers to $x_{t_k}$ in the manuscript to avoid any confusion.

• [Partial Derivatives with Respect to $x_{t_k}$]: We have clarified that the partial derivatives with respect to $x_{t_k}$ refer to the Jacobian matrix and explicitly defined this notation for clarity in the revised manuscript.

• [Table 1 and Empirical Confirmation]: While Lu et al. (2022b) demonstrated that SDE-DPM-2 generates better samples than SDE-DPM through image generation examples, our Table 1 provides additional support for this claim by comparing specific metrics (e.g., FID scores). This further substantiates the improved experimental performance of SDE-DPM-2 over SDE-DPM, offering a more detailed empirical validation of its effectiveness. We run each method multiple times and give the result in table 1.

Response to Questions:

• [VP and VE]: The main distinction between the two types of forward processes lies in the behavior of the initial error term. In the VP-type diffusion model, the initial error term decays exponentially fast, contributing minimally to the overall convergence. In contrast, for the VE-type diffusion model, as outlined in equation (17), the initial error term decays at a polynomial rate, specifically $1/T$, which becomes the dominant factor compared to the discretization error term $1/T^2$. Consequently, we conduct separate convergence analyses for the VP- and VE-type diffusion models.
• [Figure 1]: We have incorporated the empirical and theoretical convergence results for RK-2 into Figure 1, employing a logarithmic scale to depict the theoretical bounds of each method. Our analysis demonstrates that, while RK-2 and SDE-DPM show comparable empirical performance, both are less efficient than SDE-DPM-2. Moreover, we identified a gap between existing theoretical results and empirical observations, as the KL divergence for each method decreases more rapidly than the theoretical bounds. This may suggest potential directions for future research to further improve the convergence rate of these methods.
• [equation 11 and 13]: Equations (11) and (13) are not identical. The drift terms in the two equations differ: Equation (11) contains $\hat{x}_{t_K}$ in the drift term, which is a constant within each time interval, whereas Equation (13) includes $x_t$ in the drift term, which is time-dependent. We clarified this distinction in the revised manuscript, as it represents the key difference between the SDE-DPM-2 and RK-2 methods.