Efficient Integrators for Diffusion Generative Models

We thank the reviewer for their insight into our work and respond to each question below:

Weakness 1 : Difficult to figure out the "main" contribution

There are too many proposed integrators, and it is difficult to figure out which one is / ones are the "main" contribution of the paper. More specifically, it is difficult to see in which situation we should prefer one integrator over the other. To me, it seems like the authors are taking a very unclear stance on their "main" sampler : the authors are proposing a whole arsenal of samplers, and picking results which happened to beat some baselines results under very specific dataset+sampler+hyper-parameter combinations.

Response: We address this question in different parts

There are too many proposed integrators, and it is difficult to figure out which one is / ones are the "main" contribution of the paper.

As highlighted in Section 4 (Notation), the main samplers that we recommend using are the Conjugate-Splitting based PSLD samplers (CSPS) and Splitting based PSLD samplers (SPS). We understand that the samplers specified in Table 2 might suggest a lot of contributions. However, we would like to emphasize that Table 2 lists our ablation samplers, the discussion of which is necessary to construct our main samplers, which build on top of these samplers discussed in Table 2.

More specifically, it is difficult to see in which situation we should prefer one integrator over the other

From Figure 5, we make the following observations. For deterministic sampling, our CSPS samplers always perform better than SPS samplers. For stochastic sampling, except CIFAR-10, the CSPS integrators always perform better than the corresponding SPS samplers. Therefore, for most practical applications, the proposed CSPS samplers should be preferred.

the authors are proposing a whole arsenal of samplers, and picking results which happened to beat some baselines results under very specific dataset+sampler+hyper-parameter combinations.

We would like to point out that we do not perform any type of cherry-picking when comparing different sampler baselines and our proposed samplers. On the contrary, we perform an extensive comparison across multiple sampling budgets for different datasets in Figure 5, which clearly illustrates when our deterministic and stochastic samplers perform comparably or outperform different baselines. Moreover, our choice of datasets is standard in the diffusion model literature [1,2,3,4]. Regarding hyperparameters, we tune the value of \lambda for CSPS samplers during inference since Theorem 2 suggests that a specific value of \lambda is unlikely to work across multiple integration step sizes h (i.e., for different sampling budgets in terms of NFEs) and is thus a core requirement for CSPS samplers. Moreover, since we evaluate other baselines based on the scores provided in their respective papers, it is safe to assume that our comparison baselines have also been optimized for best performance.

For instance, what is the rationale behind choosing Bt=λI and Bt=λ1? When should we prefer DDIM-I over DDIM-II and vice versa? While there is a theorem on DDIM-I, I don't see any theorem on DDIM-II, which could clarify differences between DDIM-I and DDIM-II.

Response: We thank the reviewer for this question. The main rationale behind choosing $B_t$ as in $\lambda$-DDIM is mostly simplicity. Our current choice of $B_t =\lambda I$ assumes that the value of $\lambda$ is independent of time $t$, which might be overly simplistic since Theorem 2 suggests that at any time $t$, the eigenvalues $\tilde{\lambda}$ must be bounded for a stable conjugate integrator, therefore, implying a time-varying $\lambda$ might be a better choice. However, we make this simplistic choice because a time-independent lambda is more straightforward to tune during inference. We had the same rationale when choosing $B_t = \lambda 1$.

While there is a theorem on DDIM-I, I don't see any theorem on DDIM-II, which could clarify differences between DDIM-I and DDIM-II.

We would like to point out that for the choice of $B_t = \lambda I$, the result in theorem 2 simplifies, resulting in Corollary-1. However, similar simplifications are not apparent for other choices of $B_t = \lambda 1$, and Theorem 2 holds in general for other choices of $B_t$. We acknowledge that further theoretical investigation of different choices of $B_t$ could clarify differences between the two formulations and could be an interesting future direction. We already state this point as an interesting direction for future work in the conclusion.

We thank the reviewer for their positive comments. We believe that the frameworks proposed in this paper for fast sampling in diffusion models can be a fruitful starting point for further research in this area.

Questions:

Do the authors have an explanation why some methods outperform CSPS-D for very small numbers of NFE? (Figure 5, up-left plot).

