text:
\section{Introduction}

Generative models based on diffusion models~(DMs) have experienced rapid developments in the past few years and show competitive sample quality compared with generative adversarial networks~(GANs)~\cite{Dhariwal2021a,RamDhaNic22,Rombach2021}, competitive negative log likelihood compared with autoregressive models in various domains and tasks~\cite{song2021maximum, kawar2021snips}.
Besides, DMs enjoy other merits such as stable and scalable training, insensitivity to hyperparameters, and mode-collapsing resiliency~\cite{song2021maximum,Nichol2021ImprovedDD}.
However, slow and expensive sampling prevents DMs from further application in more complex and higher dimension tasks.
Once trained, GANs only forward pass neural network once to generate samples, but the vanilla sampling method of DMs needs 1000 or even 4000 steps~\cite{Nichol2021ImprovedDD,HoJaiAbb20,SonSohKin20} to pull noise back to the data distribution, which means thousands of neural networks forward evaluations.
Therefore, the generation process of DMs is several orders of magnitude slower than GANs.

How to speedup sampling of DMs has received significant attention.
With the seminal work by~\citet{SonSohKin20} on the connection between stochastic differential equations (SDEs) and diffusion models, several studies investigate acceleration via solving the corresponding SDEs more efficiently~\cite{jolicoeur2021gotta,Nichol2021ImprovedDD,tachibana2021taylor}.
Recently, a remarkable work by~\citet{DocVahKre}, \textit{critically-damped Langevin diffusion~(CLD)}, proposes a velocity-augmented diffusion process and shows improved denoising efficiency.
Most methods along this line still need hundreds of steps to generate high fidelity samples.
Another promising strategy is based on \textit{probability flow}~\cite{SonSohKin20}, which is an ordinary
differential equations~(ODE) associated with DMs that share equivalent marginal with SDE.
Simple plug-in of off-the-shelf ODE solvers can already achieve significant acceleration comparing to SDEs-based methods~\cite{SonSohKin20,DocVahKre}.
One of the arguably most popular sampling methods is \textit{denoising diffusion implicit model~(DDIM)}~\cite{SonMenErm21}, which includes both deterministic and stochastic samplers, and both show tremendous improvement in sampling quality compared with previous methods when only a small number of steps is used for the generation.
Deterministic DDIM in fact reduces to probability flow in the infinitesimal step size limit~\cite{SonMenErm21,LiuRenLin22}. Meanwhile, various approaches have been proposed to accelerate DDIM~\cite{kong2021fast,watson2021learning,LiuRenLin22}.
Recently, \citet{ZhaChe22} discovered deterministic DDIM is an ODE numerical integrator based on exponential integrator, and employed exponential multistep method to further boost the efficiency of the ODE solver.
~(More related works in \cref{app:related}.)

Although significant improvements of the DDIM in sampling efficiency have been observed empirically, the understanding of the mechanism of the DDIM is still lacking.
First, why does solving probability flow ODE provide much higher sample quality than solving SDEs, when the number of steps is small?
Second, it is shown that stochastic DDIM reduces to marginal-equivalent SDE~\cite{ZhaChe22}, but its discretization scheme and mechanism of acceleration are still unclear.
Finally, can we generalize DDIMs to other DMs and achieve similar or even better acceleration results?

In this work, we conduct a comprehensive study to answer the above questions, so that we can generalize and improve DDIM.
We start with an interesting observation that the DDIM can solve corresponding SDEs/ODE {\bf exactly} without any discretization error in finite or even one step, when the training dataset consists of only one data point.
{%
For deterministic DDIM, we find that the added noise in perturbed data along the diffusion is constant along an exact solution of probability flow ODE (see \cref{prop:ddim-ode}).
}
Besides, provided only one evaluation of log density gradient~(a.k.a. score), we are already able to recover accurate score information for any datapoints, and this explains the acceleration of stochastic DDIM for SDEs (see \cref{prop:ddim-sde}).
{%
Based on this observation, together with the manifold hypothesis, we present one possible interpretation to explain why the discretization scheme used in DDIMs is effective on realistic datasets (see \cref{fig:manifold}).}
Equipped with this new interpretation, we extend DDIM to general DMs, which we coin \textit{generalized DDIM~(gDDIM)}. With only a small but delicate change of the model parameterization, gDDIM can accelerate DMs based on general diffusion processes. Specifically, we verify the effectiveness of gDDIM on CLD.

To summarize, we have made the following contributions:
1) We provide an interpretation for the DDIM and unravel the mechanism of it.
2) The interpretation not only justifies the numerical discretization of DDIMs but also provides insights on why ODE-based samplers are preferred over SDE-based samplers when NFE is low.
3) We propose gDDIM, a generalized DDIM that can accelerate a large class of DMs.
4) We show by extensive experiments that gDDIM can drastically improve sampling quality/efficiency almost for free. Specifically, when applied to CLD, gDDIM can achieve an FID score of 2.87 with only 27 steps and 2.28 with 50 steps.

\section{Background}
In this section, we provide a brief introduction to diffusion models (DMs).
Most DMs are built on two diffusion processes in continuous-time, one forward diffusion known as the noising process that drives any data distribution to a tractable distribution such as Gaussian by gradually adding noise to the data, and one backward diffusion known as the denoising process that sequentially removes noise from noised data to generate realistic samples.
The continuous-time noising and denoising processes are modeled by stochastic differential equations~(SDEs)~\cite{SarSol19}.

\begin{figure}[tbp]
    \centering
    \includegraphics[width=1.0\textwidth]{figure/fig1.pdf}
    \caption{
        Importance of ${\bm{K}}_t$ for score parameterization ${\bm{s}}_\theta({\bm{u}}, t)= -{\bm{K}}_t^{-T} {\epsilon}_\theta({\bm{u}}, t)$ and acceleration of diffusion sampling with probability flow ODE. Trajectory of probability ODE for CLD~\cite{DocVahKre} and random pixel locations~(Left). Pixel value and output of ${\epsilon}_\theta$ in ${\bm{v}}$ channel with choice ${\bm{K}}_t={\bm{L}}_t$~\cite{DocVahKre} along the trajectory~(Mid).
        Output of ${\epsilon}_\theta$ in ${\bm{x}}, {\bm{v}}$ channels with choice ${\bm{K}}_t={\bm{R}}_t$~(Right). The smooth network output enables large stepsize and thus sampling acceleration.
    }
    \label{fig:fig1}
\end{figure}

In particular, the forward diffusion is a linear SDE with state ${\bm{u}}(t)\in\mathbb{R}^D$
\begin{equation}\label{eq:sde}
    d{\bm{u}} = {\bm{F}}_t {\bm{u}} dt + {\bm{G}}_t d{\bm{w}}, t \in [0,T]
\end{equation}
where ${\bm{F}}_t, {\bm{G}}_t \in \mathbb{R}^{D\times D}$ represent the linear drift coefficient and diffusion coefficient respectively, and ${\bm{w}}$ is a standard Wiener process. When the coefficients are piece-wise continuous, \cref{eq:sde} admits a unique solution~\cite{oks13}. Denote by $p_t({\bm{u}})$ the distribution of the solutions $\{{\bm{u}}(t)\}_{0\leq t\leq T}$ (simulated trajectories) to \cref{eq:sde} at time $t$, then $p_0$ is determined by the data distribution and $p_T$ is a (approximate) Gaussian distribution. That is, the forward diffusion \cref{eq:sde} starts as a data sample and ends as a Gaussian random variable. This can be achieved with properly chosen coefficients ${\bm{F}}_t, {\bm{G}}_t$.
Thanks to linearity of \cref{eq:sde}, the transition probability $p_{st}({\bm{u}}(t) | {\bm{u}}(s))$ from ${\bm{u}}(s)$ to ${\bm{u}}(t)$ is a Gaussian distribution. For convenience, denote $p_{0t}({\bm{u}}(t)|{\bm{u}}(0))$ by ${\mathcal{N}}(\mu_t {\bm{u}}(0), \Sigma_t)$ where $\mu_t, \Sigma_t \in \mathbb{R}^{D \times D}$.

The backward process from ${\bm{u}}(T)$ to ${\bm{u}}(0)$ of \cref{eq:sde} is the denoising process.
Remarkably, it can be characterized by the backward SDE simulated in reverse-time direction~\cite{SonSohKin2020,Anderson1982}
\begin{equation}\label{eq:r-sde}
    d {\bm{u}} = [ {\bm{F}}_t {\bm{u}} dt - {\bm{G}}_t {\bm{G}}_t^T \nabla \log p_t({\bm{u}})]dt + {\bm{G}}_t d{\bm{w}},
\end{equation}
where ${\bm{w}}$ denotes a standard Wiener process running backward in time. Here $\nabla \log p_t ({\bm{u}})$ is known as the score function. When \cref{eq:r-sde} is initialized with ${\bm{u}}(T) \sim p_T$, the distribution of the simulated trajectories coincides with that of the forward diffusion \cref{eq:sde}. Thus, ${\bm{u}}(0)$ of these trajectories are
unbiased samples from $p_0$; the backward diffusion \cref{eq:r-sde} is an ideal generative model.

In general, the score function $\nabla \log p_t({\bm{u}})$ is not accessible.
In diffusion-based generative models, a time-dependent network ${\bm{s}}_\theta({\bm{u}}, t)$, known as the score network, is used to fit the score $\nabla \log p_t({\bm{u}})$.
One effective approach to train ${\bm{s}}_\theta({\bm{u}}, t)$ is the denoising score matching~(DSM) technique~\cite{SonSohKin2020, HoJaiAbb20,Vin11} that seeks to minimize the DSM loss
\begin{equation}\label{eq:loss-dsm}
    \mathbb{E}_{t \sim {\mathcal{U}}[0,T]}\mathbb{E}_{{\bm{u}}(0), {\bm{u}}(t)|{\bm{u}}(0)}
    [\norm{
        \nabla \log p_{0t}({\bm{u}}(t) | {\bm{u}}(0)) - {\bm{s}}_\theta({\bm{u}}(t), t)
    }_{\Lambda_t}^2],
\end{equation}
where ${\mathcal{U}}[0,T]$ represents the uniform distribution over the interval $[0,\,T]$.
The time-dependent weight $\Lambda_t$ is chosen to balance the trade-off between sample fidelity and data likelihood of learned generative model~\cite{song2021maximum}.
It is discovered in \citet{HoJaiAbb20} that reparameterizing the score network by
    \begin{equation}\label{eq:scoreoriginal}
        {\bm{s}}_\theta({\bm{u}}, t) = - {\bm{K}}_t^{-T} {\epsilon}_\theta({\bm{u}}, t)
    \end{equation}
with ${\bm{K}}_t {\bm{K}}_t^{T} = \Sigma_t$ leads to better sampling quality. In this parameterization, the network tries to predict directly the noise added to perturb original data. Invoking the expression ${\mathcal{N}}(\mu_t {\bm{u}}(0), \Sigma_t)$ of $p_{0t}({\bm{u}}(t)|{\bm{u}}(0))$, this parameterization results in the new DSM loss
\begin{equation}\label{eq:eps-dsm-loss}
    {\mathcal{L}}(\theta) = \mathbb{E}_{t \sim {\mathcal{U}}[0,T]}\mathbb{E}_{{\bm{u}}(0) \sim p_0, {\epsilon} \sim {\mathcal{N}}(0, {\bm{I}}_D)}
    [\norm{
        {\epsilon} - {\epsilon}_\theta(\mu_t {\bm{u}}(0) + {\bm{K}}_t {\epsilon}, t)
    }_{{\bm{K}}_t^{-1}\Lambda_t{\bm{K}}_t^{-T}}^2].
\end{equation}

\textbf{Sampling:} After the score network ${\bm{s}}_\theta$ is trained, one can generate samples via the  backward SDE \cref{eq:r-sde} with a learned score, or the marginal equivalent SDE/ODE~\cite{SonSohKin20,zhang2021diffusion, ZhaChe22}
\begin{equation}\label{eq:lambda-sde}
    d{\bm{u}} = [{\bm{F}}_t {\bm{u}} - \frac{1+\lambda^2}{2} {\bm{G}}_t {\bm{G}}_t^T {\bm{s}}_\theta({\bm{u}}, t)]dt + \lambda {\bm{G}}_t d{\bm{w}},
\end{equation}
where $\lambda \geq 0$ is a free parameter.
Regardless of the value of $\lambda$, the exact solutions to \cref{eq:lambda-sde} produce unbiased samples from $p_0({\bm{u}})$ if ${\bm{s}}_\theta({\bm{u}}, t) = \nabla \log p_{t}({\bm{u}})$ for all $t,{\bm{u}}$. When $\lambda = 1$, \cref{eq:lambda-sde} reduces to reverse-time diffusion in \cref{eq:r-sde}. When $\lambda=0$, \cref{eq:lambda-sde} is known as the \textit{probability flow ODE}~\cite{SonSohKin20}
\begin{equation}\label{eq:ode}
    d{\bm{u}} = [{\bm{F}}_t {\bm{u}} - \frac{1}{2} {\bm{G}}_t {\bm{G}}_t^T {\bm{s}}_\theta({\bm{u}}, t)]dt.
\end{equation}

\textbf{DDPM and DDIM:} A popular DM is the \textit{Denoising diffusion probabilistic modeling~(DDPM)}~\cite{HoJaiAbb20}. For a given data distribution $p_{\rm data}({\bm{x}})$, DDPM has state ${\bm{u}} = {\bm{x}} \in \mathbb{R}^d$ and sets $p_0({\bm{u}}) = p_{\rm data}({\bm{x}})$.
Though originally proposed in the discrete time setting, DDPM can be viewed as a discretization of a continuous-time SDE with parameters
    \begin{equation}\label{eq:DDPM}
        {\bm{F}}_t= \frac{1}{2} \frac{d \log \alpha_t}{d t} {\bm{I}}_d,\quad {\bm{G}}_t = \sqrt{-\frac{d \log \alpha_t}{dt}} {\bm{I}}_d
    \end{equation}
