Why Diffusion Models Don’t Memorize: The Role of Implicit Dynamical Regularization in Training

Tony Bonnaire$^\dagger$, LPENS, Université PSL, Paris, [email protected]

Raphaël Urfin$^\dagger$, LPENS, Université PSL, Paris, [email protected]

Giulio Biroli, LPENS, Université PSL, Paris, [email protected]

Marc Mézard, Department of Computing Sciences, Bocconi University, Milano, [email protected]

$\dagger$Equal contribution.

Diffusion models have achieved remarkable success across a wide range of generative tasks. A key challenge is understanding the mechanisms that prevent their memorization of training data and allow generalization. In this work, we investigate the role of the training dynamics in the transition from generalization to memorization. Through extensive experiments and theoretical analysis, we identify two distinct timescales: an early time $\tau_\mathrm{gen}$ at which models begin to generate high-quality samples, and a later time $\tau_\mathrm{mem}$ beyond which memorization emerges. Crucially, we find that $\tau_\mathrm{mem}$ increases linearly with the training set size $n$, while $\tau_\mathrm{gen}$ remains constant. This creates a growing window of training times with $n$ where models generalize effectively, despite showing strong memorization if training continues beyond it. It is only when $n$ becomes larger than a model-dependent threshold that overfitting disappears at infinite training times. These findings reveal a form of implicit dynamical regularization in the training dynamics, which allow to avoid memorization even in highly overparameterized settings. Our results are supported by numerical experiments with standard U-Net architectures on realistic and synthetic datasets, and by a theoretical analysis using a tractable random features model studied in the high-dimensional limit.

Diffusion Models [DMs, 1, 2, 3, 4] achieve state-of-the-art performance in a wide variety of AI tasks such as the generation of images [5], audios [6], videos [7], and scientific data [8, 9]. This class of generative models, inspired by out-of-equilibrium thermodynamics [1], corresponds to a two-stage process: the first one, called forward, gradually adds noise to a data, whereas the second one, called backward, generates new data by denoising Gaussian white noise samples. In DMs, the reverse process typically involves solving a stochastic differential equation (SDE) with a force field called score. However, it is also possible to define a deterministic transport through an ordinary differential equation (ODE), treating the score as a velocity field, an approach that is for instance followed in flow matching [10].

Understanding the generalization properties of score-based generative methods is a central issue in machine learning, and a particularly important question is how memorization of the training set is avoided in practice. A model without regularization achieving zero training loss only learns the empirical score, and is bound to reproduce samples of the training dataset at the end of the backward process. This memorization regime [11, 12] is empirically observed when the training set is small and disappears when it increases beyond a model-dependent threshold [13]. Understanding the mechanisms controlling this change of regimes from memorization to generalization is a central challenge for both theory and applications. Model regularization and inductive biases imposed by the network architecture were shown to play a role [14, 15], as well as a dynamical regularization due to the finiteness of the learning rate [16]. However, the regime shift described above is consistently observed even in models where all these regularization mechanisms are present. This suggests that the core mechanism behind the transition from memorization to generalization lies elsewhere. In this work, we demonstrate -- first through numerical experiments, and then via the theoretical analysis of a simplified model -- that this transition is driven by an implicit dynamical bias towards generalizing solutions emerging in the training, which allows to avoid the memorization phase.

Contributions and theoretical picture. We investigate the dynamics of score learning using gradient descent, both numerically and analytically, and study the generation properties of the score depending on the time $\tau$ at which the training is stopped. The theoretical picture built from our results and combining several findings from the recent literature is illustrated in Figure 1. The two main parameters are the size of the training set $n$ and the expressivity of the class of score functions on which one trains the model, characterized by a number of parameters $p$; when both $n$ and $p$ are large one can identify three main regimes. Given $p$, if $n$ is larger than $n^*(p)$ (which depends on the training set and on the class of scores), the score model is not expressive enough to represent the empirical score associated to $n$ data, and instead provides a smooth interpolation, approximately independent of the training set. In this regime, even with a very large training time $\tau\to\infty$, memorization does not occur because the model is regularized by its architecture and the finite number of parameters. When $n<n^*(p)$ the model is expressive enough to memorize, and two timescales emerge during training: one, $\tau_\mathrm{gen}$, is the minimum training time required to achieve high-quality data generation; the second, $\tau_\mathrm{mem}>\tau_\mathrm{gen}$, signals when further training induces memorization, and causes the model to increasingly reproduce the training samples (left panel). The first timescale, $\tau_\mathrm{gen}$, is found independent of $n$, whereas the second, $\tau_\mathrm{mem}$, grows approximately linearly with $n$, thus opening a large window of training times during which the model generalizes if early stopped when $\tau \in [\tau_\mathrm{gen}, \tau_\mathrm{mem}]$.

Our results shows that implicit dynamical regularization in training plays a crucial role in score-based generative models, substantially enlarging the generalization regime (see right panel of Figure 1), and hence allowing to avoid memorization even in highly overparameterized settings. We find that the key mechanism behind the widening gap between $\tau_\mathrm{gen}$ and $\tau_\mathrm{mem}$ is the irregularity of the empirical score at low noise level and large $n$. In this regime the models used to approximate the score provide a smooth interpolation that remains stable for a long period of training times and closely approximates the population score, a behavior likely rooted in the spectral bias of neural networks [18]. Only at very long training times do the dynamics converge to the low lying minimum corresponding to the empirical score, leading to memorization (as illustrated in the one-dimensional examples in the left panel of Figure 1).