Response: We thank the reviewer for this question. For a very small number of NFEs, baselines like DPM-Solver, and DEIS outperform our proposed solvers. To intuitively understand this observation, we highlight that we evaluate our proposed samplers on the PSLD diffusion trained with the epsilon prediction objective. As noted in [1], for augmented space models like CLD, the dynamics in the epsilon space are oscillatory in nature, which prohibits taking larger step sizes with these models. On the other hand, baseline samplers in Figure 5, like DPM-Solver and DEIS, are evaluated on VP-SDE, which does not suffer from this problem and thus allows for larger step size. To quantify this intuition, we generalize DEIS and apply it to the PSLD diffusion ODE. The results in Table 1 below show that DEIS applied to PSLD performs significantly worse in comparison to the proposed samplers in this work across all timesteps.

Have the authors measured the total integration time? If the generation speed is our main objective, then such results should be provided as well.

Response: We assume the reviewer means the integration time for computing the coefficients required for ODE integration. Given a set of sampling timepoints ${t_i}$, since coefficients are shared between all samples, we only need to perform this process once at the start of sampling. Empirically, for our largest budget of 100 time points, numerical integration for computing coefficients takes around 20s on our setup. We have added these results in Appendix B.6

How is m0 defined during the training of 'Conjugate Symplectic Euler' and 'Conjugate Velocity Verlet'? That is, what is the generated mϵ. In my opinion, the paper would be improved if corresponding training algorithms to Algorithms 2 and 3 would be added. On a related note, does the transformation along the position dimensions xt remain a diffeomorphism in the augmented setting? If not, this could have negative implications about density estimation, and mode coverage in generation.

Response: We assume that by mentioning the training of “Conjugate Symplectic Euler” and “Conjugate Velocity Verlet”, the reviewer intends the training of the PSLD diffusion model since Conjugate Symplectic Euler and Conjugate Velocity Verlet are sampling algorithms and do not require any training. We hope this clarifies any potential misunderstanding. Next, $m_0$ i.e., the auxiliary variables in PSLD, are sampled from a Gaussian distribution $N(0, MI)$ where $M$ is the mass matrix as defined in PSLD. We encourage the reviewer to check the PSLD paper [2] for more details about $m_0$ and its training algorithm.

Regarding comments about the diffeomorphism of the transformation in the position space, we would like to point out that as long as the transformation in the augmented space itself is a diffeomorphism, we can invert the dynamics at any time to reliably recover the augmented space $[x_t, m_t]$ and then discard $m_t$ to recover $x_t$. Therefore, ensuring that the transformation remains a diffeomorphism in the $x_t$ dimension becomes moot since the diffusion happens in the joint space of $x_t$ and $m_t$.

Additional Results

Table 1: Comparison between different ODE-solvers for the CIFAR-10 dataset at various sampling budgets. DEIS was generalized using the Conjugate Integrator framework and applied to PSLD reverse diffusion ODE. DEIS (q=k) represents DEIS with polynomial extrapolation (where k represents the polynomial order), which is known to improve performance for low NFEs [5]. CSPS-D (our proposed deterministic sampler in this work) outperforms DEIS by a large margin for all sampling budgets.

References:

[1] gDDIM: Generalized denoising diffusion implicit models

[2] A Complete Recipe for Diffusion Generative Models

[3] Elucidating the Design Space of Diffusion-Based Generative Models

We thank the reviewer for their insight into our work and respond to each question below:

Combinations between datasets and methods are not consistent. For example in Figure 5, EDM is not tested on CelebA-64. Furthermore, these methods must be tested in more complex distributions such as ImageNet.

Response:

Combinations between datasets and methods are not consistent. For example in Figure 5, EDM is not tested on CelebA-64.

Response: In Figure 5 (in the main text), we report FID scores for different baseline samplers based on the results in their respective papers or by using pre-trained checkpoints from their official code implementations. Consequently, for a given dataset, we do not include results for a baseline if the corresponding pre-trained diffusion checkpoint is missing or the paper lacks results for that dataset.

Furthermore, these methods must be tested in more complex distributions such as ImageNet.

Since we demonstrate the proposed samplers in the context of PSLD (a recently proposed state-of-the-art diffusion model), due to the absence of pre-trained models for this type of diffusion model, we would need to re-train PSLD models on large datasets like ImageNet. Training large-scale diffusion models is typically computationally prohibitive, and thus in this work, we only consider small-scale datasets for which training fits within our compute limits or datasets for which pre-trained PSLD checkpoints are available. For instance, for training on ImageNet at 64x64 resolution, the authors in [3], train the diffusion model for two weeks on a setup with 32 GPUs which is computationally expensive.

