# Diffusion models for Gaussian distributions: Exact solutions and Wasserstein
errors

Emile Pierret
Institut Denis Poisson
Université d’Orléans
Orléans, France
<EMAIL_ADDRESS>
&Bruno Galerne
Institut Denis Poisson
Université d’Orléans
Orléans, France
Institut universitaire de France (IUF)
<EMAIL_ADDRESS>

###### 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 can derive the analytical solutions of the forward and backward
SDEs as well as the associated flow ODE. This provides exact expressions for
various Wasserstein errors which enable us to compare the influence of each
error type for any sampling scheme, thus allowing to monitor 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.

(a) Initialization error along the integration time

(b) Truncation error for different truncation time $\varepsilon$

Figure 1:  Wasserstein errors for the diffusion models associated with the
CIFAR-10 Gaussian. Left: Evolution of the Wasserstein distance between $p_{t}$
and the distributions associated with the continuous SDE, the continuous flow
ODE and four discrete sampling schemes with standard $\mathcal{N}_{0}$
initialization, either stochastic (Euler-Maruyama (EM) and Exponential
Integrator (EI)) or deterministic (Euler and Heun). While the continuous SDE
is less sensible than the continuous ODE (as proved by Proposition 4), the
initialization error impacts all discrete schemes with a comparable order of
magnitude. Heun’s method has the lowest error and is very close to the
theoretical ODE, except for the last step that is usually discarded when using
time truncation. Right: Wasserstein errors due to time truncation for various
truncation times $\varepsilon$. Using time truncation increases the error for
all the methods except Heun’s scheme due to instability near the origin.
Interestingly, for the standard practice truncation time
$\varepsilon=10^{-3}$, all numerical schemes have a comparable error close to
their continuous counterparts.

##  1 Introduction

Over the last five years, diffusion models have proven to be a highly
efficient and reliable framework for generative modeling [25, 14, 24, 26, 7,
15]. First introduced as a discrete process, Denoising Diffusion Probabilistic
Models (DDPM) [14] can be studied as a reversal of a continuous Stochastic
Differential Equation (SDE) [26]. 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 [26]. 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 generate
distribution is? There are four main source of errors to study for deriving
theoretical guarantees for diffusion models. First, the initialization error
is induced when approximating the marginal distribution at the end of the
forward process by a standard Gaussian distribution. Then, the discretization
error comes from the resolution of the SDE or the ODE by a numerical method.
The truncation error occurs because the backward time integration is stopped
at a small time $\varepsilon>0$ to avoid numerical instabilities due to ill-
defined score function near the origin. Finally 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 [3, 10, 15] 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 [19, 6, 4, 18, 1], eventually under an additional manifold assumption [5,
28, 2]. 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 [13].

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 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 is
Gaussian:

  * •

We give the exact solutions for both the backward SDE and probability flow
ODE.

  * •

We fully describe the Gaussian processes that occurs 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 a thorough experimental analysis 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 [15].

While our theoretical analysis relies on an exact 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. Finally, 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.

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

This preliminary section follows the seminal work of Song _et al._ [26] 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_{\mathrm{data}}$ over $\mathbb{R}^{d}$, the
forward diffusion process is the following variance preserving SDE

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