The theoretical picture described above is based on our numerical and analytical results, and builds up on previous works, in particular numerical analysis characterizing the memorization--generalization transition [19, 20], analytical works on memorization of DMs [17, 14, 13], and studies on the spectral bias of deep neural networks [18]. Our numerical experiments use a class of scores based on a realistic U-Net [21] trained on downscaled images of the CelebA dataset [22]. By varying $n$ and $p$, we measure the evolution of the sample quality (through FID) and the fraction of memorization during learning, which support the theoretical scenario presented in Figure 1. Additional experimental results on synthetic data are provided in Supplemental Material (SM, Sects. Appendix A and Appendix B). On the analytical side, we focus on a class of scores constructed from random features and simplified models of data, following [17]. In this setting, the timescales of training dynamics correspond directly to the inverse eigenvalues of the random feature correlation matrix. Leveraging tools from random matrix theory, we compute the spectrum in the limit of large datasets, high-dimensional data, and overparameterized models. This analysis reveals, in a fully tractable way, how the theoretical picture of Figure 1 emerges within the random feature framework.

Related works. - The memorization transition in DMs has been the subject of several recent empirical investigations [23, 24, 25] which have demonstrated that state-of-the-art image DMs -- including Stable Diffusion and DALL·E -- can reproduce a non-negligible portion of their training data, indicating a form of memorization. Several additional works [19, 20] examined how this phenomenon is influenced by factors such as data distribution, model configuration, and training procedure, and provide a strong basis for the numerical part of our work.

Setting: generative diffusion and score learning. Standard DMs define a transport from a target distribution $P_0$ in $\mathbb{R}^d$ to a Gaussian white noise $\mathcal{N}(0, \bm{I}_d)$ through a forward process defined as an Ornstein-Uhlenbeck (OU) stochastic differential equation (SDE):

where $\mathrm{d} \bm{\mathrm{B}}(t)$ is square root of two times a Wiener process. Generation is performed by time-reversing the SDE Equation 1 using the score function $\bm{\mathrm{s}}(\bm{\mathrm{x}}, t) = \nabla_{\bm{\mathrm{x}}} \log P_t(\bm{\mathrm{x}})$,

where $P_t(\bm{\mathrm{x}})$ is the probability density at time $t$ along the forward process, and the noise $\mathrm{d} \bm{\mathrm{B}}(t)$ is also the square root of two times a Wiener process. As shown in the seminal works [35, 36], $\bm{\mathrm{s}}(\bm{\mathrm{x}}, t)$ can be obtained by minimizing the score matching loss

where $\Delta_t=1-e^{-2t}$. In practice, the optimization problem is restricted to a parametrized class of functions $\bm{\mathrm{s}}_{\bm{\theta}}(\bm{\mathrm{x}}(t), t)$ defined, for example, by a neural network with parameters $\bm{\theta}$. The expectation over $\bm{\mathrm{x}}$ is replaced by the empirical average over the training set ($n$ iid samples $\bm{\mathrm{x}}^\nu$ drawn from $P_0$),

where $\bm{\mathrm{x}}^\nu_t(\bm{\xi})=e^{-t} \bm{\mathrm{x}}^\nu+\sqrt{\Delta_t} \bm{\xi}$. The loss in (Equation 4) can be minimized with standard optimizers, such as stochastic gradient descent [SGD, 37] or Adam [38]. In practice, a single model conditioned on the diffusion time $t$ is trained by integrating (Equation 4) over time [39]. The solution of the minimization of Equation 4 is the so-called empirical score (e.g. [12, 11]), defined as $\bm{\mathrm{s}}_{\mathrm{emp}}(\bm{\mathrm{x}}, t) = \nabla_{\bm{\mathrm{x}}} \log P_t^\mathrm{emp}(\bm{\mathrm{x}})$, with

This solution is known to inevitably recreate samples of the training set at the end of the generative process (i.e., it perfectly memorizes), unless $n$ grows exponentially with the dimension $d$ [12]. However, this is not the case in many practical applications where memorization is only observed for relatively small values of $n$, and disappears well before $n$ becomes exponentially large in $d$. The empirical minimization performed in practice, within a given class of models and a given minimization procedure, does not drive the optimization to the global minimum of Equation 4, but instead to a smoother estimate of the score that is independent of the training set with good generalization properties [13], as the global minimum of Equation 3 would do. Understanding how it is possible, and in particular the role played by the training dynamics to avoid memorization, is the central aim of the present work.

Data & architecture. We conduct our experiments on the CelebA face dataset [22], which we convert to grayscale downsampled images of size $d=32\times32$, and vary the training set size $n$ from 128 up to 32768. Our score model has a U-Net architecture [21] with three resolution levels and a base channel width of $W$ with multipliers 1, 2 and 3 respectively. All our networks are DDPMs [2] trained to predict the injected noise at diffusion time $t$ using SGD with momentum at fixed batch size $\min(n, 512)$. The models are all conditioned on $t$, i.e. a single model approximates the score at all times, and make use of a standard sinusoidal position embedding [40] that is added to the features of each resolution. More details about the numerical setup can be found in SM (Appendix A).

Evaluation metrics. To study the transition from generalization to memorization during training, we monitor the loss Equation 4 during training using a fixed diffusion time $t = 0.01$. At various numbers of SGD updates $\tau$, we compute the loss on all $n$ training examples (training loss) and on a held-out test set of 2048 images (test loss). To characterize the score obtained after a training time $\tau$, we assess the originality and quality of samples by generating 10K samples using a DDIM accelerated sampling [41]. We compute (i) the Fréchet-Inception Distance [FID, 42] against 10K test samples which we use to identify the generalization time $\tau_\mathrm{gen}$; and (ii) the fraction of memorized generated samples $f_\mathrm{mem}(\tau)$ granting access to $\tau_\mathrm{mem}$, the memorization time. Following previous numerical studies [20, 19], a generated sample $\bm{\mathrm{x}}_\tau$ is considered memorized if