The training of each pretrained model and the size of the model is not specified in the main paper. Are these models identical and trained for the same amount of time?

Response: We thank the reviewer for pointing this out. For most competitive baseline results presented in Figure 5 (like DEIS, DPM-Solver, etc.), which use the VP-SDE diffusion baseline, the size of the pre-trained diffusion model is around 108M parameters, while the size of the pre-trained PSLD model used to evaluate our proposed sampler is 97M parameters. Both models are trained for approximately the same number of steps for CIFAR-10 (around 800k-1M steps). Therefore, our comparisons are fair in terms of sample quality and time for each score function evaluation. We have added these details in Appendix E.2.

In Collorary 1, the authors show a connection between stability and the parameter λ (which should evolve with time?), however it is not clear how this hyperparrameter is chosen. Furthermore, Corollary 1 does not explain the behaviour of λ-DDIM-II, which performs best.

Response: We would like to point out that the main motivation behind including Theorem 2 is to give intuition about when projecting to the space $\hat{z}_t$ (which is parameterized by the design choice $B_t$) is a good idea. Furthermore, since we employ simplistic choices of $B_t = \lambda I$, we include Corollary 1 to provide further intuition since the expression for $\bar{\Lambda}$ simplifies. With this motivation in mind, we address the reviewer’s question in two parts:

In Collorary 1, the authors show a connection between stability and the parameter λ (which should evolve with time?), however it is not clear how this hyperparrameter is chosen

In this work, the parameter $\lambda$ is chosen at inference using grid search. Specifically, given a step size h, we tune $\lambda$ using grid search to optimize for best sample quality. In principle, while $\lambda$ can be time-dependent, we use a rather simplistic choice by setting $\lambda$ to be a constant scalar across all timesteps. While this might seem overly simplistic (due to the result in Corollary-1, which suggests that $\lambda$ would be rather time-dependent), it works well empirically and significantly reduces the effort required in tuning $\lambda$ during inference.

Furthermore, Corollary 1 does not explain the behaviour of λ-DDIM-II, which performs best.

As stated above, choosing $B_t = \lambda I$ simplifies the expression for $\bar{\Lambda}$. However, the same is not true for the choice of $B_t = \lambda 1$, and thus it is difficult to argue about the better performance of $\lambda$-DDIM-II over $\lambda$-DDIM-I. Therefore, further investigation into different theoretical aspects of Conjugate Integrators remains an exciting direction for future work. We already make this point more explicit in our future directions section at the end of the main text.

Improvements seem not to be good enough, especially considering speed. Large datasets and models are also not considered.

Response: We address this comment in two parts:

Improvements seem not to be good enough, especially considering speed

We respectfully disagree and would like to re-highlight some important empirical contributions:

In terms of absolute comparisons with state-of-the-art samplers:

1. For deterministic samplers proposed in this work, for sampling budgets between 50 and 100 steps, we outperform existing samplers like DEIS and DPM-Solver across all datasets considered in this work. We find that our proposed sampler for PSLD achieves a FID score of 2.11 in 100 steps (compared to 2.57 for the best-performing baseline DEIS) for CIFAR-10 (See Table 1 and Fig. 5 in the main text). We observe similar results for other datasets like CelebA and AFHQv2.
2. For Stochastic Samplers proposed in this work, we again outperform state-of-the-art solvers like EDM [6], where our best stochastic sampler achieves a FID score of 2.36 in 100 function evaluations (as compared to 2.63 for the baseline).

In terms of generalizing existing samplers like DPM-Solver, DEIS, etc., and applying them to the PSLD diffusion, our proposed sampler (CSPS-D) outperforms these methods by a large margin, as illustrated in Table 1 (below in the response).

Large datasets and models are also not considered

Since we demonstrate the proposed samplers in the context of PSLD (a recently proposed state-of-the-art diffusion model), due to the absence of pre-trained models for this type of diffusion model, we would need to re-train PSLD models on large datasets like ImageNet. Training large-scale diffusion models is typically computationally prohibitive, and thus in this work, we only consider small-scale datasets for which training fits within our compute budget or datasets for which pre-trained PSLD checkpoints are available. For instance, for training on ImageNet at 64x64 resolution, the authors in [3] train the diffusion model for two weeks on a setup with 32 GPUs which is computationally expensive. Moreover, we believe that the advantages of our sampling framework are already well-demonstrated on smaller-scale datasets considered in this work.

