We thank the reviewer for their careful reading and insightful feedback. Below, we address each of the questions in turn:

Q1: Yes indeed, at line 94, $s_\theta$ denotes the score function and $D_\theta$ denotes the denoiser. We will add an explicit statement in the camera-ready manuscript to avoid any ambiguity.

Q2: Both forms in Eq.3 and Eq.6 are correct. Since $\Sigma$ is the data covariance, it’s real symmetric and positive (semi-)definite, thus it admits the eigen-decomposition, $\Sigma=U\Lambda U^\top$, i.e. the principal component of data. It follows that, $\Sigma+\sigma^2 I=U(\sigma^2+\Lambda)U^\top;\ (\Sigma+\sigma^2 I)^{-1}=U(\sigma^2+\Lambda)^{-1}U^\top$.

Thus all three matrices, $\Sigma$, $\Sigma+\sigma^2 I$ and $(\Sigma+\sigma^2 I)^{-1}$ share the same eigenbasis, hence commute. So Eq.6 and Eq.3 are equivalent to each other.

Q3: Thanks for noticing this differences!

At first glance one might hope to carry out the learning dynamics analysis in the Fourier basis, since convolution is diagonal there. However, when the filter support K is strictly smaller than the signal length $N$, the DFT coefficients of the truncated filter $\gamma_\sigma=\sqrt{N}Fw_\sigma$ lie in only a K-dimensional subspace of $\mathbb{C}^N$. In other words, the Fourier coefficients of $w_\sigma$ are no longer independent but constrained to the column space of the corresponding DFT submatrix, turning the gradient dynamics into a constrained optimization problem in Fourier space.

By contrast, in the spatial domain the $K$ raw filter parameters remain independent. This allows us to leverage the well-known properties of circulant and shift matrices and obtain a much cleaner analysis (see Appendix G.4.1 and Proposition 5.4). Intuitively, each local convolutional filter “sees” only the statistics of its $K$-pixel patch and uses those statistics to linearly predict the denoised value at the center pixel.

We will expand our discussion of this technical choice in the camera-ready version to make the motivation and technical benefits clearer.

We appreciate the reviewer’s suggestions and will incorporate these clarifications into the final manuscript.

We thank the reviewer for their thorough summary and thoughtful feedback!

W1: With respect, we find this critique not entirely accurate. First, while our theoretical analysis assumes a linear, independent denoiser, it imposes no constraint on the data distribution. Second, on the empirical side (Experiment 1, Fig. 9), we show that fully connected (MLP‐based) diffusion models quantitatively exhibit the predicted inverse‐variance spectral bias and scaling coefficient, despite violating both linearity and independence assumptions. This behavior qualitatively persists across multiple natural image datasets (Experiment 2, Fig. 11). These results indicate that, although our theory relies on simplified assumptions, these do not limit its applicability to real world diffusion models with fully connected structure. However, convolutional architectures do alter the observed learning dynamics. We believe a more targeted analysis of convolutional structures could explain these differences.

W2: We acknowledge that, due to page constraints, the treatment of deeper linear and convolutional cases in the main text is relatively brief. However, the full derivations and analyses are provided in Appendices F and G. In theory, the N-fold speed-up is rigorously established in Proposition 5.3. We recognize, however, that linking this result directly to the faster empirical convergence of practical U-Nets is more challenging—isolating a single architectural component responsible for this speed-up is nontrivial.

Additionally, we performed more extensive simulation of the distributional learning dynamics for patch convolutional network. We found generally the power law scaling is attenuated, and the specific scaling exponent is strongly affected by the patch size: smaller filter size leads to more simultaneous emergence of modes and closer to 0 scaling coefficients. We think these new numerical results can better link to our empirical observations of attenuated scaling in practical UNet training. The details of the new results are presented in the section at the end.

W3: Thank you for the suggestion. Although we originally presented Equation 7 inline given its straightforward derivation, we agree that it merits its own formal proposition alongside Proposition 5.1 in the camera-ready version.

Q1: We appreciate the reviewer’s scrutiny of Appendix Figures 17–20. In those figures, we quantify the distribution of generated patches and track how their covariances evolve along the principal components of training patches. While these plots hint at a patch-space spectral bias, the effect is relatively small and less consistent than in the MLP setting or as predicted by our theory. Given this variability, we refrain from making strong claims about spectral bias in practical U-Net learning at the patch level. We believe that a refined theory incorporating additional details of the U-Net architecture may better capture these empirical curves.