where $\bm{\mathrm{a}}^{\mu_1}$ and $\bm{\mathrm{a}}^{\mu_2}$ are the nearest and second nearest neighbors of $\bm{\mathrm{x}}_\tau$ in the training set in the $L_2$ sense. In what follows, we choose to work with $k=1/3$ [20, 19], but we checked that varying $k$ to $1/2$ or $1/4$ does not impact the claims about the scaling. Error bars in the figures correspond to twice the standard deviation over 5 different test sets for FIDs, and 5 noise realizations for $\mathcal{L}_\mathrm{train}$ and $\mathcal{L}_\mathrm{test}$. For $f_\mathrm{mem}$, we report the 95% CIs on the mean evaluated with 1, 000 bootstrap samples.

Role of training set size on the learning dynamics. At fixed model capacity ($p=4\times 10^6$, base width $W=32$), we investigate how the training set size $n$ impacts the previous metrics. In the left panel of Figure 2, we first report the FID (solid lines) and $f_\mathrm{mem}(\tau)$ (dashed lines) for various $n$. All trainings dynamics exhibit two phases. First, the FID quickly decreases to reach a minimum value on a timescale $\tau_\mathrm{gen}$ ($\approx100$ K) that does not depend on $n$. In the right panel, the generated samples at $\tau=100$ K clearly differ from their nearest neighbors in the training set, indicating that the model generalizes correctly. Beyond this time, the FID remains flat. $f_\mathrm{mem}(\tau)$ is zero until a later time $\tau_\mathrm{mem}$ after which it increases, clearly signaling the entrance into a memorization regime, as illustrated by the generated samples in the right-most panel of Figure 2, very close to their nearest neighbors. Both the transition time $\tau_\mathrm{mem}$ and the value of the final fraction $f_\mathrm{mem}(\tau_\mathrm{max})$ (with $\tau_\mathrm{max}$ being one to four million SGD steps) vary with $n$. The inset plot shows the normalized memorization fraction $f_\mathrm{mem}(\tau)/f_\mathrm{mem}(\tau_\mathrm{max})$ against the rescaled time $\tau/n$, making all curves collapse and increase at around $\tau/n \approx 300$, showing that $\tau_\mathrm{mem} \propto n$, and demonstrating the existence of a generalization window for $\tau \in \left[\tau_\mathrm{gen}, \tau_\mathrm{mem}\right]$ that widens linearly with $n$, as illustrated in the left panel of Figure 1.

As highlighted in the introduction, memorization in DMs is ultimately driven by the overfitting of the empirical score $\bm{\mathrm{s}}_\mathrm{mem}(\bm{\mathrm{x}}, t)$. The evolution of $\mathcal{L}_\mathrm{train}(\tau)$ and $\mathcal{L}_\mathrm{test}(\tau)$ at fixed $t=0.01$ are shown in the middle panel of Figure 2 for $n$ ranging from 512 to 32768. Initially, the two losses remain nearly indistinguishable, indicating that the learned score $\bm{\mathrm{s}}_{\bm{\theta}}(\bm{\mathrm{x}}, t)$ does not depend on the training set. Beyond a critical time, $\mathcal{L}_\mathrm{train}$ continues to decrease while $\mathcal{L}_\mathrm{test}$ increases, leading to a nonzero generalization loss whose magnitude depends on $n$. As $n$ increases, this critical time also increases and, eventually, the training and test loss gap shrinks: for $n=32768$, the test loss remains close to the training loss, even after 11 million SGD steps. The inset shows the evolution of both losses with $\tau/n$, demonstrating that the overfitting time scales linearly with the training set size $n$, just like $\tau_\mathrm{mem}$ identified in the left panel. Moreover, there is a consistent lag between the overfitting time and $\tau_\mathrm{mem}$ at fixed $n$, reflecting the additional training required for the model to overfit the empirical score sufficiently to reproduce the training samples, and therefore to impact the memorization fraction.

Memorization is not due to data repetition. We must stress that this delayed memorization with $n$ is not due to the mere repetition of training samples, as a first intuition could suggest. In SM Sects. Appendix A and Appendix B, we show that full-batch updates still yield $\tau_\mathrm{mem}\propto n$. In other words, even if at fixed $\tau$ all models have processed each sample equally often, larger $n$ consistently postpone memorization. This confirms that memorization in DMs is driven by a fundamental $n$-dependent change in the loss landscape -- not by a sample repetition during training.