for a decreasing scalar function $\alpha_t$ satisfying $\alpha_0 = 1, \alpha_T = 0$.

The sampling scheme proposed in DDPM is inefficient; it requires hundreds or even thousands of steps, and thus number of score function evaluations~(NFEs), to generate realistic samples.
A more efficient alternative is the \textit{Denoising diffusion implicit modeling~(DDIM)} proposed in~\citet{SonMenErm21}.
It has a similar training loss \cref{eq:eps-dsm-loss} as DDPM but a different sampling scheme
\begin{equation}\label{eq:ddim}
    {\bm{x}}({t_{i-1}}) = \sqrt{\frac{\alpha_{t_{i-1}}}{\alpha_{t_{i}}}}{\bm{x}}({t_i}) +
    (
        \sqrt{1 - \alpha_{t_{i-1}} - \sigma_{t_i}^2} -
        \sqrt{1 - \alpha_{t_i}}
        \sqrt{\frac{\alpha_{t_{i-1}}}{\alpha_{t_{i}}}}
    )
    {\epsilon}_{\theta}({\bm{x}}({t_i}), t_i)
    + \sigma_{t_i} {\epsilon},
\end{equation}
where $\{\sigma_{t_i}\}$ are hyperparameters and ${\epsilon} \sim {\mathcal{N}}(0,{\bm{I}}_d)$. DDIM can generate reasonable samples within 50 NFEs.
For the special case where $\sigma_{t_i}=0$, it is recently discovered in \citet{ZhaChe22} that \cref{eq:ddim} coincides with the numerical solution to \cref{eq:ode} using an advanced discretization scheme known as the \textit{exponential integrator~(EI)} that utilizes the semi-linear structure of \cref{eq:ode}.

\textbf{CLD:} In \citet{DocVahKre} the authors propose \textit{critically-dampled Langevin diffusion~(CLD)}, a DM based on an augmented diffusion with an auxiliary velocity term. More specifically, the state of the diffusion in CLD is of the form ${\bm{u}}(t) = [{\bm{x}}(t), {\bm{v}}(t)] \in \mathbb{R}^{2d}$ with velocity variable ${\bm{v}}(t) \in \mathbb{R}^d$. The CLD employs the forward diffusion \cref{eq:sde} with coefficients
\begin{equation}\label{eq:cld-fg}
    {\bm{F}}_t :=
        \begin{bmatrix}
        0     & \beta M^{-1} \\
        \beta & -\Gamma \beta M^{-1}
        \end{bmatrix}
     \otimes {\bm{I}}_d,
    \quad
    {\bm{G}}_t :=
        \begin{bmatrix}
            0     & 0 \\
            0 & -\Gamma \beta M^{-1}
            \end{bmatrix}
     \otimes {\bm{I}}_d.
\end{equation}
Here $\Gamma>0, \beta>0, M>0$ are hyperparameters.
Compared with most other DMs such as DDPM that inject noise to the data state ${\bm{x}}$ directly, the CLD introduces noise to the data state ${\bm{x}}$ through the coupling between ${\bm{v}}$ and ${\bm{x}}$ as the noise only affects the velocity component ${\bm{v}}$ directly. It is argued that this property improves the performance of DMs \citep{DocVahKre}. It is observed in \citet{DocVahKre} that CLD has more efficient noising and denoising processes and can generate high-quality samples with less NFEs compared with DDPM.  See \cref{sec:CLD} for more discussions on CLD.

\section{Revisit DDIM: Gap between exact solution and numerical solution}\label{sec:DDIM}
The complexity of sampling from a DM is proportional to the NFEs used to numerically solve \cref{eq:lambda-sde}.
To establish a sampling algorithm with a small NFEs, we ask the bold question:

\centerline{
\textit{Can we generate samples accurately from a DM with finite steps if the score function is precise?}
}

To gain some insights to this problem, we start with the simplest scenario where the training dataset consists of only one data point ${\bm{x}}_0$.
It turns out that accurate sampling from diffusion models on this toy example is not that easy, even if the exact score function is accessible.
Most well-known numerical methods for \cref{eq:lambda-sde}, such as Euler integrator and Runge Kutta~(RK) for ODE, Euler-Maruyama~(EM) for SDE, are accompanied by discretization error and cannot recover the single data point in the train set unless infinite number of steps are used.
Indeed, in general, the exact solution to \cref{eq:lambda-sde} requires the score information along the entire continuous-time solution trajectory while any numerical method can only evaluate the score function finite times.

Surprisingly, as we show below,
DDIMs can recover the single data point in this toy example in finite steps, in fact, in one step.
We first consider the deterministic DDIM, that is, \cref{eq:ddim} with $\sigma_{t_i} =0$. In view of \cref{eq:DDPM}, the score network \cref{eq:scoreoriginal} is ${\bm{s}}_\theta({\bm{u}},t) = -\frac{{\epsilon}_\theta({\bm{u}}, t)}{\sqrt{1-\alpha_t}}$. To differentiate between the learned score and the real score, denote the ground truth version of ${\epsilon}_\theta$ by ${\epsilon}_{\rm GT}$. In our toy example, the following property holds for ${\epsilon}_{\rm GT}$.
\begin{prop}\label{prop:ddim-ode}
Assume $p_0({\bm{u}})$ is a Dirac distribution. Let ${\bm{u}}(t)$ be an arbitrary solution to the probability flow ODE \cref{eq:ode} with coefficient \cref{eq:DDPM} and the ground truth score, then ${\epsilon}_{\rm GT}({\bm{u}}(t), t) = -\sqrt{1 - \alpha_t} \nabla \log p_t({\bm{u}}(t))$ remains constant along ${\bm{u}}(t)$.
\end{prop}
\looseness=-1 We remark that even though ${\epsilon}_{\rm GT}({\bm{u}}(t), t)$ remains constant along an exact solution, the score $\nabla \log p_t({\bm{u}}(t))$ is time varying. This underscores the advantage of the parameterization ${\epsilon}_\theta$ over ${\bm{s}}_\theta$.
Inspired by \cref{prop:ddim-ode}, we devise a sampling algorithm as follows that can recover the exact data point in one step for our toy example. This algorithm turns out to coincide with the determinstic DDIM.
\begin{prop}\label{prop:ddim-ode-update-step}
With the parameterization ${\bm{s}}_\theta({\bm{u}},\tau) = -\frac{{\epsilon}_\theta({\bm{u}}, \tau)}{\sqrt{1-\alpha_\tau}}$ and the approximation ${\epsilon}_\theta({\bm{u}}, \tau) \approx {\epsilon}_\theta ({\bm{u}}(t), t)$ for $\tau \in [t-\Delta t, t]$, the solution to the probability flow ODE \cref{eq:ode} with coefficient \cref{eq:DDPM} is
    \begin{equation}\label{eq:ODEiter}
        {\bm{u}}(t-\Delta t) = \sqrt{
            \frac{\alpha_{t - \Delta t}}{\alpha_t}
         }{\bm{u}}(t) +
         (
             \sqrt{1 - \alpha_{t-\Delta t}} - \sqrt{1 - \alpha_{t}}
             \sqrt{
            \frac{\alpha_{t - \Delta t}}{\alpha_t}
         }
        ){\epsilon}_\theta ({\bm{u}}(t), t),
    \end{equation}
which coincides with deterministic DDIM.
\end{prop}
When ${\epsilon}_\theta = {\epsilon}_{\rm GT}$ as is the case in our toy example, there is no approximation error in \cref{prop:ddim-ode-update-step} and \cref{eq:ODEiter} is precise. This implies that deterministic DDIM can recover the training data in one step in our example. The update \cref{eq:ODEiter} corresponds to a numerical method known as the exponential integrator to the probability flow ODE \cref{eq:ode} with coefficient \cref{eq:DDPM} and parameterization ${\bm{s}}_\theta({\bm{u}},\tau) = -\frac{{\epsilon}_\theta({\bm{u}}, \tau)}{\sqrt{1-\alpha_\tau}}$. This strategy is used and developed recently in \citet{ZhaChe22}. \cref{prop:ddim-ode} provides the sights on why such a strategy should work.

The above discussions however do not hold for stochastic cases where $\lambda > 0$ in \cref{eq:lambda-sde} and $\sigma_{t_i} > 0$ in \cref{eq:ddim}.
Since the solutions to~\cref{eq:lambda-sde} from $t=T$ to $t=0$ are stochastic, neither $\nabla \log p_t({\bm{u}})$ nor ${\epsilon}_{\rm GT}({\bm{u}}, t)$ remains constant over time; both are affected by the stochastic noise. The denoising SDE \cref{eq:lambda-sde} is more challenging compared with the probability ODE since it injects additional noise to ${\bm{u}}(t)$.
The score information needs to remove not only noise presented in ${\bm{u}}(T)$ but also injected noise along the diffusion. In general, one evaluation of ${\epsilon}_\theta({\bm{u}}, t)$ can only provide the information to remove noise in current state ${\bm{u}}$; it cannot predict the injected noise. Can we do better? The answer is affirmative on our toy dataset. It turns out that $\nabla \log p_t({\bm{u}}(t))$ along any stochastic trajectory $\{{\bm{u}}(t)\}$ is determined by $\nabla \log p_T({\bm{u}}(T))$.
\begin{prop}\label{prop:ddim-sde}
Assume SDE coefficients \cref{eq:DDPM} and that $p_0({\bm{u}})$ is a Dirac distribution. Given any evaluation of the score function $\nabla \log p_s({\bm{u}}(s))$, one can recover $\nabla \log p_t({\bm{u}})$ for any $t, {\bm{u}}$ as
    \begin{equation}\label{eq:ddim-approx}
        \nabla \log p_t({\bm{u}}) = \frac{1 - \alpha_s}{1 - \alpha_t}
    \sqrt{\frac{\alpha_t}{\alpha_s}} \nabla \log p_s({\bm{u}}(s)) - \frac{1}{1 - \alpha_t}({\bm{u}} - \sqrt{\frac{\alpha_t}{\alpha_s}} {\bm{u}}(s)).
    \end{equation}
\end{prop}
The major difference between \cref{prop:ddim-sde} and \cref{prop:ddim-ode} is that \cref{eq:ddim-approx} retains the dependence of the score over the state ${\bm{u}}$. This dependence is important in canceling the injected noise in the denoising SDE \cref{eq:lambda-sde}.
This approximation \cref{eq:ddim-approx} turns out to lead to a numerical scheme for \cref{eq:lambda-sde} that coincide with the stochastic DDIM.

\begin{thm}\label{thm:ddims}
Given the parameterization ${\bm{s}}_\theta({\bm{u}},\tau) = -\frac{{\epsilon}_\theta({\bm{u}}, \tau)}{\sqrt{1-\alpha_\tau}}$ and the approximation ${\bm{s}}_\theta({\bm{u}}, \tau) \approx \frac{1 - \alpha_t}{1 - \alpha_\tau}
    \sqrt{\frac{\alpha_\tau}{\alpha_t}} {\bm{s}}_\theta({\bm{u}}(t),t) - \frac{1}{1 - \alpha_\tau}({\bm{u}} - \sqrt{\frac{\alpha_\tau}{\alpha_t}} {\bm{u}}(t))$ for $\tau \in [t - \Delta t, t]$, the exact solution ${\bm{u}}({t-\Delta t})$ to~\cref{eq:lambda-sde} with coefficient \cref{eq:DDPM} is
    \begin{equation}\label{eq:lambda-ddim}
        {\bm{u}}({t - \Delta t}) \sim
        {\mathcal{N}}(
            \sqrt{\frac{\alpha_{t - \Delta t}}{\alpha_t}} {\bm{u}}(t) +
            \left[
                - \sqrt{\frac{\alpha_{t - \Delta t}}{\alpha_t}} \sqrt{1-\alpha_t} +
                \sqrt{1-\alpha_{t - \Delta t} - \sigma_t^2}
            \right] {\epsilon}_\theta({\bm{u}}(t), t),
            \sigma_t^2 {\bm{I}}_d
        )
    \end{equation}
    with $\sigma_t = (1 -\alpha_{t - \Delta t})\left[
        1 - \left(
            \frac{1-\alpha_{t - \Delta t}}{1-\alpha_t}
        \right)^{\lambda^2}
        \left(
            \frac{\alpha_t}{\alpha_{t - \Delta t}}
        \right)^{\lambda^2}
    \right]$, which is the same as the DDIM \cref{eq:ddim}.
\end{thm}

Note that \cref{thm:ddims} with $\lambda=0$ agrees with \cref{prop:ddim-ode-update-step}; both reproduce the determinstic DDIM but with different derivations.

\textbf{Dirac $p_0({\bm{u}})$ and the manifold hypothesis:}
While \cref{prop:ddim-ode} and \cref{prop:ddim-sde} require the strong assumption that the data distribution is a Dirac, DDIMs work very effectively on realistic datasets, which may contain millions of datapoints~\cite{Nichol2021a}. Here we present one possible interpretation based on the manifold hypothesis~\cite{Roweis2000NonlinearDR}.

It is believed that real-world data lie on low-dimensional manifold~\cite{Tenenbaum2000AGG} embedded in a high-dimensional space and the data points are well separated in high-dimensional data space. For example, realistic images are scattered in pixel space and the distance between every two images can be very large if measured in pixel difference even if they are similar perceptually. To model this property, we consider a dataset consisting of $M$ datapoints $\{ {\bm{u}}^{(m)}\}_{m=1}^M$. The exact score is \vspace{-0.06cm}
\begin{equation*}
    \nabla \log p_t({\bm{u}}) = \sum_{m} w_m \nabla \log p_{0t}({\bm{u}} | {\bm{u}}^{(m)}), \quad w_m = \frac{p_{0t}({\bm{u}} | {\bm{u}}^{(m)})}{\sum_{m} p_{0t}({\bm{u}} | {\bm{u}}^{(m)})},
    \vspace{-0.06cm}