L1: We appreciate the reviewer drawing our attention to the concurrent work (arXiv:2505.17638). We value their careful experiments, analytical theory, and practical insights. But we regard our contributions as complementary. Their theory assumes data drawn from a white Gaussian distribution and leverages the Gaussian equivalence principle of random features to derive the feature covariance spectrum, naturally explaining the separation of timescales between generalization and memorization. However, it does not account for nontrivial data covariance structure or its interaction with sampling dynamics—precisely the focus of our paper.

In summary, there are multiple tractable models for analyzing learning dynamics: our work focuses on linear and deep linear networks, while the concurrent work explores random feature models. The linear setting renders the probability‐flow ODE sampling dynamics analytically tractable, enabling us to track the spectral evolution of the learned distribution. In contrast, random feature models facilitate discussions of feature number scaling. A promising future direction is to unify these approaches by analyzing random feature models with non-white data covariance, yielding a richer understanding of memorization and generalization as functions of the data spectrum. We will incorporate this comparison and discussion in the final version of our manuscript.

---

Numerical simulation of the patch convolution learning dynamics

We took the FFHQ32 dataset as our example, and used the statistics of patches cropped from it to compute the optimal linear filter at each noise scale $w_\sigma^*$. Then we used proposition 5.4 to obtain the learned filter at different training time, then solve the PF-ODE in the spatial domain to get samples. We computed their variance along eigenmodes and the scaling laws.

During numerical simulation, we found there are two key factors affecting the observed power law relation between image eigenvalues and convergence time, namely, the patch size of the filter $P$ and the initialization of the filter weights.

If filter weights are initialized as identity without scaling, the initial sample variance will be way too high, and the exponent will be around ~ $-0.4$ — already attenuated than the linear case solution. One practically relevant alternative is to initialize the filter weight as identity scaled by $\sigma_{data}^2 / (\sigma^2 + \sigma_{data}^2)$, inspired by the coefficient of skip connection $c_{skip}$ in EDM preconditioning scheme.

If so, we’d observe modes with both increasing and decreasing variances. The specific scaling relation is affected by the filter patch size (summarized in following table). We found that when the patch size become smaller, scaling exponent of decreasing mode approach zero. from around -0.4 for P=15 to around -0.21 for P=3 (3x3 convolution), and the convergence time generally increases with smaller patch sizes. For the mode with increasing variance, a smaller patch size also attenuated the exponent to $0.0$ (P=7) or even flip the law i.e. $0.48$ for $P=3$. Intuitively, we think, 3x3 convolution ties many modes together in the Fourier space, so during generation, many modes would be learned at the same time.

In summary, even we no longer have analytical prediction of the scaling coefficients, we can still obtain prediction through numerical simulation. It shows that the patch size of the convolution could significantly affect the scaling relation. We hypothesize that the attenuated or flat scaling relation for practical UNet trained on natural images could partially be related to its use of small 3x3 convolution.

Tab.S1 Spectral convergence time computed for patch linear convolutional networks | FFHQ32

We will add the corresponding figures to Appendix upon acceptance. If requested, we are more than happy to present the figures of these scaling laws or code generating them for reviewers to view them.

We appreciate the reviewer’s positive assessment and detailed comment on our paper! Here are the response to the questions.

Q1: The key dynamical features are the same between EDM and DSM formulations, i.e. gradient flowing along eigenmodes of data covariance or Fourier modes (convolutional case), since the two formalisms are equivalent up to scaling the state $x$ and time / noise scale (specifically if we denote $x_{t}(\epsilon\_{t})=\sqrt{\bar{\alpha}\_{t}}x\_{0}+\sqrt{1-\bar{\alpha}\_{t}}\epsilon\_{t}$ in DDPM formalism, then we can substitute our $x\_{t}\mapsto\frac{x\_{t}}{\sqrt{\bar{\alpha}\_{t}}}, \sigma\_{t}\mapsto\frac{\sqrt{1-\bar{\alpha}\_{t}}}{\sqrt{\bar{\alpha}\_{t}}}$ and translate our solution to DDPM formulation.).

The difference exists in the convergence speed and the asymptotic value of the weights at each noise scale, and (potentially) weighting of each noise scales. In the appendix, we provided a comprehensive table (Tab. 2, App. C.4) of the asymptotic value $w^*$ and convergence time scale $\tau^*$ for different variants of diffusion objective and their corresponding dynamics, including EDM, the classic DSM (eps objective), and recent flow matching.

We use the EDM formulation since it largely simplifies the notations in derivation: it uses a single symbol for noise scale, and the solutions we obtain can be translated to other formulations without much pain.

Q2: Random feature model is indeed an interesting and tractable next step, and we were planning to carry similar analysis for them.