Effect of the model capacity. To study more precisely the role of the model capacity on the memorization--generalization transition, we vary the number of parameters $p$ by changing the U-Nets base width $W \in \{8, 16, 32, 48, 64\}$, resulting in a total of $p\in\{0.26, 1, 4, 9, 16\}\times10^6$ parameters. In the left panel of Figure 3, we plot both the FID (top row) and the normalized memorization fraction (bottom row) as functions of $\tau$ for several width $W$ and training set sizes $n$. Panels A and C demonstrate that higher-capacity networks (larger $W$) achieve high-quality generation and begin to memorize earlier than smaller ones. Panels B and D show that the two characteristic timescales simply scale as $\tau_\mathrm{gen} \propto W^{-1}$ and $\tau_\mathrm{mem} \propto nW^{-1}$. In particular, this implies that, for $W>8$, the critical training set size $n_\mathrm{gm}(p)$ at which $\tau_\mathrm{mem}=\tau_\mathrm{gen}$ is approximately independent of $p$ (at least on the limited values of $p$ we focused on). When $n>n_\mathrm{gm}(p)$, the interval $\left[\tau_\mathrm{gen}, \tau_\mathrm{mem}\right]$ opens up, so that early stopping within this window yields high quality samples without memorization. In the right panel of Figure 3, we display this boundary (solid line) in the $(n, p)$ plane by fixing the training time to $\tau=\tau_\mathrm{gen}$, that we identify numerically using the collapse of all FIDs at around $W\tau_\mathrm{gen}\approx 3\times 10^6$ (see panel B), and computing the smallest $n$ such that $f_\mathrm{mem}(\tau)=0$. The resulting solid curve delineates two regimes: below the curve, memorization already starts at $\tau_\mathrm{gen}$; above the curve, the models generalize perfectly under early stopping. We repeat this experiment for $\tau=3\tau_\mathrm{gen}$ and $\tau=8\tau_\mathrm{gen}$, showing saturation to larger and larger $p$ as $\tau$ increases. Eventually, for $\tau \to \infty$, we expect these successive boundaries to converge to the architectural regularization threshold $n^\star(p)$, i.e. the point beyond which the network avoids memorization because it is not expressive enough, as found in [17] and highlighted in the right panel of Figure 1. In order to estimate $n^\star(p)$, we measure for a given $\tau$ the largest $n(\tau)$ yielding $f_\mathrm{mem}\approx0$. The curve $n(\tau)$ approaches $n^\star(p)$ for large $\tau$. We therefore estimate $n^\star(p)$ by measuring the asymptotic values of $n(\tau)$, which in practice is reached already at $\tau=\tau_\mathrm{max}=2$ M updates for the values of $W$ we focus on.

Notations. We use bold symbols for vectors and matrices. The $L^2$ norm of a vector $\bm{\mathrm{x}}$ is denoted by $\lVert \bm{\mathrm{x}} \rVert = (\sum_i \bm{\mathrm{x}}_i^2)^{1/2}$. We write $f = \mathcal{O}(g)$ to mean that in the limit $n, p \to \infty$, there exists a constant $C$ such that $\lvert f \rvert \leq C \lvert g \rvert$. Setting. We study analytically a model introduced in [17], where the data lie in $d$ dimensions. We parametrize the score with a Random Features Neural Network [RFNN, 31]

An RFNN, illustrated in Figure 4 (left), is a two-layer neural-network whose first layer weights ($\bm{\mathrm{W}}\in\mathbb{R}^{p\times d}$) are drawn from a Gaussian distribution and remain frozen while the second layer weights ($\bm{\mathrm{A}} \in \mathbb{R}^{d\times p}$) are learned during training. This model has already served as theoretical framework for studying several behaviors of deep neural network such as the double descent phenomenon [43, 44]. $\sigma$ is an element-wise non-linear activation function. We consider a training set of $n$ iid samples $\bm{\mathrm{x}}^\nu \sim P_{\bm{\mathrm{x}}}$ for $\nu = 1, \ldots, n$ and we focus on the high-dimensional limit $d, p, n\rightarrow\infty$ with the ratios $\psi_p=p/d, \psi_n=n/d$ kept fixed. We study the training dynamics associated to the minimization of the empirical score matching loss defined in (Equation 4) at a fixed diffusion time $t$. This is a simplification compared to practical methods, which use a single model for all $t$. It has been already studied in previous theoretical works [28, 17]. The loss (Equation 4) is rescaled by a factor $1/d$ in order to ensure a finite limit at large $d$. We also study the evolution of the test loss evaluated on test points and the distance to the exact score $\bm{\mathrm{s}}(\bm{\mathrm{x}})=\nabla\log P_{\bm{\mathrm{x}}}$,

where the expectations $\mathbb{E}_{\bm{\mathrm{x}}, \bm{\xi}}$ are computed over $\bm{\mathrm{x}}\sim P_{\bm{\mathrm{x}}}$ and $\bm{\xi} \sim \mathcal{N}(0, \bm{I}_d)$. The generalization loss, defined as $\mathcal{L}_\mathrm{gen} = \mathcal{L}_\mathrm{test} - \mathcal{L}_\mathrm{train}$, indicates the degree of overfitting in the model while the distance to the exact score $\mathcal{E}_\mathrm{score}$ measures the quality of the generation as it is an upper bound on the Kullback–Leibler divergence between the target and generated distributions [45, 46]. The weights $\bm{\mathrm{A}}$ are updated via gradient descent

where $\eta$ is the learning rate. In the high-dimensional limit, as the learning rate $\eta \to 0$, and after rescaling time as $\tau = k\eta / d^2$, the discrete-time dynamics converges to the following continuous-time gradient flow:

with

Assumptions. For our analytical results to hold, we make the following mathematical assumptions which are standard when studying Random Features [47, 48, 49] namely (i) the activation function $\sigma$ admits a Hermite polynomial expansion $\sigma(x)=\sum_{s=0}^\infty\frac{\alpha_s}{s!}He_s(x)$; and (ii) the data distribution $P_{\bm{\mathrm{x}}}$ has sub-Gaussian tails and a covariance $\boldsymbol{\Sigma}=\mathbb{E}_{P_{\bm{\mathrm{x}}}}[\bm{\mathrm{x}} \bm{\mathrm{x}}^T]$ with bounded spectrum. We assume that the empirical distribution of eigenvalues of $\boldsymbol{\Sigma}$ converges weakly in the high dimensional limit to a deterministic density $\rho_{\boldsymbol{\Sigma}}(\lambda)$ and that $\operatorname{Tr}(\boldsymbol{\Sigma})/d$ converges to a finite limit (for a more precise mathematical statement see SM Appendix C.3). Moreover, we make additional assumptions that are not essential to the proofs but which simplify the analysis: (iii) the activation function $\sigma$ verifies $\mu_0=\mathbb{E}_z[\sigma(z)]=0$; and (iv) the second layer $\bm{\mathrm{A}}$ is initialized with zero weights $\bm{\mathrm{A}}(\tau=0)=0$. In numerical applications, unless specified, we use $\sigma(z)=\tanh(z)$ and $P_{\bm{\mathrm{x}}}=\mathcal{N}(0, \bm{I}_d)$.

