<figure><img src="image_1.png" /></figure>

HAL
open science

# Dynamical regimes of diffusion models

Giulio Biroli, Tony Bonnaire, Valentin de Bortoli, Marc Mezard

► To cite this version:

Giulio Biroli, Tony Bonnaire, Valentin de Bortoli, Marc Mezard. Dynamical regimes of diffusion models. Nature Communications, 2024, 15 (1), pp.9957. 10.1038/s41467-024-54281-3. hal-04967666

HAL Id: hal-04967666

https://hal.science/hal-04967666v1

Submitted on 26 Feb 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Dynamical Regimes of Diffusion Models

Giulio Biroli*¹, Tony Bonnaire†¹, Valentin de Bortoli‡², and Marc Mézard§³

¹Laboratoire de Physique de l'Ecole Normale Supérieure, ENS, Université PSL, CNRS, Sorbonne Université, Université Paris Cité, F-75005 Paris, France

²Computer Science Department, ENS, CNRS, PSL University

³Department of Computing Sciences, Bocconi University, Milano, Italy

Abstract

Using statistical physics methods, we study generative diffusion models in the regime where the dimension of space and the number of data are large, and the score function has been trained optimally. Our analysis reveals three distinct dynamical regimes during the backward generative diffusion process. The generative dynamics, starting from pure noise, encounters first a 'speciation' transition where the gross structure of data is unraveled, through a mechanism similar to symmetry breaking in phase transitions. It is followed at later time by a 'collapse' transition where the trajectories of the dynamics become attracted to one of the memorized data points, through a mechanism which is similar to the condensation in a glass phase. For any dataset, the speciation time can be found from a spectral analysis of the correlation matrix, and the collapse time can be found from the estimation of an 'excess entropy' in the data. The dependence of the collapse time on the dimension and number of data provides a thorough characterization of the curse of dimensionality for diffusion models. Analytical solutions for simple models like high-dimensional Gaussian mixtures substantiate these findings and provide a theoretical framework, while extensions to more complex scenarios and numerical validations with real datasets confirm the theoretical predictions.

Keywords: Diffusion Models | Generative diffusion | Phase Transitions | Memorization | Collapse

Introduction

Machine learning has recently witnessed thrilling ad-
vancements, especially in the realm of generative mod-
els. At the forefront of this progress are diffusion models
(DMs), which have emerged as powerful tools for mod-
eling complex data distributions and generating new
realistic samples. They have become the state of the
art in generating images, videos, audio or 3d scenes
[43, 45, 46, 22, 48, 42, 3, 37]. Although the practical
success of generative diffusion models is widely rec-
ognized, their full theoretical understanding remains an
open challenge. Rigorous results assessing their conver-
gence on finite dimensional data have been obtained in
[18, 30, 17, 5, 15, 13]. However, realistic data live in high
dimensional spaces, where interpolation between data-
points should face the curse of dimensionality [20]. A
thorough understanding of how generative models es-
cape this curse is still lacking. This requires approaches
able to take into account that the number and the di-
mension of the data are simultaneously very large. In
this work we face this challenge using statistical physics
methods which have been developed to study proba-
bility distributions, disordered systems and stochastic
processes in high-dimensions [33, 12, 8].

Diffusion models work in two stages. In the forward
diffusion one starts from a data point (e.g., an image)
and gradually adds noise to it, until the image has be-

come pure noise. In the backward process one gradually
denoises the image using a diffusion in an effective force
field (the "score") which is learnt leveraging techniques
from score matching [26, 47] and deep neural networks.
In this study, we focus on diffusion models which are
efficient enough to learn the exact empirical score, i.e.
the one obtained by noising the empirical distribution of
data. In practical implementations, this should happen
when one performs a long training of a strongly over-
parameterized deep network to learn the score, in the
situation when the number of data is not too large.

Within this framework we develop a theoretical ap-
proach able to characterize the dynamics of diffusion
models in the simultaneous limit of large dimensions
and large dataset. We show that the backward generative
diffusion process consists of three subsequent dynamical
regimes. The first one is basically pure Brownian motion.
In the second one the backward trajectory finds one of
the main classes of the data (for instance if the data con-
sists of images of horses and cars, a given trajectory will
specialize towards one of these two categories). In the
third regime, the diffusion ‘collapses’ onto one of the ex-
amples of the dataset: a given trajectory commits to the
attraction basin of one of the data points, and the back-
ward evolution brings it back to that exact data point.
Since diffusion models are defined as the time reversal
of a forward noising process, the generative process has
to collapse on the training dataset under the exact empir-
ical score assumption [9, 23]. We show, by performing
a thorough analysis of the curse of dimensionality for
diffusion models, that this memorization can be avoided

<footnote>* giulio.biroli@ens.fr</footnote>

<footnote>†tony.bonnaire@ens.fr</footnote>

<footnote>‡vdebortoli@google.com</footnote>

<footnote>§marc.mezard@unibocconi.it</footnote>
at finite times only if the size of the dataset is exponentially large in the dimension. An alternative, which is the one used in practice, is to rely on regularization and approximate learning of the score, departing from its exact form. Understanding this crucial aspect of generative diffusion is a key open challenge [27, 49, 16] for which analyzing what happens when the exact empirical score is used represents a first step.

Separating these three dynamical regimes, we identify two characteristic cross-over times. The speciation time $t_S$ is the transition from pure noise to the commitment of a trajectory towards one category. The collapse time $t_C$ is the time where the backward trajectory falls into the attractor of one given data point. We provide mathematical tools to predict these times in terms of structure of data. We will first study simple models such as high-dimensional Gaussian mixtures, where we obtain a full analytical solution and hence a very detailed understanding. Within this setting, we show that in the simultaneous limit of large number and dimension of the data, the speciation and collapse cross-overs become sharp thresholds. Interestingly, both of them are related to phase transitions studied in physics. We then extend our results to more general settings and discuss the key role played by the dimensionality of data and the number of samples. Finally, we perform numerical experiments, and confront the theory to real data such as CIFAR-10, ImageNet, and LSUN, showing that our main findings hold in realistic cases. We conclude highlighting the consequences and the guidelines offered by our results, and discussing the next research steps, in particular how to go beyond the “exact empirical score” framework.

## Main Results and Setting

We focus on cases in which the data can be organized in distinct classes. For simplicity, we consider below two classes, identified in the spectrum of the covariance matrix of the data: this spectrum is assumed to display a single large eigenvalue along the principal component which we will denote $\Lambda$. This is a simplifying assumption; our analysis can be extended to more than two classes, and subclasses within classes. The data consists of $n$ data points $\vec{a} \in \mathbb{R}^d$. We assume that there exists a “true” underlying distribution $P_0(\vec{a})$ from which data are drawn, and we denote by $P_0^e(\vec{a}) = \sum_{\mu=1}^n \delta(\vec{a} - \vec{a}_\mu)/n$, the empirical distribution of the data. The components of $\vec{a}$ are normalized to be finite<fnref n="1" /> for large $d$. This implies in particular that the expectation of $|\vec{x}|^2$ grows linearly with $d$.

There exist many variants of diffusion models which are basically equivalent. We focus here on the diffusion process which consists in $d$ independent Ornstein-

Uhlenbeck Langevin equations,

$$d\vec{x}(t) = -\vec{x}(t)dt + d\vec{B}(t), \quad (1)$$

where $d\vec{B}(t)$ is square root of two times the standard Wiener process (a.k.a., Brownian motion) in $\mathbb{R}^d$. The “exact empirical score” is given by $F_i(\vec{x}, t) = \partial \log P_t^e(\vec{x})/\partial x_i$ where $P_t^e(\vec{x})$ is the noisy empirical distribution at time $t$ due to the process in (1)

$$P_t^e(\vec{x}) = \int d\vec{a} \, P_0^e(\vec{a}) \frac{1}{\sqrt{2\pi\Delta_t^d}} \exp\left(-\frac{1}{2} \frac{(\vec{x} - \vec{a}e^{-t})^2}{\Delta_t}\right). \quad (2)$$

This is the convolution of the empirical distribution of the data, $P_0^e(\vec{a})$, with a Gaussian law of variance $\Delta_t = (1 - e^{-2t})$. At long times $P_t^e(\vec{x})$ is a Gaussian distribution with zero mean and covariance equal to the identity. In DMs the generation of new data is obtained by time-reversing this process using the backward dynamics

$$-dy_i(t) = y_i dt + 2F_i(y, t)dt + d\xi_i(t), \quad (3)$$

where the noise $d\xi_i(t)$ has the same distribution as in the forward process.

Our main contribution is the characterization of three dynamical regimes in the time-reversed process in the limit of large number of data and large dimension. In regime I, at the beginning of the backward process, the random dynamical trajectories generated by (3) have not committed to a particular class of data. They have roughly the same probability to end up in one of the two classes. Figure 1 illustrates this behavior by showing that trajectories corresponding to different final classes are within a single bundle (or tube). In regime II, instead, the dynamical trajectories have committed to a particular class. In this case, the trajectory will remain in the same class until the end of the backward process. During regime II, the backward dynamics generates the features needed to produce samples in a given class, but the fate of the trajectory in terms of class is sealed. In analogy with evolutionary dynamics, we call “speciation” the cross-over between regimes I and II, which has also important connections with the concept of symmetry breaking in physics [7, 39]. As we shall show, the speciation cross-over takes place on a time-scale $t_S$ defined by

$$\Lambda e^{-2t_S} = 1, \quad \text{Speciation: I} \to \text{II} \quad (4)$$

where $\Lambda$ is the eigenvalue of the principal component of the covariance matrix of the data. Note that in the high-dimensional limit, if $\Lambda$ diverges with $d$ (typically one would expect $\Lambda \propto d$), then $t_S$ diverges logarithmically with $d$, and the speciation cross-over becomes a phase transition over time-scales of order $t_S^2$. We find that in the regimes I and II the DM generalizes (if the number

<footnote n="1">¹More precisely, we assume that the moments $\int d\vec{a} P_0(\vec{a}) a_i^p$ remain of order one for $d \to \infty$ and finite $p$.</footnote>

<footnote n="2">²Here, and in what follows, we measure time from the beginning of the forward process, i.e. large $t$ corresponds to the beginning of the backward process.</footnote>
<figure><img src="image_3.png" /><figcaption>**Figure 1:** Illustration of the three regimes of the backward dynamics through an example corresponding to a Gaussian mixture in two dimensions. Trajectories are colored white and blue according to their class at the end of the backward dynamics. In regime I, blue and white trajectories are fluctuating within the same bundle and $\vec{x}$ is similar to white noise. At the speciation time $t_S$, the ensembles of blue and white trajectories divide and head towards the distribution associated to their respective class. Regime II is where the generative process constructs an $\vec{x}$ which resembles to one element of the class (e.g., a seashore in the illustration) without being linked to any data of the training set. At the collapse time $t_C$, trajectories start to be attracted by the data point on which they collapse at $t=0$. Regime III corresponds to memorization, whereas in regime I and II, the diffusion model truly generalizes. The images on the right and on the left are illustrations obtained from our ImageNet numerical experiment (notice the collapse on the panda and seashore from the training set at $t=0$).</figcaption></figure>

of data is large enough), which means that the empirical distribution at time $t$, $P_t^e(\vec{x})$ is basically the same as the true $P_t(\vec{x})$, which is the convolution of $P_0$ and a Gaussian of variance $\Delta_t$. Therefore the distributions obtained by noising the desired $P_0$ and the empirical $P_0^e$ are identical. In regime III, the situation is completely different. There is no generalization; instead, the DM displays memorization of the training set. The probability distribution $P_t^e(\vec{x})$ no longer reflects $P_t(\vec{x})$. It actually decomposes in separated lumps around the points of the training set, and a given trajectory is committed to the attractor of the original data point, that is reached at $t=0$. Figure 1 illustrates this behavior. We call "collapse" the crossover between regimes II and III. As we shall show, the collapse takes place on a time-scale $t_C$ defined by

$$s(t_C) = s^{sep}(t_C), \quad \text{Collapse: II} \to \text{III} \quad (5)$$