As a matter of fact, some very recent work has carried out learning dynamics analysis for diffusion models focusing on the random feature model ([1] arxiv 2505.17638). The general structure of the learning dynamics remains very similar: leveraging the Gaussian equivalence principle, the key factor determining the dynamics is the covariance of the random feature of the first layer. Tools from statistical mechanics and random matrix theory can be leveraged to calculate the limiting spectrum of these covariances, which determines the convergence speed of eigenmodes. Allbeit interesting, the theoretical analysis in this work [1], was performed for the white gaussian data case. Generalization to non-white structured covariance case is likely tractable through random matrix tools, e.g. free probability.

In our camera ready version, we will add reference to this paper as concurrent work and discuss the connection and difference from our work, for readers who are interested.

One major challenge for the random feature set up is that the sampling dynamics become harder to solve using current technique (Lemma 4.1), since the score function are nonlinear in the state, the PF-ODE becomes a nonlinear dynamic system. So we can no longer track the evolution of generated distributions as easily as we did in the paper.

Q3: Thanks for mentioning the linking error, the Table on L130 should be Tab.2 and the section should be App. C.4. We will fix the typo in our final version.

Q4: In the main paper, proposition 5.4 provided the gradient flow solution of linear convolutional model with local filter, which is basically mode by mode convergence along patch eigenmodes. To obtain the sampled distribution, we’d need to translate these learned local filter to Fourier basis and solve the PF-ODE on the Fourier basis mode by mode. This part can be done by numerical integration. For any specific dataset, we can also instantiate the learned filter as an actual linear convolutional filter and solve the PF-ODE numerically in the spatial domain.

Numerical simulation of the patch convolution learning dynamics Specifically, we took the FFHQ32 dataset as our example, and used the statistics of patches cropped from it to compute the optimal linear filter at each noise scale $w_\sigma^*$. Then we used proposition 5.4 to obtain the learned filter at different training time, then solve the PF-ODE in the spatial domain to get samples. We computed their variance along eigenmodes and the scaling laws.

During numerical simulation, we found there are two key factors affecting the observed power law relation between image eigenvalues and convergence time, namely, the patch size of the filter $P$ and the initialization of the filter weights.

If filter weights are initialized as identity without scaling, the initial sample variance will be way too high, and the exponent will be around ~ $-0.4$ — already attenuated than the linear case solution. One practically relevant alternative is to initialize the filter weight as identity scaled by $\sigma_{data}^2 / (\sigma^2 + \sigma_{data}^2)$, inspired by the coefficient of skip connection $c_{skip}$ in EDM preconditioning scheme.

If so, we’d observe modes with both increasing and decreasing variances. The specific scaling relation is affected by the filter patch size (summarized in following table). We found that when the patch size become smaller, scaling exponent of decreasing mode approach zero. from around $-0.4$ for P=15 to around $-0.21$ for P=3 (3x3 convolution), and the convergence time generally increases with smaller patch sizes. For the mode with increasing variance, a smaller patch size also attenuated the exponent to $0.0$ (P=7) or even flip the law i.e. $0.48$ for $P=3$. Intuitively, we think, 3x3 convolution ties many modes together in the Fourier space, so during generation, many modes would be learned at the same time.

In summary, even we no longer have analytical prediction of the scaling coefficients, we can still obtain prediction through numerical simulation. It shows that the patch size of the convolution could significantly affect the scaling relation. We hypothesize that the attenuated or flat scaling relation for practical UNet trained on natural images could partially be related to its use of small 3x3 convolution.

Tab.S1 Spectral convergence time computed for patch linear convolutional networks | FFHQ32

We will add the corresponding figures to Appendix upon acceptance. If requested, we are more than happy to present the figures of these scaling laws or code generating them for reviewers to view them.

We appreciate the reviewer’s positive assessment of our work and the insightful comment regarding the attenuation of spectral bias in more realistic settings.

We have devoted additional effort and discussion to address this gap between theory and practice, through more extensive simulation of theory and ablation of Unet experiments.

On the theoretical side, we showed that the claimed spectral bias is a general effect intrinsic to the denoising score-matching objective, analogous to ridge-regression objectives. We demonstrated this through theoretical analysis of linear, deep linear, linear-convolution, and patch-convolution networks. After submission, we have conducted additional numerical simulations of the distributional learning dynamics for patch-convolution networks (Prop. 5.4). The new results are presented in the new section below. Briefly, in general the scaling relation for linear patch convolution network is weaker, with exponent around -0.40, this shows that the patch convolution already help learning high frequency patterns more efficiently. With proper initialization of the weights, we can obtain learning dynamics closer to reality; further, the exponent strongly depend on the patch size of the filter. Smaller filter leads to generally slower convergence of lower eigenmodes and smaller scaling exponent ~$-0.21$. So this new results shows that, the observed weaker scaling relation in UNet could be related to the small filter size (3x3) of the network, and modulated by the preconditioning and initialization.