Minor Comments: We thank the reviewer for pointing out these issues, which we will fix in our final revision.

Questions:

how the dpm-solver will formulate under your settings and when will your method outperform it?

We already provide a response to this question in point 5 of our response above.

Additional Results

Table 1: Comparison between different ODE-solvers for the CIFAR-10 dataset at various sampling budgets. DEIS was generalized using the Conjugate Integrator framework and applied to PSLD reverse diffusion ODE. DEIS (q=k) represents DEIS with polynomial extrapolation (where k represents the polynomial order), which is known to improve performance for low NFEs [5]. CSPS-D (our proposed deterministic sampler in this work) outperforms DEIS by a large margin for all sampling budgets.

References

[1] High-Resolution Image Synthesis with Latent Diffusion Models, Rombach et al.

[2] Photorealistic Text-to-Image Diffusion Models with Deep Language Understanding

[3] Elucidating the Design Space of Diffusion-Based Generative Models

[4] Score-Based Generative Modeling through Stochastic Differential Equations

[5] A Complete Recipe for Diffusion Generative Models

[6] Fast Sampling of Diffusion Models with Exponential Integrator

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

[8] Denoising Diffusion Implicit Models

This work, also reminds me of the DPM-Solvers. DPM-Solver could be viewed as a special case, it also projects the integral to another space, the lambda space in fact. It uses a prediction-correction (PC) method to achieve much higher accuracy in the projected space. So I care about two things: first how the dpm-solver will formulate under your settings and when will your method outperform it? Second, can your method benefit from PC? While this work does not compare with DPM-Solver, considering the huge influences of and similarity with DPM-Solver, it could be better if the authors could compare with it.

Response: We address this question in multiple parts.

DPM-Solver as a particular case of Conjugate Integrators

We thank the reviewer for their intuitive insight. During the response period, we established a theoretical connection of Conjugate Integrators with DPM-Solver [1]. More specifically, we find that the DPM Solver sampler for the position space-only diffusion models considered in [1] is a particular case of conjugate integrators with the choice of matrix $B_t = 0$. We have updated Proposition 2 in the main text with proof in Appendix B.3 to reflect the same.

When does our method outperform DPM-Solver?

Our method can outperform DPM-Solver (or, more broadly, other exponential integrator-based methods like DEIS) in the following use cases.

1. Exponential-based integrators do not generalize automatically to diffusion models in augmented spaces: Firstly, with the perspective presented in DPM Solver [1] and other exponential-based integrators like DEIS [5], it is not apparent how to generalize these samplers to a broader class of diffusion models i.e. augmented diffusion models (where diffusion is performed in the data space concatenated with some auxiliary variables like in CLD [2] or PSLD [3]) as considered in this work. Theoretically, one of our main contributions is a framework that generalizes ideas in fast diffusion sampling like DDIM [4], DPM-Solver [1], and DEIS [5] under a single unified framework (see Propositions 1,2 in the main text). The benefit of such a framework is that it allows the construction of analogous samplers for different types of diffusion models like CLD and PSLD (which could be helpful in different contexts) and not just position space-only diffusion models (like Stable Diffusion). Moreover, a principled framework like ours allows us to extend beyond exponential integrator-based samplers. The resulting samplers can generalize better to other types of diffusion models.

   To illustrate our point, empirically, we compare the FID scores of DEIS (generalized to PSLD by setting $B_t=0$) with “\lambda-DEIS” (obtained by setting $B_t=\lambda 1$) applied to the PSLD reverse diffusion in Table 1 (in this response) for the CIFAR-10 dataset. Our extension of DEIS i.e. $\lambda$-DEIS, outperforms the former by a large margin. We expect similar results for DPM-Solver since DEIS and DPM-Solver are exponential integrator-based methods and exhibit similar empirical performance in practice [5]. These results demonstrate that the exponential integrator-based framework does not generalize to novel diffusions automatically, therefore justifying the need for a broader framework as proposed in this work._

2. DPM Solver saturates in sample quality for higher budgets. When comparing with DPM Solver (applied to VP-SDE, which is a non-augmented diffusion model), we find that exponential integrator-based methods converge faster but to a sub-optimal FID. In contrast, our proposed sampler (though applied to PSLD) eventually converges to better overall FID scores for slightly larger budgets (FID=2.11 for our sampler vs 2.59 for DPM-Solver). Please refer to the top-left subfigure in Figure 5 in the main text, where we compare with DPM-Solver and other sampler baselines. Therefore, our method is better suited for applications with a slightly larger computing budget at inference.