where $s(t) = -\frac{1}{d} \int d\vec{x} P_t(\vec{x}) \log P_t(\vec{x})$ is the Shannon entropy per variable of the distribution at time $t$, and $s^{sep}(t) = \frac{\log n}{d} + \frac{1}{2} + \frac{1}{2}\log(2\pi\Delta_t)$ is its counterpart for a mixture of $n$ Gaussian distributions with variance $\Delta_t$ and centered on well separated points. We shall argue

below that the criterion (5) is actually valid beyond the "exact empirical score" hypothesis, and therefore provides a way to characterize the collapse (or its absence) in practical applications.

For $t \gg 1$, the distribution $P_t(\vec{x})$ becomes a $d$-dimensional Gaussian, with entropy $\frac{d}{2} + \frac{d}{2}\log(2\pi)$. In consequence, the difference between the two entropies, $f(t) = s^{sep}(t) - s(t)$, that we call "excess entropy density", equals $\frac{\log n}{d}$. On the other hand, for $t \to 0$, the entropy $s(t)$ goes to the one of $P_0(\vec{x})$, while $\log(2\pi\Delta_t) \to -\infty$. Therefore the excess entropy density goes to $-\infty$. When running the backward process, starting from large times where $f(t) = \frac{\log n}{d} > 0$, this excess entropy density decreases and crosses 0 at the collapse time $t_C$. Thus, by monitoring the time-dependence of $f(t)$, one can pinpoint when the collapse takes place. In numerical studies where $P_0$ is not known, one can approximate $f(t)$ as $f^e(t) = s^{sep}(t) - s^e(t)$, where $s^e(t)$ is the entropy per variable of the empirical distribution $P_t^e(\vec{x})$. This is a good approximation for $t \ge t_C$, where the two distributions $P_t^e$ and $P_t$ are identical and $f^e(t) = f(t)$. Instead, as we shall show, for $t < t_C$ the empirical excess entropy density vanishes while $f(t) < 0$. The parameter $\alpha = \frac{\log n}{d}$ plays a key role in the collapse, as it allows to tune the value of $t_C$. A very small $\alpha$, implies a very small excess entropy density for $t \gg 1$. In this case the collapse takes place at the very beginning of the backward dynamics. One needs $\alpha \sim O(1)$ in order to have $t_C \sim O(1)$. To diminish $t_C$ (shrink regime III) and reduce the collapse, one has to increase $\alpha$. These findings characterize a curse of dimensionality which is different but related to the one arising in supervised learning: in order to avoid memorization, the number of data has to increase exponentially with the dimension $d$. In practice, this unwanted phenomenon is avoided by using an approximate score function which is smoother than the exact one, together with a large enough dataset. We will come back to this point later in the discussion.

Our characterization of the backward dynamics is obtained using statistical physics methods developed to study phase transitions, disordered systems, and out-of-equilibrium dynamics in physics. We provide a brief introduction in the SI Appendix. The connection is not only methodological; in fact, our analysis unveils interesting relationships between the transition described above and phenomena intensively studied in physics in the last decades. In particular, the collapse transition is mutatis mutandis a glass transition in which the low energy glass states correspond to the training data. The results we obtain in Eqs. (4,5), which are at the level of rigor of theoretical physics, provide guidelines and testable predictions for realistic applications. Numerical experiments on several subsets of realistic datasets (MNIST, CIFAR-10, ImageNet, and LSUN) confirm their validity.
Analysis of Gaussian mixtures

An instructive simple example to study the backward dynamics in the large *d* and *n* limit is the Gaussian mixture model. In this case, the initial law $P_0(\vec{a})$ is the superposition of two Gaussian clusters of equal weight, which we take for simplicity with means $\pm \vec{m}$, and the same variance $\sigma^2$. We shall assume that $|\vec{m}|^2 = d\tilde{\mu}^2$, with $\sigma$ and $\tilde{\mu}$ of order 1.

**Speciation time.** In order to show the existence of regimes I and II, and compute the speciation time, we focus on the following protocol which consists in cloning trajectories. We consider a backward trajectory starting at time $t_f \gg 1$ from a point $\vec{x}_f$ drawn from a random Gaussian distribution where all components are independent with mean zero and unit variance. This trajectory evolves backward in time, through the backward process until time $t < t_f$. At this time the trajectory has reached the point $\vec{y}(t)$, at which cloning takes place. One generates for $\tau < t$ two clones, starting from the same $\vec{x}_1(t) = \vec{x}_2(t) = \vec{y}(t)$, and evolving as independent trajectories $\vec{x}_1(\tau)$ and $\vec{x}_2(\tau)$, i.e. with independent thermal noises. We compute the probability $\phi(t)$ that the two trajectories ending in $\vec{x}_1(0)$ and $\vec{x}_2(0)$ are in the same class. Defining $P(\vec{x}_1, 0|\vec{y}, t)$ as the probability that the backward process ends in $\vec{x}_1$, given that it was in $\vec{y}$ at time $t$, the joint probability of finding the trajectory in $\vec{y}$ at time $t$ and the two clones in $\vec{x}_1$ and $\vec{x}_2$ at time 0 is obtained as $Q(\vec{x}_1, \vec{x}_2, \vec{y}, t) = P(\vec{x}_1, 0|\vec{y}, t)P(\vec{x}_2, 0|\vec{y}, t)P(\vec{y}, t)$. Then $\phi(t)$ is the integral of $Q$ over $\vec{x}_1, \vec{x}_2, \vec{y}$ with the constraint $(\vec{x}_1 \cdot \vec{m})(\vec{x}_2 \cdot \vec{m}) > 0$. This simplifies into a one-dimensional integral (see SI Appendix):

$$ \phi(t) = \frac{1}{2} \int_{-\infty}^{+\infty} dy \frac{G(y, me^{-t}, \Gamma_t)^2 + G(y, -me^{-t}, \Gamma_t)^2}{G(y, me^{-t}, \Gamma_t) + G(y, -me^{-t}, \Gamma_t)}, \quad (6) $$

where $G(y, u, v)$ is a Gaussian probability density function for the real variable $y$, with mean $u$ and variance $v$, and $m = |\vec{m}| = \tilde{\mu}\sqrt{d}$, $\Gamma_t = \Delta_t + \sigma^2 e^{-2t}$. The probability $\phi(t)$ that the two clones end up in the same cluster is a decreasing function of $t$, going from $\phi(0) = 1$ to $\phi(\infty) = 1/2$. In the large $d$ limit, the scaling variable controlling the change of $\phi$ is $\tilde{\mu}\sqrt{d}e^{-t}$ which can be rewritten as $\tilde{\mu}e^{-(t-t_S)}$ by using $t_S = (1/2) \log d$. This explicitly shows that speciation takes place at the timescale $t_S$ on a window of time of order one. As detailed in the SI Appendix, this expression for $t_S$ coincides with the one obtained from the general criterion (5). We show in Fig. 2 the analytical result from (6) and direct numerical results obtained for increasingly larger dimensions. This comparison shows that our analysis is accurate already for moderately large dimensions. In the limit of infinite $d$, the analytical curve in Fig. 2 suddenly jumps from one to zero at $t/t_S = 1$, corresponding to a symmetry breaking phase transition (or a threshold phenomenon) on the time-scale $t_S$. In the numerics, following finite size scaling theory [38], we have defined the speciation

<figure><img src="image_4.png" /><figcaption>**Figure 2:** Speciation in Gaussian mixtures: Evolution of $\phi(t)$ as a function of $t/t_S$ for several values of $d$ at fixed $\tilde{\mu} = 1$ and $\sigma = 1$. The solid line corresponds to the evaluation of (6) while the dots are obtained by sampling 10 000 clone trajectories. The vertical (resp. horizontal) dashed line corresponds to $t/t_S = 1$ (resp. $\phi(t) = 0.775$). Error bars correspond to thrice standard error.</figcaption></figure>

time as the crossing point of the curves for different $d$, which corresponds approximatively to $\phi(t_S) = 0.775$ and indeed is of the order $t_S = (1/2) \log d$ for $d \to \infty$. As it happens in mean-field theories of phase transitions [35], the large dimensional limit allows to obtain a useful limiting process. In our case, this leads to a full characterization of the asymptotic backward dynamics. At its beginning, i.e. in regime I, the overlap with the centers of the Gaussian model, $\pm \vec{m} \cdot \vec{x}(t)$, is of order $\sqrt{d}$. Defining $q(t) = \vec{m} \cdot \vec{x}(t)/\sqrt{d}$, one can obtain a closed stochastic Langevin equation on $q$ in an inverted potential $V(q, t)$ (see SI Appendix),

$$ -dq = -\frac{\partial V(q, t)}{\partial q} dt + d\eta(t), \quad (7) $$

where $\eta(t)$ is square root of two times a Brownian motion, and

$$ V(q, t) = \frac{1}{2} q^2 - 2\tilde{\mu}^2 \log \cosh(qe^{-t}\sqrt{d}). \quad (8) $$

At large $d$, this potential is quadratic at times $t \gg t_S = (1/2) \log d$, and it develops a double well structure, with a very large barrier, when $t \ll t_S = (1/2) \log d$. The trajectories of $q(t)$ are subjected to a force that drives them toward plus and minus infinity. The barrier between positive and negative values of $q$ becomes so large that trajectories commit to a definite sign of $q$: this is how the symmetry breaking takes place dynamically at the timescale $t_S$, in agreement with the previous cloning results. Regime II corresponds to the scaling limit $q \to \infty$,
where $\vec{m} \cdot \vec{x}(t)$ becomes of order $d$. In this regime the rescaled overlap $\vec{m} \cdot \vec{x}(t)/d$ concentrates, and its sign depends on the set of trajectories one is focusing on. Moreover, the stochastic dynamics of the $x_i$ corresponds to the backward dynamics for a single Gaussian centered in $\pm \vec{m}$. This shows that the dynamics generalizes, see SI Appendix (and also [21] for similar results).

**Collapse time.** In order to study the collapse, we consider the probability distribution $P_t^e(\vec{x})$ given in (2) around a point $\vec{x}$ which has been obtained in the forward process starting at $t=0$ from $\vec{x} = \vec{a}_1$. Our aim is to establish whether the score obtained from $P_t^e(\vec{x})$ imposes a force that pushes trajectories toward $\vec{a}_1$ in the backward process, corresponding to memorization, or instead allows for generalization. We shall consider that both $n$ and $d$ go to infinity, keeping $\alpha = \log n/d$ fixed (it will be clear from the analysis that this is the correct ratio for the asymptotic analysis).

The vector $\vec{x}$ we consider is equal to $\vec{a}_1 e^{-t} + \vec{z} \sqrt{\Delta_t}$ where $\vec{z}$ has iid Gaussian components with zero mean and unit variance. The probability can be written as $P_t^e(\vec{x}) = [Z_1 + Z_{2...n}] / \sqrt{2\pi\Delta_t}^d$ where $Z_1 = e^{-\frac{1}{2}(\vec{x}-\vec{a}_1e^{-t})^2/(2\Delta_t)} = e^{-\frac{\vec{z}^2}{2}}$ and