\end{equation*}
which can be interpreted as a weighted sum of $M$ score functions associated with Dirac distributions.
This is illustrated in \cref{fig:manifold}. In the red color region where the weights $\{w_m\}$ are dominated by one specific data ${\bm{u}}^{(m^*)}$ and thus
$\nabla \log p_t({\bm{u}}) \approx \nabla \log p_{0t}({\bm{u}} | {\bm{u}}^{(m^*)})$.
The phenomenon caused by manifold hypothesis partially justifies the Dirac distribution assumption in \cref{prop:ddim-ode} and \cref{prop:ddim-sde}, and explains the effectiveness of DDIMs.

\begin{figure}
    \begin{minipage}[c]{0.65\textwidth}
        \includegraphics[width=\textwidth]{figure/mix_gaussian.png}
    \end{minipage}\hfill
    \begin{minipage}[c]{0.34\textwidth}
      \caption{
        Manifold hypothesis and Dirac distribution assumption. We model images dataset as a mixture of well separated Dirac distribution and visualize diffusion process. Curves in {\red{\textbf{red}}} indicate high density area spanned by $p_{0t}({\bm{u}}(t)|{\bm{u}}(0))$ by different mode and region surrounded by them indicates the phase when $p_t({\bm{u}})$ is dominated by one mode while region surrounded by {\blue{\textbf{blue}}} one is for the mixed phase. In the red region, we can approximate data score by score of dominated Dirac distribution. Such property can help accelerate sampling.
      }
    \label{fig:manifold}
    \end{minipage}
    \vspace{-0.6cm}
  \end{figure}
\section{Generalize and improve DDIM}
\label{sec:gDDIM}
The DDIM is specifically designed for DDPMs. Can we generalize it to other DMs?
With the insights in~\cref{prop:ddim-ode,prop:ddim-sde}, it turns out that with a careful chosen ${\bm{K}}_\tau$, we can generalize DDIMs to any DMs with general drift and diffusion.
We coin the resulted algorithm the~\textit{Generalized DDIM~(gDDIM)}.
Concretely, we apply the gDDIM to CLD to illustrate our approach.

\subsection{Deterministic gDDIM with~\cref{prop:ddim-ode}}

\textbf{Toy dataset}:
Motivated by \cref{prop:ddim-ode}, we ask whether there exists an ${\epsilon}_{\rm GT}$ that remains constant along a solution to the probability flow ODE \cref{eq:ode}.
We start with a special case with data distribution $p_0({\bm{u}}) = {\mathcal{N}}({\bm{u}}_0, \Sigma_0)$.
It turns out that any solution to \cref{eq:ode} is of the form
    \begin{equation}\label{eq:usolution}
        {\bm{u}}(t) = \Psi(t, 0){\bm{u}}_0 + {\bm{R}}_t {\epsilon}
    \end{equation}
with a constant ${\epsilon}$ and a time-varying ${\bm{R}}_t \in \mathbb{R}^{D\times D}$ that satisfies ${\bm{R}}_0 {\bm{R}}_0^T= \Sigma_0$ and
\begin{equation}\label{eq:r}
    \frac{d {\bm{R}}_t}{d {\bm{t}}} = ({\bm{F}}_t + \frac{1}{2}{\bm{G}}_t {\bm{G}}_t^T \Sigma_t^{-1}){\bm{R}}_t.
\end{equation}
Here $\Psi(t, s)$ is the transition matrix associated with ${\bm{F}}_\tau$; it is the solution to $\frac{\partial \Psi(t,s)}{\partial t}={\bm{F}}_t \Psi(t,s), \Psi(s,s)={\bm{I}}_D$.
Interestingly, ${\bm{R}}_t$ satisfies ${\bm{R}}_t {\bm{R}}_t^T = \Sigma_t$ like ${\bm{K}}_t$ in \cref{eq:scoreoriginal}. Based on \cref{eq:usolution} and \cref{eq:r}, we extend \cref{prop:ddim-ode} to more general DMs.
\begin{prop}\label{prop:f-ddde}
Assume the data distribution $p_0({\bm{u}})$ is ${\mathcal{N}}({\bm{u}}_0, \Sigma_0)$.
    Let ${\bm{u}}(t)$ be an arbitrary solution to the probability flow ODE \cref{eq:ode} with the ground truth score, then ${\epsilon}_{\rm GT}({\bm{u}}(t), t): = -{\bm{R}}^T_\tau \nabla \log p_t({\bm{u}}(t))$ remains constant along ${\bm{u}}(t)$.
\end{prop}
A direct consequence of \cref{prop:f-ddde} is that we can conduct accurate sampling in one step in the toy example since we can recover the score along any simulated trajectory given its value at $t=T$, if ${\bm{K}}_t$ in \cref{eq:scoreoriginal} is set to be ${\bm{R}}_t$.
This choice ${\bm{K}}_t = {\bm{R}}_t$ will make a huge difference in sampling quality as we will show later.
The fact provides guidance to design efficient sampling scheme for realistic data.

\textbf{Realistic dataset}: As the accurate score is not available for realistic dataset, we need to use learned score ${\bm{s}}_\theta({\bm{u}}, t)$ for sampling. With our new parameterization ${\epsilon}_\theta({\bm{u}}, t) = -{\bm{R}}_t^{T} {\bm{s}}_\theta({\bm{u}}, t)$ and the approximation $\tilde{{\epsilon}}_\theta({\bm{u}}, \tau) = {\epsilon}_\theta({\bm{u}}(t), t)$ for $\tau \in [t - \Delta t, t]$, we reach the update step for deterministic gDDIM by solving probability flow with approximator $\tilde{{\epsilon}}_\theta({\bm{u}}, \tau)$ exactly as
\begin{equation}\label{eq:iddim-ode-eps}
    {\bm{u}}(t - \Delta t) = \Psi(t - \Delta t, t) {\bm{u}}(t)
    +
    [
        \int_{t}^{t - \Delta t} \frac{1}{2} \Psi(t - \Delta t, \tau) {\bm{G}}_\tau {\bm{G}}^T_\tau{\bm{R}}_\tau^{-T}
    ]{\epsilon}_\theta({\bm{u}}(t), t).
\end{equation}

\textbf{Multistep predictor-corrector for ODE}:
Inspired by \citet{ZhaChe22}, we further boost the sampling efficiency of gDDIM by combining \cref{eq:iddim-ode-eps} with multistep methods~\cite{hochbruck2010exponential,ZhaChe22,LiuRenLin22}.
We derive multistep \textit{predictor-corrector} methods to reduce number of step while retaining accuracy~\cite{press2007numerical,sauer2005numerical}.
Empirically, we found that using more NFEs in predictor leads to better performance when the total NFE is small.
Thus, we only present multistep predictor for deterministic gDDIM. We include the proof and multistep corrector in~\cref{app:derivation}.
For time discretization grid $\{t_i \}_{i=0}^N$ where $t_0 = 0, t_N = T$, the $q$-th step predictor from $t_{i}$ to $t_{i-1}$ in term of ${\epsilon}_\theta$ parameterization reads
\begin{subequations}\label{eq:extrapolation}
\begin{align}
    {\bm{u}}(t_{i-1}) &= \Psi(t_{i-1}, t_i) {\bm{u}}(t_i) + \sum_{j=0}^{q-1}[{\textnormal{C}}_{ij} {\epsilon}_\theta({\bm{u}}(t_{i+j}), t_{i+j})], \label{eq:gddim-predictor-update} \\
    {\textnormal{C}}_{ij} &= \!\int_{t_i}^{t_{i-1}}  \frac{1}{2} \Psi(t_{i-1}, \tau)
    {\bm{G}}_\tau {\bm{G}}_\tau^T{\bm{R}}_\tau^{-T}
    \prod_{k \neq j} [\frac{\tau - t_{i+k}}{t_{i+j} - t_{i+k}}]
    d\tau. \label{eq:gddim-predictor-coef}
\end{align}
\end{subequations}

\subsection{Stochastic gDDIM with~\cref{prop:ddim-sde}}
Following the same spirits, we can generalize \cref{prop:ddim-sde}
\begin{prop}\label{prop:iddim-sde}
    Assume the data distribution $p_0({\bm{u}})$ is ${\mathcal{N}}({\bm{u}}_0, \Sigma_0)$.
    Given any evaluation of the score function $\nabla \log p_s({\bm{u}}(s))$, one can recover $\nabla \log p_t({\bm{u}})$ for any $t, {\bm{u}}$ as
    \begin{align}\label{eq:iddim-sde-approx-gt}
        \nabla \log p_t({\bm{u}}) = \Sigma_t^{-1} \Psi(t, s) \Sigma_s \nabla \log p_s({\bm{u}}(s)) - \Sigma_t^{-1}[{\bm{u}} - \Psi(t, s) {\bm{u}}(s)].
    \end{align}
\end{prop}
\cref{prop:iddim-sde} is not surprising; in our example, the score has closed form.
\cref{eq:iddim-sde-approx-gt} provides an accurate score estimation for our toy dataset; it also serves as score approximator in realistic data.

\textbf{Realistic dataset}:
Based on \cref{eq:iddim-sde-approx-gt}, with the parameterization ${\bm{s}}_\theta({\bm{u}}, \tau) = -{\bm{R}}_\tau^{-T} {\epsilon}_\theta({\bm{u}}, \tau)$, we propose the following gDDIM approximator $\tilde{{\bm{s}}}_\theta({\bm{u}}, \tau)$ for ${\epsilon}_\theta({\bm{u}}, \tau)$
\begin{equation}\label{eq:iddim-sde-approx}
    \tilde{{\epsilon}}_\theta({\bm{u}}, \tau) = {\bm{R}}_\tau^{-1} \Psi(\tau, s) {\bm{R}}_{s} {\epsilon}_\theta({\bm{u}}(s), s)+ {\bm{R}}_\tau^{-1}[{\bm{u}} - \Psi(\tau, s){\bm{u}}(s)].
\end{equation}
\begin{prop}\label{prop:s-gddim-update-step}
With the parameterization ${\epsilon}_\theta({\bm{u}}, t) = -{\bm{R}}_t^{T} {\bm{s}}_\theta({\bm{u}}, t)$ and the approximator $\tilde{{\epsilon}}_\theta({\bm{u}}, t)$ in~\cref{eq:iddim-sde-approx}, the solution to \cref{eq:lambda-sde} satisfies
    \begin{equation}\label{eq:iddim-sde-udpate-step}
    {\bm{u}}(t) \sim {\mathcal{N}}(\Psi(t, s){\bm{u}}(s) + [\hat{\Psi}(t, s) - \Psi(t,s)]{\bm{R}}_s{\epsilon}_\theta({\bm{u}}(s),s), {\bm{P}}_{st}),
    \end{equation}
    where $\hat{\Psi}(t, s)$ is the transition matrix associated with $\hat{{\bm{F}}}_\tau:= {\bm{F}}_\tau + \frac{1+\lambda^2}{2}{\bm{G}}_\tau {\bm{G}}_\tau^T\Sigma_\tau^{-1} $ and the covariance matrix ${\bm{P}}_{st}$ solves
\begin{equation}
    \frac{d {\bm{P}}_{s\tau}}{d\tau} = \hat{{\bm{F}}}_\tau {\bm{P}}_{s\tau} + {\bm{P}}_{s\tau} \hat{{\bm{F}}}_\tau^T + \lambda^2 {\bm{G}}_\tau {\bm{G}}_\tau^T, \quad {\bm{P}}_{ss} = 0.
\end{equation}
\end{prop}
Our stochastic gDDIM then uses \cref{eq:iddim-sde-udpate-step} for update.
Interestingly, though look quite different between stochastic gDDIM and deterministic gDDIM, there exists a connection when $\lambda=0$,
\begin{prop}\label{prop:s-gddim-update-lambda0}
    \cref{eq:iddim-sde-udpate-step} in stochastic gDDIM reduces to~\cref{eq:iddim-ode-eps} in deterministic gDDIM when $\lambda=0$.
\end{prop}

\subsection{Training and accelerating CLD}\label{sec:CLD}
We illustrate the gDDIM that is applicable to any DMs with a case study over CLD.
Denote by $ {\epsilon}_\theta({\bm{u}}, t) = [{\epsilon}_\theta({\bm{u}}, t; {\bm{x}}),{\epsilon}_\theta({\bm{u}}, t; {\bm{v}})]$ for corresponding model parameterization.
The authors of ~\citet{DocVahKre} originally propose the parameterization ${\bm{s}}_\theta({\bm{u}}, t) = -{\bm{L}}_t^{-T} {\epsilon}_\theta({\bm{u}}, t)$ where $\Sigma_t = {\bm{L}}_t {\bm{L}}_t^T$ is the Cholesky decomposition of the covariance matrix of $p_{0t}({\bm{u}}(t) | {\bm{x}}(0))$.
Built on DSM~\cref{eq:loss-dsm}, they propose \textit{hybrid score matching~(HSM)} that is claimed to be advantageous~\cite{DocVahKre}. It uses the loss
\begin{equation}\label{eq:cld-hsm-eps}
    \mathbb{E}_{t \sim {\mathcal{U}}[0,T]}\mathbb{E}_{{\bm{x}}(0), {\bm{u}}(t)|{\bm{x}}(0)}
    [\norm{
        {\epsilon} - {\epsilon}_\theta(\mu_t({\bm{x}}_0) + {\bm{L}}_t {\epsilon}, t)
    }_{{\bm{L}}_t^{-1} \Lambda_t {\bm{L}}_t^{-T}}^2].
\end{equation}
With a similar derivation~\cite{DocVahKre}, we obtain the HSM loss with our new score parameterization ${\bm{s}}_\theta({\bm{u}}, t) = -{\bm{L}}_t^{-T} {\epsilon}_\theta({\bm{u}}, t)$
as
\begin{equation}\label{eq:r-cld-hsm-eps}
    \mathbb{E}_{t \sim {\mathcal{U}}[0,T]}\mathbb{E}_{{\bm{x}}(0), {\bm{u}}(t)|{\bm{x}}(0)}
    [\norm{
        {\epsilon} - {\epsilon}_\theta(\mu_t({\bm{x}}_0) + {\bm{R}}_t {\epsilon}, t)
    }_{{\bm{R}}_t^{-1} \Lambda_t {\bm{R}}_t^{-T}}^2].
