1 Introduction↩︎

Diffusion models are a popular class of deep generative models [1]–[3]. Their key insight is to reframe generation as the reversal of a noising process, characterized by a score function. Estimation of this score function turns out to be equivalent to a denoising problem that can be solved by training a neural network on a simple regression-like objective [4]–[6]. Despite their conceptual simplicity, they are capable of achieving state-of-the-art performance in a wide range of applications [7]–[12], and appear in many of the most widely deployed generative AI systems [13].

The widespread success of diffusion models naturally raises questions about the mechanisms underpinning their strengths, motivating a growing body of theoretical work on diffusion models [14]–[17]. Of particular interest is the question of generalization [18]–[20]: understanding when a score network trained on finite data faithfully represents the underlying distribution, rather than memorizing the training set. This question is further sharpened by practical concerns around data privacy [21], making a principled understanding of the generalization-memorization tradeoff both scientifically and practically important.

Alongside these developments, an extensive body of theory has sought to explain how overparameterized neural networks generalize. These models famously challenge the classical bias–variance tradeoff by generalizing well despite perfectly memorizing noisy training data [22], [23]—a phenomenon known as benign overfitting [24]. This is explained by theoretical analysis showing that generalization is maintained by a form of self-induced regularization that occurs in the overparameterized regime. Furthermore, when taken as a function of model complexity, the test error exhibits a double descent [25]: it first decreases, peaks near the interpolation threshold, and then decreases again as model size grows (see 1, left). This behavior was first identified empirically in high-dimensional problems [26] and then theoretically grounded through the analysis of random feature models [27]. Together, these theories form a now widely internalized picture: in highly overparameterized models, fitting the training data perfectly is not only compatible with generalization but, in the right regime, actively beneficial.

One of the benefits of the diffusion model framework is that it reuses much of the previously existing machinery of deep learning, using existing neural network architectures, training with standard optimization algorithms (e.g., ADAM [28]), and fitting by minimizing a regression-like objective [2], [5]. As a result, it is tempting to conclude that the previously mentioned deep learning phenomena are at work, and benign overfitting and double descent carry over more or less intact. In this work, we show that such an assumption is largely incorrect. While several recent works have hinted that diffusion models may be incapable of benign overfitting [20], [29], [30] and other recent works have explored the relationship between over-training and memorization [19], [31], the claim that benign overfitting and double descent do not occur in diffusion models is yet to receive rigorous empirical or theoretical investigation.

For diffusion models, we show that in most practical settings, it is fundamentally impossible to have both the train and test losses simultaneously small, invalidating benign overfitting. We further show that, instead of a double descent curve, the excess risk forms a U-shaped curve with increasing complexity, akin to the classical bias–variance tradeoff. This is illustrated in the random feature network experiment in 1 where the excess risk forms a U-shape that tapers off as the width of the network grows large. In 2, we see how this phenomenon affects high-dimensional image generation problems based on a U-Net architecture [1], [32], showing that as the number of U-Net features is increased, the over-trained diffusion model goes from low quality generation, to generalization, and then towards exactly memorizing the train set. Through a combination of fundamental theoretical results, simplified models, and high-dimensional experiments, we explore this phenomenon, the mechanisms underpinning it and how diffusion models generalize in spite of it.

Concretely, we present the following contributions:

• In 3, we derive fundamental negative results showing that the population and empirical score matching loss cannot be simultaneously small in practical settings, unless the sample size grows exponentially with data dimension, regardless of the score network architecture.

• In 4, through the analysis of a linear random feature model, we provide a precise understanding of the interplay between model complexity, overfitting and noise scale. By comparing with a classical regression setup, we illustrate the key differences that prevent benign overfitting from occurring in score matching.

• In 5, we continue the analysis of 4, identifying key components of the score matching framework that can act as implicit regularization and help prevent overfitting in practice. More precisely, we identify the role of time-smoothness of the score and discuss the role of early stopping during training.

• We validate our theoretical findings by experiments conducted in high-dimensional image generation settings.

Notations. For probability measures \(\mu\) and \(\nu\) in \({\mathbb{R}^d}\), if \(\mu\) is absolutely continuous with respect to \(\nu\), denoted by \(\mu \ll \nu\), we denote \(\mathrm{d}\mu/\mathrm{d}\nu\) its density. If \(\mu \ll \nu\), the Kullback-Leibler (KL) divergence is \(\mathrm{KL}\left(\mu || \nu \right) = \int \log (\mathrm{d}\mu / \mathrm{d}\nu) \mathrm{d}\mu\), the relative Fisher information is \(\mathscr{J}(\mu || \nu) = \int \Vert \nabla \log (\mathrm{d}\mu / \mathrm{d}\nu) \Vert^2 \mathrm{d}\mu\), the Fisher information is \(\mathscr{J}(\nu) = \int \Vert \nabla \log(\mathrm{d}\nu / \mathrm{d}x) \Vert^2 \mathrm{d}\nu\). For \(a,b \in \mathbb{R}\), \(a \wedge b = \min(a,b)\), \(a \vee b = \max(a,b)\), and \(a_+ = a \vee 0\).

2 Background↩︎

By learning to reverse a noising process, diffusion models approximate a distribution \(\nu\) on \(\mathbb{R}^d\) from a dataset sampled from \(\nu\). We briefly review this construction and fix the notation used in the paper.