Emergence of the two timescales during training. We first show in Figure 5 that the behavior of training and test losses in the RF model mirrors the one found in realistic cases in Section 2, with a separation of timescales $\tau_\mathrm{gen}$ and $\tau_\mathrm{mem}$ which increases with $n$. Equation (Equation 10) is linear in $\bm{\mathrm{A}}$ and hence it can be solved exactly (see SM). The timescales of the training dynamics are given by the inverse eigenvalues of the $p\times p$ matrix $\Delta_t \bm{\mathrm{U}}/\psi_p$. Building on the Gaussian Equivalence Principle [GEP, 51, 48, 52] and the theory of linear pencils [53], [17] ([17]) derive a coupled system of equations characterizing the Stieltjes transform of the eigenvalue density $\rho(\lambda)$ of $\bm{\mathrm{U}}$ for isotropic Gaussian data that lie in a $D$-dimensional subspace with $D\le d$ and $D=\mathcal{O}(d)$. We offer an alternative derivation presented in SM for general variance using the replica method [54] -- a heuristic method from the statistical physics of disordered systems -- yielding the more compact formulation for obtaining the spectrum stated in Theorem 1. Before stating the theorem, we introduce

where $\sigma_{\bm{\mathrm{x}}}^2=\frac{\operatorname{Tr}(\boldsymbol{\Sigma})}{d}$, $\Gamma_t=e^{-2t}\sigma_{\bm{\mathrm{x}}}^2+\Delta_t=1+e^{-2t}(\sigma_{\bm{\mathrm{x}}}^2-1)$ and the expectation is over the $u, v, w$ random variables which are independent standard Gaussian $\mathcal{N}(0, 1)$.

Let $q(z)=\frac{1}{p} \operatorname{Tr}(\bm{\mathrm{U}}-z \bm{I}_p)^{-1}$, $r(z)=\frac{1}{p} \operatorname{Tr}(\boldsymbol{\Sigma}^{1/2} \bm{\mathrm{W}}^T(\bm{\mathrm{U}}-z \bm{I}_p)^{-1} \bm{\mathrm{W}} \boldsymbol{\Sigma}^{1/2})$ and $s(z)=\frac{1}{p} \operatorname{Tr}(\bm{\mathrm{W}}^T(\bm{\mathrm{U}}-z \bm{I}_p)^{-1} \bm{\mathrm{W}})$, with $z\in\mathbb{C}$. Let

Then $q(z), r(z)$ and $s(z)$ satisfy the following set of three equations:

The eigenvalue distribution of $\bm{\mathrm{U}}$, $\rho(\lambda)$, can then be obtained using the Sokhotski–Plemelj inversion formula $\rho(\lambda)=\underset{\varepsilon \rightarrow0^+}{\lim}\frac{1}{\pi}\operatorname{Im}q(\lambda+i\varepsilon)$.

We now focus on the asymptotic regime $\psi_p, \psi_n \gg 1$, typical for strongly over‑parameterized models trained on large data sets. In this limit, the spectrum of $\bm{\mathrm{U}}$ can be described analytically by the following Theorem 2.

Let $\rho$ denote the spectral density of $\bm{\mathrm{U}}$.

where $\rho_1$ is an atomless measure with support

and $\rho_2$ coincides with the asymptotic eigenvalue bulk density of the population covariance $\tilde{\bm{\mathrm{U}}}=\mathbb{E}_{\bm{\mathrm{X}}}[\bm{\mathrm{U}}]$; $\rho_2$ is independent of $\psi_n$ and its support is on the scale $\psi_p$.

The eigenvectors associated with $\delta(\lambda-{s_t^2})$ leave both training and test losses unchanged and are therefore irrelevant. In the limit $\psi_p\gg \psi_n$, the supports of $\rho_1$ and $\rho_2$ are respectively on the scales $\psi_p/\psi_n$ and $\psi_p$, i.e. they are well separated.

The proofs of both theorems are shown in SM (Appendix C). We recall that training timescales are directly related to eigenvalues $\lambda$ via the relation $\tau^{-1}= \psi_p / \Delta_t \lambda_\mathrm{\min}$ . Theorem 2 therefore demonstrates the emergence of the two training timescales $\tau_{\mathrm{mem}}$ and $\tau_{\mathrm{gen}}$ in the overparametrized regime of the RFNN model. They are respectively associated to the measures $\rho_1$ and $\rho_2$, which are well separated in regime I, for $\psi_p\gg\psi_n\gg1$, as shown in Figure 4.

Generalization: The timescale $\tau_{\mathrm{gen}}$ on which the first relaxation takes place is associated to the formation of the generalization regime. It is related to the bulk $\rho_2$ and is or order $1/\Delta_t$. This regime only depends on the population covariance $\boldsymbol{\Sigma}$ of the data and is independent of the specific realization of the dataset. On this timescale, which is of order one, both the training $\mathcal{L}_{\mathrm{train}}$ and test $\mathcal{L}_{\mathrm{test}}$ losses decrease. The generalization loss $\mathcal{L}_\mathrm{gen} = \mathcal{L}_{\mathrm{test}}-\mathcal{L}_{\mathrm{train}}$ is zero, and $\mathcal{E}_\mathrm{score}$ tends to a value that we find to scale as $\mathcal{O}(\psi_n^{-\eta})$ with $\eta\simeq0.59$ numerically (see Figure 5).