\end{equation}

Though~\cref{eq:cld-hsm-eps,eq:r-cld-hsm-eps} look similar, we cannot directly use pretrained model provided in~\citet{DocVahKre} for gDDIM.
Due to the lower triangular structure of ${\bm{L}}_t$ and the special ${\bm{G}}_t$, the solution to \cref{eq:lambda-sde} only relies on ${\epsilon}_\theta({\bm{u}}, t; {\bm{v}})$ and thus only ${\epsilon}_\theta({\bm{u}}, t; {\bm{v}})$ is learned in \citet{DocVahKre} via a special choice of $\Lambda_t$. In contrast, in our new parametrization, both ${\epsilon}_\theta({\bm{u}}, t; {\bm{x}})$ and ${\epsilon}_\theta({\bm{u}}, t; {\bm{v}})$ are needed to solve \cref{eq:lambda-sde}.
To train score model for gDDIM, we set ${\bm{R}}_t^{-1} \Lambda_t {\bm{R}}_t^{-T} = {\bm{I}}$ for simplicity, similar to the choice made in~\citet{HoJaiAbb20}.
Our weight choice has reasonable performance and we leave improvement possibilities, such as mixed score~\cite{DocVahKre}, better $\Lambda_t$ weights~\cite{song2021maximum}, for future work.
Though we requires a different training scheme of score model compared with~\citet{DocVahKre}, the modifications on the training pipeline and extra costs are almost ignorable. We include more discussions in~\cref{app:exp-impl}.
\section{Experiments}\label{sec:exp}
We validate the effectiveness of gDDIM on CLD.
We design experiments to answer the following questions.
How to verify~\cref{prop:f-ddde,prop:iddim-sde} empirically?
Can gDDIM improve sampling efficiency compared with existing works?
What differences does the choice of $\lambda$ and ${\bm{K}}_t$ make?
We present sampling results on CIFAR10.~(Implementation details and more experiments are included in~\cref{app:exp-impl})

\textbf{Choice of ${\bm{K}}_t$: ${\bm{L}}_t$ vs ${\bm{R}}_t$.} A key of gDDIM is the special choice ${\bm{K}}_t={\bm{R}}_t$ which is obtained via solving~\cref{eq:r}.
In CLD,
${\bm{L}}_t$~\cite{DocVahKre} does not obey~\cref{eq:r} and ${\bm{L}}_t \neq {\bm{R}}_t$.
As it is shown in~\cref{fig:fig1}, on real datasets with trained score model, we randomly pick pixel locations and check the pixel value and ${\epsilon}_\theta$ output along the solutions to the probability flow ODE produced by high resolution ODE solver.
With the choice ${\bm{K}}_t = {\bm{L}}_t$, ${\epsilon}^{({\bm{L}})}_\theta({\bm{u}}, t;{\bm{v}})$ suffers from oscillation like ${\bm{x}}$ value along time. However, ${\epsilon}^{({\bm{R}})}_\theta({\bm{u}}, t)$ is much more flat. We further compare samples generated by ${\bm{L}}_t$ and ${\bm{R}}_t$ parameterizaiton in \cref{tab:r_vs_l}.

\vspace{-0.2cm}
\begin{tabular}{cc}
\centering
\begin{minipage}[t]{.31\linewidth}
    \centering
    \captionof{table}{${\bm{L}}_t$ vs ${\bm{R}}_t$}\label{tab:r_vs_l}
    \begin{tabular}{lllll}
         \hline
          & \multicolumn{4}{c}{FID at different NFE}  \\
        ${\bm{K}}_t$ & 20 & 30 & 40 & 50 \\
        \hline
        ${\bm{L}}_t$ & 364  & 157  & 4.02 & 3.24 \\
        ${\bm{R}}_t$ & 3.90 & 2.66 & 2.39 & 2.28 \\
        \hline
    \end{tabular}
\end{minipage} &
\begin{minipage}[t]{.68\linewidth}
    \centering
    \captionof{table}{$\lambda$ and integrators choice with NFE=50}\label{tab:stochastic}
    \begin{tabular}{lllllll}
         \hline
          \multicolumn{7}{c}{FID at different $\lambda$}  \\
        Method & 0.0   & 0.1  & 0.3   & 0.5 & 0.7 & 1.0 \\
        \hline
        gDDIM    &  5.38 &  5.6 & 13.12 & 32  & 42 & 48\\
        EM       &  346  & 169  & 141   & 94  & 41 & 57\\
        \hline
    \end{tabular}
\end{minipage}
\end{tabular}

\textbf{Choice of $\lambda$: ODE vs SDE.} We further conduct comparison study with different $\lambda$ values. Note that polynomial extrapolation in \cref{eq:extrapolation} is not used here even when $\lambda=0$. As it is shown in~\cref{tab:stochastic}, increasing $\lambda$ deteriorates the sample quality, demonstrating our claim that determinstic DDIM has better performance than its stochastic counterpart when a small NFE is used.
We also find stochastic gDDIM significantly outperforms EM integrator.
\begin{table}[htbp]
    \caption{Unconditional CIFAR10 generative performance}
    \centering\begin{tabular}{llll}
        \hline
        Class & Model & NFE~($\downarrow$) & FID~($\downarrow$) \\
        \hline
        \multirow{5}{5em}{\textbf{CLD}}
        & our CLD-SGM~(gDDIM) & 50 & 2.28 \\
        & our CLD-SGM~(SDE, EM) & 2000 & 2.43 \\
        & our CLD-SGM~(Prob.Flow, RK45) & 157 & 3.16 \\
        \cline{2-4}
        & CLD-SGM~(Prob.Flow, RK45) by \citet{DocVahKre} & 147 & 2.71 \\
        & CLD-SGM~(SDE, EM) by \citet{DocVahKre} & 2000 & 2.23 \\
        \hline
        \multirow{3}{8em}{\textbf{Score}}
        & NCSN by \citet{song2019generative}& 1000 & 25.3 \\
        & DDPM by \citet{HoJaiAbb20}      & 1000 & 3.17 \\
        & D3PMs by \citet{austin2021structured}     & 1000 & 7.34 \\
        & DDPM++, VP by \citet{SonSohKin2020} & 2000 & 2.41 \\
        & DDPM++, sub-VP by \citet{SonSohKin2020} & 2000 & 2.41 \\
        & NCSN++ by \citet{SonSohKin2020} & 2000 & 2.20 \\
        & Likelihood SDE by \citet{song2021maximum} & - & 2.87 \\
        & Improved DDPM by \citet{Nichol2021ImprovedDD} & 4000 & 2.90 \\
        & UDM by \citet{kim2021score}                                & 2000 & 2.33 \\
        \hline
        \multirow{3}{8em}{\textbf{Score+Others}}
        & LSGM by \citet{vahdat2021score} & 138 & 2.10 \\
        & LSGM-100M by \citet{vahdat2021score} & 131 & 4.60 \\
        & Diffusion GAN by \citet{xiao2021tackling}&  4  & 3.75 \\
        \hline
    \end{tabular}
    \label{tab:cifar10}
    \vspace{-0.2cm}
\end{table}
\begin{table}[htbp]
    \caption{Efficiency comparison among training-free samplers, FID / NFE on CIFAR10
    }
    \centering\begin{tabular}{lllllll}
        \hline
        DDIM~\cite{SonMenErm21} &   FastDPM~\cite{kong2021fast}   &  DEIS~\cite{ZhaChe22} & PNDM~\cite{LiuRenLin22} & ADDIM~\cite{BaoLiZhu22}         & gDDIM    \\
        \hline
        13.36 / 10 &   9.90 / 10 &  \textbf{4.17 / 10}   &  7.05/19   &  14.0/10   & 13.41 / 10\\
        6.87  / 20 &   5.05 / 20 &  \textbf{3.31 / 20}   &  4.61/29   &  5.81/20   & 3.90 / 20\\
        4.47  / 50 &   3.32 / 50 &  2.99 / 50   &  3.68/59   &  4.04/50       & \textbf{2.87 / 27} \\
        4.16 / 100 &   2.86 /100 &  -           &  3.53/109  &  3.55/100      & \textbf{2.28 / 50} \\
        \hline
    \end{tabular}
    \label{tab:cifar10-low}
\end{table}
\section{Conclusion and Broader Impact}
\label{sec:con}

\textbf{Contribution}:
The more structural knowledge we take advantage, the more efficient we are.
In this work, we provide a clean interpretation of DDIMs based on the manifold hypothesis and the sparsity property on realistic dataset. This new perspective unboxes the numerical discretization used in DDIM and explains the advantage of ODE-based sampler over SDE-based when NFE is small.
Based on this interpretation, we extend DDIMs to general diffusion models.
The new algorithm, gDDIM, only conducts a tiny but elegant modification on the parameterization of score model and improves in sampling efficiency drastically.
We conduct extensive experiments to validate the effectiveness of our new sampling algorithm and achieve encouraging results when apply gDDIM in CLD models.

\textbf{Limitation}:
There are several promising future directions.
First, though gDDIM is designed for general DMs, we only verify it on CLD model in addition to DDPM.
gDDIM-based CLD shows better image sampling quality under $50$ NFE but still cannot produce satisfying results when NFE is extreme small, e.g. 10 or 5 steps. It is beneficial to explore more efficient diffusion process for different datasets, in which we believe gDDIM will play a important role in designing sampling algorithm.
Second, more investigation are needed to design efficient sampling algorithm by exploiting more structural knowledge in DMs.
The structural knowledge can originate from different sources such as different modality of datasets, and mathematical structure presented in specific diffusion processes. An interesting question for further improving CLD stochastic sampling is how to take advantage of the fact that gDDIM has a conservative Hamitonian component, in addition to noise, friction, and forcing from the score.

\newpage

\section{Proofs}
\label{app:derivation}
Since the proposed gDDIM is a generalization of DDIM, the results regarding gDDIM in~\cref{sec:gDDIM} are generalizations of those in \cref{sec:DDIM}.
In particular,
    \cref{prop:f-ddde} generalizes~\cref{prop:ddim-ode},
    \cref{eq:iddim-ode-eps} generalizes~\cref{prop:ddim-ode-update-step}, and
    \cref{prop:iddim-sde} generalizes~\cref{prop:ddim-sde}.
Thus, for the sake of simplicity, we mainly present proofs for gDDIM in~\cref{sec:gDDIM}.

\subsection{Deterministic gDDIM}

\subsubsection{Proof of \cref{eq:r}}

Since we assume data distribution $p_0({\bm{u}}) = {\mathcal{N}}({\bm{u}}_0, \Sigma_0)$, the score has closed form
\begin{equation}\label{eq:gt-score-gaussian}
    \nabla \log p_t({\bm{u}}) = -\Sigma_t^{-1}({\bm{u}} - \Psi(t,0){\bm{u}}_0).
\end{equation}

To make sure our construction~\cref{eq:usolution} is a solution to the probability flow ODE, we examine the condition for ${\bm{R}}_t$.
The LHS of the probability flow ODE is
\begin{align} \label{eq:r-lhs}
    d {\bm{u}} &= d [\Psi(t, 0){\bm{u}}_0 + {\bm{R}}_t {\epsilon}] \nonumber \\
          &= \dot{\Psi}(t, 0){\bm{u}}_0 dt + \dot{{\bm{R}}_t} {\epsilon} dt \nonumber \\
          &= [{\bm{F}}_t\Psi(t, 0){\bm{u}}_0 + \dot{{\bm{R}}_t} {\epsilon}] dt.
\end{align}
The RHS of the probability flow ODE is
\begin{align} \label{eq:r-rhs}
    [{\bm{F}}_t {\bm{u}} - \frac{1}{2} {\bm{G}}_t {\bm{G}}_t^T \nabla \log p_t({\bm{u}})] dt &= [{\bm{F}}_t\Psi(t, 0){\bm{u}}_0 + {\bm{F}}_t {\bm{R}}_t {\epsilon} + \frac{1}{2} {\bm{G}}_t {\bm{G}}_t^T {\bm{R}}_t^{-T} {\epsilon}]dt \nonumber \\
    &= [{\bm{F}}_t\Psi(t, 0){\bm{u}}_0 + {\bm{F}}_t{\bm{R}}_t{\epsilon}+ \frac{1}{2} {\bm{G}}_t {\bm{G}}_t^T {\bm{R}}_t^{-T} {\bm{R}}_t^{-1} {\bm{R}}_t {\epsilon}]dt,
\end{align}
where the first equality is due to $\nabla \log p_t({\bm{u}}) = -{\bm{R}}_t^{-T} {\epsilon}$.

Since \cref{eq:r-lhs,eq:r-rhs} holds for each ${\epsilon}$, we establish
\begin{align}
    \dot{{\bm{R}}}_t &= ({\bm{F}}_t+ \frac{1}{2} {\bm{G}}_t {\bm{G}}_t^T {\bm{R}}_t^{-T} {\bm{R}}_t^{-1}){\bm{R}}_t \\
                &= ({\bm{F}}_t+ \frac{1}{2} {\bm{G}}_t {\bm{G}}_t^T \Sigma_t^{-1}){\bm{R}}_t.
\end{align}

\subsubsection{Proof of \cref{prop:f-ddde}}

Similar to the proof of \cref{eq:r}, over a solution $\{{\bm{u}}(t), t\}$ to the probability flow ODE, ${\bm{R}}^{-1}({\bm{u}}(t) - \Psi(t,0){\bm{u}}_0)$ is constant. Furthermore, by~\cref{eq:gt-score-gaussian},
\begin{align}\label{eq:gt-score-gaussian-help1}
    \nabla \log p_t({\bm{u}}(t)) & = -\Sigma_t^{-1}({\bm{u}}(t) - \Psi(t,0){\bm{u}}_0) \nonumber \\
    & = -{\bm{R}}_t^{-T}{\bm{R}}_t^{-1}({\bm{u}}(t) - \Psi(t,0){\bm{u}}_0)