$$Z_{2...n} = \sum_{\mu=2}^{n} e^{-\frac{1}{2}[(\vec{x}-\vec{a}_{\mu}e^{-t})^2/(2\Delta_t),} \quad (9)$$

In the large $d$ limit (to exponential accuracy), $Z_1 \approx e^{-d/2}$. The computation of $Z_{2...n}$ is instead tricky. Even though it is a sum of $n$ uncorrelated random contributions, standard concentration methods, e.g. central limit theorem, do not apply as each term of the sum corresponds to the exponential of a random variable scaling as $\log n$ [4]. Statistical physics tools developed to study spin-glasses provide the method to solve this problem, see SI Appendix. In fact, given $\vec{x}$, $Z_{2...n}$ is the partition function of a system with $n-1$ independent 'random energy levels' $E^{\mu} = \frac{1}{2}[(\vec{a}_1 - \vec{a}_{\mu})e^{-t} + \vec{z}\sqrt{\Delta_t}]^2/(2\Delta_t)$ at thermal equilibrium at a temperature 1 (the randomness comes from $\vec{a}_{\mu}$). This is some elaboration of the 'Random Energy Model' which was introduced originally in [19] as a simple model of glass transition. A similar problem [31] was studied recently in the related context of dense associative memories, using large-deviations and replica methods (this connection was also noticed in [2]). Using a similar approach we can show (see SI Appendix) that in the large $d$ and $n$ limit the distribution of $(1/d) \log Z_{2...n}$ concentrates on a value $\psi_+(t)$ which is an increasing function of the time $t$ and of $\alpha = \lim_{d,n \to \infty} (1/d) \log n$. The analytical computation shows the existence of a collapse time $t_C$ which separates two time regimes:

* Regime III: at small times, $t < t_C$, $\psi_+(t) < -(1/2)$ and the probability $P_t^e(\vec{x})$ is dominated by the term $Z_1 = e^{-(\vec{x}-\vec{a}_1e^{-t})^2/(2\Delta_t)}$. When used in the backward diffusion, this gives a score that attracts $\vec{x}$ towards $\vec{a}_1$ at short times. With probability one, the backward trajectory, starting at time $t$ from

$\vec{x} = \vec{a}_1 e^{-t} + \vec{z} \sqrt{\Delta_t}$, collapses at the end of the backward process on data point $\vec{a}_1$. In this regime the associated Random Energy Model is in a glass phase, which precisely corresponds to memorization.

* Regimes I and II: at large times, $t > t_C$, $\psi_+(t) > -(1/2)$ and the probability $P_t^e(\vec{x})$ is dominated by the term $Z_{2...n}$. This is the regime which is not collapsed and corresponds to generalization (regime II) or Brownian motion (regime I): in a typical point $\vec{x}(t)$, drawn from the population distribution $P_t(\vec{x})$, the empirical distribution $P_t^e(\vec{x}, t)$ is equal to $P_t(\vec{x}, t)$ at leading order in $d$. In this regime the associated Random Energy Model is in the liquid (or high-temperature) phase.

The collapse time reads

$$t_C = \frac{1}{2} \log \left( 1 + \frac{\sigma^2}{n^{2/d} - 1} \right). \quad (10)$$

This equation makes manifest the curse of dimensionality: $t_C$ is of order one only when the number of data is exponential in the dimension. One needs $\log n/d \gg 1$ to push $t_C$ to zero and avoid the collapse. One can check that this equation for $t_C$ coincides with the one obtained from the general criterion (5). In Fig. 3 we plot the empirical excess entropy density $f^e(t)$ from which one deduces $t_C$ as the largest time at which $f^e(t) = 0$. The numerical results compare well with the analytical prediction, and confirms that the time-scale for collapse is well captured by our approach, even at moderately large values of $n$. When $n, d \to \infty$ at fixed $\alpha = (1/d) \log n$, our analysis shows that regime III takes place on timescale of order one, and hence after the speciation transition, as illustrated in Fig. 1.

## Generalization to realistic datasets

In the case of Gaussian mixtures we could thoroughly characterize the backward dynamics for large $n$ and $d$, using the knowledge of the initial distribution $P_0$. We now present a more general approach to the speciation and collapse transitions, for cases in which $P_0$ is not known or too complex to be analyzed exactly, but instead a dataset drawn from $P_0$ is available. In particular, we present arguments to establish criteria from Eqs. (4,5) which can be directly applied to realistic datasets.

The speciation transition can be analyzed using the covariance matrix, $C_0$, of the initial data. Our assumption that data can be organized in two distinct and very different classes translates, in terms of $C_0$, in the existence of a strong principal component with eigenvalue $\Lambda$. The speciation transition can be generically understood in terms of the forward process: it corresponds to the time-scale at which the noise added to the data blurs the principal component, and hence the connection to a given class. On this timescale, the trajectories coming from different classes in the forward process coalesce within the same
<figure><img src="image_5.png" /><figcaption>**Figure 3:** Collapse in Gaussian mixtures: Evolution of the excess entropy density $f^e(t)/\alpha$ as a function of time $t$ for several values of $d$, at fixed $n = 20 000$. The solid lines are the theoretical predictions while the dots show the results of the numerical evaluation approximating the entropy from the dataset. The vertical dashed lines represent the collapse time $t_C$ predicted analytically for Gaussian mixtures given in (10). Error bars correspond to thrice the standard error.</figcaption></figure>

bundle, as illustrated in Fig. 1. By time-reversal symmetry, it is therefore also on this time that the trajectories in the backward process separate and commit to the different classes. By evaluating the covariance matrix of the noised data $x_i(t)$, one finds

$$C(t) = C_0e^{-2t} + \Delta_t \mathbf{I}. \qquad (11)$$

The speciation time can be found comparing the two contributions on the RHS in the direction of the principal component of $C_0$: $\Lambda e^{-2t}$ vs $\Delta_t$. The first one is associated to the fluctuations between different classes, whereas the second corresponds to the broadening of $P_t(\vec{x})$ due to the noise. When the latter becomes of the same order of the former, the noise blurs the identification of trajectories in different classes. In consequence, one finds the criterion $\Lambda e^{-2t} \sim \Delta_t$, which shows that the variable controlling the speciation transition is indeed $\Lambda e^{-2t}$, as we have shown for the Gaussian mixture model. Since for large $d$ one expects a large $\Lambda$, the asymptotic window over which the speciation transition takes place is when $\Lambda e^{-2t}$ is of order one. This leads to the general result stated in (4)$^3$. The argument above is obtained studying the time at which the structure of the data, present at $t = 0$, is blurred by the noise. One can also tackle the problem starting by the other end, i.e. large times, and studying $P_t(\vec{x})$ perturbatively in $e^{-t}$, which is a small parameter

at the beginning of the backward process. By performing an expansion of (2) in $e^{-t}$, one finds that at leading order the distribution $P_t(\vec{x})$ is a multivariate Gaussian with covariance $\frac{1}{\Delta_t} [\mathbf{I} - \frac{e^{-2t}}{\Delta_t} C_0]$. The expansion described above, and detailed in the SI, is similar to the Landau expansion of phase transitions [10], in which symmetry breaking can be understood in terms of the instability of the quadratic part of $\log P_t(\vec{x})$. The speciation transition is like a symmetry breaking phase transition in $P_t(\vec{x})$ [7, 39], and this instability happens for $\Lambda e^{-2t} = \Delta_t \simeq 1$. The general result (4) can thus be obtained from two different complementary perspectives. We expect it to hold in a broad set of cases.

We now turn to the collapse transition. Generalization can be understood as the regime in which, at leading order in $d$, $P_t^e(\vec{x})$ coincides with its population counterpart $P_t(\vec{x})$ on typical noised data $\vec{x}(t)$ (i.e., drawn from $P_t(\vec{x})$). This implies that at time $t$ the training set has been forgotten. At large $d$, and for a given data $\vec{a}$, the random vectors $\vec{a}e^{-t} + \vec{z}\sqrt{\Delta_t}$, drawn from $P_t^e(\vec{x})$, lie with probability one on the ball$^4$ of radius $\sqrt{d\Delta_t}$ around $\vec{a}e^{-t}$. At small time $t$, one can therefore envision the set $M^e$ typically covered by the empirical distribution $P_t^e$ as the union of the non-overlapping balls centered around the data vectors $\vec{a}_\mu e^{-t}$. In this regime, the set $M$ typically covered by the population distribution $P_t(\vec{x})$ is clearly different. It is independent of the training set and not as singular for $t = 0$. In consequence, in this regime $P_t^e(\vec{x})$ is collapsed on the data and the diffusion model is in the memorization phase (regime III). By increasing $t$, the balls progressively grow and at a certain time, $t_C$, they cover the set $M$. Beyond this point, one expects that the empirical and population distribution coincide on $M$ (at leading order in $d$). This picture suggests a volume argument to identify the collapse: one finds the time at which the volume of $M^e$ coincides with the one of $M$. The key point is that in large $d$, the volume $V_P$ covered by typical configurations associated to a given distribution $P$ scales as $V_P = e^S$ where $S = - \int d\vec{x}P(\vec{x}) \log P(\vec{x})$ is the Shannon entropy, which scales like $d$. The set $M^e$ corresponds to $n$ d-dimensional Gaussian distributions with mean zero and covariance matrix given by the identity times $\Delta_t$. For $t \le t_C$, the volume of $M^e$ is the one of $n$ non-overlapping Gaussians, which therefore reads: $V_{M^e} = ne^{S_G}$ where $S_G = (d/2)(1 + \log(2\pi\Delta_t))$ is the entropy of one of the Gaussian distributions. On the other hand, the volume of $M$ is given by $e^{ds(t)}$. By requiring the equality of the two volumes, $V_{M^e} = V_M$, one finds the general result for $t_C$ given in (5)$^5$. This criterion can also be obtained from the Random Energy Model method, generalizing our derivation for the Gaussian Mixture model. Note that the arguments above can also be applied to cases in which the score is learned

<footnote>$^3$We have decided to associate the speciation transition with the time at which the variable $\Lambda e^{-2t}$ takes the value one, but this is a convention and the choice of another number of order one would work too when $\Lambda$ is large.</footnote>

<footnote>$^4$Actually they are located close to the circumference, as this is the region with the largest volume in large $d$.</footnote>

<footnote>$^5$Note that all the identities used in the volume argument are correct up to corrections exponentially small in $d$. For this reason, (5) is expressed in terms of intensive quantities, and hence valid up to vanishing corrections for $d \to \infty$.</footnote>
<table><caption>**Table 1:** Datasets for speciation experiments</caption><thead><tr><th>Dataset name</th><th>n</th><th>d</th><th>Λ</th><th>t<sub>S</sub></th></tr></thead><tbody><tr><td>1. MNIST</td><td>10 000</td><td>1024</td><td>7.66</td><td>1.02</td></tr><tr><td>2. CIFAR</td><td>3000</td><td>3072</td><td>16.72</td><td>1.41</td></tr><tr><td>3. ImageNet16</td><td>2000</td><td>768</td><td>3.05</td><td>0.56</td></tr><tr><td>4. ImageNet32</td><td>2000</td><td>3072</td><td>12.11</td><td>1.25</td></tr><tr><td>5. LSUN</td><td>40 000</td><td>12 288</td><td>60.52</td><td>2.05</td></tr></tbody></table>

by approximate models. In this case the distribution, $P_t(\vec{x})$ and its associated entropy density $s(t)$, have to be replaced by their model-dependent counterparts. The arguments above offer another way, besides the exact computation by the Random Energy Model method, to illustrate the relationship between the collapse and the glass transition studied in physics [6]. In fact, the ways in which $P_t^e(\vec{x})$ covers the space of $\vec{x}$ before and after the collapse is the exact counterpart of how the Boltzmann law covers the configuration space before and after the ideal glass transition [6]. In the glass phase the Boltzmann law is formed by lumps centered around amorphous optimal configurations (ideal glasses). Whereas in the liquid phase the Boltzmann law is spread over all configurations. In the case studied here, the elements of the training set are the counterparts of the ideal glass configurations.

Having shown the generality of our criteria for speciation and collapse, we now test them on realistic data.

## Numerical experiments

In the following we focus on various image datasets and we learn the score function using a state-of-the-art neural network with a finite number of samples $n$. We assume that, by using a heavily over-parameterized model, the resulting estimate $\hat{\mathcal{F}}$ of the score is close to the exact one.

**Settings and datasets.** We train a diffusion model as first described in [25], and corresponding to (1), with a time horizon $T = 5.05$ and a linear schedule for the variance. The denoiser is a U-Net [40] with an architecture similar to other traditional implementations of DMs [25, 44, 39]. Our training sets are constructed by focusing on data divided in two classes and extracted from the image datasets: MNIST [29], CIFAR-10 [28], downsampled ImageNet [14], and LSUN [50]. For example, in the case of ImageNet, the two classes we use are panda and seashore (as illustrated in Fig. 1). The large variety of datasets allows us to explore different kinds of images and several values of $n$ and $d$, see Table 1. We refer to the SI Appendix for details on the processing of these datasets, the denoiser, and how the numerical experiments presented below are carried out.

**Speciation time.** To numerically extract the speciation time and compare it to the theoretical result given by (4), we use the same cloning procedure introduced for Gaussian mixtures. We generate two independent

<figure><img src="image_6.png" /><figcaption>**Figure 4:** Speciation in realistic datasets: Evolution of $\phi(t)$, the probability that the two clones end up in the same class, as a function of $t/t_S$ for several image datasets. The values of $t_S$ are the theoretical prediction for the speciation time obtained using (4) and listed in Table 1. The dashed horizontal line indicate $\phi(t) = 0.775$ and the errorbars correspond to thrice the standard error.</figcaption></figure>

trajectories at time $t$ of the backward process, and numerically estimate the probability $\phi(t)$ that they end up in the same class by averaging over many initial conditions and backward processes. In order to recognize the class of $\vec{x}_1(0)$ and $\vec{x}_2(0)$, we use a classifier with a ResNet-18 architecture [24], which has a test accuracy larger than 95% on the datasets we focus on. The corresponding results are shown in Fig. 4. The time axis has been rescaled by the theoretically predicted $t_S = (1/2) \log \Lambda$ from (4). The values of $\Lambda$ and $t_S$ are listed in Table 1. Figure 4 shows that indeed the speciation phenomenon is at play in realistic DMs: $\phi(t)$ goes from 0.5 when $t \gg t_S$ to one when $t \ll t_S$, in a way which is qualitatively analogous to Fig. 2. When rescaled with respect to $t_S$, the time-dependence exhibits remarkable similarity across vastly different image sets, suggesting evidence of the common underlying phenomenon of speciation. Moreover, our prediction for the speciation time, $t_S$, captures well the behavior found in these numerical experiments.

**Collapse time.** We focus on the datasets ImageNet16, ImageNet32, and LSUN. In order to be able to study thoroughly the collapse phenomenon, we have to keep the number of training data small ($n = 200$ for LSUN and $n = 2000$ for ImageNet). Otherwise, the model is not expressive enough to represent the singular behavior of the exact score at small times. To estimate the collapse time numerically, and compare it to the theoretical result from (5), we follow two distinct strategies. First, we use again the cloning procedure, but now we estimate the probability $\phi_C(t)$ that the two cloned trajectories collapse onto the same data of the training set at time zero.
The top left panel of Fig. 5 shows the evolution of $\phi_C(t)$. The cross-over time where $\phi_C(t)$ goes from zero to one provides a first numerical estimate of $t_C$. An alternative estimation is found by tracking, during the backward process, the index of the closest neighbor in the dataset, noted $\mu_*(t) \in \{1, \dots, n\}$. For a single trajectory, we estimate the collapse time $\hat{t}_C$ as the last time during the backward process at which this index changes, meaning that for $t < \hat{t}_C$, the nearest neighbor of $\vec{x}(t)$ is always the same. The top right panel of Fig. 5 represents the distribution of $\hat{t}_C$ measured in this way on the LSUN dataset for illustration. By computing the average of this distribution for each dataset, we obtain the estimates $\overline{\hat{t}}_C$ of the collapse time. We indicated them by vertical dashed lines in all panels of Fig. 5. The two different approaches to estimate the collapse time agree well, as seen in the top left panel where $\overline{\hat{t}}_C$ crosses $\phi_C(t)$ at around 0.60 for all datasets. We can now test the theoretical criterion for the collapse time given by (5), or equivalently by the largest time at which the empirical excess entropy density $f^e(t)$ is equal to zero, see the bottom panel of Fig. 5. The shape of $f^e(t)/\alpha$ is remarkably similar to the one obtained analytically for the Gaussian Model (Fig. 3). It clearly shows evidences of a transition. The theoretical estimate, $f^e(t_C) = 0$, also compares very well with the numerical ones $\overline{\hat{t}}_C$.

In conclusion, the numerical experiments presented in this section demonstrate the presence of speciation and collapse in realistic image datasets, and validate our theory, in particular the criteria from Eqs. (4) and (5) to identify the times at which these phenomena takes place.

<figure><img src="image_7.png" /><figcaption>**Figure 5:** Collapse in realistic datasets (ImageNet16, ImageNet32 and LSUN): (Top Left) Evolution of $\phi_C(t)$, the probability that two cloned trajectories collapse on the same data of the training set at time zero. (Top Right) Histograms of $\hat{t}_c$ derived from the last-changing indices $\mu_*$ on 4000 generated samples for the LSUN dataset trained with $n = 200$. (Bottom) Evolution of the empirical excess entropy $f^e(t)/\alpha$. In all panels, the colored vertical dashed lines indicate the average of $\overline{\hat{t}}_C$. The errorbars correspond to thrice the standard error.</figcaption></figure>

## Discussion

In this work we have analyzed the backward dynamics of diffusion models assumed to be efficient enough to learn the exact empirical score. We have shown that for large number of data $n$, large dimension $d$, and in absence of regularization, the backward dynamics displays three different regimes. We have characterized the cross-overs between them, dubbed *speciation* and *collapse*, which become true transitions in the large $n, d$ limit. Interestingly, both of them have physical counterparts in the theory of phase transitions. Speciation is a symmetry-breaking transition of $P_t(\vec{x})$ at which the most prominent classes are generated [39, 7]. Collapse corresponds to a glass transition at which $P_t(\vec{x})$ fragments in an ensemble of lumps centered around the training data. Our approach characterizes the time at which speciation and collapse take place in terms of structure of data. The speciation time is determined by the eigenvalue of the principal component of the covariance matrix of the data, whereas the collapse time is governed by the entropy of the noised data distribution.

Although we have focused on a specific diffusion model, our results hold for the large variety of DMs which are based on inverting the noising process [1]. These meth-

ods use different procedures to implement the denoising process $P_\infty(\vec{x}) \rightarrow P_t(\vec{x}) \rightarrow P_0(\vec{x})$. However, as we have shown, the phenomena of speciation and collapse, can be understood and analyzed focusing on the properties of the distribution $P_t(\vec{x})$ alone. In consequence, since all these different DMs lead to the same backward evolution of $P_t(\vec{x})$, the criteria from Eqs. (4) and (5) for speciation and collapse hold in general, whether the generative process is deterministic and flow-based or stochastic and diffusion-based.

Our work was done within the “exact empirical score” hypothesis. However, the phenomena we analyzed are relevant also when the score is learned approximately, as shown by our numerical experiments. Moving beyond the “exact empirical score” hypothesis opens up multiple avenues for further research. On the one hand, it would be interesting to analyze simple models of data and score in the limit of large number of data, dimension and parameters, as initiated in [7, 16]. On the other hand, it would be important to develop a quantitative study of the role of regularization on the phenomena we presented in this work. As discussed in [7], using a good model of the score is a key ingredient to generate data in the right proportions, i.e. reproducing the correct
weights of the classes. In consequence, regularization methods could be detrimental in regime I and for spe- ciation: by preventing the score to be close to the exact one they could lead to data with correct features (e.g., realistic images) but wrong proportions with respect to the training set. On the contrary, regularization is benefi- cial to avoid the collapse. As we have shown, diffusion models are “cursed”: one needs an exponential number of data in *d* to avoid the collapse on datapoints from the training set. Regularization allows to circumvent this problem in practice if the number of data is large enough [27] (see also SI Appendix). The volume argument we develop in this work can be applied beyond the “exact empirical score” hypothesis, and could offer a way to analyze quantitatively how the collapse depends on *n*, *d* and the capacity of the model used to learn the score.

Finally, our results also provide suggestions to im-
prove and understand procedures used in practical ap-
plications. In particular, taking into account in practical
implementations the existence of three regimes of the
backward dynamics could lead to better performances.

Acknowledgments

The authors thank F. Bach, S. Mallat, A. Montanari for
helpful discussions. GB acknowledges support from
the French government under the management of the
Agence Nationale de la Recherche as part of the “In-
vestissements d’avenir” program, reference ANR-19-
P3IA0001 (PRAIRIE 3IA Institute). MM acknowledges
the support of the PNRR-PE-AI FAIR project funded by
the NextGeneration EU program.

References

[1] Michael S Albergo, Nicholas M Boffi, and Eric Vanden-Eijnden. Stochastic interpolants: Auni-fying framework for flows and diffusions. *arXiv preprint arXiv:2303.08797*, 2023.

[2] Luca Ambrogioni. The statistical thermodynam-
ics of generative diffusion models. arXiv preprint
arXiv:2310.17467, 2023.

[3] Omer Bar-Tal, Hila Chefer, Omer Tov, Charles Herrmann, Roni Paiss, Shiran Zada, Ariel Ephrat, Junhwa Hur, Yuanzhen Li, Tomer Michaeli, et al. Lumiere: A space-time diffusion model for video generation. arXiv preprint arXiv:2401.12945, 2024.

[4] Gérard Ben Arous, Leonid V Bogachev, and Stanislav A Molchanov. Limit theorems for sums of random exponentials. *Probability theory and related fields*, 132:579–612, 2005.

[5] Joe Benton, Valentin De Bortoli, Arnaud Doucet, and George Deligiannidis. Linear convergence bounds for diffusion models via stochastic localization. *arXiv preprint arXiv:2308.03686*, 2023.

[6] Ludovic Berthier and Giulio Biroli. Theoretical perspective on the glass transition and amorphous materials. *Reviews of modern physics*, 83(2):587, 2011.

[7] Giulio Biroli and Marc Mézard. Generative diffusion in very large dimensions. *J. Stat. Mech.*, page 093402, 2023.

[8] Tony Bonnaire, Davide Ghio, Kamesh Krishna-murthy, Francesca Mignacco, Atsushi Yamamura, and Giulio Biroli. High-dimensional non-convex landscapes and gradient descent dynamics. arXiv preprint arXiv:2308.03754, 2023.

[9] Patrick Cattiaux, Giovanni Conforti, Ivan Gentil, and Christian Léonard. Time reversal of diffusion processes under a finite entropy condition. arXiv preprint arXiv:2104.07708, 2021.

[10] Paul M Chaikin, Tom C Lubensky, and Thomas A Witten. *Principles of condensed matter physics*, volume 10. Cambridge university press Cambridge, 1995.

[11] Paul M Chaikin, Tom C Lubensky, and Thomas A Witten. *Principles of condensed matter physics*, volume 10. Cambridge university press Cambridge, 1995.

[12] Patrick Charbonneau, Enzo Marinari, Giorgio Parisi, Federico Ricci-tersenghi, Gabriele Sicuro, Francesco Zamponi, and Marc Mezard. *Spin Glass Theory and Far Beyond: Replica Symmetry Breaking after 40 Years*. World Scientific, 2023.

[13] Sitan Chen, 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.

[14] Patryk Chrabaszcz, Ilya Loshchilov, and Frank Huter. A downsampled variant of imagenet as an alternative to the cifar datasets, 2017.

[15] Giovanni Conforti, Alain Durmus, and Marta Gentiloni Silveri. Score diffusion models without early stopping: finite fisher information is all you need. *arXiv preprint arXiv:2308.12240*, 2023.

[16] Hugo Cui, Florent Krzakala, Eric Vanden-Eijnden, and Lenka Zdeborová. Analysis of learning a flow-based generative model from limited sample complexity. arXiv preprint arXiv:2310.03575, 2023.

[17] Valentin De Bortoli. Convergence of denoising dif-
fusion models under the manifold hypothesis. arXiv
preprint arXiv:2208.05314, 2022.

[18] Valentin De Bortoli, James Thornton, Jeremy Heng, and Arnaud Doucet. Diffusion schrödinger bridge with applications to score-based generative modeling. *Advances in Neural Information Processing Systems*, 34:17695–17709, 2021.
[19] Bernard Derrida. Random-energy model: An exactly solvable model of disordered systems. *Physical Review B*, 24(5):2613, 1981.

[20] David L Donoho et al. High-dimensional data analysis: The curses and blessings of dimensionality. *AMS math challenges lecture*, 1(2000):32, 2000.

[21] Davide Ghio, Yatin Dandi, Florent Krzakala, and Lenka Zdeborová. Sampling with flows, diffusion and autoregressive neural networks: A spin-glass perspective. *arXiv preprint arXiv:2308.14085*, 2023.

[22] Florentin Guth, Simon Coste, Valentin De Bortoli, and Stephane Mallat. Wavelet score-based generative modeling, 2022.

[23] Ulrich G Haussmann and Etienne Pardoux. Time reversal of diffusions. *The Annals of Probability*, pages 1188–1205, 1986.

[24] Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun. Deep residual learning for image recognition. In *2016 IEEE Conference on Computer Vision and Pattern Recognition (CVPR)*, pages 770–778, 2016. doi: 10.1109/CVPR.2016.90.

[25] Jonathan Ho, Ajay Jain, and Pieter Abbeel. Denoising diffusion probabilistic models. *Advances in Neural Information Processing Systems*, 2020.

[26] Aapo Hyvärinen and Peter Dayan. Estimation of non-normalized statistical models by score matching. *Journal of Machine Learning Research*, 6(4), 2005.

[27] Zahra Kadkhodaie, Florentin Guth, Eero P Simoncelli, and Stéphane Mallat. Generalization in diffusion models arises from geometry-adaptive harmonic representation. *arXiv preprint arXiv:2310.02557*, 2023.

[28] Alex Krizhevsky, Geoffrey Hinton, et al. Learning multiple layers of features from tiny images. *Toronto, ON, Canada*, 2009.

[29] Y. Lecun, L. Bottou, Y. Bengio, and P. Haffner. Gradient-based learning applied to document recognition. *Proceedings of the IEEE*, 86(11):2278–2324, 1998. doi: 10.1109/5.726791.

[30] Holden Lee, Jianfeng Lu, and Yixin Tan. Convergence for score-based generative modeling with polynomial complexity. *arXiv preprint arXiv:2206.06227*, 2022.

[31] Carlo Lucibello and Marc Mézard. The exponential capacity of dense associative memories. *Phys. Rev. Lett*, 132:077301, 2024.

[32] Calvin Luo. Understanding diffusion models: A unified perspective. *arXiv preprint arXiv:2208.11970*, 2022.

[33] Marc Mezard and Andrea Montanari. *Information, physics, and computation*. Oxford University Press, 2009.

[34] Marc Mézard, Giorgio Parisi, and Miguel Angel Virasoro. *Spin glass theory and beyond: An Introduction to the Replica Method and Its Applications*, volume 9. World Scientific Publishing Company, 1987.

[35] Manfred Opper and David Saad. *Advanced mean field methods: Theory and practice*. MIT press, 2001.

[36] Adam Paszke, Sam Gross, Francisco Massa, Adam Lerer, James Bradbury, Gregory Chanan, Trevor Killeen, Zeming Lin, Natalia Gimelshein, Luca Antiga, Alban Desmaison, Andreas Kopf, Edward Yang, Zachary DeVito, Martin Raison, Alykhan Tejani, Sasank Chilamkurthy, Benoit Steiner, Lu Fang, Junjie Bai, and Soumith Chintala. Pytorch: An imperative style, high-performance deep learning library. In *Advances in Neural Information Processing Systems* 32, pages 8024–8035. Curran Associates, Inc., 2019.

[37] Ben Poole, Ajay Jain, Jonathan T Barron, and Ben Mildenhall. Dreamfusion: Text-to-3d using 2d diffusion. *arXiv preprint arXiv:2209.14988*, 2022.

[38] Vladimir Privman. *Finite size scaling and numerical simulation of statistical systems*. World Scientific, 1990.

[39] Gabriel Raya and Luca Ambrogioni. Spontaneous symmetry breaking in generative diffusion models. *arXiv preprint arXiv:2305.19693*, 2023.

[40] Olaf Ronneberger, Philipp Fischer, and Thomas Brox. U-net: Convolutional networks for biomedical image segmentation. In *Medical Image Computing and Computer-Assisted Intervention–MICCAI 2015: 18th International Conference, Munich, Germany, October 5-9, 2015, Proceedings, Part III 18*, pages 234–241. Springer, 2015.

[41] David Ruelle. A mathematical reformulation of derrida's rem and grem. *Communications in Mathematical Physics*, 108:225–239, 1987.

[42] Chitwan Saharia, William Chan, Saurabh Saxena, Lala Li, Jay Whang, Emily L Denton, Kamyar Ghasemipour, Raphael Gontijo Lopes, Burcu Karagol Ayan, Tim Salimans, et al. Photorealistic text-to-image diffusion models with deep language understanding. *Advances in Neural Information Processing Systems*, 35:36479–36494, 2022.

[43] Jascha Sohl-Dickstein, Eric Weiss, Niru Maheswaranathan, and Surya Ganguli. Deep unsupervised learning using nonequilibrium thermodynamics. In *International Conference on Machine Learning*, 2015.
[44] Jiaming Song, Chenlin Meng, and Stefano Ermon. Denoising diffusion implicit models. *arXiv preprint arXiv:2010.02502*, 2020.

[45] Yang Song and Stefano Ermon. Generative modeling by estimating gradients of the data distribution. *Advances in Neural Information Processing Systems*, 2019.

[46] 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.

[47] Pascal Vincent. A connection between score matching and denoising autoencoders. *Neural computation*, 23(7):1661–1674, 2011.

[48] Ling Yang, Zhilong Zhang, Yang Song, Shenda Hong, Runsheng Xu, Yue Zhao, Yingxia Shao, Wentao Zhang, Bin Cui, and Ming-Hsuan Yang. Diffusion models: A comprehensive survey of methods and applications. *arXiv preprint arXiv:2209.00796*, 2022.

[49] TaeHo Yoon, Joo Young Choi, Sehyun Kwon, and Ernest K Ryu. Diffusion probabilistic models generalize when they fail to memorize. In *ICML 2023 Workshop on Structured Probabilistic Inference {\&} Generative Modeling*, 2023.

[50] Fisher Yu, Ari Seff, Yinda Zhang, Shuran Song, Thomas Funkhouser, and Jianxiong Xiao. Lsun: Construction of a large-scale image dataset using deep learning with humans in the loop, 2016.

[51] Jean Zinn-Justin. *Quantum field theory and critical phenomena*, volume 171. Oxford university press, 2021.
# Supplementary Information (SI) Appendix
## Dynamical Regimes of Diffusion Models

Giulio Biroli, Tony Bonnaire, Valentin de Bortoli, and Marc Mézard

In this Appendix we first present more details on the analytical frameworks we used to study the backward dynamics of DMs, and then provide additional information on the numerical experiments. As already stressed in the main text, we need specific methods to study the limit of large dimension and large number of data. Such methods have been developed in statistical physics to study dynamics and thermodynamics of a large number of degrees of freedom. We will refer to original articles and, when needed, provide a short-introduction to the methods.

## A Landau-type expansion for estimating the speciation time

The speciation transition is similar to the symmetry breaking phenomenon [39, 7] taking place at thermodynamic phase transitions, for instance when a ferromagnetic system develops a non-zero magnetization at low temperature. One way to study this phenomenon is to construct perturbatively the free energy as a function of the field [11]. In our case, we proceed in a similar way by obtaining $\log P_t(\vec{x})$ in a perturbative expansion in $e^{-t}$ valid at large times. This approach is justified since speciation takes place at large times.

We rewrite the probability to be at $\vec{x}$ at time $t$ as

$$ 
\begin{aligned}
P_t(\vec{x}) &= \int d\vec{a} P_0(\vec{a}) \frac{1}{\sqrt{2\pi\Delta_t^d}} \exp \left( -\frac{1}{2} \frac{(\vec{x} - \vec{a}e^{-t})^2}{\Delta_t} \right) \\
&= \frac{1}{\sqrt{2\pi\Delta_t^d}} \exp \left( -\frac{1}{2} \frac{\vec{x}^2}{\Delta_t} + g(\vec{x}) \right)
\end{aligned} 
\quad (12) $$

where $g(\vec{x})$, defined as

$$ g(\vec{x}) = \log \int d\vec{a} P_0(\vec{a}) \exp \left( -\frac{1}{2} \frac{\vec{a}^2 e^{-2t}}{\Delta_t} \right) \exp \left( \frac{e^{-t} \vec{x} \cdot \vec{a}}{\Delta_t} \right) \quad (13) $$

can be interpreted in a field-theoretic (or probabilistic) approach, as a generative function for connected correlations of the variables $\vec{a}$ [51]. Expanding this function at large times, one finds

$$ g(\vec{x}) = \frac{e^{-t}}{\Delta_t} \sum_{i=1}^{d} x_i \langle a_i \rangle + \frac{1}{2} \frac{e^{-2t}}{\Delta_t^2} \sum_{i,j=1}^{d} x_i x_j [\langle a_i a_j \rangle - \langle a_i \rangle \langle a_j \rangle] + O((xe^{-t})^3) \quad (14) $$

where the brackets $\langle \cdot \rangle$ denote the expectation value with respect to the effective distribution $P_0(\vec{a})e^{-\vec{a}^2 e^{-2t}/(2\Delta_t)}$. Using this expansion, one finds at large times:

$$ \log P_t(\vec{x}) = C + \frac{e^{-t}}{\Delta_t} \sum_{i=1}^{d} x_i \langle a_i \rangle - \frac{1}{2\Delta_t} \sum_{i,j=1}^{d} x_i M_{ij} x_j + O((xe^{-t})^3) \quad (15) $$

where $C$ is an $\vec{x}$-independent term and

$$ M_{ij} = \delta_{ij} - e^{-2t} [\langle a_i a_j \rangle - \langle a_i \rangle \langle a_j \rangle] \quad (16) $$

The curvature of $\log P_t(\vec{x})$ depends on the spectrum of the matrix $M$. At large times $M$ is close to the identity, all its eigenvalues are positive. The shape changes qualitatively at the largest time $t_S$ such that the largest eigenvalue of $M$ crosses zero. The speciation time is characterized by a change of curvature of the effective potential $-\log P_t(\vec{x})$.

Notice that the effective measure used for computing the correlations of $a_i$ variables appearing in the definition of $M$ can be substituted at large times by the original measure $P_0(\vec{a})$. Therefore $M \simeq \mathbb{I} - e^{-2t}\hat{C}_0$ where $\hat{C}_0$ is the covariance of the distribution $P_0$, which can be estimated as the covariance of the data. This leads to the (second) general criterion discussed in the main text for the speciation transition.
# B Gaussian mixtures: speciation time and asymptotic stochastic processes in the high-dimensional limit

We consider the case when $P_0(\vec{a})$ is the superposition of two Gaussian clusters of equal weight, with means $\pm \vec{m}$, and the same variance $\sigma^2$. We ensure that the two Gaussians are well separated in the limit of large dimensions, by assuming that $|\vec{m}|^2 = d\tilde{\mu}^2$, with $\sigma$ and $\tilde{\mu}$ of order 1. In the following we present the detail of the analysis of the backward dynamics in the limit of large dimension. In the study of regimes I and II, we shall assume that $P_t^e(\vec{x})$ coincides with its population counterpart $P_t(\vec{x})$, and use the latter. We justify this assumption in the analysis of regime III, by showing that for $\alpha = \frac{\log n}{d}$ finite, the speciation time, and hence regime I and II, take place before the collapse, i.e. when $P_t^e(\vec{x})$ coincides with its population counterpart $P_t(\vec{x})$.

## B.1 Spectral estimate of $t_S$

The speciation time computed from the eigenvalue criterion of the previous section is easily computed in this case. The matrix $M$ is given by $M_{ij} = (1 - \sigma^2 e^{-2t})\delta_{ij} - e^{-2t}m_i m_j$ and its largest eigenvalue is $(1 - \sigma^2 e^{-2t} - d\tilde{\mu}^2 e^{-2t})$. We get therefore in the large $d$ limit $t_S = \frac{1}{2} \log(d\tilde{\mu}^2)$ which up to subleading corrections identifies the speciation timescale as

$$t_S = \frac{1}{2} \log(d).$$

## B.2 Asymptotic stochastic process in regime I and symmetry breaking

In the large dimensional limit, we can provide a full analytic study of the dynamics in regime I, i.e. at the beginning of the backward process following [7]).