Can the proposed samplers be combined with PC?

Yes, our samplers can be extended to incorporate PC-based sampling as well since special cases like DPM-Solver and DEIS can automatically incorporate correction steps.

The stability analysis seems to avoid the key points? I would like to know a two-phase conclusion: first, when will matrix A_t stability have an inverse, this is the key to the stable conjuncture; second, under the condition that A_t is inversible and stable, when will eq 7 and 8 as ODE solvers be stable? Theorem 2 seems like a naive application of ODE stability conditions, and gets rid of the key part. A_t is an integral of multiple components, so we may not directly assume it to be stably invertible. But you also do not assume it stable anyway (minimum eigenvalue larger than some positive constant). It would also be helpful if the authors could explain in what reality settings those conditions can be satisfied.

Response: We thank the reviewer for pointing this issue out. We address the question in two parts.

When will A_t be invertible:

Theoretically, since the expression for $A_t$ in Eqn. 6 consists of multiple terms in the integral with the matrices $B_t$, $F_t$, and $C_\text{skip}(t)$ having no specific properties, it's non-trivial to conclude if arbitrary choices of these matrices can ensure invertibility on $A_t$. An alternate (but less elegant) fix would be to update our mapping instead to ensure an invertible transformation from $z_t$ to $\hat{z}_t$ by setting $\hat{z}_t = (A_t + \delta I) z_t$ where $\delta > 0$ is a small constant to ensure non zero eigenvalues at any time t.

However, empirically, when using the PSLD diffusion model (as considered in this work), we do not notice any instabilities because of this issue across our experiments with multiple datasets. We have highlighted this aspect in Section 3.1 and Appendix B.5 in the paper.

When will the ODE solvers derived in Eqn. 8 be stable:

Theorem 2 states when the ODE integrator in Eqn. 8 is stable, provided the mapping $A_t$ is invertible. We have updated Theorem 2 to reflect the assumption on the invertibility of $A_t$.

Theorem 3 gives error analysis based on h, I have concerns about its error with respect to the original t. After the projection is to solve the original problem, measuring the error in a projected space could be less meaningful. This is the same problem in the DPM-Solver paper, they only analyze errors when lambda is stably transformed from t, but in reality is often not the case. But DPM-Solver performs great generally so I am not criticizing it. Just point out that error analysis with respect to t variable is preferred for typical ode numerical analysis.

Response: We think there might be some misunderstanding here.

Firstly, we would like to clarify that Theorem 3 presents the error analysis for Naive Splitting Integrators (specifically for the Naive Velocity Verlet (NVV) sampler). As presented in Section 3.2, Splitting Integrators are independent of the Conjugate Integrators presented in Section 3.1 and, thus, do not involve a projection step. Therefore, the truncation error in Theorem 3 is measured in the original rather than the projected space. Furthermore, we don’t combine Conjugate Integrators with Splitting Integrators until Section 3.3. To clarify further, in Section 3.3, we numerically solve individual splitting updates using conjugate integrators.

Secondly, the error analysis in Theorem 3 is based on $h$, which is the step size during numerical integration. We would like to point out that error analysis based on the integration step size is commonplace in the numerical analysis of ODEs. We hope this clarifies any misunderstandings about the proposed samplers.

We thank the reviewer for their insight into our work and respond to each question below:

I don't know whether I deviate too much from the author's experiences. In my experience, in most cases, 25-step DPM-solver sampling could produce good enough sampling results for large models like StableDiffusions on large diverse image distribution; and less than 15-step DPM-solver sampling could produce good enough results for small diffusion models in small data domain. Please tell me if you observe a different phenomenon, and I will be happy to test it. Getting back to the subject, are the 100 and 50 NFEs considered by this paper a bit redundant?

Response: For non-augmented diffusion models (i.e., without auxiliary variables) like Stable Diffusion, indeed, DPM-Solver can produce plausible samples in 15-25 steps. However,