\end{align}
\cref{eq:gt-score-gaussian-help1} implies that $-{\bm{R}}_t^{T} \nabla \log p_t({\bm{u}}(t))$ is a constant and invariant with respect to $t$.

\subsubsection{Proof of \cref{eq:ODEiter,eq:iddim-ode-eps}}

We derive \cref{eq:iddim-ode-eps} first.

The update step is based on the approximation $\tilde{{\epsilon}}_\theta({\bm{u}}, \tau) = {\epsilon}_\theta({\bm{u}}(t), t)$ for $\tau \in [t-\Delta t, t]$.
The resulted ODE with $\tilde{{\epsilon}}_\theta$ reads
\begin{equation}
    \dot{{\bm{u}}} = {\bm{F}}_\tau {\bm{u}} + \frac{1}{2} {\bm{G}}_\tau {\bm{G}}_\tau^T {\bm{R}}_\tau^{-1} {\epsilon}_\theta({\bm{u}}(t),t),
\end{equation}
which is a linear ODE. The closed-form solution reads
\begin{equation}
    {\bm{u}}(t_{i-1}) = \Psi(t_{i-1}, t_{i}) {\bm{u}}(t)
    +
    [
        \int_{t_i}^{t_{i-1}} \frac{1}{2} \Psi(t_{i-1}, \tau) {\bm{G}}_\tau {\bm{G}}^T_\tau{\bm{R}}_\tau^{-T}
    ]{\epsilon}_\theta({\bm{u}}(t), t),
\end{equation}
where $\Psi(t,s)$ is the transition matrix associated ${\bm{F}}_t$, that is, $\Psi$ satisfies
\begin{equation}
    \frac{d \Psi(t,s)}{dt} = {\bm{F}}_t \Psi(t,s) \quad \Psi(s,s) = {\bm{I}}_D.
\end{equation}

When the DM is specified to be DDPM, we derive~\cref{eq:ODEiter} based on~\cref{eq:iddim-ode-eps} by expanding the coefficients in~\cref{eq:iddim-ode-eps} explicitly as
\begin{align*}
    \Psi(t,s) &= \sqrt{\frac{\alpha_t}{\alpha_s}}, \nonumber \\
    \Psi(t-\Delta t, t) &= \sqrt{\frac{\alpha_{t-\Delta t}}{\alpha_t}}, \\
    \int_t^{t-\Delta t} \frac{1}{2} \Psi(t-\Delta t, \tau) {\bm{G}}_\tau {\bm{G}}_\tau^{T} {\bm{R}}_{\tau}^{-1} d \tau
    &= \int_t^{t-\Delta t} -\frac{1}{2}
    \sqrt{\frac{\alpha_{t-\Delta t}}{\alpha_\tau}}
    \frac{d \log \alpha_\tau}{d \tau}
    \frac{1}{\sqrt{1 - \alpha_\tau}} d\tau \nonumber \\
    &= \sqrt{\alpha_{t-\Delta t}} \sqrt{\frac{1-\alpha_\tau}{\alpha_\tau}} \Bigg|_{\alpha_{t}}^{\alpha_{t-\Delta t}}  \nonumber \\
    &= \sqrt{1 - \alpha_{t-\Delta t}} - \sqrt{1 - \alpha_{t}}
    \sqrt{
   \frac{\alpha_{t-\Delta t}}{\alpha_t}
    }.
\end{align*}

\subsubsection{Proof of Multistep Predictor-Corrector}

Our Multistep Predictor-Corrector method slightly extends the traditional linear multistep Predictor-Corrector method to incorporate the semilinear structure in the probability flow ODE with an exponential integrator~\cite{press2007numerical,hochbruck2010exponential}.

\textbf{Predictor:}

For~\cref{eq:ode}, the key insight of the multistep predictor is to use existing function evaluations ${\epsilon}_\theta({\bm{u}}(t_{i}), t_{i}), {\epsilon}_\theta({\bm{u}}(t_{i+1}), t_{i+1}), \cdots, {\epsilon}_\theta({\bm{u}}(t_{i+q-1}), t_{i+q-1})$ and their timestamps $t_{i}, t_{i+1}, \cdots, t_{i+q-1}$ to fit a $q-1$ order polynomial ${}^{p}{\epsilon}_\theta(t)$ to approximate ${\epsilon}_\theta({\bm{u}}(\tau),\tau)$. With this approximator $\tilde{{\epsilon}}_\theta({\bm{u}}, \tau) = {}^{p}{\epsilon}_\theta(\tau)$ for $\tau \in [t_{i-1}, t_{i}]$, the multistep predictor step is obtained by solving
\begin{align}\label{eq:multistep-predictor-key}
    \frac{d {\bm{u}}}{dt} &= {\bm{F}}_t {\bm{u}} + \frac{1}{2} {\bm{G}}_\tau {\bm{G}}_\tau^T {\bm{R}}_\tau^{-T} \tilde{{\epsilon}}_\theta({\bm{u}}, \tau) \nonumber \\
    & = {\bm{F}}_t {\bm{u}} + \frac{1}{2} {\bm{G}}_\tau {\bm{G}}_\tau^T {\bm{R}}_\tau^{-T} {}^{p}{\epsilon}_\theta(\tau),
\end{align}
which is a linear ODE. The solution to~\cref{eq:multistep-predictor-key} satisfies
\begin{equation}\label{eq:multistep-predictor-sol}
    {\bm{u}}(t_{i-1}) = \Psi(t_{i-1}, t_{i}) {\bm{u}}(t_i) +
    \int_{t_i}^{t_{i-1}} \frac{1}{2}{\bm{G}}_\tau {\bm{G}}_\tau^T {\bm{R}}_\tau^{-T} {}^{p}{\epsilon}_\theta(\tau) d\tau.
\end{equation}
Based on Lagrange formula, we can write ${}^{p}{\epsilon}_\theta(\tau)$ as
\begin{equation}\label{eq:poly}
    {}^{p}{\epsilon}_\theta(\tau) = \sum_{j=0}^{q-1} [\prod_{k \neq j}\frac{\tau - t_{i+k}}{t_{i+j} - t_{i+k}}] {\epsilon}_\theta({\bm{u}}_{t_{i+j}}, t_{i+j}).
\end{equation}
Plugging~\cref{eq:poly} into~\cref{eq:multistep-predictor-sol}, we obtain
\begin{align}
    {\bm{u}}(t_{i-1}) &= \Psi(t_{i-1}, t_i) {\bm{u}}(t_i) + \sum_{j=0}^{q-1}[{}^{p}{\textnormal{C}}^{(q)}_{ij} {\epsilon}_\theta({\bm{u}}(t_{i+j}), t_{i+j})], \label{eq:app-gddim-predictor-update} \\
    {}^{p}{\textnormal{C}}^{(q)}_{ij} &= \!\int_{t_i}^{t_{i-1}}  \frac{1}{2} \Psi(t_{i-1}, \tau)
    {\bm{G}}_\tau {\bm{G}}_\tau^T{\bm{R}}_\tau^{-T}
    \prod_{k \neq j,k=0}^{q-1} [\frac{\tau - t_{i+k}}{t_{i+j} - t_{i+k}}]
    d\tau, \label{eq:app-gddim-predictor-coef}
\end{align}
which are ~\cref{eq:gddim-predictor-update,eq:gddim-predictor-coef}.
Here we use ${}^{p}{\textnormal{C}}^{(q)}_{ij}$ to emphasis these are constants used in the $q$-step predictor. The 1-step predictor reduces to~\cref{eq:iddim-ode-eps}.

\textbf{Corrector:}

Compared with the explicit scheme for the multistep predictor, the multistep corrector  behaves like an implicit method~\cite{press2007numerical}. Instead of constructing ${}^{p}{\epsilon}_\theta(\tau)$ to extrapolate model output for $\tau \in [t_{i-1}, t_{i}]$ as in the predictor, the $q$ step corrector aims to find ${}^{c}{\epsilon}_\theta(\tau)$ to interpolate ${\epsilon}_\theta({\bm{u}}(t_{i-1}), t_{i-1}), {\epsilon}_\theta({\bm{u}}(t_{i}), t_{i}), {\epsilon}_\theta({\bm{u}}(t_{i+1}), t_{i+1}), \cdots, {\epsilon}_\theta({\bm{u}}(t_{i+q-2}), t_{i+q-2})$ and their timestamps $t_{i-1}, t_{i}, t_{i+1}, \cdots, t_{i+q-2}$. Thus, ${\bm{u}}(t_{i-1})$ is obtained by solving
\begin{equation}
    {\bm{u}}(t_{i-1}) = \Psi(t_{i-1}, t_{i}) {\bm{u}}(t_{i}) +
    \int_{t_i}^{t_{i-1}} \frac{1}{2}{\bm{G}}_\tau {\bm{G}}_\tau^T {\bm{R}}_\tau^{-T} {}^{c}{\epsilon}_\theta(\tau) d\tau.
\end{equation}
Since ${}^{c}{\epsilon}_\theta(\tau)$ is defined implicitly, it is not easy to find ${}^{c}{\epsilon}_\theta(\tau), {\bm{u}}(t_{i-1})$. Instead, practitioners bypass the difficulties by interpolating ${{\epsilon}}_\theta(\bar{{\bm{u}}}(t_{i-1}), t_{i-1}), {\epsilon}_\theta(\bar{{\bm{u}}}(t_{i}), t_{i}), {\epsilon}_\theta(\bar{{\bm{u}}}(t_{i+1}), t_{i+1}), \cdots, {\epsilon}_\theta(\bar{{\bm{u}}}(t_{i+q-2}), t_{i+q-2})$ where $\bar{{\bm{u}}}(t_{i-1})$ is obtained by the predictor in~\cref{eq:multistep-predictor-sol} and $\bar{{\bm{u}}}(t_{i})={\bm{u}}({t_{i}}),\bar{{\bm{u}}}(t_{i+1})={\bm{u}}({t_{i+1}}),\cdots,\bar{{\bm{u}}}(t_{i+q-2})={\bm{u}}({t_{i+q-2}})$. Hence, we derive the update step for corrector based on
\begin{equation}\label{eq:multistep-corrector-sol}
    {\bm{u}}(t_{i-1}) = \Psi(t_{i-1}, t_{i}) {\bm{u}}(t_{i}) +
    \int_{t_i}^{t_{i-1}} \frac{1}{2}{\bm{G}}_\tau {\bm{G}}_\tau^T {\bm{R}}_\tau^{-T} {}^{c}{\epsilon}_\theta(\tau) d\tau,
\end{equation}
where ${}^{c}{\epsilon}_\theta(\tau)$ is defined as
\begin{equation}\label{eq:corrector-poly}
    {}^{c}{\epsilon}_\theta(\tau) = \sum_{j=-1}^{q-2} [\prod_{k \neq j}\frac{\tau - t_{i+k}}{t_{i+j} - t_{i+k}}] {\epsilon}_\theta(\bar{{\bm{u}}}_{t_{i+j}}, t_{i+j}).
\end{equation}

Plugging \cref{eq:corrector-poly} into \cref{eq:multistep-corrector-sol}, we reach the update step for the corrector
\begin{align}
    {\bm{u}}(t_{i-1}) &= \Psi(t_{i-1}, t_i) {\bm{u}}(t_i) + \sum_{j=-1}^{q-2}[{}^{c}{\textnormal{C}}^{(q)}_{ij} {\epsilon}_\theta(\bar{{\bm{u}}}(t_{i+j}), t_{i+j})], \label{eq:app-gddim-corrector-update} \\
    {}^{c}{\textnormal{C}}^{(q)}_{ij} &= \!\int_{t_i}^{t_{i-1}}  \frac{1}{2} \Psi(t_{i-1}, \tau)
    {\bm{G}}_\tau {\bm{G}}_\tau^T{\bm{R}}_\tau^{-T}
    \prod_{k \neq j,k=-1}^{q-2} [\frac{\tau - t_{i+k}}{t_{i+j} - t_{i+k}}]
    d\tau. \label{eq:app-gddim-corrector-coef}
\end{align}
We use ${}^{c}{\textnormal{C}}^{(q)}_{ij}$ to emphasis constants used in the $q$-step corrector.

\textbf{Exponential multistep Predictor-Corrector:}

Here we present the Exponential multistep Predictor-Corrector algorithm. Specifically, we employ one $q$-step corrector update step after an update step of the $q$-step predictor. The interested reader can easily extend the idea to employ multiple update steps of corrector or different number of steps for the predictor and the corrector. We note coefficients ${}^{p}{\textnormal{C}}, {}^{c}{\textnormal{C}}$ can be calculated using high resolution ODE solver once and used everywhere.

\begin{algorithm}[H]
    \caption{Exponential multistep Predictor-Corrector}
    \begin{algorithmic}
    \label{alg:empc}
    \STATE \textbf{Input}: Timestamps $\{t_i\}_{i=0}^N$, step order $q$, coefficients for predictor update ${}^{p}{\textnormal{C}}$, coefficients for corrector update ${}^{c}{\textnormal{C}}$
    \STATE \textbf{Instantiate}: ${\bm{u}}({t_N}) \sim p_T({\bm{u}})$
    \FOR{$i$ in $N, N-1, \cdots, 1$}
        \STATE \commet{predictor update step}
        \STATE $q_{\text{cur}} = \min(q, N-i+1)$ \commet{handle warming start, use lower order multistep method}
        \STATE ${\bm{u}}_{t_{i-1}} \leftarrow $ Simulate \cref{eq:app-gddim-predictor-update} with $q_{\text{cur}}$-step predictor
        \STATE \commet{corrector update step}
        \STATE $q_{\text{cur}} = \min(q, N-i+2)$ \commet{handle warming start, use lower order multistep method}
        \STATE $\bar{{\bm{u}}}(t_{i-1}),\bar{{\bm{u}}}(t_{i}), \cdots, \bar{{\bm{u}}}(t_{i+q_{\text{cur}}-1})  \leftarrow $
        ${{\bm{u}}}(t_{i-1}),{{\bm{u}}}(t_{i}), \cdots, {{\bm{u}}}(t_{i+q_{\text{cur}}-1})$
        \STATE ${\bm{u}}_{t_{i-1}} \leftarrow $ Simulate \cref{eq:app-gddim-corrector-update} with $q_{\text{cur}}$-step corrector
    \ENDFOR
    \end{algorithmic}