As shown in Sect. C of this Appendix, when one studies the dynamics at times $t > t_C$ the empirical distribution at time $t$, $P_t^e(\vec{x})$ is well approximated by $P_t(\vec{x})$ which is the convolution of the initial distribution $P_0$ (a mixture of Gaussians centered at $\pm \vec{m}$) and the diffusion kernel proportional to $e^{-(\vec{x}-\vec{a}e^{-t})^2/2}$. The explicit expression is thus

$$P_t(\vec{x}) = \frac{1}{2\sqrt{2\pi\Gamma_t^d}} \left[ e^{-(\vec{x}-\vec{m}e^{-t})^2/(2\Gamma_t)} + e^{-(\vec{x}+\vec{m}e^{-t})^2/(2\Gamma_t)} \right] \quad (17)$$

where $\Gamma_t = \sigma^2 e^{-2t} + \Delta_t$ goes to 1 at large times. The log of this probability is

$$\log P_t(\vec{x}) = -\frac{\vec{x}^2}{2\Gamma_t} + \log \cosh \left( \vec{x} \cdot \vec{m} \frac{e^{-t}}{\Gamma_t} \right) \quad (18)$$

and hence the score reads

$$S_i(\vec{x}) = -\frac{x_i}{\Gamma_t} + m_i \frac{e^{-t}}{\Gamma_t} \tanh \left( \vec{x} \cdot \vec{m} \frac{e^{-t}}{\Gamma_t} \right) \quad (19)$$

