Diffusion models for Gaussian distributions: Exact solutions and Wasserstein errors

## Abstract

Abstract Diffusion or score-based models recently showed high performance in image generation.
They rely on a forward and a backward stochastic differential equations (SDE). The sampling of a data distribution is achieved by solving numerically the backward SDE or its associated flow ODE.
Studying the convergence of these models necessitates to control four different types of error: the initialization error, the truncation error, the discretization and the score approximation.
In this paper, we study theoretically the behavior of diffusion models and their numerical implementation when the data distribution is Gaussian.
In this restricted framework where the score function is a linear operator, we derive the analytical solutions of the backward SDE and the probability flow ODE.
We prove that these solutions and their discretizations are all Gaussian processes, which allows us to compute exact Wasserstein errors induced by each error type for any sampling scheme.
Monitoring convergence directly in the data space instead of relying on Inception features, our experiments show that the recommended numerical schemes from the diffusion models literature are also the best sampling schemes for Gaussian distributions.

## 1 Introduction

Over the last five years, diffusion models have proven to be a highly efficient and reliable framework for generative modeling [, , , , , ].
First introduced as a discrete process, Denoising Diffusion Probabilistic Models (DDPM) [] can be studied as a reversal of a continuous Stochastic Differential Equation (SDE) [].
A forward SDE progressively transforms the initial data distribution by adding more and more noise as time grows.
Then, the reversal of this process, called backward SDE, allows us to approximately sample the data distribution starting from Gaussian white noise.
Moreover, the SDE is associated with an Ordinary Differential Equations (ODE) called probability flow [].
This flow preserves the same marginal distributions as the backward SDE and provides another way to sample the score-based generative model.

An important issue about diffusion models is the theoretical guarantees of convergence of the model:
How close to the data distribution the generated distribution is?
There are four main sources of errors to study for deriving theoretical guarantees for diffusion models: (a) the initialization error is induced when approximating the marginal distribution at the end of the forward process by a standard Gaussian distribution. (b) The discretization error comes from the resolution of the SDE or the ODE by a numerical method. (c) The truncation error occurs because the backward time integration is stopped at a small time ${\epsilon}>0$ to avoid numerical instabilities due to ill-defined score function near the origin. (d) The score approximation error accounts for the mismatch between the ideal score function and the one given by the network trained using denoising score-matching.

Despite these numerous sources of errors,
a lot of numerical and theoretical research has been led to assess the generative capacity of diffusion models.
Several articles [, ]
provide strong experimental studies for the choices of sampling parameters.
On the theoretical side, several works derive upper bounds on the 1-Wasserstein or TV distance between the data and the model distributions by making assumptions on the $L^{2}$-error between the ideal and learned score functions and on the compacity of the support of the data [, , , , , ], eventually under an additional manifold assumption [, , ].
Yet, on one hand, to the best of our knowledge, the derived theoretical bounds mostly rely on worst case scenario and are not tight enough to explain the practical efficiency of diffusion models.
On the other hand, numerical considerations mostly rely on Inception feature distributions through the FID metric [].

Ideally, given a data distribution of interest, one would like to have an adapted estimation of the discrepancy between the data and the diffusion model samples, thus enabling adaptive hyperparameter selection for the sampling procedure.
As a first step towards reaching this goal, in the present work we study diffusion models applied to Gaussian data distributions.
While this setting has a priori no practical interest, since simulating Gaussian variates does not require a diffusion model, it provides a large parametric family of distributions for which the errors involved in diffusion model sampling can be completely understood.

When restricting the data distribution to be Gaussian, the resulting score function is a simple linear operator. Exploiting this specificity allows us to derive the following contributions under the assumption that the data distribution is Gaussian:

- •
We give the exact solutions for both the backward SDE and the probability flow ODE.
- •
We fully describe the Gaussian processes that occur when using classical sampling discretization schemes.
- •
We derive exact 2-Wasserstein errors for the corresponding sample distributions and are able to assert for the influence of each error type on these errors, as illustrated by Figure 1.

Our theoretical study allows for an analytical evaluation of any numerical sampler, either stochastic or deterministic.
In particular, it confirms the strength of best practice scheme such as Heun’s method for the ODE flow [].
We provide our source code that can be applied to any Gaussian data distribution of interest and gives insight to calibrate parameters of a diffusion sampling algorithm (see supplementary material).

While our theoretical analysis relies on an exactly known score function,
we conduct additional experiments to assess for the influence of the score approximation error.
Surprisingly, in the context of texture synthesis, we show that with a score neural network trained for modeling a specific Gaussian micro-texture a stochastic Euler-Maruyama sampler is more faithful to the data distribution than Heun’s method, thus highlighting the importance of the score approximation error in practical situations.

Plan of the paper: First, we recall in Section 2 the continuous framework for SDE-based diffusion models.
Section 3 presents our main theoretical results detailing the exact backward SDE and probability flow ODE solutions when supposing the data distribution to be Gaussian.
Section 4 gives explicit Wasserstein error formulas when sampling the corresponding processes,
yielding to an ablation study for comparing the influence of each error type on several sampling schemes.
In Section 5, we study numerically a special case of Gaussian distribution for texture synthesis in order to evaluate the influence of the score approximation error occurring with a standard network architecture. Finally, we address discussion and limitations of our framework in Section 6.

Figure: (a) Initialization error along the integration time
Refer to caption: x1.png

## 2 Preliminaries: Score-based models through diffusion SDEs

This preliminary section follows the seminal work of Song *et al.* [] and introduces specific notation to differentiate the exact backward process and the generative backward process obtained when starting from a white noise.
Given a target distribution $p_{\rm{data}}$ over $\mathbb{R}^{d}$, the forward diffusion process is the following variance preserving SDE

$$ $d{\bm{x}}_{t}=-\beta_{t}{\bm{x}}_{t}dt+\sqrt{2\beta_{t}}d\bm{w}_{t},\quad 0 \leq t\leq T,\quad{\bm{x}}_{0}\sim p_{\rm{data}}$ (1) $$

where $\left(\bm{w}_{t}\right)_{t\geq 0}$ is a $d$-dimensional Brownian motion and $\beta$ is a positive weight function.
The distribution $p_{\rm{data}}$ is noised progressively and the function $\beta$ is the variance of the added noise by time unit.
We denote by $p_{t}$ the density of $\left({\bm{x}}_{t}\right)$ for $t>0$ since $p_{\rm{data}}$ can be supported on a lower-dimensional manifold [].
The SDE is designed so that $p_{T}$ is close to the Gaussian standard distribution that we denote
$\mathcal{N}_{0}$ in whole paper.
Under some assumptions on the distribution $p_{\rm{data}}$ [], the backward process $\left({\bm{x}}_{T-t}\right)_{0\leq t\leq T}$ verifies the backward SDE

$$ $d{\bm{y}}_{t}=\beta_{T-t}({\bm{y}}_{t}+2\nabla_{\bm{y}}\log p_{T-t}({\bm{y}}_{ t}))dt+\sqrt{2\beta_{T-t}}d\bm{w}_{t},\quad 0\leq t<T,\quad{\bm{y}}_{0}\sim p_ {T}.$ (2) $$

The objective is now to solve this reverse equation to sample ${\bm{y}}_{T}\sim p_{\mathrm{data}}$.
However, the distribution $p_{T}$ is in general not known, and image(^1^11Although we may refer to data as images, our analysis is fully general and applies to any vector-valued diffusion model.) generation is achieved by sampling

$$ $d{\tilde{\bm{y}}}_{t}=\beta_{T-t}({\tilde{\bm{y}}}_{t}+2\nabla_{\bm{y}}\log p_ {T-t}({\tilde{\bm{y}}}_{t}))dt+\sqrt{2\beta_{T-t}}d\bm{w}_{t},\quad 0\leq t<T, \quad{\tilde{\bm{y}}}_{0}\sim\mathcal{N}_{0}.$ (3) $$

Note that approximating $p_{T}$ by $\mathcal{N}_{0}$ for the initialization ${\bm{y}}_{0}$ makes that the solution of the SDE of Equation 3 is not exactly the target distribution $p_{\rm{data}}$.
An alternative way to approximately sample $p_{\rm{data}}$ is to use that every diffusion process is associated with a deterministic process whose trajectories share the same marginal probability densities $\left(p_{t}\right)_{0<t\leq T}$ as the SDE [].
The deterministic process associated with Equation 2 is

$$ $d{\bm{x}}_{t}=\left[-\beta_{t}{\bm{x}}_{t}-\beta_{t}\nabla_{{\bm{x}}}\log p_{t }({\bm{x}}_{t})\right]dt,\quad 0<t\leq T,\quad{\bm{x}}_{0}\sim p_{\mathrm{data }}.$ (4) $$

This ODE can be solved in reverse-time to sample ${\bm{x}}_{0}$ from ${\bm{x}}_{T}\sim p_{T}$. Given $({\bm{x}}_{t})_{0\leq t\leq T}$ solution of Equation 4, $({\bm{x}}_{T-t})_{0\leq t\leq T}$ is solution of

$$ $d{\bm{y}}_{t}=\left[\beta_{T-t}{\bm{y}}_{t}+\beta_{T-t}\nabla_{{\bm{y}}}\log p _{T-t}({\bm{y}}_{t})\right]dt,\quad 0\leq t<T.$ (5) $$

Again, in practice, the ODE which is considered to achieve image generation is

$$ $d{\widehat{\bm{y}}}_{t}=\left[\beta_{T-t}{\widehat{\bm{y}}}_{t}+\beta_{T-t} \nabla_{{\widehat{\bm{y}}}}\log p_{T-t}({\widehat{\bm{y}}}_{t})\right]dt,\quad 0 \leq t<T,\quad{\widehat{\bm{y}}}_{0}\sim\mathcal{N}_{0},$ (6) $$

where $p_{T}$ is replaced by $\mathcal{N}_{0}$.
As a consequence of this approximation, the property of conservation of the marginals $\left(p_{t}\right)_{0\leq t\leq T}$ does not occur.
We denote by $\left({\tilde{q}}_{t}\right)_{0\leq t\leq T}$, respectively $\left({\widehat{q}}_{t}\right)_{0\leq t\leq T}$, the marginals of $\left({\tilde{\bm{y}}}_{t}\right)_{0\leq t\leq T}$ and $\left({\widehat{\bm{y}}}_{t}\right)_{0\leq t\leq T}$ and ${\tilde{p}}_{t}={\tilde{q}}_{T-t}$, ${\widehat{p}}_{t}={\widehat{q}}_{T-t}$ the marginals of $\left({\tilde{\bm{y}}}_{T-t}\right)_{0\leq t\leq T}$ and $\left({\widehat{\bm{y}}}_{T-t}\right)_{0\leq t\leq T}$ such that ${\tilde{p}}_{t}$ and ${\widehat{p}}_{t}$ are approximations of $p_{t}$.

## 3 Exact SDE and ODE solutions

Our approach relies on deriving explicit solutions to the various SDE and ODE.
We begin with the forward SDE in full generality obtained in applying the variation of constants (see the proof in Appendix B.1). This resolution also provides an ODE verified by the covariance matrix of ${\bm{x}}_{t}$, that we denote ${\bm{\Sigma}}_{t}=\mathrm{Cov}({\bm{x}}_{t})$.

###### Proposition 1 (Solution of the forward SDE) .

The strong solution of Equation 1 can be written as:

$$ ${\bm{x}}_{t}=e^{-B_{t}}{\bm{x}}_{0}+{\bm{\eta}}_{t},\quad 0\leq t\leq T,$ (7) $$

where $B_{t}=\int_{0}^{t}\beta_{s}ds$ and ${\bm{\eta}}_{t}=e^{-B_{t}}\int_{0}^{t}e^{B_{s}}\sqrt{2\beta_{s}}d\bm{w}_{s}$ is a Gaussian process independent of ${\bm{x}}_{0}$ whose covariance matrix is $(1-e^{-2B_{t}})\bm{I}$. Consequently, the covariance matrix ${\bm{\Sigma}}_{t}$ of ${\bm{x}}_{t}$ is

$$ ${\bm{\Sigma}}_{t}=e^{-2B_{t}}{\bm{\Sigma}}+(1-e^{-2B_{t}})\bm{I}.$ (8) $$

where ${\bm{\Sigma}}$ is the covariance matrix of ${\bm{x}}_{0}\sim p_{\rm{data}}$.
Futhermore, ${\bm{\Sigma}}_{t}$ is invertible for $t>0$ and verifies the matrix-valued ODE

$$ $d{\bm{\Sigma}}_{t}=2\beta_{t}(\bm{I}-{\bm{\Sigma}}_{t})dt,\quad 0<t\leq T.$ (9) $$

For a general data distribution $p_{\rm{data}}$, solving the backward SDE in infeasible, the main reason being that the expression of the score function to integrate is unknown.
To circumvent this obstacle, we now suppose that the data distribution is Gaussian.

###### Assumption 1 (Gaussian assumption) .

$p_{\rm{data}}$ is a centered Gaussian distribution $\mathcal{N}({\bm{0}},{\bm{\Sigma}})$.

Note that ${\bm{\Sigma}}$ may be non-invertible and thus $p_{\rm{data}}$ supported on a strict subspace of $\mathbb{R}^{d}$, a special case of manifold hypothesis.
Consequently, the matrix ${\bm{\Sigma}}_{t}$ is in general only invertible for $t>0$.
Under Gaussian assumption, $({\bm{x}}_{t})$ is a Gaussian process with marginal distribution $p_{t}=\mathcal{N}({\bm{0}},{\bm{\Sigma}}_{t})$ and consequently the score is the linear function

$$ $\nabla\log p_{t}({\bm{x}})=-{\bm{\Sigma}}_{t}^{-1}{\bm{x}},\quad 0<t\leq T.$ (10) $$

Note that the linearity of the diffusion score characterizes Gaussian distributions as detailed by Proposition 5 in Appendix A.

The cornerstone of our work is that under Gaussian assumption we can derive an exact solution of the backward SDE, without supposing that the initial condition is Gaussian.

###### Proposition 2 (Solution of the backward SDE under Gaussian assumption) .

Under Gaussian assumption, the strong solution to the SDE of Equation 2:

$$ $d{\bm{y}}_{t}=\beta_{T-t}({\bm{y}}_{t}+2\nabla_{\bm{y}}\log p_{T-t}({\bm{y}}_{ t}))dt+\sqrt{2\beta_{T-t}}d\bm{w}_{t},\quad 0\leq t<T$ (11) $$

with ${\bm{y}}_{0}$ following any initial distribution can be written as:

$$ ${\bm{y}}_{t}=e^{-(B_{T}-B_{T-t})}{\bm{\Sigma}}_{T-t}{\bm{\Sigma}}^{-1}_{T}{\bm {y}}_{0}+{\bm{\xi}}_{t},\quad 0\leq t\leq T$ (12) $$