1. Firstly, Stable Diffusion [1] is position space only, i.e., there are no auxiliary variables in its design. In this work, we consider the sampler design for a broader class of diffusion models, and it's not apparent whether DPM-Solver works effectively in such cases. In this work, we propose Conjugate Integrators, which allow us to bring samplers like DEIS [6], DPM-Solver [7] (see our response to your point 5), and DDIM [8] in a common framework. This further allows us to extend these samplers in a principled way to a more broad class of diffusion models like PSLD. This is illustrated in Table 1 below, where we first extend DEIS (which is an exponential-based integrator like DPM-Solver) to PSLD and show that it does not generalize well to diffusions in augmented spaces like PSLD even for NFE=50 and NFE=100, which justifies our use of 100 and 50 NFEs in this work.
2. Secondly, for data-space-only diffusions like VP-SDE, it is commonly observed that state-of-the-art ODE-based samplers like DPM-Solver and DEIS converge quickly in fewer steps but saturate to a sub-optimal FID for even smaller datasets like CIFAR-10 (see the Top-left subfigure in Figure 5) and therefore might not be suitable for applications which might allow for a larger sampling budget. On the contrary, though applied to the PSLD diffusion, our proposed sampler in Figure 5 converges to a better FID score of 2.11 on CIFAR-10 (vs. 2.59 for DPM-Solver). Therefore, when evaluating diffusion models, we think it is more reasonable to evaluate sampler quality on a range of sampling budgets (visualized as sampling speed in NFE vs sample quality in FID).

Overall, we think the sampling budgets of 50 and 100 NFEs are not redundant for the results presented in this paper.

In eq 1, the authors actually simplify the original diffusion process to linear ones, and replace the f(z_t) function to linear map F_t z_t, so do please kindly mention that in the corresponding place. Otherwise, it may cause confusion. I care about the scope of this simplification, could it be applied to generation diffusion processes? Especially those in large-scale datasets and models? For example,

Response: Thanks for pointing this out. We have added a brief note in the background specifying that we primarily consider diffusion processes with affine drift, i.e., $f(z_t) = F_t z_t$. To the best of our knowledge,

1. The assumption of an affine drift is an essential component for obtaining tractable perturbation kernels $p(x_t|x_0)$ needed for training using denoising score matching.
2. Secondly, the assumption works well empirically and is employed in most state-of-the-art diffusion models like EDM [3], VP-SDE [4], and PSLD [5]. Therefore, diffusion models with affine drifts are widely adopted for most state-of-the-art diffusion models (including large-scale diffusion models like Stable Diffusion [1], Imagen [2], etc.)

Thank you for the response.

I have an immediate question about the stability issue of eq 70.

Notice that you have ABA^-1 in the integral. The ode will be very unstable (huge numerical error) if A has a small eigenvalues, adding a small delta will not help it, as it still has an error larger than O(1/sigma).

By the way, noting eq 7, if A^-1 or ABA^-1 is infinite, then this sde does not have a solution as it as infinite Lipschitz on z.

We thank the reviewer for the follow-up questions. There are two aspects to this issue: Theoretical and empirical.

Theoretically, for the transformation $A_t + \delta I$, the minimum eigenvalue will be $1/(\lambda_\text{min} + \delta)$. Thus, with a small enough delta, we can ensure that the overall transformation matrix $A_t + \delta I$ is invertible (on a machine with finite numerical precision). Moreover, since our framework does not require additional training during inference, $\delta$ could even be a hyperparameter tuned during inference (whether this is even required is an empirical question, which we address below) to obtain a stable inverse.

Given that our transformation is invertible (using a $\delta >0$ parameter to stabilize eigenvalues), regarding whether the ODE in Eqn. 7 has a solution, or whether the integrator in Eqn. 8 has significant discretization errors due to small eigenvalues of $A_t^{-1}$, we would like to highlight that both of these cases depend not only on the invertibility of $A_t$ but also on the design parameter $B_t$ which is essential in the design of conjugate integrators. For instance, for $B_t = \lambda I$ (equivalent to $\lambda$-DDIM-I in this work), the expression $A_tB_tA_t^{-1}$ simplifies to $\lambda I$ thus bypassing all problems stated by the reviewer. Thus, the overall stability of both the ODE in Eqn. 7 and our conjugate integrator in Eqn. 8 depend on $B_t$ and not just on the invertibility of $A_t$. A more formal analysis of the stability of the ODE in Eqn. 7 would likely involve analysis of the eigenvalues of the jacobian of the vector field in Eqn. 7 at an equilibrium point. We acknowledge that this could be an interesting direction for future work.