On the empirical side, MLPs trained on Gaussian or image data exhibit power-law scaling with an exponent of about –0.4 on natural images and around –1 for Gaussian data. However, under the same setup, UNets learn all high-variance modes at roughly the same time and likewise for low-variance modes: variance-increasing modes still converge faster than variance-decreasing modes, so spectral bias persists, but the scaling exponent is effectively zero.

Our current hypothesis is that certain modern UNet design choices attenuate or circumvent the standard spectral bias, enabling more efficient convergence to the natural-image spectrum. To test this, we systematically ablated various UNet features on FFHQ32—specifically, multi-level/multi-resolution design, depth, and channel width. (Without multi-level design and skip connections, a UNet reduces to a ResNet.) Surprisingly, we found that channel width is the key factor affecting the spectral signature of learning. Reducing channels from 128 to 8 or 16 lengthened convergence times and revealed a more pronounced spectral bias—consistent with the theory for patch-linear convolution networks. When channel width is held constant, UNets with 1–4 levels exhibit the same near-zero exponent; similarly, ResNets of varying depth share this signature.

These findings highlight a limitation of existing theory: it typically assumes convolutional networks with few channels. Extending analyses to multi- or infinite-channel limits—perhaps via Convolutional Neural Tangent Kernel (CNTK) methods (e.g. arXiv:2203.09255)—could better capture the learning dynamics of practical high-dimensional UNets.

In summary, understanding the evolution of the sample distribution during training of practical diffusion models remains challenging, and we do not claim to solve this problem in one paper. Nonetheless, our contributions are twofold: (1) solving tractable cases of diffusion learning dynamics and demonstrating their relevance to deep nonlinear MLPs, which can be relevant in non-image, e.g. tabular settings; and (2) revealing qualitative differences in practical UNet and ResNet cases. We believe that future theoretical work must incorporate design elements such as local convolution and multi-channel structure—otherwise, generic setups (e.g. two-layer MLPs, random features, linear models) common in diffusion theory (e.g. arXiv:2502.00336; 2505.24769; 2505.17638) will likely fail to account for these spectral-learning phenomena. We hope this message resonates with the theory community.

---

Numerical simulation of the patch convolution learning dynamics We took the FFHQ32 dataset as our example, and used the statistics of patches cropped from it to compute the optimal linear patch filter at each noise scale $w_\sigma^*$. Then we used proposition 5.4 to obtain the learned filter at different training time, then solve the PF-ODE in the spatial domain to get samples. We computed their variance along eigenmodes and the scaling laws.

During numerical simulation, we found there are two key factors affecting the observed power law relation between image eigenvalues and convergence time, namely, the patch size of the filter $P$ and the initialization of the filter weights.

If filter weights are initialized as identity without scaling, the initial sample variance will be way too high, and the exponent will be around ~ $-0.4$ — already attenuated than the linear case solution. One practically relevant alternative is to initialize the filter weight as identity scaled by $\sigma_{data}^2 / (\sigma^2 + \sigma_{data}^2)$, inspired by the coefficient of skip connection $c_{skip}$ in EDM preconditioning scheme.

If so, we’d observe modes with both increasing and decreasing variances. The specific scaling relation is affected by the filter patch size (summarized in following table). We found that when the patch size become smaller, scaling exponent of decreasing mode approach zero. from around -0.4 for P=15 to around -0.21 for P=3 (3x3 convolution), and the convergence time generally increases with smaller patch sizes. For the mode with increasing variance, a smaller patch size also attenuated the exponent to $0.0$ (P=7) or even flip the law i.e. $0.48$ for $P=3$. Intuitively, we think, 3x3 convolution ties many modes together in the Fourier space, so during generation, many modes would be learned at the same time.

In summary, even we no longer have analytical prediction of the scaling coefficients, we can still obtain prediction through numerical simulation. It shows that the patch size of the convolution could significantly affect the scaling relation. We hypothesize that the attenuated or flat scaling relation for practical UNet trained on natural images could partially be related to its use of small 3x3 convolution.

Tab.S1 Spectral convergence time computed for patch linear convolutional networks | FFHQ32

We will add the corresponding figures to Appendix upon acceptance. If requested, we are more than happy to present the figures of these scaling laws or code generating them for reviewers to view them.