where ${\bm{\xi}}_{t}=e^{-(B_{T}-B_{T-t})}{\bm{\Sigma}}_{T-t}\int_{0}^{t}{\bm{\Sigma}
}^{-1}_{T-s}e^{-(B_{T}-B_{T-s})}\sqrt{2\beta_{T-s}}d\bm{w}_{s}$ is a Gaussian process with covariance matrix $\mathrm{Cov}({\bm{\xi}}_{t})={\bm{\Sigma}}_{T-t}-e^{-2(B_{T}-B_{T-t})}{\bm{
\Sigma}}^{2}_{T-t}{\bm{\Sigma}}^{-1}_{T}$. Finally:

$$ $\mathrm{Cov}({\bm{y}}_{t})={\bm{\Sigma}}_{T-t}+e^{-2(B_{T}-B_{T-t})}{\bm{ \Sigma}}_{T-t}^{2}{\bm{\Sigma}}^{-1}_{T}\left({\bm{\Sigma}}^{-1}_{T-t}\mathrm{ Cov}({\bm{y}}_{0}){\bm{\Sigma}}^{-1}_{T}{\bm{\Sigma}}_{T-t}-\bm{I}\right),$ (13) $$

and in particular, if $\mathrm{Cov}({\bm{y}}_{0})$ and ${\bm{\Sigma}}$ commute,

$$ $\mathrm{Cov}({\bm{y}}_{t})={\bm{\Sigma}}_{T-t}+e^{-2(B_{T}-B_{T-t})}{\bm{ \Sigma}}_{T-t}^{2}{\bm{\Sigma}}_{T}^{-2}\left[\mathrm{Cov}({\bm{y}}_{0})-{\bm{ \Sigma}}_{T}\right].$ (14) $$

While not as straightforward as the forward case, the proof also relies on applying the variation of constants and is given in Appendix B.2.
Note that if ${\bm{y}}_{0}$ is correctly initialized at $p_{T}$, $y_{T-t}\sim p_{t}$ at each time $0\leq t\leq T$.
As shown by the following proposition (proved in Appendix B.3),
the flow ODE also has an explicit solution under Gaussian assumption which is related to optimal transport (OT).

###### Proposition 3 (Solution of the ODE probability flow under Gaussian assumption) .

The solution to the reverse-time probability flow ODE of Equation 5:

$$ $d{\bm{y}}_{t}=\left[\beta_{T-t}{\bm{y}}_{t}+\beta_{T-t}\nabla_{{\bm{y}}}\log p _{T-t}({\bm{y}}_{t})\right]dt,\quad 0\leq t<T$ (15) $$

for any $\bm{y}_{0}$ is:

$$ ${\bm{y}}_{t}={\bm{\Sigma}}^{-1/2}_{T}{\bm{\Sigma}}^{1/2}_{T-t}{\bm{y}}_{0}, \quad 0\leq t\leq T,$ (16) $$

which is the application of the OT map between $p_{T}$ and $p_{T-t}$ to the initial condition $\bm{y}_{0}$. Consequently, the covariance matrix $\mathrm{Cov}({\bm{y}}_{t})$ verifies

$$ $\mathrm{Cov}({\bm{y}}_{t})={\bm{\Sigma}}^{-1/2}_{T}{\bm{\Sigma}}^{1/2}_{T-t} \mathrm{Cov}({\bm{y}}_{0}){\bm{\Sigma}}^{1/2}_{T-t}{\bm{\Sigma}}^{-1/2}_{T}, \quad 0\leq t\leq T,$ (17) $$

and in particular, if $\mathrm{Cov}({\bm{y}}_{0})$ and ${\bm{\Sigma}}$ commute,

$$ $\mathrm{Cov}({\bm{y}}_{t})={\bm{\Sigma}}_{T-t}{\bm{\Sigma}}^{-1}_{T}\mathrm{ Cov}({\bm{y}}_{0}),\quad 0\leq t\leq T.$ (18) $$

Here we must highlight a subtle issue: Whatever the initial distribution of ${\bm{y}}_{0}$ is, the ODE solution consists in applying the OT map between $p_{T}$ and $p_{T-t}$ at time $t$.
If ${\bm{y}}_{0}$ follows $p_{T}$, then ${\bm{y}}_{T-t}\sim p_{t}$ at each time $0\leq t\leq T$.
But since in practice one cannot truly sample $p_{T}$ and uses ${\bm{y}}_{0}\sim\mathcal{N}_{0}$ instead,
the resulting flow is not an OT flow (even though it involves an OT mapping) and the distribution of ${\bm{y}}_{T}$ differs from $p_{\rm{data}}$.

###### Proposition 1 (Solution of the forward SDE) .

###### Assumption 1 (Gaussian assumption) .

###### Proposition 2 (Solution of the backward SDE under Gaussian assumption) .

###### Proposition 3 (Solution of the ODE probability flow under Gaussian assumption) .

#### Links with related work.

Some parts of the previous propositions have been stated in previous work.

The expression of the SDE solution of Proposition 1 is given without proof in [] and the ODE verified by the variance is given in [] (citing []) but it is generalized here to the full covariance matrix (Equation 9).

To the best of our knowledge Proposition 2 is new and is the cornerstone for our analytical and numerical study.
Gaussian mixtures have been studied in the context of diffusion models [, , ] since they also provide an explicit analytical score.
However, solving exactly the backward SDE is not feasible for Gaussian mixtures as far as we know.

The relation between OT and probability flow ODE has been discussed in [, ]. [] show that, in general, the flow ODE solution is not an OT between $p_{\rm{data}}$ and $\mathcal{N}_{0}$ at infinite time $T\rightarrow+\infty$, thus contradicting a conjecture of [].
Yet, they briefly discuss the Gaussian case as special case for which the conjecture is valid. Indeed, [] derive the solution of the flow ODE under Gaussian assumption at infinite time horizon [, Appendix B].
More recently, an expression of the solution of the flow ODE relying on the eigendecomposition of the covariance matrix of the data in Gaussian case is given in [] assuming ${\bm{y}}_{0}\sim\mathcal{N}_{0}$.
None of these works discuss the mismatch between the OT map and the initialization of ${\bm{y}}_{0}$.
Our Proposition 3 highlights that the generated process is not an OT flow due to the initialization error.

## 4 Exact Wasserstein errors

The specificity of the Gaussian case allows us to study precisely the different types of error with the expression of the explicit solution of the backward SDE. In what follows, we designate by Wasserstein distance the 2-Wasserstein distance which is known in closed forms when applied to Gaussian distributions [].
For two centered Gaussians $\mathcal{N}({\bm{0}},{\bm{\Sigma}}_{1})$ and $\mathcal{N}({\bm{0}},{\bm{\Sigma}}_{2})$ such that ${\bm{\Sigma}}_{1},{\bm{\Sigma}}_{2}$ are simultaneously diagonalizable with respective eigenvalues $\left(\lambda_{i,1}\right)_{1\leq i\leq d},\left(\lambda_{i,2}\right)_{1\leq i
\leq d},$

$$ $\operatorname{\text{$\textbf{W}_{2}$}}2(\mathcal{N}({\bm{0}},{\bm{\Sigma}}_{1} ),\mathcal{N}({\bm{0}},{\bm{\Sigma}}_{2}))^{2}=\sum_{1\leq i\leq d}(\sqrt{ \lambda_{i,1}}-\sqrt{\lambda_{i,2}})^{2}$ (19) $$

as used in []. In the literature, the quality of the diffusion models is measured with FID []
which is the $\operatorname{\text{$\textbf{W}_{2}$}}2$-error between Gaussians fitted to the Inception features [] of two discrete datasets.
Here we use the $\operatorname{\text{$\textbf{W}_{2}$}}2$-errors directly in data space, which is more informative and allows us to provide theoretical $\operatorname{\text{$\textbf{W}_{2}$}}2$-errors.
To illustrate our theoretical results, we consider the CIFAR-10 Gaussian distribution, that is, the Gaussian distribution such that ${\bm{\Sigma}}$ is the empirical covariance of the CIFAR-10 dataset.
As shown in Appendix C, images produced by this model are not interesting due to a lack of structure, but the corresponding covariance has the advantage of reflecting the complexity of real data.

#### The initialization error.

As discussed in Sections 2 and 3,
the marginals of
both generative processes ${\tilde{\bm{y}}}$ and ${\widehat{\bm{y}}}$ following respectively Equation 6 and Equation 3
slightly differs from $p_{t}$ due to their common white noise initial condition.
This implies an error that we call the initialization error. The distance between $\left({\tilde{p}}_{t}\right)_{0\leq t\leq T}$, $\left({\widehat{p}}_{t}\right)_{0\leq t\leq T}$ and $\left(p_{t}\right)_{0\leq t\leq T}$ can be explicitly studied in the Gaussian case with the following proposition (proved in Appendix B.4).

###### Proposition 4 (Marginals of the generative processes under Gaussian assumption) .

Under Gaussian assumption,
$\left({\tilde{\bm{y}}}_{t}\right)_{0\leq t\leq T}$ and $\left({\widehat{\bm{y}}}_{t}\right)_{0\leq t\leq T}$
are Gaussian processes.
At each time $t$, ${\tilde{p}}_{t}$ is the Gaussian distribution $\mathcal{N}({\bm{0}},{\tilde{{\bm{\Sigma}}}}_{t})$ with ${\tilde{{\bm{\Sigma}}}}_{t}={\bm{\Sigma}}_{t}+e^{-2(B_{T}-B_{t})}{\bm{\Sigma}}
^{2}_{t}{\bm{\Sigma}}_{T}^{-1}({\bm{\Sigma}}^{-1}_{T}-\bm{I})$ and ${\widehat{p}}_{t}$ is the Gaussian distribution $\mathcal{N}({\bm{0}},{\widehat{{\bm{\Sigma}}}}_{t})$ with ${\widehat{{\bm{\Sigma}}}}_{t}={\bm{\Sigma}}^{-1}_{T}{\bm{\Sigma}}_{t}$.
For all $0\leq t\leq T$, the three covariance matrices ${\bm{\Sigma}}_{t}$, ${\tilde{{\bm{\Sigma}}}}_{t}$ and ${\widehat{{\bm{\Sigma}}}}_{t}$ share the same range.
Furthermore, for all $0\leq t\leq T$,

$$ $\operatorname{\text{$\textbf{W}_{2}$}}2({\tilde{p}}_{t},p_{t})\leq \operatorname{\text{$\textbf{W}_{2}$}}2({\widehat{p}}_{t},p_{t})$ (20) $$

which shows that at each time $0\leq t\leq T$ and in particular for $t=0$ which corresponds to the desired outputs of the sampler, the SDE sampler is a better sampler than the ODE sampler when the exact score is konwn.

In practice the initialization error for the SDE and ODE samplers may vary by several orders of magnitude, as shown for the CIFAR-10 example in Figure 1.(a) (solid lines) which illustrates Equation 20.

###### Proposition 4 (Marginals of the generative processes under Gaussian assumption) .

#### The discretization error.

The implementation of the SDE and the ODE implies to choose a discrete numerical scheme. We propose to study four different schemes presented in Table 1. The classical Euler-Maruyama (EM) is used in [] and the exponential integrator (EI) in [] to sample from the SDE of Equation 3.
The Euler method is the simplest ODE solver and Heun’s scheme is recommended in [] to model the ODE of Equation 6. Under Gaussian assumption, the eigenvalues of the covariance matrices can be computed numerically recursively for each scheme to evaluate the Wasserstein distance.
Indeed, under Gaussian assumption, the score is a linear operator and the discrete schemes lead to linear operations described in Table 1.
Then, a Gaussian initialization for ${\bm{y}}_{0}$ provides a sequence of centered Gaussian processes $(y^{\Delta,\cdot}_{k})_{k}$ and if ${\bm{y}}_{0}$ follows $p_{T}$ or $\mathcal{N}_{0}$, the covariance matrix $\mathrm{Cov}(y^{\Delta,\cdot}_{k})$
admit the same eigenvectors as ${\bm{\Sigma}}$ and we can use Equation 19 to compute Wasserstein distances. Let us illustrate the computation of the eigenvalues with the EM scheme.
Denoting $\left(\lambda_{i,t}\right)_{1\leq i\leq d}$ the eigenvalues of ${\bm{\Sigma}}_{t}$ and $\big{(}\lambda_{i,k}^{\Delta,\text{EM}}\big{)}_{1\leq i\leq d}$ the eigenvalues of the covariance matrix of the Euler-Maruyama discretization of the SDE at the $k$th step, $1\leq k\leq N-1$, the relation verified by these eigenvalues is

$$ $\lambda_{i,k+1}^{\Delta,\text{EM}}=\big{(}1+\Delta_{t}\beta_{T-t_{k}}(1-\tfrac {2}{\lambda_{i,T-t_{k}}})\big{)}^{2}\lambda_{i,k}^{\Delta,\text{EM}}+2\Delta_{ t}\beta_{T-t_{k}},1\leq i\leq d,0\leq k\leq N-2$ (21) $$