At the beginning of the backward process (for $t \gg t_S$), the $x_i$ are i.i.d. Gaussian variables with unit variance. In consequence the overlap with the centers of the Gaussian model, $\vec{m} \cdot \vec{x}(t)$, scales as $\sqrt{d}$ in the large dimensional limit. Strictly speaking, regime I is defined as the scaling regime (or the time-window) in which $\vec{m} \cdot \vec{x}(t)$ keeps this scaling with $d$. Introducing the quantity $q(t) = \frac{1}{\sqrt{d}} \vec{m} \cdot \vec{x}(t)$ and using the notation $t_S = (1/2) \log d$, one can therefore write the backward equation on each $x_i$ as:

$$-dx_i = x_i + 2 \left( -\frac{x_i}{\Gamma_t} + m_i \frac{e^{-t}}{\Gamma_t} \tanh \left( q(t) \frac{e^{-(t-t_S)}}{\Gamma_t} \right) \right) dt + d\eta_i(t) \quad (20)$$

where $d\eta_i(t)$ is square root of two times the Brownian motion. This shows that in regime I, on the timescale $t_S$, each $x_i$ satisfies a Langevin equation in which the interactions with all the other variables is through the fluctuating rescaled overlap $q(t)$. It is a kind of dynamical mean-field equation [8].