Noising process. The noising process is given by an SDE of Ornstein–Uhlenbeck-type in \({\mathbb{R}^d}\), \[\label{eq:forward-process-sde} \mathrm{d}{X_t}= -\kappa {X_t}\mathrm{d}t + \sqrt{2} \mathrm{d}\mathrm{B}_t\;, \qquad {X}_0 \sim \nu\;,\tag{1}\] where \(\kappa \geq 0\) and \(\mathrm{B}_t\) is a standard \(d\)-dimensional Brownian motion. We assume that \(\mathbb{E}_{X\sim\nu} [ \Vert X \Vert^2 ] < \infty\). It is well-known that \({X_t}|{X}_0 \sim \mathrm{N}(\alpha_t {X}_0, \sigma_t^2 {\mathrm{I}}_d)\), where the mean and variance are given by \[\label{eq:cond95score95form} \alpha_t = e^{- \kappa t} \;, \qquad \sigma^2_t = \kappa^{-1} (1 - \alpha_t^2) \;,\tag{2}\] with the convention \(\sigma_t^2 = 2t\) when \(\kappa=0\). Fixing a time horizon \(T > 0\) and writing \({Y}_t := {X}_{T-t}\), the classical result of [33], [34] identifies the law of the time-reversed process as the solution of \[\mathrm{d}{Y}_t = \kappa {Y}_t \mathrm{d}t + 2 \nabla \log {p}_{T-t}({Y}_t) \mathrm{d}t + \sqrt{2} \mathrm{d}\mathrm{B}_t\;, \qquad {Y}_0 \sim {p}_T \;,\] where \(p_t\) denotes the density of \(X_t\). Thus, if we could simulate \({Y}_t\) starting from \({Y}_0 \sim {p}_T\), we would obtain samples \(Y_T \sim \nu\). The only unknown quantity is the score function \(\nabla \log {p}_t\), which must be learned from data. Once we have an approximate score function \(\hat{s}:{\mathbb{R}^d}\times [0,T] \to {\mathbb{R}^d}\), we generate samples by numerically simulating \[\mathrm{d}\widehat{{Y}}_t = \kappa \widehat{{Y}}_t \mathrm{d}t + 2 \hat{s}(\widehat{{Y}}_t, t) \mathrm{d}t + \sqrt{2} \mathrm{d}\mathrm{B}_t\;, \qquad {Y}_0 \sim q_0 \;.\] For \(\kappa > 0\), we set \(q_0=\mathrm{N}(0, \kappa^{-1} {\mathrm{I}}_d)\) as it approximates \(p_T\) for large \(T\), and for \(\kappa = 0\), we set \(q_0 = \mathrm{N}(0, 2 T {\mathrm{I}}_d)\). We denote by \(q_t\) the probability density of \(\widehat{{Y}}_t\) for \(t \in [0, T]\).

Score matching. The score function is approximated using a deep neural network. Ideally, the goal is to minimize the (population) explicit score matching (ESM) loss, \[\begin{align} {\mathscr{L}_{\mathrm{ESM}}}({s}, t) := \mathbb{E}\big [ \big \| {s}(t, {X_t}) - \nabla \log p_t({X_t}) \big \|^2 \big ]\;, \qquad {\mathscr{L}_{\mathrm{ESM}}}({s}, \varpi) = \int {\mathscr{L}_{\mathrm{ESM}}}({s}, t) \, \mathrm{d}\varpi(t)\;, \end{align}\] where \(\varpi\) is some non-negative measure on \([0, T]\). While the time-weighting \(\varpi\) can be chosen freely, a canonical choice is the evidence lower bound (ELBO) weighting, \(\varpi_{\operatorname{ELBO}}^\epsilon (\mathrm{d}t) = \mathbb{1}_{[\epsilon, T]}(t)\), for \(\epsilon > 0\). Under this weighting, the explicit score matching loss upper bounds the marginal KL divergence with \(\mathrm{KL}\left({p}_\epsilon || {{q}}_{T-\epsilon} \right) \leq {\mathscr{L}_{\mathrm{ESM}}}(s, \varpi_{\operatorname{ELBO}}^\epsilon) + \mathrm{KL}\left({p}_T || q_0 \right)\). Since the ESM loss depends on the intractable score \(\nabla \log p_t\), a tractable alternative is considered: the (population) denoising score matching (DSM) loss, \[{\mathscr{L}_{\mathrm{DSM}}}(s, \varpi) := \int \mathbb{E}[\ell_t(s, X_0)] \mathrm{d}\varpi(t)\;, \qquad \ell_t(s, X_0) := \mathbb{E}\big [ \big \| s(t, {X_t}) - \nabla \log p_{t\vert 0}({X_t}|X_0) \big \|^2 \big | X_0 \big ]\;,\] where \(p_{t|0}(\cdot|\cdot)\) is the density of \(X_t\) given \(X_0\). This loss is equal to the ESM loss up to a constant whilst only using the simpler conditional score \(\nabla \log p_{t\vert 0}\), turning the learning problem into a denoising problem [5], [6]. It is then estimated using a dataset \(S = \{Z_i\}_{i=1}^n \sim \nu^{\otimes n}\), leading to the empirical DSM loss, \[\begin{align} \label{eq:emp95dsm} {\widehat{\mathscr{L}}_{\mathrm{DSM}}}(s, \varpi) := \frac{1}{n} \sum_{i=1}^n \int \ell_t (s, Z_i) \, \mathrm{d}\varpi(t) \;, \qquad n \geq 1\;. \end{align}\tag{3}\] We also define the empirical ESM loss, related to the above by an additive data-dependent constant, \[\begin{align} \label{eq:emp95esm} {\widehat{\mathscr{L}}_{\mathrm{ESM}}}({s}, t) := \mathbb{E}\big [ \big \| {s}(t, {\widehat{X}_t}) - \nabla \log \hat{p}_t({\widehat{X}_t}) \big \|^2 \big | S \big ]\;, \qquad {\widehat{\mathscr{L}}_{\mathrm{ESM}}}({s}, \varpi) = \int {\widehat{\mathscr{L}}_{\mathrm{ESM}}}({s}, t) \, \mathrm{d}\varpi(t)\;, \end{align}\tag{4}\] where \(\mathrm{d}{\widehat{X}_t}= -\kappa {\widehat{X}_t}\mathrm{d}t + \sqrt{2} \mathrm{d}\mathrm{B}_t,\) \(\require{upgreek} {\widehat{X}}_0 \sim n^{-1} \sum_{i=1}^n \updelta_{Z_i}\) and for \(t>0\), \(\hat{p}_t\) is the density of \({\widehat{X}_t}\).