Empirically, we find that for the PSLD ODE, NFE=100 (i.e., 100 discretized time points) with $\delta=0$, $B_t=\lambda 1$ (i.e., $\lambda$-DDIM-II), the minimum eigenvalue of $A_t$ across all time points is around 1.0089 (which is near $t=\epsilon=1e-3$). Therefore, in practice, we do not observe such instabilities and the inverse of the transformation can be easily computed.

We will make some of these points more explicit in our final revision. Please feel free to revert back if you have any additional questions.

Thank you for your follow-up. From your response, we sense a couple of misconceptions (please feel free to correct us if we misunderstood something in your response).

1. Firstly, Theorem 3 states the error analysis for Splitting Integrators. We once again highlight that the error analysis in Theorem 3 is strictly applicable to Splitting Integrators and not Conjugate Integrators (which involve a projection step) and, therefore, separate from this discussion. (Perhaps you meant Theorem 2?)
2. We specifically deal with ODEs in Conjugate Integrators in this work, not SDEs. Therefore, Ito's calculus is not applicable here, so there is no notion of a strong/weak solution. (Though extension to SDEs can be a potential direction for future work, which we also highlight in the conclusion.)

With these points in mind, let's consider the following aspects:

Existence of a solution: This is a good point. Below, we examine the boundedness of different terms in the ODE in Eqn. 7 using matrix norms (Frobenius norms in particular denoted as ||.|| in subsequent discussion):

1. Boundedness of $A_t B_t A_t^{-1}$: It is straightforward to see that $||A_t B_t A_t^{-1}|| \leq \kappa_t ||B_t||$ where $\kappa_t = ||A_t|| ||A_t^{-1}||$ is the time-dependent condition number. Thus, for this matrix to be bounded, we need $\kappa_t ||B_t||$ to be bounded. This implies that the upper bound on the norm of $A_t B_t A_t^{-1}$ is not only dependent on $\kappa_t$ but also on $||B_t||$. Thus, this regularity condition provides a guideline for the end-user to choose $B_t$ appropriately. We will add this regularity condition to the theorem.

2. Boundedness of $d\Phi_t/dt$: Similarly, it can be shown that $||d\Phi_t/dt|| \leq 1/2 ||A_t|| ||G_t||^2 ||C_\text{out}||$. Since the norm of the diffusion coefficient $G_t$ is bounded, the boundedness of the norm of $||d\Phi/dt||$ primarily depends on the norms $||A_t||$ and $||C_\text{out}||$. Thus, our second regularity condition is that the norm of $A_t$ and $C_\text{out}$ must be bounded.

3. Boundedness of $\epsilon_{\theta}(z_t, t)$: This is our last regularity condition which usually holds for modern neural network design.

Therefore, for our vector field in Eqn. 7 to be bounded, we have three regularity conditions that we highlight above. We thank the reviewer for pointing out these issues and will include the proposed regularity conditions in the extended statement for our theorem. It is worth noting that we have tried to be generic in our treatment of the proposed regularity conditions. The exact magnitude of these bounds will depend on the specific choices of design matrices like $F_t$, $G_t$, $C_\text{skip}$, $C_\text{out}$ and $B_t$.

Regarding overall integration time, we already provide estimates for computing coefficients $\Phi_t$ in Appendix B.6 (see response to Reviewer eNRC). We will include the total wall clock running time for different NFE budgets in our subsequent revision. Additionally, we would like to highlight that empirically, in our experiments, we always set $\delta=0$.

Regarding the $(A_t + \delta)$ issue, as already shown above, the norm of $||A_tB_tA_t^{-1}||$ depends on $\kappa_t ||B_t||$. Thus, a time-dependent $\delta$ might not be required since we control the norm $||B_t||$ during sampling anyway by tuning $B_t$ as illustrated in this paper. This automatically helps us control the bound on $||A_tB_tA_t^{-1}||$ and allows us to use a small $\delta$ independent of time.

Regarding stable invertibility of A_t: As pointed out by the reviewer, we will mention this condition more explicitly in the revised version of our theorem.

We hope our response provides more theoretical intuition and clarifies your concerns.

References

[1] Score-Based Generative Modeling through Stochastic Differential Equations

[2] Score-Based Generative Modeling with Critically-Damped Langevin Diffusion

[3] A Complete Recipe for Diffusion Generative Models

[4] Elucidating the Design Space of Diffusion-Based Generative Models

[5] Denoising Diffusion Implicit Models

[6] Fast Sampling of Diffusion Models with Exponential Integrator

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