In order to obtain the equation on $q(t)$, one can sum the equation on all $x_i$. One therefore finds that this projection of the trajectory point $\vec{x}(t)$ on the direction linking the centers of the two Gaussians in the mixture satisfies the closed backward stochastic differential equation cited in the main text:

$$-dq = -\frac{\partial V(q, t)}{\partial q} dt + d\eta(t) \quad (21)$$
where $d\eta(t)$ is square root of two times the Brownian motion, and the potential $V(q, t)$ reads

$$V(q, t) = \frac{1}{2}q^2 - 2\tilde{\mu}^2 \log \cosh \left( qe^{-t}\sqrt{d} \right) \qquad (22)$$

Clearly this potential is quadratic at times $t \gg (1/2) \log d$, and it develops a double well structure when $t \ll (1/2) \log d$. Equation (7) is a Langevin equation in an inverted potential $V(q, t)$. As explained in the main text, the resulting dynamics leads to the symmetry breaking in which trajectories commit to one of the two classes. The scaling variable controlling such phenomenon is $e^{-\sqrt{t}d}$, as also found by the other means explained in the main text and this Appendix.

Note that the time where the curvature at $q = 0$ vanishes, $\frac{\partial^2 V}{\partial q^2}(q = 0, t^*) = 0$, is

$$t^* = \frac{1}{2} \log(2d\tilde{\mu}^2). \qquad (23)$$

In the main text we identify the speciation time with $t_S = \frac{1}{2} \log(d\tilde{\mu}^2)$, or equivalently $e^{-2t_S}d\tilde{\mu}^2 = 1$ but any other choice of constant different from one, e.g. $1/2$ as above, would be equivalent in the large $d$ limit (it gives subleading correction to $t_S$ in the large $d$ limit).

## B.3 Analytic computation of the speciation time defined from cloning

We now consider the criterion for speciation used in numerical experiments, namely the probability $\phi(t)$ that two trajectories, cloned at time $t$, end at time 0 in the same class. For gaussian mixtures, this can be computed as follows.

$$P(t) = \left[ \int_{\vec{x}_1 \cdot \vec{m} > 0, \vec{x}_2 \cdot \vec{m} > 0} + \int_{\vec{x}_1 \cdot \vec{m} < 0, \vec{x}_2 \cdot \vec{m} < 0} \right] d\vec{x}_1 d\vec{x}_2 d\vec{y} \left( P(\vec{x}_1, 0|\vec{y}, t) P(\vec{x}_2, 0|\vec{y}, t) P(\vec{y}, t) \right) \qquad (24)$$

where $P(\vec{x}_1, 0|\vec{y}, t)$ is the probability that the backward process ends in $\vec{x}_1$ knowing that it was in $\vec{y}$ at time $t$. This conditional probability can be rewritten as $P(\vec{x}_1, 0|\vec{y}, t) = P(\vec{x}_1, 0, \vec{y}, t)/P(\vec{y}, t)$. We can now use the forward process to compute all the probabilities involved, using $P(\vec{x}_1, 0, \vec{y}, t) = P(\vec{y}, t|\vec{x}_1, 0)P(\vec{x}_1, 0)$:

$$P(\vec{x}_1, 0, \vec{y}, t) = P_0(\vec{x}_1) \frac{e^{-\frac{(\vec{y}-\vec{x}_1 e^{-t})^2}{2\Delta_t}}}{\sqrt{2\pi \Delta_t^d}}$$

hence leading to:

$$P(\vec{x}_1, 0|\vec{y}, t) = P_0(\vec{x}_1) \frac{e^{-\frac{(\vec{y}-\vec{x}_1 e^{-t})^2}{2\Delta_t}}}{\sqrt{2\pi\Delta_t^d}} \frac{1}{P(\vec{y}, t)} \qquad (25)$$

We now insert this expression into the equation for $P(t)$ and integrate over $\vec{x}_1$ and $\vec{x}_2$. The integrals can be done by noticing that this is equivalent to drawing both $x_1$ and $x_2$ from one Gaussian with probability $1/4$ and from the other Gaussian with probability $1/4$. One thus finds:

$$\phi(t) = \int \frac{d\vec{y}}{2\sqrt{2\pi(\sigma^2 e^{-2t} + \Delta_t)^d}} \frac{e^{-\frac{2}{2(\sigma^2 e^{-2t} + \Delta_t)}(\vec{y}-\vec{x}_1 e^{-t})^2}}{e^{-\frac{2}{2(\sigma^2 e^{-2t} + \Delta_t)}(\vec{y}-\vec{x}_1 e^{-t})^2} + e^{-\frac{2}{2(\sigma^2 e^{-2t} + \Delta_t)}(\vec{y}+\vec{\mu}e^{-t})^2}} \qquad (26)$$

All the integrals on the components orthogonal to $\vec{m}$ factorize and can be done: they give one. One ends up with a one-dimensional integral equal to:

$$\phi(t) = \frac{1}{2} \int dy \frac{G(y, me^{-t}, \Gamma_t)^2 + G(y, -me^{-t}, \Gamma_t)^2}{G(y, me^{-t}, \Gamma_t) + G(y, -me^{-t}, \Gamma_t)}, \qquad (27)$$

where $G(y, u, v)$ is a Gaussian probability density function for the real variable $y$, with mean $u$ and variance $v$, $m = |\vec{m}| = \tilde{\mu}\sqrt{d}$, $\Gamma_t = \Delta_t + \sigma^2 e^{-2t}$. This one-dimensional integral is easily done numerically. For large $d$, the probability $\phi(t)$ that the two clones end up in the same cluster is a decreasing function of $t$, going from $\phi(0) = 1$ to $\phi(\infty) = 1/2$, and the speciation time for this Gaussian mixture model, estimated from $\phi(t_s) = .775$ is of order $t_s \simeq (1/2) \log d$ in the limit of large dimensions.
## B.4 Asymptotic stochastic process in regime II

At the end of regime I the overlap of $\vec{x}(t)$ with the centers of the Gaussian model diverge: for some trajectories $q(t)$ goes to plus infinity, for others to minus infinity. This corresponds to a change of scaling regime. In fact, beyond regime I, the overlap $\vec{x}(t) \cdot \vec{m}$ scales proportionally to $d$.

After the speciation transition, trajectories are committed to a given center. Conditioning on trajectories which correspond to the center $+\vec{m}$, and using that $\vec{x}(t) \cdot \vec{m} \to +\infty$, the score (19) simplifies to:

$$S_i^+(\vec{x}) = -\frac{x_i}{\Gamma_t} + m_i \frac{e^{-t}}{\Gamma_t} \qquad (28)$$

This is the score of the backward process of a single Gaussian centered in $+\vec{m}$. The resulting equation on $x_i$ is therefore:

$$-dx_i = \left( -\frac{x_i}{\Gamma_t} + m_i \frac{e^{-t}}{\Gamma_t} \right) dt + d\eta_i(t) \qquad (29)$$

where $d\eta_i(t)$ is square root of two times a Brownian motion. This equation is the one of the backward process of the single Gaussian centered in $+\vec{m}$, and therefore guarantees that all trajectories evolving with such equation will generate the single Gaussian centered in $+\vec{m}$.

We have focused on the trajectories committed to the center $+\vec{m}$. One can proceed analogously for the ones committed to the center $-\vec{m}$. The results above still holds with $m_i$ replaced by $-m_i$.

We conclude this section on regime I and II by stressing that we expect that our analysis can be generalized to other models. The Curie-Weiss already studied in [7] is an example. The dynamics in regime II has also been studied in [21] using a model for the score.

# C Gaussian Mixtures: collapse time

## C.1 Setting and methods

In order to study the structure of $P_t^e(\vec{x})$ around a point $\vec{x} = \vec{a}^1 e^{-t} + \vec{z} \sqrt{\Delta_t}$, where $\vec{z}$ has i.i.d. Gaussian distributed components with mean zero and variance one, we start from the representation $P_t^e(\vec{x}) = [Z_1 + Z_{2...n}] / \sqrt{2\pi\Delta_t^d}$ where $Z_1 = e^{-(\vec{x}-\vec{a}^1 e^{-t})^2/(2\Delta_t)} = e^{-(\vec{z}^2)/2}$ and

$$Z_{2...n} = \sum_{\mu=2}^{n} e^{-(\vec{x}-\vec{a}^{\mu} e^{-t})^2/(2\Delta_t)} \qquad (30)$$

We now study $P_t(\vec{x})$ in the limit of large $d$, large $n$, keeping $\alpha = \log n/d$ fixed. The first piece behaves as $Z_1 \simeq e^{-d/2}$. In the second piece, $Z_{2...n}$, we assume, without loss of generality, that the first $n_+$ data points were sampled from the Gaussian with mean $\vec{m}$, and the last $n_- = n - n_+$ ones were sampled from the Gaussian with mean $-\vec{m}$. In the large $n$ limit, because of the law of large numbers $n_+/n$ goes to 1/2. We can write $Z_{2...n} = Z_+ + Z_-$, where

$$Z_+ = \sum_{\mu=2}^{n_+} e^{-(\vec{x}-\vec{a}^\mu e^{-t})^2/(2\Delta_t)} \qquad (31)$$

$$Z_- = \sum_{\mu=n_++1}^{n} e^{-(\vec{x}-\vec{a}^\mu e^{-t})^2/(2\Delta_t)} \qquad (32)$$

$Z_1$, $Z_+$ and $Z_-$ give the contribution to $P_t(\vec{x})$ from, respectively, the data point $\mu = 1$ (which is supposed to have been generated from the gaussian with mean $\vec{m}$), the other $n_+ - 1$ data points generated from this same gaussian with mean $\vec{m}$, and the $n_-$ data points generated from the gaussian with mean $-\vec{m}$. The two quantities $Z_\pm$ are partition functions and we shall show that, for typical values of $\vec{x} = \vec{a}^1 e^{-t} + \vec{z} \sqrt{\Delta_t}$, the free energies $(1/d) \log Z_\pm$ concentrates in the large $d$ limit around a value $\psi_\pm(t)$. Furthermore, $\psi_- < \psi_+$. Therefore in the large $d$ limit the collapse transition is identified as the time $t_c$ such that $\psi_+(t_c) = -1/2$.