\end{algorithm}

\subsection{Stochastic gDDIM}

\subsubsection{Proof of \cref{prop:iddim-sde}}

Assuming that the data distribution $p_0({\bm{u}})$ is ${\mathcal{N}}({\bm{u}}_0, \Sigma_0)$ with a given $\Sigma_0$, we can derive the mean and covariance of $p_t({\bm{u}})$ as
\begin{align}
    \mu_t &= \Psi(t, 0) \\
    \frac{d \Sigma_t}{d t} &= {\bm{F}}_t \Sigma_t + \Sigma_t {\bm{F}}_t^{T} + {\bm{G}}_t {\bm{G}}_t^T.
\end{align}
Therefore, the ground truth score reads
\begin{equation}\label{eq:app-gt-gaussian-score}
    \nabla \log p_t({\bm{u}}) = -\Sigma_t^{-1} ({\bm{u}} - \Psi(t,0){\bm{u}}_0).
\end{equation}

We assume $\Sigma_0$ is given but ${\bm{u}}_0$ is unknown. Fortunately, ${\bm{u}}_0$ can be inferred via one score evaluation as follows.
Given evaluation $\nabla \log p_s({\bm{u}}(s))$, we can recover ${\bm{u}}_0$ as
\begin{equation}\label{eq:app-recover-vu0}
    {\bm{u}}_0 = \Psi(0, s)[\Sigma_s \nabla \log p_s({\bm{u}}(s)) + {\bm{u}}(s)].
\end{equation}
Plugging \cref{eq:app-recover-vu0} and $\Psi(t,s)=\Psi(t,0)\Psi(0,s)$ into~\cref{eq:app-gt-gaussian-score}, we recover~\cref{eq:iddim-sde-approx-gt}.

\subsubsection{Proof of \cref{prop:s-gddim-update-step}}
With the approximator $\tilde{{\epsilon}}_\theta({\bm{u}}, t)$ defined in~\cref{eq:iddim-sde-approx} for $\tau \in [s, t]$, \cref{eq:lambda-sde} can be reformulated as
\begin{align}\label{eq:s-gddim-sde}
    d {\bm{u}} &= {\bm{F}}_\tau {\bm{u}} d\tau + \frac{1+\lambda^2}{2} {\bm{G}}_\tau {\bm{G}}_\tau^T {\bm{R}}_\tau^{-T} \tilde{{\epsilon}}_\theta({\bm{u}}, \tau)dt + \lambda {\bm{G}}_\tau d{\bm{w}} \nonumber \\
    & = ({\bm{F}}_\tau + \frac{1+\lambda^2}{2} {\bm{G}}_\tau {\bm{G}}_\tau^T {\bm{R}}_\tau^{-T} {\bm{R}}_\tau^{-1}) {\bm{u}} dt \nonumber \\
    & \quad + \frac{1+\lambda^2}{2} {\bm{G}}_\tau {\bm{G}}_\tau^T {\bm{R}}_\tau^{-T} {\bm{R}}_\tau^{-1}\Psi(\tau, s)({\bm{R}}_s {\epsilon}_\theta({\bm{u}}(s), s) - {\bm{u}}(s))dt + \lambda {\bm{G}}_\tau d{\bm{w}}.
\end{align}
Define $\hat{{\bm{F}}}_\tau := {\bm{F}}_\tau + \frac{1+\lambda^2}{2} {\bm{G}}_\tau {\bm{G}}_\tau^T {\bm{R}}_\tau^{-T} {\bm{R}}_\tau^{-1} = {\bm{F}}_\tau + \frac{1+\lambda^2}{2} {\bm{G}}_\tau {\bm{G}}_\tau^T \Sigma^{-1}_\tau$,
and denote by $\hat{\Psi}(t,s)$ the transition matrix associated with it.
Clearly, \cref{eq:s-gddim-sde} is a linear differential equation on ${\bm{u}}$, and the conditional probability $\hat{p}_{st}({\bm{u}}(t)|{\bm{u}}(s))$ associated with it is a Gaussian distribution.

Applying ~\citet[Eq~(6.6,6.7)]{SarSol19}, we obtain the exact expressions
\begin{align}\label{eq:s-ddim-mean-sol}
    \text{Mean}\, =\,\, &\hat{\Psi}(t,s){\bm{u}}(s)
    - [\int_s^t \hat{\Psi}(t, \tau)\frac{1+\lambda^2}{2} {\bm{G}}_\tau {\bm{G}}_\tau^T \Sigma^{-1}_\tau \Psi(\tau, s)]d\tau\, {\bm{u}}(s) \nonumber \\
    & +[\int_s^t \hat{\Psi}(t, \tau)\frac{1+\lambda^2}{2} {\bm{G}}_\tau {\bm{G}}_\tau^T \Sigma^{-1}_\tau \Psi(\tau, s)]d\tau \,{\bm{R}}_s {\epsilon}_\theta({\bm{u}}(s), s)
\end{align}
for the mean of $\hat{p}_{st}({\bm{u}}(t)|{\bm{u}}(s))$. Its covariance ${\bm{P}}_{s\tau}$ satisfies
\begin{equation}
    \frac{d {\bm{P}}_{s\tau}}{d \tau} = \hat{{\bm{F}}_\tau} {\bm{P}}_{s\tau}
    + {\bm{P}}_{s\tau} \hat{{\bm{F}}_\tau}^{T} + \lambda^2 {\bm{G}}_\tau {\bm{G}}_\tau^T, \quad {\bm{P}}_{ss} = 0.
\end{equation}

\cref{eq:s-ddim-mean-sol} has a closed form expression with the help of the following lemma.
\begin{lemma}\label{lemma:s-gddim}
    $$\int_s^t \hat{\Psi}(t, \tau)\frac{1+\lambda^2}{2} {\bm{G}}_\tau {\bm{G}}_\tau^T \Sigma^{-1}_\tau \Psi(\tau, s) = \hat{\Psi}(t, s) - {\Psi}(t, s)$$
\end{lemma}
\begin{proof}
   For a fixed $s$, we define ${\bm{N}}(t) = \int_s^t \hat{\Psi}(t, \tau)\frac{1+\lambda^2}{2} {\bm{G}}_\tau {\bm{G}}_\tau^T \Sigma^{-1}_\tau \Psi(\tau, s)$ and ${\bm{M}}(t) = \hat{\Psi}(t, s) - {\Psi}(t, s)$. It follows that
    \begin{align}
        \frac{d {\bm{N}}(\tau)}{d\tau} &= \hat{{\bm{F}}}_\tau {\bm{N}}(\tau) + \frac{1+\lambda^2}{2} {\bm{G}}_\tau {\bm{G}}_\tau^T \Sigma^{-1}_\tau \Psi(\tau, s) \nonumber \\
        & = {\bm{F}}_\tau {\bm{N}}(\tau) + \frac{1+\lambda^2}{2} {\bm{G}}_\tau {\bm{G}}_\tau^T \Sigma^{-1}_\tau [{\bm{N}}(\tau) + \Psi(\tau, s)] \\
        \frac{d {\bm{M}}(\tau)}{d\tau} &= \hat{{\bm{F}}}_\tau \hat{\Psi}(\tau, s) - {\bm{F}}_\tau \Psi(\tau, s) \nonumber \\
        &= {\bm{F}}_\tau {\bm{M}}(\tau) + \frac{1+\lambda^2}{2} {\bm{G}}_\tau {\bm{G}}_\tau^T \Sigma^{-1}_\tau \hat{\Psi}(\tau, s).
    \end{align}
Define ${\bm{E}}(t)={\bm{N}}(t) - {\bm{M}}(t)$, then
    \begin{equation}
         \frac{d {\bm{E}}(t)}{dt} = ({\bm{F}}_t + \frac{1+\lambda^2}{2} {\bm{G}}_t {\bm{G}}_t^T \Sigma_t^{-1}){\bm{E}}(t).
    \end{equation}
On the other hand, ${\bm{N}}(s) = {\bm{M}}(s) = 0$ which implies ${\bm{E}}(s)=0$. We thus conclude ${\bm{E}}(t)=0$ and ${\bm{N}}(t) = {\bm{M}}(t)$.
\end{proof}
Using~\cref{lemma:s-gddim}, we simplify~\cref{eq:s-ddim-mean-sol} to
\begin{equation}
    \Psi(t, s){\bm{u}}(s) + [\hat{\Psi}(t, s) - \Psi(t,s)]{\bm{R}}_s{\epsilon}_\theta({\bm{u}}(s),s),
\end{equation}
which is the mean in~\cref{eq:iddim-sde-udpate-step}.

\subsubsection{Proof of \cref{thm:ddims}}

We restate conclusion presented in \cref{thm:ddims}. The exact solution ${\bm{u}}({t-\Delta t})$ to~\cref{eq:lambda-sde} with coefficient \cref{eq:DDPM} is
\begin{equation}
    {\bm{u}}({t - \Delta t}) \sim
    {\mathcal{N}}(
        \sqrt{\frac{\alpha_{t - \Delta t}}{\alpha_t}} {\bm{u}}(t) +
        \left[
            - \sqrt{\frac{\alpha_{t - \Delta t}}{\alpha_t}} \sqrt{1-\alpha_t} +
            \sqrt{1-\alpha_{t - \Delta t} - \sigma_t^2}
        \right] {\epsilon}_\theta({\bm{u}}(t), t),
        \sigma_t^2 {\bm{I}}_d
    )
\end{equation}
with $\sigma^2_t = (1 -\alpha_{t - \Delta t})\left[
    1 - \left(
        \frac{1-\alpha_{t - \Delta t}}{1-\alpha_t}
    \right)^{\lambda^2}
    \left(
        \frac{\alpha_t}{\alpha_{t - \Delta t}}
    \right)^{\lambda^2}
\right]$, which is the same as the DDIM \cref{eq:ddim}.

\cref{thm:ddims} is a concrete application of~\cref{eq:iddim-sde-udpate-step} when the DM is a DDPM and ${\bm{F}}_\tau, {\bm{G}}_\tau$ are set to~\cref{eq:DDPM}. Thanks to the special form of ${\bm{F}}_\tau$, $\Psi$ has the expression
\begin{equation}
    \Psi(t, s) = \sqrt{\frac{\alpha_t}{\alpha_s}} {\bm{I}}_d,
\end{equation}
and $\hat{\Psi}$ satisfies
\begin{align}
    \log \hat{\Psi}(t, s) &= \int_s^t [
        \frac{1}{2} \frac{d \log \alpha_\tau}{d \tau}
        - \frac{1 + \lambda^2}{2}\frac{d \log \alpha_\tau}{d \tau} \frac{1}{1 - \alpha_\tau}
    ] d\tau \\
    \hat{\Psi}(t, s) & = \left(
        \frac{1 - \alpha_t}{1-\alpha_s}
    \right)^{\frac{1+\lambda^2}{2}}
    \left(
        \frac{\alpha_s}{\alpha_t}
    \right)^{\frac{\lambda^2}{2}}.
\end{align}

\textbf{Mean:}
Based on~\cref{eq:iddim-sde-udpate-step}, we obtain the mean of $\hat{p}_{st}$ as
\begin{align}
    & \sqrt{\frac{\alpha_t}{\alpha_s}} {\bm{u}}(s) +
    \left[
        - \sqrt{\frac{\alpha_t}{\alpha_s}} \sqrt{1-\alpha_s} +
        \left(
            \frac{1 - \alpha_t}{1-\alpha_s}
        \right)^{\frac{1+\lambda^2}{2}}
        \left(
            \frac{\alpha_s}{\alpha_t}
        \right)^{\frac{\lambda^2}{2}}\sqrt{1-\alpha_s}
    \right]
    {\epsilon}_\theta({\bm{u}}(s), s) \\
    & = \sqrt{\frac{\alpha_t}{\alpha_s}} {\bm{u}}(s) +
    \left[
        - \sqrt{\frac{\alpha_t}{\alpha_s}} \sqrt{1-\alpha_s} +
        \sqrt{
            (1 - \alpha_t)
            \left(
                \frac{1-\alpha_{t - \Delta t}}{1-\alpha_t}
            \right)^{\lambda^2}
            \left(
                \frac{\alpha_t}{\alpha_{t - \Delta t}}
            \right)^{\lambda^2}
        }
    \right]
    {\epsilon}_\theta({\bm{u}}(s), s) \\
    & = \sqrt{\frac{\alpha_t}{\alpha_s}} {\bm{u}}(s) +
    \left[
        - \sqrt{\frac{\alpha_t}{\alpha_s}} \sqrt{1-\alpha_s} +
        \sqrt{
            1 - \alpha_t - \sigma_s^2
        }
    \right]
    {\epsilon}_\theta({\bm{u}}(s), s),
\end{align}
where
\begin{equation}
    \sigma_s^2 = (1 -\alpha_{t})\left[
        1 - \left(
            \frac{1-\alpha_{t}}{1-\alpha_s}
        \right)^{\lambda^2}
        \left(
            \frac{\alpha_s}{\alpha_{t}}
        \right)^{\lambda^2}
    \right].
\end{equation}
Setting $(s, t) \leftarrow (t, t-\Delta t$), we arrive at the mean update in~\cref{eq:lambda-ddim}.

\textbf{Covariance:}
It follows from
\[
    \frac{d {\bm{P}}_{s\tau}}{d\tau} =
    2[
        \frac{d \log \alpha_\tau}{2 d \tau} -
        \frac{1+\lambda^2}{2}
        \frac{d \log \alpha_\tau}{d \tau}
        \frac{1}{1 - \alpha_\tau}
    ] {\bm{P}}_{s\tau}
    - \lambda^2 \frac{d \log \alpha_\tau}{d \tau} {\bm{I}}_d,~ {\bm{P}}_{ss} = 0
    \]