where $\left(\boldsymbol{w}_{t}\right)_{t\geq 0}$ is a $d$-dimension Brownian
motion and $\beta$ is a positive weight function. The distribution
$p_{\mathrm{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({\boldsymbol{x}}_{t}\right)$ for $t>0$ since $p_{\mathrm{data}}$ can be
supported on a lower-dimensional manifold [5]. 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_{\mathrm{data}}$ [21], the backward process
$\left({\boldsymbol{x}}_{T-t}\right)_{0\leq t\leq T}$ verifies the backward
SDE

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

The objective is now to solve this reverse equation to sample
${\boldsymbol{y}}_{T}\sim p_{\mathrm{data}}$. However, the distribution
$p_{T}$ is in general not known, and image111Although 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{\boldsymbol{y}}}_{t}=\beta_{T-t}({\tilde{\boldsymbol{y}}}_{t}+2\nabla\log
p_{T-t}({\tilde{\boldsymbol{y}}}_{t}))dt+\sqrt{2\beta_{T-t}}d\boldsymbol{w}_{t},\quad
0\leq t<T,\quad{\tilde{\boldsymbol{y}}}_{0}\sim\mathcal{N}_{0}.$ (3)

Note that approximating $p_{T}$ by $\mathcal{N}_{0}$ for the initialization
${\boldsymbol{y}}_{0}$ makes that the solution of the SDE (3) is not exactly
the target distribution $p_{\mathrm{data}}$. An alternative way to
approximately sample $p_{\mathrm{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 [26]. The deterministic process associated with Equation (2) is

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

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

$d{\boldsymbol{y}}_{t}=\left[\beta_{T-t}{\boldsymbol{y}}_{t}+\beta_{T-t}\nabla_{{\boldsymbol{y}}}\log
p_{T-t}({\boldsymbol{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{\boldsymbol{y}}}_{t}=\left[\beta_{T-t}{\widehat{\boldsymbol{y}}}_{t}+\beta_{T-t}\nabla_{{\widehat{\boldsymbol{y}}}}\log
p_{T-t}({\widehat{\boldsymbol{y}}}_{t})\right]dt,\quad 0\leq
t<T,\quad{\widehat{\boldsymbol{y}}}_{T}\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{\boldsymbol{y}}}_{t}\right)_{0\leq t\leq T}$ and
$\left({\widehat{\boldsymbol{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{\boldsymbol{y}}}_{T-t}\right)_{0\leq t\leq T}$
and $\left({\widehat{\boldsymbol{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 ${\boldsymbol{x}}_{t}$,
that we denote
${\boldsymbol{\Sigma}}_{t}=\operatorname{Cov}({\boldsymbol{x}}_{t})$.

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

The strong solution of Equation (1) can be written as:

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

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

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

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

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

For a general data distribution $p_{\mathrm{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_{\mathrm{data}}$ is a centered Gaussian distribution
$\mathcal{N}({\boldsymbol{0}},{\boldsymbol{\Sigma}})$.

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

$\nabla\log
p_{t}({\boldsymbol{x}})=-{\boldsymbol{\Sigma}}_{t}^{-1}{\boldsymbol{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 Equation (2) can be written
as:

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

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

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

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

$\operatorname{Cov}({\boldsymbol{y}}_{t})={\boldsymbol{\Sigma}}_{T-t}+e^{-2(B_{T}-B_{T-t})}{\boldsymbol{\Sigma}}_{T-t}^{2}{\boldsymbol{\Sigma}}_{T}^{-1}\left[{\boldsymbol{\Sigma}}^{-1}_{T}\operatorname{Cov}({\boldsymbol{y}}_{0})-\boldsymbol{I}\right]$
(13)

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. As shown by
the following proposition (proved in Appendix B.3), the flow ODE also has an
explicit solution under Gaussian assumption which related to optimal
transport.

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

The solution to the probability flow ODE (4) under Gaussian assumption
corresponds to the optimal transport map between $p_{T}$ and
$p_{\mathrm{data}}$. More precisely, for any ${\boldsymbol{y}}_{0}$,

${\boldsymbol{y}}_{t}={\boldsymbol{\Sigma}}^{-1/2}_{T}{\boldsymbol{\Sigma}}^{1/2}_{T-t}{\boldsymbol{y}}_{0},\quad
0\leq t\leq T,$

is the solution of the reverse-time ODE (5). Consequently, the covariance
matrix $\operatorname{Cov}({\boldsymbol{y}}_{t})$ verifies

$\operatorname{Cov}({\boldsymbol{y}}_{t})={\boldsymbol{\Sigma}}^{-1/2}_{T}{\boldsymbol{\Sigma}}^{1/2}_{T-t}\operatorname{Cov}({\boldsymbol{y}}_{0}){\boldsymbol{\Sigma}}^{1/2}_{T-t}{\boldsymbol{\Sigma}}^{-1/2}_{T},\quad
0\leq t\leq T,$ (14)

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

$\operatorname{Cov}({\boldsymbol{y}}_{t})={\boldsymbol{\Sigma}}^{-1}_{T}{\boldsymbol{\Sigma}}_{T-t}\operatorname{Cov}({\boldsymbol{y}}_{0}),\quad
0\leq t\leq T.$ (15)

Here we must highlight a subtle issue: Whatever the initial distribution is,
the ODE solution consists in applying the optimal transport map between
$p_{T}$ and $p_{\mathrm{data}}$. Since in practice one cannot truly sample
$p_{T}$, the resulting flow is not an optimal transport flow and the
distribution of ${\boldsymbol{y}}_{T}$ differs from $p_{\mathrm{data}}$.

#### 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 [12]
and the ODE verified by the variance is given in [26] (citing [22]) but it is
generalized here to the full covariance matrix (Equation (9)). To the best of
our knowledge, Proposition 2 is new and one the most important contribution of
the paper. Similar to our approach, Gaussian mixtures have been studied in the
context of diffusion model [29, 30, 23] since they also provide an explicit
analytical score. However, solving the backward SDE is not feasible for
Gaussian mixtures as far as we know. The ODE [17, 16] can be interpreted in
the infinite time as an optimal transport (OT) between the prior distribution
and the Gaussian standard distribution. The relation between optimal transport
and probability flow ODE (also called Fokker-Planck ODE) has been discussed in
[17, 16] in the asymptotic case where $T\mapsto+\infty$. Both these papers
discussed the Gaussian case, but our Proposition 3 highlights that the
generated process is not an optimal transport flow due to the intialiazation
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[8]. For two centered Gaussians
$\mathcal{N}({\boldsymbol{0}},{\boldsymbol{\Sigma}}_{1})$ and
$\mathcal{N}({\boldsymbol{0}},{\boldsymbol{\Sigma}}_{2})$ such that
${\boldsymbol{\Sigma}}_{1},{\boldsymbol{\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}({\boldsymbol{0}},{\boldsymbol{\Sigma}}_{1}),\mathcal{N}({\boldsymbol{0}},{\boldsymbol{\Sigma}}_{2}))^{2}=\sum_{1\leq
i\leq d}(\sqrt{\lambda_{i,1}}-\sqrt{\lambda_{i,2}})^{2}.$ (16)

In the literature, the quality of the diffusion models is measured with FID
[13] which is the $\operatorname{\text{$\textbf{W}_{2}$}}2$-error between
Gaussians fitted to the Inception features [27] 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 ${\boldsymbol{\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{\boldsymbol{y}}}$ and ${\widehat{\boldsymbol{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({\tilde{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{\boldsymbol{y}}}_{t}\right)_{0\leq
t\leq T}$ and $\left({\widehat{\boldsymbol{y}}}_{t}\right)_{0\leq t\leq T}$
are Gaussian processes. At each time $t$, ${\tilde{p}}_{t}$ is the Gaussian
distribution
$\mathcal{N}({\boldsymbol{0}},{\tilde{{\boldsymbol{\Sigma}}}}_{t})$ with
${\tilde{{\boldsymbol{\Sigma}}}}_{t}={\boldsymbol{\Sigma}}_{t}+e^{-2(B_{T}-B_{t})}{\boldsymbol{\Sigma}}^{2}_{t}{\boldsymbol{\Sigma}}_{T}^{-1}({\boldsymbol{\Sigma}}^{-1}_{T}-\boldsymbol{I})$
and ${\widehat{p}}_{t}$ is the Gaussian distribution
$\mathcal{N}({\boldsymbol{0}},{\widehat{{\boldsymbol{\Sigma}}}}_{t})$ with
${\widehat{{\boldsymbol{\Sigma}}}}_{t}={\boldsymbol{\Sigma}}^{-1}_{T}{\boldsymbol{\Sigma}}_{t}$.
For all $0\leq t\leq T$, the three covariance matrices
${\boldsymbol{\Sigma}}_{t}$, ${\tilde{{\boldsymbol{\Sigma}}}}_{t}$ and
${\widehat{{\boldsymbol{\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})$
(17)

which shows for $t=0$ that 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).

#### 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 [26] and the exponential
integrator (EI) in [5] to sample from the SDE (3). The Euler method is the
simplest ODE solver and Heun’s scheme is recommended in [15] to model the ODE
(6). Under Gaussian assumption, the eigenvalues of the covariance matrix can
be computed numerically recursively for each scheme to evaluate the
Wasserstein distance. More precisely, all the covariance matrices are
diagonalizable in the diagonalization basis of ${\boldsymbol{\Sigma}}$, and
thus Equation (16) is valid. For example, denoting
$\left(\lambda^{t}_{i}\right)_{1\leq i\leq d}$ the eigenvalues of
${\boldsymbol{\Sigma}}_{t}$ and $\big{(}\lambda^{EM,k}_{i}\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^{\text{EM},k+1}_{i}=\big{(}1+\Delta_{t}\beta_{T-t_{k}}(1-\tfrac{2}{\lambda^{T-t_{k}}_{i}})\big{)}^{2}\lambda^{\text{EM},k}_{i}+2\Delta_{t}\beta_{T-t_{k}},1\leq
i\leq d,0\leq k\leq N-2$ (18)

with initialization $\lambda^{0,\text{EM}}_{i}=1,1\leq i\leq 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.

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

#### The truncation error.

As discussed in [26], 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 ${\boldsymbol{\Sigma}}$ is not
invertible. This approximation is called the truncation error. As a
consequence, even without error initialization, 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 5.(b) for both continuous and numerical solutions. For the
standard practice truncation time $\varepsilon=10^{-3}$ [26, 15], 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.

|  | 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-4 |  4.77 |  4.77 |  0.65 |  0.65 |  0.31 |  0.31 |  0.15 |  0.16
$\varepsilon=10^{-5}$ |  2.5E-3 |  2.6E-3 |  4.77 |  4.77 |  0.65 |  0.65 |  0.31 |  0.31 |  0.16 |  0.16
$\varepsilon=10^{-3}$ |  0.17 |  0.17 |  4.67 |  4.67 |  0.69 |  0.69 |  0.39 |  0.39 |  0.27 |  0.27
$\varepsilon=10^{-2}$ |  1.35 |  1.35 |  4.56 |  4.56 |  1.69 |  1.69 |  1.50 |  1.50 |  1.42 |  1.42
EI | $\varepsilon=0$ |  0 |  6.7E-4 |  2.81 |  2.81 |  0.57 |  0.57 |  0.30 |  0.30 |  0.16 |  0.16
$\varepsilon=10^{-5}$ |  2.5E-3 |  2.6E-3 |  2.81 |  2.81 |  0.57 |  0.57 |  0.30 |  0.30 |  0.16 |  0.16
$\varepsilon=10^{-3}$ |  0.17 |  0.17 |  2.91 |  2.91 |  0.66 |  0.66 |  0.41 |  0.41 |  0.28 |  0.28
$\varepsilon=10^{-2}$ |  1.35 |  1.35 |  3.93 |  3.93 |  1.76 |  1.76 |  1.55 |  1.55 |  1.45 |  1.45
Euler | $\varepsilon=0$ |  0 |  0.07 |  1.72 |  1.78 |  0.38 |  0.44 |  0.19 |  0.26 |  0.10 |  0.17
$\varepsilon=10^{-5}$ |  2.5E-3 |  0.07 |  1.72 |  1.78 |  0.38 |  0.44 |  0.20 |  0.26 |  0.10 |  0.17
$\varepsilon=10^{-3}$ |  0.17 |  0.19 |  1.72 |  1.78 |  0.42 |  0.48 |  0.27 |  0.32 |  0.21 |  0.25
$\varepsilon=10^{-2}$ |  1.35 |  1.36 |  2.21 |  2.25 |  1.41 |  1.43 |  1.37 |  1.38 |  1.36 |  1.37
Heun | $\varepsilon=0$ |  0 |  0.07 |  7.09 |  7.09 |  0.72 |  0.73 |  0.21 |  0.22 |  0.05 |  0.09
$\varepsilon=10^{-5}$ |  2.5E-3 |  0.07 |  6.48 |  6.48 |  0.64 |  0.65 |  0.18 |  0.20 |  0.05 |  0.09
$\varepsilon=10^{-3}$ |  0.17 |  0.19 |  0.56 |  0.57 |  0.13 |  0.15 |  0.16 |  0.18 |  0.17 |  0.19
$\varepsilon=10^{-2}$ |  1.35 |  1.36 |  1.37 |  1.38 |  1.35 |  1.36 |  1.35 |  1.36 |  1.35 |  1.36
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).

#### 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^{-2}$. 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^{-5}$. As [15], our conclusions lead to the choice of Heun’s
scheme as the go-to method.

(a) Initialization error at final time

(b) Truncation error at final time for $\varepsilon=10^{-3}$

Figure 2:  Eigenvalue contribution to the Wasserstein error. The magnitude of
the Wasserstein error is influenced by the eigenvalues of the covariance of
the Gaussian distribution. Left: Contribution to the Wasserstein error for the
continuous equations and the discretization schemes with standard
initialization $\mathcal{N}_{0}$. Right: Same plot when using a truncation
time $\varepsilon=10^{-3}$. All schemes use $N=1000$ steps. While we prove
that the continuous SDE is always better than the continuous ODE (Proposition
4), it is not the same for the discrete schemes. With a truncation time
$\varepsilon=10^{-3}$ (b), Heun’s method is nearly as good as the continuous
ODE solution for all eigenvalues, which shows it is well-adapted to any
Gaussian distribution.

#### 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
${\boldsymbol{\Sigma}}$. Indeed, as seen in Equation (16), 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. 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 ${\boldsymbol{\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.

##  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 [5]), 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 [11] associated with an RGB texture
$\boldsymbol{u}\in\mathbb{R}^{3\times{M\times N}}$ which is defined as the
stationary Gaussian distribution that has covariance equal the autocorrelation
of $\boldsymbol{u}$. More precisely, this distribution is sampled using
convolution with a white Gaussian noise [11]: Denoting $m\in\mathbb{R}^{3}$
the chanelwise mean of $\boldsymbol{u}$ and
$\mathbf{t}_{c}=\frac{1}{\sqrt{MN}}(\boldsymbol{u}_{c}-m_{c})$, $1\leq c\leq
3$, its associated _texton_ , for $\boldsymbol{w}\sim\mathcal{N}_{0}$ of size
$M\times N$ the chanelwise convolution
${\boldsymbol{x}}=m+\mathbf{t}\star\boldsymbol{w}\in\mathbb{R}^{3\times{M\times
N}}$ follows $\operatorname{ADSN}(\boldsymbol{u})$. This distribution is the
Gaussian $\mathcal{N}(m,{\boldsymbol{\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}({\boldsymbol{0}},{\boldsymbol{\Sigma}})$.
The matrix ${\boldsymbol{\Sigma}}$ is a well-known convolution matrix [9], its
eigenvectors and associated eigenvalues can be computed in the Fourier domain,
as done in Appendix E.2. ${\boldsymbol{\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$ with multiplicity $2MN$ and we can conduct the same analysis as
before (see Appendix D). To evaluate if a set of $N_{\text{samples}}$ sampled
images is close to the ADSN distribution $p_{\mathrm{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}({\boldsymbol{0}},\boldsymbol{\Gamma})$ such that
$\boldsymbol{\Gamma}$ admits the same eigenvectors as ${\boldsymbol{\Sigma}}$,
we compute

$\operatorname{\text{$\textbf{W}^{\text{\tiny emp.}}_{2}$}}(p^{\text{\tiny
emp.}},p_{\mathrm{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}}$
(23)

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

Original $\boldsymbol{u}$ |  |  |  |  |  |
---|---|---|---|---|---|---
ADSN |  |  |  |  |  |
$p_{\theta}^{\text{EM}}$ |  |  |  |  |  |
$p_{\theta}^{\text{Heun}}$ |  |  |  |  |  |
Figure 3:  Texture samples generated with the learned score. First row:
original image $\boldsymbol{u}$ and its DFT modulus (for all DFT modulus we
display the sum of the DFT modulus of the three color chanels and apply a
logarithmic contrast change). Second row: three samples of
$\operatorname{ADSN}(\boldsymbol{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.

#### Learning the score function.

We train the network using the code222Code available at
https://github.com/yang-song/score_sde_pytorch associated with the paper [26].
We choose the architecture of DDPM, which is a U-Net described in [14], 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 [26]. $\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_{\mathrm{data}}=\operatorname{ADSN}(\boldsymbol{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 (23). 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
${\boldsymbol{\Sigma}}$ of $\operatorname{ADSN}(\boldsymbol{u})$. We complete
this dedicated empirical measure with the standard FID. This 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 consider as some initialization error
ocuring at each step.

| Exact score distribution | Learned score distribution
---|---|---
$p$ |  $\operatorname{\text{$\textbf{W}_{2}$}}2$($p$,$p_{\mathrm{data}}$) $\downarrow$ |  $\operatorname{\text{$\textbf{W}^{\text{\tiny emp.}}_{2}$}}$($p^{\text{\tiny emp.}}$,$p_{\mathrm{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
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 4. 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).

##  6 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. 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 [15].

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.

Acknowledgements: The authors acknowledge the support of the project MISTIC
(ANR-19-CE40-005).

## References

  * [1] Joe Benton, Valentin De Bortoli, Arnaud Doucet, and George Deligiannidis.  Nearly $d$-linear convergence bounds for diffusion models via stochastic localization.  In The Twelfth International Conference on Learning Representations, 2024.
  * [2] Minshuo Chen, Kaixuan Huang, Tuo Zhao, and Mengdi Wang.  Score approximation, estimation and distribution recovery of diffusion models on low-dimensional data.  In Andreas Krause, Emma Brunskill, Kyunghyun Cho, Barbara Engelhardt, Sivan Sabato, and Jonathan Scarlett, editors, Proceedings of the 40th International Conference on Machine Learning, volume 202 of Proceedings of Machine Learning Research, pages 4672–4712. PMLR, 23–29 Jul 2023.
  * [3] Sitan Chen, Sinho Chewi, Holden Lee, Yuanzhi Li, Jianfeng Lu, and Adil Salim.  The probability flow ODE is provably fast.  In Thirty-seventh Conference on Neural Information Processing Systems, 2023.
  * [4] Sitan Chen, Sinho Chewi, Jerry Li, Yuanzhi Li, Adil Salim, and Anru Zhang.  Sampling is as easy as learning the score: theory for diffusion models with minimal data assumptions.  In The Eleventh International Conference on Learning Representations, 2023.
  * [5] Valentin De Bortoli.  Convergence of denoising diffusion models under the manifold hypothesis.  Transactions on Machine Learning Research, 2022.
  * [6] Valentin De Bortoli, James Thornton, Jeremy Heng, and Arnaud Doucet.  Diffusion schrödinger bridge with applications to score-based generative modeling.  In Advances in Neural Information Processing Systems, volume 34, pages 17695–17709. Curran Associates, Inc., 2021.
  * [7] Prafulla Dhariwal and Alexander Nichol.  Diffusion models beat GANs on image synthesis.  In Advances in Neural Information Processing Systems, volume 34, pages 8780–8794. Curran Associates, Inc., 2021.
  * [8] D.C Dowson and B.V Landau.  The fréchet distance between multivariate normal distributions.  Journal of Multivariate Analysis, 12(3):450–455, 1982.
  * [9] Sira Ferradans, Gui-Song Xia, Gabriel Peyré, and Jean-François Aujol.  Static and dynamic texture mixing using optimal transport.  In Scale Space and Variational Methods in Computer Vision, 2013.
  * [10] Giulio Franzese, Simone Rossi, Lixuan Yang, Alessandro Finamore, Dario Rossi, Maurizio Filippone, and Pietro Michiardi.  How much is enough? a study on diffusion times in score-based generative models.  Entropy, 25(4), 2023.
  * [11] Bruno Galerne, Yann Gousseau, and Jean-Michel Morel.  Random Phase Textures: Theory and Synthesis.  IEEE Transactions on Image Processing, 20(1):257–267, 2011.
  * [12] Xuefeng Gao and Lingjiong Zhu.  Convergence analysis for general probability flow odes of diffusion models in wasserstein distances.  ArXiv, abs/2401.17958, 2024.
  * [13] Martin Heusel, Hubert Ramsauer, Thomas Unterthiner, Bernhard Nessler, and Sepp Hochreiter.  Gans trained by a two time-scale update rule converge to a local Nash equilibrium.  In I. Guyon, U. Von Luxburg, S. Bengio, H. Wallach, R. Fergus, S. Vishwanathan, and R. Garnett, editors, Advances in Neural Information Processing Systems, volume 30. Curran Associates, Inc., 2017.
  * [14] Jonathan Ho, Ajay Jain, and Pieter Abbeel.  Denoising diffusion probabilistic models.  In Advances in Neural Information Processing Systems, volume 33, pages 6840–6851. Curran Associates, Inc., 2020.
  * [15] Tero Karras, Miika Aittala, Timo Aila, and Samuli Laine.  Elucidating the design space of diffusion-based generative models.  In Proc. NeurIPS, 2022.
  * [16] Valentin Khrulkov, Gleb Ryzhakov, Andrei Chertkov, and Ivan Oseledets.  Understanding DDPM latent codes through optimal transport.  In The Eleventh International Conference on Learning Representations, 2023.
  * [17] Hugo Lavenant and Filippo Santambrogio.  The flow map of the fokker–planck equation does not provide optimal transport.  Applied Mathematics Letters, 133:108225, 2022.
  * [18] Holden Lee, Jianfeng Lu, and Yixin Tan.  Convergence of score-based generative modeling for general data distributions.  In NeurIPS 2022 Workshop on Score-Based Methods, 2022.
  * [19] Holden Lee, Jianfeng Lu, and Yixin Tan.  Convergence for score-based generative modeling with polynomial complexity.  In Proceedings of the 36th International Conference on Neural Information Processing Systems, NIPS ’22, Red Hook, NY, USA, 2024. Curran Associates Inc.
  * [20] B. Øksendal.  Stochastic Differential Equations: An Introduction with Applications.  Universitext. Springer Berlin Heidelberg, 2010.
  * [21] E. Pardoux.  Grossissement d’une filtration et retournement du temps d’une diffusion.  In Jacques Azéma and Marc Yor, editors, Séminaire de Probabilités XX 1984/85, pages 48–55, Berlin, Heidelberg, 1986. Springer Berlin Heidelberg.
  * [22] S. Särkkä and A. Solin.  Applied Stochastic Differential Equations.  Institute of Mathematical Statistics Textbooks. Cambridge University Press, 2019.
  * [23] Kulin Shah, Sitan Chen, and Adam Klivans.  Learning mixtures of gaussians using the DDPM objective.  In Thirty-seventh Conference on Neural Information Processing Systems, 2023.
  * [24] Jiaming Song, Chenlin Meng, and Stefano Ermon.  Denoising diffusion implicit models.  In International Conference on Learning Representations, 2021.
  * [25] Yang Song and Stefano Ermon.  Generative modeling by estimating gradients of the data distribution.  In H. Wallach, H. Larochelle, A. Beygelzimer, F. d’ Alché-Buc, E. Fox, and R. Garnett, editors, Advances in Neural Information Processing Systems, volume 32. Curran Associates, Inc., 2019.
  * [26] Yang Song, Jascha Sohl-Dickstein, Diederik P Kingma, Abhishek Kumar, Stefano Ermon, and Ben Poole.  Score-based generative modeling through stochastic differential equations.  In International Conference on Learning Representations, 2021.
  * [27] Christian Szegedy, Vincent Vanhoucke, Sergey Ioffe, Jon Shlens, and Zbigniew Wojna.  Rethinking the inception architecture for computer vision.  In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition (CVPR), June 2016.
  * [28] Li Kevin Wenliang and Ben Moran.  Score-based generative model learn manifold-like structures with constrained mixing.  In NeurIPS 2022 Workshop on Score-Based Methods, 2022.
  * [29] Martin Zach, Erich Kobler, Antonin Chambolle, and Thomas Pock.  Product of gaussian mixture diffusion models.  Journal of Mathematical Imaging and Vision, Mar 2024.
  * [30] Martin Zach, Thomas Pock, Erich Kobler, and Antonin Chambolle.  Explicit diffusion of gaussian mixture model based image priors.  In Luca Calatroni, Marco Donatelli, Serena Morigi, Marco Prato, and Matteo Santacesaria, editors, Scale Space and Variational Methods in Computer Vision, pages 3–15, Cham, 2023. Springer International Publishing.

##  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:

  1. (i)

${\boldsymbol{x}}_{0}\sim\mathcal{N}({\boldsymbol{0}},{\boldsymbol{\Sigma}})$
for some covariance ${\boldsymbol{\Sigma}}$.

  2. (ii)

$\forall t>0,\nabla_{x}\log p_{t}({\boldsymbol{x}})$ is linear w.r.t
${\boldsymbol{x}}$.

  3. (iii)

$\exists t>0,\nabla_{x}\log p_{t}({\boldsymbol{x}})$ is linear w.r.t
${\boldsymbol{x}}$.

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

###### Proof.

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

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

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

##  Appendix B Proofs of Section 3

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

We aim at solving:

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

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

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

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

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

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

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

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

$\text{Var}\left(\int_{0}^{t}e^{B_{s}}\sqrt{2\beta_{s}}d\boldsymbol{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$
(28)

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

$\operatorname{Cov}\left({\boldsymbol{\eta}}_{t}\right)=e^{-2B_{t}}(e^{2B_{t}}-1)\boldsymbol{I}=\left(1-e^{-2B_{t}}\right)\boldsymbol{I}.$
(29)

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

And,

$d{\boldsymbol{\Sigma}}_{t}=-2\beta_{t}e^{-2B_{t}}({\boldsymbol{\Sigma}}-\boldsymbol{I})dt=\left[-2\beta_{t}{\boldsymbol{\Sigma}}_{t}+2\beta_{t}\boldsymbol{I}\right]dt=-(2\beta_{t}{\boldsymbol{\Sigma}}_{t}-2\beta_{t}{\boldsymbol{\Sigma}}^{-1}_{T}){\boldsymbol{\Sigma}}_{t}dt.$
(30)

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

We aim at solving

$d{\boldsymbol{y}}_{t}=\beta_{T-t}({\boldsymbol{y}}_{t}+2\nabla\log
p_{T-t}({\boldsymbol{y}}_{t}))dt+\sqrt{2\beta_{T-t}}d\boldsymbol{w}_{t},\quad
0\leq t\leq T$ (31)

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

$\displaystyle d{\boldsymbol{z}}_{t}$
$\displaystyle=e^{C_{t}}{\boldsymbol{\Sigma}}^{-1}_{T-t}d{\boldsymbol{y}}_{t}-e^{C_{t}}d[{\boldsymbol{\Sigma}}^{-1}]_{T-t}{\boldsymbol{y}}_{t}dt+\beta_{T-t}{\boldsymbol{z}}_{t}dt$
(32)
$\displaystyle=\left[{\boldsymbol{\Sigma}}^{-1}_{T-t}e^{C_{t}}\beta_{T-t}({\boldsymbol{y}}_{t}-2{\boldsymbol{\Sigma}}^{-1}_{T-t}{\boldsymbol{y}}_{t})-\beta_{T-t}{\boldsymbol{z}}_{t}+2\beta_{T-t}{\boldsymbol{\Sigma}}^{-1}_{T-t}{\boldsymbol{z}}_{t}\right]dt+\sqrt{2\beta_{T-t}}e^{C_{t}}{\boldsymbol{\Sigma}}^{-1}_{T-t}d\boldsymbol{w}_{t}$
(33) (using Equation (9)) (34)
$\displaystyle=\beta_{T-t}(1-2{\boldsymbol{\Sigma}}^{-1}_{T-t}){\boldsymbol{z}}_{t}dt-\beta_{T-t}{\boldsymbol{z}}_{t}dt+2\beta_{T-t}{\boldsymbol{\Sigma}}^{-1}_{T-t}{\boldsymbol{z}}_{t}dt+e^{C_{t}}\sqrt{2\beta_{T-t}}{\boldsymbol{\Sigma}}^{-1}_{T-t}d\boldsymbol{w}_{t}$
(35)
$\displaystyle=\sqrt{2\beta_{T-t}}e^{C_{t}}{\boldsymbol{\Sigma}}^{-1}_{T-t}d\boldsymbol{w}_{t}.$
(36)

Consequently,

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

And,

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

Finally,

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

By the multidimensional Itô’s isometry,

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

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

$\displaystyle dA_{s}$
$\displaystyle=2\beta_{T-s}A_{s}ds-e^{2C_{s}}d\left[{\boldsymbol{\Sigma}}^{-1}\right]_{T-s}$
(42)
$\displaystyle=2\beta_{T-s}A_{s}ds+e^{2C_{s}}\left[-2\beta_{T-s}{\boldsymbol{\Sigma}}^{-1}_{T-s}+2\beta_{T-s}{\boldsymbol{\Sigma}}^{-2}_{T-s}\right]ds$
(43) $\displaystyle=2e^{2C_{s}}\beta_{T-s}{\boldsymbol{\Sigma}}^{-2}_{T-s}ds.$
(44)

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

Finally,
$\operatorname{Cov}({\boldsymbol{\xi}}_{t})={\boldsymbol{\Sigma}}_{T-t}^{2}\left({\boldsymbol{\Sigma}}^{-1}_{T-t}-e^{-2C_{t}}{\boldsymbol{\Sigma}}^{-1}_{T}\right)={\boldsymbol{\Sigma}}_{T-t}-e^{-2C_{t}}{\boldsymbol{\Sigma}}_{T-t}^{2}{\boldsymbol{\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}$
(46)

that provides

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

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

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

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

As done in [16], the matrix ${\boldsymbol{\Sigma}}^{1/2}_{t}$ admits a
derivative which is
$d\left[{\boldsymbol{\Sigma}}^{1/2}\right]_{t}=\frac{1}{2}d{\boldsymbol{\Sigma}}_{t}{\boldsymbol{\Sigma}}^{-1/2}_{t}$
because it is diagonalisable. Let us check that
${\boldsymbol{y}}_{t}={\boldsymbol{\Sigma}}^{-1/2}_{T}{\boldsymbol{\Sigma}}_{T-t}^{1/2}{\boldsymbol{y}}_{0}$
is solution of the ODE (4),

$\displaystyle d{\boldsymbol{y}}_{t}$
$\displaystyle=-{\boldsymbol{\Sigma}}^{-1/2}_{T}\frac{1}{2}d{\boldsymbol{\Sigma}}_{T-t}{\boldsymbol{\Sigma}}^{-1/2}_{t}{\boldsymbol{y}}_{0}dt$
(49)
$\displaystyle={\boldsymbol{\Sigma}}^{-1/2}_{T}\left[\beta_{T-t}{\boldsymbol{\Sigma}}_{T-t}-2\beta_{T-t}\boldsymbol{I}\right]{\boldsymbol{\Sigma}}^{-1/2}_{T-t}{\boldsymbol{y}}_{0}dt\quad(\text{using
Equation\eqref{eq:ode_covariance}})$ (50)
$\displaystyle=\left[\beta_{T-t}-2\beta_{T-t}{\boldsymbol{\Sigma}}^{-1}_{T-t}\right]{\boldsymbol{y}}_{t}dt$
(51)
$\displaystyle=\left[\beta_{T-t}+2\beta_{T-t}\nabla_{{\boldsymbol{y}}}\log
p_{T-t}({\boldsymbol{y}}_{t})\right]{\boldsymbol{y}}_{t}dt.$ (52)

Finally,

$\operatorname{Cov}({\boldsymbol{y}}_{t})={\boldsymbol{\Sigma}}^{-1/2}_{T}{\boldsymbol{\Sigma}}_{T-t}^{1/2}\operatorname{Cov}({\boldsymbol{y}}_{0}){\boldsymbol{\Sigma}}^{-1/2}_{T-t}{\boldsymbol{\Sigma}}_{T}^{1/2}.$
(53)

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

$\operatorname{Cov}({\boldsymbol{y}}_{t})={\boldsymbol{\Sigma}}^{-1}_{T}{\boldsymbol{\Sigma}}_{T-t}\operatorname{Cov}({\boldsymbol{y}}_{0}).$
(54)

###  B.4 Proof of Proposition 4

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

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

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

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

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

$\sqrt{\lambda^{t}_{i}}-\sqrt{\tilde{\lambda}^{t}_{i}}=\sqrt{\lambda^{t}_{i}}\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).$
(57)

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

$\displaystyle\sqrt{\lambda^{t}_{i}}-\sqrt{\widehat{\lambda}^{t}_{i}}=\sqrt{\lambda^{t}_{i}}\left(1-\sqrt{\frac{1}{\lambda_{i}^{T}}}\right)$
(58)
$\displaystyle=\sqrt{\lambda^{t}_{i}}\left(1-\sqrt{1+\left(\frac{1}{\lambda_{i}^{T}}-1\right)}\right).$
(59)

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}}$
(60)
$\displaystyle=e^{-2(B_{T}-B_{t})}\frac{e^{-2B_{t}}(\lambda_{i}-1)+1}{e^{-2B_{T}}(\lambda_{i}-1)+1}$
(61)
$\displaystyle=\frac{(\lambda_{i}-1)+e^{2B_{t}}}{(\lambda_{i}-1)+e^{2B_{T}}}$
(62) $\displaystyle<1.$ (63)

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.$
(64)

Thus,

$\displaystyle\left|\sqrt{\lambda^{t}_{i}}-\sqrt{\tilde{\lambda}^{t}_{i}}\right|$
$\displaystyle=\sqrt{\tilde{\lambda}^{t}_{i}}-\sqrt{\lambda^{t}_{i}}$ (65)
$\displaystyle=\sqrt{\lambda^{t}_{i}}\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)$
(66)
$\displaystyle<\sqrt{\lambda^{t}_{i}}\left(\sqrt{1+\left(\frac{1}{\lambda_{i}^{T}}-1\right)}-1\right)$
(67) $\displaystyle=\sqrt{\widehat{\lambda}^{t}_{i}}-\sqrt{\lambda^{t}_{i}}$
(68)
$\displaystyle=\left|\sqrt{\lambda^{t}_{i}}-\sqrt{\widehat{\lambda}^{t}_{i}}\right|.$
(69)

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

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

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

In this case, for $1\leq i\leq d$,
$\widehat{\lambda}^{t}_{i}=\tilde{\lambda}^{t}_{i}=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.$
(70)

Finally,

$\displaystyle\left|\sqrt{\lambda^{t}_{i}}-\sqrt{\tilde{\lambda}^{t}_{i}}\right|$
$\displaystyle=\sqrt{\lambda^{t}_{i}}-\sqrt{\tilde{\lambda}^{t}_{i}}$ (71)
$\displaystyle=\sqrt{\lambda^{t}_{i}}\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)$
(72)
$\displaystyle<\sqrt{\lambda^{t}_{i}}\left(1-\sqrt{1+\left(\frac{1}{\lambda_{i}^{T}}-1\right)}\right)$
(73) $\displaystyle=\sqrt{\lambda^{t}_{i}}-\sqrt{\widehat{\lambda}^{t}_{i}}$
(74)
$\displaystyle=\left|\sqrt{\lambda^{t}_{i}}-\sqrt{\widehat{\lambda}^{t}_{i}}\right|.$
(75)

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}).$
(76)

##  Appendix C Gaussian CIFAR-10 samples

The Gaussian CIFAR-10 produces unstructured images. A grid of samples is
presented in Figure 4.

Figure 4:  CIFAR-10 Gaussian samples. Samples are generated from the Gaussian
distribution fitting the CIFAR-10 dataset.

##  Appendix D 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 4.

(a) Initialization error along the integration time.

(b) Truncation error for different truncation time $\varepsilon$.

Figure 5:  Wasserstein errors for the diffusion models associated with the Gaussian microtextures. Left: Evolution of the Wasserstein distance between $p_{t}$ and the distributions associated with the continuous SDE, the continuous flow ODE and four discrete sampling schemes with standard $\mathcal{N}_{0}$ initialization, either stochastic (Euler-Maruyama (EM) and Exponential Integrator (EI)) or deterministic (Euler and Heun). While the continuous SDE is less sensible than the continuous ODE (as proved by Proposition 4), the initialization error impacts all discrete schemes. Heun’s method has the lowest error and is very close to the theoretical ODE, except for the last step that is usually discarded when using time truncation. Right: Wasserstein errors due to time truncation for various truncation times $\varepsilon$. Heun’s scheme is not defined without truncation time due to the zero eigenvalue. Interestingly, for the standard practice truncation time $\varepsilon=10^{-3}$, all numerical schemes have a comparable error close to their continuous counterparts. |  | 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-6 | 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^{-3}$ | 3.84 | 3.84 | 51.92 | 51.92 | 10.55 | 10.55 | 6.80 | 6.80 | 5.16 | 5.16
$\varepsilon=10^{-2}$ | 17.09 | 17.09 | 48.24 | 48.24 | 20.39 | 20.39 | 18.57 | 18.57 | 17.79 | 17.79
EI | $\varepsilon=0$ | 0 | 5.2E-6 | 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^{-3}$ | 3.84 | 3.84 | 31.94 | 31.94 | 9.74 | 9.74 | 6.76 | 6.76 | 5.24 | 5.24
$\varepsilon=10^{-2}$ | 17.09 | 17.09 | 41.49 | 41.49 | 21.02 | 21.02 | 18.95 | 18.95 | 17.99 | 17.99
Euler | $\varepsilon=0$ | 0 | 6.4E-3 | 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^{-3}$ | 3.84 | 3.84 | 6.79 | 6.79 | 4.85 | 4.85 | 4.41 | 4.41 | 4.14 | 4.14
$\varepsilon=10^{-2}$ | 17.09 | 17.09 | 18.52 | 18.52 | 17.35 | 17.35 | 17.22 | 17.22 | 17.15 | 17.15
Heun | $\varepsilon=0$ | 0 | 6.4E-3 | - | - | - | - | - | - | - | -
$\varepsilon=10^{-5}$ | 0.36 | 0.36 | 2.4E+3 | 2.4E+3 | 3.0E+2 | 3.0E+2 | 1.1E+2 | 1.1E+2 | 40.00 | 40.00
$\varepsilon=10^{-3}$ | 3.84 | 3.84 | 15.42 | 15.42 | 2.25 | 2.25 | 3.40 | 3.40 | 3.73 | 3.73
$\varepsilon=10^{-2}$ | 17.09 | 17.09 | 16.59 | 16.59 | 17.07 | 17.07 | 17.09 | 17.09 | 17.09 | 17.09
Table 4:  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.

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

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

For a given image $\boldsymbol{v}\in\mathbb{R}^{3\times{M\times N}}$, we
define the DFT of $\boldsymbol{v}$,
$\widehat{\boldsymbol{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{\boldsymbol{v}}_{c,\xi}=\sum_{x\in{M\times
N}}\boldsymbol{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$ (77)

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

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

we have:

$\widehat{\boldsymbol{u}\star\boldsymbol{v}}=\widehat{\boldsymbol{u}}\odot\widehat{\boldsymbol{v}},$
(79)

where $\odot$ is the componentwise product.

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

Let $\boldsymbol{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}(\boldsymbol{u})$ is the Gaussian distribution of
$\boldsymbol{X}=\mathbf{t}\star\boldsymbol{w}$ such that:

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

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

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

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

$\widehat{\boldsymbol{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}}$ (82)

we have,

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

Furthermore, the family $\left(\boldsymbol{u}^{1,\xi}\right)_{\xi\in{M\times
N}}$ is orthogonal. Thus, the eigenvalues of ${\boldsymbol{\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$ with multiplicity $2MN$.

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

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

We have

$\displaystyle{\boldsymbol{\Sigma}}\boldsymbol{u}^{2,\xi}$
$\displaystyle=0.\boldsymbol{u}^{2,\xi}$ (92)
$\displaystyle{\boldsymbol{\Sigma}}\boldsymbol{u}^{3,\xi}$
$\displaystyle=0.\boldsymbol{u}^{3,\xi}.$ (93)

Then, applying the orthonomalization of Gram-Schmidt on each tuple
$(\boldsymbol{u}^{1,\xi},\boldsymbol{u}^{2,\xi},\boldsymbol{u}^{3,\xi})_{\xi\in\mathbb{R}^{M\times
N}}$, we obtain an orthonormal basis in the Fourier domain
$(\boldsymbol{v}^{1,\xi},\boldsymbol{v}^{2,\xi},\boldsymbol{v}^{3,\xi})_{\xi\in\mathbb{R}^{M\times
N}}$ of eigenvectors of ${\boldsymbol{\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{\boldsymbol{v}}}^{j_{1},\xi_{1}}\right)^{T}\widehat{\boldsymbol{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{\boldsymbol{v}}}^{j_{1},\xi_{1}}_{x_{1}}\widehat{\boldsymbol{v}}^{j_{2},\xi_{2}}_{x_{2}}$
(94) $\displaystyle=\mathbf{1}_{\begin{subarray}{c}j_{1}=j_{2}\\\
\xi_{1}=\xi_{2}\end{subarray}}$ (95)

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

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

In particular,

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

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

Let consider a Gaussian distribution
$\mathcal{N}({\boldsymbol{0}},\boldsymbol{\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}}$,

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

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

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

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

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

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

$\displaystyle\left(\overline{\widehat{\boldsymbol{v}}}^{j,\xi}\right)^{T}\widehat{\boldsymbol{Y}}$
$\displaystyle=\sqrt{\lambda_{1}^{\xi}}\left[\overline{\widehat{\boldsymbol{w}}}^{T}\widehat{\boldsymbol{v}}^{j,\xi}\right]\left\|\widehat{\boldsymbol{v}}^{j,\xi}\right\|^{2}=\sqrt{\lambda_{1}^{\xi}}\left[\overline{\widehat{\boldsymbol{w}}}^{T}\widehat{\boldsymbol{v}}^{j,\xi}\right]$
(102)
$\displaystyle\left|\left(\overline{\widehat{\boldsymbol{v}}}^{j,\xi}\right)^{T}\widehat{\boldsymbol{Y}}\right|^{2}$
$\displaystyle=\lambda_{j}^{\xi}\left|\overline{\widehat{\boldsymbol{w}}}^{T}\widehat{\boldsymbol{v}}^{j,\xi}\right|^{2}$
(103)
$\displaystyle\mathbb{E}\left[\left|\left(\overline{\widehat{\boldsymbol{v}}}^{j,\xi}\right)^{T}\widehat{\boldsymbol{Y}}\right|^{2}\right]$
$\displaystyle=\lambda_{j}^{\xi}\mathbb{E}\left[\left|\overline{\widehat{\boldsymbol{w}}}^{T}\widehat{\boldsymbol{v}}^{j,\xi}\right|^{2}\right]$
(104)
$\displaystyle\mathbb{E}\left[\left|\overline{\widehat{\boldsymbol{w}}}^{T}\widehat{\boldsymbol{v}}^{j,\xi}\right|^{2}\right]$
$\displaystyle=\sum_{1\leq c_{1},c_{2}\leq
3}\mathbb{E}\left[\overline{\widehat{\boldsymbol{w}}}_{c_{1}}(\xi)\widehat{\boldsymbol{w}}_{c_{2}}(\xi)\right]\widehat{\boldsymbol{v}}_{c_{1}}^{j,\xi}(\xi)\overline{\widehat{\boldsymbol{v}}}_{c_{2}}(\xi)\text{
by Equation\leavevmode\nobreak\ \eqref{eq:v_inner} }$ (105)
$\displaystyle=\sum_{1\leq c\leq
3}\mathbb{E}\left[\left|\widehat{\boldsymbol{w}}_{c}(\xi)\right|^{2}\right]\left|\widehat{\boldsymbol{v}}_{c}^{j,\xi}(\xi)\right|^{2}\text{
because the chanels are inependent}$ (106) $\displaystyle=3MN\sum_{1\leq c\leq
3}\left|\widehat{\boldsymbol{v}}_{c}^{j,\xi}(\xi)\right|^{2}\text{ because
$\mathbb{E}\left[\left|\widehat{\boldsymbol{w}}_{c}(\xi)\right|^{2}\right]=MN$}$
(107) $\displaystyle=3MN\text{ by Equation\leavevmode\nobreak\
\eqref{eq:v_norm}. }$ (108)

Finally,

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

Finally, for a given sampling $\left(\boldsymbol{Y}_{k}\right)_{1\leq k\leq
N_{\text{samples}}}$ following the distribution
$\mathcal{N}({\boldsymbol{0}},\boldsymbol{\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{\boldsymbol{v}}}^{j,\xi}\right)^{T}\widehat{\boldsymbol{Y}_{k}}\right|^{2}.$
(110)

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

$\operatorname{\text{$\textbf{W}^{\text{\tiny
emp.}}_{2}$}}(\mathcal{N}^{\text{emp.}}({\boldsymbol{0}},\boldsymbol{\Gamma}),\operatorname{ADSN}(\boldsymbol{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)}$
(111)

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.