As already pointed out in the main text, in order to obtain $Z_\pm$, standard concentration methods, e.g. central limit theorem, do not apply as each term of the sum corresponds to the exponential of a random variable scaling as $\log n$ [4]. $Z_\pm$ correspond to the partition function of a system with $n-1$ independent 'random energy levels'. This is some elaboration of the 'Random Energy Model' which was introduced originally in [19] as a simple model of
glass transition, and studied in the limit we are interested in, i.e. $n, d$ going to infinity with $\alpha = \frac{\log n}{d}$ fixed. This model has played a central role in the physics of glassy systems [34]. It was first solved with methods of theoretical physics [19]. It was then thoroughly studied by probabilistic rigorous methods [41].

For a mathematical introduction that make connections with physics and computer science see the book [33]. For our purposes, it is also useful the recent work [31] where a similar problem was addressed in the context of dense associative memories.

## C.2 Computation of $\psi_+$

Given a time $t$, $Z_+$ is the partition function of a system with $n_+ - 1$ independent 'energy levels' $E^\mu = (\vec{x} - \vec{a}^\mu e^{-t})^2/(2\Delta_t)$ at thermal equilibrium at a temperature 1. All these energies are independent and identically distributed, and they are of order $d$. Writing $\varepsilon^\mu = E^\mu/d$, we denote by $\rho_t(\varepsilon)$ the probability density function from which these reduced energies are drawn independently. This density $\rho_t(\varepsilon)$ satisfies a large deviation principle with a rate function $f_t(\varepsilon)$:

$$ \rho_t(\varepsilon) = e^{-d f_t(\varepsilon)} \qquad (33) $$

We define $g_t(\lambda)$ the Legendre transform $g_t(\lambda) = \max_{\varepsilon}[-f_t(\varepsilon) - \lambda\varepsilon]$ and use

$$ e^{dg_t(\lambda)} = \int d\varepsilon \, e^{-d[f_t(\varepsilon) + \lambda\varepsilon]} = \mathbb{E}_+ e^{-\lambda(\vec{x} - \vec{a}e^{-t})^2/2\Delta_t} \qquad (34) $$

where $\mathbb{E}_+$ is an expectation on $\vec{a}$ drawn from the gaussian measure of mean $\vec{m}$ and variance $\sigma^2$. The right-hand side of (34) is the convolution of two gaussians which is easily computed. Using the fact that $(\vec{x} - \vec{m}e^{-t})^2/d$ concentrates at large $d$ around $\sigma_t^2 + \Delta_t$, where $\sigma_t^2 = \sigma^2 e^{-2t}$, we get:

$$ g_t(\lambda) = \frac{1}{2} \log \frac{\Delta_t}{\Delta_t + \lambda\sigma_t^2} - \frac{1}{2} \frac{\lambda(\Delta_t + \sigma_t^2)}{\Delta_t + \lambda\sigma_t^2} \qquad (35) $$

The Legendre transform can be done analytically. The maximum of $g_t(\lambda) + \lambda\varepsilon$ is obtained at

$$ \lambda^*(\varepsilon) = \frac{\sigma_t^2 - 4\Delta_t\varepsilon + \sqrt{\sigma_t^4 + 8\Delta_t^2\varepsilon + 8\Delta_t\sigma_t^2\varepsilon}}{4\sigma_t^2\varepsilon} \qquad (36) $$

The inverse transform is obtained from the maximum of $-f_t(\varepsilon) - \lambda\varepsilon$ which is obtained at

$$ \varepsilon^*(\lambda) = -\frac{dg_t}{d\lambda} = \frac{\Delta_t^2 + 2\Delta_t\sigma_t^2 + \lambda\sigma_t^2}{2(\Delta_t + \lambda\sigma_t^2)^2} \qquad (37) $$

and the rate function is

$$ f_t(\varepsilon) = \frac{\Delta_t}{2\sigma_t^2(\sigma_t^2 + A)} [-8\Delta_t\varepsilon + (1 - 6\varepsilon)\sigma_t^2 + (1 + 2\varepsilon A)] - \frac{1}{2} \log \frac{4\Delta_t\varepsilon}{\sigma_t^2 + A} \qquad (38) $$

where

$$ A = \sqrt{\sigma_t^4 + 8\varepsilon\Delta_t(\Delta_t + \sigma_t^2)} \qquad (39) $$

The partition function is $Z_+ = \sum_{\mu=2}^{n_+} e^{-d\varepsilon^\mu}$. The annealed approximation to this partition function is $Z_+^{ann} = n_+ \int d\varepsilon e^{-d(\varepsilon+f_t(\varepsilon))}$. Using Laplace's method, we see that this integral is dominated by $\varepsilon = \varepsilon^*(\lambda = 1) = 1/2$, and gives

$$ \psi_+^{ann} = \lim_{d \to \infty} \frac{1}{d} \log Z_+^{ann} = \alpha - \frac{1}{2} - f_t\left(\frac{1}{2}\right) = \alpha + g_t(1) = \alpha + \frac{1}{2} \log \frac{\Delta_t}{\Delta_t + \sigma_t^2} - \frac{1}{2} \qquad (40) $$

The methods for analyzing random energy models are standard (see [19, 33]). For a presentation close to the developments which we use here, see [31], in particular its Appendix A. Using either with first and second moment methods, or with replicas, we find that the annealed free energy is exact, namely $\psi_+ = \psi_+^{ann}$ when $t < t_{cond}$, but at $t > t_{cond}$ a condensation phenomenon occurs, the partition function is dominated by the levels with lowest energy, $\varepsilon_{min}$ given by the smallest root of $\alpha = \hat{f}_t(\varepsilon)$. The condensation time $t_{cond}$ is thus characterized by the fact that
$\epsilon_{min} = 1/2$, which gives:

$$ \alpha = f_{t_{cond}}(1/2) \qquad (41) $$

So the final expression for the free energy of this random energy model is:

$$ \begin{aligned} \psi_+ = \lim_{d \to \infty} \frac{1}{d} \log Z_+ &= \alpha + \frac{1}{2} \log \frac{\Delta_t}{\Delta_t + \sigma_t^2} - \frac{1}{2} \quad \text{if } t < t_{cond} \\ &= -\frac{1}{2} \quad \text{if } t > t_{cond} \end{aligned} \qquad (42) $$

and its condensation transition is

$$ t_{cond} = \frac{1}{2} \log \left[ 1 + \frac{\sigma^2}{n^{2/d} - 1} \right] \qquad (43) $$

### C.3 Computation of $\psi_-$

Given $t$, the partition function $Z_-$ is again a partition function of a random energy model, which can be studied following exactly the same steps as for $Z_+$. The computation of $g_t(\lambda)$ is slightly different. Now, $(\vec{x} - \vec{m}e^{-t})^2/d$ concentrates at large $d$ around $4(\vec{m}^2/d)e^{-2t} + \sigma_t^2 + \Delta_t$. This gives:

$$ g_t(\lambda) = \frac{1}{2} \log \frac{\Delta_t}{\Delta_t + \lambda\sigma_t^2} - \frac{1}{2} \frac{\lambda[M_t + \Delta_t + \sigma_t^2]}{\Delta_t + \lambda\sigma_t^2} \qquad (44) $$

where $M_t = 4(\vec{m}^2/d)e^{-2t}$. We now list the changes in the subsequent formulations, keeping the same notations as in the previous section. The maximum of $g_t(\lambda) + \lambda\epsilon$ is obtained at

$$ \lambda^*(\epsilon) = \frac{\sigma_t^2 - 4\Delta_t\epsilon + \sqrt{\sigma_t^4 + 8\Delta_t^2\epsilon + 8\Delta_t\sigma_t^2\epsilon + 8\Delta_t M_t\epsilon}}{4\sigma_t^2\epsilon} \qquad (45) $$

The inverse transform is obtained from the maximum of $-f_t(\epsilon) - \lambda\epsilon$ which is at

$$ \varepsilon^*(\lambda) = -\frac{dg_t}{d\lambda} = \frac{\Delta_t^2 + 2\Delta_t\sigma_t^2 + \lambda\sigma_t^2 + M_t\Delta_t}{2(\Delta_t + \lambda\sigma_t^2)^2} \qquad (46) $$

and the rate function is

$$ f_t(\epsilon) = \frac{\Delta_t [-8\Delta_t\epsilon + (1-6\epsilon)\sigma_t^2 + (1+2\epsilon A)] + M_t [\sigma_t^2 - 8\epsilon\Delta_t + A]}{2\sigma_t^2(\sigma_t^2 + A)} - \frac{1}{2} \log \frac{4\Delta_t\epsilon}{\sigma_t^2 + A} \qquad (47) $$

where

$$ A = \sqrt{\sigma_t^4 + 8\epsilon\Delta_t(\Delta_t + \sigma_t^2) + 8\Delta_t M_t\epsilon} \qquad (48) $$

The annealed approximation to the partition function is $Z_-^{ann} = n_- \int d\epsilon e^{-d(\epsilon+f_t(\epsilon))}$. This integral is dominated by

$$ \epsilon = \epsilon^*(\lambda = 1) = \frac{1}{2} \left[ 1 + \frac{\Delta_t M_t}{(\Delta_t + \sigma_t^2)^2} \right] \qquad (49) $$

and gives

$$ \psi_-^{ann} = \lim_{d \to \infty} \frac{1}{d} \log Z_-^{ann} = \alpha + g_t(1) = \alpha + \frac{1}{2} \log \frac{\Delta_t}{\Delta_t + \sigma_t^2} - \frac{1}{2} - \frac{1}{2} \frac{M_t}{\Delta_t + \sigma_t^2} \qquad (50) $$

The annealed free energy is exact when $t < t_{cond}$. The condensation takes place at $t > t_{cond}$, then the partition function is dominated by the levels with lowest energy, $\epsilon_{min}$ given by the smallest root of $\alpha = \hat{f}_t(\epsilon)$. The
condensation time $t_{cond}$ is thus characterized by the fact that $\epsilon_{min} = \epsilon^*(\lambda = 1)$, which gives:

$$ \alpha = f_{t_{cond}} \left( \epsilon = \frac{1}{2} \left[ 1 + \frac{\Delta_t M_t}{(\Delta_t + \sigma_t^2)^2} \right] \right) \qquad (51) $$

So the final expression for the free energy of this random energy model is:

$$ \begin{aligned} \psi_{-} = \lim_{d \to \infty} \frac{1}{d} \log Z_{-} &= \alpha + \frac{1}{2} \log \frac{\Delta_t}{\Delta_t + \sigma_t^2} - \frac{1}{2} - \frac{1}{2} \frac{M_t}{\Delta_t + \sigma_t^2} && \text{if } t < t_{cond} \\ &= -\frac{1}{2} \left[ 1 + \frac{\Delta_t M_t}{(\Delta_t + \sigma_t^2)^2} \right] && \text{if } t > t_{cond} \end{aligned} \qquad (52) $$

and its condensation transition is

$$ t_{cond} = \frac{1}{2} \log \left[ 1 + \frac{\sigma^2}{n^{2/d} - 1} \right] \qquad (53) $$

It is easy to see that at any time $\psi_{-} < \psi_{+}$.

## C.4 Final expression

Going back to the representation $P_t^e(\vec{x}) = [Z_1 + Z_{2...n}] / \sqrt{2\pi\Delta_t^d}$, where $Z_1 = e^{-(\vec{x}-\vec{a}^1e^{-t})^2/(2\Delta_t)} = e^{-(\vec{z}^2)/2}$, it is made of three contributions:

1. $Z_1$ which behaves at large $d$ as $e^{-d/2}$

2. $Z_+$ which behaves at large $d$ as $e^{-d\psi_+}$

3. $Z_-$ which behaves at large $d$ as $e^{-d\psi_-}$

In the limit of large dimensions $d \to \infty$, the contribution from $Z_-$ is irrelevant because $\psi_- < \psi_+$. The comparison of the two other terms depends on whether $\psi_+$ is larger or smaller than $-1/2$. We thus find two time regimes separated by a collapse transition time $t_C$.

