(Continuous Grönwall Lemma) Let $I=[a,b]$ denote a compact interval of the
real line with $a<b$. Let $\alpha,\beta,u$ be continuous real-valued functions
defined on $I$. Assume $\beta$ is non-negative, $\alpha$ is non-decreasing and
$u$ satisfies the integral inequality

$u(t)\leq\alpha(t)+\int_{a}^{t}\beta(s)u(s)\hskip
1.69998pt\mathrm{d}s,\qquad\forall t\in I.$

Then

$u(t)\leq\alpha(t)\exp\left(\int_{a}^{t}\beta(s)\hskip
1.69998pt\mathrm{d}s\right),\qquad\forall t\in I.$

######  Lemma B.11.

Assume the following forward SDE is satisfied
$dX_{t}=F(t,X_{t})dt+G(t)dW_{t},t\in[0,T]$, where $T>0$, $F(t,x)$ is Lipschitz
with respect to $(t,x)$, $G$ is continuous and $X_{0}$ is an integrable random
variable. Then, there exists $C>0$ such that for all $s<t\in[0,T]$ with
$t-s\leq 1$,

$\mathbb{E}[|X_{t}-X_{s}|]\leq C\sqrt{t-s}.$

###### Proof.

We take $s=0$ and apply the triangular inequality:

$|X_{t}-X_{0}|\leq\int_{0}^{t}|F(u,X_{u})-F(0,X_{0}))|\mathrm{d}u+|\int_{0}^{t}G(u)dW_{u}|+t|F(0,X_{0})|.$

Setting $u(t)=\mathbb{E}[|X_{t}-X_{0}|]$ and taking the expectation, we
deduce:

$u(t)\leq
K\frac{t^{2}}{2}+t\mathbb{E}[|F(0,X_{0})|]+\mathbb{E}[(\int_{0}^{t}G(u)dW_{u})^{2}]^{\frac{1}{2}}+K\int_{0}^{t}u(s)ds,$

where $K$ is a positive constant. Note that
$\mathbb{E}[(\int_{0}^{t}G(u)dW_{u})^{2}]^{\frac{1}{2}}=[\int_{0}^{t}G^{2}(s)ds]^{\frac{1}{2}}$
by the Itô isometry property which is less than $C\sqrt{t}$ since $G$ is
continuous. Thus, we have proved that:

$u(t)\leq\alpha(t)+K\int_{0}^{t}u(s)ds,$

where $\alpha(t)=K\frac{t^{2}}{2}+t\mathbb{E}[|F(0,X_{0})|]+C\sqrt{t}$ is non-
decreasing. By Lemma B.10: $u(t)\leq\alpha(t)\exp(Kt)\leq\alpha(t)\exp(KT)$.
Since $\alpha(t)\leq C\sqrt{t}$ for $t\in[0,1]$, the lemma holds. ∎

######  Lemma B.12.

(Discrete Grönwall Lemma) Consider a real number sequence $(u_{n})_{n}$ such
that

$u_{n+1}\leq a_{n+1}u_{n}+b_{n+1},n\geq 0$

where $(a_{n})$ and $(b_{n})_{n}$ are two given sequences such that $b_{n}$ is
positive. Then

$u_{n}\leq A_{n}u_{0}+\sum_{k=1}^{n}A_{k,n}b_{k},$

where $A_{n}=\prod_{k=1}^{n}a_{k}$ and
$A_{k,n}=A_{n}/A_{k}=\prod_{i=k+1}^{n}a_{i}$.

###  B.3 Proof of Corollary 4.3

####  B.3.1 Convergence of SEEDS-2:

Let $\\{\mathbf{y}_{t_{i}}\\}_{i}$ be the discrete stochastic process defined
as follows:

$\mathbf{y}_{t_{0}}\leftarrow\mathbf{x}_{T},\mathbf{y}_{t_{i}}\leftarrow\frac{\alpha_{t_{i}}}{\alpha_{t_{i-1}}}\mathbf{y}_{t_{i-1}}-2\alpha_{t_{i}}\int_{\lambda_{t_{i-1}}}^{\lambda_{t_{i}}}e^{-u}\bm{\epsilon}(\bm{u}_{i},s_{i})\mathrm{d}u-\sqrt{2}\alpha_{t_{i}}\int_{\lambda_{t_{i-1}}}^{\lambda_{t_{i}}}e^{-s}\mathrm{d}\bm{\bar{\omega}}(s),$

with $s_{i}\leftarrow
t_{\lambda}\left(\lambda_{t_{i-1}}+\frac{h_{i}}{2}\right)$ and

$\bm{u}_{i}\leftarrow\frac{\alpha_{s_{i}}}{\alpha_{t_{i-1}}}\mathbf{y}_{t_{i-1}}-2\sigma_{s_{i}}\left(e^{\frac{h_{i}}{2}}-1\right)\bm{\epsilon}(\mathbf{y}_{t_{i-1}},t_{i-1})-\sqrt{2}\alpha_{s_{i}}\int_{\lambda_{t_{i-1}}}^{\lambda_{s_{i}}}e^{-s}\mathrm{d}\bm{\bar{\omega}}(s),$

Then $\\{\mathbf{y}_{t_{i}}\\}_{i}$ has the same distribution as
$\\{\widetilde{\mathbf{x}}_{t_{i}}\\}_{i}$ in Algorithm 3. We can compute the
difference:

$\mathbf{y}_{t_{i}}-\mathbf{x}_{t_{i}}=\frac{\alpha_{t_{i}}}{\alpha_{t_{i-1}}}(\mathbf{y}_{t_{i-1}}-\mathbf{x}_{t_{i-1}})+\Gamma;\Gamma=\Gamma_{1}+\Gamma_{2}$

and:

$\displaystyle\Gamma_{1}$ $\displaystyle=$
$\displaystyle-2\alpha_{t_{i}}\int_{\lambda_{t_{i-1}}}^{\lambda_{t_{i}}}e^{-u}\left[\bm{\epsilon}(\bm{u}_{i},s_{i})-\bm{\epsilon}(x_{s_{i}},s_{i})\right]\mathrm{d}u$
$\displaystyle\Gamma_{2}$ $\displaystyle=$
$\displaystyle-2\alpha_{t_{i}}\int_{\lambda_{t_{i-1}}}^{\lambda_{t_{i}}}e^{-u}\left[\bm{\epsilon}(x_{s_{i}},s_{i})-\widehat{\bm{\epsilon}}(\widehat{\mathbf{x}}_{u},u)\right]\mathrm{d}u$

Similarly to the case $k=1$, $\mathbb{E}[|\Gamma_{2}|]\leq Ch^{2}$. Note that
$\mathbb{E}[|\Gamma_{1}|]\leq Ch\Delta_{i}$ with
$\Delta_{i}=\mathbb{E}[|\bm{u}_{i}-\mathbf{x}_{s_{i}}|]$. Introduce:
$U_{i}=\mathbb{E}[|\mathbf{y}_{t_{i}}-\mathbf{x}_{t_{i}}|]$. We will bound
$\Delta_{i}$ by a function of $U_{i-1}$. For this, recall:

$\bm{u}_{i}=\frac{\alpha_{s_{i}}}{\alpha_{t_{i-1}}}\mathbf{y}_{t_{i-1}}-2\alpha_{s_{i}}\int_{\lambda_{t_{i-1}}}^{\lambda_{s_{i}}}e^{-\lambda}\bm{\epsilon}(\mathbf{y}_{t_{i-1}},t_{i-1})\mathrm{d}\lambda-\sqrt{2}\alpha_{s_{i}}\int_{\lambda_{t_{i-1}}}^{\lambda_{s_{i}}}e^{-\lambda}\mathrm{d}\bm{\bar{\omega}}(\lambda),$

and write the difference:

$\bm{u}_{i}-\mathbf{x}_{s_{i}}=\frac{\alpha_{s_{i}}}{\alpha_{t_{i-1}}}(\mathbf{y}_{t_{i-1}}-\mathbf{x}_{t_{i-1}})+2\alpha_{t_{i}}\int_{\lambda_{t_{i-1}}}^{\lambda_{t_{i}}}e^{-u}(\widehat{\bm{\epsilon}}(\widehat{\mathbf{x}}_{u},u)-\bm{\epsilon}_{\theta}(\mathbf{y}_{t_{i-1}},\lambda_{t_{i-1}}))\mathrm{d}\lambda.$

By the triangular inequality:

$|\widehat{\bm{\epsilon}}(\widehat{\mathbf{x}}_{u},u)-\bm{\epsilon}(\mathbf{y}_{t_{i-1}},\lambda_{t_{i-1}})|\leq
C(|\mathbf{x}_{t_{\lambda}(u)}-\mathbf{x}_{t_{i-1}}|+|\mathbf{x}_{t_{i-1}}-\mathbf{y}_{t_{i-1}}|+|t_{\lambda}(u)-t_{i-1}|),$

and so:
$\mathbb{E}[|\widehat{\bm{\epsilon}}(\widehat{\mathbf{x}}_{u},u)-\bm{\epsilon}(\mathbf{y}_{t_{i-1}},\lambda_{t_{i-1}})|]\leq
C(\sqrt{h}+U_{i-1}+h)$. Finally:

$\Delta_{i}\leq\frac{\alpha_{s_{i}}}{\alpha_{t_{i-1}}}U_{i-1}+Ch\left(\sqrt{h}+U_{i-1}+h\right)\leq
U_{i-1}+Ch\left(\sqrt{h}+U_{i-1}\right),$

and $\mathbb{E}[|\Gamma_{1}|]\leq Ch^{5/2}+C(h+h^{2})U_{i-1}$. This gives the
bound:

$U_{i}\leq Ch^{2}+(\frac{\alpha_{t_{i}}}{\alpha_{t_{i-1}}}+Ch+Ch^{2})U_{i-1}.$

Now using Theorem B.6, the proof can now be finished following as in proof in
Section B.2.4.

####  B.3.2 Convergence of SEEDS-3:

Continuing with the same notations as before, we will prove by analogy that:

$U_{i}\leq
Ch^{2}+(\frac{\alpha_{t_{i}}}{\alpha_{t_{i-1}}}+Ch+Ch^{2}+Ch^{3})U_{i-1},$
(49)

so that, by Theorem B.6, we obtain the desired result.

Let $\mathbf{y}_{t_{i}}=\widetilde{\mathbf{x}}_{t_{i}}$. We have:

$y_{t_{i}}=\frac{\alpha_{t_{i}}}{\alpha_{t_{i-1}}}y_{t_{i-1}}-2\alpha_{t_{i}}\int_{\lambda_{t_{i-1}}}^{\lambda_{t_{i}}}e^{-\lambda}\bm{\epsilon}(y_{t_{i-1}},t_{i-1})d\lambda$
$-\frac{2\alpha_{t_{i}}}{h_{i}r_{2}}\int_{\lambda_{t_{i-1}}}^{\lambda_{t_{i}}}e^{-\lambda}(\lambda-\lambda_{t_{i-1}})(\bm{\epsilon}(u_{2i},s_{2i})-\bm{\epsilon}(y_{t_{i-1}},t_{i-1}))d\lambda+\text{Noise},$

and
$x_{t_{i}}=\frac{\alpha_{t_{i}}}{\alpha_{t_{i-1}}}x_{t_{i-1}}-2\alpha_{t_{i}}\int_{\lambda_{t_{i-1}}}^{\lambda_{t_{i}}}e^{-\lambda}\widehat{\bm{\epsilon}}(\widehat{\mathbf{x}}_{\lambda},\lambda)\mathrm{d}\lambda+\text{Noise}$,
where Noise is the same in both equations. From the proof for $k=1$, it
suffices to bound:

$\Gamma=\frac{1}{h_{i}}\int_{\lambda_{t_{i-1}}}^{\lambda_{t_{i}}}e^{-\lambda}(\lambda-\lambda_{t_{i-1}})(\bm{\epsilon}(u_{2i},s_{2i})-\bm{\epsilon}(y_{t_{i-1}},t_{i-1}))d\lambda,$

in the $L^{1}$ norm. By the Lipschitz property of $\bm{\epsilon}$:

$\displaystyle\mathbb{E}[|\Gamma|]$ $\displaystyle\leq$ $\displaystyle
Ch(\mathbb{E}[|u_{2i}-x_{s_{2i}}|]+\mathbb{E}[|x_{t_{i-1}}-x_{s_{2i}}|]+\mathbb{E}[|y_{t_{i-1}}-x_{t_{i-1}}|])+Ch^{2}$
$\displaystyle\leq$ $\displaystyle
Ch(\mathbb{E}[|u_{2i}-x_{s_{2i}}|]+\sqrt{h}+U_{i-1})+Ch^{2}.$

Now, let us bound $\mathbb{E}[|u_{2i}-x_{s_{2i}}|]$ and for this, write

$u_{2i}=\frac{\alpha_{s_{2i}}}{\alpha_{t_{i-1}}}y_{t_{i-1}}-2\alpha_{s_{2i}}\int_{\lambda_{t_{i-1}}}^{\lambda_{s_{2i}}}e^{-\lambda}\bm{\epsilon}(y_{t_{i-1}},t_{i-1})\mathrm{d}\lambda$
$-2\alpha_{s_{2i}}\int_{\lambda_{t_{i-1}}}^{\lambda_{s_{2i}}}e^{-\lambda}(\lambda-\lambda_{t_{i-1}})\frac{\epsilon(u_{2i-1},s_{2i-1})-\epsilon(y_{t_{i-1}},t_{i-1})}{r_{1}h_{i}}\mathrm{d}\lambda+\text{Noise}$

and

$x_{s_{2i}}=\frac{\alpha_{s_{2i}}}{\alpha_{t_{i-1}}}x_{t_{i-1}}-2\alpha_{s_{2i}}\int_{\lambda_{t_{i-1}}}^{\lambda_{s_{2i}}}e^{-\lambda}\widehat{\bm{\epsilon}}(\widehat{\mathbf{x}}_{\lambda},\lambda)\mathrm{d}\lambda+\text{Noise},$

where Noise is the same in both equations. So

$\displaystyle u_{2i}-x_{s_{2i}}$ $\displaystyle=$
$\displaystyle\frac{\alpha_{s_{2i}}}{\alpha_{t_{i-1}}}(y_{t_{i-1}}-x_{t_{i-1}})-2\alpha_{s_{2i}}\int_{\lambda_{t_{i-1}}}^{\lambda_{s_{2i}}}e^{-\lambda}(\bm{\epsilon}(\mathbf{y}_{t_{i-1}},t_{i-1})-\widehat{\bm{\epsilon}}(\widehat{\mathbf{x}}_{\lambda},\lambda))\mathrm{d}\lambda$
$\displaystyle-2\alpha_{s_{2i}}\int_{\lambda_{t_{i-1}}}^{\lambda_{s_{2i}}}e^{-\lambda}(\lambda-\lambda_{t_{i-1}})D\mathrm{d}\lambda,$

where
$D=\frac{\epsilon(u_{2i-1},s_{2i-1})-\epsilon(y_{t_{i-1}},t_{i-1})}{r_{1}h_{i}}$.
From the convergence proof of SEEDS-1, and the Lipschitz property of
$\epsilon$, we obtain:

$\mathbb{E}[|u_{2i}-x_{s_{2i}}|]\leq
U_{i-1}+ChU_{i-1}+Ch^{2}+Ch\mathbb{E}[|u_{2i-1}-y_{t_{i-1}}|].$

Again, by the triangular inequality:

$\mathbb{E}[|u_{2i-1}-y_{t_{i-1}}|]\leq\mathbb{E}[|u_{2i-1}-x_{s_{2i-1}}|]+C\sqrt{h}+U_{i-1},$

and since