Memorization: The timescale $\tau_{\mathrm{mem}}$, on which the second stage of the dynamics takes place, is associated to overfitting and memorization. It is related to the bulk $\rho_1$, and scales as $\psi_p / \Delta_t \lambda_\mathrm{\min}$, where $\lambda_\mathrm{\min}$ is the left edge of $\rho_1$ . In the overparameterized regime $p \gg n$, $\tau_{\mathrm{mem}}$ becomes large and of order $\psi_n/\Delta_t$, thus implying a scaling of $\tau_{\mathrm{mem}}$ with $n$. On this timescale, the training loss decreases while the test loss increases, converging to their respective asymptotic values as computed in [17]. Figure 5 indeed shows that all training and test curves separate, correspondingly the generalization loss $\mathcal{L}_{\mathrm{gen}}$ increases, at a time that scales with $\psi_p /\Delta_t \lambda_{\min}$, as shown in the inset.

As $n$ increases, the asymptotic ($\tau \rightarrow \infty$) generalization loss $\mathcal{L}_{\mathrm{gen}}$ decreases, indicating a reduced overfitting. For $n>n^*(p) = p$, although some overfitting remains (i.e., $\mathcal{L}_{\mathrm{gen}} > 0$), the value of $\mathcal{L}_{\mathrm{gen}}$ is sensibly reduced, and the model is no longer expressive enough to memorize the training data, as shown in [17]. This regime corresponds to the Architectural Regularization phase in Figure 1. We show in Figure 5 (panel C) how the generalization loss $\mathcal{L}_{\mathrm{gen}}$ varies in the $(n, p)$ plane depending on the time $\tau$ at which training is stopped. In agreement with the above results, we find that the generalization--memorization transition line depends on $\tau$ and moves upward for larger values of $\tau$, similarly to the numerical results exposed in Figure 3 and the illustration in Figure 1.

We have shown that the training dynamics of neural network-based score functions display a form of implicit regularization that prevents memorization even in highly overparameterized diffusion models. Specifically, we have identified two well-separated timescales in the learning: $\tau_\mathrm{gen}$, at which models begins to generate high-quality, novel samples, and $\tau_\mathrm{mem}$, beyond which they start to memorize the training data. The gap between these timescales grows with the size of the training set, leading to a broad window where early stopped models generate novel samples of high-quality. We have demonstrated that this phenomenon happens in realistic settings, for controlled synthetic data, and in analytically tractable models. Although our analysis focuses on DMs, the underlying score‑learning mechanism we uncover is common to all score‑based generative models such as stochastic interpolants [55] or flow matching [10]; we therefore expect our results to generalize to this broader class.

Limitations and future works. - While we derived our results under SGD optimization, most DMs are trained in practice with Adam [38]. In SM Sects. Appendix A and Appendix D, we show that the two key timescales still arise using Adam, although with much fewer optimization steps. Studying how different optimizers shift these timescales would be valuable for practical usage.

The authors thank Valentin De Bortoli for initial motivating discussions on memorization--generalization transitions. RU thanks Beatrice Achilli, Jérome Garnier-Brun, Carlo Lucibello and Enrico Ventura for insightful discussions. RU is grateful to Bocconi University for its hospitality during his stay, during which part of this work was conducted. This work was performed using HPC resources from GENCI-IDRIS (Grant 2025-AD011016319). GB acknowledges support from the French government under the management of the Agence Nationale PR[AI] RIE-PSAI (ANR-23-IACL-0008). MM acknowledges the support of the PNRR-PE-AI FAIR project funded by the NextGeneration EU program. After completing this work, we became aware that A. Favero, A. Sclocchi, and M. Wyart [58] had also been investigating the memorization--generalization transition from a similar perspective.

This document provides detailed derivations and additional experiments supporting the main text (MT). In Appendix A, we give details about the numerical experiments carried out in Section 2. In Appendix B we provide additional numerical experiments on simplified score and data models. Appendix C gives formal proofs of the main theorems of Section 3. Finally, Appendix D exposes more details on the numerical experiments of Section 3.

Dataset. All numerical experiments in Section 2 of the MT use the CelebA face dataset [22]. We center-crop each RGB image to $32\times 32$ pixels and convert to grayscale images in order to accelerate the training of our Diffusion Models (DMs). To precisely control the samples seen by a model, no data augmentation is applied, and we vary the training set size $n$ in the window $\left[128, 32768\right]$. Examples of training samples are shown in the left-most block of Figure 6.

Architecture. As commonly done in DDPMs implementations [e.g., 2, 4], the network approximating the score function is a U-Net [21] made of three resolution levels, each containing two residual blocks with channel multipliers $\{1, 2, 4\}$ respectively. We apply attention to the two coarsest resolutions, and embed the diffusion time via sinusoidal position embedding [40]. The base channel width $W$ varies from $16$ to $64$ depending on the experiment, resulting in a total of $1$ to $16$ million trainable parameters.

Time reparameterization. Compared to the framework presented in the MT, the DDPMs we train make use of a time reparameterization of the forward and backward processes with a variance schedule $\{\beta_{t'}\}_{t'=1}^T$, where $T$ is the time horizon given as a number of steps, fixed to 1000 in our experiments. The variance is evolving linearly from $\beta_1 = 10^{-4}$ to $\beta_{1000} = 2 \times 10^{-2}$. A sample at time $t'$, denoted $\bm{\mathrm{x}}(t')$, can be expressed from $\bm{\mathrm{x}}(0)$ as the following interpolation