3 Fundamental limitations of overfitting in diffusion models↩︎

Benign overfitting occurs when good generalization performance is achieved despite interpolating training data. In this section, we derive impossibility results that show that this cannot happen in diffusion models without an exceedingly large sample size. Unlike [sec:fine-grainded-linear-analysis,sec:preventing-overfitting], the results presented here hold regardless of the score network architecture. In this sense, they identify fundamental limitations of the score matching framework itself, rather than of any particular class of models.

3.1 Warm-up: A Fisher information lower bound↩︎

We begin with the following simple lemma, which is a generic lower bound on the sum of the population and empirical ESM losses. The proof can be found in 8.1.1.

1 shows that for both the empirical and population score matching losses to be small, we must have the empirical and population measures \(\hat{p}_t\) and \(p_t\) close also. The measure of closeness here is the relative Fisher information, which under additional concentration assumptions, can be shown to upper bound the KL divergence (see 8.1.2). As \(\hat{p}_t\) is initialized at the empirical distribution \(\require{upgreek} n^{-1} \sum_{i=1}^n \updelta_{Z_i}\) where \((Z_1,\dots, Z_n) \sim \nu^{\otimes n}\), this suggests that benign overfitting only happens when the sample size is large enough. We make this argument precise in the next subsection.

3.2 A quantitative lower bound under time-integration↩︎

We now refine the previous result by obtaining lower bounds on the amount of data required for the empirical and population losses to be small simultaneously. The lower bound is a function of the effective dimension of the data distribution \(\nu\) which we capture through the lower Rényi dimension \(d'\), formally defined in 7.1. In cases where \(\nu\) is supported on a compact manifold, \(d'\) is exactly the manifold dimension. We also extend Lemma 1 by considering the full time-integrated loss as opposed to an individual time-slice, requiring only that the weighting \(\varpi\) has positive non-decreasing density on an interval. The bound is generic and holds for any random score function \(\hat{s}\) that can be taken as a function of training data. Given a real-valued function \(w\), we define the generalized inverse \(w^{-}(r) := \inf \{t : w(t) \geq r\}\), with the convention \(\inf \emptyset = \infty\).

The assumption on the measure \(\varpi\), encompasses several standard choices of time-weighting within the continuous training framework [2]. In the following remark, we consider several of these, obtaining closed-form bounds.