$u_{2i-1}-x_{s_{2i-1}}=\frac{\alpha_{s_{2i-1}}}{\alpha_{t_{i-1}}}(y_{t_{i-1}}-x_{t_{i-1}})-2\alpha_{s_{2i-1}}\int_{\lambda_{t_{i-1}}}^{\lambda_{s_{2i-1}}}e^{-\lambda}(\widehat{\bm{\epsilon}}(\widehat{\mathbf{x}}_{\lambda},\lambda)-\bm{\epsilon}(y_{t_{i-1}},t_{i-1})\mathrm{d}\lambda,$

we have, as before,

$\mathbb{E}[|u_{2i-1}-x_{s_{2i}}|]\leq U_{i-1}+ChU_{i-1}+Ch^{2}.$

Combining the previous inequalities leads to (49). This finishes the proof of
Corollary 4.3.

##  Appendix C Implementation Details

The SEEDS solvers used in our experiments are exactly the variants we
described in Algorithms 2, 3, 4 in the main part of the paper. In particular,
SEEDS-2 and SEEDS-3 solvers are completely determined by one-parameter
families of deterministic exponential integrators from [10] of order 2 and 3
respectively and prescribed by the following Butcher tableaux:

$\begin{array}[]{l|ll}0&&\\\ c_{2}&c_{2}\varphi_{1,2}&\\\
\hline\cr&\left(1-\frac{1}{2c_{2}}\right)\varphi_{1}&\frac{1}{2c_{2}}\varphi_{1}\end{array}\qquad\begin{array}[]{l|lll}0&0&0&0\\\
c_{2}&c_{2}\varphi_{1,2}&0&0\\\
\frac{2}{3}&\frac{2}{3}\varphi_{1,3}-\frac{4}{9c_{2}}\varphi_{2,3}&\frac{4}{9c_{2}}\varphi_{2,3}&0\\\
\hline\cr&\varphi_{1}-\frac{3}{2}\varphi_{2}&0&\frac{3}{2}\varphi_{2}\end{array}$

In all experiments we fix the parameter $c_{2}=0.5$ for SEEDS-2 and
$c_{2}=\frac{1}{3}$ for SEEDS-3. We point out that the Butcher tableau here
associated to SEEDS-3 is the result of a weakening of the stiff order
conditions so it might suffer from order reduction in the deterministic case.
We also point out that a full theory of stiff order conditions for stochastic
exponential Runge-Kutta methods for semi-linear DEs with a time-varying linear
coefficient have not been yet developed to the authors knowledge. As such,
this is only an analogy whose purpose is to clarify how our solvers relate to
well-known solvers from the literature, but such Butcher tableaux do not
rigorously reflect their convergence order.

###  C.1 Stabilization of the exponential terms

In the proposed algorithms, one subtle detail is to re-arrange all equations
in order for them to make use only of the $\texttt{expm1}(h)$ function which
computes $e^{h}-1$ with great numerical stability, specially for small values
of $h$. We have for values $r_{1}=1/3$, $r_{2}=2/3$, the following identity

$\displaystyle\sqrt{e^{2h_{i}}-e^{h_{i}}}z_{i}+\sqrt{e^{h_{i}}-1}v_{i}$
$\displaystyle=$
$\displaystyle\sqrt{e^{h_{i}}-1}\left(\sqrt{e^{h_{i}}}z_{i}+v_{i}\right)$
$\displaystyle=$
$\displaystyle\sqrt{e^{h_{i}}-1}\left(e^{\frac{h_{i}}{2}}z_{i}+v_{i}\right)$
$\displaystyle=$
$\displaystyle\sqrt{\texttt{expm1}(h_{i})}\left(\left(\texttt{expm1}\left(\frac{h_{i}}{2}\right)+1\right)z_{i}+v_{i}\right)$
$\displaystyle=$
$\displaystyle\sqrt{\texttt{expm1}(h_{i})}\left(\texttt{expm1}\left(\frac{h_{i}}{2}\right)+1\right)z_{i}+\sqrt{\texttt{expm1}(h_{i})}v_{i}.$

Now by using

$\displaystyle e^{2(r_{2}-r_{1})h_{i}}$ $\displaystyle=$ $\displaystyle
e^{2\left(\frac{2}{3}-\frac{1}{3}\right)h_{i}}=e^{2\frac{1}{3}h_{i}}=e^{r_{2}h_{i}},$

we now compute

$\displaystyle\sqrt{e^{2r_{2}h_{i}}-e^{r_{2}h_{i}}}z^{1}_{i}+\sqrt{e^{r_{2}h_{i}}-1}z^{2}_{i}$
$\displaystyle=$
$\displaystyle\sqrt{e^{r_{2}h_{i}}-1}\left(\sqrt{e^{r_{2}h_{i}}}z^{1}_{i}+z^{2}_{i}\right)$
$\displaystyle=$
$\displaystyle\sqrt{e^{r_{2}h_{i}}-1}\left(e^{\frac{h_{i}}{3}}z^{1}_{i}+z^{2}_{i}\right)$
$\displaystyle=$
$\displaystyle\sqrt{\texttt{expm1}(r_{2}h_{i})}((\texttt{expm1}(r_{1}h_{i})+1)z^{1}_{i}+z^{2}_{i})$
$\displaystyle\sqrt{e^{2h_{i}}-e^{2r_{2}h_{i}}}z^{1}_{i}$ $\displaystyle=$
$\displaystyle\sqrt{e^{3r_{2}h_{i}}-e^{2r_{2}h_{i}}}z^{1}_{i}$
$\displaystyle=$ $\displaystyle
e^{r_{2}h_{i}}\sqrt{e^{r_{2}h_{i}}-1}z^{1}_{i}$ $\displaystyle=$
$\displaystyle\sqrt{\texttt{expm1}(r_{2}h_{i})}(\texttt{expm1}(r_{2}h_{i})+1)z^{1}_{i}.$

Finally, we obtain

$\displaystyle\sqrt{e^{2h_{i}}-e^{2r_{2}h_{i}}}z^{1}_{i}+\sqrt{e^{2r_{2}h_{i}}-e^{r_{2}h_{i}}}z^{2}_{i}+\sqrt{e^{r_{2}h_{i}}-1}z^{3}_{i}$
$\displaystyle=$
$\displaystyle\sqrt{\texttt{expm1}(r_{2}h_{i})}((\texttt{expm1}(r_{2}h_{i})+1)z^{1}_{i}+(\texttt{expm1}(r_{1}h_{i})+1)z^{2}_{i}+z^{3}_{i}).$

###  C.2 Noise schedules parameterizations

The inverses of $\lambda(t)$ are given in the (VP) linear and cosine schedules
respectively by

$\displaystyle t_{\lambda}(\lambda)$ $\displaystyle=$
$\displaystyle\frac{2\log(e^{-2\lambda}+1)}{\sqrt{\beta_{\min}^{2}+2(\beta_{\max}-\beta_{\min})\log(e^{-2\lambda}+1)}+\beta_{\min}},$
$\displaystyle t_{\lambda}(\lambda)$ $\displaystyle=$
$\displaystyle\frac{2(1+s)}{\pi}\arccos\left(\text{exp}\Big{(}-\dfrac{1}{2}\log(e^{-2\lambda}+1)+\log\cos\big{(}\cfrac{\pi
s}{2(1+s)}\big{)}\Big{)}\right)-s.$

In EDM Noise Prediction case, the inverses of $\lambda(t)$, as given in (37),
with respect to the NPFO and NRSDE are respectively given by

$\displaystyle
t_{\lambda}(\lambda)=\sigma_{d}\tan(e^{-\lambda})\quad\text{and}\quad
t_{\lambda}(\lambda)=\dfrac{\sigma_{d}}{\sqrt{\frac{1}{\sigma^{2}_{d}e^{-2\lambda}}-1}}.$

###  C.3 EDM discretization

We follow Karras et al. [12] to implement the EDM discretization timesteps
$\\{t_{i}\\}_{i=0}^{M}$ as $t_{i}=\sigma^{-1}(\sigma_{i})$ such that, for
$\rho>0$,

$\displaystyle\sigma_{i<M}=\left(\sigma_{\text{max}}^{\frac{1}{\rho}}+\frac{i}{N-1}\left(\sigma_{\text{min}}^{\frac{1}{\rho}}-\sigma_{\text{max}}^{\frac{1}{\rho}}\right)\right)^{\rho}\quad\text{and}\quad\sigma_{N}=0.$

From the definition, we note that $\sigma_{0}=\sigma_{\text{max}}$ and
$\sigma_{M-1}=\sigma_{\text{min}}$, where $\sigma_{\text{min}}$ and
$\sigma_{\text{max}}$ denote the minimum and maximum noise magnitude,
respectively. We also keep default value $\rho=7$ as in [12]. However, we
figured out that when using EDM discretization with linear schedule, the noise
schedule improvement in iDDPM pre-trained models [19] would result in two
consecutive time steps of the same value, i.e. $t_{j}=t_{j+1}$ for some index
$j=1,\dots,M-1$ and large steps ($M\geq 61$). Thus, for SEEDS-3 and DPM-
Solver-3 [15], the usage of the function
$\varphi_{2}(h)=\dfrac{e^{h}-h-1}{h^{2}}$ (see Appendix D for more details)
will cause zero division error where $h=\lambda_{t_{j+1}}-\lambda_{t_{j}}=0$.
Therefore, in our implementation, we ignore the noise schedule improvement of
iDDPM [19] models when using solvers of order three.

###  C.4 Final sampling step

The sampling phase in DPMs using SEEDS follows the RSDE, which requires
gradual computing through discretization time steps $\\{t_{i}\\}_{i=0}^{M}$
and the latter goes from $t_{0}=T$ to $t_{M}=0$. In our implementation, to
avoid the logarithm of zero error at the last step, i.e.
$\log\sigma_{t_{M}}=\log(0)$, we stop the sampling phase at step $(M-1)$.
Hence, the NFE used in a run will be given as

$\displaystyle\text{NFE}=k\times(M-1),$

where $k$ represents the order of the solver. We also do not use the
“denoising” trick, i.e., ignoring the random noise at the last step and leave
it to further research.

##  Appendix D Reminders on Stochastic Exponential Integrators

Let us consider a SDE of the form

$\mathrm{d}\mathbf{x}(t)=[a(t)\mathbf{x}(t)+c(t)f(\mathbf{x}(t),t)]\mathrm{d}t+g(t)\mathrm{d}\bm{\omega}(t),$
(50)

where $a,c:[0,T]\rightarrow\mathbb{R}$ and
$g:[0,T]\rightarrow\mathbb{R}^{d\times d}$. In other words we concentrate to
high-dimensional semi-linear non autonomous SDEs with additive noise. The
objective of this section is to construct explicit stochastic exponential
derivative-free methods for the above equation following the Runge-Kutta (RK)
approach. These methods ideally should fulfill the following properties:

  1. 1.

If $f\equiv 0$, then (50) can be solved exactly;

  2. 2.

If $g\equiv 0$, then a SEEDS method for (50) identifies with an exponential RK
(ERK) method;

  3. 3.

If $a\equiv 0$, then a SEEDS method for (50) identifies with a stochastic RK
(SRK) method and if moreover $g\equiv 0$ then it identifies with a classical
RK ODE method.

Before tackling the aimed SEEDS problem let us rapidly recall elementary
constructions of RK, ERK, weak and strong SRK methods. We will not deal with
time-adaptive methods here.

###  D.1 Derivative-free exponential ODE schemes

####  D.1.1 Runge-Kutta approach

Derivative-free schemes are obtained by comparing the Itô-Taylor expansion of
the above paragraph with expressions of $\mathbf{x}(t)$ in terms of its
intermediate evaluations between $s$ and $t=s+h$ and the Taylor expansions of
such evaluations. As a simple example, in the ODE regime

$\mathrm{d}\mathbf{x}(t)=f(\mathbf{x}(t),t)\mathrm{d}t.$ (51)

As such, analytic solutions to the above equations are of the form

$\displaystyle\mathbf{x}(t)$ $\displaystyle=$
$\displaystyle\mathbf{x}(s)+\int_{s}^{t}f(\mathbf{x}(\tau),\tau)\mathrm{d}\tau.$

Now, Taylor expansion gives, up to order 2:

$\displaystyle\mathbf{x}(t)$ $\displaystyle=$
$\displaystyle\mathbf{x}(s)+h\mathbf{x}^{\prime}(s)+\frac{h^{2}}{2}\mathbf{x}^{\prime\prime}(s)+\mathcal{O}(h^{3})$
$\displaystyle=$
$\displaystyle\mathbf{x}(s)+hf(\mathbf{x}(s),s)+\frac{h^{2}}{2}(\mathbf{x}^{\prime}(s)\partial_{\mathbf{x}}f(\mathbf{x}(s),s)+\partial_{t}f(\mathbf{x}(s),s))+\mathcal{O}(h^{3})$
$\displaystyle=$
$\displaystyle\mathbf{x}(s)+hf(\mathbf{x}(s),s)+\frac{h^{2}}{2}(f(\mathbf{x}(s),s)\partial_{\mathbf{x}}f(\mathbf{x}(s),s)+\partial_{t}f(\mathbf{x}(s),s))+\mathcal{O}(h^{3}).$

A straightforward recursion yields a Taylor expansion

$\displaystyle\mathbf{x}(t)$ $\displaystyle=$
$\displaystyle\mathbf{x}(s)+\sum_{k=1}^{n}\frac{h^{k}}{k!}L_{t}^{k-1}f(\mathbf{x}(t),t)+\int_{s}^{t}\cdots\int_{s}^{\tau}L_{t}^{n}f(\mathbf{x}(\tau),\tau)\mathrm{d}\tau^{n+1},$

where we denote the generalized infinitesimal operator of the solution
$\mathbf{x}$ of (51) by

$L_{t}(\bullet)=\partial_{t}(\bullet)+f(\mathbf{x}(t),t)\cdot\partial_{\mathbf{x}}(\bullet).$

Derivative-free methods seek to get rid of the derivatives in $L_{t}$ by
computing Taylor expansions of $f$ at well-chosen points. In the explicit one-
step case, this amounts of defining

$\widehat{\mathbf{x}}(t)=\widehat{\mathbf{x}}(s)+h\sum_{k=1}^{n}\Phi(\widehat{\mathbf{x}}(s),s,h),$

where the $\Phi$ function does not contain derivatives of $f$. The general
high order case is given by well-tuned coefficients in the following scheme

$\displaystyle\widehat{\mathbf{x}}(t)$ $\displaystyle:=$
$\displaystyle\widehat{\mathbf{x}}(s)+h\sum_{i=1}^{n}\alpha_{i}f(\mathbf{x}_{i},s+c_{i}h)$
$\displaystyle\mathbf{x}_{i}$ $\displaystyle=$
$\displaystyle\widehat{\mathbf{x}}(s)+h\sum_{j=1}^{n}a_{i,j}f(\mathbf{x}_{j},s+c_{j}h).$

By denoting $\alpha=[\alpha_{1}\ldots\alpha_{n}]$, $c=[c_{1}\ldots c_{n}]^{T}$
and $A=(a_{i,j})$, these can be represented by a Butcher tableau of the form

$\begin{array}[]{c|c}c&A\\\ \hline\cr&\alpha\end{array}$

In the following sections we will present extended version of this tableau for
representing more involved numerical schemes. As such, the Euler, midpoint,
Heun and general second order explicit methods are respectively:

$\begin{array}[]{c|c}0&0\\\
\hline\cr&1\end{array}\qquad\begin{array}[]{c|cc}0&0&0\\\ 1/2&1/2&0\\\
\hline\cr&0&1\end{array}\qquad\begin{array}[]{c|cc}0&0&0\\\ 1&1&0\\\
\hline\cr&1/2&1/2\end{array}\qquad\begin{array}[]{c|cc}0&0&0\\\
c_{2}&c_{2}&0\\\ \hline\cr&1-\frac{1}{2c_{2}}&\frac{1}{2c_{2}}\end{array}$

####  D.1.2 Exponential Runge-Kutta approach

We now concentrate on a non-autonomous semi-linear ODE of the form

$\mathrm{d}\mathbf{x}(t)=[f(t)\mathbf{x}(t)+g(\mathbf{x}(t),t)]\mathrm{d}t,$
(52)

where $f(t)=\bar{f}(t)\cdot\operatorname{Id}_{d}$ and
$g(\mathbf{x}(t),t)=\bar{g}(\mathbf{x}(t),t)\cdot\operatorname{Id}_{d}$ are
diagonal time-dependent matrices with identical diagonal coefficients. In
particular, we have $[f^{(k)}(t),f^{(l)}(s)]=0$, a property we will need to
facilitate exponential matrix multiplications. As $f(s)$ is a $d$-dimensional
diagonal matrix, the fundamental matrix for (52), usually given by the Peano-
Baker series, simplifies in this case to the form

$\Phi(t;s):=e^{\int_{s}^{t}f(r)\mathrm{d}r}=\sum_{n=0}^{\infty}\frac{1}{n!}\left(\int_{s}^{t}f(r)\mathrm{d}r\right)^{n},$

which satisfies $\Phi^{\prime}(t;s)=f(t)\Phi(t;s)$ and $\Phi(s;s)=1$. Then the
exact solution for (52) is given, via the variation of constants formula, by

$\displaystyle\mathbf{x}(t)$ $\displaystyle=$
$\displaystyle\Phi(t;s)\left(\mathbf{x}(s)+\int_{s}^{t}\Phi^{-1}(\tau;s)g(\mathbf{x}(\tau),\tau)\mathrm{d}\tau\right)=\Phi(t;s)\mathbf{x}(s)+\int_{s}^{t}\Phi(t;\tau)g(\mathbf{x}(\tau),\tau)\mathrm{d}\tau.$

In light of this integral form, one can formalize a general class of

exponential $n$-stage RK methods

$\displaystyle\mathbf{x}_{i}$ $\displaystyle=$
$\displaystyle\gamma_{i}(h,f_{s})\mathbf{x}(s)+\sum_{j=1}^{n}a_{ij}(h,f_{s})g(\mathbf{x}_{j},s+c_{j}h)$
$\displaystyle\mathbf{x}(t)$ $\displaystyle=$
$\displaystyle\gamma_{0}(h,f_{s})\mathbf{x}(s)+\sum_{i=1}^{n}b_{i}(h,f_{s})g(\mathbf{x}_{i},s+c_{i}h),$

where $\gamma_{0},\gamma_{i},a_{ij},b_{j}$ are functions of the step-size $h$,
$f$and $f_{s}(r):=\int_{0}^{r}f(s+\tau)\mathrm{d}\tau$. There are two possible
approaches to create exponential integrators, namely the exponential time-
differencing (ETD) approach that uses the variation of constants formula and
makes use of the $\varphi$ functions, and the Lawson approach, also know as
integrating factor (IF), which makes a change of variables on the above SDE
thus avoiding the use of the $\varphi$ functions but computing exponential
factors step-wise.

The Lawson and the ETD approaches for Exponential Euler schemes

There are two choices one can make when computing first order approximations
of

$\displaystyle\mathbf{x}(t)$ $\displaystyle=$
$\displaystyle\Phi(t;s)\mathbf{x}(s)+\int_{s}^{t}\Phi(t;\tau)g(\mathbf{x}(\tau),\tau)\mathrm{d}\tau.$

First, by interpolating $g(\mathbf{x}(\tau),\tau)$ as $g(\mathbf{x}(s),s)$ we
obtain

$\displaystyle\mathbf{x}(t)$ $\displaystyle=$
$\displaystyle\Phi(t;s)\mathbf{x}(s)+g(\mathbf{x}(s),s)\int_{s}^{t}\Phi(t;\tau)\mathrm{d}\tau.$

Now two choices remain. Either $\int_{s}^{t}\Phi(t;\tau)\mathrm{d}\tau$ is
computed exactly or we again interpolate $\Phi(t;\tau)$ as $\Phi(t;s)$. Taking
for simplicity $f\equiv A$ to be constant and by denoting $h=t-s$, the first
case yields the ETD Euler method

$\displaystyle\widehat{\mathbf{x}}(t)$ $\displaystyle=$
$\displaystyle\Phi(t;s)\widehat{\mathbf{x}}(s)+h\gamma_{1}(f,t,s)g(\widehat{\mathbf{x}}(s),s)$
$\displaystyle=$ $\displaystyle
e^{Ah}\widehat{\mathbf{x}}(s)+h\varphi_{1}(Ah)g(\widehat{\mathbf{x}}(s),s)$
$\displaystyle=$
$\displaystyle\widehat{\mathbf{x}}(s)+h\varphi_{1}(Ah)[g(\widehat{\mathbf{x}}(s),s)-\widehat{\mathbf{x}}(s)],$

and the second yields the IF Euler (also called Lawson-Euler) method

$\displaystyle\widehat{\mathbf{x}}(t)$ $\displaystyle=$
$\displaystyle\Phi(t;s)\widehat{\mathbf{x}}(s)+h\Phi(t;s)g(\widehat{\mathbf{x}}(s),s)$
$\displaystyle=$ $\displaystyle
e^{Ah}(\widehat{\mathbf{x}}(s)+hg(\widehat{\mathbf{x}}(s),s)).$

Now, consider a solution

$\displaystyle\mathbf{x}(t)$ $\displaystyle=$
$\displaystyle\Phi(t;s)\mathbf{x}(s)+\int_{0}^{h}\Phi(h;\tau)g(\mathbf{x}(s+\tau),s+\tau)\mathrm{d}\tau.$

Exponential methods then aim to approximate the term
$g(\mathbf{x}(s+\tau),s+\tau)$ by its interpolation polynomial in certain non-
confluent quadrature nodes $c_{1},\ldots,c_{n}$.

The ETD approach

In this case, the variation of constants formula yields

$\displaystyle\mathbf{x}(s+ch)$ $\displaystyle=$ $\displaystyle
e^{\int_{s}^{s+ch}a(\tau)\mathrm{d}\tau}\mathbf{x}(s)+\int_{s}^{s+ch}e^{\int_{\tau}^{s+ch}a(r)\mathrm{d}r}f(\mathbf{x}(\tau),\tau)\mathrm{d}\tau.$
$\displaystyle=$ $\displaystyle
e^{\int_{0}^{ch}a(s+\tau)\mathrm{d}\tau}\mathbf{x}(t)+\int_{0}^{ch}e^{\int_{\tau-s}^{ch}a(s+r)\mathrm{d}r}f(\mathbf{x}(s+\tau),s+\tau)\mathrm{d}\tau.$

Now, as before, the Taylor expansion of $f$ yields

$\displaystyle f(\mathbf{x}(s+\tau),s+\tau)$ $\displaystyle=$
$\displaystyle\sum_{j=1}^{q}\frac{\tau^{j-1}}{(j-1)!}f^{(j-1)}(\mathbf{x}(s),s)$
$\displaystyle+\int_{0}^{\tau}\frac{(\tau-\tau_{1})^{q-1}}{(q-1)!}f^{(q)}(\mathbf{x}(s+\tau_{1}),s+\tau_{1})\mathrm{d}\tau_{1}.$

Recall that the $\varphi$ functions are given in an integral form as follows

$\varphi_{k+1}(t)=\int_{0}^{1}e^{(1-\delta)t}\frac{\delta^{k}}{k!}\mathrm{d}\delta,$

which satisfy $\varphi_{k}(0)=\frac{1}{k!}$. Now denote

$\displaystyle\varphi_{j}(t,a)$ $\displaystyle:=$
$\displaystyle\frac{1}{t^{j}}\int_{0}^{t}e^{\int_{\tau}^{t}a(r)\mathrm{d}r}\frac{\tau^{j-1}}{(j-1)!}\mathrm{d}\tau,\qquad
j\geqslant 1.$ $\displaystyle\varphi_{1}(h,a)$ $\displaystyle=$
$\displaystyle\frac{1}{h}\int_{0}^{h}e^{\int_{\tau}^{h}a(r)\mathrm{d}r}\mathrm{d}\tau=\int_{0}^{1}e^{h\int_{\theta}^{1}a(r)\mathrm{d}r}\mathrm{d}\theta.$

The exact solution at $s+ch$ now reads

$\displaystyle\mathbf{x}(s+ch)$ $\displaystyle=$ $\displaystyle
e^{\int_{0}^{ch}a(s+\tau)\mathrm{d}\tau}\mathbf{x}(s)+\sum_{j=1}^{q}(ch)^{j}\varphi_{j}(ch,a)f^{(j-1)}(\mathbf{x}(s),s)$
$\displaystyle+\int_{0}^{ch}e^{\int_{\tau}^{ch}a(s+r)\mathrm{d}r}\left(\int_{0}^{\tau}\frac{(\tau-\tau_{1})^{q-1}}{(q-1)!}f^{(q)}(\mathbf{x}(s+\tau_{1}),s+\tau_{1})\mathrm{d}s\right)\mathrm{d}\tau_{1}.$

Using the left endpoint rule yields

$\displaystyle\mathbf{x}(t)$ $\displaystyle=$ $\displaystyle
e^{\int_{0}^{h}a(s+\tau)\mathrm{d}\tau}\mathbf{x}(s)+f(\mathbf{x}(s),s)\int_{0}^{h}e^{\int_{\tau}^{h}a(s+r)\mathrm{d}r}\mathrm{d}\tau+\mathcal{O}(h^{2})$
$\displaystyle=$ $\displaystyle
e^{\int_{0}^{h}a(s+\tau)\mathrm{d}\tau}\mathbf{x}(s)+h\varphi_{1}(h,a)f(\mathbf{x}(s),s)+\mathcal{O}(h^{2}).$

Second order examples:

The condition $b_{2}(z)c_{2}=\varphi_{2}(z)$ implies
$b_{2}=\varphi_{2}(z)/c_{2}$ and we obtain

$\displaystyle\widehat{\mathbf{x}}$ $\displaystyle=$ $\displaystyle
e^{chA}\mathbf{x}(s)+ch\varphi_{1}(chA)f(\mathbf{x}(s),s)$
$\displaystyle\mathbf{x}(t)$ $\displaystyle=$ $\displaystyle
e^{hA}\mathbf{x}(s)+h\left(\varphi_{1}(hA)-\frac{1}{c}\varphi_{2}(hA)\right)f(\mathbf{x}(s),s)+\frac{h}{c}\varphi_{2}(hA)f(\widehat{\mathbf{x}},s+ch).$
(53)

A second method is obtained by weakening the above condition to
$b_{2}(0)c_{2}=\varphi_{2}(0)=1/2$ giving

$\displaystyle\widehat{\mathbf{x}}$ $\displaystyle=$ $\displaystyle
e^{chA}\mathbf{x}(s)+ch\varphi_{1}(chA)f(\mathbf{x}(s),s)$
$\displaystyle\mathbf{x}(t)$ $\displaystyle=$ $\displaystyle
e^{hA}\mathbf{x}(s)+h\varphi_{1}(hA)\left(1-\frac{1}{2c}\right)f(\mathbf{x}(s),s)+\frac{h}{2c}\varphi_{1}(hA)f(\widehat{\mathbf{x}},s+ch).$
(54)

In some cases this method can suffer from order reduction and not reach order
2 of convergence. Moreover, setting $c=1/2$ gives the exponential midpoint
method:

$\displaystyle\widehat{\mathbf{x}}$ $\displaystyle=$ $\displaystyle
e^{chA}\mathbf{x}(s)+ch\varphi_{1}(chA)f(\mathbf{x}(s),s)$
$\displaystyle\mathbf{x}(t)$ $\displaystyle=$ $\displaystyle
e^{hA}\mathbf{x}(s)+h\varphi_{1}(hA)f(\widehat{\mathbf{x}},s+ch).$ (55)

The above order 1 and 2 exponential methods are represented by the following
exponential Butcher tableaux

$\begin{array}[]{c|c}0&0\\\
\hline\cr&\varphi_{1}\end{array}\qquad\begin{array}[]{l|ll}0&&\\\
c_{2}&c_{2}\varphi_{1,2}&\\\
\hline\cr&\varphi_{1}-\frac{1}{c_{2}}\varphi_{2}&\frac{1}{c_{2}}\varphi_{2}\end{array}\qquad\begin{array}[]{l|ll}0&&\\\
c_{2}&c_{2}\varphi_{1,2}&\\\
\hline\cr&\left(1-\frac{1}{2c_{2}}\right)\varphi_{1}&\frac{1}{2c_{2}}\varphi_{1}\end{array}\qquad\begin{array}[]{l|ll}0&&\\\
c_{2}&c_{2}\varphi_{1,2}&\\\ \hline\cr&0&\varphi_{1}\end{array}$

We check that when formally setting $A=0$ then the first two methods are
identical and give the generic 2nd order RK method, the choice $c=1$ gives the
Heun method and the choice $c=1/2$ gives the midpoint method.

Fourth order methods

It can be shown that ERK methods need at least 5 stages to achieve order 4. By
setting formally $A=0$ these methods do not have a non-exponential counterpart
to our knowledge.

5-stage sequential We have a fourth-order ERK scheme given by

$\begin{array}[]{c|ccccc}0&&&&&\\\ \frac{1}{2}&\frac{1}{2}\varphi_{1,2}&&&&\\\
\frac{1}{2}&\frac{1}{2}\varphi_{1,3}-\varphi_{2,3}&\varphi_{2,3}&&&\\\
1&\varphi_{1,4}-2\varphi_{2,4}&\varphi_{2,4}&\varphi_{2,4}&&\\\
\frac{1}{2}&\frac{1}{2}\varphi_{1,5}-2a_{5,2}-a_{5,4}&a_{5,2}&a_{5,2}&\frac{1}{4}\varphi_{2,5}-a_{5,2}&\\\
\hline\cr&\varphi_{1}-3\varphi_{2}+4\varphi_{3}&0&0&-\varphi_{2}+4\varphi_{3}&4\varphi_{2}-8\varphi_{3}\end{array}$

with

$a_{5,2}=c_{5}\varphi_{2,5}-\varphi_{3,4}+c_{5}^{2}\varphi_{2,4}-c_{5}\varphi_{3,5}.$

This can also be represented by Rosenbrock-like Butcher tableau

$\begin{array}[]{c|ccccc}0&&&&&\\\ \frac{1}{2}&\frac{1}{2}\varphi_{1,2}&&&&\\\
\frac{1}{2}&\frac{1}{2}\varphi_{1,3}&\varphi_{2,3}&&&\\\
1&\varphi_{1}&\varphi_{2}&\varphi_{2}&&\\\
\frac{1}{2}&\frac{1}{2}\varphi_{1,5}&a_{5,2}&a_{5,2}&\frac{1}{4}\varphi_{2,5}-a_{5,2}&\\\
\hline\cr&\varphi_{1}&0&0&4\varphi_{3}-\varphi_{2}&4\varphi_{2}-8\varphi_{3}\end{array}$

with

$a_{5,2}=c_{5}\varphi_{2,5}-c_{3}\varphi_{3,3}+c_{5}^{2}\varphi_{2}-\varphi_{3}.$

This traduces into (recall that
$D_{i}=f(\mathbf{x}_{i},s+c_{i}h)-f(\mathbf{x}(s),s)$)

$\displaystyle\mathbf{x}(t)$ $\displaystyle=$
$\displaystyle\mathbf{x}(s)+h(\varphi_{1}(hA))F(\mathbf{x}(s),s)+h(4\varphi_{3}(hA)-\varphi_{2}(hA))D_{4}$
$\displaystyle+h(4\varphi_{2}(hA)-8\varphi_{3}(hA))D_{5}$
$\displaystyle\mathbf{x}_{2}$ $\displaystyle=$
$\displaystyle\mathbf{x}(s)+c_{2}h\varphi_{1}(c_{2}hA)F(\mathbf{x}(s),s)$
$\displaystyle\mathbf{x}_{3}$ $\displaystyle=$
$\displaystyle\mathbf{x}(s)+hc_{3}\varphi_{1}(c_{3}hA)F(\mathbf{x}(s),s)+h\varphi_{2}(c_{3}hA)D_{2}$
$\displaystyle\mathbf{x}_{4}$ $\displaystyle=$
$\displaystyle\mathbf{x}(s)+h\varphi_{1}(hA)F(\mathbf{x}(s),s)+h\varphi_{2}(hA)(D_{2}+D_{3})$
$\displaystyle\mathbf{x}_{5}$ $\displaystyle=$
$\displaystyle\mathbf{x}(s)+h\varphi_{1}(hA)F(\mathbf{x}(s),s)+ha_{5,2}(hA)(D_{2}+D_{3})+h(c_{5}^{2}\varphi_{2}(c_{5}hA)-a_{5,2}(hA))D_{4}.$

Inspired by this, we deduce the following fourth order Algorithm 5 specialized
for DPMs in the VP noise prediction regime.

Algorithm 5 DPM-Solver-4

1:def DPM-
Solver-4$(\bm{\epsilon}_{\theta},\tilde{\mathbf{x}}_{t_{i-1}},t_{i-1},t_{i},r=0.5)$:

2: $h_{i}\leftarrow\lambda_{t_{i}}-\lambda_{t_{i-1}}$

3: $s_{2}\leftarrow t_{\lambda}\left(\lambda_{t_{i-1}}+rh_{i}\right),\quad
s_{3}\leftarrow t_{\lambda}\left(\lambda_{t_{i-1}}+rh_{i}\right)$

4: $s_{4}\leftarrow t_{\lambda}\left(\lambda_{t_{i-1}}+h_{i}\right),\quad
s_{5}\leftarrow t_{\lambda}\left(\lambda_{t_{i-1}}+rh_{i}\right)$

5:
$\bm{k}_{1}\leftarrow\bm{\epsilon}_{\theta}(\tilde{\mathbf{x}}_{t_{i-1}},t_{i-1})$

6:
$\bm{k}_{2}\leftarrow\frac{\alpha_{s_{2}}}{\alpha_{s}}\tilde{\mathbf{x}}_{t_{i-1}}-\sigma_{s_{2}}(e^{rh}-1)\bm{k_{1}}$

7:
$\bm{k}_{3}\leftarrow\frac{\alpha_{s_{3}}}{\alpha_{s}}\tilde{\mathbf{x}}_{t_{i-1}}-\sigma_{s_{3}}(e^{rh}-1)\bm{k_{1}}-\sigma_{s_{3}}\left(4\frac{e^{rh}-1}{h}-2\right)[\epsilon_{\theta}(\bm{k_{2}},s_{2})-\bm{k_{1}}]$

8:
$\bm{k}_{4}\leftarrow\frac{\alpha_{s_{4}}}{\alpha_{s}}\tilde{\mathbf{x}}_{t_{i-1}}-\sigma_{s_{4}}(e^{h}-1)\bm{k_{1}}-\sigma_{s_{4}}\left(\frac{e^{h}-1}{h}-1\right)[\epsilon_{\theta}(\bm{k_{3}},s_{3})+\epsilon_{\theta}(\bm{k_{2}},s_{2})-2\bm{k_{1}}]$

9:
$A=\sigma_{s_{5}}(e^{rh}-1)\bm{k_{1}}-\frac{1}{4}\sigma_{s_{5}}\left(\frac{e^{h}-1}{h}-1\right)[\bm{k_{1}}+\epsilon_{\theta}(\bm{k_{2}},s_{2})+\epsilon_{\theta}(\bm{k_{3}},s_{3})]$

10:
$B=\sigma_{s_{5}}\left(\frac{e^{rh}-1}{h}-\frac{1}{2}\right)[\bm{k_{1}}+4\epsilon_{\theta}(\bm{k_{2}},s_{2})+4\epsilon_{\theta}(\bm{k_{3}},s_{3})-\epsilon_{\theta}(\bm{k_{4}},s_{4})]$

11:
$C=\sigma_{s_{5}}\left(\frac{e^{h}-1+4(e^{rh}-1)-3h}{h^{2}}-1\right)[-\bm{k_{1}}-\epsilon_{\theta}(\bm{k_{2}},s_{2})-\epsilon_{\theta}(\bm{k_{3}},s_{3})+\epsilon_{\theta}(\bm{k_{4}},s_{4})]$

12:
$\bm{k}_{5}\leftarrow\frac{\alpha_{s_{5}}}{\alpha_{t_{i-1}}}\tilde{\mathbf{x}}_{t_{i-1}}-A-B-C$

13:
$D=\sigma_{t}(e^{h}-1)\bm{k_{1}}-\sigma_{t}\left(\frac{e^{h}-1}{h}-1\right)[4\epsilon_{\theta}(\bm{k_{5}},s_{5})-\epsilon_{\theta}(\bm{k_{4}},s_{4})-3\bm{k_{1}}]$

14:
$E=\sigma_{t}\left(4\frac{e^{h}-1-h}{h^{2}}-2\right)[\bm{k_{1}}+\epsilon_{\theta}(\bm{k_{4}},s_{4})-2\epsilon_{\theta}(\bm{k_{5}},s_{5})]$

15:
$\tilde{\mathbf{x}}_{t_{i}}\leftarrow\frac{\alpha_{t}}{\alpha_{t_{i-1}}}\tilde{\mathbf{x}}_{t_{i-1}}-D-E$

16: Return $\mathbf{x}_{t_{i}}$

###  D.2 Derivative-free exponential SDE schemes

Let $\mathbf{x}_{t}$ be the path at the continuous limit $h\to 0$, and
$\\{\widehat{\mathbf{x}}_{t}\\}_{t_{0}}^{t_{M}}$ be the discretized numerical
path, computed by a numerical scheme $\mathcal{S}$ with $M=1/h$ steps of
length $h>0$. Then, $\mathcal{S}$ has

  1. 1.

strong order of convergence $\gamma$ if there is $K>0$ such that

$\mathbb{E}[|\mathbf{x}_{t_{M}}-\widehat{\mathbf{x}}_{t_{M}}|]\leq
Kh^{\gamma},$ (56)

  2. 2.

weak order of convergence $\beta$ if there is $K>0$ and a function class
$\mathcal{K}$ such that

$|\mathbb{E}[\phi(\mathbf{x}_{t_{M}})]-\mathbb{E}[\phi(\widehat{\mathbf{x}}_{t_{M}})]|\leq
Kh^{\beta},\qquad\forall\phi(\cdot)\in\mathcal{K}.$ (57)

Strong convergence is concerned with the precision of the path, while the weak
convergence is with the precision of the moments. As, for diffusion models,
the center of attention is the evolution of the probability densities rather
than that of the noising process of single data samples, weak convergence is
enough to guarantee the well-conditioning of our numerical schemes. Moreover,
when the diffusion coefficient vanishes, then both strong and weak convergence
(with the choice $\phi=\operatorname{id}$) reduce to the usual deterministic
convergence criterion for ODEs.

####  D.2.1 Strong and Weak Stochastic Runge-Kutta approach

In all what follows $\omega$ is considered a $d$-dimensional Wiener process
(with identity diffusion matrix).

Consider the following SDE

$\mathrm{d}\mathbf{x}(t)=f(\mathbf{x}(t),t)\mathrm{d}t+g(t)\mathrm{d}\omega(t),$

where $g(t)=\hat{g}(t)\cdot\operatorname{Id}_{d}$ is considered here as a
diagonal matrix with identical diagonal entries $\hat{g}(t)$. Given an initial
value independent of $\omega$, the integral form of $\mathbf{x}(t)$ is given
by

$\displaystyle\mathbf{x}(t)$ $\displaystyle=$
$\displaystyle\mathbf{x}(s)+\int_{s}^{t}f(\mathbf{x}(\tau),\tau)d\tau+\int_{s}^{t}g(\tau)\mathrm{d}\omega(\tau).$

The idea underlying stochastic numerical schemes is very similar to the one
used in the deterministic approach, that is to take expansions of the terms
inside the integrals based at the integral’s initial value and replace the
obtained derivatives that appear by interpolated approximations. A key
difference here is that as we have to consider Itô-Taylor expansions, the
infinitesimal operators are different but most importantly most of the
stochastic iterated integrals will need to be approximated in an appropriate
sense, whenever that is possible. We will develop this expansion up to triple
integrals.

####  D.2.2 Truncated Itô-Taylor expansions

Applying Itô formula to $h=f$ or $g$ yields

$\displaystyle h(\mathbf{x}(t),t)$ $\displaystyle=$ $\displaystyle
h(\mathbf{x}(s),s)+\int_{s}^{t}g(\tau)\cdot\partial_{\mathbf{x}}h(\mathbf{x}(\tau),\tau)\mathrm{d}\omega(\tau)$
$\displaystyle+\int_{s}^{t}\left(\partial_{t}h(\mathbf{x}(\tau),\tau)+f(\mathbf{x}(\tau),\tau)\cdot\partial_{\mathbf{x}}h(\mathbf{x}(\tau),\tau)+\frac{g^{2}(\tau)}{2}\partial_{\mathbf{x}^{2}}^{2}h(\mathbf{x}(\tau),\tau)\right)\mathrm{d}\tau.$

This allows us to define two differential operators $L,G$ as

$\displaystyle L_{t}$ $\displaystyle=$
$\displaystyle\partial_{t}+f(\mathbf{x}(t),t)\cdot\partial_{\mathbf{x}}+\frac{g^{2}(t)}{2}\cdot\partial^{2}_{\mathbf{x}}$
(58) $\displaystyle G_{t}$ $\displaystyle=$ $\displaystyle
g(t)\cdot\partial_{\mathbf{x}}.$ (59)

In particular $L_{t}g(t)=\partial_{t}g(t)$ and $G_{t}g(t)=0$. With this
notation, we have

$\displaystyle h(\mathbf{x}(t),t)$ $\displaystyle=$ $\displaystyle
h(\mathbf{x}(s),s)+\int_{s}^{t}L_{t}h(\mathbf{x}(\tau),\tau)\mathrm{d}\tau+\int_{s}^{t}G_{t}h(\mathbf{x}(\tau),\tau)\mathrm{d}\omega(\tau),$

so our solution reads

$\displaystyle\mathbf{x}(t)$ $\displaystyle=$
$\displaystyle\mathbf{x}(s)+\int_{s}^{t}f(\mathbf{x}(\tau_{1}),\tau_{1})\mathrm{d}\tau_{1}+\int_{s}^{t}g(\tau_{1})\mathrm{d}\omega(\tau_{1})$
$\displaystyle=$
$\displaystyle\mathbf{x}(s)+\int_{s}^{t}\left(f(\mathbf{x}(t),t)+\int_{s}^{\tau_{1}}L_{t}f(\mathbf{x}(\tau_{2}),\tau_{2})\mathrm{d}\tau_{2}+\int_{s}^{\tau_{1}}G_{t}f(\mathbf{x}(\tau_{2}),\tau_{2})d\omega(\tau_{2})\right)\mathrm{d}\tau_{1}$
$\displaystyle+\int_{s}^{t}\left(g(t)+\int_{s}^{\tau_{1}}L_{t}g(\tau_{2})\mathrm{d}\tau_{2}+\int_{s}^{\tau}G_{t}g(\tau_{2})\mathrm{d}\omega(\tau_{2})\right)\mathrm{d}\omega(\tau_{1})$
$\displaystyle=$
$\displaystyle\mathbf{x}(s)+f(\mathbf{x}(s),s)h+g(t)(\omega(t)-\omega(s))+R_{1}.$

Now $G_{t}g(\tau_{2})=0$ and

$g(s)+\int_{s}^{\tau_{1}}L_{t}g(\tau_{2})\mathrm{d}\tau_{2}=g(s)+\int_{s}^{\tau_{1}}\partial_{t}g(\tau_{2})\mathrm{d}\tau_{2}=g(s)+g(\tau_{1})-g(s)=g(\tau_{1}).$

So we have

$\int_{s}^{t}\left(g(s)+\int_{s}^{\tau_{1}}L_{t}g(\tau_{2})\mathrm{d}\tau_{2}+\int_{s}^{\tau}G_{t}g(\tau_{2})\mathrm{d}\omega(\tau_{2})\right)\mathrm{d}\omega(\tau_{1})=g(s)(\omega(t)-\omega(s)),$

and now our solution reads

$\mathbf{x}(t)=\mathbf{x}(s)+f(\mathbf{x}(s),s)h+g(s)(\omega(t)-\omega(s))+R_{1},$

where

$\displaystyle R_{1}$ $\displaystyle=$
$\displaystyle\int_{s}^{t}\int_{s}^{\tau_{1}}L_{t}f(\mathbf{x}(\tau_{2}),\tau_{2})\mathrm{d}\tau_{2}\mathrm{d}\tau_{1}+\int_{s}^{t}\int_{s}^{\tau_{1}}G_{t}f(\mathbf{x}(\tau_{2}),\tau_{2})\mathrm{d}\omega(\tau_{2})\mathrm{d}\tau_{1}.$

Now we have

$\displaystyle L_{t}f(\mathbf{x}(t),t)$ $\displaystyle=$ $\displaystyle
L_{t}f(\mathbf{x}(s),s)+\int_{s}^{t}g(\mathbf{x}(\tau),\tau)\cdot\partial_{\mathbf{x}}L_{t}f(\mathbf{x}(\tau),\tau)\mathrm{d}\omega(\tau)$
$\displaystyle+\int_{s}^{t}\left(\partial_{t}L_{t}f(\mathbf{x}(\tau),\tau)+f(\mathbf{x}(\tau),\tau)\cdot\partial_{\mathbf{x}}L_{t}f(\mathbf{x}(\tau),\tau)+\frac{g^{2}(\mathbf{x}(\tau),\tau)}{2}\partial_{\mathbf{x}^{2}}^{2}L_{t}f(\mathbf{x}(\tau),\tau)\right)\mathrm{d}\tau$
$\displaystyle=$ $\displaystyle
L_{t}f(\mathbf{x}(s),s)+\int_{s}^{t}L^{2}_{t}f(\mathbf{x}(\tau),\tau)\mathrm{d}\tau+\int_{s}^{t}G_{t}L_{t}f(\mathbf{x}(\tau),\tau)\mathrm{d}\omega(\tau),$

and

$\displaystyle G_{t}f(\mathbf{x}(t),t)$ $\displaystyle=$ $\displaystyle
G_{t}f(\mathbf{x}(s),s)+\int_{s}^{t}g(\mathbf{x}(\tau),\tau)\cdot\partial_{\mathbf{x}}G_{t}f(\mathbf{x}(\tau),\tau)\mathrm{d}\omega(\tau)$
$\displaystyle+\int_{s}^{t}\left(\partial_{t}G_{t}f(\mathbf{x}(\tau),\tau)+f(\mathbf{x}(\tau),\tau)\cdot\partial_{\mathbf{x}}G_{t}f(\mathbf{x}(\tau),\tau)+\frac{g^{2}(\mathbf{x}(\tau),\tau)}{2}\partial_{\mathbf{x}^{2}}^{2}G_{t}f(\mathbf{x}(\tau),\tau)\right)\mathrm{d}\tau$
$\displaystyle=$ $\displaystyle
G_{t}f(\mathbf{x}(s),s)+\int_{s}^{t}L_{t}G_{t}f(\mathbf{x}(\tau),\tau)\mathrm{d}\tau+\int_{s}^{t}G^{2}_{t}f(\mathbf{x}(\tau),\tau)\mathrm{d}\omega(\tau).$

Now, if we denote $\mathrm{d}\omega^{0}(\tau)=\mathrm{d}\tau$,
$\mathrm{d}\omega^{1}(\tau)=\mathrm{d}\omega(\tau)$ and

$I_{(i)}=\int_{s}^{t}\mathrm{d}\omega^{i}(\tau_{1})\qquad
I_{(i,j)}=\int_{s}^{t}\int_{s}^{\tau_{1}}\mathrm{d}\omega^{j}(\tau_{2})\mathrm{d}\omega^{i}(\tau_{1})\qquad
i,j\in\\{0,1\\},$

applying the same procedure to $R_{1}$ leads to

$\displaystyle R_{1}$ $\displaystyle=$
$\displaystyle\int_{s}^{t}\int_{s}^{\tau_{1}}L_{t}f(\mathbf{x}(\tau_{2}),\tau_{2})\mathrm{d}\tau_{2}\mathrm{d}\tau_{1}+\int_{s}^{t}\int_{s}^{\tau_{1}}G_{t}f(\mathbf{x}(\tau_{2}),\tau_{2})\mathrm{d}\omega(\tau_{2})\mathrm{d}\tau_{1}$
$\displaystyle=$
$\displaystyle\int_{s}^{t}\int_{s}^{\tau_{1}}\left[L_{t}f(\mathbf{x}(t),t)+\int_{s}^{\tau_{2}}L^{2}_{t}f(\mathbf{x}(\tau_{3}),\tau_{3})\mathrm{d}\tau_{3}+\int_{s}^{\tau_{2}}G_{t}L_{t}f(\mathbf{x}(\tau_{3}),\tau_{3})\mathrm{d}\omega(\tau_{3})\right]\mathrm{d}\tau_{2}\tau_{1}$
$\displaystyle+\int_{s}^{t}\int_{s}^{\tau_{1}}\left[G_{t}f(\mathbf{x}(t),t)+\int_{s}^{\tau_{2}}L_{t}G_{t}f(\mathbf{x}(\tau_{3}),\tau_{3})\mathrm{d}\tau+\int_{s}^{\tau_{2}}G^{2}_{t}f(\mathbf{x}(\tau_{3}),\tau_{3})\mathrm{d}\omega(\tau_{3})\right]\mathrm{d}\omega(\tau_{2})\mathrm{d}\tau_{1}$
$\displaystyle=$ $\displaystyle
L_{t}f(\mathbf{x}(t),t)\int_{s}^{t}\int_{s}^{\tau_{1}}\mathrm{d}\tau_{2}\mathrm{d}\tau_{1}+G_{t}f(\mathbf{x}(t),t)\int_{s}^{t}\int_{s}^{\tau_{1}}\mathrm{d}\omega(\tau_{2})\mathrm{d}\tau_{1}$
$\displaystyle+\int_{s}^{t}\int_{s}^{\tau_{1}}\int_{s}^{\tau_{2}}L^{2}_{t}f(\mathbf{x}(\tau_{3}),\tau_{3})\mathrm{d}\tau_{3}\mathrm{d}\tau_{2}\mathrm{d}\tau_{1}+\int_{s}^{t}\int_{s}^{\tau_{1}}\int_{s}^{\tau_{2}}G_{t}L_{t}f(\mathbf{x}(\tau_{3}),\tau_{3})\mathrm{d}\omega(\tau_{3})\mathrm{d}\tau_{2}\mathrm{d}\tau_{1}$
$\displaystyle+\int_{s}^{t}\int_{s}^{\tau_{1}}\int_{s}^{\tau_{2}}L_{t}G_{t}f(\mathbf{x}(\tau_{3}),\tau_{3})\mathrm{d}\tau_{3}\mathrm{d}\omega(\tau_{2})\mathrm{d}\tau_{1}$
$\displaystyle+\int_{s}^{t}\int_{s}^{\tau_{1}}\int_{s}^{\tau_{2}}G^{2}_{t}f(\mathbf{x}(\tau_{3}),\tau_{3})\mathrm{d}\omega(\tau_{3})\mathrm{d}\omega(\tau_{2})\mathrm{d}\tau_{1}$
$\displaystyle=$ $\displaystyle
L_{t}f(\mathbf{x}(s),s)I_{(0,0)}+G_{t}f(\mathbf{x}(s),s)I_{(0,1)}+L^{2}_{t}f(\mathbf{x}(s),s)I_{(0,0,0)}+G_{t}L_{t}f(\mathbf{x}(s),s)I_{(0,0,1)}$
$\displaystyle+L_{t}G_{t}f(\mathbf{x}(s),s)I_{(0,1,0)}+G^{2}_{t}f(\mathbf{x}(s),s)I_{(1,1,0)}+R_{2},$

and with $R_{2}$ consisting on quadruple integrals. As such, our solution now
reads

$\displaystyle\mathbf{x}(t)$ $\displaystyle=$
$\displaystyle\mathbf{x}(s)+f(\mathbf{x}(s),s)h+g(s)(\omega(t)-\omega(s))+L_{t}f(\mathbf{x}(s),s)I_{(0,0)}+G_{t}f(\mathbf{x}(s),s)I_{(0,1)}$
$\displaystyle+L^{2}_{t}f(\mathbf{x}(s),s)I_{(0,0,0)}+G_{t}L_{t}f(\mathbf{x}(s),s)I_{(0,0,1)}$
$\displaystyle+L_{t}G_{t}f(\mathbf{x}(s),s)I_{(0,1,0)}+G^{2}_{t}f(\mathbf{x}(s),s)I_{(1,1,0)}+R_{2}.$

Now, in the SDE regime, one cannot get rid of the iterated Itô integral and so
stochastic RK methods cannot be derived as simple extensions of their
deterministic counterparts. In order to continue we now take into account the
fact that the diffusion SDE has additive and diagonal noise. In this case,
both the Itô and the Stratonovich SDE coincide.

Iterated integrals

Now, for simplicity, set $t=0$. We then have $I_{(0)}=h$,
$I_{(0,0)}=\int_{0}^{h}\int_{0}^{\tau_{1}}\mathrm{d}\tau_{2}\mathrm{d}\tau_{1}=\frac{h^{2}}{2}$,
$I_{(0,0,0)}=\frac{h^{3}}{6}$. Now notice that

$\displaystyle I_{(1)}$ $\displaystyle:=$
$\displaystyle\hat{w}_{h}\sim\mathcal{N}(0,h)$ $\displaystyle I_{(0,1)}$
$\displaystyle:=$
$\displaystyle\hat{z}_{h}:=\int_{0}^{h}\int_{0}^{\tau_{1}}\mathrm{d}\omega(\tau_{2})\mathrm{d}\tau_{1}=\underset{n\rightarrow\infty}{\lim}\frac{h}{n}\sum^{n-1}_{i=0}\sum^{i}_{j=1}\epsilon_{j},\qquad\epsilon_{j}\sim\mathcal{N}\left(0,\frac{h}{n}\right).$

Additionally, $\hat{w}_{h}$ and $\hat{z}_{h}$ satisfy
$\mathbb{E}[\hat{w}_{h}^{2}]=h$,

$\displaystyle\mathbb{E}[(\hat{w}_{h}h-\hat{z}_{h})^{2}]$ $\displaystyle=$
$\displaystyle\mathbb{E}\left[\left(\int_{0}^{h}\tau\mathrm{d}\omega(\tau)\right)^{2}\right]=\frac{1}{3}h^{3}$
$\displaystyle\mathbb{E}[\hat{w}_{h}\hat{z}_{h}]$ $\displaystyle=$
$\displaystyle\mathbb{E}\left[\hat{w}_{h}\int_{0}^{h}\tau\mathrm{d}\omega(\tau)\right]=\mathbb{E}\left[\int_{0}^{h}\tau\mathrm{d}\tau\right]=\frac{1}{2}h^{2}$
$\displaystyle\mathbb{E}[\hat{z}_{h}^{2}]$ $\displaystyle=$
$\displaystyle\mathbb{E}[(\hat{w}_{h}h-\hat{z}_{h})^{2}-h^{2}\hat{w}_{h}^{2}+2h\hat{w}_{h}\hat{z}_{h}]=\frac{1}{3}h^{3}.$

####  D.2.3 Integral approximations

Weak Approximations

When crafting weak stochastic approximations to SDEs one may replace multiple
Itô integrals by other random variables satisfying the corresponding moment
conditions. We will denote $\hat{I}_{\alpha}$ the approximation of
$I_{\alpha}$ for $\alpha$ a multi-index following [13, Corollary 5.12.1]. Of
course, the deterministic integrals $I_{(0,\ldots,0)}$ need not to be
approximated.

First order approximations

The random variable $\hat{I}_{(1)}$ must satisfy for some constant $K$:

$|\mathbb{E}[\hat{I}_{(1)}]|+|\mathbb{E}[(\hat{I}_{(1)})^{3}]|+|\mathbb{E}[(\hat{I}_{(1)})^{2}-h]|\leqslant
Kh^{2}$

Two possible choices for $\hat{I}_{(1)}$ are either
$\hat{I}_{(1)}\sim\mathcal{N}(0,h)$ or $\hat{I}_{(1)}$ is a two-pointed
distributed discrete random variable with

$\mathbb{P}\left[\hat{I}_{(1)}=\pm\sqrt{h}\right]=\frac{1}{2}.$

Second order approximations

The random variable $\hat{I}_{(1)}$ must satisfy for some constant $K$:

$|\mathbb{E}[\hat{I}_{(1)}]|+|\mathbb{E}[(\hat{I}_{(1)})^{3}]|+|\mathbb{E}[(\hat{I}_{(1)})^{5}]|+|\mathbb{E}[(\hat{I}_{(1)})^{2}-h]|+|\mathbb{E}[(\hat{I}_{(1)})^{4}-3h^{2}]|\leqslant
Kh^{2}.$

Two possible choices for $\hat{I}_{(1)}$ are either
$\hat{I}_{(1)}\sim\mathcal{N}(0,h)$ or $\hat{I}_{(1)}$ is a three-pointed
distributed random variable with

$\mathbb{P}\left[\hat{I}_{(1)}=\pm\sqrt{3h}\right]=\frac{1}{6},\qquad\mathbb{P}[\hat{I}_{(1)}=0]=\frac{2}{3}.$

The rest follows from the above calculations:

$\hat{I}_{(0,1)}=\frac{1}{2}h\hat{I}_{(1)},\qquad i=0,1.$

Third order approximations

One can choose $\hat{I}_{(1)}\sim\mathcal{N}(0,h)$,
$\hat{I}_{(0,1)}\sim\mathcal{N}\left(0,\frac{1}{3}h^{3}\right)$ satisfying
$\mathbb{E}[\hat{I}_{(1)}\hat{I}_{(0,1)}]=\frac{1}{2}h^{2}$.

Then, one can deduce the following:

$\displaystyle\hat{I}_{(1,0,0)}=\hat{I}_{(0,1,0)}=\hat{I}_{(0,0,1)}$
$\displaystyle=$ $\displaystyle\frac{1}{6}h^{2}\hat{I}_{(1)}$
$\displaystyle\hat{I}_{(0,1,1)}=\hat{I}_{(1,0,1)}=\hat{I}_{(1,1,0)}$
$\displaystyle=$ $\displaystyle\frac{1}{6}h(\hat{I}_{(1)}^{2}-h).$

Thus, we can write the solution weak approximation as

$\displaystyle\mathbf{x}(t)$ $\displaystyle=$
$\displaystyle\mathbf{x}(s)+f(\mathbf{x}(s),s)h+g(s)\hat{I}_{(1)}+L_{t}f(\mathbf{x}(s),s)\frac{h^{2}}{2}+R_{2}$
(60)
$\displaystyle+G_{t}f(\mathbf{x}(s),s)\hat{I}_{(0,1)}+L^{2}_{t}f(\mathbf{x}(s),s)\frac{h^{3}}{6}+G_{t}L_{t}f(\mathbf{x}(s),s)\frac{1}{6}h^{2}\hat{I}_{(1)}$
$\displaystyle+L_{t}G_{t}f(\mathbf{x}(s),s)\frac{1}{6}h^{2}\hat{I}_{(1)}+G^{2}_{t}f(\mathbf{x}(s),s)\frac{1}{6}h(\hat{I}_{(1)}^{2}-h).$

An example of such a pair
$(\hat{I}_{(1)},\hat{I}_{(0,1)})=(\hat{w}_{h},\hat{z}_{h})$ can be easily
obtained as follows

$\displaystyle\left[\begin{array}[]{c}\hat{w}_{h}\\\
\hat{z}_{h}\end{array}\right]$ $\displaystyle=$
$\displaystyle\left[\begin{array}[]{cc}\sqrt{h}&0\\\
\frac{h\sqrt{h}}{2}&\frac{h\sqrt{h}}{2\sqrt{3}}\end{array}\right]\left[\begin{array}[]{c}u_{1}\\\
u_{2}\end{array}\right],\qquad
u_{1},u_{2}\overset{\text{i.i.d.}}{\sim}\mathcal{N}(0,1).$

Indeed, for such a pair we have

$\displaystyle\mathbb{E}\left[\left[\begin{array}[]{c}\hat{w}_{h}\\\
\hat{z}_{h}\end{array}\right]\left[\begin{array}[]{cc}\hat{w}_{h}&\hat{z}_{h}\end{array}\right]\right]$
$\displaystyle=$ $\displaystyle\left[\begin{array}[]{cc}h&\frac{h^{2}}{2}\\\
\frac{h^{2}}{2}&\frac{h^{3}}{3}\end{array}\right].$

In light of the above expression, the truncated Taylor expansions we refer in
the main part of the paper consists on the consideration of only the
coefficients in $L_{t}^{k}$. The only noise noise contribution we will
consider corresponds to $g(s)\hat{I}_{(1)}$.

##  Appendix E Experiment Details

We evaluate the Fréchet inception distance (FID) after generating 50K samples
with each solver, and compare with the statistics of real-data. In our
experiments we make use of the code from [12] for continuously trained models
as well as their reference FID stats111https://nvlabs-fi-
cdn.nvidia.com/edm/fid-refs/ and that of [15] for discretely trained models.

All the experiments of SEEDS for continuous-time models are parameterized
within the EDM framework with the discretization of type EDM, linear schedule,
and scaling none as described on [12] in noise prediction mode unless
explicitly stated. We use the SEEDS-3 method that has 3 NFEs per step and
fixed step-size and report FID scores at NFEs divisible by 3.

We leverage the explicit Langevin-like “churn” trick using in [12] to add or
remove noise in the sampling phase. Specifically, [12] uses 4 hyper-parameters
$S_{\text{churn}},S_{\text{min}},S_{\text{max}}$ and $S_{\text{noise}}$ in
which $S_{\text{churn}}$ controls the overall amount of stochasticity added
before giving the input to the SEEDS-3 method when the noise level (or time
step in EDM configuration) $t_{i}$ is contained in the noise interval
$[S_{\text{min}},S_{\text{max}}]$. It means that the EDM proposed sampler is
stochastic under some conditions of those hyper-parameters and deterministic
otherwise, while our method is completely stochastic. In our experiments, we
set $S_{\text{churn}}=0$ except for ImageNet-64 EDM optimized model. We
noticed that using the additional stochasticity indeed helps to improve the
image quality as in Fig. 1 (c). Moreover, setting $S_{\text{noise}}$ slightly
above 1 might correct the errors in earlier steps more effectively as
indicated in [12].

###  E.1 Pre-trained model specifications

For producing the CIFAR-10 time-continuous results in Table 1, we use the VP
DDPM++ continuous architecture. These models are publicly available in
conditional222https://nvlabs-fi-
cdn.nvidia.com/edm/pretrained/baseline/baseline-cifar10-32x32-cond-vp.pkl and
unconditional333https://nvlabs-fi-
cdn.nvidia.com/edm/pretrained/baseline/baseline-cifar10-32x32-uncond-vp.pkl
versions and were directly derived from [24] under the Apache 2.0 license. On
the unconditional mode (Figure 1 (a-b)), we first generate the FID curves of 3
types of DPM-Solver (with orders 1, 2 and 3) using the updates from their
official implementation444https://github.com/LuChengTHU/dpm-solver in noise
prediction mode. Taking profit of the tuning advancements proposed by [12], we
used a linear noise schedule with $\beta_{\text{d}}=19.1$ and
$\beta_{\min}=0.1$ that slightly differs from the original parameters from
[24] but were proven beneficial. We set the end time of sampling $\varepsilon$
to $1$e-$4$ as recommended by [15, Appendix D.2]. The values of all benchmark
models for Figure 1 (b) were taken directly from tables provided by [15].

In our FFHQ-64 experiments, we employ the unconditional VP
pretrained555https://nvlabs-fi-
cdn.nvidia.com/edm/pretrained/baseline/baseline-ffhq-64x64-uncond-vp.pkl model
provided by [12].

For the CelebA-64 experiments, we use the pre-trained VP unconditional model
whose checkpoint666https://drive.google.com/file/d/1R_H-
fJYXSH79wfSKs9D-fuKQVan5L-GR/view?usp=sharing is provided by [23]. We use the
Type-1 discretization proposed in [15] to ensure compatibility of our method
with the prescribed trained steps of such model.

For ImageNet-64, we both use the baseline and the optimized pre-trained models
given in [12]. We note that the baseline is trained on the iDDPM class of
model [19], which actually uses different preconditioning and thus the change
of variables compared to the optimized model. The Figure 1 (c) was obtained
using the EDM-preconditioned checkpoint777https://nvlabs-fi-
cdn.nvidia.com/edm/pretrained/edm-imagenet-64x64-cond-adm.pkl. The added noise
settings of SEEDS-3 solver were not subject to a grid-search optimization
procedure. The chosen hyper-parameters were $S_{\operatorname{churn}}=11$,
$S_{\operatorname{noise}}=1.003$, $S_{\operatorname{min}}=0.05$, and
$S_{\operatorname{max}}=15$ but we are confident that this configuration can
be optimized to further improve our results.

###  E.2 Noise vs. Data Prediction approaches

In Appendix B of [16], the authors compare DPM-Solver2 and DPM-Solver++(2S),
which amounts on comparing in our framework the difference between the
obtained exponential integrators for the PFO in the noise and data prediction
regimes to detect exactly a coefficient on the non-linear term that is absent
in the noise prediction regime. The term they find corresponds exactly to the
difference between applying the variation of constants formula before (instead
of after) replacing the score function with the desired neural network. In
Tab. 2 we report both data and noise prediction SEEDS-3. At low NFEs the DP
approach gives better results but stabilizes in high NFEs at a FID score that
is worse than the one the NP approach reaches.

Table 2: Comparison between noise prediction $F_{\theta,t}$ and data
prediction $D_{\theta,t}$ modes of SEEDS-3 on CIFAR-10 (VP uncond. cont.).
Sampling method / NFE | 9 | 30 | 60 | 90 | 150 | 165 | 180
---|---|---|---|---|---|---|---
SEEDS-3 data prediction | 60.75 | 22.42 | 12.47 | 2.95 | 2.51 | 2.54 | 2.55
SEEDS-3 noise prediction | 471.29 | 288.20 | 33.92 | 3.76 | 2.40 | 2.39 | 2.48

###  E.3 Low vs. High stage Solvers

Similar to DPM-Solver [15], the FID scores in Tab. 3 and 4 show that at low
NFEs, higher stage methods performs more poorly while at higher NFEs, DPM-
Solver-3 and SEEDS-3 are better than their counterparts with 1 and 2 stages.

Table 3: FID comparison between SEEDS (Ours) and DPM-Solver for low NFEs on CIFAR-10 VP uncond. discrete. We recomputed the DPM-Solver score using the "non-deep" model while [15] reports results for the "deep" architecture. The symbol ⋆ is used when using 1 NFE more and † when using 1 NFE less because the given NFE cannot be divided by 2 or 3. This corresponds to Figure 1 (a). Method / NFE | 10 | 12 | 15 | 20 | 30 | 40 | 50 | 100
---|---|---|---|---|---|---|---|---
DPM-Solver-1 | 22.90 | 17.73 | 13.36 | 9.78 | 6.87 | 5.77 | 5.17 | 4.22
DPM-Solver-2 | 12.22 | 6.52 |  ⋆4.55 | – | 3.75 | 3.68 | 3.64 | 3.60
DPM-Solver-3 |  †66.92 | 9.72 | 5.32 |  ⋆3.83 | 3.66 | – |  ⋆3.61 |  †3.58
SEEDS-1 | 303.48 | 239.79 | 279.84 | 192.68 | 84.78 | 45.26 | 28.18 | 8.24
SEEDS-2 | 481.09 | 473.48 |  ⋆430.98 | 305.88 | 223.01 | 51.41 | 11.10 | 3.19
SEEDS-3 |  †483.04 | 482.19 | 479.63 |  ⋆462.61 | 280.48 |  †247.44 |  ⋆62.62 |  †3.53

Table 4: FID comparison between SEEDS (Ours) and DPM-Solver for high NFEs on CIFAR-10 VP uncond. discrete. We recomputed the DPM-Solver score using the "non-deep" model while [15] reports results for the "deep" architecture. The symbol ⋆ is used when using 1 NFE more because the given NFE cannot be divided by 2 or 3. This corresponds to Figure 1 (a). Method / NFE | 150 | 200 | 300 | 510
---|---|---|---|---
DPM-Solver-3 | 3.59 |  ⋆3.58 | – | 3.58
SEEDS-1 | – | 4.07 | 3.40 | –
SEEDS-3 | 3.12 |  ⋆3.08 | 3.14 | 3.24

###  E.4 Deterministic vs. Stochastic Solver

Deterministic solvers as DPM-Solver [15] are fast and well-adapted to
applications in which speed is the most concern. As shown in [15], Table 3 and
4, DPM-Solver converges to a local minimum at early steps and cannot be
improved in large NFEs. Moreover, preconditioned deterministic solver in EDM
gives optimal quality on unconditional CIFAR-10 [12, Fig. 5 (b)]. However, for
more complicated datasets as ImageNet64, the stochasticity indeed helps
improve the samples quality [12, Fig. 5 (c)]. We can consider the random noise
as a corrector that approaches a better local or even global minimum. At
$S_{\text{churn}}=0$, our SEEDS-3 with stochasticity gives lower FID score
than EDM deterministic Heun (see Fig. 1 (c)). SEEDS-3 also reaches the best
quality prior to the number steps needed in Euler Maruyama and other solvers
as in Table 5.

Table 5: FID comparison between SEEDS-3 (Ours) and other solvers on CIFAR-10
VP uncond. discrete. The symbol ⋆ is used when using 1 NFE more and † when
using 1 NFE less because the given NFE cannot be divided by 2 or 3. This
corresponds to Figure 1 (b).
Method / NFE | 10 | 12 | 15 | 20 | 50 | 200 | 1000
---|---|---|---|---|---|---|---
Euler-Maruyama | 278.67 | 246.29 | 197.63 | 137.34 | 32.63 | 4.03 | 3.16
Analytic DDPM | 35.03 | 27.69 | 20.82 | 15.35 | 7.34 | 4.11 | 3.84
Analytic DDIM | 14.74 | 11.68 | 9.16 | 7.20 | 4.28 | 3.60 | 3.86
DDIM | 13.58 | 11.02 | 8.92 | 6.94 | 4.73 | 4.07 | 3.95
DPM-Solver-3 |  †6.92 | 9.72 | 5.32 |  ⋆3.83 |  ⋆3.61 |  ⋆3.58 | –
SEEDS-3 |  †483.04 | 482.19 | 479.63 |  ⋆462.61 |  ⋆62.62 |  ⋆3.08 | –

###  E.5 Hardware configuration

During the experiments, we used three Linux-based servers with 60GB memory
each and 4 GPUs NVIDIA V100 32GB, 4 GPUs NVIDIA V100 16GB, and 2 GPUs NVIDIA
V100 32GB, respectively. Table 6 shows the detail of the configuration
utilized for each experiment. We noted that when using the 4 GPUs
configuration, the FID results were slightly lower (around 2%), even after
using a stacked fixed random seed. We run the experiments multiple times and
reported the minimum FID value each time.

Table 6: Details of GPUs utilized during the experiments. Experiment | Model | Number | GPU Size
---|---|---|---
CIFAR-10 continuous | Nvidia V100 | 4 | 16 GB
CIFAR-10 discrete | Nvidia V100 | 2 | 32 GB
FFHQ64 | Nvidia V100 | 4 | 16 GB
CelebA64 | Nvidia A100 | 2 | 16 GB
ImageNet64 | Nvidia V100 | 4 | 32 GB

###  E.6 Licences

Pre-trained models:

  * •

CIFAR-10 models by [24]: Apache V2.0 license

  * •

FFHQ-64 model by [12]: Creative Commons Attribution-NonCommercial-ShareAlike
4.0 International License.

  * •

CelebA-64 model by [23]: Apache V2.0 license

  * •

ImageNet-64 model by [12]: Apache V2.0 license

  * •

Inception-v3 model by [25]: Apache V2.0 license

###  E.7 Supporting samples

In this subsection, we report the image grids supporting our claims in the
Experiments section.

NP

DP

NFE=9  NFE=21 NFE=30 NFE=90 NFE=129 NFE=165

Figure 3: Samples on CIFAR-10 from low to high NFEs by SEEDS-3 in Noise
Prediction (NP) and Data Prediction (DP) modes, using conditional VP
continuous baseline model [12].

NP

DP

NFE=9  NFE=21 NFE=30 NFE=90 NFE=129 NFE=165

Figure 4: Samples on CIFAR-10 from low to high NFEs by SEEDS-3 in Noise
Prediction (NP) and Data Prediction (DP) modes, using unconditional VP
continuous baseline model [12].

SEEDS-3

DPM-Solver-3

NFE=9  NFE=21 NFE=30 NFE=60  NFE=90 NFE=150

Figure 5: Samples on CelebA-64 from low to high NFEs by SEEDS-3 and DPM-
Solver-3, using pre-trained model from [23].

NP

DP

NFE=9  NFE=21 NFE=30 NFE=60  NFE=90 NFE=150

Figure 6: Samples on FFHQ-64 from low to high NFEs by SEEDS-3 in Noise
Prediction (NP) and Data Prediction (DP) modes, using unconditional VP
continuous baseline model [12].

NP

DP

NFE=30  NFE=60 NFE=90 NFE=150  NFE=210 NFE=270

Figure 7: Samples on ImageNet-64 from low to high NFEs in Noise Prediction
(NP) and Data Prediction (DP) modes, using conditional EDM optimized model
[12].

NP

DP

NFE=30  NFE=60 NFE=90 NFE=150  NFE=210 NFE=270

Figure 8: Samples on ImageNet-64 from low to high NFEs in Noise Prediction
(NP) and Data Prediction (DP) modes, using conditional EDM baseline model
[12].

Figure 9: Example of samples on CelebA-64 generated by SEEDS-3 in 90 NFEs,
using pre-trained model from [23].

(a) DPM-Solver-3 at 44 NFE.

(b) SEEDS-3 at 201 NFE.

Figure 10: Visual sample quality comparison between DPM-Solver-3 and SEEDS-3
using their optimal settings and unconditional CIFAR-10 discrete model [24].

(a) DPM-Solver-3 at 80 NFE.

(b) EDM at 80 NFE.

(c) SEEDS-3 at 150 NFE.

Figure 11: Visual sample quality comparison between DPM-Solver-3, EDM and
SEEDS-3 using their optimal settings and unconditional VP FFHQ-64 continuous
model [12].

(a) DPM-Solver-3 at 30 NFE.

(b) EDM at 90 NFE.

(c) SEEDS-3 at 150 NFE.

Figure 12: Visual sample quality comparison between DPM-Solver-3, EDM and
SEEDS-3 using their optimal settings and conditional VP CIFAR-10 baseline
continuous model [12].  Figure 13: Example of samples on LSUN-Bedroom-256 by
SEEDS-3 in 201 NFEs, using pre-trained model from [5].