where $\overline{\alpha}(t') = \prod_{s=1}^{t'} (1-\beta_s)$, and $\bm{\xi}$ is a standard and centered Gaussian noise. This is a reparameterization of the Ornstein-Uhlenbeck process from Equation 1 defined through time $t$ in the MT, with

Training. All DMs are trained with Stochastic Gradient Descent (SGD) at fixed learning rate $\eta = 0.01$, fixed momentum $\beta=0.95$ and batch size $B=\min(n, 512)$. We focus on SGD to facilitate the analysis of time scaling, avoiding problems that may cause alternative adaptive optimization schemes like Adam [38]. We train each model for at least 2M SGD steps, sometimes more for large values of $n$ displaying memorization only later. We do not employ exponential moving average or learning-rate warm-up.

Generation. To accelerate sampling while preserving FID, we employ the DDIM sampler of [41] ([41]) which replaces the Markovian reverse SDE with a deterministic, non-Markovian update. Given a trained denoiser $\bm{\xi}_{\bm{\theta}}(\bm{\mathrm{x}}_t, t)$, we iterate for $t=T', \ldots, 1$

with $T'=200$. During training, we generate at 40 milestones a set of 10, 000 samples to assess generalization and memorization. Examples of samples obtained from a model trained on $n=1024$ samples with base width $W=32$ are shown in the middle and right blocks from Figure 6 for two training times, $\tau=190$ K and $\tau=1.62$ M. At $\tau=190$ K the model generalizes ($f_\mathrm{mem}=0\%$) and achieve a test FID of 35.1. After too much training, memorization sets in and, by $\tau=1.62$ M steps, nearly half the generated samples reproduce training images ($f_\mathrm{mem}=47.2\%$).

Statistical evaluation. FIDs [42] are computed¹ using 10, 000 generated samples and 10, 000 test samples, averaged over 5 independent runs with disjoint test sets. Error bars in the MT denote twice the standard deviation. Training and test losses are estimated similarly over 5 repeated evaluation on $n$ training samples and 2048 test samples, and give negligible confidence intervals. For the memorization fraction $f_\mathrm{mem}(\tau)$, we report the standard error on the mean obtained via bootstrap resampling of the 10, 000 generated samples. We also verified that the scaling in the memorization time $\tau_\mathrm{mem}$ is insensitive to the choice of the threshold $k$ used to define $f_\mathrm{mem}$ in Equation 6 by testing larger and lower values.

Computing resources. Most trainings were performed on Nvidia H100 GPUs (80GB of memory). A typical run of 2M steps takes approximately 50 hours on two GPUs and vary with the model size (defined through its base width $W$). In total, we train 18 distinct models for the several $n, W$ configurations of the MT. The longest training ($n=32768$ and $W=32$ in Figure 2) ran for 11M steps. The generation of $10, 000$ samples over 40 training times takes around an additional hour per model on the same hardware support.

All the experiments in the MT use a fixed batch size $B=512$, and in Sect Section 2 we emphasize that the observed $\mathcal{O}(n)$ scaling of $\tau_\mathrm{mem}$ cannot be explained by repetition over training samples. To validate this statement, the left panel of Figure 7 shows FID and memorization fraction curves when we train the models with full-batch updates ($B=n$) for $n\in\left[128, 512\right]$. At any fixed $\tau$, every sample has been seen exactly $\tau$ times. Yet $\tau_\mathrm{mem}$ continues to grow linearly with $n$, as shown in the inset. This demonstrates that larger datasets reshape the loss landscape -- requiring proportionally more updates to overfit -- rather than simply increasing memorization through repeated exposure of training samples.

We conclude this section by repeating our analysis at fixed $W=64$ using the Adam optimizer [38] instead of SGD with momentum. The learning rate is $\eta = 1\times10^{-4}$, gradient averages take values $(\beta_1, \beta_2)=(0.9, 0.999)$, and batch size $B=\min(512, n)$. We keep all other settings and evaluation metrics as above. As shown in the right panel of Figure 7, Adam yields the same two-phase training dynamics with first a generalization regime with $f_\mathrm{mem}=0$ and good performances (small FID), and later a memorization phase at $\tau_\mathrm{mem}\propto n$, as shown in the inset. The only difference is that both $\tau_\mathrm{gen}$ and $\tau_\mathrm{mem}$ occur after much fewer steps compared to SGD. This also points out that the emergence of the two well-separated timescales and their scaling is a fundamental property of the loss landscape.

The aim of this section is to show our results hold for other data distributions than natural images, and alternative score model that U-Net architectures.

Data distribution. We focus on data iid sampled from a $d$-dimensional Gaussian Mixture Model (GMM) made of two balanced Gaussians centered on $\pm \bm{\mu}$ with unit covariance, i.e.,

In what follows, we choose to work with $\bm{\mu} = \bm{1}_d$, with $\bm{1}_d = \left[1, \ldots, 1\right] ^{\mkern-1.5mu\mathsf{T}} \in \mathbb{R}^d$. In this controlled setup, the generalization score can be computed analytically from $\mathbb{P}_0$ and reads