We highlight that the lower bound is increasing exponentially in the intrinsic dimension and depends on \(\epsilon\), supporting both the idea that manifold distributed data is beneficial to diffusion models and that early stopping in the reverse process can benefit generalization [15], [20], [35], [36]. Since many popular implementations of diffusion models take \(\epsilon \in \left\{ 10^{-3}, 10^{-5} \right\}\), the expectation that the amount of data scales with \(\epsilon^{-d'/2}\) is unreasonably strong. Thus, in most practical settings, diffusion models will not exhibit the property of benign overfitting.

The proof, presented in 8.1.3, uses a similar approach as 1, showing that for empirical and population losses to be small, \(\hat{p}_t\) and \(p_t\) must be close. Instead of the Fisher information, we derive a connection to the total variation via an argument based on Girsanov’s theorem.

4 Fine-grained analysis on linear RFNNs↩︎

Section 3 establishes at a general level that benign overfitting cannot occur in diffusion models in most practical settings. However, it does not reveal how overfitting manifests, nor how it depends on the interplay between model size, data dimension, and noise scale. To get a finer understanding, we now analyze a concrete and tractable model class: two-layer linear random feature networks [37], similar to the linear activation case of the models considered in [19], [38], [39]. Despite the simplicity of the linear setting, it has been central to the theoretical analysis of benign overfitting and double descent [25], [40], suggesting that it is rich enough to capture overfitting in diffusion models.

Setup. We study \(2\)-layer linear score networks of the form \(s_{W,A} : x \mapsto A W x\), where \(W \in \mathbb{R}^{p \times d}\) is fixed and \(A \in \mathbb{R}^{d \times p}\) is a trainable parameter. We fix a time \(t \in (0,T]\) and consider the minimization of the empirical risk defined by \(A \mapsto \sigma_t^2 {\widehat{\mathscr{L}}_{\mathrm{DSM}}}(s_{W,A}, t)\). We denote by \(\mathbf{\Sigma}\) the covariance matrix of the data distribution \(\nu\) and define the empirical covariance matrix \(\widehat{\mathbf{\Sigma}}:= n^{-1} \sum_{i=1}^n Z_i Z_i^T\) along with \(\mathbf{\Sigma}_t := \alpha_t^2 \mathbf{\Sigma}+ \sigma_t^2 {\mathrm{I}}_d\) and \(\widehat{\mathbf{\Sigma}}_t := \alpha_t^2 \widehat{\mathbf{\Sigma}} + \sigma_t^2 {\mathrm{I}}_d\).

In the following lemma, we characterize the DSM and ESM loss associated with the empirical risk minimization problem. The proof can be found in 8.2.

For \(p\geq d\), the ESM and DSM loss at \(\widehat{A}\) are independent of \(p\). By the proof of 2, this is because the empirical risk minimization (ERM) problem becomes equivalent to the ERM associated with a \(d\)-dimensional linear model (independent of \(p\)), similar to [30]. See 8.2 for details. We specify our setup further in the next assumption.

The next proposition is a precise characterization of overfitting in a proportional asymptotic regime.

Proposition 1 demonstrates how the population loss is a U-shaped function of model complexity, captured by the relative model width, \(\psi_p\). When \(\psi_p > \psi_d\), we observe that the empirical risk of ?? is of order \(\mathcal{O}_{} \left( \psi_d^{-1} \right)\) and can, therefore, become arbitrarily small and interpolate training data when the data dimension is large. In the same regime, we observe that the population score matching loss is of order \(\sigma_t^{-4}\) and therefore explodes as \(t\to 0^+\). This can also be observed clearly in 3, displaying the exact asymptotics of the risks, obtained in the proofs in 8.2. It demonstrates that interpolation and generalization are mutually incompatible. We also study the dependence on \(t\) in 3, showing that the rate of explosion of the excess risk as \(t \to 0^+\) depends on the model complexity.

Comparison with regression. To understand how these results differ from traditional regression, we consider a multi-output regression problem with the same data and architecture. Given a target \(B_\star \in \mathbb{R}^{d \times d}\) and noiseless labels \(Y_i := B_\star Z_i\), we minimize the empirical risk \(A \mapsto n^{-1} \sum_{i=1}^n \left\Vert s_{W, A}(Z_i) - Y_i \right\Vert^2 + \lambda\left\Vert AW \right\Vert_{\mathrm{F}}^2\). We focus on the overparameterized regime \(p > d > n\) and \(\lambda = 0\), where minimizers \(A\) satisfy \((AW - B_\star) \widehat{\mathbf{\Sigma}} = 0\). Since \(\widehat{\mathbf{\Sigma}}\) has rank \(n < d\) a.s., this gives a \(d(d-n)\)-dimensional affine subspace of solutions for \(AW\), contrasting sharply with 2, where score matching produces a unique value of \(A W = -\widehat{\mathbf{\Sigma}}_t^{-1}\) even without regularization.

The excess risk for the regression problem is given directly by \(\| (AW - B_\star) \mathbf{\Sigma}^{1/2} \|_{\mathrm{F}}^2\). The expression has a similar form as the score matching ESM loss in Lemma 2 once replacing \(B_\star\) with \(-\mathbf{\Sigma}_t^{-1}\), and \(\mathbf{\Sigma}\) with \(\mathbf{\Sigma}_t\) in the weighting. To resolve the under-determination, we consider the ridgeless limit \(\lambda \to 0^+\), recovering the minimum-norm interpolator. This projects \(AW\) on the orthogonal of \(\mathrm{ker}(\widehat{\mathbf{\Sigma}})\), leading to the excess risk \(\| B_\star \, \Pi_{\mathrm{ker}(\widehat{\mathbf{\Sigma}})} \mathbf{\Sigma}^{1/2} \|_{\mathrm{F}}^2\), where \(\Pi_{\mathrm{ker}(\widehat{\mathbf{\Sigma}})}\) is the column-wise projection onto the kernel of \(\widehat{\mathbf{\Sigma}}\). Thus, generalization can occur in the overparameterized regime due to two factors: (i) the minimum-norm bias projects \(AW\) on the orthogonal of \(\mathrm{ker}(\widehat{\mathbf{\Sigma}})\); (ii) the residual is small when \(B_\star\) aligns with \(\widehat{\mathbf{\Sigma}}^{1/2}\) within the kernel of \(\widehat{\mathbf{\Sigma}}\). This is made precise in the benign overfitting literature that identifies a form of self-induced regularization resulting from alignment [24], [40].

The picture is fundamentally different for score matching. First, since the empirical problem exactly determines the solution \(AW = - \widehat{\mathbf{\Sigma}}_t^{-1}\), there can be no minimum-norm bias that projects \(AW\) orthogonally to \(\mathrm{ker}(\widehat{\mathbf{\Sigma}})\). Moreover, whereas \(B_\star\) could reasonably align with \(\mathbf{\Sigma}^{1/2}\) on the null space in such a way that makes their product small, it is less likely that \(-\mathbf{\Sigma}_t^{-1}\) aligns in the same way with \(\mathbf{\Sigma}_t^{1/2}\). Indeed, if \(\mathbf{\Sigma}\) is not full rank, we have \(\| \mathbf{\Sigma}_t^{-1} \Pi_{\mathrm{ker}(\widehat{\mathbf{\Sigma}})} \mathbf{\Sigma}_t^{1/2} \|_{\mathrm{F}}^2 \geq (d-n)\sigma_t^{-2}\) which explodes as \(t \to 0\). This is made precise in the proof of 1, leading to the first term of ?? , corresponding to the contribution of the null space of \(\widehat{\mathbf{\Sigma}}\) to the ESM loss.

5 Preventing overfitting in diffusion models↩︎

The absence of benign overfitting in diffusion models suggests that regularization is required to achieve good generalization in the overparameterized regime. In this section, we show that the score matching framework contains forms of implicit regularization that help prevent overfitting.

5.1 Implicit regularization due to time smoothness↩︎

In most practical applications of diffusion models, a single neural network is trained across multiple noise levels. Whilst primarily done for efficiency, the resulting time-smoothness of the score function has been shown empirically to benefit generalization in diffusion models [41]. Here, we study the impact of time-smoothness of the score estimator in the setting of 4 by considering the extreme case where it is time-independent.

Setup. We parameterize the score network as in 4 with \(p > d\) and we consider a time weighting \(\varpi\) which is a probability distribution on \([0,T]\) such that we can define the probability measure \(\pi(\mathrm{d}t) \propto \sigma_t^{-2} \varpi(\mathrm{d}t)\). For any distribution \(\omega \in \left\{ \pi, \varpi \right\}\) on \([0,T]\), we write \[\begin{align} \alpha_\omega^2 := \int \alpha_t^2 \mathrm{d}\omega(t)\;, \quad \sigma_\omega^2 := \int \sigma_t^2 \mathrm{d}\omega(t)\;, \quad \mathbf{\Sigma}_\omega := \int \mathbf{\Sigma}_t \mathrm{d}\omega(t)\;, \quad \widehat{\mathbf{\Sigma}}_\omega := \int \widehat{\mathbf{\Sigma}}_t \mathrm{d}\omega(t)\;. \end{align}\] In the next lemma, we characterize the empirical risk minimization in this setting.

Next, we make this observation precise by characterizing the asymptotic score matching losses.

2 shows that time-smoothness effectively prevents the explosion of the population loss observed in 1. It also shows that the limiting empirical risk cannot be arbitrarily small as soon as \(\varpi\) is not a Dirac distribution (by similar calculations, we reach the same conclusion for \({\widehat{\mathscr{L}}_{\mathrm{ESM}}}(s_{W,\widehat{A}}, \varpi)\)). This contrasts with the fixed noise scale case of 1, even when \(\psi_d\) is large, suggesting that the time-smoothness prevents overfitting. Moreover, while in 1 the asymptotic ESM loss could explode as \(t \to 0^+\), in the present setting the ESM loss stays bounded, showing that the time-smoothness acts as a form of implicit regularization, as hinted in 2.

In the proof of 2 (see 8.3), we derive exact asymptotics for the population and empirical losses, which we report in the right figure in 3 (green curves). Using the asymptotics of the previous section, we compare it against an overfit linear score network that has independent weights at each time, corresponding to the completely non-time-smooth case (purple curves). We observe that time-smoothness increases the empirical risk while decreasing the population loss, showing that it is indeed acting as an implicit regularization mechanism.

We verify this in a practical image generation setting similar to 2, but considering a time-independent DDPM U-Net model trained across a range \([t, t+r]\) for varying \(r\). We evaluate train and test error at \(t\) to see how changing \(r\) affects the error. We present the results in [fig:no95time95exp] where we find that increasing \(r\) increases the train error at \(t\) and therefore is indeed preventing overfitting. Interestingly, the test error first decreases and then increases for large \(r\), suggesting a tradeoff on \(r\).

5.2 Implicit regularization via early stopping↩︎

Next, we study the impact of early stopping during the training of the score network. To this end, we fix a time \(t \in [0,T]\) and study the dynamics of the parameter \(A\) under gradient flow, \[\begin{align} \label{eq:gradient-flow-dynamics} \frac{\mathrm{d}}{\mathrm{d}\tau}A(\tau) = -p^{-1}\nabla_A (\sigma_t^2 {\widehat{\mathscr{L}}_{\mathrm{DSM}}}(s_{W,A(\tau)}, t)) \;, \quad \tau \geq 0 \;, \end{align}\tag{7}\] where \(W \in \mathbb{R}^{p \times d}\) is as in 4, and the \(p^{-1}\) factor is ensuring a finite limit as \(p \to \infty\) [19]. Without loss of generality, we assume that \(A(0) = 0\). We show in 8.4 that in the limit \(p \to \infty\) with \((n,d)\) fixed, the matrix \(A(\tau) W\) converges to a finite limit following the gradient flow of a linear model (as shown in 8.4). Moreover, when \(\tau \to \infty\), \(A(\tau) W \to - \widehat{\mathbf{\Sigma}}^{-1}_t\), which leads to the same score matching bounds as in 1 for \(p > d\). In the next proposition, we show that for small enough values of \(\tau\), we avoid the explosion of the test loss in 1.

This result shows that early stopping prevents the explosion of the test loss as \(t \to 0\) given that the early stopping \(\tau\) is of order \(\sigma_t^{-2}\). This suggests that, despite not being benign, overfitting and memorization appear later in training, allowing early stopping to effectively prevent it, which aligns with recent results [19], [30]. The advantages of early stopping can also be observed in high-dimensional experiments: in [fig:train95unet], we report the evolution of the test and train errors for a U-Net across iterations, observing that there is a unique optimal stopping time where the ESM loss is controlled as the number of parameters grows.

6 Related works & Conclusion↩︎

Random features analysis & linear models. While we focus on quantifying overfitting, we also derive exact asymptotics of the losses in the proof of 1, complementing prior analyses [30], [38], [39]. As discussed in 4, the two-layer linear RFNNs analysis recovers linear diffusion models as a particular case, extending the regularized linear model analysis of [30]. In particular, we recover the conclusion that overfitting is driven by the null space. Recent work also studied the learning dynamics of diffusion models [30], [42], including RFNN settings [18], [19]. 5.2 supplements these approaches by quantifying the role of early stopping in mitigating overfitting in linear networks. Overall, we provide a fine-grained characterization of overfitting in linear RFNNs and, crucially, incorporate the time-smoothness of the models with time-integrated losses (5.1), which is beyond the scope of prior work.

Diffusion model generalization. In addition to the aforementioned works, there has been a lot of recent interest in the generalization properties of diffusion models [18], [43], [44]. In particular, minimax statistical rates have been obtained [12], [14], [15], the relationship with data geometry has been explored [31], [35], [45] and, more recently, algorithm-dependent generalization bounds were proposed for diffusion models [20], [29]. Our work offers a new but complementary perspective, by suggesting that the mechanisms causing generalization in diffusion models are fundamentally different from those in classical settings.

Statistical lower bounds. Recent work analyzed memorization through the development of lower bounds [46]. In particular, 1 can be compared with [47], who showed that \(\mathscr{J}(\hat{p}_t || p_t) = {\widehat{\mathscr{L}}_{\mathrm{DSM}}}(\nabla \log p_t, t) - {\widehat{\mathscr{L}}_{\mathrm{DSM}}}(\nabla \log \hat{p}_t, t)\). This is of a different nature than the left-hand side of 5 , which is the sum of the empirical and population ESM losses for a generic score. That being said, the lower bounds on \(\mathscr{J}(\hat{p}_t || p_t)\) derived in [47] for mixture data could easily be combined with 1. Our results in 3.2 differ in that they consider a very generic setting and exploit the intrinsic dimension of the data in statistical guarantees. Note that 5 does not require any restrictive assumptions on \(\nu\), such as the log-Sobolev inequality used in [48].

6.0.0.1 Conclusion.

In this work, we showed that the phenomenon of benign overfitting, known to occur in many deep learning settings, does not occur in diffusion models. We first derive a fundamental limitation of score matching, showing that the empirical and population losses cannot be small without an exceedingly large training set. Then, through a combination of high-dimensional experiments and a fine-grained analysis in a simplified setting, we showed that the double descent phenomenon of deep learning does not occur and that the test error follows a U-shaped curve with respect to model complexity. Finally, by incorporating the time-smoothness of the score and early-stopped gradient dynamics into our experiments and linear network analysis, we demonstrated that these key components of score estimation can help prevent overfitting.

Acknowledgements↩︎

U.S. is partially supported by the French government under the management of Agence Nationale de la Recherche as part of the “Investissements d’avenir” program, reference ANR-19-P3IA-0001 (PRAIRIE 3IA Institute). T.F., B.D., and U.S. are supported by the European Research Council Starting Grant DYNASTY – 101039676. A.D. is supported by the France 2030 program with the reference ANR-25-PEIA-0001 (THEOREM project). A.D. is funded by the European Union (ERC-2022-SYG-OCEAN-101071601). Views and opinions expressed are however those of the author(s) only and do not necessarily reflect those of the European Union or the European Research Council Executive Agency. Neither the European Union nor the granting authority can be held responsible for them. A.D. is supported by Hi! Paris and Agence Nationale de la Recherche (Grant 11-LABX-0047). This work received government funding administered by the National Research Agency (ANR) under the France 2030 program "Hi! PARIS", grant number ANR-23-IACL-0005. The authors would also like to thank George Deligiannidis, Giovanni Conforti and Francis Bach for inspiring discussions. This paper was partially motivated by a mid-seminar argument with Francis and the SIERRA team initiated by a statement similar to the title of this work.

7 Additional Technical Background↩︎

7.1 Intrinsic dimensions↩︎

Our main results in 3 are expressed in terms of the Rényi dimension of a Borel measure \(\mu\) on \({\mathbb{R}^d}\). This notion was first developed by [49] and then rigorously treated in [50], [51]. It is a standard notion in the theory of chaotic dynamical systems. We provide below the definition.

7.2 Random covariance matrices↩︎

We recall here some basic random matrix theory results related to random covariance matrices and the Marchenko-Pastur theorem. We invite the reader to consult [54] for additional details.

Given a symmetric matrix \(\mathbf{\Sigma}\in \mathbb{R}^{d \times d}\), we denote its spectrum by \(\mathrm{Spec}(\mathbf{\Sigma})\) (where the eigenvalues are counted with their multiplicity). The spectral measure of \(\Sigma\) is defined by \[\require{upgreek} \begin{align} \label{eq:spectral-measure} \widehat{\mu}_{\mathbf{\Sigma}} := \frac{1}{d} \sum_{\lambda \in \mathrm{Spec}(\mathbf{\Sigma})} \updelta_\lambda \;. \end{align}\tag{8}\] We recall below the celebrated Marchenko-Pastur theorem [55]. To this end, we first define the Marchenko-Pastur distribution below.

The Marchenko-Pastur’s theorem is given below.

8 Omitted Proofs↩︎

In this section, we present the omitted proofs of the results from [sec:fundamental-limitations,sec:fine-grainded-linear-analysis,sec:preventing-overfitting].

8.1 Omitted proofs of 3↩︎

8.1.1 Proof of 1.↩︎

We present below the proof of 1.

8.1.2 Lower bound under log-Sobolev inequalities↩︎

In this subsection, we present some additional results to complement 1. Our goal is to show that, under the assumption that the data distribution satisfies a logarithmic Sobolev inequality, we can obtain a similar lower bound to 1, where the relative Fisher information term can be replaced by other divergences, such as the total variation distance.

8.1.2.1 Conventions for probabilistic inequalities.

We define the total variation distance as \[\begin{align} \mathrm{TV}(\mu,\nu) := \sup_A |\mu(A) - \nu(A)|\;. \end{align}\] We say that a Borel probability distribution \(\nu\) satisfies the log-Sobolev inequality with constant \(\rho\) if for all differentiable \(f \in \mathrm{L}^1(\nu)\), we have \[\begin{align} \operatorname{Ent}_\nu(f) \leq \frac{\rho}{2} \int \frac{\left\Vert \nabla f \right\Vert^2}{f} \mathrm{d}\nu\;. \end{align}\] Equipped with these definitions, we have the following proposition.

8.1.3 Proof of 1↩︎

Before giving the proof of 1, we provide some technical lemmas below. Recall that we denote by \(p_{t|0}(\cdot|x)\) the conditional density of \(X_t\) given \(X_0 = x\), defined by 2 .

The proof of 1 is based on the following proposition, providing a lower bound on the sample size, provided that two KL divergences associated to the population and empirical forward processes are both smaller than an absolute constant. Below, we use the following convention and notation for the total variation distance: \[\begin{align} \label{eq:tv95distance} \left\Vert \mu - \nu \right\Vert_{\mathrm{TV}} := \sup_A |\mu(A) - \nu(A)| \;, \end{align}\tag{9}\] for probability measure \(\mu\) and \(\nu\). We also recall that we denote by \(p_t\) the probability density of \(X_t\) following 1 and by \(\hat{p}_t\) the probability density of the process \(\mathrm{d}{\widehat{X}_t}= -\kappa {\widehat{X}_t}\mathrm{d}t + \sqrt{2} \mathrm{d}\mathrm{B}_t\) initialized at \(\require{upgreek} {\widehat{X}}_0 \sim \frac{1}{n} \sum_{i=1}^n \updelta_{Z_i}\), for \(t > 0\), where \(S := (Z_1, \dots,Z_n) \sim \nu^{\otimes n}\) is the dataset.

We can now present the proof of 1.

8.2 Omitted proofs of 4↩︎

In this section, we present the proofs of 2 and 1.

8.2.1 Proof of 2↩︎

We present below the proof of 2.

8.2.2 Proof of 1↩︎

We prove 1 as a corollary of the following more general result, which gives the exact asymptotic behavior of the score matching loss and the empirical risk.

We present below the proof of 1.

8.3 Omitted proofs of 5.1↩︎

In this section, we present the proofs of 5.1, related to the impact of the time integration in the score matching loss.

8.3.1 Proof of 3↩︎

8.3.2 Proof of 2↩︎

We provide below the proof of 2. We first provide, in the next proposition, the exact asymptotics of the score matching losses in the case of time-integrated losses.

We present below the proof of 2.

8.4 Proof of 5.2↩︎

We present below the proof of 3.

9 Experimental details↩︎

In this section, we provide some additional details on the experiments presented in the main text.

9.1 Random features experiments↩︎

In this short subsection, we quickly describe the empirical setup for the diffusion random feature experiment presented in 1.

9.1.0.1 Regression experiments.

In the regression experiment presented in 1, we use a classical random feature regression setting. We consider a data distribution of the form \[\require{upgreek} \begin{align} y_i = \langle \upbeta_\star, x_i \rangle + \sigma \epsilon_i \;, \quad x_i \sim \mathrm{N}(0, {\mathrm{I}}_d)\;, \quad \epsilon \sim \mathrm{N}(0, 1) \;, \end{align}\] where \(\epsilon_i\) and \(x_i\) are independent and we sample \(n\) i.i.d. pairs \((x_i,y_i)\) for \(i \in \left\{ 1,\dots,n \right\}\). We consider a model of the form \[\begin{align} f_{A,W}(x) = \frac{A}{\sqrt{p}} \varrho \left( \frac{W}{\sqrt{d}} x \right) \;, \end{align}\] where \(A \in \mathbb{R}^{d \times p}\), \(W \in \mathbb{R}^{p \times d}\) is a random matrix with i.i.d. standard Gaussian entries, and \(\varrho (x) := \max(x,0)\) is the ReLU activation.

We consider the minimization of the empirical risk \[\begin{align} \widehat{A}\in \mathop{\mathrm{arg\,min}}_A \left\{ \frac{1}{n} \sum_{i=1}^n (y_i - f_{A,W}(x_i))^2 + \frac{\lambda}{\sqrt{p}} \left\Vert A \right\Vert_{\mathrm{F}}^2 \right\} \;, \end{align}\] with \(\lambda\) a ridge regularization parameter.

Then, we evaluate the test (using a validation set) and train risk evaluated at the empirical risk minimizer \(\widehat{A}\) and report it in 1.

Hyperparameters details. In 1, we use \(n=10^2\), \(d=20\), \(N_g = 10\), \(\lambda = 10^{-4}\), and \(p \in [5, 4 \cdot 10^2]\).

9.1.0.2 Diffusion experiments.

We consider the following family of random features score networks, which is similar to existing works [19], [38], for \(x \in {\mathbb{R}^d}\) \[\begin{align} f_{A,W}(x) = \frac{A}{\sqrt{p}} \varrho \left( \frac{W}{\sqrt{d}} x \right) \;, \end{align}\] where \(A \in \mathbb{R}^{d \times p}\), \(W \in \mathbb{R}^{p \times d}\) is a random matrix with i.i.d. standard Gaussian entries, and \(\varrho (x) := \max(x,0)\) is the ReLU activation. For the data distribution, we take the data distribution to by \(\nu = \mathrm{N}(0, d^{-1}{\mathrm{I}}_d)\).

Given a fixed \(t > 0\), a random \(W\) and a dataset \(S^{(n)} := (Z_1,\dots,Z_n) \sim \nu^{\otimes n}\), we consider the minimization of the regularized empirical risk, defined by \[\begin{align} A \mapsto \frac{1}{n} \sum_{i=1}^n \mathbb{E}_{\Xi} \left[ \left\Vert \sigma_t f_{A,W} (\alpha_t Z_i + \sigma_t \Xi) + \Xi \right\Vert^2 \right] + \frac{\lambda}{\sqrt{p}} \left\Vert A \right\Vert_{\mathrm{F}}^2 \;, \end{align}\] with \(\Xi \sim \mathrm{N}(0,{\mathrm{I}}_d)\) and \(\lambda\) a ridge regularization parameter. An important difference with the regression case above is that, in order to make the empirical risk minimization tractable, we use a Monte Carlo approximation of the expectation over \(\Xi \sim \mathrm{N}(0,{\mathrm{I}}_d)\). Given an integer \(N_g \in \mathbb{N}^\star\) and \((\xi_{ij})_{1 \leq i \leq n,~1 \leq j \leq N_g} \sim \mathrm{N}(0,{\mathrm{I}}_d)^{\otimes (n N_g)}\), we therefore compute \[\begin{align} \widehat{A}\in \mathop{\mathrm{arg\,min}}_A \left\{ \frac{1}{n} \sum_{i=1}^n \sum_{j=1}^{N_g} {\left\Vert \sigma_t f_{A,W} (\alpha_t Z_i + \sigma_t \xi_{ij}) + \xi_{ij} \right\Vert^2} + \frac{\lambda}{\sqrt{p}} \left\Vert A \right\Vert_{\mathrm{F}}^2 \right\} \end{align}\] by noting that it can be written as the solution of a linear regression problem. This procedure is similar to [39].

After estimating \(A\), we compute the empirical and population denoising score matching loss, where the latter is estimated using a validation set. The hyperparameters details are provided below.

Hyperparameters details. In 1, we use \(\kappa = 1\), \(n=2 \cdot 10^4\), \(t = 10^{-1}\), \(d=3\), \(N_g = 10\), \(\lambda = 10^{-4}\), and \(p \in [2, 10^4]\).

9.2 High-dimensional experiments↩︎

Our experiments are performed using an implementation of the DDPM model from [2]. We use the configuration titledvp.ddpm.cifar10_continuous which implements the DDPM model of [1] but for the continuous-time variance preserving setting. The architecture is a U-Net [32] with the encoder and decoder each consisting of four resolution levels \((32\times32, 16\times16, 8\times8, 4\times4)\), with two residual blocks per level and utilizes self-attention at a resolution of \(16\times16\). The model is conditioned on time, with the timestep encoded via a sinusoidal embedding and injected into each residual block.

For figures 2 and [fig:train95unet] we modify the number of features (NF) parameter (denoted nf in the configuration file and taking a default value of \(128\)). This parameter defines the base channel width used throughout the model which is multiplied by \(m = 1, 2, 2, 2\) to produce the actual channel widths for each of the four resolution levels. The number of parameters in the model scales quadratically with channel width. We consider NF\(\in \{2, 4, 8, 16, 32, 64, 128\}\) to see how the number of parameters affects overparameterization. Some technical modifications to the code were required to support smaller values of NF.

For all high-dimensional experiments we use CIFAR-10 but restrict to the first \(1,000\) examples to make overfitting easier to achieve. Evaluation is performed on the full \(10,000\) held-out examples. We train the model with the default setup in the codebase (ADAM with \(10\%\) dropout and EMA decay \(0.9999\)) and we train for \(1,000,000\) iterations, which is larger than typical to better guarantee overfitting. In the experiments of Figure 2 we choose the checkpoint with smallest train loss. The samples generated on the right-hand side of Figure 2 uses the DDIM implementation in the codebase to generate samples with a fixed seed across NF values. We perform log-likelihood calculations using the bpd implementation in the codebase.

For the experiment in Figure [fig:no95time95exp], we consider the DDPM network but we restrict the access of the neural network to the time input. However, we do still scale the network properly according to the time input. We then train the model as usual but restrict the time variable to the range \([t, t + r]\) for \(t=0.1\) and \(r \in [0.0001, 0.00032, 0.001, 0.0032, 0.01, 0.032, 0.1, 0.5]\). We then evaluate the train and test error at \(t\) to obtain the plot.

9.3 Compute resources↩︎

Our experiments were run on an Amazon EC2 G6e (g6e.xlarge) instance which has a single NVIDIA L40S Tensor Core GPU.

9.4 Licenses↩︎

Codebases:

• Score-Based Generative Modeling through Stochastic Differential Equations [2]: Apache License 2.0.

Datasets:

• CIFAR-10 [59]: MIT license.

References↩︎

---

1. We refer to Theorem 2.4 of [54] for a discussion on the different conditions required by the Marchenko-Pastur theorem.↩︎