with initialization $\lambda_{i,0}^{\Delta,\text{EM}}=\left\{\begin{array}[]{ll}1&\mbox{if ${\bm{y}
}_{T}$ is initialized at $\mathcal{N}_{0}$ }\\
\lambda_{i,T}&\mbox{if ${\bm{y}}_{T}$ is initialized at $p_{T}$ }\\
\end{array}\right.\quad 1\leq i\leq d.$. More detailed computations for EM and formulas for other schemes are presented in Appendix D.
For each scheme, we recursively compute the eigenvalues at each time discretization and present the observed Wasserstein distance in Figure 1.(a).
We can observe that Heun’s method provide the lower Wasserstein distance, followed by EM, EI and the Euler scheme. Note that the discrete schemes does not preserve the range of the covariance matrix, contrary to the continuous formulas. This explains the fact that the Wasserstein distance increases at the final step.

**Table 1: Stochastic and deterministic discretization schemes. EM and EI disctretize the backward SDE of Equation 3, Euler and Heun schemes discretize of the probability flow ODE of Equation 6 with a regular time schedule $\left(t_{k}\right)_{0\leq k\leq N}$ with stepsize $\Delta_{t}=\tfrac{T}{N}$**
| SDE schemes | Euler-Maruyama (EM) | $\left\{\begin{array}[]{rl}{\tilde{\bm{y}}^{\Delta,\text{EM}}}_{0}&\sim\mathcal<br>{N}_{0}\\<br>{\tilde{\bm{y}}^{\Delta,\text{EM}}}_{k+1}&={\tilde{\bm{y}}^{\Delta,\text{EM}}}<br>_{k}+\Delta_{t}\beta_{T-t_{k}}\left({\tilde{\bm{y}}^{\Delta,\text{EM}}}_{k}-2{<br>\bm{\Sigma}}_{T-t_{k}}^{-1}{\tilde{\bm{y}}^{\Delta,\text{EM}}}_{k}\right)+<br>\sqrt{2\Delta_{t}\beta_{T-t_{k}}}{\bm{z}}_{k},\leavevmode\nobreak\ {\bm{z}}_{k<br>}\sim\mathcal{N}_{0}\end{array}\right.$<br>(22) |
| --- | --- | --- |
| Exponential integrator (EI) | $\begin{array}[]{l}\left\{\begin{array}[]{rl}{\tilde{\bm{y}}^{\Delta,\text{EI}}<br>}_{0}&\sim\mathcal{N}_{0}\\<br>{\tilde{\bm{y}}^{\Delta,\text{EI}}}_{k+1}&={\tilde{\bm{y}}^{\Delta,\text{EI}}}<br>_{k}+\gamma_{1,k}\left({\tilde{\bm{y}}^{\Delta,\text{EI}}}_{k}-2{\bm{\Sigma}}^<br>{-1}_{T-t_{k}}{\tilde{\bm{y}}^{\Delta,\text{EI}}}_{k}\right)+\sqrt{2\gamma_{2,<br>k}}z_{k},\leavevmode\nobreak\ {\bm{z}}_{k}\sim\mathcal{N}_{0}\end{array}\right<br>.\\<br>\text{where}\leavevmode\nobreak\ \gamma_{1,k}=\exp(B_{T-t_{k}}-B_{T-t_{k+1}})-<br>1\leavevmode\nobreak\ \text{and}\leavevmode\nobreak\ \gamma_{2,k}=\frac{1}{2}(<br>\exp(2B_{T-t_{k}}-2B_{T-t_{k+1}})-1)\end{array}$<br>(23) |  |
| ODE schemes | Explicit Euler | $\left\{\begin{array}[]{rl}{\widehat{\bm{y}}}^{\Delta,\text{Euler}}_{0}&\sim<br>\mathcal{N}_{0}\\<br>{\widehat{\bm{y}}}^{\Delta,\text{Euler}}_{k+1}&={\widehat{\bm{y}}}^{\Delta,<br>\text{Euler}}_{k}+\Delta_{t}f(t_{k},{\widehat{\bm{y}}}^{\Delta,\text{Euler}}_{<br>k})\quad\text{with}\leavevmode\nobreak\ f(t,{\bm{y}})=\beta_{T-t}{\bm{y}}-<br>\beta_{T-t}{\bm{\Sigma}}_{T-t}^{-1}{\bm{y}}\end{array}\right.\quad\quad\quad\quad\thinspace$<br>(24) |
| Heun’s method | $\left\{\begin{array}[]{rl}{\widehat{\bm{y}}}^{\Delta,\text{Heun}}_{0}&\sim<br>\mathcal{N}_{0}\\<br>{\widehat{\bm{y}}}^{\Delta,\text{Heun}}_{k+1/2}&={\widehat{\bm{y}}}^{\Delta,<br>\text{Heun}}_{k}+\Delta_{t}f(t_{k},{\widehat{\bm{y}}}^{\Delta,\text{Heun}}_{k}<br>)\quad\text{with}\leavevmode\nobreak\ f(t,{\bm{y}})=\beta_{T-t}{\bm{y}}-\beta_<br>{T-t}{\bm{\Sigma}}_{T-t}^{-1}{\bm{y}}\\<br>{\widehat{\bm{y}}}^{\Delta,\text{Heun}}_{k+1}&={\widehat{\bm{y}}}^{\Delta,<br>\text{Heun}}_{k}+\frac{\Delta_{t}}{2}\left(f(t_{k},{\widehat{\bm{y}}}^{\Delta,<br>\text{Heun}}_{k})+f(t_{k+1},{\widehat{\bm{y}}}^{\Delta,\text{Heun}}_{k+1/2})<br>\right)\end{array}\right.\quad\quad\quad\thinspace\thinspace$<br>(25) |  |

#### The truncation error.

As discussed in [], it is preferable to study the backward process on $[\varepsilon,T]$ instead of $[0,T]$ because the score is a priori not defined for $t=0$, which occurs in our case if ${\bm{\Sigma}}$ is not invertible.
This approximation is called the truncation error. As a consequence, even without initialization error, the backward process leads to $p_{\varepsilon}$ and not $p_{0}$.
Under Gaussian assumption, it is possible to explicit this error with the expression given in Proposition 3 and 2 as done in Figure 1.(b) for both continuous and numerical solutions.
For the standard practice truncation time $\varepsilon=10^{-3}$ [, ], all numerical schemes have an error close to the corresponding continuous solution.
Using a lower $\varepsilon$ value is only relevant for the continuous SDE solution.

**Table 2: Ablation study of Wasserstein errors for the CIFAR-10 Gaussian. For a given discretization scheme, the table presents the Wasserstein distance associated with the truncation error for different values of $\varepsilon$. The columns $p_{T}$ and $\mathcal{N}_{0}$ show the influence of the initialization error. The continuous column corresponds to the continuous SDE or ODE linked with the scheme (identical values for EM, EI and Euler, Heun) and a given number of integration steps $N$.**
|  |  | Continuous | $N=50$ | $N=250$ | $N=500$ | $N=1000$ |  |  |  |  |  |
| --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- |
|  |  | $p_{T}$ | $\mathcal{N}_{0}$ | $p_{T}$ | $\mathcal{N}_{0}$ | $p_{T}$ | $\mathcal{N}_{0}$ | $p_{T}$ | $\mathcal{N}_{0}$ | $p_{T}$ | $\mathcal{N}_{0}$ |
| EM | $\varepsilon=0$ | 0 | 6.7E-04 | 4.78 | 4.78 | 0.65 | 0.66 | 0.32 | 0.32 | 0.16 | 0.16 |
| $\varepsilon=10^{-5}$ | 4.1E-03 | 4.2E-03 | 4.77 | 4.77 | 0.66 | 0.66 | 0.32 | 0.32 | 0.16 | 0.16 |  |
| $\varepsilon=10^{-4}$ | 0.03 | 0.03 | 4.76 | 4.76 | 0.66 | 0.66 | 0.32 | 0.32 | 0.17 | 0.17 |  |
| $\varepsilon=10^{-3}$ | 0.18 | 0.18 | 4.68 | 4.68 | 0.70 | 0.70 | 0.40 | 0.40 | 0.27 | 0.27 |  |
| EI | $\varepsilon=0$ | 0 | 6.7E-04 | 2.81 | 2.81 | 0.57 | 0.57 | 0.30 | 0.30 | 0.16 | 0.16 |
| $\varepsilon=10^{-5}$ | 4.1E-03 | 4.2E-03 | 2.81 | 2.81 | 0.57 | 0.57 | 0.30 | 0.30 | 0.16 | 0.16 |  |
| $\varepsilon=10^{-4}$ | 0.03 | 0.03 | 2.82 | 2.82 | 0.58 | 0.58 | 0.31 | 0.31 | 0.17 | 0.17 |  |
| $\varepsilon=10^{-3}$ | 0.18 | 0.18 | 2.91 | 2.91 | 0.67 | 0.67 | 0.41 | 0.41 | 0.29 | 0.29 |  |
| Euler | $\varepsilon=0$ | 0 | 0.07 | 1.72 | 1.78 | 0.38 | 0.44 | 0.20 | 0.26 | 0.10 | 0.17 |
| $\varepsilon=10^{-5}$ | 4.1E-03 | 0.07 | 1.72 | 1.78 | 0.38 | 0.44 | 0.20 | 0.26 | 0.10 | 0.17 |  |
| $\varepsilon=10^{-4}$ | 0.03 | 0.08 | 1.72 | 1.78 | 0.38 | 0.44 | 0.20 | 0.26 | 0.11 | 0.17 |  |
| $\varepsilon=10^{-3}$ | 0.18 | 0.19 | 1.73 | 1.79 | 0.42 | 0.48 | 0.27 | 0.32 | 0.21 | 0.25 |  |
| Heun | $\varepsilon=0$ | 0 | 0.07 | - | - | - | - | - | - | - | - |
| $\varepsilon=10^{-5}$ | 4.1E-03 | 0.07 | 23.42 | 23.42 | 2.86 | 2.87 | 1.05 | 1.06 | 0.37 | 0.38 |  |
| $\varepsilon=10^{-4}$ | 0.03 | 0.08 | 4.68 | 4.68 | 0.43 | 0.44 | 0.12 | 0.14 | 0.03 | 0.08 |  |
| $\varepsilon=10^{-3}$ | 0.18 | 0.19 | 0.58 | 0.59 | 0.13 | 0.15 | 0.16 | 0.18 | 0.17 | 0.19 |  |

#### Ablation study.

We propose in Table 2 an ablation study to monitor the magnitude of each error and their accumulation for various sampling schemes for the CIFAR-10 example.
In accordance with Proposition 4,
the initialization error influences the ODE schemes, while SDE schemes are not affected.
Schemes having a sufficient number of steps are not sensitive to the truncation error for $\varepsilon<10^{-3}$, except Heun’s scheme which is unstable near to origin.
The discretization error is the more important approximation but it becomes very low for a sufficient number of steps. The lower Wasserstein error is provided by Heun’s method with $1000$ steps, $\varepsilon=10^{-4}$.
As [], our conclusions lead to the choice of Heun’s scheme as the go-to method.

Figure: (a) Initialization error at final time
Refer to caption: x3.png

#### Influence of eigenvalues.

The above observations and conclusions are observed on the CIFAR-10 Gaussian.
However, in general, they depend on the eigenvalues of the covariance matrix ${\bm{\Sigma}}$.
Indeed, as seen in Equation 19, the Wasserstein distance is separable and each eigenvalue contributes to increase it.
In Figure 2, we evaluate the contribution of each eigenvalue by plotting $\lambda\mapsto|\sqrt{\lambda}-\sqrt{\lambda^{\text{scheme}}}|$ for each scheme.
Figure 2.(a) demonstrates that for the continuous equations, the error increases with the eigenvalues except for a strong decrease for $\lambda=1$.
Besides, as proved in the proof of Proposition 4 (see Appendix B.4), the error for the SDE is always lower than the error for the ODE and we can observe how tight is Equation 20.
Unfortunately, once discretized the stochastic schemes are not as good as the continuous solutions.
The EI scheme is the more stable along the range of eigenvalues but in the end it is in general more costly than the others in terms of Wasserstein error.
Without truncation time, Heun’s method fails for low eigenvalues because ${\bm{\Sigma}}$ is not stably invertible.
However, as seen in Figure 2.(b), with a truncation time $\varepsilon=10^{-3}$, Heun’s method is very close to the continuous ODE solution.
This shows that for any Gaussian distribution Heun’s method introduces nearly no additional discretization error, making this scheme the one to favor in practice.
Our code allows for the evaluation of any covariance matrix and the computation of Figure 1 and Table 2 (provided the eigenvalues can be computed, see the supplementary).

## 5 Numerical study of the score approximation

So far our theoretical and numerical study has been conducted under the hypothesis that the score function is known, thus discarding the evaluation of the score approximation.
In practice, for general data distribution, the score function is parameterized by a neural network trained using denoising score-matching.
This learned score function is not perfect and while theoretical studies assume the network to be close to the theoretical one (with uniform or adaptative bounds, see the discussion in []),
such an hypothesis is hard to check in practice, especially in our non compact setting.
Thus, we propose in this section to train a diffusion models on a Gaussian distribution and evaluate numerically the impact of the score approximation.

#### The Gaussian ADSN distribution for microtextures.

So far our running example was the CIFAR-10 Gaussian but we will now turn to another example that produces visually interesting images, namely Gaussian micro-textures.
We consider the asymptotic discrete spot noise (ADSN) distribution [] associated with an RGB texture $\bm{u}\in\mathbb{R}^{3\times{M\times N}}$ which is defined as the stationary Gaussian distribution that has covariance equal the autocorrelation of $\bm{u}$.
More precisely, this distribution is sampled using convolution with a white Gaussian noise []:
Denoting $m\in\mathbb{R}^{3}$ the channelwise mean of $\bm{u}$ and $\mathbf{t}_{c}=\frac{1}{\sqrt{MN}}(\bm{u}_{c}-m_{c})$, $1\leq c\leq 3$, its associated *texton*, for $\bm{w}\sim\mathcal{N}_{0}$ of size $M\times N$
the channelwise convolution ${\bm{x}}=m+\mathbf{t}\star\bm{w}\in\mathbb{R}^{3\times{M\times N}}$ follows $\operatorname{ADSN}(\bm{u})$. This distribution is the Gaussian $\mathcal{N}(m,{\bm{\Sigma}})$.
To deal with zero mean Gaussian, adding the mean $m$ is considered as a post-processing to visualize samples and we study $\mathcal{N}({\bm{0}},{\bm{\Sigma}})$.
The matrix ${\bm{\Sigma}}$ is a well-known convolution matrix [], its eigenvectors and associated eigenvalues can be computed in the Fourier domain, as done in Appendix F.2.
${\bm{\Sigma}}$ admits the eigenvalues $\lambda^{\xi,\operatorname{ADSN}}_{1}=|\widehat{\mathbf{t}}_{1}|^{2}(\xi)+|
\widehat{\mathbf{t}}_{2}|^{2}(\xi)+|\widehat{\mathbf{t}}_{3}|^{2}(\xi),\xi\in
\mathbb{R}^{M\times N}$ and 0 0 with multiplicity $2MN$ and we can conduct the same analysis as before (see Appendix E).
To evaluate if a set of $N_{\text{samples}}$ sampled images is close to the ADSN distribution $p_{\rm{data}}$, we evaluate a problem-specific empirical Wasserstein distance: Supposing that the $N_{\text{samples}}$ are drawn from a Gaussian distribution $p^{\text{\tiny emp.}}=\mathcal{N}({\bm{0}},\bm{\Gamma})$ such that $\bm{\Gamma}$ admits the same eigenvectors as ${\bm{\Sigma}}$, we compute

$$ $\operatorname{\text{$\textbf{W}^{\text{\tiny emp.}}_{2}$}}(p^{\text{\tiny emp. }},p_{\rm{data}})=\sqrt{\sum_{\xi\in\mathbb{R}^{3{M\times N}}}\left(\sqrt{ \lambda^{\xi,\text{emp.}}_{1}}-\sqrt{\lambda^{\xi,\operatorname{ADSN}}_{1}} \right)^{2}+\lambda^{\xi,\text{emp.}}_{2}+\lambda^{\xi,\text{emp.}}_{3}}$ (26) $$

where $(\lambda^{\xi,\text{emp.}}_{i})_{\xi\in\mathbb{R}^{{M\times N}},1\leq i\leq 3}$ are estimators of the eigenvalues of $\bm{\Gamma}$ given in Appendix F.3.

Figure: Figure 3: Texture samples generated with the learned score. First row: original image $\bm{u}$ and its DFT modulus (for all DFT modulus we display the sum of the DFT modulus of the three color channels and apply a logarithmic contrast change). Second row: three samples of $\operatorname{ADSN}(\bm{u})$ with their associated DFT moduli. Third and fourth row: Samples generated with the learned score with EM and Heun’s discretization schemes and their associated DFT moduli. While both schemes use the same learned score function, the generation with Heun’s scheme can produce out-of-distribution samples, as seen with the third sample.
Refer to caption: extracted/5986824/Figures/Samples/ADSN/tissu_256.png

#### Learning the score function.

We train the network using the code(^2^22Code available at [https://github.com/yang-song/score_sde_pytorch](https://github.com/yang-song/score_sde_pytorch)) associated with the paper []. We choose the architecture of DDPM, which is a U-Net described in [], with the parameters proposed for the dataset CelebaHQ256 to deal with the 256$\times$256 ADSN model associated with the top-left image of Figure 3. We use the training procedure corresponding to DDPM cont. in [].
$\beta$ is linear from 0.05 to 10 with $T=1$.
We train over 1.3M iterations, and we generate at each iteration a new batch of ADSN samples.
We implement the stochastic EM and derministic Heun schemes replacing the exact score by its learned version with $N=1000$ steps and a trunction time $\varepsilon=10^{-3}$. We name $p_{\theta}^{\text{EM}}$ and $p_{\theta}^{\text{Heun}}$, the corresponding distributions and present samples in Figure 3.
Both distributions accumulate the four error types.

#### Evaluation of the score approximation.

It is not possible to compute theoretically the Wasserstein distance between $p_{\rm{data}}=\operatorname{ADSN}(\bm{u})$ and $p_{\theta}^{\text{EM}},p_{\theta}^{\text{Heun}}$ due to the non-linearity of the learned score. To compute an empirical Wasserstein error between it, we use Equation 26.
Let us precise that this approximation underestimates the real Wasserstein distance since it wrongly assumes that the distributions $p_{\theta}^{\text{EM}}$, $p_{\theta}^{\text{Heun}}$ are Gaussian with a covariance matrix diagonalizable in the same basis than the covariance matrix ${\bm{\Sigma}}$ of $\operatorname{ADSN}(\bm{u})$.
We complete this dedicated empirical measure with the standard FID.
These metrics are reported in Table 3 where for theoretical distributions that are fast to sample we add the standard deviations computed on $25$ different $50k$-samplings.
For this Gaussian distribution, the score approximation is by far the most impactful source of error.
We observe that the stochastic EM sampling is more resilient to score approximation than the deterministic Heun’s scheme, resulting in out-of-distribution samples (Figure 3).
We may explain this behavior by recalling the results of Proposition 4 that shows that SDE solutions are less sensitive to initialization errors than ODE.
Indeed, adding noise at each iteration tends to mitigate the accumulated errors, and score approximation may be considered as some initialization error ocurring at each step.

**Table 3: Numerical evaluation of the score approximation for a Gaussian microtexture model. For two schemes, the EM discretization of the backward SDE and Heun’s method associated with the flow ODE, the table shows the Wasserstein distance and FID for theoretical and learned distributions. The theoretical $\operatorname{\text{$\textbf{W}_{2}$}}2$ value is computed with explicit formulas, as done in Table 5. The FID and empirical $\operatorname{\text{$\textbf{W}_{2}$}}2$ w.r.t the theoretical distribution are computed on $25$ samplings of 50K images while only one sampling of 50K images is drawn for the parametric distributions (to limit computation time).**
|  | Exact score distribution | Learned score distribution |  |  |  |
| --- | --- | --- | --- | --- | --- |
| $p$ | $\operatorname{\text{$\textbf{W}_{2}$}}2$($p$,$p_{\rm{data}}$) $\downarrow$ | $\operatorname{\text{$\textbf{W}^{\text{\tiny emp.}}_{2}$}}$($p^{\text{\tiny emp.}}$,$p_{\rm{data}}$) $\downarrow$ | FID($p^{\text{\tiny emp.}}$,$p_{\mathrm{data}}^{\text{\tiny emp.}}$) $\downarrow$ | $\operatorname{\text{$\textbf{W}^{\text{\tiny emp.}}_{2}$}}$($p^{\text{\tiny emp.}}_{\theta}$,$p_{\mathrm{data}}^{\text{\tiny emp.}}$) $\downarrow$ | FID($p^{\text{\tiny emp.}}_{\theta}$,$p_{\mathrm{data}}^{\text{\tiny emp.}}$) $\downarrow$ |
| EM | 5.16 | 5.1630$\pm$7E-5 | 0.0891$\pm$8E-4 | 15.6 | 01.02 |
| Heun | 3.73 | 3.7323$\pm$2E-4 | 0.0447$\pm$6E-4 | 56.7 | 19.48 |

## 6 Discussion and limitations

The main limitation of our work is that our theoretical and numerical results are limited to Gaussian distributions.
Resorting to diffusion models for sampling Gaussian distributions is not necessary in practice, rather we use Gaussian distributions as a test case family to provide insight on diffusion models.

A natural extension of this work is to compute error types for more complex distributions (e.g multimodal) such as Gaussian mixtures models (GMM).
However, generalizing our results for these more complex distributions one faces two main difficulties. First, to the best of our knowledge, we are unable to derive exact solutions to the backward SDE or the flow ODE under GMM assumption, even though the score has a known analytical expression [, , ].
Another key feature of this study is to evaluate exactly the Wasserstein error by using Equation 19, strongly relying on the Gaussian assumption.
A closed-form of the Wasserstein distance between two GMMs is not known,
leading to alternative distance definitions for such models [].
Hence, to compare the distributions generated in practice with exact solutions of time continuous equations under GMM assumption, as we do for the Gaussian case, one should solve two open theoretical problems.

## 7 Conclusion

By restricting the analysis of diffusion models to the specific case of Gaussian distributions, we were able to derive exact solutions for both the backward SDE and its associated probability flow ODE.
We demonstrate that regarding the initialization error, the SDE sampler is more resilient than the ODE sampler for Gaussian distributions.
Additionally, we characterized the discrete Gaussian processes arising when discretizing these equations.
This allowed us to provide exact Wasserstein errors for the initialization error, the discretization error, and the truncation error as well as any of their combinations.
This theoretical analysis led to conclude that Heun’s scheme is the best numerical solution, in accordance with empirical previous work [].

To conclude our work we conducted an empirical analysis with a learned score function using standard architecture which showed that the score approximation error may be the most important one in practice.
This suggests that assessing the quality of learned score functions is an important research direction for future work.

## Appendix A Characterization of Gaussian distributions through diffusion models

The following proposition shows that our Gaussian assumption occurs if and only if the score function is linear.

###### Proposition 5 .

The three following propositions are equivalent:

- (i)
${\bm{x}}_{0}\sim\mathcal{N}({\bm{0}},{\bm{\Sigma}})$ for some covariance ${\bm{\Sigma}}$.
- (ii)
$\forall t>0,\nabla_{x}\log p_{t}({\bm{x}})$ is linear w.r.t ${\bm{x}}$.
- (iii)
$\exists t>0,\nabla_{x}\log p_{t}({\bm{x}})$ is linear w.r.t ${\bm{x}}$.

In this case, for $t>0$, $\nabla_{{\bm{x}}}\log p_{t}({\bm{x}})=-{\bm{\Sigma}}^{-1}_{t}{\bm{x}}$, with ${\bm{\Sigma}}_{t}$ defined in Proposition 1.

###### Proof.

$(ii)\Rightarrow(iii)$ is clear.

If $(i)$, for $t>0$, $p_{t}({\bm{x}})=C_{t}\exp\left(-\frac{1}{2}{\bm{x}}^{T}{\bm{\Sigma}}^{-1}_{t}{
\bm{x}}\right)$. Consequently, $\nabla_{{\bm{x}}}\log p_{t}({\bm{x}})=-{\bm{\Sigma}}^{-1}_{t}{\bm{x}}$ and $(i)\Rightarrow(ii)$

If $(iii)$, there exists $A$ such that $\nabla_{{\bm{x}}}\log p_{t}({\bm{x}})=A{\bm{x}}$. Consequently, $p_{t}({\bm{x}})=C_{t}\exp(-\frac{1}{2}{\bm{x}}^{T}A{\bm{x}})$ and ${\bm{x}}_{t}$ is Gaussian. This provides that ${\bm{x}}_{0}=e^{B_{t}}{\bm{x}}_{t}-{\bm{\eta}}_{t}$ is Gaussian and $(iii)\Rightarrow(i)$.
∎

###### Proposition 5 .

###### Proof.

## Appendix B Proofs of Section 3

### B.1 Proposition 1 : Solution of the forward SDE

We aim at solving:

$$ $d{\bm{x}}_{t}=-\beta_{t}{\bm{x}}_{t}dt+\sqrt{2\beta_{t}}d\bm{w}_{t},\quad{\bm{ x}}_{0}\sim p_{\rm{data}}.$ (27) $$

By considering ${\bm{z}}_{t}=e^{B_{t}}{\bm{x}}_{t}$ where $B_{t}=\int_{0}^{t}\beta_{s}ds$,

$$ $d{\bm{z}}_{t}=\beta_{t}e^{B_{t}}{\bm{x}}_{t}+e^{B_{t}}d{\bm{x}}_{t}=\beta_{t}e ^{B_{t}}{\bm{x}}_{t}+e^{B_{t}}(-\beta_{t}{\bm{x}}_{t}dt+\sqrt{2\beta_{t}}d\bm{ w}_{t})=\sqrt{2\beta_{t}}e^{B_{t}}d\bm{w}_{t}.$ (28) $$

Consequently, for $0\leq t\leq T$,

$$ ${\bm{z}}_{t}={\bm{z}}_{0}+\int_{0}^{t}\sqrt{2\beta_{s}}e^{B_{s}}d\bm{w}_{s},{ \bm{z}}_{0}=e^{B_{0}}{\bm{x}}_{0}={\bm{x}}_{0}$ (29) $$

and for $0\leq t\leq T$,

$$ ${\bm{x}}_{t}=e^{-B_{t}}{\bm{z}}_{t}=e^{-B_{t}}{\bm{x}}_{0}+e^{-B_{t}}\int_{0}^ {t}e^{B_{s}}\sqrt{2\beta_{s}}d\bm{w}_{s}=e^{-B_{t}}{\bm{x}}_{0}+{\bm{\eta}}_{t}.$ (30) $$

By Itô’s isometry (see e.g []),

$$ $\text{Var}\left(\int_{0}^{t}e^{B_{s}}\sqrt{2\beta_{s}}d\bm{w}_{s}\right)=\int_ {0}^{t}2\beta_{s}e^{2B_{s}}ds=[e^{2B_{s}}]_{0}^{t}=e^{2B_{t}}-e^{2B_{0}}=e^{2B _{t}}-1$ (31) $$

which provides the covariance matrix of ${\bm{\eta}}_{t}$:

$$ $\mathrm{Cov}\left({\bm{\eta}}_{t}\right)=e^{-2B_{t}}(e^{2B_{t}}-1)\bm{I}=\left (1-e^{-2B_{t}}\right)\bm{I}.$ (32) $$

Because ${\bm{x}}_{0}$ and ${\bm{\eta}}_{t}$ are independent, ${\bm{\Sigma}}_{t}=e^{-2B_{t}}{\bm{\Sigma}}+\left(1-e^{-2B_{t}}\right)\bm{I}$.

And,

$$ $d{\bm{\Sigma}}_{t}=-2\beta_{t}e^{-2B_{t}}({\bm{\Sigma}}-\bm{I})dt=-2\beta_{t} \left[{\bm{\Sigma}}_{t}-\bm{I}\right]dt$ (33) $$

### B.2 Proposition 2 : Solution of the backward SDE under Gaussian assumption

We aim at solving

$$ $d{\bm{y}}_{t}=\beta_{T-t}({\bm{y}}_{t}-2{\bm{\Sigma}}^{-1}_{T-t}{\bm{y}}_{t}) dt+\sqrt{2\beta_{T-t}}d\bm{w}_{t},\quad 0\leq t\leq T$ (34) $$

Denoting $C_{t}=\int_{0}^{t}\beta_{T-s}ds$, by considering ${\bm{z}}_{t}={\bm{\Sigma}}^{-1}_{T-t}e^{C_{t}}{\bm{y}}_{t}$,

$$ $\displaystyle d{\bm{z}}_{t}$ $\displaystyle=e^{C_{t}}{\bm{\Sigma}}^{-1}_{T-t}d{\bm{y}}_{t}+e^{C_{t}}d[{\bm{ \Sigma}}^{-1}_{T-t}]{\bm{y}}_{t}+\beta_{T-t}{\bm{z}}_{t}dt$ (35) $\displaystyle=e^{C_{t}}{\bm{\Sigma}}^{-1}_{T-t}d{\bm{y}}_{t}-e^{C_{t}}{\bm{ \Sigma}}^{-1}_{T-t}d[{\bm{\Sigma}}_{T-t}]{\bm{\Sigma}}^{-1}_{T-t}{\bm{y}}_{t}+ \beta_{T-t}{\bm{z}}_{t}dt\text{ by derivative of the inverse matrix}$ (36) $\displaystyle=e^{C_{t}}{\bm{\Sigma}}^{-1}_{T-t}\left[\beta_{T-t}({\bm{y}}_{t}- 2{\bm{\Sigma}}^{-1}_{T-t}{\bm{y}}_{t})dt+\sqrt{2\beta_{T-t}}d\bm{w}_{t}\right] -2\beta_{T-t}e^{C_{t}}{\bm{\Sigma}}^{-1}_{T-t}\left[{\bm{\Sigma}}_{T-t}-\bm{I} \right]{\bm{\Sigma}}^{-1}_{T-t}{\bm{y}}_{t}dt+\beta_{T-t}{\bm{z}}_{t}dt$ (37) (using Equation 9 ) (38) $\displaystyle=\left[{\bm{\Sigma}}^{-1}_{T-t}e^{C_{t}}\beta_{T-t}({\bm{y}}_{t}- 2{\bm{\Sigma}}^{-1}_{T-t}{\bm{y}}_{t})-\beta_{T-t}{\bm{z}}_{t}+2\beta_{T-t}{ \bm{\Sigma}}^{-1}_{T-t}{\bm{z}}_{t}\right]dt+\sqrt{2\beta_{T-t}}e^{C_{t}}{\bm{ \Sigma}}^{-1}_{T-t}d\bm{w}_{t}$ (39) $\displaystyle=\beta_{T-t}(\bm{I}-2{\bm{\Sigma}}^{-1}_{T-t}){\bm{z}}_{t}dt- \beta_{T-t}{\bm{z}}_{t}dt+2\beta_{T-t}{\bm{\Sigma}}^{-1}_{T-t}{\bm{z}}_{t}dt+e ^{C_{t}}\sqrt{2\beta_{T-t}}{\bm{\Sigma}}^{-1}_{T-t}d\bm{w}_{t}$ (40) $\displaystyle=\sqrt{2\beta_{T-t}}e^{C_{t}}{\bm{\Sigma}}^{-1}_{T-t}d\bm{w}_{t}.$ (41) $$

Consequently,

$$ ${\bm{z}}_{t}={\bm{z}}_{0}+\int_{0}^{t}\sqrt{2\beta_{T-s}}e^{C_{s}}{\bm{\Sigma} }^{-1}_{T-s}d\bm{w}_{s}={\bm{\Sigma}}^{-1}_{T}{\bm{y}}_{0}+\int_{0}^{t}\sqrt{2 \beta_{T-s}}e^{C_{s}}{\bm{\Sigma}}^{-1}_{T-s}d\bm{w}_{s}.$ (43) $$

And,

$$ ${\bm{y}}_{t}=e^{-C_{t}}{\bm{\Sigma}}_{T-t}{\bm{z}}_{t}=e^{-C_{t}}{\bm{\Sigma}} _{T-t}{\bm{\Sigma}}^{-1}_{T}{\bm{y}}_{0}+e^{-C_{t}}{\bm{\Sigma}}_{T-t}\int_{0} ^{t}{\bm{\Sigma}}^{-1}_{T-s}e^{C_{s}}\sqrt{2\beta_{T-s}}d\bm{w}_{s}.$ (44) $$

Finally,

$$ ${\bm{y}}_{t}=e^{-C_{t}}{\bm{\Sigma}}_{T-t}{\bm{\Sigma}}^{-1}_{T}{\bm{y}}_{0}+{ \bm{\xi}}_{t}\quad\text{with}\quad{\bm{\xi}}_{t}=e^{-C_{t}}{\bm{\Sigma}}_{T-t} \int_{0}^{t}{\bm{\Sigma}}^{-1}_{T-s}e^{C_{s}}\sqrt{2\beta_{T-s}}d\bm{w}_{s}.$ (45) $$

By the multidimensional Itô’s isometry,

$$ $\mathrm{Cov}(\int_{0}^{t}{\bm{\Sigma}}^{-1}_{T-s}e^{C_{s}}\sqrt{2\beta_{T-s}}d \bm{w}_{s})=2\int_{0}^{t}e^{2C_{s}}\beta_{T-s}{\bm{\Sigma}}^{-2}_{T-s}ds.$ (46) $$

Now, remark that for $A_{s}=e^{2C_{s}}{\bm{\Sigma}}^{-1}_{T-s}$,

$$ $\displaystyle dA_{s}$ $\displaystyle=2\beta_{T-s}A_{s}ds+e^{2C_{s}}d\left[{\bm{\Sigma}}_{T-s}^{-1}\right]$ (47) $\displaystyle=2\beta_{T-s}A_{s}ds-2\beta_{T-s}e^{2C_{s}}\left[\bm{I}-{\bm{ \Sigma}}^{-1}_{T-s}\right]{\bm{\Sigma}}^{-1}_{T-s}ds\text{ (using \lx@cref{ creftypecap~refnum}{eq:ode_covariance})}$ (48) $\displaystyle=2e^{2C_{s}}\beta_{T-s}{\bm{\Sigma}}^{-2}_{T-s}ds.$ (49) $$

$$ $\mathrm{Cov}\left(\int_{0}^{t}{\bm{\Sigma}}^{-1}_{T-s}e^{C_{s}}\sqrt{\beta_{T- s}}d\bm{w}_{s}\right)=\int_{0}^{t}dA_{s}=\left[A_{s}\right]_{0}^{t}=e^{2C_{t}} {\bm{\Sigma}}^{-1}_{T-t}-{\bm{\Sigma}}^{-1}_{T}.$ (50) $$

Finally, $\mathrm{Cov}({\bm{\xi}}_{t})={\bm{\Sigma}}_{T-t}^{2}\left({\bm{\Sigma}}^{-1}_{
T-t}-e^{-2C_{t}}{\bm{\Sigma}}^{-1}_{T}\right)={\bm{\Sigma}}_{T-t}-e^{-2C_{t}}{
\bm{\Sigma}}_{T-t}^{2}{\bm{\Sigma}}^{-1}_{T}$

We have the final formula considering:

$$ $C_{t}=\int_{0}^{t}\beta_{T-s}ds=\int_{T-t}^{T}\beta_{x}dx=\int_{0}^{T}\beta_{x }dx-\int_{0}^{T-t}\beta_{x}dx=B_{T}-B_{T-t}$ (51) $$

that provides

$$ $\mathrm{Cov}({\bm{y}}_{t})={\bm{\Sigma}}_{T-t}+e^{-2(B_{T}-B_{T-t})}{\bm{ \Sigma}}_{T-t}^{2}{\bm{\Sigma}}^{-1}_{T}\left({\bm{\Sigma}}^{-1}_{T-t}\mathrm{ Cov}({\bm{y}}_{0}){\bm{\Sigma}}^{-1}_{T}{\bm{\Sigma}}_{T-t}-\bm{I}\right).$ (52) $$

In particular, if $\mathrm{Cov}({\bm{y}}_{0})$ and ${\bm{\Sigma}}$ commute,

$$ $\mathrm{Cov}({\bm{y}}_{t})={\bm{\Sigma}}_{T-t}+e^{-2(B_{T}-B_{T-t})}{\bm{ \Sigma}}_{T-t}^{2}{\bm{\Sigma}}^{-1}_{T}\left({\bm{\Sigma}}^{-1}_{T}\mathrm{ Cov}({\bm{y}}_{0})-\bm{I}\right).$ (53) $$

### B.3 Proposition 3 : Solution of the ODE probability flow under Gaussian assumption

As done in [], the matrix ${\bm{\Sigma}}^{1/2}_{t}$ admits a derivative which is $d\left[{\bm{\Sigma}}_{t}^{1/2}\right]=\frac{1}{2}d{\bm{\Sigma}}_{t}{\bm{\Sigma
}}^{-1/2}_{t}$ because it is diagonalisable. Let us check that

$$ ${\bm{y}}_{t}={\bm{\Sigma}}^{-1/2}_{T}{\bm{\Sigma}}_{T-t}^{1/2}{\bm{y}}_{0}$ (54) $$

is solution of the ODE of Equation 5:

$$ $\displaystyle d{\bm{y}}_{t}$ $\displaystyle=-{\bm{\Sigma}}^{-1/2}_{T}\frac{1}{2}d{\bm{\Sigma}}_{T-t}{\bm{ \Sigma}}^{-1/2}_{T-t}{\bm{y}}_{0}$ (55) $\displaystyle={\bm{\Sigma}}^{-1/2}_{T}\left[\beta_{T-t}{\bm{\Sigma}}_{T-t}- \beta_{T-t}\bm{I}\right]{\bm{\Sigma}}^{-1/2}_{T-t}{\bm{y}}_{0}dt\quad(\text{ using \lx@cref{creftypecap~refnum}{eq:ode_covariance}})$ (56) $\displaystyle=\left[\beta_{T-t}{\bm{\Sigma}}_{T-t}-\beta_{T-t}\bm{I}\right]{ \bm{\Sigma}}^{-1}_{T-t}{\bm{\Sigma}}^{-1/2}_{T}{\bm{\Sigma}}^{1/2}_{T-t}{\bm{y }}_{0}dt\quad(\text{by commutativity})$ (57) $\displaystyle=\left[\beta_{T-t}-\beta_{T-t}{\bm{\Sigma}}^{-1}_{T-t}\right]{\bm {y}}_{t}dt$ (58) $\displaystyle=\left[\beta_{T-t}+\beta_{T-t}\nabla_{{\bm{y}}}\log p_{T-t}({\bm{ y}}_{t})\right]{\bm{y}}_{t}dt.$ (59) $$

Let us discuss the link between this solution and OT.
The formula of OT map between two centered Gaussian distributions $\mathcal{N}({\bm{0}},{\bm{\Sigma}}_{1})$ and $\mathcal{N}({\bm{0}},{\bm{\Sigma}}_{2})$ is well known. In [], the authors give the linear map (affine when the distributions are not centered) $\bm{T}:\bm{X}\mapsto\bm{A}\bm{X}$ with

$$ $\bm{A}={\bm{\Sigma}}_{1}^{-\frac{1}{2}}\left({\bm{\Sigma}}_{1}^{\frac{1}{2}}{ \bm{\Sigma}}_{2}{\bm{\Sigma}}_{1}^{\frac{1}{2}}\right)^{\frac{1}{2}}{\bm{ \Sigma}}_{1}^{-\frac{1}{2}}.$ (60) $$

When ${\bm{\Sigma}}_{1}$ and ${\bm{\Sigma}}_{2}$ commute, this expression simplifies to:

$$ $\bm{A}={\bm{\Sigma}}_{1}^{-1/2}{\bm{\Sigma}}_{2}^{1/2}.$ (61) $$

We showed that the solution (Equation 54) of the backward probability flow in the finite interval $[0,t]$, with $0\leq t\leq T$, corresponds to applying to the initial point ${\bm{y}}_{0}$ the linear map

$$ $\bm{A}={\bm{\Sigma}}_{T}^{-\frac{1}{2}}{\bm{\Sigma}}_{T-t}^{\frac{1}{2}},$ (62) $$

that is, the OT map between $p_{T}=\mathcal{N}({\bm{0}},{\bm{\Sigma}}_{T})$ and
$p_{T-t}=\mathcal{N}({\bm{0}},{\bm{\Sigma}}_{T-t})$.

Let us now derive the covariance matrix of the solution, which characterises a Gaussian distribution.

$$ $\mathrm{Cov}({\bm{y}}_{t})={\bm{\Sigma}}^{-1/2}_{T}{\bm{\Sigma}}_{T-t}^{1/2} \mathrm{Cov}({\bm{y}}_{0}){\bm{\Sigma}}^{-1/2}_{T-t}{\bm{\Sigma}}_{T}^{1/2}.$ (63) $$

In particular, if $\mathrm{Cov}({\bm{y}}_{0})$ and ${\bm{\Sigma}}$ commute,

$$ $\mathrm{Cov}({\bm{y}}_{t})={\bm{\Sigma}}^{-1}_{T}{\bm{\Sigma}}_{T-t}\mathrm{ Cov}({\bm{y}}_{0}).$ (64) $$

### B.4 Proof of Proposition 4

For $0\leq t\leq T$, denoting $\left(\lambda_{i,t}\right)_{1\leq i\leq d}$ the eigenvalues of ${\bm{\Sigma}}_{t}$, the eigenvalues of ${\tilde{{\bm{\Sigma}}}}_{t}=\mathrm{Cov}({\tilde{\bm{y}}}_{T-t})$ are

$$ $\tilde{\lambda}_{i,t}=\lambda_{i,t}+e^{-2(B_{T}-B_{t})}\lambda_{i,t}^{2}\frac{ 1}{\lambda_{i,T}}\left(\frac{1}{\lambda_{i,T}}-1\right),\quad i=1,\ldots,d.$ (65) $$

and the eigenvalues of ${\widehat{{\bm{\Sigma}}}}_{t}=\mathrm{Cov}({\widehat{\bm{y}}}_{T-t})$ are

$$ $\widehat{\lambda}_{i,t}=\frac{\lambda_{i,t}}{\lambda_{i,T}},\quad 1\leq i\leq d.$ (66) $$

Consequently, $\operatorname{\text{$\textbf{W}_{2}$}}2(p_{t},{\tilde{p}}_{t})$ is the sum of the squares of all:

$$ $\sqrt{\lambda_{i,t}}-\sqrt{\tilde{\lambda}_{i,t}}=\sqrt{\lambda_{i,t}}\left(1- \sqrt{1+e^{-2(B_{T}-B_{t})}\lambda_{i,t}\frac{1}{\lambda_{i,T}}\left(\frac{1}{ \lambda_{i,T}}-1\right)}\right).$ (67) $$

Similarly, $\operatorname{\text{$\textbf{W}_{2}$}}2(p_{t},{\widehat{p}}_{t})$ is the sum of the squares of all:

$$ $\displaystyle\sqrt{\lambda_{i,t}}-\sqrt{\widehat{\lambda}_{i,t}}=\sqrt{\lambda _{i,t}}\left(1-\sqrt{\frac{1}{\lambda_{i,T}}}\right)$ (68) $\displaystyle=\sqrt{\lambda_{i,t}}\left(1-\sqrt{1+\left(\frac{1}{\lambda_{i,T} }-1\right)}\right).$ (69) $$

Let us now compare individually these differences.

$$ $\displaystyle\frac{e^{-2(B_{T}-B_{t})}\lambda_{i,t}\frac{1}{\lambda_{i,T}} \left(\frac{1}{\lambda_{i,T}}-1\right)}{\frac{1}{\lambda_{i,T}}-1}$ $\displaystyle=e^{-2(B_{T}-B_{t})}\frac{\lambda_{i,t}}{\lambda_{i,T}}$ (70) $\displaystyle=e^{-2(B_{T}-B_{t})}\frac{e^{-2B_{t}}(\lambda_{i}-1)+1}{e^{-2B_{T }}(\lambda_{i}-1)+1}$ (71) $\displaystyle=\frac{(\lambda_{i}-1)+e^{2B_{t}}}{(\lambda_{i}-1)+e^{2B_{T}}}$ (72) $\displaystyle<1.$ (73) $$

Case 1: $0<\lambda_{i}<1$ and $t>0$

In this case, $\lambda_{i,T}<1$ and:

$$ $0<e^{-2(B_{T}-B_{t})}\lambda_{i,t}\frac{1}{\lambda_{i,T}}\left(\frac{1}{ \lambda_{i,T}}-1\right)<\frac{1}{\lambda_{i,T}}-1.$ (74) $$

Thus,

$$ $\displaystyle\left|\sqrt{\lambda_{i,t}}-\sqrt{\tilde{\lambda}_{i,t}}\right|$ $\displaystyle=\sqrt{\tilde{\lambda}_{i,t}}-\sqrt{\lambda_{i,t}}$ (75) $\displaystyle=\sqrt{\lambda_{i,t}}\left(\sqrt{1+e^{-2(B_{T}-B_{t})}\lambda_{i} ^{t}\frac{1}{\lambda_{i,T}}\left(\frac{1}{\lambda_{i,T}}-1\right)}-1\right)$ (76) $\displaystyle<\sqrt{\lambda_{i,t}}\left(\sqrt{1+\left(\frac{1}{\lambda_{i,T}}- 1\right)}-1\right)$ (77) $\displaystyle=\sqrt{\widehat{\lambda}_{i,t}}-\sqrt{\lambda_{i,t}}$ (78) $\displaystyle=\left|\sqrt{\lambda_{i,t}}-\sqrt{\widehat{\lambda}_{i,t}}\right|.$ (79) $$

Case 2: $\lambda_{i}=0$ and $t=0$.

In this case, for $1\leq i\leq d$, $\widehat{\lambda}_{i,T}=\tilde{\lambda}_{i,T}=0$.

Case 3: $\lambda_{i}=1$.

In this case, for $1\leq i\leq d$, $\widehat{\lambda}_{i,t}=\tilde{\lambda}_{i,t}=1$.

Case 4: $1<\lambda_{i}$.

In this case, $\lambda_{i,T}\geq 1$, and
$\frac{e^{-2(B_{T}-B_{t})}\lambda_{i,t}\frac{1}{\lambda_{i,T}}\left(\frac{1}{
\lambda_{i,T}}-1\right)}{\frac{1}{\lambda_{i,T}}-1}=e^{-2(B_{T}-B_{t})}\frac{
\lambda_{i,t}}{\lambda_{i,T}}<1$ provides

$$ $e^{-2(B_{T}-B_{t})}\lambda_{i,t}\frac{1}{\lambda_{i,T}}\left(\frac{1}{\lambda_ {i,T}}-1\right)>\frac{1}{\lambda_{i,T}}-1.$ (80) $$

Finally,

$$ $\displaystyle\left|\sqrt{\lambda_{i,t}}-\sqrt{\tilde{\lambda}_{i,t}}\right|$ $\displaystyle=\sqrt{\lambda_{i,t}}-\sqrt{\tilde{\lambda}_{i,t}}$ (81) $\displaystyle=\sqrt{\lambda_{i,t}}\left(1-\sqrt{1+e^{-2(B_{T}-B_{t})}\lambda_{ i,T}\frac{1}{\lambda_{i,T}}\left(\frac{1}{\lambda_{i,T}}-1\right)}\right)$ (82) $\displaystyle<\sqrt{\lambda_{i,t}}\left(1-\sqrt{1+\left(\frac{1}{\lambda_{i,T} }-1\right)}\right)$ (83) $\displaystyle=\sqrt{\lambda_{i,t}}-\sqrt{\widehat{\lambda}_{i,t}}$ (84) $\displaystyle=\left|\sqrt{\lambda_{i,t}}-\sqrt{\widehat{\lambda}_{i,t}}\right|.$ (85) $$

This case study provides:

$$ $\operatorname{\text{$\textbf{W}_{2}$}}2({\tilde{p}}_{t},p_{t})\leq \operatorname{\text{$\textbf{W}_{2}$}}2({\widehat{p}}_{t},p_{t}).$ (86) $$

## Appendix C Gaussian CIFAR-10 samples

The Gaussian CIFAR-10 produces unstructured images. A grid of samples is presented in Figure 4. To sample from this Gaussian, the empirical covariance matrix of size $\mathbb{R}^{(3\times 32\times 32)\times(3\times 32\times 32)}$is computed and then the SVD decomposition to extract a square root matrix (see source code).

Figure: Figure 4: CIFAR-10 Gaussian samples. Samples are generated from the Gaussian distribution fitting the CIFAR-10 dataset.
Refer to caption: extracted/5986824/Figures/Gaussian_cifar10/cifar10_samples.png

## Appendix D Computation of the 2-Wasserstein distances for numerical schemes

The 2-Wasserstein errors can be computed by using Equation 19 recalled here:

$$ $\operatorname{\text{$\textbf{W}_{2}$}}2(\mathcal{N}({\bm{0}},{\bm{\Sigma}}_{1} ),\mathcal{N}({\bm{0}},{\bm{\Sigma}}_{2}))^{2}=\sum_{1\leq i\leq d}(\sqrt{ \lambda_{i,1}}-\sqrt{\lambda_{i,2}})^{2}.$ ( 19 ) $$

for two centered Gaussians $\mathcal{N}({\bm{0}},{\bm{\Sigma}}_{1})$ and $\mathcal{N}({\bm{0}},{\bm{\Sigma}}_{2})$ such that ${\bm{\Sigma}}_{1},{\bm{\Sigma}}_{2}$ are simultaneously diagonalizable with respective eigenvalues $\left(\lambda_{i,1}\right)_{1\leq i\leq d},\left(\lambda_{i,2}\right)_{1\leq i
\leq d}$. We aim at computing these errors between the Gaussian process folowing $\left(p_{t}\right)_{0\leq t\leq T}$ and the processes induced by the discretization schemes. Table 1 shows that each discretization scheme leads to a discrete time Gaussian process whose covariance matrix is diagonalizable in the basis of ${\bm{\Sigma}}$ when intialize them with either $\mathcal{N}_{0}$ (usual sampling) or $p_{T}$ (no initialization error). Let detail this point for the Euler-Maruyama (EM) scheme. Let denote $\left(\bm{v}_{i}\right)_{1\leq i\leq d}$ a basis of eigenvectors of ${\bm{\Sigma}}$ and its associated eigenvalues $\left(\lambda_{i}\right)_{1\leq i\leq d}$. For $0\leq t\leq T$, ${\bm{\Sigma}}_{t}=e^{-2B_{t}}{\bm{\Sigma}}+(1-e^{-2B_{t}})\bm{I}$ and for all $1\leq i\leq d$,

$$ ${\bm{\Sigma}}_{t}\bm{v}_{i}=\left(e^{-2B_{t}}\lambda_{i}+(1-e^{-2B_{t}})\right )\bm{v}_{i}$ (87) $$

Consequently ${\bm{\Sigma}}_{t}$ admits the same eigenvectors than ${\bm{\Sigma}}$ with associated eigenvalues $\left(\lambda_{i,t}\right)_{1\leq i\leq d}=\left(e^{-2B_{t}}\lambda_{i}+(1-e^{
-2B_{t}})\right)_{1\leq i\leq d}$. Let study the covariance matrix of the EM process.
Let denote $\left({\bm{\Sigma}}_{k}^{\text{$\Delta$,EM}}\right)_{\begin{subarray}{c}1\leq i
\leq d,0\leq k\leq N\end{subarray}}$ the covariance matrix of the Gaussian process generated by the EM scheme at each step and $\left(\lambda_{i,k}^{\text{$\Delta$,EM}}\right)_{\begin{subarray}{c}1\leq i
\leq d,0\leq k\leq N\end{subarray}}$ its eigenvalues. First,

$$ ${\bm{\Sigma}}_{0}^{\text{$\Delta$,EM}}=\left\{\begin{array}[]{ll}\bm{I}&\mbox{ if ${\bm{y}}_{T}$ is initialized at $\mathcal{N}_{0}$ }\\ {\bm{\Sigma}}_{T}&\mbox{if ${\bm{y}}_{T}$ is initialized at $p_{T}$ }\\ \end{array}\right..$ (88) $$

And consequently,

$$ $\lambda_{i,0}^{\text{$\Delta$,EM}}=\left\{\begin{array}[]{ll}1&\mbox{if ${\bm{ y}}_{T}$ is initialized at $\mathcal{N}_{0}$ }\\ e^{-2B_{T}}\lambda_{i}+(1-e^{-2B_{T}})&\mbox{if ${\bm{y}}_{T}$ is initialized at $p_{T}$ }\\ \end{array}\right.\quad 1\leq i\leq d.$ (89) $$

Then, by Table 1,

$$ $\displaystyle{\tilde{\bm{y}}^{\Delta,\text{EM}}}_{1}$ $\displaystyle=\left(\bm{I}+\Delta_{t}\beta_{T-t_{0}}\left(\bm{I}-2{\bm{\Sigma} }_{T-t_{0}}^{-1}\right)\right){\tilde{\bm{y}}^{\Delta,\text{EM}}}_{0}+\sqrt{2 \Delta_{t}\beta_{T-t_{0}}}{\bm{z}}_{0},\leavevmode\nobreak\ {\bm{z}}_{0}\sim \mathcal{N}_{0}$ (90) $$

and

$$ $\displaystyle{\bm{\Sigma}}_{1}^{\text{$\Delta$,EM}}$ $\displaystyle=\left(\bm{I}+\Delta_{t}\beta_{T-t_{0}}\left(\bm{I}-2{\bm{\Sigma} }_{T-t_{0}}^{-1}\right)\right){\bm{\Sigma}}_{0}^{\text{$\Delta$,EM}}\left(\bm{ I}+\Delta_{t}\beta_{T-t_{0}}\left(\bm{I}-2{\bm{\Sigma}}_{T-t_{0}}^{-1}\right) \right)^{T}+2\Delta_{t}\beta_{T-t_{0}}\bm{I}$ (91) $\displaystyle=\left(\bm{I}+\Delta_{t}\beta_{T-t_{0}}\left(\bm{I}-2{\bm{\Sigma} }_{T-t_{0}}^{-1}\right)\right)^{2}{\bm{\Sigma}}_{0}^{\text{$\Delta$,EM}}+2 \Delta_{t}\beta_{T-t_{0}}\bm{I}\text{ because ${\bm{\Sigma}}$ and ${\bm{\Sigma }}_{0}^{\text{$\Delta$,EM}}$ commute}.$ (92) $$

Let $1\leq i\leq d$,

$$ $\displaystyle{\bm{\Sigma}}_{1}^{\text{$\Delta$,EM}}\bm{v}_{i}$ $\displaystyle=\left[\left(1+\Delta_{t}\beta_{T-t_{0}}\left(\bm{I}-\frac{2}{ \lambda_{i,T-t_{0}}}\right)\right)^{2}\lambda_{i,0}^{\text{$\Delta$,EM}}+2 \Delta_{t}\beta_{T-t_{0}}\right]\bm{v}_{i}$ (93) $$

Consequently, $(\bm{v}_{i})_{1\leq i\leq d}$ is also a basis of eigenvectors of ${\bm{\Sigma}}_{1}^{\text{$\Delta$,EM}}$ and

$$ $\lambda_{i,1}^{\text{$\Delta$,EM}}=\left(1+\Delta_{t}\beta_{T-t_{0}}\left(\bm{ I}-\frac{2}{\lambda_{i,T-t_{0}}}\right)\right)^{2}\lambda_{i,0}^{\text{$\Delta $,EM}}+2\Delta_{t}\beta_{T-t_{0}},\quad 1\leq i\leq d.$ (94) $$

Thus, we can obtain the eigenvalues $\left(\lambda_{i,k}^{\text{$\Delta$,EM}}\right)_{\begin{subarray}{c}1\leq i
\leq d,0\leq k\leq N\end{subarray}}$ at each time and plot at each time

$$ $\sqrt{\sum_{1\leq i\leq d}\left(\sqrt{\lambda_{i,T-t_{k}}}-\sqrt{\lambda_{i,k} ^{\text{$\Delta$,EM}}}\right)},\quad 1\leq k\leq N$ (95) $$

as done in Figure 1. These computations can be led for the different schemes, as presented in Table 4.

**Table 4: Recursive form of the eigenvalues of the covariance matrix associated with the Gaussian process generated by the different schemes for a regular time schedule $\left(t_{k}\right)_{0\leq k\leq N}$ with steps $\Delta_{t}=\tfrac{T}{N}$.**
| SDE schemes | EM | $\lambda_{i,k+1}^{\Delta,\text{EM}}=\left(1+\Delta_{t}\beta_{T-t_{k}}\left(1-<br>\tfrac{2}{\lambda_{i,T-t_{k}}}\right)\right)^{2}\lambda_{i,k}^{\Delta,\text{EM<br>}}+2\Delta_{t}\beta_{T-t_{k}},\quad 1\leq i\leq d,0\leq k\leq N-1$ |
| --- | --- | --- |
| EI | $\lambda_{i,k+1}^{\Delta,\text{EI}}=\left(1+\gamma_{1,k}\left(1-\tfrac{2}{<br>\lambda_{i,T-t_{k}}}\right)\right)^{2}\lambda_{i,k}^{\Delta,\text{EI}}+2\gamma<br>_{2,k}\quad 1\leq i\leq d,0\leq k\leq N-1$ |  |
| ODE schemes | Euler | $\lambda_{i,k+1}^{\Delta,\text{Euler}}=\left(1+\Delta_{t}\beta_{T-t_{k}}\left(1<br>-\tfrac{1}{\lambda_{i,T-t_{k}}}\right)\right)^{2}\lambda_{i}^{\text{Euler},k}<br>\quad 1\leq i\leq d,0\leq k\leq N-1$ |
| Heun | $\lambda_{i,k+1}^{\Delta,\text{Heun}}=\left(1+\frac{\Delta_{t}}{2}\beta_{T-t_{k<br>}}\left(1-\tfrac{1}{\lambda_{i,T-t_{k}}}\right)+\frac{\Delta_{t}}{2}\beta_{T-t<br>_{k+1}}\left(1-\tfrac{1}{\lambda_{i,T-t_{k+1}}}\right)\left(1+\Delta_{t}\beta_<br>{T-t_{k}}\left(1-\tfrac{1}{\lambda_{i,T-t_{k}}}\right)\right)\right)\lambda_{i<br>,k}^{\Delta,\text{Heun}}$ |  |
|  |  | $\quad 1\leq i\leq d,0\leq k\leq N-1$ |

## Appendix E Theoretical Wasserstein distance for the ADSN model

As done for the Gaussian CIFAR-10, the Wasserstein errors can be computed for the ADSN model as shown in Figure 5 and Table 5.

Figure: (a) Initialization error along the integration time.
Refer to caption: x5.png

**Table 5: Ablation study of Wasserstein errors for the Gaussian microtextures. For a given discretization scheme, the table presents the Wasserstein distance associated with the truncation error for different values of $\varepsilon$. The columns $p_{T}$ and $\mathcal{N}_{0}$ show the influence of the initialization error. The continuous column corresponds to the continuous SDE or ODE linked with the scheme (identical values for EM, EI and Euler, Heun). Note that the Heun scheme is not defined without truncation time due to the zero eigenvalue.**
|  |  | Continuous | $N=50$ | $N=250$ | $N=500$ | $N=1000$ |  |  |  |  |  |
| --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- |
|  |  | $p_{T}$ | $\mathcal{N}_{0}$ | $p_{T}$ | $\mathcal{N}_{0}$ | $p_{T}$ | $\mathcal{N}_{0}$ | $p_{T}$ | $\mathcal{N}_{0}$ | $p_{T}$ | $\mathcal{N}_{0}$ |
| EM | $\varepsilon=0$ | 0 | 5.2E-06 | 53.37 | 53.37 | 10.58 | 10.58 | 6.27 | 6.27 | 4.02 | 4.02 |
| $\varepsilon=10^{-5}$ | 0.36 | 0.36 | 53.35 | 53.35 | 10.57 | 10.57 | 6.26 | 6.26 | 4.02 | 4.02 |  |
| $\varepsilon=10^{-4}$ | 1.15 | 1.15 | 53.21 | 53.21 | 10.53 | 10.53 | 6.25 | 6.25 | 4.03 | 4.03 |  |
| $\varepsilon=10^{-3}$ | 3.84 | 3.84 | 51.92 | 51.92 | 10.55 | 10.55 | 6.80 | 6.80 | 5.16 | 5.16 |  |
| EI | $\varepsilon=0$ | 0 | 5.2E-06 | 30.91 | 30.91 | 8.85 | 8.85 | 5.71 | 5.71 | 3.84 | 3.84 |
| $\varepsilon=10^{-5}$ | 0.36 | 0.36 | 30.92 | 30.92 | 8.85 | 8.85 | 5.72 | 5.72 | 3.84 | 3.84 |  |
| $\varepsilon=10^{-4}$ | 1.15 | 1.15 | 31.01 | 31.01 | 8.92 | 8.92 | 5.78 | 5.78 | 3.90 | 3.90 |  |
| $\varepsilon=10^{-3}$ | 3.84 | 3.84 | 31.94 | 31.94 | 9.74 | 9.74 | 6.76 | 6.76 | 5.24 | 5.24 |  |
| Euler | $\varepsilon=0$ | 0 | 6.4E-03 | 5.69 | 5.70 | 3.27 | 3.27 | 2.50 | 2.51 | 1.87 | 1.87 |
| $\varepsilon=10^{-5}$ | 0.36 | 0.36 | 5.70 | 5.71 | 3.28 | 3.28 | 2.53 | 2.53 | 1.90 | 1.90 |  |
| $\varepsilon=10^{-4}$ | 1.15 | 1.15 | 5.80 | 5.80 | 3.43 | 3.43 | 2.72 | 2.72 | 2.15 | 2.15 |  |
| $\varepsilon=10^{-3}$ | 3.84 | 3.84 | 6.79 | 6.79 | 4.85 | 4.85 | 4.41 | 4.41 | 4.14 | 4.14 |  |
| Heun | $\varepsilon=0$ | 0 | 6.4E-03 | - | - | - | - | - | - | - | - |
| $\varepsilon=10^{-5}$ | 0.36 | 0.36 | 2.4E+03 | 2.4E+03 | 3.0E+02 | 3.0E+02 | 1.1E+02 | 1.1E+02 | 40.00 | 40.00 |  |
| $\varepsilon=10^{-4}$ | 1.15 | 1.15 | 2.3E+02 | 2.3E+02 | 26.34 | 26.34 | 8.54 | 8.54 | 2.01 | 2.01 |  |
| $\varepsilon=10^{-3}$ | 3.84 | 3.84 | 15.42 | 15.42 | 2.25 | 2.25 | 3.40 | 3.40 | 3.73 | 3.73 |  |

## Appendix F Study of the covariance matrix of the ADSN distribution

### F.1 Reminders on the Discrete Fourier Transform (DFT)

For a given image $\bm{v}\in\mathbb{R}^{3\times{M\times N}}$, we define the DFT of $\bm{v}$, $\widehat{\bm{v}}\in\mathbb{R}^{3\times{M\times N}}$ such that for $1\leq c\leq 3$,$\xi\in\mathbb{R}^{{M\times N}}$

$$ $\widehat{\bm{v}}_{c,\xi}=\sum_{x\in{M\times N}}\bm{v}_{c,x}\exp({-\frac{2i\pi x _{1}\xi_{1}}{M}})\exp(-\frac{2i\pi x_{2}\xi_{2}}{N}),\quad i^{2}=-1$ (96) $$

where $\widehat{\bm{v}}_{c,\xi}$ is the value of $\widehat{\bm{v}}$ at coordinate $\xi$ of the k-th channel of $\widehat{\bm{v}}$. For $\bm{u}\in\mathbb{R}^{3{M\times N}}$, by defining $\bm{u}\star\bm{v}$ the periodic convolution such that for $1\leq c\leq 3,x\in\mathbb{R}^{{M\times N}}$:

$$ $\left(\bm{u}\star\bm{v}\right)_{c,x}=\sum_{y\in{M\times N}}\bm{u}_{c,x-y}\bm{v }_{c,y}$ (97) $$

we have:

$$ $\widehat{\bm{u}\star\bm{v}}=\widehat{\bm{u}}\odot\widehat{\bm{v}},$ (98) $$

where $\odot$ is the componentwise product.

### F.2 Eigenvectors of the covariance matrix of the ADSN distribution

Let $\bm{u}\in\mathbb{R}^{3\times{M\times N}}$ and its associated texton $\mathbf{t}\in\mathbb{R}^{3\times{M\times N}}$. The distribution $\operatorname{ADSN}(\bm{u})$ is the Gaussian distribution of $\bm{X}=\mathbf{t}\star\bm{w}$ such that:

$$ $\bm{X}_{i}=\mathbf{t}_{i}\star\bm{w}\in\mathbb{R}^{{M\times N}},1\leq i\leq 3, \bm{w}\sim\mathcal{N}_{0}$ (99) $$

Consequently, denoting ${\bm{\Sigma}}$ the covariance of $\operatorname{ADSN}(\bm{u})$, for $\bm{v}\in\mathbb{R}^{3{M\times N}}$,

$$ $\widehat{{\bm{\Sigma}}\bm{v}}_{i}=\widehat{\mathbf{t}}_{i}\overline{\widehat{ \mathbf{t}}_{1}}\widehat{\bm{v}}_{1}+\widehat{\mathbf{t}}_{i}\overline{ \widehat{\mathbf{t}}_{2}}\widehat{\bm{v}}_{2}+\widehat{\mathbf{t}}_{i} \overline{\widehat{\mathbf{t}}_{3}}\widehat{\bm{v}}_{3}=\widehat{\mathbf{t}}_{ i}\left(\overline{\widehat{\mathbf{t}}_{1}}\widehat{\bm{v}}_{1}+\overline{ \widehat{\mathbf{t}}_{2}}\widehat{\bm{v}}_{2}+\overline{\widehat{\mathbf{t}}_{ 3}}\widehat{\bm{v}}_{3}\right)$ (100) $$

This equation proves that the kernel of ${\bm{\Sigma}}$ contains the kernel of $\bm{v}\in\mathbb{R}^{3\times{M\times N}}\mapsto\overline{\widehat{\mathbf{t}}_
{1}}\widehat{\bm{v}}_{1}+\overline{\widehat{\mathbf{t}}_{2}}\widehat{\bm{v}}_{
2}+\overline{\widehat{\mathbf{t}}_{3}}\widehat{\bm{v}}_{3}\in\mathbb{R}^{{M
\times N}}$ which has a dimension greater than $2MN$. Consequently, 0 0 is eigenvalue of ${\bm{\Sigma}}$ with multiplicity greater than $2MN$. Furthermore, for $\xi\in\mathbb{R}^{{M\times N}}$, denoting $\bm{u}^{1,\xi}$ such that:

$$ $\widehat{\bm{u}}_{i}^{1,\xi}(\omega)=\mathbf{1}_{\omega=\xi}\widehat{\mathbf{t }}_{i}(\omega),1\leq i\leq 3,\omega\in\mathbb{R}^{{M\times N}}$ (101) $$

we have,

$$ ${\bm{\Sigma}}\bm{u}^{1,\xi}=(|\widehat{\mathbf{t}}_{1}(\xi)|^{2}+|\widehat{ \mathbf{t}}_{2}(\xi)|^{2}+|\widehat{\mathbf{t}}_{3}(\xi)|^{2})\bm{u}^{1,\xi}.$ (102) $$

Furthermore, the family $\left(\bm{u}^{1,\xi}\right)_{\xi\in{M\times N}}$ is orthogonal. Thus, the eigenvalues of ${\bm{\Sigma}}$ are $\left(|\widehat{\mathbf{t}}_{1}(\xi)|^{2}+|\widehat{\mathbf{t}}_{2}(\xi)|^{2}+
|\widehat{\mathbf{t}}_{3}(\xi)|^{2}\right)_{\xi\in{M\times N}}$ and 0 0 with multiplicity $2MN$.

For $\xi\in\mathbb{R}^{M\times N}$, we denote $\bm{u}^{2,\xi}$,$\bm{u}^{3,\xi}$ such that for $\omega\in\mathbb{R}^{{M\times N}}$:

$$ $\displaystyle\left\{\begin{array}[]{ll}\widehat{\bm{u}}_{1}^{2,\xi}(\omega)&=- \mathbf{1}_{\omega=\xi}\overline{\widehat{\mathbf{t}}}_{3}(\omega)\\ \widehat{\bm{u}}_{2}^{2,\xi}(\omega)&=0\\ \widehat{\bm{u}}_{3}^{2,\xi}(\omega)&=\mathbf{1}_{\omega=\xi}\overline{ \widehat{\mathbf{t}}}_{1}(\omega)\\ \end{array}\right.$ (106) $\displaystyle\left\{\begin{array}[]{ll}\widehat{\bm{u}}_{1}^{3,\xi}(\omega)&=0 \\ \widehat{\bm{u}}_{2}^{3,\xi}(\omega)&=-\mathbf{1}_{\omega=\xi}\overline{ \widehat{\mathbf{t}}}_{3}(\omega)\\ \widehat{\bm{u}}_{3}^{3,\xi}(\omega)&=\mathbf{1}_{\omega=\xi}\overline{ \widehat{\mathbf{t}}}_{2}(\omega)\\ \end{array}\right.$ (110) $$

We have

$$ $\displaystyle{\bm{\Sigma}}\bm{u}^{2,\xi}$ $\displaystyle=0.\bm{u}^{2,\xi}$ (111) $\displaystyle{\bm{\Sigma}}\bm{u}^{3,\xi}$ $\displaystyle=0.\bm{u}^{3,\xi}.$ (112) $$

Then, applying the orthonomalization of Gram-Schmidt on each tuple $(\bm{u}^{1,\xi},\bm{u}^{2,\xi},\bm{u}^{3,\xi})_{\xi\in\mathbb{R}^{M\times N}}$, we obtain an orthonormal basis in the Fourier domain $(\bm{v}^{1,\xi},\bm{v}^{2,\xi},\bm{v}^{3,\xi})_{\xi\in\mathbb{R}^{M\times N}}$ of eigenvectors of ${\bm{\Sigma}}$. More precisely, for $\xi_{1},\xi_{2}\in\mathbb{R}^{{M\times N}}$, $1\leq j_{1},j_{2}\leq 3$,

$$ $\displaystyle\left(\overline{\widehat{\bm{v}}}^{j_{1},\xi_{1}}\right)^{T} \widehat{\bm{v}}^{j_{2},\xi_{2}}$ $\displaystyle=\sum_{\begin{subarray}{c}x_{1}\in{M\times N}\\ x_{2}\in{M\times N}\end{subarray}}\overline{\widehat{\bm{v}}}^{j_{1},\xi_{1}}_ {x_{1}}\widehat{\bm{v}}^{j_{2},\xi_{2}}_{x_{2}}$ (113) $\displaystyle=\mathbf{1}_{\begin{subarray}{c}j_{1}=j_{2}\\ \xi_{1}=\xi_{2}\end{subarray}}$ (114) $$

which is applying the square root of ${\bm{\Sigma}}$ to the white Gaussian noise $\bm{w}$. Furthermore, we can ensure that for $\xi\neq\omega\in\mathbb{R}^{{M\times N}},1\leq j\leq 3$, $\widehat{\bm{v}}^{j,\xi}(\omega)=0$ such that only the frequency $\xi$ is active in the Fourier transform of ${\bm{v}}^{j,\xi}$. Consequently, for $\bm{w}\in\mathbb{R}^{3{M\times N}}$,

$$ $\overline{\widehat{\bm{w}}}^{T}{\bm{v}}^{j,\xi}=\sum_{1\leq i\leq 3}\overline{ \widehat{\bm{w}}}_{i}(\xi)\widehat{\bm{v}}^{j,\xi}_{i}(\xi).$ (115) $$

In particular,

$$ $\left(\widehat{\overline{\bm{v}}}^{j,\xi}\right)^{T}=\|\widehat{\bm{v}}^{j,\xi }\|^{2}=\sum_{1\leq i\leq 3}\left|{\bm{v}}^{j,\xi}_{i}(\xi)\right|^{2}=1.$ (116) $$

### F.3 Computation of the empirical Wasserstein error in the ADSN covariance diagonalization basis

Let consider a Gaussian distribution $\mathcal{N}({\bm{0}},\bm{\Gamma})$ such that there exists $(\lambda_{1}^{\xi},\lambda_{2}^{\xi},\lambda_{3}^{\xi})_{\xi\in\mathbb{R}^{M
\times N}}$ such that for all $\xi\in\mathbb{R}^{{M\times N}}$,

$$ $\bm{\Gamma}\bm{v}^{j,\xi}=\lambda_{j}^{\xi}\bm{v}^{j,\xi},\quad 1\leq j\leq 3.$ (117) $$

Let $\bm{w}\sim\mathcal{N}_{0}\in\mathbb{R}^{3{M\times N}}$, $(\bm{v}^{1,\xi},\bm{v}^{2,\xi},\bm{v}^{3,\xi})_{\xi\in\mathbb{R}^{M\times N}}$ is an orthonormal basis in the Fourier domain such that:

$$ $\displaystyle\widehat{\bm{w}}$ $\displaystyle=\sum_{\xi\in\mathbb{R}^{{M\times N}}}\left(\left[\overline{ \widehat{\bm{w}}}^{T}\widehat{\bm{v}}^{1,\xi}\right]\widehat{\bm{v}}^{1,\xi}+ \left[\overline{\widehat{\bm{w}}}^{T}\widehat{\bm{v}}^{2,\xi}\right]\widehat{ \bm{v}}^{2,\xi}+\left[\overline{\widehat{\bm{w}}}^{T}\widehat{\bm{v}}^{3,\xi} \right]\widehat{\bm{v}}^{3,\xi}\right)$ (118) $$

A sample drawn from $\mathcal{N}({\bm{0}},\bm{\Gamma})$ has the same distribution as $\bm{Y}$ given by

$$ $\widehat{\bm{Y}}=\sum_{\xi\in\mathbb{R}^{{M\times N}}}\sqrt{\lambda_{1}^{\xi}} \left[\overline{\widehat{\bm{w}}}^{T}\widehat{\bm{v}}^{1,\xi}\right]\widehat{ \bm{v}}^{1,\xi}+\sum_{\xi\in\mathbb{R}^{{M\times N}}}\sqrt{\lambda_{2}^{\xi}} \left[\overline{\widehat{\bm{w}}}^{T}\widehat{\bm{v}}^{2,\xi}\right]\widehat{ \bm{v}}^{2,\xi}+\sum_{\xi\in\mathbb{R}^{{M\times N}}}\sqrt{\lambda_{3}^{\xi}} \left[\overline{\widehat{\bm{w}}}^{T}\widehat{\bm{v}}^{3,\xi}\right]\widehat{ \bm{v}}^{3,\xi}.$ (120) $$

Note that the three channels of $\bm{w}$ are independent. Furthermore, for $1\leq j\leq 3$

$$ $\displaystyle\left(\overline{\widehat{\bm{v}}}^{j,\xi}\right)^{T}\widehat{\bm{ Y}}$ $\displaystyle=\sqrt{\lambda_{1}^{\xi}}\left[\overline{\widehat{\bm{w}}}^{T} \widehat{\bm{v}}^{j,\xi}\right]\left\|\widehat{\bm{v}}^{j,\xi}\right\|^{2}= \sqrt{\lambda_{1}^{\xi}}\left[\overline{\widehat{\bm{w}}}^{T}\widehat{\bm{v}}^ {j,\xi}\right]$ (121) $\displaystyle\left|\left(\overline{\widehat{\bm{v}}}^{j,\xi}\right)^{T} \widehat{\bm{Y}}\right|^{2}$ $\displaystyle=\lambda_{j}^{\xi}\left|\overline{\widehat{\bm{w}}}^{T}\widehat{ \bm{v}}^{j,\xi}\right|^{2}$ (122) $\displaystyle\mathbb{E}\left[\left|\left(\overline{\widehat{\bm{v}}}^{j,\xi} \right)^{T}\widehat{\bm{Y}}\right|^{2}\right]$ $\displaystyle=\lambda_{j}^{\xi}\mathbb{E}\left[\left|\overline{\widehat{\bm{w} }}^{T}\widehat{\bm{v}}^{j,\xi}\right|^{2}\right]$ (123) $\displaystyle\mathbb{E}\left[\left|\overline{\widehat{\bm{w}}}^{T}\widehat{\bm {v}}^{j,\xi}\right|^{2}\right]$ $\displaystyle=\sum_{1\leq c_{1},c_{2}\leq 3}\mathbb{E}\left[\overline{\widehat {\bm{w}}}_{c_{1}}(\xi)\widehat{\bm{w}}_{c_{2}}(\xi)\right]\widehat{\bm{v}}_{c_ {1}}^{j,\xi}(\xi)\overline{\widehat{\bm{v}}}_{c_{2}}(\xi)\text{ by \lx@cref{ creftypecap~refnum}{eq:v_inner} }$ (124) $\displaystyle=\sum_{1\leq c\leq 3}\mathbb{E}\left[\left|\widehat{\bm{w}}_{c}( \xi)\right|^{2}\right]\left|\widehat{\bm{v}}_{c}^{j,\xi}(\xi)\right|^{2}\text{ because the channels are inependent}$ (125) $\displaystyle=3MN\sum_{1\leq c\leq 3}\left|\widehat{\bm{v}}_{c}^{j,\xi}(\xi) \right|^{2}\text{ because $\mathbb{E}\left[\left|\widehat{\bm{w}}_{c}(\xi) \right|^{2}\right]=MN$}$ (126) $\displaystyle=3MN\text{ by \lx@cref{creftypecap~refnum}{eq:v_norm}. }$ (127) $$

Finally,

$$ $\displaystyle\mathbb{E}\left[\left|\left(\overline{\widehat{\bm{v}}}^{j,\xi} \right)^{T}\widehat{\bm{Y}}\right|^{2}\right]=3MN\lambda_{1}^{\xi}$ (128) $$

Finally, for a given sampling $\left(\bm{Y}_{k}\right)_{1\leq k\leq N_{\text{samples}}}$ following the distribution $\mathcal{N}({\bm{0}},\bm{\Gamma})$, an estimator of $\lambda_{j}^{\xi}$ is:

$$ $\lambda_{j}^{\xi,\text{emp.}}=\frac{1}{3N_{\text{samples}}MN}\sum_{k=1}^{N_{ \text{samples}}}\left|\left(\overline{\widehat{\bm{v}}}^{j,\xi}\right)^{T} \widehat{\bm{Y}_{k}}\right|^{2}.$ (129) $$

The empirical Wasserstein distance between the Gaussian distribution $\mathcal{N}({\bm{0}},\Gamma)$ and the ADSN model with texton $\mathbf{t}$ is:

$$ $\operatorname{\text{$\textbf{W}^{\text{\tiny emp.}}_{2}$}}(\mathcal{N}^{\text{ emp.}}({\bm{0}},\bm{\Gamma}),\operatorname{ADSN}(\bm{u}))=\sqrt{\sum_{\xi\in \mathbb{R}^{{M\times N}}}\left(\left(\sqrt{\lambda_{1}^{\xi,\text{emp.}}}- \sqrt{\lambda_{1}^{\xi,\operatorname{ADSN}}}\right)^{2}+\lambda_{2}^{\xi,\text {emp.}}+\lambda_{3}^{\xi,\text{emp.}}\right)}$ (130) $$

with $\lambda_{1}^{\xi,\operatorname{ADSN}}=|\widehat{\mathbf{t}}_{1}(\xi)|^{2}+|
\widehat{\mathbf{t}}_{2}(\xi)|^{2}+|\widehat{\mathbf{t}}_{3}(\xi)|^{2}$ for $\xi\in\mathbb{R}^{{M\times N}}$.

Furthermore, the computations can be vectorized by componentwise products in the Fourier domain.

## Peer Reviews

confidence: 4
contribution: 2
initial_rating: 3
presentation: 2
questions: 1. Can the authors add appropriate citations for proposition 1, for instance eq 6.19 and 6.20 in [SARKKA & SOLIN 2019]. 
2. The ADSN data samples are out of distribution for the classifier used to compute FID scores. Could the authors use a different metric? for instance MMD or other metrics/divergences. 

As mentioned in the weaknesses section, the findings of this paper are not clearly situated within the wider literature, that is what this paper is adding beyond what is known already. For instance, 

1. Error due to discretization is a widely studied problem in numerical sampling. Heun’s is a second-order integration scheme leading to slower accumulation of error compared to Euler’s method. Would the Gaussian data analysis agree with the fact that  a higher order scheme such as RK4 provide a more accurate integration of the PF-ODE?
2. In line 502-503, the authors mention that “the score approximation error is by far the most impactful source of error”. However, it is already known that score estimation error can lead to an error in sampling, for instance in [Chen et al 2023], [De Bortoli et al 2021], etc. 
3. Can the Gaussian data analysi provide any theoretical justification or insight into the robustness of SDE sampling over ODE sampling besides empirical evidence?

**References**

[Chen et al 2023] Chen, Sitan, Sinho Chewi, Jerry Li, Yuanzhi Li, Adil Salim, and Anru R. Zhang. "Sampling is as easy as learning the score: theory for diffusion models with minimal data assumptions." *arXiv preprint arXiv:2209.11215* (2022).

[De Bortoli et al 2021] De Bortoli, V., Thornton, J., Heng, J. and Doucet, A., 2021. Diffusion schrödinger bridge with applications to score-based generative modeling. Advances in Neural Information Processing Systems, 34, pp.17695-17709.
soundness: 2
strengths: In section 4, particular table 2, the authors provide a thorough analysis of the impact of the:

1. initialization error 
2. truncation error
3. discretization error 

for commonly used samplers such as EM, Heun and Euler in the centered Gaussian case.
summary: The authors define a centered Gaussian $N(0, \Sigma)$ as the data distribution and assume that the inference SDE, or forward process, is VP-SDE. VP-SDE is a linear diffusion process and it’s transition kernel, a Gaussian distribution, can be computed in closed-form. This allows the authors to

1. Compute the exact score $\nabla \log q(x_t)$ since the process $x_t$ for all time $t$ is Gaussian
2. Having access to a closed-form for the score allows one to compute exact solutions to the backward SDE (prop 2) and the probability flow ODE (prop 3), which then allows one to analyze the effect of:
    1. impact of initialization error in the ODE and SDE generative models
    2. impact of discretization error
    3. impact of truncation error
    4. influence of the eigenvalues of the covariance matrix

Using the exact derivations, the authors then analyze a Gaussian data distribution with the covariance computed from samples from the CIFAR10 dataset. The authors analyze the $W_2$ error of different samplers on this dataset and present results in table 2. 

In their next experiment, the authors learn the score function using a U-Net and see if their findings translate to when the score function architecture is nonlinear in a setting where the exact score is available as well.
weaknesses: While the authors have developed a simple framework for analyzing the different sources of error while sampling, it is limited to a Gaussian data distribution. It is not made clear how these insights generalize beyond Gaussian data and fits into the wider literature. See questions.

The score function is available in closed-form when the data distribution is a mixture of Gaussians, and as the authors mention the Wasserstein error is not computable in closed form. However, empirical metrics/divergences such as MMD, kernel Stein discrepancies, etc are available. 

Citations:

1. Proposition 1 has also been derived in appendix A in [Albergo 2023], where the authors derive the score for a number of different forward diffusion processes, flows, etc. Can the authors cite them in the main paper?
2. Similarly, for proposition 3, a solution to the PF-ODE has been derived in appendix B in [Khrulkov 2022]. Can the authors cite them in the main paper?

---

confidence: 3
contribution: 1
initial_rating: 3
presentation: 3
questions: I know that considering non-gaussian case is out of scope of this paper, but I would like to know the followings.
1. In the ablation study, the authors reported some influences based on the classified error categories (initialization/discretization/truncation errors), but it is not clear where the gaussian assumption is being used in the consideration of this error.
2. In Figure 2, as the authors commented, all errors around $\lambda=1$ drop suddenly. Why does it happen? In addition, similar sudden drops happen some other $\lambda$s. Are there any explanations on why it happens?
3. In the same page, the authors leave a comment "... for any Gaussian distribution Heun's method introduces nearly no additional discretization error, ..." Is this significantly influenced by the gaussian assumption, or can this be stated regardless of the gaussian distribution?
4. In page 9, to my understanding, the authors approximate W2 between $p_\theta$ and $p_{data}$ by replacing $p^{emp} \to p_\theta$ in eq (26). If so, how do you get $\lambda^{\xi, emp} = \lambda^{\xi, \theta}$ ? It would be related to Appendix F.3, but I'm not sure how to get $\Gamma$ in this case.
5. Besides these questions, if the authors have any insights indicating that the current results can be commonly applied outside the gaussian assumption, please make it clear.
soundness: 3
strengths: - clarity(+)
    - The paper is well written and easy to follow. The theory part is easy to understand.
    - The numerical experiments are somewhat informative.
summary: Theory:
- It is well known that the W2-distance between two gaussian distributions is exactly computable. If we consider diffusion process described by a linear drift term, it is also well known that the corresponding integral kernel is also gaussian, which means if we consider **the initial distribution as gaussian (gaussian assumption)**, the remaining all marginal distributions at diffusion time t are also gaussian distributions as commented above of eq (10). The degrees of freedom of the distribution is described by covariance matrix $\Sigma_t$.
- The authors also provide solutions of backward SDE (Proposition 2) and backward probability flow ODE (Proposition 3) under the gaussian assumption, which provide exact expressions of each covariance matrix and covariance on $t=T$ with forward time expression.
- By considering normal distribution at $t=T$ with forward time expression, the above expressions provide covariance at time t with backward SDE and ODE, which provide an inequality of W2-distance (Proposition 4).
- In addition, the authors classify errors as:
    - the initialization error, which measures discrepancy between p_T and model's latent distribution,
    - the discretization error, which measures errors caused by numerical integral of SDE or ODE,
    - the truncation error, which measures an error caused by integral truncation at t=$\epsilon>0$.

Numerical experiments:
- They conducted some experiments to measure these errors in W2-sense (Table 2) and reported (each $\Sigma_t$'s) eigenvalue dependencies of each error with "gaussianized" CIFAR-10 data (Figure 2).
- They studied effects of score approximations by models (section 5).
weaknesses: - significance (-?)
    - Thanks to the gaussian assumption, everything is very clear. However, at the same time, it makes the all results a little weak. I am not sure about how these exact results help in understanding the errors in actual learning outcomes. I summarized some questions about it below.

---

confidence: 3
contribution: 3
initial_rating: 8
presentation: 3
questions: - I was a bit surprised that closed form solutions for the forward and backwards SDEs in the Gaussian case have not been presented in the literature earlier on. My guess would have been that they could be found in SDE textbooks. The authors provide some newer references to related results at the end of section 3. My question would be if you have thoroughly checked earlier SDE literature (before diffusion models became popular in the ML world) if these results have been covered previously?
soundness: 3
strengths: - the paper presents an interesting description and investigation of the sources of errors in diffusion models
- though restricted to specific cases, I found the discussions and results interesting. I believe the findings can be useful for practical application of diffusion models
- the paper is well written and nicely presented
summary: The paper presents a theoretical investigation of the different sources of errors in diffusion models in the case of Gaussian distributions. I.e., the setting is on purpose limited to this simple case in order to be able to derive closed form solutions to the forward and backwards diffusions and to provide a direct error investigation.
weaknesses: - one can of course criticize the Gaussian assumption as being limited and not reflecting the complexity of real world data. I still believe the paper was an interesting read, even though it covers a specific case

---

confidence: 5
contribution: 3
initial_rating: 6
presentation: 3
questions: 1.The paper compares several numerical sampling schemes, but has it covered all the optimal possible schemes?
2.For complex data distributions, analytic solutions may be difficult to obtain, which may affect the computational efficiency of the model. Does the paper consider the issue of computational efficiency in practical applications?
soundness: 3
strengths: 1.The paper derives exact analytical solutions for the forward and inverse stochastic differential equations (SDEs), as well as the probabilistic flow ordinary differential equations (ODEs), within the framework of Gaussian distributions. 
2.The initialization error, truncation error, discretization error, and  score approximation error  are systematically analyzed. 
3.Different numerical sampling schemes are compared through numerical experiments, and in particular, the effectiveness of Heun's method for Gaussian distributions is demonstrated.
summary: This paper studies theoretically the behavior of diffusion models and their numerical implementation when the data distribution is Gaussian. In this restricted framework where the score function is a linear operator, we derive the analytical solutions of the backward SDE and the probability flow ODE. We prove that these solutions and their discretizations are all Gaussian processes, which allows us to compute
exact Wasserstein errors induced by each error type for any sampling scheme.
weaknesses: 1.While the paper provides an in-depth theoretical analysis, it may lack validation on diverse real-world datasets.
2.Although it is possible to derive exact solutions under the assumption of a Gaussian distribution, obtaining such analytical solutions for more complex data distributions may be challenging, which could limit the generalization of the theory. 
3.The paper points out that score approximation errors may be significant in practical applications but does not propose an effective solution for minimizing these errors.

Meta-review: The paper provided exact solutions to the diffusion model under the Gaussian distribution, and the associated error in Wasserstein distance due to the tractable form of Gaussian. However, the main drawback of the paper is that diffusion for the Gaussian distribution is the least interesting distribution from both theoretical and practical perspectives, since the main fight ground for diffusion models is for multi-modal distributions. Therefore, this paper has limited novelty and is of limited interest to the ICLR readership.

Decision: Reject