Score model. The denoise $\bm{\xi}_{\bm{\theta}}(\bm{\mathrm{x}}_t, t)$ is implemented as a lightweight residual multi-layer neural network (see Figure 8): an input layer projecting $\mathbb{R}^d \to \mathbb{R}^W$, followed by three identical residual blocks and an output layer projecting back to $\mathbb{R}^d$. Each block consists of two fully connected layers of width $W$, a skip connection, and a layer normalization. We encode the diffusion time $t$ via sinusoidal position embedding and add it to the first feature of each block. The total number of parameter in the network is $p(d, W) = W(2d + 13) + d + 6W^2$. For $d=8$, and $W=128$, the reference setting of this section, this yields $p=102, 024$ trainable parameters.

Training and computing resources. Unless otherwise specified, we train every model of this section with SGD at fixed learning rate $\eta = 6\times 10^{—3}$ and momentum $\beta=0.95$ using full-batch updates $B=n$ for $n\in\{128, 256, 512, 1024, 2048, 4096\}$, running for up to 4M updates. All experiments are executed on an Nvidia RTX 2080 Ti, with the largest $n=4096$ requiring around 10 hours to complete.

Generalization and memorization metrics. In addition to the memorization fraction $f_\mathrm{mem}(\tau)$, we exploit this controlled setting where we know the true data distribution $\mathbb{P}_0$ to directly measure how closely it matches the generated distribution $\mathbb{P}_{\bm{\theta}}$ via the Kullback-Leibler (KL) divergence

The cross-entropy term $\mathbb{E}_{\mathbb{P}_{\bm{\theta}}}\left[\log\mathbb{P}_0\right]$ is easy to estimate using Monte Carlo,

where $\{\tilde{\bm{\mathrm{x}}}_\mu\}_{\mu=1}^N$ are $N=10, 000$ samples drawn from the model with parameters $\bm{\theta}(\tau)$ at training time $\tau$. Estimating the negative entropy term $\mathbb{E}_{\mathbb{P}_{\bm{\theta}}}\left[\log\mathbb{P}_{\bm{\theta}}\right]$ is more challenging, since DMs only give access to the score function $\bm{\mathrm{s}}_{\bm{\theta}}(\bm{\mathrm{x}}, t) = \nabla_{\bm{\mathrm{x}}} \log \mathbb{P}_{\bm{\theta}}(\bm{\mathrm{x}})$ and not the underlying probability distribution $\mathbb{P}_{\bm{\theta}}$. We can however employ time integration to express it as a function of the score only,

with

This expression assumes that the model learns an accurate representation of the score function. It is noteworthy to mention that samples are generated using standard Euler-Maruyama discretization of the backward process Equation 2 of the MT over $T=1000$ timesteps.

In Figure 9, the left panel shows how the KL divergence and memorization fraction evolve with training time $\tau$ for different training set sizes $n$ at fixed width $W=128$, while the right panel fixes $n=2048$ and varies $W$. In both cases, we observe two distinct phases. First, the KL divergence decreases to near zero on a timescale $\tau_\mathrm{gen}$ independent of $n$ during which the model fully generalizes ($f_\mathrm{mem}=0$). Beyond $\tau_\mathrm{gen}$, both $D_\mathrm{KL}(\mathbb{P}_{\bm{\theta}}|\mathbb{P}_0)$ and $f_\mathrm{mem}$ begin to rise at a time $\tau_\mathrm{mem}$ that scales linearly with $n$, as highlighted by the inset of the left panel. In contrast, $\tau_\mathrm{mem}$ scales with $W^{-1}$, as shown in the inset of the right panel. While the precise dependence of $\tau_\mathrm{gen}$ with $W$ remains inconclusive in this setting and require a more careful analysis, these results on the GMM mirror the main findings of the MT: the training dynamics of DMs unfolds first in a generalization phase and only later -- at $\tau_\mathrm{mem} \propto n/W$ -- memorization begins.

Conditional generation aims to sample from distributions of the form $p(\bm{\mathrm{x}} | \bm{\mathrm{y}})$, where $\bm{\mathrm{y}}$ denotes a conditioning variable such as a class label, a text embedding, or any other contextual information. DMs can naturally be extended to this setting using for instance classifier-free guidance [56]. Although conditioning often improves sample quality in practice, memorization effects have also been observed in models trained conditionally [25, 59, 60]. Our analysis does not rely on the model being unconditional since these variables typically enter the model as additional inputs and we expect our result to hold in this setting as well. To investigate it, we train a classifier-free guidance model to generate sample from our Gaussian mixture conditionally on the class label, and compute the memorized fraction as a function of $\tau$ that we report in Figure 10. In the inset, when rescaling the training time by $n$, the curves for $n \in \{256, 512, 1024\}$ all collapse perfectly, confirming that the phenomenon persists in the conditional setting. For more complex datasets, $\tau_\mathrm{mem}$ and $\tau_\mathrm{gen}$ may in fact depend on the conditioning variable and intermediate regimes could exist where certain classes have already entered the generalization (or memorization) phase while others have not yet.

In the following we provide the mathematical arguments and the proofs for the statement in the MT. The section using the replica method is not mathematically rigorous but uses a well established method of theoretical physics, which has been already shown to provide correct results in several cases. The final result is rigorous, since it can alternatively be obtained from the rigorous free random matrix approach of [17], as shown in Appendix C.4.

We recall here the notations used throughout Section 3 of the MT and Appendix C of the SM.

Unless specified, all the expectation values are taken for standard Gaussian variables. We will denote

the matrix whose columns are the data point vectors and likewise we decompose $\bm{\mathrm{W}}$ as

where $\bm{\mathrm{W}}_i \in \mathbb{R}^{ d}$ denotes the $i$ th row of $\bm{\mathrm{W}}$.

We recall the definitions of the matrices $\bm{\mathrm{U}}$ and $\bm{\mathrm{V}}$