that
\[
    {\bm{P}}_{st} = (1 -\alpha_t)\left[
        1 - \left(
            \frac{1-\alpha_t}{1-\alpha_s}
        \right)^{\lambda^2}
        \left(
            \frac{\alpha_s}{\alpha_t}
        \right)^{\lambda^2}
    \right].
\]
Setting $(s, t) \leftarrow (t, t-\Delta t)$, we recover the covariance in~\cref{eq:lambda-ddim}.

\subsubsection{Proof of \cref{prop:s-gddim-update-lambda0}}

When $\lambda=0$, the update step in~\cref{eq:iddim-sde-udpate-step} from $s$ to $t$ reads
\begin{align}\label{eq:lambda0-s-gddim}
    {\bm{u}}(t) &= \Psi(t, s){\bm{u}}(s) + [\hat{\Psi}(t, s) - \Psi(t, s)] {\bm{R}}_s {\epsilon}_\theta({\bm{u}}(s), s).
\end{align}
Meanwhile, the update step in~\cref{eq:iddim-ode-eps} from $s$ to $t$ is
\begin{equation}\label{eq:appx-d-gddim}
    {\bm{u}}(t) = \Psi(t, s){\bm{u}}(s) + [
        \int_{s}^t \frac{1}{2}\Psi(t, \tau) {\bm{G}}_\tau {\bm{G}}_\tau^{T}{\bm{R}}_\tau^{-T} d\tau
    ] {\epsilon}_\theta({\bm{u}}(s), s).
\end{equation}

\cref{eq:appx-d-gddim,eq:lambda0-s-gddim} are equivalent once we have the following lemma.
\begin{lemma}
    When $\lambda=0$, $$\int_{s}^t \frac{1}{2}\Psi(t, \tau) {\bm{G}}_\tau {\bm{G}}_\tau^{T}{\bm{R}}_\tau^{-T} d\tau= [\hat{\Psi}(t, s) - \Psi(t, s)] {\bm{R}}_s.$$
\end{lemma}
\begin{proof}
    We introduce two new functions
    \begin{align}
        {\bm{N}}(t) &:= \int_{s}^t \frac{1}{2}\Psi(t, \tau) {\bm{G}}_\tau {\bm{G}}_\tau^{T}{\bm{R}}_\tau^{-T} d\tau \\
        {\bm{M}}(t) &:= [\hat{\Psi}(t, s) - \Psi(t, s)] {\bm{R}}_s.
    \end{align}
    First, ${\bm{N}}(s) = {\bm{M}}(s) = 0$.
    Second, they satisfy
    \begin{align}
        \frac{d {\bm{N}}(t)}{d t} &= {\bm{F}}_t {\bm{N}}(t) + \frac{1}{2} {\bm{G}}_t {\bm{G}}_t^{T} {\bm{R}}_t^{-T} \\
        \frac{d {\bm{M}}(t)}{d t} &= [\hat{{\bm{F}}}_t \hat{\Psi}(t, s) - {\bm{F}}_t \Psi(t, s)] {\bm{R}}_s \\
        &={\bm{F}}_t {\bm{M}}(t) + \frac{1}{2}{\bm{G}}_t {\bm{G}}_t^{T} \Sigma_t^{-1} \hat{\Psi}(t, s) {\bm{R}}_s \label{eq:lambda0-s-gddim-m}.
    \end{align}
Note $\hat{\Psi}$ and ${\bm{R}}$ satisfy the same linear differential equation as
    \begin{equation}
        \frac{d \hat{\Psi}(t,s)}{d t} = [{\bm{F}}_t + \frac{1}{2}{\bm{G}}_t {\bm{G}}_t^T \Sigma_t^{-1}] \hat{\Psi}(t,s),
        \quad
        \frac{d {\bm{R}}_t}{dt} = [{\bm{F}}_t + \frac{1}{2}{\bm{G}}_t {\bm{G}}_t^T \Sigma_t^{-1}] {\bm{R}}_t.
    \end{equation}
It is a standard result in linear system theory (see \citet[Eq(2.34)]{SarSol19}) that $\hat{\Psi}(t,s) = {\bm{R}}_t {\bm{R}}^{-1}_s$. Plugging it and ${\bm{R}}_t {\bm{R}}_t^T = \Sigma_t$ into \cref{eq:lambda0-s-gddim-m} yields
    \begin{equation}
        \frac{d {\bm{M}}(t)}{d t} = {\bm{F}}_t {\bm{M}}(t) + \frac{1}{2} {\bm{G}}_t {\bm{G}}_t^{T} {\bm{R}}_t^{-T}.
    \end{equation}
Define ${\bm{E}}(t) = {\bm{N}}(t) - {\bm{M}}(t)$, then it satisfies
    \begin{equation}
        {\bm{E}}(s) = 0 \quad \frac{d{\bm{E}}(t)}{d t} = {\bm{F}}_t {\bm{E}}(t),
    \end{equation}
    which clearly implies that ${\bm{E}}(t)=0$. Thus, ${\bm{N}}(t) = {\bm{M}}(t)$.
\end{proof}
\section{More related works}
\label{app:related}

Learning generative models with DMs via score matching has received tremendous attention recently~\cite{sohl2015deep,lyu2012interpretation, song2019generative,SonSohKin20,HoJaiAbb20,Nichol2021ImprovedDD}.
However, the sampling efficiency of DMs is still not satisfying.
\citet{jolicoeur2021gotta} introduced adaptive solver for SDEs associated with DMs for the task of image generation.
\citet{SonMenErm21} modified the forward noising process into a non-Markov process without changing the training objective function. The authors then proposed a family of samplers, including deterministic DDIM and stochastic DDIM, based on the modifications. Both of the samplers demonstrate significant improvements over previous samplers. There are variants of the DDIM that aim to further improve the sampling quality and efficiency.
\citet{BaoLiZhu22} improved the DDIM by optimizing the reverse variance in DMs.
\citet{WatChaHo22} generalized the DDIM in DDPM with learned update coefficients, which is trained by minimizing an external perceptual loss.
\citet{Nichol2021ImprovedDD} tuned the variance of schedule of DDPM.
\citet{LiuRenLin22} found that the DDIM is a pseudo numerical method and proposed a pseudo linear multi-step method for it.
\citet{ZhaChe22} discovered that DDIMs are numerical integrators for marginal-equivalent SDEs, and the deterministic DDIM is actually an exponential integrator for the probability flow ODE.
They further utilized exponential multistep methods to boost sampling performance for VPSDE.

Another promising approach to accelerate diffusion model is distillation for the probability flow ODE.
\citet{luhman2021knowledge} proposed to learn the map from noise to data in a teacher-student fashion, where supervised signals are provided by simulating the deterministic DDIM.
The final student network distilled is able to generate samples with reasonable quality within one step.
\citet{Salimans2022} proposed a new progressive distillation approach to improve the training efficiency and stability. This distillation approach relies on solving the probability flow ODE and needs extra training procedure. Since we generalize and improve the DDIM in this work, it will be beneficial to combine this distillation method with our algorithm for better performance in the future.

Recently, numerical structures of DMs have received more and more attentions; they play important roles in efficient sampling methods. \citet{DocVahKre} designed Symmetric Splitting CLD Sampler~(SSCS) that takes advantage of Hamiltonian structure of the CLD and demonstrated advantages over the naive Euler-Maruyama method. \citet{ZhaChe22} utilized the semilinear structure presented in DMs and showed that the exponential integrator gave much better sampling quality than Euler method.
Further tailoring these integrators to account for the stiff property of ODEs~\cite{hochbruck2010exponential, Whalen2015} is a promising direction in fast sampling for DMs.

\section{More experiment details}
\label{app:exp-impl}

We present practical implementation of gDDIM and its application to CLD. We include training details for ${\bm{R}}_t$-based CLD proposed in this work and discuss necessary calculation overhead for executing gDDIM.
More experiments are conducted to verify the effectiveness of gDDIM compared with other sampling algorithms.
We report image sampling performance over average of 3 runs with different random seeds.

\subsection{Training CLD}
Our training pipeline, model architectures and hyperparameters are similar to those in~\citet{DocVahKre}. The main differences are in the choice of ${\bm{K}}_t$ and loss weights ${{\bm{K}}_t^{-1}\Lambda_t{\bm{K}}_t^{-T}}$.

We change from ${\bm{K}}_t={\bm{L}}_t$ to ${\bm{K}}_t={\bm{R}}_t$. Unlike ${\bm{L}}_t$ which has triangular structure and closed form expression, we rely on high accuracy numerical solver to solve ${\bm{R}}_t$.
The triangular structure of ${\bm{L}}_t$ and sparse pattern of ${\bm{G}}_t$ for CLD in \cref{eq:cld-fg} also have an impact on the training loss function of the score model.
Due to special structure of ${\bm{G}}_t$,
we only need to learn ${\bm{s}}_\theta({\bm{u}}, t;{\bm{v}})$ signals present in the velocity channel.
When ${\bm{K}}_t = {\bm{L}}_t$, ${\bm{s}}_\theta({\bm{u}}, t) = -{\bm{L}}^{-T}_t {\epsilon}^{({\bm{L}}_t)}_\theta({\bm{u}}, t)$ with ${\bm{L}}^{-T}_t$ being upper triangular. Thus we only need to train ${\epsilon}^{({\bm{L}}_t)}_\theta({\bm{u}}, t;{\bm{v}})$ to recover ${\bm{s}}_\theta({\bm{u}}, t;{\bm{v}})$.
In contrast, ${\bm{R}}_t$ does not share the triangular structure as ${\bm{L}}_t$;
both ${\epsilon}^{({\bm{R}}_t)}_\theta({\bm{u}}, t;{\bm{x}}),{\epsilon}^{({\bm{R}}_t)}_\theta({\bm{u}}, t;{\bm{v}})$ are needed to recover ${\bm{s}}_\theta({\bm{u}}, t;{\bm{v}})$.
Therefore,~\citet{DocVahKre} sets loss weights
\begin{equation}
    {\bm{L}}_t^{-1} \Lambda_t {\bm{L}}_t^{-T} = \begin{bmatrix}
        0     & 0 \\
        0     & 1
        \end{bmatrix}
     \otimes {\bm{I}}_d,
\end{equation}
while we choose
\begin{equation}
    {\bm{R}}_t^{-1} \Lambda_t {\bm{R}}_t^{-T} = \begin{bmatrix}
        1     & 0 \\
        0     & 1
        \end{bmatrix}
     \otimes {\bm{I}}_d.
\end{equation}

As a result, we need to double the channel in the output layer in our new parameterization associated with ${\bm{R}}_t$, though the increased number of parameters in last layer is negligible compared with other parts of
diffusion models.

We include the model architectures and hyperparameters in \cref{tab:model}. In additional to the standard size model on CIFAR10, we also train a smaller model for CELEBA to show the efficacy and advantages of gDDIM.
\begin{table}[htbp]
    \caption{Model architectures and hyperparameters}
    \label{tab:model}
    \centering
    \begin{tabular}{lll}
        \hline
        Hyperparameter                  & CIFAR10                   &        CELEBA \\
        \hline
        Model                           &                           &               \\
        EMA rate                        &  0.9999                   &         0.999 \\
        \# of ResBlock per resolution   &  8                        &         2     \\
        Normalization                   &  Group Normalization      & Group Normalization \\
        Progressive input               &  Residual                 &    None \\
        Progressive combine             &  Sum                      &    N/A  \\
        Finite Impluse Response         &  Enabled                  &   Disabled \\
        Embedding type                  & Fourier                   &   Positional\\
        \# of arameters                 & $\approx$ 108M            & $\approx$ 62M \\
        \hline
        Training                        &                           &           \\
        \# of iterations                & 1m                      &   150k    \\
        Optimizer                       & Adam                      &   Adam    \\
        Learning rate                   & 2$\times 10^{-4}$         & 2$\times 10^{-4}$ \\
        Gradient norm clipping          & 1.0                       & 1.0       \\
        Dropout                         & 0.1                       & 0.1       \\
        Batch size per GPU              & 32                        & 32 \\
        GPUs                            & 4 A6000                   & 4 A6000 \\
        Trainig time                    &  $\approx$ 78h         & $\approx$ 16h   \\
        \hline
    \end{tabular}
\end{table}

\subsection{Calculation of constant coefficients}

In gDDIM, many coefficients cannot be obtained in closed-form.
Here we present our approach to obtain them numerically. Those constant coefficients can be divided into two categories, solutions to ODEs and definite integrals. We remark that these coefficients only need to be calculated once and then can be used everywhere. For CLD,
each of these coefficient corresponds to a
$2\times2$ matrix. The calculation of all these coefficients can be done within 1 min.

\textbf{Type I: Solving ODEs}
The problem appears when we need to evaluate ${\bm{R}}_t$ in \cref{eq:r} and $\hat{\Psi}(t,s)$ in
\begin{equation}
    \frac{d{\hat{\Psi}(t,s)}}{dt} = \hat{{\bm{F}}}_t \hat{\Psi}(t,s), \quad \hat{\Psi}(s,s) = {\bm{I}}.
\end{equation}
Across our experiments, we use RK4 with a step size $10^{-6}$ to calculate the ODE solutions. For $\hat{\Psi}(t,s)$, we only need to calculate $\hat{\Psi}(t,0)$ because $\hat{\Psi}(t,s) = \hat{\Psi}(t,0) [\hat{\Psi}(s,0)]^{-1}$. In CLD, ${\bm{F}}_t, {\bm{G}}_t, {\bm{R}}_t, \Sigma_t$ can be simplified to a $2\times 2$ matrix; solving the ODE with a small step size is extremely fast. We note $\Psi(t,s)$ and $\Sigma_t$ admit close-form formula~\cite{DocVahKre}.
Since the output of numerical solvers are discrete in time, we employ a linear interpolation to handle query in continuous time.
Since ${\bm{R}}_t, \hat{\Psi}$ are determined by the forward SDE in DMs, the numerical results can be shared.
In stochastic gDDIM~\cref{eq:iddim-sde-udpate-step}, we apply the same techniques to solve ${\bm{P}}_{st}$.