* at large times, $t > t_C$, $\psi_+ > -(1/2)$ and the probability $P_t^e(\vec{x})$ is dominated by the term $Z_+$. This is the regime which is not collapsed.

* at small times, $t < t_C$, $\psi_+ < -(1/2)$ and the probability $P_t^e(\vec{x})$ is dominated by the term $Z_1 = e^{-(\vec{x}-\vec{a}^1e^{-t})^2/(2\Delta_t)}$. When used in the backward diffusion, this gives a score that attracts $\vec{x}$ towards $\vec{a}^1$ at short times.

We have thus shown that, in the backward process, trajectories which are at time $t$ at typical points of the form $\vec{x} = \vec{a}^1e^{-t} + \vec{z}\sqrt{\Delta_t}$ are attracted towards $\vec{a}^1$ if $t < t_C$. This identifies $t_C$ as the collapse time.

Looking at the explicit expression (42) of $\psi_+$, we see that the collapse time $t_C$ is identical to the condensation time of the random energy model, $t_{cond}$. This result, derived exactly in this case of gaussian mixtures, is actually a general result, as discussed in the main paper. It could also be shown using the REM approach, provided the distribution of energies satisfies a large deviation principle.

# D Details on the numerical experiments

## D.1 Gaussian mixtures

**Speciation experiment.** In the case of a mixture of two Gaussian clusters centered on $\pm \vec{m} \in \mathbb{R}^d$ with variance $\sigma^2$, the score function $\vec{F}(\vec{x}, t)$ can be analytically expressed as

$$ \vec{F}(\vec{x}(t), t) = \vec{m} \frac{e^{-t}}{\Gamma_t} \tanh \left( \frac{e^{-t}}{\Gamma_t} \vec{x}(t) \cdot \vec{m} \right) - \frac{\vec{x}(t)}{\Gamma_t}, \qquad (54) $$

where $\Gamma_t = \Delta_t + \sigma^2 e^{-2t}$, with $\Delta_t = 1 - e^{-2t}$. We can then discretize the stochastic differential equation associated to the backward process as

$$ \vec{x}(t+1) = \vec{x}(t) + \eta [\vec{x}(t) + 2\mathcal{F}(\vec{x}(t), t)] + \xi \sqrt{2\eta}, \qquad (55) $$
where $\xi \sim \mathcal{N}(0, I)$, $\vec{x}_T \sim \mathcal{N}(0, I)$ and $\eta = T/L$, with $T = 10$ the time horizon and $L = 1000$ the number of discrete steps. Two clones $\vec{x}_1(t)$ and $\vec{x}_2(t)$ share the same trajectory until a given time below which they have independant noise realizations $\xi$. At the end of the backward diffusion, one can then recover the class of the clone $i$ by projecting it onto $\vec{m}$ as sign $(\vec{m} \cdot \vec{x}_i(0))$. In the experiment corresponding to Fig. 2 of the main text, we use $\vec{m} = [1, \dots, 1]$, $\sigma^2 = 1$, and each point is obtained by averaging over 10 000 initial conditions.

**Collapse experiment.** The main component of the collapse experiment is the excess entropy density $f(t)$ which vanishes for $t \le t_C$. When dealing with Gaussian mixtures, $f(t)$ can be derived exactly and we additionally establish a strategy to approximate it numerically as follows. At a given time $t$, we draw $n'$ samples $\vec{x}^{(\nu)}(t) = e^{-t}a^{(\mu)} + \sqrt{1 - e^{-2t}}\xi^{(\nu)}$, where $\mu$ is chosen uniformly in $\{1, \dots, n\}$. The entropy $s(t) = -\int d\vec{x} P_t(\vec{x}) \log P_t(\vec{x})$ can then be approximated as the empirical average over the $n'$ samples, leading to the empirical estimate $f^e(t)$ of the excess entropy, using

$$s^e(t) = -\frac{1}{n'd} \sum_{\nu=1}^{n'} \log P_t^e(\vec{x}_t^{(\nu)}), \qquad (56)$$

where $P_t^e(\vec{x})$ is given as usual in terms of the original dataset by

$$P_t^e(\vec{x}) = \frac{1}{n} \sum_{\mu=1}^{n} \frac{1}{\sqrt{2\pi\Delta_t^d}} \exp\left(-\frac{(\vec{x} - \vec{a}^\mu e^{-t})^2}{2\Delta_t}\right). \qquad (57)$$

This approximation is valid for $t > t_C$, where $P_t(\vec{x}) \approx P_t^e(\vec{x})$. As shown in Fig. 2 of the main text where we fix $n = 20000$ and $n' = 250000$ for several $d$, this estimate fits quite remarkably the analytical curve. Such a procedure therefore allows for the numerical derivation of the collapse $t_C$ from the dataset only, under the assumption that one uses the exact empirical score.

## D.2 Realistic datasets

**Denoising Diffusion Probabilistic Models.** In the second part of the paper, we learn the score by training a Denoising Diffusion Probabilistic Model (DDPM) in discrete time, as introduced by [25]. In this context, the forward process has a variance schedule $\{\beta_t^{t'}/t'\}_{t'=1}^L$, where $L$ is the time horizon given as a number of steps, fixed to 1000 in our experiments. In our case, the variance is evolving linearly from $\beta_1 = 10^{-4}$ to $\beta_{1000} = 2 \times 10^{-2}$. A sample at timestep $t'$, denoted $\vec{x}(t')$ can therefore be expressed readily from its initial state, $\vec{x}(0) = \vec{a}$, as

$$\vec{x}(t') = \sqrt{\bar{\alpha}(t')}\vec{a} + \sqrt{1 - \bar{\alpha}(t')}\xi(t'), \qquad (58)$$

where $\bar{\alpha}(t') = \prod_{s=1}^{t'}(1 - \beta_s)$ and $\xi$ is a standard and centered Gaussian noise. This equation in fact corresponds to the discretization of the Ornstein-Uhlenbeck Eq. (1) given in the main text under the following reparameterization of the timestep $t'$,

$$t = -\frac{1}{2} \log (\bar{\alpha}(t')), \qquad (59)$$

where $t$ is the time as defined in the main text. We then train a neural network to learn $\xi_\theta(\vec{x}(t'), t')$, the noise at time $t'$ by iteratively optimizing the loss function

$$\mathcal{L}(\theta) = \|\xi(t') - \xi_{\theta}(\vec{x}(t'), t')\|^2. \qquad (60)$$

Once learned, the denoiser can then be used to generate a new sample as

$$\vec{x}(t' - 1) = \frac{1}{\sqrt{1 - \beta_{t'}}} \vec{x}(t') - \frac{\beta_{t'}}{\sqrt{(1 - \beta_{t'}) (1 - \bar{\alpha}_{t'})}} \xi_{\theta}(\vec{x}(t'), t') + \sqrt{\beta_{t'}} \vec{z}, \qquad (61)$$

where $t'$ runs from $T$ to one, $\vec{x}(T) \sim \mathcal{N}(0, I)$, and $\vec{z} \sim \mathcal{N}(0, I)\delta(t' > 1)$, which, in the continuous limit, is equivalent to the score-based diffusion studied in the main text [32], where $-\xi_\theta(\vec{x}(t'), t') / \sqrt{1 - \bar{\alpha}_{t'}}$ corresponds to the score.
**Datasets and preprocessing.** Our experiments are based on five datasets: MNIST [29], CIFAR [28], ImageNet16/32 [14], and LSUN [50]. For each dataset, we focus only on two well-distinct and balanced classes that are: zero and seven for MNIST ($n = 10000$), horses and cars for CIFAR ($n = 3000$), lesser pandas and seashores for ImageNet16 and 32 ($n = 2000$), and churches and conference rooms for LSUN ($n = 40000$ or $n = 200$). In the case of MNIST, the images are zero-padded from their original size to 32 × 32. For LSUN, images are center-cropped at size 256 × 256 before being resized to 64 × 64. All the datasets but MNIST have three color channels that are kept for the training, and the data are centered in each channel. The different datasets therefore allow to span several values of $d$ ranging from 1024 (for MNIST) to 12288 for LSUN. The values of $\Lambda$ given in the main text corresponding to the largest eigenvalue of the covariance matrix are computed from the first channel only.

**Architecture and training.** For all the datasets, the denoiser $\xi_{\theta}(\tilde{x}(t'), t')$ has a U-Net architecture similar to previous applications of DDPMs [25, 44]. More precisely, it implements four resolution levels, each having two convolutional residual blocks with group normalization, and self attention is applied to the two intermediary levels. The time is added into each block through a standard sinusoidal position embedding, resulting in a total of 25.7 million parameters. The models are then trained using ADAM optimizer with a fixed learning rate of $10^{-4}$ for 350k steps with batches of size 128, except for LSUN for which we reduced it to 64. All the $n$ images were used to train the models, without any data augmentation. It is worth mentioning that the results exposed in the main text do not depend on the formal architecture of the denoiser, nor the optimization scheme, as long as the network is able to learn a score close to the true one (or equivalently, a sufficiently good denoiser). In Figs. 6 to 10, we show several samples generated from our six trained models.

**Speciation experiment.** To estimate numerically the speciation time from the probability $\phi(t)$ that the two clones end up in the same class, we need to classify the generated samples. To do so, we train the PyTorch [36] implementation of the ResNet-18 architecture [24] with pre-trained weights on ImageNet. One model is trained for each of the dataset using all the $n$ samples. The resulting classifiers yield 95% test accuracy at worst on LSUN, and more than 99% for the other datasets.

**Collapse experiment.** To study the collapse, we devise two methods. First, we exploit the cloned trajectories but we now focus on the indices of the nearest neighbors in the training set at the end of the backward process. Using a small number of training data ($n = 2000$ for ImageNet and $n = 200$ for LSUN) ensures the backward dynamics to collapse onto one of the training datapoint. When $t < t_C$, the two clones have the same nearest neighbor, with a very small Euclidean distance, showing that the collapse has indeed taken place. This approach allows us to compute $\phi_C(t)$, the probability that the two clones collapse on the same datapoint. We illustrate this phenomenon on the LSUN dataset trained with either $n = 200$ or $n = 40000$ in Fig. 11. For a random batch of generated samples, we show the four nearest neighbors $\{a^{\mu_1}, a^{\mu_2}, a^{\mu_3}, a^{\mu_4}\}$ (ordered by increasing distance) in the training set. In the case of a model trained with a small number of datapoints, the generated images collapse onto $\mu_1$, as shown in Fig. 11a. The model being fixed, one solution to avoid the collapse is to increase $n$, resulting in the generation of new samples, as illustrated in Fig. 11b. In a second approach, we exploit the indices of nearest neighbors in the training set *during* the backward process. At any times $t < t_C$, the index of the nearest neighbor $\mu_1(t)$ is therefore fixed. At $t = 0$, the distance $d_{\mu_1}$ to the nearest neighbor is vanishing, as seen before. The average of this distribution for $n = 4000$ generated samples yields the estimate $\hat{t}_C$. The entropy appearing in $f^e(t)$ is computed the same way as in Gaussian mixtures (see Sect. D.1), making use of the training sets and the forward process, with $n' = 400000$.
<figure><img src="image_8.png" /><figcaption>Figure 6: MNIST one and eight generated samples with $n = 10\,000$.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Figure 7: CIFAR horses and cars generated samples with $n = 3000$.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Figure 8: ImageNet16 lesser pandas and seashores generated samples with $n = 2000$.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Figure 9: ImageNet32 lesser pandas and seashores generated samples with $n = 2000$.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Figure 10: LSUN churches and conference rooms generated samples with (a) $n = 200$ or (b) $n = 40\,000$.</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>**Figure 11:** Illustration of the collapse on the LSUN dataset trained with (a) *n* = 200 or (b) *n* = 40 000 datapoints. Each row corresponds to a generated sample shown in the left-most column. The other columns display the first four nearest neighbors in the training set, denoted as *a*<sup>μ1</sup>, *a*<sup>μ2</sup>, *a*<sup>μ3</sup>, and *a*<sup>μ4</sup>. Below each of them is indicated the *L*<sub>2</sub> norm *d* to the generated sample.</figcaption></figure>