\textbf{Type II: Definite integrals}

The problem appears in the derivation of update step in \cref{eq:iddim-ode-eps,eq:gddim-predictor-coef,eq:app-gddim-corrector-coef}, which require coefficients such as ${}^{p}{\textnormal{C}}^{(q)}_{ij}, {}^{c}{\textnormal{C}}^{(q)}_{ij}$. We use step size $10^{-5}$ for the integration from $t$ to $t - \Delta t$.
The integrand can be efficiently evaluated in parallel using GPUs. Again, the coefficients, such as ${}^{p}{\textnormal{C}}^{(q)}_{ij}, {}^{c}{\textnormal{C}}^{(q)}_{ij}$, are calculated once and used afterwards if we need sample another batch with the same time discretization.

\subsection{More experiments on the choice of score parameterization}

Here we present more experiments details and more experiments regarding \cref{prop:ddim-ode} and differences between the two parameterizations involving ${\bm{R}}_t$ and ${\bm{L}}_t$.

\textbf{Toy experiments: }

\begin{figure}[htbp]
    \centering
    \includegraphics[width=0.5\linewidth]{figure/ddpm_1d.pdf}
    \caption{Trajectory and ${\epsilon}$ of Probability Flow solution in DDPM}
    \label{fig:app_ddim_1d}
\end{figure}

\begin{figure}[htbp]
    \centering
    \includegraphics[width=\linewidth]{figure/cld_1d.pdf}
    \caption{Trajectory and ${\epsilon}$ of Probability Flow solution in CLD}
    \label{fig:app_gddim_1d}
\end{figure}

Here we present more empirical results to demonstrate the advantage of proper ${\bm{K}}_t$.

In VPSDE, the advantage of DDIM has been verified in various models and datasets~\cite{SonMenErm21}. To empirically verify~\cref{prop:ddim-ode}, we present one toy example where the data distribution is a mixture of two one dimension Gaussian distributions. In this toy example, the ground truth $\nabla \log p_t({\bm{u}})$ is known. As it shown in~\cref{fig:app_ddim_1d}, along the solution to probability flow ODE, the score parameterization $-{\bm{R}}^{T}_t \log p_t({\bm{x}})$ enjoys better smoothing property.

In CLD, we present a similar study. As the covariance in CLD is no longer diagonal, we find the difference of ${\bm{L}}_t$ parameterization and ${\bm{R}}_t$ parameterization is large in~\cref{fig:app_gddim_1d}.

\textbf{Image experiments:}

\begin{table}
    \centering
    \begin{tabular}{ll|llll}
        \hline
         && \multicolumn{3}{c}{FID / IS at different NFE}  \\
       $q$ & ${\bm{K}}_t$ & 20 & 30 & 40 & 50 \\
       \hline
       \multirow{2}{*}{0}
           & ${\bm{L}}_t$ & 461.72 / 1.18   & 441.05 / 1.21   & 244.26 / 2.39   & 120.07 / 5.25 \\
           & ${\bm{R}}_t$ & 17.94 / 8.47   & 9.84 / 9.00   & 6.87 / 9.19   & 5.46 / 9.31 \\
       \hline
       \multirow{2}{*}{1}
           & ${\bm{L}}_t$ & 368.38 / 1.32   & 232.39 / 2.43   & 99.45 / 4.92   & 57.06 / 6.70 \\
           & ${\bm{R}}_t$ & 8.14 / 9.10   & 4.35 / 9.46   & 3.32 / 9.58   & 2.92 / 9.65 \\
       \hline
       \multirow{2}{*}{2}
           & ${\bm{L}}_t$ & 463.33 / 1.17   & 166.90 / 3.56   & 4.12 / 9.25   & 3.31 / 9.38 \\
           & ${\bm{R}}_t$ & \textbf{4.54 / 9.42}   & \textbf{2.66 / 9.64}   & \textbf{2.39 / 9.74}   & 2.32 / 9.75 \\
       \hline
       \multirow{2}{*}{3}
           & ${\bm{L}}_t$ & 464.36 / 1.17   & 463.45 / 1.17   & 463.32 / 1.17   & 240.45 / 2.29 \\
           & ${\bm{R}}_t$ & 341.70 / 1.42   & 292.31 / 1.70   & 13.27 / 10.15   & \textbf{2.28 / 9.77} \\
       \hline
    \end{tabular}
    \caption{More experiments on CIFAR10}
    \label{tab:app_r_vs_l_img_cifar}
\end{table}

\begin{table}
    \centering
    \begin{tabular}{ll|llll}
        \hline
         && \multicolumn{3}{c}{FID at different NFE}  \\
       $q$ & ${\bm{K}}_t$ & 20 & 30 & 40 & 50 \\
       \hline
       \multirow{2}{*}{0}
           & ${\bm{L}}_t$ & 446.56   & 430.59   & 434.79   & 379.73 \\
           & ${\bm{R}}_t$ & 37.72   & 13.65   & 12.51   & 8.96 \\
       \hline
       \multirow{2}{*}{1}
           & ${\bm{L}}_t$ & 261.90   & 123.49   & 105.75   & 88.31   \\
           & ${\bm{R}}_t$ & 12.68   & 7.78   & 5.93   & 5.11 \\
       \hline
       \multirow{2}{*}{2}
           & ${\bm{L}}_t$ & 446.74   & 277.28   & 6.18   & 5.48  \\
           & ${\bm{R}}_t$ & \textbf{6.86}   & \textbf{5.67}   & \textbf{4.62}   & 4.19 \\
       \hline
       \multirow{2}{*}{3}
           & ${\bm{L}}_t$ & 449.21   & 440.84   & 443.91   & 286.22 \\
           & ${\bm{R}}_t$ & 386.14   & 349.48   & 20.14   & \textbf{3.85} \\
       \hline
    \end{tabular}
    \caption{More experiments on CELEBA}
    \label{tab:app_r_vs_l_img_celeba}
\end{table}

We present more empirical results regarding the comparison between ${\bm{L}}_t$ and ${\bm{R}}_t$. Note that we use exponential integrator for the parametrization ${\bm{L}}_t$ as well, similar to DEIS~\cite{ZhaChe22}. We vary the polynomial order $q$ in multistep methods~\cite{ZhaChe22} and test sampling performance on CIFAR10 and CELEBA. In both datasets, we generate 50$k$ images and calcualte their FID. As shown in~\cref{tab:app_r_vs_l_img_cifar,tab:app_r_vs_l_img_celeba}, ${\bm{R}}_t$ has significant advantages, especially when NFE is small.

We also find that multistep method with a large $q$ can harm sampling performance when NFE is small. This is reasonable; the method with larger $q$ assumes the nonlinear function is smooth in a large domain and may rely on outdated information for approximation, which may worsen the accuracy.

\subsection{More experiments on the choice of $\lambda$}

\begin{figure}[htbp]
    \centering
    \includegraphics[width=1.0\textwidth]{figure/ddpm_lambda_traj.pdf}
    \caption{Sampling with various $\lambda$ and accurate score on a Toy example. Trajectories are shown from $T$ to $0$.}
    \label{fig:oned_lambda_toy_example}
\end{figure}

To study the effects of $\lambda$, we visualize the trajectories generated with various $\lambda$ but the same random seeds in~\cref{fig:oned_lambda_toy_example} on our toy example. Clearly, trajectories with smaller $\lambda$ have better smoothing property while trajectories with large $\lambda$ contain much more randomness. From the fast sampling perspective, trajectories with more stochasticity are much harder to predict with small NFE compared with smooth trajectories.

We include more qualitative results on the choice of $\lambda$ and comparison between the Euler-Maruyama (EM) method and the gDDIM in~\cref{fig:app_cifar_lambda,fig:app_celeba_lambda}. Clearly, when NFE is small, increasing $\lambda$ has a negative effect on the sampling quality of gDDIM. We hypothesize that $\lambda=0$ already generates high-fidelity samples and additional noise may harm the sampling performance. With a fixed number of function evaluations, information derived from score network fails to remove the injected noise as we increase $\lambda$. On the other hand, we find that the EM method shows slightly better quality as we increase $\lambda$. We hypothesize that the ODE or SDEs with small $\lambda$ has more oscillations than SDEs with large $\lambda$. It is known that the EM method has a very bad performance for oscillations systems and suffers from large discretization error~\cite{press2007numerical}. From previous experiments, we find that ODE in CLD is highly oscillated.

We also find both methods perform worse than Symmetric Splitting CLD Sampler~(SSCS)~\cite{DocVahKre} when $\lambda=1$.
The improvement by utilizing Hamiltonian structure and SDEs structure is significant. This  encourages further exploration that incorporates Hamiltonian structure into gDDIM in the future. Nevertheless, we also remark that SSCS with $\lambda=1.0$ performs much worse than gDDIM with $\lambda=0$.

\subsection{More comparisons}

We also compare the performance of the CLD model we trained with that claimed in~\citet{DocVahKre} in~\cref{tab:app_comp}. We find that our trained model performs worse than \citet{DocVahKre} when a blackbox ODE solver or EM sampling scheme with large NFE are used.
There may be two reasons. First, with similar size model, our training scheme not only needs to fit $\nabla_{{\bm{v}}} \log p_t({\bm{u}})$, but also $\nabla_{{\bm{x}}} \log p_t({\bm{u}})$, while~\citet{DocVahKre} can allocate all representation resources of neural network to $\nabla_{{\bm{v}}} \log p_t({\bm{u}})$. Another factor is the mixed score trick on parameterization, which is shown empirically have a boost in model performance~\cite{DocVahKre} but we do not include it in our training.

We also compare our algorithm with more accelerating sampling methods in~\cref{tab:app_comp}. gDDIM has achieved the best sampling acceleration results among training-free methods, but it still cannot compete with some distillation-based acceleration methods.

\begin{table}[htbp]
    \caption{More comparison on CIFAR10}
    \centering
    \begin{tabular}{llll}
        \hline
        Class & Model & NFE~($\downarrow$) & FID~($\downarrow$) \\
        \hline
        \multirow{5}{5em}{\textbf{CLD}}
        & our CLD-SGM~(gDDIM)      & 50   & 2.28 \\
        & our CLD-SGM~(SDE, EM)               & 2000 & 2.43 \\
        & our CLD-SGM~(Prob.Flow, RK45)       & 157  & 3.16 \\
        \cline{2-4}
        & CLD-SGM~(Prob.Flow, RK45) by \citet{DocVahKre} & 147 & 2.71 \\
        & CLD-SGM~(SDE, EM) by \citet{DocVahKre} & 2000 & 2.23 \\
        \hline
        \multirow{3}{8em}{\textbf{Training-free Accelerated Score}}
        & DDIM by \citet{SonMenErm21} & 100 & 4.16 \\
        & Gotta Go Fast by \citet{JolLiPic21}    & 151  & 2.73 \\
        & Analytic-DPM by \citet{BaoLiZhu22}  & 100 & 3.55 \\
        & FastDPM by \citet{KonPin21}         & 100 & 2.86 \\
        & PNDM by \citet{LiuRenLin22} & 100 & 3.53 \\
        & DEIS by \citet{ZhaChe22}      & 50  & 2.99 \\
        \hline
        \multirow{3}{8em}{\textbf{Training-needed Accelerated Score}}
        & DDSS by \citet{WatChaHo22} & 25 & 4.25 \\
        & Progressive Distillation by \citet{Salimans2022} & 4 & 3.0 \\
        & Knowledge distillation by \citet{luhman2021knowledge} & 1 & 9.36 \\
        \hline
    \end{tabular}
    \label{tab:app_comp}
    \vspace{-0.2cm}
\end{table}

\subsection{Code licenses}

We implemented gDDIM and related algorithms in Jax. We have used code from a number of sources in \cref{tab:license}.
\begin{table}[h]
    \centering
    \begin{tabular}{lll}
    \hline
     URL & Citation & License \\
    \hline
    \url{https://github.com/yang-song/score_sde} & \cite{SonSohKin2020} &  Apache License 2.0  \\
    \url{https://github.com/nv-tlabs/CLD-SGM} & \cite{DocVahKre} &  NVIDIA Source Code License \\
    \url{https://github.com/qsh-zh/deis} & \cite{ZhaChe22} & Unknown \\
    \hline
    \end{tabular}
    \caption{Code License}
    \label{tab:license}
\end{table}

\begin{figure}[htbp]
    \centering
    \includegraphics[width=0.80\textwidth]{figure/cifar_L.pdf}
    \includegraphics[width=0.80\textwidth]{figure/cifar_R.pdf}
    \caption{Comparison between ${\bm{L}}$~(Upper) and ${\bm{R}}$~(Lower) with exponential integrator on CIFAR10.}
\end{figure}

\begin{figure}[htbp]
    \centering
    \includegraphics[width=0.80\textwidth]{figure/celeba_L.pdf}
    \includegraphics[width=0.80\textwidth]{figure/celeba_R.pdf}
    \caption{Comparison between ${\bm{L}}$~(Upper) and ${\bm{R}}$~(Lower) with exponential integrator on CELEBA.}
\end{figure}

\begin{figure}[htbp]
    \centering
    \includegraphics[width=0.80\textwidth]{figure/cifar_em.pdf}
    \includegraphics[width=0.80\textwidth]{figure/cifar_sdeis.pdf}
    \caption{Comparison between EM~(Upper) and gDDIM~(Lower) on CIFAR10.}
    \label{fig:app_cifar_lambda}
\end{figure}

\begin{figure}[htbp]
    \centering
    \includegraphics[width=0.80\textwidth]{figure/celeba_em.pdf}
    \includegraphics[width=0.80\textwidth]{figure/celeba_sdeis.pdf}
    \caption{Comparison between EM~(Upper) and gDDIM~(Lower) on CELEBA.}
    \label{fig:app_celeba_lambda}
\end{figure}