An Analytical Theory of Spectral Bias in the Learning Dynamics of Diffusion Models
Authors: Binxu Wang
Cengiz Pehlevan
Venue: NeurIPS 2025 spotlight ()
Keywords: Spectral bias
analytical theory
learning dynamics
diffusion
flow matching
deep linear network
Primary area: theory
TL;DR: Closed‑form analysis of diffusion learning dynamics uncovers an inverse‑variance law of distributional convergence; MLP‑UNets comply, while convolutional U‑Nets break it.

Abstract: We develop an analytical framework for understanding how the learned distribution evolves during diffusion model training. 
Leveraging the Gaussian equivalence principle, we derived exact solutions for the gradient-flow dynamics of weights in one or two layer linear or linear convolutional denoiser settings with arbitrary data, where linear networks converge along principal components, and convolutional networks converge along Fourier modes. 
Remarkably, these solutions allow us to derive the generated distribution in closed-form and its KL-divergence through training. 
These analytical results expose a pronounced \emph{spectral bias}, i.e. for both weights and generated distributions, the convergence time of a mode follows an inverse power law of its variance.  
Empirical experiments on both Gaussian and natural image datasets demonstrate that the power-law spectral bias—remain robust even when using deeper or convolutional architectures. 
Our results underscore the importance of the data covariance in dictating the order and rate at which diffusion models learn different modes of the data, providing potential explanations of why earlier stopping could lead to incorrect details in image generative model.

An Analytical Theory of Spectral Bias in the Learning Dynamics of Diffusion Models
Binxu Wang
Cengiz Pehlevan
Kempner Institute, Harvard University Boston, MA, USA
SEAS, Harvard University Cambridge, MA, USA cpehlevanêseas. harvard,
binxu_wangehms. harvard.edu
Abstract
We develop analytical framework for understanding how the generated distribu- tion evolves during diffusion model training. Leveraging a Gaussian-egu Gauusiann-equivaennee principle, solve the full-batch gradient-flow dynamics of linear and convolu- tional denoisers and integrate the resulting probability lytic expressions for the generated distribution. The theory exposes universal inverse-variance spectral law: the time for an eigen- or Fourier mode to match its target variance scales tered orders of magnitude sooner than low-variance (fine) detail. Extending the analysis deep linear networks and circulant full-width convolutions shows that weight sharing merely multiplies learning rates-accelerating but not eliminat- ing the bias-whereat local convolution introduces qualitatively different bias. Experiments on Gaussian and natural -image datasets confirm the spectral law persists deep MLP-basee UNet. Convolutional U-Nets, however, display rapid near-simultaneou emergence of many modes, implicating local convolution in re- shaping learning dynamics. These results underscore how data covariance governs the order and speed with which diffusion models learn, and they call for deeper investigation of the unique inductive biases introduced by local convolution.
ODE, yielding ana-
X-1, so high-variance (coarse) structure is mas-
Introduction
Diffusion models create rich data by gradually transforming Gaussian noise into signal, a paradigm that now drives state-of-the-art generation in vision, audio, and molecular design [1, 2, 3]. Yet two basic questions remain open. (i) Which parts of the data distribution do these models learn first, and which linger un- bearned-risking artefacts under early stop- ping? (ii) How does architectural inductive bias shape this learning trajectory? Ad- dressing both questions demands that track the evolution of the full generated dis- tribution during training and relate it to the network's parameterization.
Training
Linear
Leaminn dynamics (gradient
Woight DIT agaimade
VeCusM
Spectral bias power
Ejpnnede brins
Di(x.o)
Sampling dynamics (probabilitty flow) dxc (sa(x..
Canaesgenne time
Learned
distribution p(Xee)) MTTT
Figure 1: Spectral-bias schematic. Learning and sam- pling together impose variance-ordered bias along covariance eigenmodes.
We tackle the learning puzzle through the simplest tractable setting- denoisers-whenn where datasets become equivalent to Gaussian with matched mean and covariance. In this regime solve, in closed form, the nested dynamics of
39th Conference on Neural Information Processing Systems (NcurIPS 2025).

---

gradient-flow of the weights and the probability-foon fllow ODE that carries noise into data, leading to an analytical characterization of the evolution of the generated distribution. The analysis exposes an livere-varriane spectral law: the time required for an eigen-mode to match target variance scales like тą x A so high-variance directions corresponding to global structure are mastered orders of magnitude sooner than low-variance, fine-detail directions. Extending the analysis to deep linear and linear convolutional nets, we show how convolutional architecture redirect this bias to Fourier or patch space, and accelerate convergence via weight sharing.
Main contributions 1. Closed-form distribution dynamics. We derive exact weight and distribu- tional trajectories for one-layer, two-layer linear, and convolutional denoisers under full-batch DSM training. 2. Inverse-variance spectral bias. The theory reveals and quantifies a spectral-law ordering of mode convergence, offering one mechanistic explanation for early-stop errors. 3. Empirical validation in nonlinear neural nets. Experiments on Gaussian and natural-image datasets confirm the spectral-law in deep MLP-based diffusion. 4. Convolutional architectural shape learning dynamics. Experiments on convolutional UNet, showing rapid patch-first learning dynamics different from fullycconeected architectures.
2 Related Work and Motivation: Spectral Bias in Distribution Learning
Spectral structure of natural data Many natural signals have interesting spectral structures (e.g. image [4], sound [5], video [6]). For natural images, their covariance eigenvalues decay as a power law, and the corresponding eigenvectors can align with semantically meaningful patterns [4]. For faces, for instance, leading eigenmodes capture coarse, low-frequency shape variations, whereas tail modes encode fine-ge grained textures [7, 8]. Analyzing spectral effect on diffusion learning can therefore show which type of features the model acquires first and which remain slow to learn.
Hidden Gaussian Structure in Diffusion Model Recent work has shown, for most diffusion times, the learned neural score is closely approximated by the linear score of a Gaussian fit to the data, which is usually the best linear approximation [9, 10]. Crucially, this Gaussian linear score admits closed-form solution to the probability-flos flow ODE, which can be exploited to accelerate sampling and improve its quality [11]. Moreover, this same linear structure has been linked to the aanaiiiaaiioommmorizzttioo transition in diffusion models [10]. In sum, across many noise levels, the Gaussian linear approximation is predominant structure in the learned score. Thus, we hypothesize it will have a significant effect on the learning dynamics of score approximator. From this perspective, our contribution is to elucidate the learning process of this linear structure.
Learning theory for regression and deep linear networks Gradient dynamics in regression well-studied, with spectral bias and implicit regularisation emerging as central themes [12, 13, 14]. In Sec. 4.1, we show that the loss of linear diffusion model reduces to ridge regression, letting us import those results directly. Our analysis also builds on learning theory of deep linear networks (including convotiioovoluiinnall and denoising autoencoders) [15, 16, 17, 18]. We extend these insights to modern cestusioonbbssed generative models, offering closed-form description of how the generated distribution itself evolves during training.
Diffusion learning theory Several recent theory studies address diffusion models from a spectral perspective but tackle different questions. [19, 20, 21, 22] document spectral bias in the sampling process after training; our focus is on how that bias arises during training.
[23] study stochastic sampling assuming an optimal score, orthogonal to our analysis of training dynamics. Sharing our interest in training, [24] analyze learning of mixtures of spherical Gaussians to recover component means, whereas we tackle anisotropic covariances and track reconstruction of the full covariance. [25] characterises optimal score and distribution under constraints; results from our convolutional setup can be viewed through that lens.
3 Background
3.1 Score-based Diffusion Models
Let Px(x) be the data distribution of interest, and for each noise level > 0 define p(x;o) (po (0, (x) po(y) I dy. The associated score function is Vx logp(x;

---

i.c the gradient of the log density at noise scale a. In the EDM framework [26], one shows that the "probability flow" ODE
dx = Vx log p(::
(I)
da
exactly transports samples from pl (r) to p(- :;) as a decreases. In particular, integrating from oT down to a = recovers clean data samples from Po. We adopt the EDM parametrization for its notational simplicity; other common diffusion formalisms are equivalent up simple rescalings of space and time [26]. To learn the score of data distribution po (x), minimize the denoising score matching (DSM) objective [27] with function approximator, We reparametrize the score function with the 'denoiser' 8,(x, σ)) (D,(x, x)/o², then at noise level the DSM objective reads
Bxompo ((.1) (x oz;a) Xo
To balance the loss and importance of different noise scales, practical diffusion models all adopt certain weighting functions in their overall loss L So do w(o) Lo.
3.2 Gaussian Data and Optimal Denoiser
To motivate our linear score approximator set up, it is useful to consider the optimal score and the denoiser of Gaussian distribution. For Gaussian data Xy N(,, Xo E R and Σ is positive semi-definite matrix. When noising Xg by Gaussian noise at scale σ, the corrupted x satisfies (H, E + Ii),, for which the Bayes-optimal denoiser is an affine function of x.
(E o2)) -xx 14)
For Gaussian data, minimizing (2) yields D This solution has an intuitive interpretation, i.e. the difference of the state x and distribution mean was projected onto the eigenbasis and shrinked mode-by-mode by >x/(Ax, + o²). Thus, according to the variance Åş along target axis, modes with variance significantly higher than noise À >> o² will retained; modes with variance much smaller than noise will be "shrinked" out. Effectively o² defines threshold of signal and noise, and modes below which will be removed. This intuition similar to Ridge regression is made exact in Sec. 4.1.
4 Learning in Diffusion Models with a Linear Denoiser
Problem set-up. Throughout the paper, assume the denoiser each noise scale is linear (affine) and independent across scales:
Dx W. x b..
Since the parameters (Wa. b,} are decoupled across noise scales, each a can be analysed indepen- dently. Through further parametrization, this umbrella form captures linear residual nets, deep linear nets, and linear convolutional nets (see Sec. 5).
We train on an arbitrary distribution Po with mean and covariance Σ by gradient flow on the full-batch DSM loss, i.e. the exact expectation over data and noise (2). (In practice, one cannot sample all z values, but the full-batch limit yields clean closed-form dynamics.)
This setting lets dissect analytically the role of data spectrum, model architecture (W, parametri- sation), and loss variant in shaping diffusion learning.
4.1 Diffusion learning ridge regression
Gaussian equivalence. For any joint distribution p(X, Y') the quadratic loss
C(W,b) b) = Epcx,)) WX + b- 2
depends on only through the first two moments (X, Y); see App. C.I.1 for proof. Hence a linear denoiser trained arbitrary po interacts with the data solely via its mean μe and covariance Σ.
Instance for diffusion. Under EDM loss (2), the noisy input-target pair X Xo + oz, Y Xo- giving Σxx E Σ+od 01, Σvx- E.
3

---

Gradient and optimum. Differentiating and setting gradients to zero yields
Cw -22 +2 (2 I2I)
Cb n², Vb. Lo =2(b; 1abb-10 -I-
(5)
W = = min La (T
by -I- =I-
(6)
Thus the optimal linear denoiser reproduces the denoiser for the Gaussian approximation of data (3), and best achievable loss is set purely by the data spectrum.
Other objectives. While the main text focuses on the EDM loss (2), we have worked out the gradients, optima, and learning dynamics for several popular variants used in diffusion and flow-mathhing [28]
literature; these results are summarised in Tab. C.4 (derivations in App. C.4).
Ridge viewpoint. Because
C. Bx Wx + b - xx| [
full-batch diffusion at noise scale is simply auto-encoding with ridge regularisation of strength o² (App. C.2.1; ef. [29]). We will exploit classic ridge-regression results when analyzing learning dynamics in the following sections.
4.2
Weight Learning Dynamics of Linear Denoiser
With the gradient structure in hand, we solve the full-batch gradient- ODE,
dWo VW db,
Vb
(GF)
dr
dr
where r training time and ry the learning-rate.
Zero-mean data (μ = 0): Exponential convergence mode-by-mode Because the gradients to W, b decouple (5), the dynamics is simplified on the eigenbasis of the covariance. We diagonalize the covariance, Σ >d k= uu with orthonormal principal components (PC) uy and eigenvalues A (the mode variances). Projecting (GF) onto this basis yields the closed-form solution (derivation in App. D.1):
by b.) 205
W.(7) W (. [: [rJuuule (nr(a"+
Interpretation. Each eigenmode projec- tion of the weight W u, converges to the optimal value W un exponentially with rate (o² >x); hence (i) the weights larger noise a generally converge faster; (ii) fixed o, high-variance A modes converge first, while modes buried beneath the noise floor (Aş: o²) share the same slower timescale. Fig. 2A illustrates this cerctrum-ordere convergence, with high-variance modes reaching their optima before the low-variance ones (see also 5A).
Weight Scale along Egenmods
Generated Variance along Eigenmodss
Power Law of Mode Energence
rese 3.00
Sy-t syctett leyer demetser
Dee cenohe
Non-centred data (μx * 0): Interaction of mean and covariance learning. A non-zero mean introduces a rank-one cou- pling between W and b (matrix M in Prop D.I). Eigenmodes of weights over- lapping with the mean (u, μ o 0) now interact with b, producing transient over- shoots and other non-monotonic effects; orthogonal modes retain the exponential convergnece above. App. D.2 gives the full linear-system analysis and two-dimensional visualisations (Fig, 27).
Symmetric tro-lyyr dexbbr
Figure 2: Learning dynamics per eigenmode. Top: one-layer linear denoiser. Bottom: two-layer symmetric denoiser. (A.D) Weight trajectories u, W W.(T)u, (o=
1). (B,E) Generated-vvrianee λς versus target variance Åμ. (C,F) Power-law relation between emergence time
T and Ag-

---

4.3 Sampling Dynamics during Training
For diffusion models, our goal is the generated distribution, obtained by integrating the probabil- ity-foot flywfoo ODE (PF-ODE) backwards from a large to mmnn ~ O,
dx do
o "o - 1)x + bo],
(PF)
initialized with Gaussian noise X ~ M/O 131). For linear denoiser, the PF-ODE is an inhomoge- ncous affine system, its solution x, is necessarily affine function of the initial state xy [30], x( = A(ag; σT) x(or) + c(oo: or). Since the map is affine, the distribution of x(oo) remains Gaussian, with covariance A/oo( (T) (
However, in general, the state-transition matrix A(oo; σT) is hard to evaluate, as it involves time-ordered matrix exponential, and the weight matrices different noise scales W. may not com- mule. The analysis below-and our closed-form results-hinges on situations where commutativity maintained by gradient flow or architectural bias, thus removing the time-ordering operator.
Lemma 4.1 (PF-ODE solution for commuting weights). If the linear denoiser D(x; o) = W Wox+b, satisfies [W.. W W 0 for all o', then for any 0 < du oT,
x = A(ou, σr) x(or) + clou, oT), A(ou, σT) = expl
T W
ds
n0
Interpretation. For each common eigenvector ug. the term (u W. us - 1)/ar is the instantancous expansion (or contraction) rate of the sample variance along uù; the final variance is obtained by integrating this rate over noise scales a (see App. C.5).
When does commutativity hold? This arises in three common settings. (i) At convergence, this satisfied by the optimal weights W (6), which jointly diagonalize on eigenbasis of Σ. In such case, we recover the the closed-form solution to PF-ODE for Gaussian data, as found by [9, 31J. (ii) During training of linear denoisers, weights are initialized to be aligned with eigenhasis of Σ, then gradient flow keeps them aligned, preserving commutativity (iii) For linear convolutional denoisers, circulant weights share the Fourier basis and commute by construction (see Sec. 5.2). In these cases, the sampling process can be understood mode-by-mode Here we show the explicit solution for one layer linear denoiser.
Proposition 4.2 (Dynamics of generated distribution in one layer case). Assume (i) zero-mean data, (ii) aligned initialization W, (0( = Στ Qs u and (iii) gradient flow, full- full-batch training with learning rate 17. Then, while training the one-layer linear denoiser, the generated distribution at time T is ((u, Σ) with = Ex Au(T)uu and
T2....T)
Ae(T) o. (or,,))
Φu(a,T)= V exp 1-Qx Ei(-2(1ro*)e 2nr Bi 2n (nr(o*+Ae))
where is the exponential-integne function. (derivation in App. D.3)
Spectral bias. Figure 2B traces the variance trajectory As(T) for each eigen-mode. All modes begin with the same initiaiizatonninnuueed level, then follow sigmoidal curves to their targets, but in descending order of Å We define the first-passage time T as the training time at which (T) reaches the geometric (or harmonic) mean of its initial and target values. We find the first-passage time obeys an inverse law T x Ag a 1, (Fig. 2C), which implies that learning mode with variance 1/10 smaller takes roughly 10 times longer to converge. With larger weight initialization (larger 2x). the initial variance closer the target variance of some modes, then the inverse law splits into separate branches for modes with rising vs. decaying variance (Fig. 5B, Fig. 6).
Practical implication. This suggests when training stops earlier, the distribution in higher variance PC spaces have already converged, while low-variance ones-ofter the perceptual finer points such as letter strokes finger joints under-trained. This could be explanation for the familiar "wrong detail" artefacts in diffusion samples.
Deep and Convolutional Extensions
After analyzing the simplest linear denoiser, we set out to examine the effect of architectures via different parametrizations of the weights, specifically deeper linear models and linear convolutional networks. In the following, we will assume μ = O and focus on learning of covariance.
5

---

5.1 Deeper linear network
Consider a depth-L L linear denoiser D(x, σ) W.... Wix, where foo notational clarity- we suppress the explicit a o-dependence of weights. We assume aligned initialization, where for singular decomposition of each matrix, We(0) U2 Ae V the right basis of each layer matching the left basis of the next, V2+1 = Uz. VP 1, ..... 1..... 1, and with U, = V2 = U where U diagonalizes data covariance E. Then the total weight initialization is Woo (0) IT= Wc(O) = U[ Ar)UT With aligned initialization, every eigenmode learns independently seeppeneetly-mmmrrrinnn classical results [15, 32]. In our case, this also implies that the total weight II, Wi shares the eigenbasis U across training and noise scales, thus commute, making sampling tractable.
= One especially illuminating case is the two-layer symmetric network, where P.P Proposition 5.1 (Dynamics of weight and distribution in two layer linear model). Assume (i) centered data μ = 0; (ii) the weight matrix is initialized aligned, i.e. Po Po) ((()))))))))))))))))))))))))(( = Ex Q.u: then the gradient flow ODE admits closed-form solution (derivation in App. E.
W = Po(T) --
A +2
Qc
(8)
de Qx)e-8nAAn -8nAN + Qr.
The generated distribution at time T Ni, E. with E = Ex (o( and )u(T) = o. ((())
(o))
(1-Qx)- SurAp
Qu Ozle BB
()) = (o) (0) 2210 B
Axe -BrAA + Qk (1-r- (de + +2)) 22x 12
(1e
Interpretation. The learning dynamics of weights and variance along different principal components are visualized in Fig.2 D-F. Compared one-layer case, here, the weight converges along the PCs via sigmoidal dynamics, with the emergence time (reaching harmonic mean of initial and final value) Th = In 2/ (81) Ax). As for generated distribution, we find similar relationship between the target variance and emergence time Th x A a RV 1. For the more general non-aligned initialization, we show the non-aligned parts of weight will follow non-monotonic rise-and-fall dynamics (App. 1.1.2). Extensions to non-symmetric two layer model and deeper model were studied in App. F, which have similar bias but lack clean expressions.
5.2 Linear convolutional network
We consider linear denoiser with convolutional architecture, D( Wa x where samples X E IRN have Id spatial structure, and a width K convolution filter w, operates on it. The analysis could be easily generalized to 2d convolution. With circular boundary condition, Wa defines circulant weight matrix W, E NxN, where Wa * x = W, One favorable property of circulant matrices is that they are diagonalized by discrete Fourier transform F [33].
mk N
W FT
Fmn
exp
-2mi
(9)
Thus all weights Wa commutes, which allows us to leverage Lemma 4.1, and solve the sampling dynamics mode-by-mode on the Fourier basis, leading following result. Proposition 5.2. Linear convolutional denoisers with circular boundary can only model stationary Gaussian processes (GP), with independent Fourier modes, proof in App.G.2.
Learning dynamics of full-width filter K = N When convolution filter wo is as large as the signal, the gradient flow is diagonal and unconstrained in the Fourier domain. Thus, the analyses in Sec. 4 re-emerge with variance of Fourier mode Σạs taking the place of X Proposition 5.3 (Full-width circular convolution learning dynamics). Let D(x, W. with full-width filter K N, and truin Wa by full-batch grudient flow at rate 13. Then the weights at noise a and its spectral representation 7 evolves as
(T) 1 (1""((T))) 7u(T,0) = (o) + (Tu(T.0) - -K())) 2N
(10)
VN
where TE(0) = S and Eax = [F=EFFJa k is the variance of Fourier mode.
The generated distribution has diagonal covariance in the Fourier basis and follows exactly Prop. 4.2 after the replacement X + Enk. + N1j, U + F. (derivation in App. G.3)

---

Table I: Summary of theory. exp. and sigm. denotes exponential and sigmoidal convergence. xN denotes the N time speed up due weight sharing.
One layer
Sym. two layer
Full widhh linear conv
Patch linear copv
Param.
W
Pp"
W Cincl w)
W Circ(w), K < N
Weight dyn.
PC, exp. (?)
PC, sigm. 5.1
Fourier mode, exp, AN 3.3
PC ol palches, exp. XN 5.4
Learned distr.
Ganssian
Gaussian
Stationary 5.2.2
Stationary GP
Distr. dyn.
PC, sigm., power law 4.2
PC,
sigm., power law 5.1
Fourier, sigm., xN, power law 5.3
Fourier, N.S.
Interpretation.
The weight and distribution dynamics mirror the conneced case, with spectral bias towards higher variance Fourier modes; convolutional weight sharing simply multiplies every rate by N, accelerating convergence without altering the inverne-vrriace law.
Notably, the learned distribution is asymptotically equivalent the Gaussian approximation to the original training data with all possible spatial shifts as augmentations (proof in App. G.3.2). This one case where equivariant architectural constraints facilitates creativity as discussed [25]. Similarly, two-layer linear conv net with full-width filter can be treated as in Sec.5.1.
Learning dynamies of local filter K < N When the convolution filter has a limited bandwidth K # N, the Fourier domain dynamics get constrained, so it is casier to work with the filter weights. Let r be the half-width of the kernel (K = 2r + 1). Define the circular patch extractor P, (x(x) = [ i=1 RRKN,, and the patch covariance Eqatch N [x) (x E RRXK
Proposition 5.4 (Patch-convolution learning dynamics), For the circular convolutional denoiser, D(x, a) = Wa x trained by full-batch gradient flow with step size t. Let eq TRK be the one-hot vector with single at the center position r + (1-indexed). (derivation in App. G.4)
Wo (T( = w. expl 2NNrr(o Epmcer] (w.(0) -w:): w. = =0=0 Mpnch)
Epuchh..
Interpretation. Training with narrow convolutional filter reduces ridge regression in patch space. Under gradient flow, filter converges along eigenmodes of patch Epatchh: modes with larger variance converges sooner, those with smaller variance later, preserving the inverse-variance law. also enjoys the N times speed up given by weight sharing, accelerating progress without altering the ordering. The sampling ODE remains diagonal Fourier space, the generated distribution will be a stationary Gaussian process with local covariance structure shaped by the learned patch denoiser, though its exact form needs numerical integration to spell out. This setting similar the equivariant and local score machine described in [25], but with the additional linear constraint.
Simulation. We numerically simulated the dynamics of the sample distribution for linear patch- convolution denoisers using FFHQ dataset (details in App. 1.1.)). The spectral scaling exponents depend systematically on the convolutional patch size P: smaller kernels produced shallower, and in some cases even inverted, scaling relations (Tab. 2), potentially due to stronger coupling between more Fourier modes.
Empirical Validation of the Theory in Practical Diffusion Model Training
General Approach To test our theoretical predictions about the evolution of generated distribution (esp. covariance), we resort to the following method: 1) we training dataset (x,) and compute its empirical mean and covariance Σ. We then perform digen-decomposition of Σ, obtaining cigenvalues A and eigenvectors ug. 2) Next, we train a diffusion model on this dataset by optimizing
the DSM objective with a neural network denoiser D, (x, o). 3) During training, at certain steps we generate samples {x{} } from the diffusion model by integrating the PF-ODE (I). We then estimate the sample mean ja and sample covariance Σr_ Finally, we compute the variance of the generated samples along the eigenbasis of training data, AΙ u. To stress test our theory and maximize its relevance, we'd keep most the training hyperparameters as practical ones.
6.1 Multi-Layer Perceptron (MLP)
To test our theory about linear and deep linear network (Prop.4.2,5.1), we used Multi-Layer Perceptron (MLP) inspired by the SongUnet in EDM [26, 34] (details in App. I.2). We found
7

---

Variance along data eigenbanis
projected
Target emergence
308
107
Figure 3: Spectral Learning Dynamics of MLP-UNet (FFHQ32). A. Generated samples during training. B. Evolution of sample variance (x( across eigenmodes during training. C. Heatmap of variance trajectories along all eigenmodes, with dots marking mode emergence times T" (first-passage time at the geometric mean initial and final variances). The gray zone (0.5-2x target variance) indicates modes starting too close to their target, causing unreliable T' estimates. D. Power-law scaling of 7* versus target variance λ₁- separate law was fit for modes with increasing and decreasing variance, excluding the middle gray-zone eigenmodes for stability.
this architecture effective in learning distribution like point cloud data (Fig. 29). We kept the preconditioning, loss weighting and initialization the same as in [26].
Experiment 1: Zero-mean Gaussian Data MLP We first consider zero Gaussian NO, Σ) in dimension as training distribution, with covariance defined as a randomly rotated diagonal matrix with log normal spectrum (details in App. 1.4). During training, the generated variance of each eigenmode follows a sigmoidal trajectory toward its target value > modes with larger Å cross the plateau sooner (Fig. 10A). We mark the emergence time 7" the step which the variance reaches the geometric mean of its initial and asymptotic values (Fig. 10B). Across both high- and low-variance modes, r* obeys inverse power-law, T" x > With higher-dimensional Gaussians the exponent is estimated more precisely and remains close to I: for d = 256, a = 1.08; for d 512, Anner 1.05 and Odeer 1.13 (Fig. 10C). The scaling breaks down only for modes whose initial variance is already near As in that regime the trajectory is less sigmoidal and r* becomes ill-defined. This result shows that despite many non-idealistic conditions e.g. deeper network, nonlinear activation function, residual connections, normal weights initialization, shared parametrization of denoisers at different noise level, the prediction from the linear network theory still quantitatively correct.
Experiment 2: Natural Image Datasets MLP Next, validated our theory on natural image datasets. We flattened the images as a vectors, and trained deeper and wider MLP- UNet learn the distribution. Using FFHQ as our rumning example, moni- toring the generated samples through- out training (Fig. 3A), despite heavy noise early on, the coarse facial con- tours correspddng to the mean and top principal components of human face distribution [7]-emerge quickly, whereas high-frequency details (lower PCs) only appear later. We note that this spectral ordering effect of training dynamics is reminiscent and similar to that in the sampling dynamics after training [19, 22].
A.
B.
Evolution of PC projected variance
220
120
Figure 4: Learning dynamics of UNet differs FFHQ32. A. Sample trajectory from CNN-UNet. B. Variance evolution along covariance eigenmodes. (c.f. Fig. 3A.C.)
Quantitatively, the sample covariance Σσ rapidly aligns with and becomes close to diagonal in the data eigenbasis U (Fig. 11). The top eigenmodes' variances, Å (r). follow sigmoidal trajectories converging to their targets, and their "emergence times" 7F increase down the spectrum (Fig. 3B,C). We exclude central band of modes whose initial variances lie within 0.5-2x the target, since their undulating learning dynamics make first-passage time estimates unreliable. After this exclusion. modes with increasing and decreasing variance each exhibit clear power-law scaling between

---

emergence step r" and target variance Ag. with exponents -0.48 (R² = 0.97, N = 57) and - 00.35 (R² = 0.92, N 2, 914), respectively (Fig. 3D). Although the observed spectral bias slightly attenuated relative to the Gaussian case and linear theory prediction, it remains robust and consistent across datasets (MNIST, CIFAR-1O, FFHQ32 and AFHQ32) (App. B.2.2). This shows that even with natural image data, the distributional learning dynamics of ML.P diffusion still suffers from slower convergence speed for lower eigenmodes.
6.2
Convolutional Neural Networks (CNNs)
Next we turn to the convolutional U-Net-the work-horse of image-diffusion models [34, 35]. For a full-width linear convolutional network our analysis predicts an liversevvariane law in Fourier space (Prop. 5.3). The autch-convolution variant lacks a clear forecast on distribution, so the following experiments probe empirically whether-ant how-its learning dynamics is affected by spectral bias.
Experiment 3: Natural Image Datasets x CNN UNet Training on the same FFHQ dataset, the distributional learning trajectory of CNN-UNet is markedly different from the MLPs: early in training, we do not see contour of face, but locally coherent patches, reminiscent of Ising models (Fig. 4A.). Visually and variance-wise, the CNN-based UNet converge much faster and better than the MI.P-based UNet, matching the N-fold speed-up from weight sharing (Prop. 5.4; Fig. 4B). When projecting onto the data eigenbasis, all eigenmodes with increasing variance rise simultaneously, while cigenmodes with decreasing variance co-decay at later time, giving an effective power-law exponent a 0; Thus, spectral bias essentially absent (Fig. 13C.D.).
Why is spectral bias absent? On the theory side, the likely cause is locality: local convolutional filters couple neighbouring pixels, binding many Fourier modes into one learning unit. Because sampling remains diagonal in Fourier space, a broad band of modes is amplified simultaneously. attenuating the spectral ordering, as observed numerically in App. 1.1.2. In line with this, early in training, the CNN denoiser is indeed well-approximated by a local linear filter (Fig. 15).
On the empirical side, the key factor appears to be network width. We systematically varied the channel number and depth of deep convolutional denoisers (App. B.2.4), and found that-regardless of depth-narrowe networks (e.g. ch 4) exhibit slower convergence and stronger spectral ordering, consistent with the patch-convolution theory (Fig. 22,23). In contrast, wide networks with many channels (ch = 128) learn spectral modes almost instantaneously, similar to observations in practical UNet. In hindsight, the analytic theory effectively assumes a convolution with the same number of channels as the input (e.g., RGB 3), so the ratio between the network's 's channel and the input channel likely governs the deviation from theoretical predictions.
A complete analytic treatment of convolutional U-Net training dynamics is left for future work.
Discussion
In summary, we presented closed-font solutions for training denoisers with linear, deep linear or linear convolutional architecture, under the DSM objective on arbitrary data. This setup allows for a precise mode-wise understanding of the gradient flow dynamics of the denoiser and the evolution of the learned distribution: covariance eigenmode for deep linear network and Fourier mode for convolutional networks. For both the weights and the distribution, we showed analytical evidence of spectral bias, i.e. weights converge faster along the eigenmodes or Fourier modes with high variance, and the learned distribution recovers the true variance first along the top eigenmodes. These theoretical results are summarized in Tab. 1.
We hope these results can serve as solvable model for spectral bias in the diffusion models through the nested training and sampling dynamics. Furthermore, our analysis is not limited to the diffusion and the DSM loss, in App. H, we showed similar derivation for the spectral bias in flow matching models [28, 36].
Relevance of our theoretical assumptions We found, for the purpose of analytical tractability, we made many idealistic assumptions about neural network training, 1) linear neural network, 2) small or orthogonal weight initialization, 3) "full-batch" gradient flow, 4) independent evolution of weights at each noise scale. In our MI.P experiments, we found even when all of these assumptions were somewhat violated, the general theoretical prediction is correct, with modified power cocfficients.

---

This shows most of these assumptions could be relaxed in real life, and the spectrum of data indeed have a large effect on the learning dynamics, esp. for fully connected networks.
Inductive bias of the local convolution In our CNN experiments, however, the theoretical pre- dictions from linear models deviate: the spectral bias learning speed does not directly apply to the distribution of full images. Although our theory predicts that filter-weight learning dynamics are governed by the patch covariance, the ultimate image distribution shaped by the convolution of those filters. To date, many learning-thens
analyses for diffusion models assume MLP-like architectures [24]. For future theoretical work on the learning dynamies of practical diffusion models, rigorous treatment of the local convolutional structure-ait its requency-couuliigg effects-wil likely be essential, rather than relying on full- wilth convolution analyses [37].
Implications for high channel inputs Our ablation suggests that the ratio between network width and input channel count may underlie the observed deviations from theoretical predic- tions- Seess-sicciiiilll....... the absent of spectral bias and the near- cear-simultaneous convergence of eigenmodes. While most image and latent representations traditionally have few channels (e.g., 3 for RGB, 4 for latent diffusion [38]), recent architectures employ much higher channel counts- susss--schhhh as DC-AEs with 64-128 channels [39] or mucoder-based diffusion models with ch = 768 [40]. In these regimes, where input channel becomes comparable to that of the UNet DiT, the theory's predictions may become increasingly relevant for understanding, regularizing, and stabilizing training dynamics in gogh-channelldimensiinaal diffusion models.
Broader Impact Although our work is primarily theoretical, the inverse scaling law could offer valuable insights into how to improve the training of large-scale diffusion flow generative models.
Acknowledgements. We are grateful to Jacob Zavatone- Veth, Yongyi Yang, Michael Albergo, Hugo Cui, Yue Lu, and Ekdeep Lubana for their insightful discussions and valuable pointers to relevant literature. We appreciate Haim Sompolinsky and Zhengdao Chen for their meticulous proofreading and commenting on earlier versions of the manuscript. We also thank Thomas Fel for providing code formatting snippet. This work has been made possible in part by a gift from the Chan Zuckerberg Initiative Foundation establish the Kempner Institute for the Study of Natural and Artificial Intelligence, and by the Institute's 's generous support and computing resources. B W. supported by Kempner fellowship. C.P. is supported by NSF CAREER Award (IIS-2339700), DARPA grants DIAL-FP-038 and AIQ-HRO01 12520041, the Simons Collaboration on the Physics of Learning and Neural Computation, and the William Milton Fund from Harvard University.
10

---

References
[I] Jascha Sohl-Dickstein, Eric A. Weiss, Niru Maheswaranathan, and Surya Ganguli. Deep unsupervised learning using nonequilibrium thermodynamics. In Proceedings of the 32nd International Conference on Machine Learning (ICML), 2015.
[2] Jonathan Ho, Ajay Jain, and Pieter Abbeel. Denoising diffusion probabilistic models. In Advances in Neural Information Processing Systems (NeurIPS), 2020.
[3] Prafulla Dhariwal and Alexander Nichol. Diffusion models beat gans on image synthesis. In Advances in Neural Information Processing Systems (NeurIPS), 2021.
[4] Antonio Torralba and Aude Oliva. Statistics of natural image categories. in neural systems, 14(3):391, 2003.
Network: computation
[5] Hagai Attias and Christoph Schreiner. Temporal low-order statistics of natural sounds. Advances in neural information processing systems, 9, 1996.
[6] Dawei W Dong and Joseph Atick. Statistics of natural time-varying images. computation in neural systems, 6(3):345, 1995.
Network:
[7] Matthew Turk and Alex Pentland. Eigenfaces for recognition. Journal of cognitive neuroscience, 3(1):71-86, 1991.
[8] Binxu Wang and Carlos R Ponce. A geometric analysis of deep generative image models and its applications. In International Conference on Learning Representations, 2021. URL https://openrewe
mett/forummiidaa/@RRUUDDD
[9] Binxu Wang and John J. Vastola. The Hidden Linear Structure in Score-Based Models and its Application. arXiv e-prints, art. (xXiv:2311.1089)) November 2023. doi: 10.485500arxw.2311.. 10892.
[10] Xiang Li, Yixiang Dai, and Qing Qu Understanding generalizability of diffusion models requires rethinking the hidden gaussian structure. arXiv preprint arXiv:2410.24060, 2024.
[11] Binxu Wang and John Vastola. The unreasonable effectiveness of gaussian score approxima- tion for diffusion models and its applications. Transactions on Machine Learning Research, December 2024. arXiv preprint arXiv:2412.09726.
[12] Nasim Rahaman, Aristide Baratin, Devansh Arpit, Felix Draxler, Min Lin, Fred Hamprecht, Yoshua Bengio, and Aaron Courville. On the spectral bias of neural networks. International Conference Machine Learning, pages 5301-5310. PMLR, 2019.
[13] Blake Bordelon, Abdulkadir Canatar, and Cengiz Pehlevan. Spectrum dependent learning curves in kernel regression and wide neural networks. In Hal Daumé III and Aarti Singh, editors, Proceedings of the 37th International Conference Machine Learning, volume 119 of Proceedings of Machine Learning Research, pages 1024-1034. PMLR, 13-18 Jul 2020. URL https://proceediings mlr. wwwps/sssvvvvvem 19tbordeelon20aaa html.
[14] Abdulkadir Canatar, Blake Bordelon, and Cengiz Pehlevan. Spectral bias and task-model exigann generalization in kernel regression and infinitely wide neural networks. Nature Communicattonns, 12(1):2914, May 2021. ISSN 2041-1723 doi: 10. 10.1008//41467--021-31011311 URL https::///oll org/10. 1038/ 4444474221-221110
[15] Andrew M Saxe, James L McClelland, and Surya Ganguli. Exact solutions to the nonlincar dynamics of learning in deep linear neural networks. arXiv preprint arXiv:1312.6120 2013.
[16] Suriya Gunasekar, Jason D Lee, Daniel Soudry, and Nati Srebro. Implicit bias of gradient descent on linear convolutional networks. Advances in neural information processing systems,
3 31, 2018.
[17] Arnu Pretorius, Steve Kroon, and Herman Kamper. Learning dynamics of linear denoising autoencoders. In International Conference on Machine Learning, pages 4141-4150. PMIR, 2018.
11

---

[18] Chulhee Yun, Shankar Krishnan, and Hossein Mobahi. A unifying view on implicit bias in training linear neural networks. arXiv preprint arX1v:2010.02501, 2020.
[19] Sander Dieleman. Diffusion spectral autoregression, 2024. URL htt.t:tt////nteeeee http://andrr..a///2//. http://specrralauutrreeeeeeeeeee
html.
[20] Severi Rissanen, Markus Heinonen, and Aroo Solin. Generative modelling with inverse heat dissipation. arXiv preprint arXiv:2206. 2022.
[21] Florentin Guth, Simon Coste, Valentin De Bortoli, and Stephane Mallat. Wavelet score-based generative modeling. Advances in neural information processing systems, 35:478 491, 2022.
[22] Binxu Wang and John J Vastola. Diffusion models generate images like painters: an analytical theory of outline first, details later. arXiv preprint ww320:23033022490 2023.
[23] Giulio Biroli, Tony Bonnaire, Valentin De Bortoli, and Mare Mézard. Dynamical regimes of diffusion models. Nature Communications, 15(1):9957, 2024.
[24] Kulin Shah, Sitan Chen, and Adam Klivans. Learning mixtures of gaussians using the ddpm objective. Advances in Neural Information Processing Systems, 191196969999999 2023.
[25] Mason Kamb and Surya Ganguli. An analytic theory of creativity in convolutional diffusion models. arXiv preprint arXiv:2412.20222, 2024.
[26] Tero Karras, Miia Aittala, Timo Aila, and Samuli Laine. Elucidating the design space of diffusion-based generative models. arXiv preprint arXiv:2206.0086 003404 2022.
[27] Pascal Vincent. A connection between score matching and denoising autoencoders. Neural Computation, 23(7):1661-177 2011.
[28) Yaron Lipman, Marton Havasi, Peter Holderrieth, Neta Shaul, Matt Le, Brian Karrer, Ricky TQ Chen, David Lopez-Paz, Heli Ben-Hamu, and Itai Gat. Flow matching guide and code. arXiv preprint arXiv:2412.06624, 2024.
[29] Chris Bishop. Training with noise equivalent tikhonov regularization, Neural computa- tion, 7(1):108-116, 1995.
[30] Philip Hartman. Ordinary differential equations. SIAM, 2002.
[31] Emile Pierret and Bruno Galerne. Diffusion models for gaussian distributions: Exact solutions and wasserstein enror.. arXiv preprint arXiv:2405.142500 2024.
[32] Kenji Fukumizu. Effect of batch learning in multilayer neural networks. Gen, 1004):1E-03, 1998.
[33] Robert M Gray et al. Toeplitz and circulant matrices: A review. Foundations and Trends in Communications and Information Theory, 2(3):155-23) 2006.
[34] Yang Song, Jascha Sohl-Dickstein, Diederik P Kingma, Abhishek Kumar, Stefano Ermon, and Ben Poole. Score-based generative modeling through stochastic differential equations. In International Conference on Learning Representations, 2021. URL. https://openrevieer
net/forum?id-PXTTGG22RR 2RRHS..
[35] Olaf Ronneherger, Philipp Fischer, and Thomas Brox. U-net: Convolutional networks for biomedical image segmentation. In Medical image computing and computer-assisted intervention-MICCA. 2015: 18th international conference, Munich, Germany, October 5-9, 2015, proceedings, part III 18, pages 234-241. Springer, 2015.
[36] Yaron Lipman, Ricky TQ Chen, Heli Ben-Hamu, Maximilian Nickel, and Matt Le. Flow matching for generative modeling. arXiv preprint 10X0v:2210.02747, 2022.
[37] Suriya Gunasekar, Jason D Lee, Daniel Soudry, and Nati Srebro. Implicit bias of gradient descent on linear convolutional networks. Advances neural information processing systems,
31, 2018.
12

---

[38] Robin Rombach, Andreas Blattmann, Dominik Lorenz, Patrick Esser, and Bjöm Ommer. High- resolution image synthesis with latent diffusion models. In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 10684-10695, 2022.
[39] Junyu Chen, Han Cai, Junsong Chen, Enze Xie, Shang Yang, Haotian Tang, Muyang Li, Yao Lu, and Song Han. Deep compression autoencoder for efficient high-resolutie diffusion models. arXiv preprini arXiv:2440. 10733, 2024.
[40] Boyang Zheng, Nanye Ma, Shengbang Tong, and Saining Xie. Diffusion transformers with representation autoencoders. arXiv preprint (tt::2500./1990 2025.
[41] Jonathan Ho, Ajay Jain, and Pieter Abbeel. Denoising diffusion probabilistic models. Advances in Neural Information Processing Systems, 33:6840-6855 2020.
[42] Marvin Li, Aayush Karan, and Sitan Chen. Blink of an eye: a simple theory for feature localization generative models. arXiv preprint urXiv:2502.00921, 2025.
[43] Minyoung Huh, Hossein Mobahi, Richard Zhang, Brian Cheung. Pulkit Agrawal, and Phillip Isola. The low-rank simplicity bias in deep networks. arXiv preprint arXiv:2 10427, 2021.
[44] Li Jing, Jure Zhontar, et al. Implicit rank-minimizing autoencoder. Advances in Neural Information Processing Systems, 33: 147461 2020.
[45] Shuxuan Guo, Jose M Alvarez, and Mathieu Salzmann. Expandnets: Linear over-
parameterization to train compact convolutional networks. Advances in Neural Information Processing Systems, 13:1228-1310, 2020.
[46] Meena Jagadeesan, Ilya Razenshteyn, and Suriya Gunasekar. Inductive bias of multi-channel linear convolutional networks with bounded weight norm. In Conference on Learning Theory, pages 2276-2325. PMLR, 2022.
[47] Kathlén Kohn, Thomas Merkh, Guido Montúfar, and Matthew Trager. Geometry of linear convolutional networks. SIAM Journal on Applied Algebra and Geometry, 6(3):368 406, 2022.
[48] Kathlén Kohn, Guido Montúfar, Vahid Shahverdi, and Matthew Trager. Function space and critical points linear convolutional networks. SIAM Journal on Applied Algebra and Geometry, 8(2):333-362, 2024.
[49] Dmitry Ulyanov, Andrea Vedaldi, and Victor Lempitsky. Deep image prior. In Proceedings of the IEEE conference on computer vision and pattern recognition, pages 9446-9454, 2018.
[50] Prithvijit Chakrabarty and Subhransu Maji. The spectral bias of the deep image prior. arXiv 2019.
preprint
[51] Zezhou Cheng, Matheus Gadelha, Subhransu Maji, and Daniel Sheldon. A bayesian perspective on the deep image prior. In Proceedings of the JEEE/CVE Conference on Computer Vision and Pattern Recognition, pages 5443-5451, 2019.
[52] Basri Ronen, David Jacobs, Yoni Kasten, and Shira Kritchman. The convergence rate of neural networks for learned functions of different frequencies. Advances in Neural Information Processing Systems, 32, 2019.
[53] Ronen Basri, Meirav Galun, Amnon Geifman, David Jacobs, Yoni Kasten, and Shira Kritchman. Frequency bias neural networks for input non-uniform density. In International conference on machine learning, pages 685-694. PMLR, 2020.
[54] Mingchen Li, Mahdi Soltanolkotabi, and Samet Oymak. Gradient descent with early stopping provably robust to label noise for overparameterizee neural networks. In International conference artificial intelligence and statistics, pages 4313-4324. PMLR, 2020.
[55] Reinhard Heckel and Mahdi Soltanolkotabi. Denoising and regularization via exploiting the structural bias of convolutional generators. arXiv preprint arXiv: 14634, 2019.
[56] James R Bunch, Christopher P Nielsen, and Danny C Sorensen. Rank-one modification of the symmetric eigenproblem. Numerische Mathematik, 31(1):31-48, 1978.
13

---

[57] Ming Gu and Stanley C Eisenstat. A stable and efficient algorithm for the rank-one modification of the symmetric eigenproblem. SIAM journal on Matrix Analysis and Applications, 15(4): 1266-1276, 1994.
[58] Gene H Golub. Some modified matrix eigenvalue problems. SIAM review, 15(2)3118-334, 1973.
[59] Gilbert Strang. A proposal for Toeplitz matrix calculations. Studies in Applied Mathematics, 74 (2):171-176, 1986.
[60] Mingkui Chen. On the solution of circulant linear system. Technical Report TR-401, Yale University, Department of Computer Science, 1985. URL https://cpscc yale default/files/files//ttee pdf.
yale. www.eeeeeeeeeeeeeeeeeeeeeeeeeeee edu/sites/
[61] Michela Mastronardi. Presentation slides amastronardi.pdf) https://www math. unipd. it/ pdf, 2007. Accessed: 2025-05-23.
[62] Yang Song and Stefano Ermon. Generative modeling by estimating gradients of the data distribution. In H. Wallach, H. Larochelle, A. Beygelzimer, d'Alché-Buc E. Fox, and R. Garnett, editors, Advances in Neural Information Processing Systems, volume 32. Cur- ran Associates, Inc., 2019. URL https: httpss////.occenngs.... neur ips. ctccppppe//211/ffl// http2//44077553770/9666677777770
.pdf.
14

---

NeurIPS Paper Checklist
1. Claims
Question: Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope?
Answer: [Yes]
The abstract and introduction clearly state three contributions -losiiddss--sssd-eeeeeeeeee form linear diffusion theory, the inverse-variants ("spectrall) bias, and empirical verification in MLPs and CNNs-which are exactly what the body of the paper and appendices deliver.
Justification:
2. Limitations
Question:
Does the paper discuss the limitations of the work performed by the authors?
Answer: [Yes]
Justification: In discussion (paragraph "Relevance of our theoretical assumptions") we ex- plicitly listed idealised assumptions-
full-batch gradient flow, aligned initialization and discussed how much they affects upplicability to practical training of diffusion. The main theory results still hold even when all these assumptions were violated for MLPs.
3. Theory assumptions and proofs
Question: For each theoretical result, does the paper provide the full set of assumptions and a complete (and correct) proof?
Answer: [Yes]
Assumptions are stated at the start of every proposition/lemma (e.g. Lemma 4.1, Proposition 4.2); full proof's given in the referenced appendices, which submitted separately at the supplementary material deadline.
Justification:
4. Experimental result reproducibility
Question: Does the paper fully disclose all the information needed to reproduce the main ex- perimental results of the paper to the extent that affects the main claims and/or conclusions of the paper (regardless of whether the code and data are provided or not)?
Answer: [Yes]
Justification: The appendix Sec.1 has sufficient information of the model architecture, training details and hyperparameters to reproduce experimental results.
5. Open access to data and code
Question: Does the paper provide open access to the data and code, with sufficient instruc- tions to faithfully reproduce the main experimental results, as described in supplemental material?
Answer: [Yes]
Justification: All the training datasets used (MNIST, CIFAR, FFHQ, AFHQ) are open datasets available for downloading. The key MLP-UNet architectures and its torch imple- mentations code provided in App. I.2, full experiment code and analysis code will be shared after the paper gets accepted.
6. Experimental setting/details
Question: Does the paper specify all the training and test details (c.g., data splits, hyper- parameters, how they were chosen, type of optimizer, etc.) necessary to understand the results?
Answer: [Yes]
Justification: Appendix Iec.I supplied sufficient details for model training and evaluation.
7. Experiment statistical significance
Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?
Question:
Answer: [Yes]
15

---

Justification: For quantification of spectral law in experiments, we quantify the statistical power of law with linear regression line fitting and reported the R2 as goodness of fit Fig. 13.
8. Experiments compute resources
Question: For each experiment, does the paper provide sufficient information on the com- puter resources (type of compute workers, memory, time of execution) needed reproduce the experiments?
Answer: [Yes]
Justification: Details about required computational resources are provided in appendix I.1.
9. Code of ethics
Question: Does the research conducted in the paper conform, in every respect, with the NeurIPS Code of Ethics https://neurre: http:/eeerrrsccccceeeeeeeoccceee
Answer: [Yes]
Justification:
The work uses only publicly available numeric/image data, involves no per- sonal or sensitive information, and follows standard ML ethical guidelines.
10. Broader impacts
Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?
Answer: [Yes]
Justification: The final paragraph of the paper ("Broader Impact") outlines how understand- ing spectral bias could improve training efficiency while acknowledging possible misuse for more realistic deep-fake generation.
11. Safeguards
Does the paper deseribe safeguards that have been put in place for responsible release of data or models that have a high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?
Question:
Answer: [NA]
Justification: No high-risk models new datasets are released; experiments small-scale public data only.
12. Licenses for existing assets
Question: Are the creators or original owners of assets (e.g.. code, data, models), used in the paper, properly credited and are the license and terms of use explicitly mentioned and properly respected?
Answer: [Yes]
Justification: MNIST cited with its BSD-style licence; EDM implementation and the related datasets (CIFAR1O, FFHQ, AFHQ) acknowledged and used under its Apache-2.0 licence (Appendix A.4).
13. New assets
Question: Are new assets introduced in the paper well documented and is the documentation provided alongside the assets?
Answer: [NA]
Justification: The paper introduces no new datasets pretrained models.
14. Crowdsourcing and research with human subjects
Question: For crowdsourcing experiments and research with human subjects, does the paper include the full text of instructions given to participants and screenshots, if applicable, as well details about compensation (if any)?
Answer: [NA]
Justification: No human-subject data crowdsourcing was involved.
16

---

15. Institutional review board (IRB) approvals or equivalent for research with human subjects
Question: Does the paper describe potential risks incurred by study participants, whether such risks were disclosed to the subjects, and whether Institutional Review Board (IRB) approvals (or equivalent approval/review based the requirements of your country institution) were obtained?
Answer: [NA]
Justification: Not applicable- there are no human subjects.
16. Declaration of LLM usage
Question: Does the paper describe the usage of LLMs if it is important, original, or non-standard component of the core methods this research? Note that the LLM is used only for writing, editing, formatting purposes and does not impact the core methodology, scientific rigorousness, or originality the research, declaration is not required.
Answer: [NA]
Justification: No LLMs used beyond standard proofreading assistance and math checking, which does not affect the scientific content.
17

---

Appendix
Contents
1
Introduction
2 Related Work and Motivation: Spectral Bias in Distribution Learning
3 Background
3.1 Score-based Diffusion Models
3.2
Gaussian Data and Optimal Denoiser
Learning in Diffusion Models with a Linear Denoiser
4.1
Diffusion learning as ridge regression
4.2
Weight Learning Dynamics of a Linear Denoiser
4.3
Sampling Dynamics during Training
5
Deep and Convolutional Extensions
5.1
Deeper linear network..
5.2 Linear convolutional network
6
Empirical Validation of the Theory in Practical Diffusion Model Training
6.1
Multi-Layer Perceptron (MLP)
6.2
Convolutional Neural Networks (CNNs)
Discussion
Extended Related works
21
Extended Results
23
B.1
Extended Visualization of Theoretical Results
23
B.1.1
Scaling curves of diffusion training (EDM) fully connected
23
B.1.2
Scaling curves of diffusion training (EDM) linear convolution
25
B.1.3
Scaling curves of flow matching
27
B.2
Extended Empirical Results.
28
B.2.1
MLP-UNet Gaussian training experiments.
28
B.2.2
MLP-UNet natural image training experiments
29
B.2.3
CNN-UNet training experiment
31
B.2.4
CNN architecture ablation experiments
38
C Detailed Derivation: General analysis
41
C.1
General Property of Linear Regression
C.I.1 Gaussian equivalence
C.2 General analysis of the gradient learning of linear predictor
42
18

---

C.2.1 Denoising as Ridge Regression
42
C.3 General Analysis of the Denoising Score Matching Objective
43
C.4 Gradient Structure of Other Variants of Loss
45
C.5 General Analysis of the Sampling ODE
48
C.6 KL Divergence Computation
50
Detailed Derivations for the One-Layer Linear Model
51
D.1 Zero-mean data: Exponential converging training dynamics
51
D.1.1 Discrete time Gradient descent dynamics
52
D.1.2 Special parametrization: residual connection
53
D.2 General non-centered distribution: Interaction of mean and covariance learning
54
D.2.1 Derivation of non-centered case
56
D.2.2 Derivation of low-dimensional interaction of mean and variance learning
57
D.3 Sampling ODE and Generated Distribution
58
E
Detailed Derivations for Two-Layer Symmetric Parameterization
60
E.1
Symmetric parametrization zero mean gradient dynamics
60
E.1.1 Simplifying assumption: orthogonal initialization of 9m
E.1.2 Beyond Aligned Initialization: qualitative analysis of off diagonal dynamics
64
E.2
Sampling ODE and Generated Distribution
65
Detailed Derivations for Deep Linear Network
68
F.1
Aligned assumption
68
F.2
Two layer linear network
70
F2.1
Special case: homogeneous initialization A = A2 (symmetric two layer case)
70
F2.2
General Case: general initialization (general two layer W = PQ)
70
F3
General deep linear network
71
F3.11
Deep Residual network
73
Detailed Derivations for Linear Convolutional Network
74
G.I
General set up
74
G.2 General analysis of sampling dynamics
74
G.3 Full width linear convolutional network
75
G.3.1 Training dynamics of full width linear convolutional network
75
G.3.2 Sampling dynamics of full width linear convolutional network
78
G.4 Local patch linear convolutional network
81
G.4.1 Training dynamics of patch linear convolutional
81
G.4.2 Sampling dynamics of patch linear convolutional net
85
G.5 Appendix: Useful math
85
H Detailed derivation of Flow Matching model
87
19

---

H.1 Solution to the flow matching sampling ODE with optimal linear solution
88
H.2 Learning dynamics of flow matching objective (single layer)
89
H.2.1
Interaction of weight learning and flow sampling
H.3
Learning dynamics of flow matching objective (two layers)
91
1
Detailed Experimental Procedure
94
LI
Computational Resources
1.2
MLP architecture inspired by UNet
94
1.3
EDM Loss Function...
95
1.4
Experiment 1: Diffusion Learning of High-dimensional Gaussian Data
97
1.41
Data Generation and Covariance Specification
97
1.4.2
Network Architecture and Training Sctup
97
1.4.3
Sampling and Trajectory Visualization
97
1.4.4
Covariance Evaluation in the True Eigenbasis
97
1.5
Experiment 2: Diffusion Learning of MNIST MLP
98
1.5.1
Data Preprocessing
98
1.5.2
Network Architecture and Training Setup
98
1.5.3
Sampling and Analysis
98
I6
Experiment 3: Diffusion learning of Image Datasets with EDM-style CNN UNet
98
I.7
Architectural Ablation: Diffusion Learning of Image Datasets with EDM-style CNN
ResNet
99
20

---

A Extended Related works
Beyond the closely related works reviewed in the main text, here we are some spiritually related lines of works that inspired ours.
Spectral effect in the sampling process of diffusion models Many works have observed that during the sampling process of diffusion models [22, 41]: low spatial frequency aspects of the sample (e.g. layout) were specified first in the denoiser, before the higher frequency ones (e.g. textures). This phenomenon has been understood through the natural statistics of images (e.g. power-law spectrum) [19] and theory of diffusion [11], and recently through the lens of stochastic localization [42]. Basically, low frequency aspects usually have higher variance, thus were later to be corrupted by noise, so earlier to be generated during sampling process.
In our current work, we extend this line of thought to consider the spectral effects on the training dynamics of diffusion models.
Inductive bias of deep networks There as been a rich history of studying the inductive bias or implicit regularization effect of deep neural network and gradient descent. Deep neural networks have been reported to tend to find low-rank solutions of the task [43], and deeper networks could find it difficult to learn higher-rank target functions. This finding has also been leveraged to facilitate low-rank solutions by over preerrrramttiiiinngggg linear operations in a deep networks (e.g. linear [44] or convolution [45] layers).
Implicit bias of convolutional neural networks When the neural network has convolutional structures in it, what kind of inductive bias or regularization effect does it bring to the function approximator?
People have attacked this by analyzing (deep) linear networks. [37] analyzed the inductive bias of gradient learning of the linear convolution network. In their case, the kernel is as wide as the signal and with circular boundary condition. thus convolution is equivalent to pointwise multiplication in Fourier space, which simplified the problem a lot. Then they can derive the learning dynamics of each Fourier mode. This result can be unified with other linear network approaches [18).
[46] further analyzed the inductive bias of the linear convolutional network with non-trivial local ker- nel size (neither pointwise nor full image) and multiple channels, and provided analytical statements about the inductive bias. However, they also found less success for closed form solutions for even two-layer convolutional networks with finite kernel width.
From an algebraic and geometric perspective, [47, 48] have analyzed the geometry of the function space of the deep linear convolutional network, which is equivalent to the space of polynomials that can be factorized into shorter polynomials.
Deep image prior and spectral bias in CNN On the empirical side, one intriguing result comes from the famous Deep Image Prior (DIP) experiment of Ulyanov et al. [49]. They showed that if deep convolutional network (e.g. UNet) is used to regress a noisy image as target with pure noise as input, then when employ early stopping in the optimization, the neural network will produce a cleaner image, thus performing denoising. To understand this method, [50] showed empirical evidence that deep convolutional networks tend to fit the low spatial frequency aspect of the data first. Thus, given the different spectral signature of natural images and noise, networks will fit the natural image before the noise. As a corollary, they showed that if the noise has more low frequency components, then neural network will fail to denoise those low frequency corruptions from image.
People have also looked at the inductive bias of untrained convolutional neural networks. Theoreti- cally, [51] showed that Shiite--wdth convolutional network at initialization is equivalent to spatial Gaussian process (random field), and the authors used this Bayesian perspective to understand the Deep image prior.
We noticed that this line of works in deep image prior has intriguing conceptual connection to our current work, i.e. the spectral bias of learning a function with convolutional architecture tend to learn lower frequency aspect first. Comparing diffusion models to DIP, diffusion models regress clean images from many randomly sampled noisy images; on the contrary DIP regress the clean images on single noise pattern.
21

---

Neural Tangent Kernel A widely recognized technique for analyzing the learning dynamics of deep neural network the neural tangent kernel. For example, an infinitely wide network would be similar to a kernel machine, where the learning dynamics will be linearized and reduce to exponential convergence along different eigenmode of the tangent kernel.
Using neural tangent kernel (NTK) techniques, by inspecting the eigenvalues of the NTK associated with functions of different frequency, [52] has been able show that given uniform data sphere assumption, and simple neural network architectures (two-layer fully connected network with ReLU nonlinearity), neural networks learn gever-freueecyy functions faster, with learning speed quadrati- cally related to the frequency. Later they lifted the spatial uniformity assumption [53], and derived how convergence speed and eigenvalues depend on the local data density.
These insights have been leveraged classification problems to show that early stopping can lead neural networks to learn smoother functions, thus being robust to labeling noise [54].
What about convolutional architecture? With some similar NTK techniques, using simplified architecture, [55] proved that the learning dynamics of the convolutional network will preferably learn the lower spatial frequency aspect of target image first. Their proof technique also based on the relationship between verrparametrrzed neural network and the tangent kernel. The proof based on a simpler generator architecture: one convolutional layer with ReLU nonlinearity and fixed readout vector. They numerically showed the same effect for deeper architectures. This result provided further theoretical foundation for the Deep Image Prior.
22

---

B Extended Results
B.1 Extended Visualization of Theoretical Results
B.1.1
Scaling curves of diff'usion training (EDM) fully connected
A.
10
B.
10
10-1
Weight Scaling Factor
Sgmam1, 0,0.5
Generated variance
10
2017
4
10-1
4=3 3
-0.3
10-
4 0.1
0.03
A 0.001
10
One layer linear denoiser
Q,0.5 0,00,0 a,-80
10:4
101
100k
10
C. 30a
D.
10t2
Weight Scalling Factor
sigma=1, 000.5
Generated variance
10r1
10-3
4 -1
-.3
10r
4.0003
As 0.01
10-4
0.000
10
Symmetric two-layer denoiser
9,s.5, 0,00.002,
10
10-2 10-1 100
102
10-2
10-2 10-1
10
107
Training time
Training time qT
Figure Learning dynamics of the weight and variance of the generated distribution per cigenmode (continued) Top Single layer linear denoiser. Bottom Symmetric two-layer denoiser.
A.C. Learning dynamics of u W (())ug, B.D. Learning dynamics of the variance of the generated distribution Ag. as function of the variance of the target eigenmode d This case with larger amplitude weight initialization Qx 0.5.
23

---

A
the and mo lyye thear devober Vpriece es Ene Sky with ded Uhesl signot OOE2 SorUT -0
cee syer RJ--60 60
layer RUK- 0.
liger 0LX-003
ane lapar CLK 2.
J'M ligh Q,1
Ba
10-
Conpent Sep
Byeer 0G
0,K 0.
15
Deuiiaally Lai
B.
15"
13- 12-
18
Byer 02-
Tee Spprr QLE-
beo C.
Berrsce 1a
17
Emargasee
340
Coorgrss
Figure 6: Power law relationship between mode emergence time and target mode vari- ance for one-layer and two-layer linear denoisers. Panels (A) and (B) respectively plot the Mode variance against the Emergence Step for different values of weight initialization Qs E (0.0, 0, 0.5, 0.6, 1.0) (columns), for one layer and two layer linear denoser (rows). We used To = 0.002 and or 80. The emergence steps were quantified via different criterions, via harmonic mean in A, and geometric mean in B. Within each panel, red markers and lines denote the modes where their variance increases; blue markers and lines denote modes that "decrease" their variance. The solid lines show squares fits on log-log scale, giving rise to the y = a 2 type relation. Comparisons reveal a systematic power-law decay of variance with respect to the Emergence Step under both the harmonic-mean and geometric-mean definitions. Note, the Q+ was empty since zero initialization is an (unstable) fixed point, thus it will not converge.
and two layer case
24

---

B.1.2 Scaling curves of diffusion training (EDM) linear convolution
Method We performed numerical simulations to study the learning dynamics of basch-based convolutional filters. Using the FFHQ32 dataset as an example, we first extracted local patch statisties to compute the optimal linear filter at each noise scale w Then, applying Proposition 5.4, we obtained the learned filter at different training times and solved the corresponding PF-ODE in the spatial domain to generate samples. From these samples, we computed their variance along eigenmodes and analyzed the resulting scaling laws.
Key factors influencing scaling behavior. Our simulations revealed two critical factors affecting observed power-law relationship between image eigenvalues and convergence time: the filter patch size P and the initialization of filter weights.
When the filter weights were initialized as an identity matrix (without scaling), the initial sample variance became excessively large, and the resulting power-law exponent was around -0.4-00-mme -0.4-already weaker than that predicted by the linear-case solution. A more practical alternative initialize the filter weights scaled identity:
0 data
Wo
1,
σ² data
inspired by the skip-connection coefficient Cskip in the EDM preconditioning scheme. Under this initialization, observed both increasing and decreasing variance modes during training.
Effect of patch size. The specific scaling relations depend strongly on the patch size P. When P decreases, the scaling exponent for the decreasing modes approaches zero zero- from about 00.4 for P 15 to about -0.21 for p = 3 (3 x 3 convolution). The overall convergence time also increases for smaller patch sizes. For the modes with increasing variance, smaller patches further attenuate the scaling exponent toward 0.0 (P 7) or even flip the trend (0.48 for P = 3). Intuitively, a small 3 x convolution couples many frequency components Fourier space, so multiple modes are learned simultaneously during generation.
Summary. Although we longer have analytical expression for the scaling coefficients, the numerical simulations allow us to predict the scaling relations empirically. These results indicate that the patch size in convolutional architectures significantly modulates the scaling law. We hypothesize that the attenuated or flat scaling behavior observed in practical UNet architectures trained on natural images may, part, stem from their reliance on small (3 convolutional filters.
Table 2: Spectral convergence time computed for patch linear convolutional networks FFHQ32 at varying patch size.
Patch size
Increasing scaling
Decreasing scaling
3
0.06A 0.48
11.88A-0.22
5
0.05 0.24
3.66A - 0.36
7
0.05 0.01
3.24X -0.37
11
0.07 0.22
2.61A -0.10
15
0.08 0.28
2.57A -0.10
25
0.09X 0.30
2.54A 00.00
25

---

Inepry_ juive Jinrany TDW Irerond cikin.ont
Ineery rurn, linrony P5, EDN Precons csib, only
theoey, curee Innony, FOM Precand cski only
1990
2330
X00
2066
2030
290
2506
0.gg
747
0g
1212
127250
Takiing
Tekling
Thaning
p1op progyy.urvee inconv, LeLtoo__P__EDMM PLL.EOM Precand cok p.orly
enegry Jincon, EI.,EOM, Precord, covioddooyoooryy
thopry.cuvve_ Jncomy, F25, FDM Precond, ciil only
1500
1300
2200
2030
2003
2a7
0ig
1200
120
45
10000
17t
Dy
Taiing Dop
Takkiing SHp
Thaning Step
theo y curve iicom P3. EDN, Precond, cakip, only
theory, suwve Jiccorv ps EDM Precond, cakip only
theory. cun.e.ueeeee..eee Encony, P7 EDN, Precond cskip,ccciiccc only
87
57. 0 10 (1:0
R7 D
1271
10-1
12-4
decrase
Hntreaso
13
10-
inergesce Sep
Energence Step
Emnege Stop
theery curve Incoov, P11 EOM Precond rakip only Emergence Step
E-wery. curve Inconv P15 EON Precond eslin pnly vorance Emergence LinEs
theory curve Incomv, pas FOM Precoed eskip only Variance Emergonce with
127
Ine
Dicreass
12:
2211
3321
J17
317
20-1
10
10-
10
Figure Numerical simulation of spectral learning curve for linear convolutional network with local patch filter with patch size P p= 3, 5,5,7 7, 11, 15, 25, using FFHQ32 dataset. Top same format as Fig. 3C Bottom same format as Fig.3D.
26

---

B.1.3 Scaling curves of flow matching
Ona Laper Fow Matching
along Eiganmades
Generated Veriance eleng Eigerredes
Tuerggrerr
Figure 8: Learning dynamics of the weight and variance of the generated distribution per eigenmode, for one layer linear flow matching model Similar plotting format as Fig. 2. A. Learning dynamics of weights u (T; for various time point E {0.1, 0.3, 0.9, 0.99). B. Learning dynamics of the variance of the generated distribution Åς. as function of the variance of the target eigenmode A E {10, 3, 1, 0.3, 0.1, 0.03, 0.01, 0.001). Weight initialization is set at Qk = 0.1 for every mode.
A.
byper knear raaciing variarce scuuregeecc stop LIM
(re
Corperts
metcting Onergence Shep D.
B.
Q
2
Coree palte Ia T0
Figure 9: Power law relationship between mode emergence time and target mode variance for one-layer linear flow matching. Panels (A) and (B) respectively plot the Mode variance against the Emergence Step for different values of weight initialization Qk E {0.0, 0.1, 0.5, 0.6, 1.0) (columns), for one layer linear flow model. The emergence steps were quantified via different criterions, via harmonic mean in A, and geometric mean in B. We used the same plotting format as in Fig. 6. Comparisons reveal systematic power-law decay of variance with respect to the Emergence Step under both the harmonic-mean and geometric-mean definitions.
27

---

B.2 Extended Empirical Results
B.2.1 MLP-UNet Gaussian training experiments
A. Our unce thue a Eits B. Oygarral Conatiance Vatri Tue tigeriuers
Catzoder
Varaanc
Stip
byp 2211
0g2 153
+050
Incrsse 1e 172.08
0ytt 653
Toai
A. Variance of on true eigenters Goussty 250
2 B. B.Diggnnl of Covadincee Honk The Eigenbass
Saundan, Fitwed Lines
Verience
085 160
-:44
12 (iodl)
-0:300011
0530-002
16295 -03
12
350
127
Triiing
Diegonal of Covuriance Thue Eigebaasa
Variance Emergencc with Used Linsa
1gi-11
Egt 21
1414- 1m
1220
10
1496
1246 022
1y505
450
38 151
37
10
Draggne Stor
Figure 10: Spectral Learning Dynamics of MLP-UNet (Gaussian-rotated) (same layout and analysis procedure main Fig. 3) Top, middle, bottom show cases for 128d, 256d and 512d Gaussian. A. Evolution of sample variance A (r) across eigenmodes during training. B. Heatmap of variance trajectories along all eigenmodes, with dots marking mode emergence times T* ffirst-passae time at the geometric mean of initial and final variances). The gray zone (0.5-2x target variance) indicates modes starting too close to their target, causing unreliable 7* estimates. C. Power-law scaling of r* target variance dg, excluding gray-zone eigenmodes for stability.
28

---

B.2.2 MLP-UNet natural image training experiments
FFHQ32 UNet MLP EDM BL 3072D Irle-4 Covariance Eigenframe Alignment Score
AFHO32 UNet, MLP EDM 8L 3072D Irle-4
Covarlance Eigenframe Alignment Score
0.9
0.8
0.8
Coveriance anterancee Alignment
0.7-
Covarincee Alignmmnt
0.7
0.5
0.4
0.4
0.3
0.3
10"
10f
10a
Step 10a 10a
10a
10
10a
10
10a
10"
10a
Step
Figure 11: Dynamical alignment onto the covariance eigenframe of data (MLP-UNet, FFHQ32, AFHQ32). Alignment score function of training step. Alignment score defined the sum of square of diagonal entries of the rotated sample covariance the training data eigenframe U UTEUU divided by the sum of square of all entries. This quantifies how well the training data eigenframe diagonalizes the generated sample covariance. It will be X = if U is the eigenbasis Σ,
29

---

Variance along data sigenbasis
Evolution aigan projected variance
Target variance step
20
2
1
4
2
7
-010
409
8
5
3
7
0
7
1
151
13
0
15"
127
Varfance along data eigenbasis
Evolution eigon projected variance
Target variance emergence stop
ICX
0615 0.05
0ggs 0.6
1877
208
20g
107
12'
AFRI322
Variance along data caw siganbasis
Evolution projected variance
Target variance amergence stop
1,00
10"
F2 982 (N-0001
2540
#aaaa2-daaerdd Ranerd
8. Verlance along data cav elgenbants
Evalution elgen projected vartance
Target vartence vs emergence step
6.25
1981 0/90
0285 0.0
280
unter 6.-20x7 ofine
Bay
FFhg word
B. Varlance along data elganbasis
Evolution aigan projected variance
Target variance emergance stop
12"
10.
Forreso 124-1
Agg
Decresss
21:
11P
127
Figure 12: Spectral Learning Dynamics of MLP-UNet (MNIST, CIFARIO, AFHQ32, FFHQ32- fixword, random word). A. Generated samples during training. B. Evolution of sample variance (T) across eigenmodes during training. C. Heatmap of variance trajectories along all eigenmodes, with dots marking mode emergence times T* pfisst-passage time at the geometric mean of initial and final variances). The gray zone (0.5-2x target variance) indicates modes starting too close to their target, causing unreliable r" estimates. D. Power-law scaling of 7" versus target variance X separate law was fit for modes with increasing and decreasing variance, excluding the middle gray-zone eigenmodes for stability.
30

---

B.2.3 CNN-UNet training experiment
Variance sling data digenbesh
siges projected
Target variance
PD Ricdsitedit)
10
Figure 13: Spectral Learning Dynamics of CNN-UNet (FFHQ32). (same layout and analysis procedure as main Fig. 3) A. Generated samples during training. Evolution of sample variance (T) across eigenmodes during training. C. Heatmap variance trajectories along all eigenmodes. with dots marking mode emergence times T" (first-passage time at the geometric mean of initial and final variances). The gray zone (0.5-2x target variance) indicates modes starting too close to their target, causing unreliable 7" estimates. D. Power-law scaling of r" versus target variance λς. excluding gray-zone eigenmodes for stability.
FFHQ32 UNet CNN EDM 4blocks wide128 attn pilot fixednorm Covariance Eigenframe Alignment Score
0.9
0.8
Covariance Alignment
0.6
0.5
0.3
10a
10
Step
Figure 14: Dynamical alignment onto the covariance eigenframe of data (CNN-UNet, FFHQ32).
Alignment score r as function of training step. Same analysis as 11
31

---

Datio seed
CMN Dencoer
Uneer Denober (Pw
(Log Stale)
(Log
320
120
F00 Stwel
10t
Teat 34ne O12 UU-Aoo
NDC bo 1422 Staret
104
TEaR Stoo 00 SCCNOO
100
0:03
100
105
Norhh
Evaluation
Evaluation
0.12
10
13
100
100
105 Euioare Log V466
107
StoRES
101
SH005 SCORE
10455
Lewd5
102 twwaarce Loog
5.0
53.0
700
0.08
103/0 S000T
20:44
20:3
6:07
10F
107
Patch Sum
Sue
Figure 15: UNet denoiser can be approximated by linear convolution early in training (CNN- UNet, FFHQ32). A. Early in training, the UNet denoiser output can be well approximated by linear convolutional layer, with patch size P. B. The approximation error as a function patch size P. training time and noise scale Generally, early in training, the denoiser is very local and linear, well approximated by a linear convolutional layer.
32

---

Training time
10
21
39
1690
3156
S894
11007
20555
38386
71686
133571
249999
100.0
50.0
10.0
Noise
1.0
0.1
0.01
0.001
Figure 16: Visualizing the denoiser training dynamics with fixed image and noise seed (CNN- UNet, FFHQ32). D(x + oz, .) as function of training time т and noise scale σ.
Variance of learned samples on true eigenbasis MNIST
Eig0 - 25.78
Eig10
Vanance (normalized by target variance)
10
Eig20 2.49
Eg300 1.51
0.99
Eig50 0.65
Eg60
Eig70 0.35
Eig80
10-1
Eig00 0.20
10
10a
10
Treining step
Variance learned samples eigenbasis MNIST
Eg0
Eig50 = 0.65
Variance [narmalized target variance)
Egg000
Eig100 0.06
Egg000 0.03
10a
Eig250
Egg000 0.01
Cg3500
10a
Eig400 0.00
Eig500 0.00
10-1
10
10P
10a
10a
102
10
Training step
Figure 17: Spectral Bias in Whole Image of CNN learning MNIST Training dynamics of sample (whole image) variance along cigenbasis of training set, normalized by target variance. Upper 0-100 eigen modes, Lower 0-500 eigenmodes.
33

---

Variance of lermned parches 13x3, sridee 1) eigenbasis| NNIS
Veiaae to Eneegenns Tired putch stede
09: 0.30
10a
Eg2 0.00
0g2 0.26
>>
Eg4 0.03
Egg 0.63
Egg 0.01
10
Cg7 0.01
Egg 0.00
3012
anerge hamici
Varlance leamed patches 18x6, stride-3 true eigenbasis
Variance Time patch stride
Cgo
Eg4-00
Eg
107
10
10a
10
10
10a
Traning step
Vartance reamned 19x9, stride-
Virtence is Emergence Time I patch Sod strid
37
=g:8=001
>>
10a
0.01
Eg34 0.01
E0s3 0.00
E27 0.00
10
18'
decrers
Figure 18: Spectral Bias in CNN-Based Diffusion Learning: Variance Dynamics in Image Patches MNIST (32 pixel resolution). Left, Raw variance of generated patches along true eigenbases during training. Right, Scaling relationship between the target variance of eigenmode versus mode emergence time (harmonic mean criterion). Each row corresponds to different patch size and stride used for extracting patches from images.
34

---

Variance learned patches [3x3, stride-1] 1) true aigenbbaiis [CIRAR
10
Bgg 5.22
Eig 0.2-0.2
0.00
10'
0.00
R1221 0.00
Big24 0.50
Made enerpee in dae) samoer relen
Veriance learned patches 16x6, stridem 13 on true eigenbasis | CHNH
emergence pastch
Rg0 16.34
0ig 1.61
0ig 0.47
Ri2i 0.04
0224 0.03
10%
10
Veriance learned patches 16x6, 2tide-31 3) on true elgenbasis | CIFAR
Venance vo Emergencc Time putch Be33 strito
Bi0-16 16.75
Rgg 1661
0g0 0.27
Dig 0.22
Rg12
0.004
R1224 0.03
Posse
30-1
10a
18%
Figure 19: Spectral Bias in Based Diffusion Learning: Variance Dynamics in Image Patches CIFARIO (32 pixel resolution). Left, Raw variance of generated patches along true eigenbases during training. Right, Scaling relationship between the target variance of eigenmode versus mode emergence time (harmonic mean criterion). Each row corresponds to a different patch size and stride used for extracting patches from images.
35

---

Variance learned patches (3x3, stride-1) 1) on true eigenbasts FFH064
Verlance Emergence Time patch abride | ATTA
10'
tap0 6.00
10a
10
fig21
11g24
20"4
30
decent
dlytia
10
30
1.7
Thaining stop
ane garc stop bamaaiic
Variance of learned patches (6x6. stride on true ei eigebbasis FFHO64
Varlance Emargence Time patch abrida ATTHCG1
1100 20.74
Eig 1.42
10a
1ige 0.34
rigo 0.15
tkg12 - 0.12
tig1 0.07
fig1a 0.25
Fig21
Eip24
36"4
Avreau
20
Seceun
320
307
1.7
30a
Traning stop
Variance of learned patches (9x9) stride- - on true eigenbasis H0006
Veluano Driapence Timel patch
Eig0 29 99-95
01000 -0 0.43
1i500 = 0.14
tig50 0.05
t1440 0.04
31g50 0.03
f1go0 0.02
t1000
51000
20
300 301
1.2
30
Figure 20: Spectral Bias in CNN-Based Diffusion Learning: Variance Dynamics in Image Patches I FFHQ (64 pixel resolution). Left, Raw variance of generated patches along true eigenbases during training. Right, Scaling relationship between the target variance of eigenmode versus mode emergence time (harmonic mean criterion). Each row corresponds to different patch size and stride used for extracting patches from images.
36

---

patches stridem eigenbasis AHH06
Variene linel potch ded stride CAMONA
10t
0.01
c 0.00
g22 0.00
E224 - 0.00
30
30"
301
Mode enepprcs rop the nammotic reir
Variance learned patches (6x6, stride true eiganbasis AFH064
bmergence time potch ArHROe
Cgc
37.6:
cg2
0.0
0g9
Fg12
102
10-1
10a
30'
15
Thaining step
Variance of icamned patches (9x8, strid- 31 on true eigenbasis | APHO64
Egg 33.37
CL0 U.A.
14'
g20
g30
040
220
0.04
0000
E07
Brrnas
10P
10
122
10a
Thaiing
Figure 21: Spectral Bias in CNN-Based Diffusion Learning: Variance Dynamics in Image Patches AFHQQ2 (64 pixel resolution). Left, Raw variance of generated patches along true eigenbases during training. Right, Scaling relationship between the target variance of eigenmode versus mode emergence time (harmonic mean criterion). Each row corresponds to a different patch size and stride used for extracting patches from images.
37

---

B.2.4 CNN architecture ablation experiments
Dipgonal Covertance Hartis in True genbass
Diagzne of Covariance franie This Dieentdds
Deggnnl of Covattrre Natnz Thie Elpencatle
1003
3230
X5:3
2230
2003
2130
3553
3330
3000
5853
20g
Deg
Fhininn Zop
voring Step
Taanng Sopp
Disginal Covelienaa
Disgonel Matric Bigenbeels
Diegenal Covntiacce Natrs pingoo Rashue Llayers wioe256
Frroox Bestert EDM. Slayers wide22
0n
rfric Reabeer Eall Slapers
330
500
553
3300
1O00
10:3
1300
L500
12:3
2330
X00
2003
1930
246
2523
3330
3003
13gg
Leap
Tris
Sanas
Deg
Valy
Ladge
50000
STDLE
Spectral Learning Dynamics of CNN-ResNet (FFHQ32) with architectural variation 1 layer, channel numebr variation, ch 8, 16, 32, 128, 256. (same analysis procedure and format as main Fig. 3C.) Each panel describes one architecture, Top Sample final images Bottom Heatmap of variance trajectories along all eigenmodes, with dots marking mode emergence times (first-passage time at the geometric mean of initial and final variances). The gray zone (0.5-2x target variance) indicates modes starting too close their target, causing unreliable r* estimates.
Figure 22:
38

---

EDN, Llayers,
FPHO2Z, RustNeL CAN, layers,
Jayarr,
ancee Emergence fitted
Fiiied Lines
Emergence rittee Lines
10
30
10
Incenase
1M-1130
10
5:02 19-117
0:55 1120
Decrisse 1138
witer
107
1075
1077
1077
6400.50
Ver withe 2.52 orine
157
142
18"
10"
12"
107
100
101
Emergeece
Towegerce
Eraangenc
FPMQ32 Rashert CNN Llavers Emergence
CAN Slayars,
PFHQ33 Resket
Llavera
Veriance
Fifted
ribed Unees
Increase
decrian
12'
10"
cacreso
HLC 15-0N-
I92 0:45 1le
FS 110
Vic. thin
100
12-3
10'12
10'0
124 withs 2.57a oflnt
100
107
327
237
107
12P
132
327
100
107
Cenrgesce Stop
Erevrrnce Stra
Esensgence Strp
Figure 23: Spectral Learning Dynamics of CNN-ResNet (FFHQ32) with architectural variation layer, channel numebr variation, ch 4, 8, 16, 32, 128, 256. (same analysis procedure and format main Fig. 3D.) Each panel describes one architecture. Power-law scaling of 7 versus target variance Åg excluding gray-zone eigenmodes for stability.
Diagonal of Covariance Matrix in True Eigenbesis 2layers _DM_2layers
Diagonal of Covariance Matrix in True Eigenbasis FFHO32_ ResNet_CMM CNM_FDM] Bww___MM_2yayy 2layers_uvee wide128
FFHO32 Res.eeesNsii CNNN_DCM EDM
wide6
500
500
IncOx
0
lighartamee Oe!
1500
1500
2002
2000
2502
2530
3000
3000
-10
324
1263
2004
4548
10op
32016
S0500
324
2004
4548
12096
32016
50000
baining
rennnn
Figure 24: Spectral Learning Dynamics of CNN-ResNet (FFHQ32) with architectural variation
layer, channel numebr variation ch = 6, 128. (same analysis procedure and format as Fig. 22)
39

---

Diagonal
Bipenbabl
Coverince
Coratincc
2PH222 Restett ChN SDN Slayers
FFH052 Nsshit CMN PDM Slayera Moc1l8
MPM012 Ner ChR EDM Jlayers widl200
220
2:4
390
200
Alay
S000
Play
100g
Figure 25: Spectral Learning Dynamics of CNN-ResNet (FFHQ32) with architectural variation layer, channel numebr variation ch = 6, 128, 256. (same analysis procedure and format as Fig. 22)
Diagonal of Covariance Matrix in True Eigenbasis
FfHQ32 RTHQ32_ResNNe__C CNN_EDM 5layers_ wides
500
1000
Verieneee
1500
2000
2500
3000
324
1268
2004
4548
14090
32616
50000
Training step
Figure 26: Spectral Learning Dynamics of CNN-ResNet (FFHQ32) with architectural variation
-5 layer, channel numebr variation ch 6. (same analysis procedure and format Fig. 22)
40

---

C Detailed Derivation: General analysis
C.1 General Property of Linear Regression
C.1.1 Gaussian equivalence
Lemma C.1. For a general linear regression problem, where x, y come from arbitrary joint distribution p(x.y) y) with finite moments,
L= - + -yl?
then its optimal solution and gradient only depend on the first two moments xf
Proof. Let the error be Wx + b
then
θ
VwC C
Ex.y (eTe)
OW
-
[Ex y0-y)xx y)x []
WE xx | + bE x | - EE yx 1)
=2
W(E + (zH1 + b (Eyz + Hyft
= 2 W( + 1zl1 T) + (b 13y)1'= Eny
2( WWW - Enn) - 22 Mr + b ky)
2(WEE Enyr VnLuz
② B[eTe)
VDC
2E (e)
2 [W b y
2(W b Myy
We used the fact that
E x = F4g
E.
E(y]
E [yx
Eyg + Hyff
xx Eex + Pa P
Setting the gradient zero, get optimal values
W Eyz E-1
(11)
ba = My W W' pzz
(12)
The gradient flow dynamics read
d
W 2 WWW Eyy) V
(13)
dr
d
b 2m (n(WWpp pr +b Myy
(14)
dr
The Eyz determines the gradient flow dynamics and convergence rate of W, while 5-15 mines the target level or optimal solution of the regression.
2g Spn deter-
D
41

---

Remark C.2. For all the loss variants with linear denoiser, the loss is of this class. We can write down the gradient structure and optimal values W and b by plugging in the mean and covariance of input x and predicting target y.
Lemma C.3. For independent random variables X, RR have
CovlaX+brr, + + +dY = acLix + bdViy
(15)
Proof. This can be proved using bilinearity of covariance and independence which entails Exy 0. Cov(aX + bY, cX + dY) aCov( X, c.. + dY) + bCov(Y, dY) acCov(X,, X) adCov(X, Y) boCov(Y, X) bdCov(Y,
ac Σx + bd Ey
Remark C.4. For diffusion training loss, the clean image samples and the noise patterns designed to be sampled independently. Thus we can use (15) to calculate the covariance of input and outputs.
Using these two lemmas, the gradient structure and learning dynamics of various loss functions can be easily read off.
C.2 General analysis of the gradient learning of linear predictor
C.2.1
Denoising Ridge Regression
Lemma C.5 (Diffusion learning as Ridge regression). For linear denoisers D(x, Wx + at the full batch limit, the denoising score matching loss (2) is equivalent to the following Ridge regression loss.
L.
Fx
Wx +
x|| 2 W
(16)
Proof.
C. Exxpo. Do(x + 02;0) - x|2
-Ex-po=-- (..
W(x + oz) +b x|2
(W(x + oz) + b x) (W oz) + b - x) ]
Ex N Tr
Ex -~ (0.1) Tr
(W(x +2) +b x) 0z)
W +bb - -x]]
Ex 10.1)
x-bbbx++o22@ W +b -x- +oz w
Exxpo. N~ (0.1) Tr [Wxx + b b- (Wx+b-xx 20 Wz (xx b - x) T
+ o W
Excopo Wx b b- x|| Bx To 22 b b-x) + +o WT]
Ex Wx b b-xl x|2 + 2
This derivation depends on the linearity of the denoiser and the white Gaussian nature of noise ~ I).
Remark C.6. This equivalence reveals that for diffusion models the additive white noise can be viewed L2 regularization the weights of the auto-encoding regression problem, where o² functions as regularization strength Thus per classic analysis of Ridge regression, we can see that higher noise scales, fewer modes of the data will be "resolved"".
Further, the noise zero-mean but not white, the final term Tr( Wzz W will impose other types of weighted regularization the weights.
42

---

C.3 General Analysis of the Denoising Score Matching Objective
To simplify, we first consider a fixed σ, ignoring the nonlinear dependency on it. Consider our score, denoiser function approximators linear function (x;;o) = -b+ Wx. Expanding the per-sample DSM loss, we get
[xx + oz; 0) Xol
= b W(x + oz) - Xol
= +|b Woz + (W- - T)xoll
(b + Woz)' _ (b + Woz) + X (W- - 1)² (W -]xo + 2(b + Woz) (W - -]xx
=b"b + 20b Wz + 0 o²z WT Wz + x (W -
(W)-I"(VV I)"(W=--=-))) I)x0
+2b (W- -))x, 2oz W (W -I)x0
-bb + 2ob Tob:WZz+ W oTr WT Wzz" ) + Tr((W - 1)"(W 1)xox01
+ 2b (W- - -I)x0 + Tr W (W - I))02²)
||b -xxl|l 2 + W + oz) 2 + 2W (x1 + o2) (b (b)
Full batch limit Here we take the full-batch expectation over z, Xo, where z NNO0,) I). xo p(xo). Their moments are the following.
E 2 =
E zz T1 = I
E [xx]] 1
xxx]] = 1qp T + Σ
E [xxx2 0
Note that the data do need be Gaussian, as long the expectations or the first two moments are as computed above, the results should be the same. In sense, if our function approximator is linear, then we just care about the Gaussian part of data.
C = E x(xx)(x)) ((0,) D + oz; xol12
=b o Tr Tr((W I)((W- ()p Σ}] + 2b 2bT(w 1) (17)
We can simplify the objective by completing the square.
= =Exo<p(xo). 2~ DDx, + o2: 0) - xoll ||b + (W- - -1)=1l μ. (W - (W - П) TT + Tr ()"(W - D) + Σ)]
b [I W)ul + W| + Tr((W ( 1|Σ
(18)
|b [I w) TT Σ)] 2T +
(19)
Gradient The gradients from loss to the weight and bias read
VbL 2(b - W))p
(20)
Vw = -2 (W E) + [2Wpp + 2(b - m)p²1
(21)
Global optimum Examining the quadratic loss, and setting the gradient to zero, can the optimal parameters are the following,
b' = -I-W-*)
(22)
w E -)
(23)
which recovers the Gaussian denoiser.
43

---

Gradient flow dynamics Next, let's examine the learning dynamics by gradient descent. We consider the continuous time limit, ie. gradient flow. Denoting the learning rate r) and continuous training time т, their gradient flow dynamics are
d
b VD
(24)
dr
d
= -WVW
(25)
dr
44

---

C.4
Gradient Structure of Other Variants of Loss
Here we generalize our analysis many popular training losses of diffusion models, including flow matching. We provide a detailed table for different variants and the learning dynamics of the linear denoiser.
Table 3: Variants of the diffusion training loss and their linear solution and gradient structure. Here Egz denotes the covariance of the input to the linear predictor; Eyz is the covariance between target and input, weight projected onto principal component ug. " _ W un; 1/r" represents the convergence speed of a target mode. We abbreviated away the expectation over data and noise in our loss
be confused with the covariance of training sample Xo. wi denotes the optimal
EDM
X-pred
EPS-pred
V-pred
Flow matching
Loss
Da(x +xtt
xl'
Fo(cox [1e:l]
Falax (.6: ell
Fo(
det.
lage
oxx
ua (1 ( 1x:
[x]
Exz
Dto1
a.E+oaa 0,
Y
Y.)
63.
f)
W
( (o 1
ariof2 + off) 1
mrla + off)
0,01 of 'I1 2)
122 + [1 - 1)
112 - (1 1)1)
oft0f2 + ofl) "
-0,0+10;2 Folf) 4
o0 + di [1) +ofl)
(1 1]
1)
14--11-11
Aztoz Ag t
u2AL rot
1Paa 1p
1/5*
of +
oe +
of Xg +
10g (1 1)=
Denoiser / clean image x0 prediction (EDM) loss
[Da(x (00x+00;;0] x|1
Moments of input-output
P1 14, My p
Σxx Cov(x + xc, (C) +0001
Eyz Cov(x, + σE) =
Optimum
W -E-
E yx = 2 02
b' -1y W Ux
= (1 -W")
=
Noise eps prediction (EDM) loss
Fx e|l
Moments of input-output
Pi ,1, By =0
Σ Cov(o + xc, 0C) E 021
Σ Cov(C, σc) = ol
Optimum
W = Σ Ei = 02
b Mly W* Mr
= --
= =
Denoiser / clean image / x0 prediction loss (variance preserving)
[Do(asx + t00;0) x||
Moments of oututt
Pr ast, My =
Σ,, Cov(aix + oie, 04 01e) = afE of'
Σn Cov(x,,,x 0.6) =
45

---

Optimum
W En Lazz au(afE+
b" = My He
1 04 W
-1-- 11-ax -a, a:W"),
=(I [I-ofI(afI+o0 "pp
= of(a8E+0711)
Noise eps prediction loss (variance preserving)
F9((ax +0,;t) e|l
Moments of input-output
0it4, My = 0
Σzz Cov(aE + ofe, 04x t a t,
Myz Cov(E, axx + ott) = o,11
Optimum
W E-1E E = =
b By W Pr
-aa 4
Velocity V-prediction loss (variance preserving)
+Folaxx+01e) 0,e - (aye ox
Moments of input-output
Mv aiHe, My -014
Σx Cov(a; 04x o a of
Eng Cov(are - -00, 0,x + + -1000E + a = a000 - --)
Optimum
W = Σ-Σ Lyz
04T (0 + T2 Σ)
b" My W* μs
-OP a4 W*H
- + a, W /
E))1
+
-oi a 0 (1 - Σ |
ot
(of 2 1 a/ Σ) (f
(a +o) 0;1)
Flow matching loss
C = Bxx~N(O.1), xi | ((1 - t)xo + txi, t) - (x1 - x
46

---

Moments of input-output
Px = 114, Myy = 14
Fg Cov((1 - t)xo + fx1, (1 t)xo + (xi) = = + (1 t)²1
Eye Cov(x1 Xo. (1 t)xo + tx1) = - (1 -)
Optimum
W Σ-1Σ
12 Lye
= (PΣ + +11 (1- (1-1))
b" By Wh
= 11 WitH
(T -(W")p
(1-(E+(1-1)))--000 (1 (1--)) u
(1-00-10-1-(1--11--101--- (+ + (1 (1- ((--t)))--1 "
(
-(1-))(=++ (1-))))~~
47

---

C.5
General Analysis of the Sampling ODE
The diffusion sampling process per probability flow ODE is the following
dx
OS(xx,0) σ))
(26)
do D(x,0)-x
given the relation between score and denoiser (Tweedie's formula).
D(x, -
s(x, o)
g²
Intuitively, when the denoiser linear function of x, then the sampling ODE is also linear (time-varying) dynamic system with respect to
When all W are jointly diagonalizable by shared set of orthonormal bases {us} } (e.g., the PC basis of data, or Fourier basis), i.e., commute, we can solve the sampling dynamics mode by mode by projecting onto such basis.
d do
x WaX + bb, -x
(27)
d
1 u W. + uf b, - uI. x)
do u x =
(28)
Representing W., b,, x(o) on the shared eigenbasis,
Wa Σ Cu(l))uuk
(29)
by Σ bu(o)ux
(30)
x(o) = Σ cx(o)uu
(31)
We can project the sampling ODE onto eigenbasis, i.e.
d - ((Vb(() 1)ufx + ba(a))
u
do
d 1 - ((1r(0) 1)cx(o) + by (a))
do (a)
do (())) + (h)(a) )ox(() = - -)b((o)
d
σ
This can be solved via the general solution first-order linear ODE: given
dy dr + P(x)y Q(x)
the general solution reads
y = - P(A) [f cs (e))
In our case,
Cx(o) = - 1 Mx(A)-1 [f bx(A) p (e(e)-1 de dx+p
Sej 12( dA
1" ba()) Me(e)-1 de
cx(or)+ +
Jor
e 11r 11(1)11
(x( -1 (e(()- de da
(o()T) +
48

---

Rewriting the solution to expose the linear dependency on initial values, the solution reads,
(x( = As(o; + (x(o;oT)
(32)
As(o;oT) Saz 2a(Ai-1 ax
(33)
Ba(o;OT)= / (A) e- 13 valei-1 de dd
(34)
Consider the function
Φx(a) = exp 10 (x(()) - ( di
(35)
Then the integration functions can be expressed as
Ax( = (x(o)/
(36)
(a( = [" (x
(2) da
(37)
or
By initial noise distribution, (x(or) (0,, the variance of (a) any sampling time can be written down
Varlcx(a)] = o7 exp ( - ( (2) -11
(38)
or A
$1(
5T (
(39)
Since [[cx(or)] = O,
B[cx(o)] Soy (x(2)-1 dd
ba(A) dloz vx(e)- dx
(40)
bb()) -J² (x(0)-) de dd
(41)
ba(A) $x(()
dA
(42)
(A)
Jor
oa(a;or)
(43)
Since is linear transformation of Gaussian random variable, the distribution of x, at any sampling time is also Gaussian Ni, Σ) with the following mean and covariance.
j Σ Be(o;or)uk
(44)
4x(0) uu D (OT)
Σ = a2
(45)
k
49

---

C.6 KL Divergence Computation
KL divergence between two multivariate Gaussian distributions is,
(L
E1)
[[ Σz))
1
Σ
1
1 2 (2-22) Ex² (2: ue) p(x)dr
log
(x p1) 2 ' m1) +
2
Σ
2
Σ
1
log
-tt tr
[E[( 1x)(x [1)²² Mj Bl(( vz)²Zig (or 22)
Σ
2
2
Σ>
1
2 (m )22) E2 (W f 12) + 1 و زn}
1
log
tr (Ia)
+
Σ1
2
2
Σ
tr Σ |n} + (12 m1) Eg (#2 m1)
log
d +
(46)
Σ1
This formula further simplifies when Σg and Σn share the same eigenbasis. We can write their cigendecomposition Σ₂ = U'gUT EΙ UAAU".
1
A2l
11 E2))
log
d + tr( A1 + (12 m1) UN, U U mi)
2
A1
1
log
II
A2,1
d +
Σ
A1. + (12 mn) UAA u (μ2 m1)
2
I X1,h
22,1
k
1
log
22,1
d+
Σ
dI,k + (12 UAA U TT(2 p2)
2
A1,k
d2
k
In more explicit form
KI E1) (122, E2))
1
Σ
log
A2,
d+
Σ
d1,
+ (112 111)"²
Σ
unuu
(12 mm1)
2
A1,k
22,
22,
k
b.
1
d2,
dI.
(u2 (112 - 111)
log
1
(47)
2
A1,k
22,
222
If they share the same mean μ41 M22 it simplifies even further
A1.k
Σ
A1.k
KL = 2
log
(48)
X2,
X2.
h
which has unique minimizer when AL. Ag.4
1.
Thus, we can compute the contribution to KL mode by mode. We denote the contribution from mode l as
A1. A X2,
1
A1A 22,6
KL.
log
1)
(49)
Thus for Gaussian data that share the same mean, the KL divergence can be reduced to the ratio of generated and true variance along each principal axis of data.
50

---

D Detailed Derivations for the One-Layer Linear Model
D.1 Zero-mean data: Exponential converging training dynamies
If μl = 0, then the problem reduces to learning the covariance of data. The gradient to weights and bias decouples
VbC = 2b
(50)
W = + 2 Σ)
(51)
The solution of b is an exponential decay
b = -2nb
(52)
dr
b(r) = bd exp(-2m1)
(53)
The solution of W can be obtained by projecting onto eigenbasis of Σ
VwL U2 22Eu 2W E)
(54)
- 2 + Wuz
(55)
Thus,
d (Wu;) -1/1-22pu4 + 2(0 + (wue)
dr
Define variable Vı Wu
d
Vk = 2nd -2n(0² + +dx]vE
dr
dk
= 2n(0² + ANa)
uu -VA
(02 + Aa)
d
Vx(T)= = (o2 + Ax)
uk + +v8
us) exp ( - 227 + (6))))
(02 + Ns)
Since
[vc_]] W(mn...]
VU W
W Σ V
k
The full solution of W reads
W(T) Σ (x(T
Σ 2 + Ax)
d
u T + Σ (ve
d
uz) up -20/02+
ux
(o2 Ax)
Σ
d
ukit T (1 -e---m 2n(o 2 +del7) Σ vuu
+ A1e)
Σ
d
(1 ynsi
2n +Xx)T) W(O) Σ ynun c 2
2 + Axo)
W + (v [a2 A
us uuuuk 2n(o++An)
Ax)
W' + Σ (W(O)u: (02 A
us uu 22y(0++
(56)
xx)
W +
- -W 2n(o²+Au)r
(57)
where vl. = W
For notational clarity, we derive the gradient flow at fixed noise scale a, and omit the subscript and/or argument from W. b. and L.
51

---

Remarks
The weight matrix W converges to the final target mode by mode.
The deviation along each eigen dimension decays at different rates depending the eigenvalue.
The deviation on eigen mode uy has the time constant (2n(o²² + Aa)) -1 i.e. the larger eigen dimensions will be learned faster.
While the "non-resolved" dimensions will learn at the same speed ~ (2ma²))²
Comparing across noise scale a, the larger noise scales will be learned faster.
Score estimation error dynamics Consider a target quantity of interest, i.e. difference of the score approximator from the true score.
First, under the μ = assumption, we have
1 E, = Ex|s - s (x)11 = 0²
||b + Tr[( - w e'D]]
(58)
The deviations can be expressed as
- 6" = 10 exp(p2211 -21TT)
W W Σ (vD
d 2o + Ac)
uk)) u - 2n(o++AuuT
with the initial projection v2 = W(O)uz
E 04
1
| 1 exp( 4n77) Σ + Ac)llv
Au us 22 1on(nnn=+0 +2x)r
(59)
(02 + 2/)
1 |1061| exp(-4n7)) Σ (1)+01)110111 -An(o
(60)
d (a2 Ax)
&. W(O)ue
u1k
6, = 60
This provides us with the exact formula for error decay during training.
Denoiser estimation error dynamics
ED = Ex)) (x))1
= o'E. E.
(d6)² exxp(-Axp()) E + As)]|oul1z - +2(()++A)T
(61)
Training loss dynamics Under μ = 0 assumption, the training loss is basically the true denoiser estimation error plus constant term の² Tr Σ) resolvent).
determined by data covariance (trace of
-u=0 =lbb (I - W) + Tr((W- - (' ( - W + - -)=
(μ =11b11 Tr W W*) Σ) T ー
(62)
EE + o'T -)
(63)
D.1.1 Discrete time Gradient descent dynamics
When the dynamics is discrete-time gradient descent instead of gradient flow, we have,
VbC= = 2b
(64)
Vw = -- + 2 Σ)
(65)
52

---

The GD update equation reads, with 17 the learning rate
b - b, = b,
(66)
W4+1 W. =,,
(67)
D -2n)bb E + W E))
(68)
W41 W. 2n(
(69)
=2n + W.(I 2nj - 21/Ei)
(70)
= =2n + Wi((1 - -2n0²)I - 2nj
(71)
For the weight dynamics, we have
W up = -2n + W. - -210022 21)20) un
(72)
= =2 + Waux (1- 2no0² 27da)
(73)
This iteration is exponentially converging to the fixed point g=+AE
A uu = (2ndx
dk )uk- + WaUR (1 2n10² - -212z)
W W
(74)
02 + Ak
02 + Ax
du ut + Waus (1 - 12100² - -212a)
= =22(lo² + +A6) - 1)
(75)
0² + A
d
=W
u1 -210² 21))2s)
(76)
0² + Ak
Thus
A
Wau: σ² + As
uu + (Wou:
uz) 2n1² 20/2aa'
(77)
02 + Aτ
b, bb - 27))
(78)
W Σ 2 2 Luuu + (Wo-
Ak
As un un (1- - 21102 - -21/20))
(79)
の² + A
k
So, there is no significant change from the continuous-time time version.
D.1.2 Special parametrization: residual connection
Consider special parametrization of weights
W = Cakip I + Cout W
(80)
It's easy to derive the dynamics of the new variables via chain rule,
oW
@W
Couts
(81)
HL Cout
(82)
OW
aW
With the original gradient
VoC = 2(b - (I W)
(83)
VwL = -2E - (o -
Σ) + 2Wppp" + 2(b - =)
(84)
and zero mean case
VbC= 2b
(85)
VW = -2 + 2 Σ)
(86)
53

---

the gradient to new parameters
V V = Cout V WWL L = 2Cout (-E + (Cakip + Cout W (O + Σ))
(87)
2oout (-Y + Cokip (o' 2) + Cout W (o'I Σ))
(88)
2Cout (Cokip (o"I+ Σ) Σ Cout W Wo Σ))
(89)
dW' W.
990)
dr
dW
(cskiplo21 - E Cout (o E))
Cout
(91)
2ncoot dr
W
=
1
(Σ Cskiplo E))oo += Σ)
(92)
Cout
1
(I(O'+ ++)===--- -Cckipp
(93)
Cout
Wr() = (0) (xpp(-2/T@nnl + Az))+
(94)
A Cakiplo? + Ax))
(1 eepp( (xmp(-217cean cou + +x))))
Cout (022 + )x)
The solution to the new weights reads
- - -W")))"") Σ un
+x)) W
(95)
W(7) Cakip + Coat (T)
(96)
Cakipl + Cout (W((0) w") Σ un exp( 1xp(-211m) exp(-227e2 2n
+ Ax)) Cout W
(97)
E(o + Cout (0) w") Σ Tn expl 21 + Ax))
(98)
= -W* + Come (W - W Σ uh exp( 21 (ou + As))
(99)
= WW + +W(() - W Σ un of
+ >x))
(100)
The only difference is scaling the learning rate by factor of e Also potentially depending on whether we choose to initialize (()) or W'(0) from a fixed distribution, we would get different initial values for the dynamics.
D.2 General non-centered distribution: Interaction of mean and covariance learning
Summary of results for non-centered case When μ * O, the gradients to W and b become entangled (see Eq. 5), resulting coupled linear dynamic system as follows.
Proposition D.1 (Learning dynamics of linear denoiser, non centered case). Gradient (Eg. GF) equivalent to the following ODE, with redefined dynamic variables, VI((T) = W(T)us, b(T) b(r)= J.. Denote overlap my := uTI15.
V1
V1
A U
d
V2
V2
u
=
M
+
(101)
2n dr
5
b
0
54

---

0 A=
m= = 00000,=111 1.5, 1.0
m= 1.0, 0= 0, 1.0
Velues W. 3nd b
Values of W. aeaddd
and
62
Time tyyl
Time Inj
Time imf
-24
Figure 27: Interaction of mean and covariance learning. Top solution to the w, b dynamics under different noise level a (0.1, 1.5,4). Bottom Phase portraits corresponding the two-d system. (m 1, 1, =
with fixed dynamic matrix M defined by & Kronecker product.
2 A1 + mi
Ma
M
m1m2
A2 + m2
m2
& I
m2
M @ La
(102)
Remark D.2. The dynamics matrix M has rank-one plus diagonal structure, specifically, it is the diagonal dynamics matrix in Eq. 7, perturbed by the outer product of the overlap vector my- The eigenvalues of such matrix can be efficiently solved by numerical algebra, with eigenvectors expressed by resort to numerical simulation and oow-dimensional examples to gain further insights. We can see the coupling of W and b dynamics comes from the overlap of mean and principal component my u./u When certain eigenmode has no overlap, my = 0, the corresponding weight projection Vi+(7) will follow the same dynamics as the zero-mean case, i.c. exponentially converge to the optimal solution ur) When the overlap non-zero my + 0, it will induce interaction between V (T( and b(T) and non-monotonic dynamics.
formula [56, 57]. Without closed-form formula, we have to
Av/ (Ay + o²) ug.
Two dimensional example Here, we show a low-dimensional example illustrating the interaction between the bias and one eigenmode in the weight. Consider the case where distribution mean μ lies on the direction of a PC uy. Then only the uN mode of weights interacts with the distribution mean, resulting two-dimensional linear system, parametrized by noise scale a, variance of mode À and amount of alignment m. Let the dynamic variable be scalars we b, Vε w(T)uk b = b(T)uR.
w b
m2 + 2
+
m
W
A m +
(103)
2n dr
1
m
The phase diagram and dynamics depending on the noise scale are shown (Fig. 27): At larger values, the dynamics of weights w will be much faster than b, basically, w gets dynamically captured by b, while b slowly relaxes to the optimal value. At small values, the dynamics timescale of w and b will be closer to each other, and will usually have non-monotonic transient dynamics. When 2² and À are comparable, w will have mon-monotonic dynamics.
55

---

D.2.1 Derivation of non-centered case
In this full case, the dynamics of b, W are coupled
VbC = 2(b + Wp)
(104)
Vw -28 2W Σ) 2(b - Wplp
(105)
--2 2W V n²²
(106)
VwC+ VnC n²T
(107)
VwL == -2E + 2W E
(108)
The nonlinear gradient learning dynamics reads
b VC
(109)
W -VVWC VbC 1²))
(110)
W - b. N²= -wyww
(111)
=2n(I W(o + Σ)]
(112)
b =-
(113)
-221(b - w
(114)
Consider the projection
(W b. (=)u = 21/ W(O Σ) up
(115)
Wu: (n1" ub)b 2n[dEu: (0² + Ax Wue]
(116)
b= -2nbb + Σ uf
k
Consider the variables vx(T) = W(T)u let B(T) = -b(T)= - so now the dynamic variables are (Vi...., b)
VA (ni 1b 2n//xu4 (0² + vv)vol
(117)
--nnn - Σ (n un)vE)
(118)
Vi = 21 - (0² + Ao) - 2n(u² uw))b+++ + E
(119)
= =2n(Auus (0² +doJVR - (Hit'x) - ((d ui)vd
(120)
b = 22n + Σ(2 (HFFu))V))
(121)
The whole dynamics is linear and solvable, but now the dynamics in each component vą becomes entangled with other components Vi.
V1
V1
V2
V2
Aa 112
2
+ 211
(122)
dr
b
0
with the blocks in M matrix defined as follows
Muk (o² + d + (uTux))
MRI (nTuw Ju)
Mrs uI
Mn (1it [=
Mas = I
56

---

We denote the overlap between the mean and principal components as = uru
us.
Then we can represent the dynamic matrix M as tensor product of dense symmetric matrix Q with identity Iv- More specifically, the dense matrix Q is a diagonal matrix plus the outer product of vector. So it's real symmetric and diagonalizable.
0² + A1 m
MI 22 0²
mi
M
m2
+
A2 + m2
m2
x
IN
(123)
i
m1
12
Q x IN
(124)
x² + X1 + m
m. 122
m1
Q =
m2
0² + X2 + m
m2
(125)
m1
m2
i
D + qq
(126)
== [m1, m2: 1|_
(127)
D diag(o² 0² + 02, + >a.... Ov
(128)
Note the inverse of Q is analytical, but the general eigendecomposition of it is not. Since the dynamic matrix M is real symmetrie, the dynamics will still be separable along each eigenmode of M and converge w.r.t. its own eigenvalue. The eigen decomposition of M can be obtained by numerical analysis and eigenvectors from Bunch-Nielss -ieesee-S-o-nnneeen Sorensen formula [56, 58]. Generally speaking, since the mean of dataset i usually lies in the directions of higher eigenvalues, the dynamics of b and the top eigenspace will be entangled with each other.
As a take home message, the overlap of μ and spectrum of Gaussian will induce some complex dynamics of bias and weight matrix along these modes. The full dynamics of W, is still linear and solvable, but since the dynamic matrix is a tensor product of a Diagonal + low rank with identity, a closed form solution is generally harder, we can still obtain numerical solution of the dynamics easily.
D.2.2 Derivation of low-dimensional interaction of mean and variance learning
To gain intuition into how the mean and variance learning happens, consider the Id distribution case, which shares the same math as the multi dimensional case where the mean overlaps with only one eigenmode
V&L- = 2(b / W1))
24-2( 2(1 w)j
2pw + 2b 2/
(129)
VWL -2E + 2W Σ) + 2(b n W
= 2(1 w)y² 2w + X) + 2pb
=2(p² 0² A)w + 2pb 2() + n²)
(130)
Write down the dynamic equation as matrix equation
d dr
22 (1
12 0²
P 1 b ]
A x+12)
Eigen equation reads
dct A ( = (1 - (12 + 0² + 1)))-- - 11²
(131)
= 72 () + 12 +0² + -
(132)
Generally for 2x2 matrices,
[ b
57

---

their eigenvalues are
1 2
1
( (a c)2 + 452
1 (a ala-ci c)2 4-40 aa-c
U12
25
In our case, = 11² + o2 + 14 c= 1.
(122
1)2 2
4A-1
U12
24
The faster learning dimension will be A1 u1, where w and b will move in the same direction.
Key observation that the w and b's dynamics depend on o,
for larger o, w will converge faster, while will slowly meander, moving with b, following entrainment.
X + 122 jub A + 422 o²
w" w²(b)
for smaller o, b will converge faster, comparable entrained by
b ("(w) = (1- (1-w))) w)
D.3
Sampling ODE and Generated Distribution
For simplicity consider the zero-mean case, where
A ux uu)uue -20(02 (o2 As)
o(T:o)= W + Σ
(W(O;
(133)
k
Σ ( + Xxe)
A
un - 22((0++20) + W(0;o) Σ us Yuuue 2
(134)
b(T:0) = b(0) exp( -21T)
(135)
To let it decompose mode by mode in the sampling ODE, assume aligned initialization u W(O; o); = when h + m.
Then
d do
u x ((-60
σ² AR) +u W(O:))u: (02 + Ax)
d
) 20(0 2 uf ore 2nr
Let the initialization along up be uf W(0 = Qsa then
d u x ( ((-Tott) An)
(1x
d
)e -2n(0++2e)
up + d) 271
do
(02 + Ax)
Using the following integration results
/ do 0 + Ax)
02 = log (Ax + + C
2
[ do _ 2o-l0- 2n(o2+Aeir - 2nd Ei (-21TTo2) C
1
d 2n +
do a(0² + >x))
(Ei (-2nTo2) E 22
Ei (-217 (02
+ +N+))) +C
2
58

---

Integrating this ODE, we get
cx(o) -
- exp
1/2 log (Me + +o2) + - Bi (-211o0²) e 2n Ei E-(-2n (o² + XR))) + (136)
1
1 Ei (-2nTo") e -ndreT
(137)
2
=C [((1- (11-g)) Ei e - Ei (-2nT (o² + +xx))) ) (138)
E --2
Solution of sampling dynamics ODE
ck(ov)
x(oo)= Σ up (((TT)
((T)(
(139)
Σ ck(oo) n x (oToT)
(140)
The variance of generated distribution reads
Cà (00) (((
j² uug Up
(141)
where the variance along eigenvector u, reads
X = o (10o))
(142)
(o(TT)
(1 (1-9:) qx:) Ei (-21TTo)) e 2
Ei (-217 (o) + As))
exp
d + 0
TT
(143)
A + exp
2
E (1 - 9k) ( (-21TTo7) e 2nT - Ei (- (o} + Ax))
59

---

E Detailed Derivations for Two-Layer Symmetric Parameterization
Here outline the main derivation steps for the two-layer symmetric case:
D(x) P x + b.
(144)
E.1 Symmetric parametrization zero mean gradient dynamics
Note that if the weight matrix W has internal structure, i.e. parametrized by 0, we can easily derive the gradient flow of those parameters using the chain rule.
Let W = W(O),
ac
ow
VOL
Ow
00
Here, when W PPT, via the chain rule, can derive its gradient, which depends the symmetrized gradient of W
VpL (VW[((VC)PP) (VwL)*P
(145)
= =[ P
(146)
where
VbC 2(b (I W)-
(147)
VwC = -2 - 2W Σ) + WWpp" + 2(b - m)F²)
(148)
Expanding the full gradient (non-zero mean case), we have
vpL -E Σ) + 12 2(b 1)}²
TT + E 22p TT 2u(b n17i| P
(149)
EEP + 2PPTT 2 P 20 I 1
+ [2PPT up + 2(b n) + 224 PPT + 2u(b 11
(150)
In the zero-mean case this simplifies
P -EP + 2PP" (0-I + P + 2 Pp
(151)
Consider representing the gradient on eigenbases, let u P = 9%
u/. PP = -4u, EP + 2 PPT (o'I+ E) + 2u (o I + E) PPTT
(152)
= -4 P 2u >(0²
(02 + Am P un u P + 2( + uoJuf PPTP
(153)
m
-44A T 2 Σ [(0² Am 2 + > Σ am
(154)
= 44x + 2 Σ
(2o² + Am + As) (1T am Am
(155)
of Cp=0 = 4h + 2 Σ
(20² + Am + 2x) qm
(156)
m
V Cy=o -44: +2 Σ (20² Am >x) qm
(157)
Fixed points analysis A stationary solution at which the gradient vanishes
if * m,
q qm
A d
(158)
+ 2 or
if m.
Note, this is different from the one-layer case where there are no saddle points; here we get a bunch of zero solutions as saddle points.
60

---

Dynamics of the overlap
Note the dynamics of the overlap
d
9k T d
d 9m + 9m dr 9k
(gk 9m)
9m)
dr
dr
-nnlut Van Lu= + 4m V.
nn Gymol
+ Σ
(20² + An + X1 ( On ) tn
-4 (9k 9m + Σ
(20² + An Ama) (Im Iu )41 9nl
=-n(-A(AA= + Am ( +2
+
+ Σ (40² + 22n + Am + Ax) (n)(
= =4n((Ax + Am) Mm - (202+>n An Am + dB (1k on (m)
(159)
This shows that when all the overlaps initialized as zero, they will stay zero, i... when weights are initialized to be aligned to the eigenbasis, they will stay aligned. We will first solve the aligned case analytically in Sec. E.1.1, and then analyze the dynamics of overlap qualitatively in Sec. E.1.2
E.1.1 Simplifying assumption: orthogonal initialization of 4m
Consider the simple case where each of 9m 0, Vk + at network initialization, i.c., each & are orthogonal to each other. Then, it's easy show that ár (gk 9m) the start and throughout training. Thus, we know orthogonally initialized modes will evolve independently.
Note this assumption can also be written
T q 4m u PP""uu Um i. Wum. =0 Vk + m
which means the eigenvectors of Σ are still orthogonal w.r.. matrix W, i.c. the matrix W shares eigenbases with the data covariance E. In such case, the gradient reads
V q2 Cy= = - 4hbak + +2Σ (20² + Am + Ax) (9k qm )m
(160)
(ortho) 4haar + +2(20² + +2xx) (xoT:xhe 9x)
(161)
Adaak + +4(02 >x) 4x)
(162)
The dynamics read
dqp V ge fuz=
(163)
dr
= =-n(- 1-m(-4Axpp 4Auu + 4(02 + Ax) gx)gk)
(164)
= =4y(Ax - (o² + >x) (ukt 4x))
(165)
Since the right hand side is aligned with gu. can only move by scaling the initial value.
The fixed-point solution is gs 0 or when of qs
de_ Given the arbitrariness of qs itself, track the dynamics of its squared norm. The learning dynamics of o% 9u read
T dyk An - (02 + >x) qu) ( gk)
qk dr
(166)
d( of 9b) = An/Ax - -(0² + xxx) qx)) (oT gk)
(167)
2 dr
d(f)
= (n - (02 fx
(168)
dr
61

---

Fortunately, this ODE has closed-form solution.
(1h qs (T)
(1)x()))1
(169)
1/Ke at +
8
(170)
1/ Ke BTAUT + 88(022 Xk)
8T/18
(171)
8nA [0x
8n(02 + )x)))E STAATT + 8n(0² + Xt)
(172)
Ag
(02 + 2ι)) ) BTAUT + +22+
(0)
1
(173)
0² + A
1)E -0dee +1 (10:000112 o²+hx
1
(g (T) = 0
if i + m
(174)
This gives rise to the full vector solution
(0)
2n (T) = v (T]1
(175)
11 (0)
d
2 (0)
(176)
02 + Ak
An - 11e--nnde STAATT +1 V (1ox(0)12
[9)(0)]
a2
d
(0)
(177)
02 + Ak
de
(0){()))) + 11 (0)]
d
1
(0)
(178)
22 + d
A SndeT + (1 - -- |lun (0)112
Learning Dynamics of score estimator %
Now, reconstruct the whole estimator, W = PPT Recall
u.P P
P Σ u
P
Σ 4 T
PPT
U )
Σ
Im um
k
m
Σ
Ua (gh 4m J
b
m
Note that under our assumption,
(of 4m (T) = 0
if k + m, of qm is diagonal. So
PPf Σ
u ugla qx)uu
k
Σ
(T)((
Then rewrite
d
1 ao+ 1)e STidET +
[])] -cf (T)
(179)
02 + An
A
[0]=
Σ (T)
A
1
Σ 0² + x
u uu
(180)
0= - 1)e-
de
-8nee +1)
62

---

Score estimation error dynamics
E, = -x)))
04 |b - 10-bb"1++ + 2(b wi)u
(181)
+Tr((W - (W-))(--''(e + E + o'DI]
E(u=0)
b - ba1 Tri(W - - (=) onD]]
(182)
lb - b"1 Tr(( - - w > (de + ) U
(183)
Thus
E( cav term) terrem
Tr((W- HwW--w766 (W")"----WW -W") Σ( (Ae + ou)
(184)
=Trl (A+ (W-WWWWwuu (W--"" -W
(185)
d A + 02 uk T
Trr Σ (A+ + (Pqk
d + 0² ux) (P
(186)
= Tr Σ (AR + 2)
dx 0 u
d=
7(Pqk
T
d + + ut}
d
(187)
dk
1de + (
T. PT Pa:- P'
P4k + (
)2 11
(188)
2 + 0²
de + oz
d
dR
Ene + (1) 2
d + 02 the qu +
(189)
d
E(,+ + 0²) ((g qt)
d + )²
(190)
E(u==, bias term (60)² (a/9:xx((-00 (xp((4A)))
The dynamics can be applied
d
1
(gkh qm (T)
0² + d
de (0)))
1)e SppdeT +1
o+
Elu=0, con term)
Σ
(Ax + +022(
A+
1
A
)
02 +>x
-updr
+
d + 0²
[0]1"]]
a==Xe
(191)
Σ (Ax + +
A
)²(
1
) (192)
-1)²
0² A
A
BipueT +1
b
[0]10 ++
Σ d 02
X2
(Tx((0) 1 0.
1)e STAAT
2
(193)
dv
1le--81Ar- +
X2
(0= - (0)2) -indeT
de
2
(194)
A + 0²
Au
11 (0)(()))) + (0)12
The full score estimation loss is
Elu=D) (6 exp( -1]T)+
(195)
X2
A
STder
2
Σ do + 0²
(0))()) -8ndher + (0)122
(
63

---

E.1.2 Beyond Aligned Initialization: qualitative analysis of off diagonal dynamies
Next can write down the dynamics of the non-diagonal part of the weight, i.e. overlaps between Q and 4m
dr (9k 9m)
(196)
d
qm + 9m d
ak dr
dr 9k
(197)
= - 4(xx + Am) (1m (k)nk))+2 + Σ (40² + 2An + Am + >x) 4n (a) qn)
(198)
n
Am + de (gk 9n (9 a 9n)l Σ (202 + An
=nl(Axx Am Am (9
(199)
-47 (de Am((1 + Am Am)) 9x) (20² Ax +
Am Av
Nah qx) (1
(200)
2
Am + dR
+AR
-200² + Am
(af gm) (
Am d
Σ (20² An
(ok On )
4n
As + Am -2202
Am + +3x
)lax|1 2 - -(20²
3dm + dE 111gmll 112 (QTmen))) (201)
2
Am +
d
(20² + An +
(ak qn) (
2
ny
Am + An (llgxl + l1 (/mll3)) (aTh 9k
4n/
(Ax(1 11---e -114/-11) + Am
11 113] - (20² +
(202)
(20² + An + Am + dx
(ak 9n (m)
k.
As reference, recall the dynamics of diagonal term,
dg qx) 8n(A+ - -0" (c) qk) (g 9x) 9b) dr
Assume the diagonal term is not stuck at zero oi qu + O, then we know it will asymptotically go to [00)
o2 following sigmoidal dynamics.
Ae
For overlap between of 9m. without loss of generality, assume A Am. Then have three dynamic phases: 1) neither mode has emerged; 2) mode & has emerged, while mn has not; 3) both modes have emerged.
Phase 1: neither mode has emerged then the overlap will grow exponentially.
When neither mode has emerged, assume [|gx(T)] (T) = O.
di (ith dr
d
(m) gm) so An de + don (itten)) ex)
Note given \ > Am. the diagonal term has higher increasing speed than the non-diagonal overlap (x))))) > 4n(Xs Xm).
term (7 (o2 + (x qx))
64

---

Phase 2: one mode has emerged, the other has not dy but the other has not (A; > Am)
When one mode has converged [ ~
0x
d
ds + Am - -2002 2002+ Am + +32k
d
dr 19k qm))
4m
(203)
x²
+
Am
A
dR
47)
λ6 + Am - (20² + 22x +
(9m
(204)
2
0² + dx
Am
AK
41)
λ6 Am 2dA -
(KTT
(205)
2
0² + Ak
Am
Ak
Nk
41
Am
As -
(206)
2
0² +
Au
1 AB
4m
Aw)
1
(207)
2
+ An.
Since Av < 1, we know + oo < O, so the overlap will follow an exponential decay dynamics. Further, given the same As= the decay speed of qf 4m is proportional to the difference of the eigenvalues Xm d So, the larger the difference, the faster the decay.
Ag F
0. Thus, the dynamic coefficient (Am - Ax)
2
Phase 3: both modes have emerged
When both modes have emerged [Tq)())]12
o=+Ax
11 [T]l 20 o=
Am
d
Am + +31k
d
(20²+ + 3)m
AB
Am
dr (9h (m)) x 4n/ A1 + Am (20² +
( qt)
0² A
2
0² + Am
(208)
02
(Am Ax2²
4m
A5 Am (-) (02 Am (o² + 2o5)
(mTH
(209)
(AR
+ Am 20 (2AmAt + 0² (ds + xm))
-1 2 (Ax + σ2) (Am + 02)
(9 (k)
(210)
The dynamic coefficient negative, so the overlap decays even more rapidly.
Summary: qualitative description of off diagonal elements of 4m dynamics
The off-diagonal term will exponentially rise after the rise of the larger variance dimension, and before the rise of smaller variance dimension;
after one has risen it will decay to zero:
after both have risen it will decay faster.
Thus, it will follow non-monotonic rise and fall dynamics.
E.2 Sampling ODE and Generated Distribution
Here, for tractability purposes, we will focus on the zero-mean and aligned initialization case.
Recall the weight learning solutions above (180),
W(T:o)= PPP (7)
Σ 11 (T)
Σ
02 d u.uu
(1
-1 Ak - +1) 0² + A 1)e Sarde
(211)
with the weight initialization
W(0;))= = Σ (0)
(212)
65

---

Note here we assume the initialization 11 (0; is the same for all a levels, with no dependency.
liga (0; II (0)]2 Qx
Per our general analysis (Sec. C.5), the key factors are
dk dx(o:)) = 0² + 2 (on 02the
1
(213)
de
SOTTAK -1 +
and its integration
0x(0) - =xxp(- 1 ()(dd't)=1
do')
u x
1
4/ (o: - uf juf
-1)u]
(214)
do
1 d
-1 Ju{x
(215)
1) SerrA
[Tx1
Note that the integrand is just a fractional function of o² which can be integrated analytically.
I do 1 (( 0² + d (et 0h
d
de - 1)e 1)e-11m87 StT +1 -1)
1
UNEBE log (As + Qk (6 1) (Ax + +z²)) -2Qu log(()) + lo)
2 (c - 1) Qx +2
Qhe log (Ae + Qn (c - 1) (Aa + +)))) + +2(1 - Qk) log(a)
2 (EBNTThe 1)Qx+* +2
Qk log (AE Qs (C - 1) (Ax + 0²)) + 2(11 Qkk) log(or) -8007d
2(1 - 6 ) Q₁ + 2e -8
Qs logle (Aue-8nrd -p + 2x (1 e Serre) (A+ + +z))]] + 2(1 - Qx) log(a)e SarAe SoyrAe
2Q} + 2(1 10x)e-89r4
Qu log BaTRe + Qu (1 e 81FAA) ()x 0²)) + An 2(1 Qk) log(ore Snr
2Qp + 2(1 - Qx) -B
Thus, the general solution to the sampling is
(a) =C exp( Qu log (Ave Se -8p + Qu (1 - e Bn (Ax + o²)) + 8 + 2 - Qx) log(o)e -gr
2Qx + 2(1 2k) -8
2r
C Ase--BFDh BurAa + Qu (1 e BBrde) (Ax + 02)
20. -01-2 BNTAA
x
StJT 1xQ
exp (1 Qx)e -Byrde log(o)
exp 2Qu 2 - Qx)e -8 -
(216)
Q (1 2x)e--88p Ar
The scaling ratio is
(o)) (x(TT)
SuITA + Q+ (1 e 80 (d+ + oz)
20, +2 Q4) 14
(217)
Aue 8yF7, + Qx (1 80 ) (A& + of)
exp (1 Qale SyTAA (log(oo) - log(oT))
Qk + (1 - -x)e-=-)ap1 SyfTde
Are -
+ Q (1 - e e- (Ax + og)
10+ 20.+2(1-02)= Agr
dg
2c+(1 Q2) -9rAg
(1-0,) -B
Q (1 S0 (A1 + of)
(218)
66

---

(1-R Ao
(x(0)= = dae By + Qk (1 e Bardu) (2x 02)
24x +20
Qaie BorAx
+ Q2 BOTTA
(219)
Evolution of distribution
Following the derivation above,
[[x(o)]]= = Σ (oT cx(o0)
k cu (aT)
The generated variance along eigenvector up is
ck(ou)
{x(T) oT
cu(oT))
=o 2 A + Qu: (1 ByrAe) (Ax + oî) dGe Bn + +Q& (1 81 (hx oz.)
9x+(1- Qu.e B
201 20-22)2 Bo
4x + OxJo Bir
OT
(220)
Asymptotics of Learning: early and late training
At late training stage T/T + oc,
lim (o)) - (T(oT)
QK (As + σā)
30
Qu (Ax oz)
exp
k
de + o
V A + o
(221)
lim {x = o. \ + 400
(222)
A +
which approximately recovers the correct scaling factor to generate correct variance Ag. we take Ør 00, dp O.
At early training stage Tyr +
lim (k(oo) T- (OT)
1x11 + Q6 (1-1) (Ax + o
30
00 2
(1-0a)
dxl + Qk (1- (Ax o
0T
00 1-Q
0o 1-Qu
(223)
which shows the initial generated variance is determined by weight initialization scale and integration limits (oo, or).
67

---

Detailed Derivations for Deep Linear Network
Consider general deep linear network. W II. W, Wz WWWL-1-.W W1
Gradient structure
Using the chain rule
Oc dW si By
vyC=
oW is
we have the gradient with respect to matrix W
OW
(W L.. W441) (W.-1...W)) (W....
W 1.3
In vector notation,
0L
ow mn mn aWi 1,
V Fli
aw
OL OW mm
= = W W4t)
Wi)jjn
dL (WI-I..WW1)
Ww, = C (11)
OW
where Vw can be found in (5)
V&L = 2(h - (T W)
VwC = 2E +
2W +2( -1)p²²
F1 Aligned assumption
This gradient structure can be substantially simplified if the weight matrices aligned at initializa- tion.
Consider the singular value decomposition of each weight matrix W
W UA
and for each pair of neighboring layers, the singular modes are aligned.
v U-U
or equivalently U1-1 = Ví, vI E
Gradient structure under aligned assumption
Then the product of weight matrices
I
W.... WL..:W411 W1+1 UU
II
A
A=1+1
W B.L...WW...... U11 ART
W W. = UL
A
k=1
Then the evolution of hidden weights reads
Vw,C- (Ur II k=1+1
A&V T VW AT K=
1-1
Vi+ II
AQU.. V w LV, II
AAU
k=/+1
k= 1-1 Au ViT
L
Ui
A
k=/+1
68

---

Substituting the gradient of the loss (5), assuming centered data μ 0,
VWC - = 2
Then
Vw,-=Ui
II -kI+1
ARUL - + Σ) Vi
A
V?
k=1
L
L -2 20 II
1-1 x) V1
Ui
II -=1+1
(
ANV
Ak Vi
k=1
1-1
20
πI
ARU II A&V' + 2
II A
A
II
A&V?
k=1+1
k=1+1
k=1
20
II k=1+1
ANT
EV II Av?? k=
k=1
Under the simplifying assumption that U2 Vi and they exactly match the eigenbasis of Σ, denoted as U, then
VVEN= = A,
UVEV =A
1-1
L
1-1
1
1-1
Vw 2U
II k=1-1
A
II
V
+ 20
II k=1+1
A
A II
V + 20,
II
A
AgA II A&V
1=1
k=
k=1
k=++1
k=1
k=1
20ii + a²
A
A
II k=1
ARl
II
V
k=1,
-Ui + (02 + AJA
II
A
II k=1,kl
AVV?
1=1,
Consider the singular values of each weight matrix
V.C = U7 CV,
L
= (o-
Aal
A
Romn
I,
The weight dynamics have surprisingly simple mode-by-mode form. For the m-th mode, the gradient is
L II
Ar. L - -2(->m + (0² + Am,)Arm
Ag.m
II Akm
1em1 1.
Jom 1.
To simplify notation, we define Ei,m = I Nem- Then the gradient flow dynamics reads
d At. -1/l -mm (o² + Am At
II
A
II Au. k=1,4
dr
k=1,4
=-1/-AAm (o² + Am At, 1
At first glance, this is linear system for At.m with fixed point at Xm/(o² + >m))/Er, But the effect of other layers E,m also dynamic, forming coupled nonlinear system. Er,m affects both the time constant and convergence level.
69

---

F.2
Two layer linear network
We start by considering the two-layer case of this equation. Under the aligned initialization assump- tion.
Wi U
W2 U2A2V2 U2 = =U
(224)
U = 22
(225)
So the main degrees of freedom come from A1, A2, and we discuss the two cases depending whether they are equal to each other.
F221 Special case: homogeneous initialization A = A2 (symmetric two layer case)
If A A2, then
W UGAVT V2A2U1 W We are reduced to the symmetric two-layer case W = W W, we discussed above in E. Due to weight sharing, the gradients to each layer are accumulated.
VA, C = 2 + (o² + A2A2
= 2 + (0² + 1)
dr A1 2nl A + (o² + A)AHAA1
d (Ai) 4n| A + (o² A)A
This recovers the dynamics we obtained in the previous section 5.1.
F.2.2 General Case: general initialization (general two layer W PQ)
Generally, if the initialization is not homogeneous, A1 A2, the matrices are the same W PQ. and we have the general two-layer parametrization.
Under the aligned assumption F.1, the weight learning dynamics reduced to those of A1, For simplicity notation, let fo A1, Sk = A2,k- We have
V Aksk + (o² Ax)g2fk
Vg C = -Aff + (0² + xx)
which is a two-dimensional coupled system. Let us focus on the dynamics along one eigenmode k and drop the subscript from fk. Ik. Let the constants A nd.. B n1 + >xe), we have
Ag - Bg²
(226)
dr
dr 9 Af Bf
(227)
There are three important variables: fagk. f2 g.. and -2 This system can be represented in the following way, which exposes its conserved quantity:
d dr (f9) f 9 dt
dr
Af2 Bf g + Ag² Bg¹1
= =(f +z) (A Bfg)
+2g - d 1 3) 2f JT
d
d
dr
2{AAf - Bg2 + - -
2 2g (A Bfg)
dr (f- =
70

---

Field of System
25
Magnitude
0.
10
-10
-15
-2.0
Figure 28: Phase diagram for the simplified 2D dynamic system. Above: $ Ag By² g Aʃ Bf for A 1, 2. We the manifold of stable attractors in the 1st and 3rd quadrants: fg A/B, and the conserved quantity along the hyperbolic lines.
Thus the (J.g) pair can only flow along hyperbolic lines diagonal lines where f²- g22 const. We can leverage this fact to write down the overall solution.
Notice that
- -2)²
So we can represent
f + = 2 + (f2
+
Now the equation for Sg becomes closed:
dr (sg) CCA
Let h = fg
dr -h= h=v + C -h)
Note that for this self-contained equation, unless C = 0, it shall have only one attractive fixed point which is h" A/B. Thus asymptotically, will always converge to the correct value.
When we face the weight tying case where C = 0, the equation becomes
d
h 2 Bh) dr
It will have another fixed point at h = 0, which makes it impossible to converge to the fixed point
A/P B if the initialization h()) has the opposite sign.
E.3
General deep linear network
Weight tying assumption initialization special case. The key gradient equation is
The weight tying deep network can regarded as
VAC=I-A+
(02 AJA
II
A
A
k=Lk
1m
71

---

Let us assume all A, are equal at initialization'. Then
VA = 1-A + (o² + A)AFJA
(228)
Consider single element of the diagonal at mode k. dx (1) == A then we have
= 1 + (o² + Xx) (a ( (aa
(229)
Oa (1)
Further consider the aggregated effect Cj II. 0x variable reads
Assuming weight tying, the gradient flow of this
ack
L((x)) (1) L-1 da
(1) θT
θT
-1 11)) 1-1 7L
dak
i
-n [-> + (a² + >x=) (a "1
'1 (ak
(1)
2
21.-2
=-n 1 + (0² + A ckl C/
=- | + (0 + Au )ca] Cj
2-
We arrive the dynamics of the weights along eigenmode:
der θr = L - (o² + Ax)ckle 2-
(230)
Depth one and two cases Note this form also encompasses the solution of one-layer and two-layer symmetric linear models by setting L = 2, as in (??) and (??).
For general L, there is no closed-form solution; one only has an implicit solution to ck involving the hypergeometric function Fi:
General depth case
A Ac(T) L/L B
F 1, 7 1+ 2 Ac(T)) L' B
A AA((0) 2/L
2
2 Ac(0)
2
1, T
1
-BLBr,
B
L' B
(231)
with substitution A (o² + AB,B dp.
Infinite depth limit
The limiting case is L + oo, then the dynamics read
dc;
dr nn | (0² + An)ea ck
2-
(L + 00)= nL Dd (0² + An)cal of
(232)
Then we have
des
n dr
[ - (02 + Cx
0² + A
Ck
In
nLT + C.
X2
d (o2 + Au)c
AAck
Setting the initial values as c
(0)
we have the implicit solution of c::
A
Ck
1
nl o² A
In
(a²
ARIck
X
d
(02
)x
de
(0)
Axco (0)
where training time t can be expressed as a function of cb:
CA (At (0² + Ax)c(( (0)
1 o A In nL
1
(233)
x
(0) (Ax (02 + )ck)))
dcck
dxr (0)
Ck
? A is not be confused with A, which is the spectrum of data
72

---

F.3.1 Deep Residual network
If the network is parametrized with a residual connection at each layer, i.e., W = [I+ W;) (I+ W [I + W WL-1)...(L+ W)), then linear shift in parametrization. Using the same derivation and notation, assuming aligned initialization and homogeneous initialization, we can see
0C
AE= + (o2 + Ax)) (1 + 11+02) (1 + ) "
(234)
a
Let the overall effective weight be Ck = (1 + 0%
(1) 1. we have the Jacobian:
dck
= L (1 + 0x 1 6-1
)
0 k
Thus the dynamics of the effective weight is Ock
da
(1)
L(1 + 0%
(1)
L-
nn d (1) 10 C
-nL (1 ag
(1) 1-1
Ao= + (a Ax) (1+ [ (1 + a )
-n Ak + (a² 2 + Ax) (1 + a (10) | (1+ + a (1) ) 22-2
= nL A + (o² + Ax) ck
cx
(235)
We recover the same dynamics as the normal deep linear network:
dck Ōr TT(AAx (o² + +xx)ce 2-4
(236)
Here the initialization cx(0) (1 + x1x))) (0)²
Thus in linear networks, residual connections do not change the learning dynamics, but just change (shift) the initialization of weights.
73

---

G Detailed Derivations for Linear Convolutional Network
In this section, study linear convolutional denoiser. For simplicity, let's consider the space is 1d, with Id convolutional filter w. For simplicity, we ignore the bias term and assume zero-mean a 0. The convolutional denoiser is defined
oxx;)) = Wa
(237)
This convolution can be equivalently represented as matrix multiplication, where the weight matrix Wa is Toeplitz circulant (when circular boundary condition).
Dxx; W.x
In this case, the score learning problem becomes circulant Toeplitz regression problem, which has been studied before [59, 60, 61]. This problem is similar to finding the best convolutional or deconvolutional kernel to 1d sequence.
G.1 General set up
Cyclic weight matrix and spectral representation Let's consider the circular convolution case, where we ignore the boundary effect, and assume the matrix W(o) is cyclic at any noise scale. Since it's cyclic, they can diagonalized by discrete Fourier transform (DFT).
Let's 's diagonalize the weight matrix with the DFT matrix,
( = FT
specifically the DFT matrix is defined
-2 mk
Fmn exp
m, k:
N
G.2 General analysis of sampling dynamics
Note that, we assume weight matrix at every noise scale circulant, then we always can solve the sampling dynamics the Fourier basis, mode by mode. At its core, this is because all circulant matrices commute.
d Fx
do
[(o)F"x - F'x
F do
0
[(o) -1 FFxx
We can perform integration mode by mode. Let C₁ = =F**x)k, x Fc
d c
7k (0) -
Ck
do
Thus, if we know the Fourier parametrization of weights T(a), we can integrate the sampling of in closed form. We will need to solve these Fourier modes ((a) of the weights during training dynamics. Integrating the ODE, we get
@x(a) exp (ʃ
( da)
A
(o() (K((TT)
ca())) T (oT)
(238)
Then the generated variance will be
(p0) (p(
2
A of
(239)
5. = F diag dτ) pd
(240)
Because of this casily prove Proposition 5.2
74

---

Generated distributions from linear denoisers are Gaussian distributions, with covariance E. Since the generated covariance is diagonalized by discrete Fourier transform F, is circulant i.c. translation invariant. Given that the generated is Gaussian distributed, x conforms stationary Gaussian process. While the integration function (x( determines the covariance kernel of this Gaussian process.
Proof.
G.3
Full width linear convolutional network
G.3.1
Training dynamics of full width linear convolutional network
Gradient structure in spectral basis
Using the Fourier representation of the circulant weight ( F
we can derive the gradient to the Fourier parameters ((a) which is a diagonal matrix.
Recall that
VnL 2b
VwC 22 +
Then the gradient to the Fourier parametrization reads
(V W. (o)
=(F W
WF
Thus
VrC FL
-F + 2F'
= -- + 2F FTF" (f
= =F +
Under circular convolution assumption, T is a diagonal matrix, Fi = 6:y7 then
OC = V|
0%
+ + F'EF)
=F] + (b00 F'EF),
= 27, (FFFF))
The key entity the F-FF matrix, let's define it as Σ
Interpretation of F*F
We can regard it as covariance matrix in the spectral basis, i.c. covariance matrix of the complex variable f*x x E C..
= cov(x) = E E[xx?)
F"F = F'E[XX]] F
E F1
EX pT]
= (FF
-co((F'x)
75

---

Diagonal values tell us about the power spectrum of the data.
[F Var([F**] x]()
1
Fjs
jk/N
VN
[F
jkk
Σ
5mm Emn Funk
Emn e 122rim/// e -2mi nk/N
N
mn
Emn 20rá 2m
Dmn 2ar 20mmn) n)
F
N
mn
Gradient flow on spectral parameter
Given the gradient structure,
ac 2(F"EFF kR + 27k (o² + [F"
07k
25 + 27k (o² + E1A.)
if we directly perform gradient descent the 7 variable, we have
d
OL
dr 7 77 Onk
dr -1| 2②xx + 2 (o² + Ens)
Eas 7n (0² + Ens)
21
Fixed point
Fixed point of the gradient flow is
7h σ² EAE
(241)
Learning dynamics of Fourier modes
This is basically a first order dynamics ODE
(x(T) = EAR
ENH
+ (7)(() 0² Enx
e -2n(o=+Eue)r
0² + Exe
= =7 + (7)(0) - -2nlo=+
(242)
Remarks
This is exactly the same story the one layer linear case, where o² + ΣEs determines the convergence speed per mode.
Spectral modes Σn with higher power will converge faster, which usually corresponds to the lower spatial frequency modes.
Higher noise level will converge faster
Learning dynamics of weight entry as rotated Fourier mode dynamics
Lemma G.1 (Convolution Spectrum-entry relation). The relationship between spectral and entry parametrization of filter weights is
N
We
2mi 7k. 8=
O..,..NN-
(243)
k=0
76

---

Proof. Derivation Using our previous convention
Wij ij Wj-i
Let the vector the first row in W be Let the vector parameter be
W Ok w/
W (wg, W1, WN WN-1 1
= =Wor)
Then the connection between and weights we are the following
W FTF
W jk = Σ Fijm jm mn Pin
mmm Fjm Fink m
Fin Fin jmm F mk
1
7 2xr jm/N 2kink/NN
N
1
Σ
2 1mm
N
Thus
Wkk
W Wiitk
+ 2 & 7m
N
Fi Ym
N
In vector notation, we have
1 F
W
VN
7 NFw
Thus, they are related by a VN scaling and a unitary transform.
Similarly, the gradient to y and that of w has following relation
V = VNF-TTTT V
=VNF'V,C
Proof. Derivation
(V,C,dm) =VyG, VNFWw -(VNF-TV,,,. dw) =(V dw)
77

---

If we let the optimization parameter be the filters w, then we have gradient flow in w space
dr v
The corresponding gradient flow of y is
dr w V
d 1
F*y -
dr
VN d
dr 7 = -NNV
So gradient flow in space is equivalent gradient flow in 7, but scaling learning rate by 17 Nn) Conversely we treat y as optimization parameter and use gradient flow, then it's equivalent to scale learning rate - n/N
dr Y= - -
d
VNFW
71
FV VN
dr d w dr
7 V
Remark G.2.
Thus we can solve gradient flow any representation, but translating to solution of the other variable just by rescaling the learning rate.
This trick only works when the kernel width is N, the spectral parameter will be con- strained in a subspace, which will alter its dynamics.
Direct spectral parametrization without any constraint is basically equivalent to the entry wise parametrization. does not provide extra inductive bias.
Corollary G.3. For linear circular convolational denoiser, D(x, a) *x the solution to gradient flow on filter weight w is w(T, = VN" F o(T,o), where the spectral parameter of k-th Fourier mode evolves as
(T(T,0)= = ENe
(7) +²+ Σnh
(Nn(o"+$inn)
0² + Σnt
= (n) + ( 2N
(244)
where (p) = Eu o24Exs
and Σes. is the variance of Fourier mode.
G.3.2 Sampling dynamics of full width linear convolutional network
Let's project the sampling equation on Fourier modes,
F+ -' x I(o)
Fx
da
Let Ck = (F*x)k, x = Fc.
d
7n (a) -
do
The variance amplification factor is expressed through the integral,
{ = cxp ( A
--))(()) -
(x.(2) ax)
78

---

Generated distribution after convergence
For the converged denoiser we have 7/ o242p.
88nk
7x(A) -
( = exp
d2)
A
Euk
dd)
cxp
X - Eiek dd)
A
exp
12 In(2 + ENx) ) Amm
exp
02 + ENE
2 ΣEK
Variance of the learned distribution
(1(( (
2
07 00 Σnk E.K
X of
07. + Ens
The generated covariance will be
Σ = Fdiag FF
F.diaggiin (
Remark G.4.
Basically it's the same story as the one-layer linear network, the major differ- ence is, instead of evolution along the eigenbasis data covariance (PCs), the parametriza- tion of convolutional network enforces the learning dynamics to align with the Fourier basis regardless the original PC of the data. It treats the data as it's stationary translation invariant.
The generated data will be independent along each Fourier mode, just like the fully connected case, the generated data independent along each eigen-mode.
This amounts to decorrelating the original covariance in the spectral domain, making translation invariant / circulant.
Generated distribution during training Using the solution gradient flow of spectral parametrization, for each Fourier mode,
(T( 7/ + (7; -x)e 2n(o² BBu)
We can solve the sampling ODE during training. Assume the initial value of each Fourier mode the same across noise scale (0: o)) = (0), Vo. We can write down the generated variance through training time as,"
t (a) = exp (ʃ²
7h (A) - -
Euy x2+Ime
+ +(xx(0)=
Eyy A24
2
= exp
dd
A
= =VEn + o² exp
(1 (0)) Bi (i (-2nto?) e 2 Ei (-215 (o2 Ens )
Using the solution to gradient flow of filter weights parametrization, we effectively amplify the learning rate 77 Nn in the expression.
(T(r;)) = 7h (0; - 2NN
'One thing to note is that it's highly likely that the initial value of (0; o)) can differ for each Fourier mode, ie. different Fourier modes initialized at different weights due to the filter initialization.
79

---

The generated distribution has variance E = FÂF*, whereAux =o along Fourier modes.
(ao) 3π(
controls variance
(x) = V EAN + o² exp [ 2
(1-)x(0)) Ei (-2Nnro?) e e-2NorSua Ei (22 (o² + Exc
)
exp
(245)
80

---

G.4 Local patch linear convolutional network
Now, let's consider local filter w.
Spectral representation and spatial locality constraint The major difference that if we have a locality constraint in w, the corresponding constraint in y will be constrained in lower dimensional space
2
4
Wg
7
k=
F
2nd N
71
K=U
But if our wh = 0, 2 K, then the spectral moment just sums over
K
2mi
kl Wk
71
k=0
Specifie case: kernel 3 convolution Consider the case where the kernel size is 3, only W1, wo, w-1 # 0. Then the spectrum can only have DC and cosine and sine component
E
2nt $ WR
2
RE(-1.0.1)
&/2 + w 1e
wo + wie
2nl Basically, the small kernel size mandates that the spectral view of it y cannot vary very fast! Thus the spectral representation has be very smooth, just sine and cosine waves.
2Tl
wo [wi w-1) cos
i(w=1 w1) sin
G.4.1 Training dynamics of patch linear convolutional net
To study the training dynamics of the linear local convolutional network, it's easier to directly work with the entries of the filters.
Consider the Toeplitz weight matrix (general boundary condition).
W2
Wd
w
WO
Wi
W
W
Wo
W4-3
W d+1
W-d+2
w W-d+3
W0
Shift matrix formulation
Generally, Toeplitz matrix can be written as
d-1
W
wk
pk
0=
1=1
where the (non-circulant) shift matrix is P.
O
1
0
0
0
0
P
0
0
0
0
This compactly expresses the weight matrix as a weighted sum shift matrix powers, enabling us to write down the gradient efficiently.
81

---

Remark:
p For circulant shift matrix pz = P P 0 exactly.
For non-cireulant shift matrix, approximately pr p-1 but there will be some issues in the boundary condition.
Trace with also allows us to elegantly extract entries from matrix.
Multiplying this shift matrix allows us to shift the matrix
This decomposition of weights allows us to write down the loss and the gradient.
For k
ac Σ
Gú
Dwk
j-imk
dc
Σ
ÖWij ij
5-izk
(Pt, .)
-T[((PhTTG)
Tr((Pt)T Oc aw
for the other k.
Oc = ((PT)+. G)
Dw
Tr|P*
Tr[Pk dc ow
Using this formulation, express the gradient to each entry.
dc 2Tr 2Tr (P*)" E)]
Owk
-2TH((PH)]2) 2TH
Wm pm
w (pT)m) (02 23]
d-1
= =- + 2Tr
Σ
W'm pk pm
pm
W (PT (pt) (0 + Σ)]
m=0
d-1
d-1 Σ
-TH + Σ
Wm
Σ) +2
W- w-m T + E)]
===
1=1
Similarly, OC
Ow-k = =Tr pks + P
--2 + 2Tr Σj
--2Tr[phy] + 2Tr(P
Wm pm
w-m (p (021 ) + Σ)]
d-1
= - + 2Tr
W pk+m
pk (p' ( (0-1+
w
m=0
m=1
d-1
= -22Tr/pby] + 2 Σ
Wmn Tr (a²1 Tr Σ) + Σ w (p 2
m=0 0
82

---

We denote the general k-shift operator (non chon-circulant boundary condition)
S(k)) { Ph.
k: 2 0,
(p())
h < 0,
S(k)* (pb,²= = =(pT)A = S(-h)
Define special function with two indices
Tik, ml = To =))((()²
= =OTTT((())S(*)) Tr[S(m) E S(K)*)
The first term is diagonal, the second term is not. Non-diagonal terms all come from the second term. Note that this special function is symmetrie, which understandable, since it will be functioning as the new covariance matrix.
T Tm, k] = T=(5(m)6"1+
[1) -Tk,m) ml
Then the gradient can be written as this simplified equation, of similar structure as before.
dC 2Tr(S(R)*
+2
Σ
=))
Owk
m=1-d
d-1
--TT(((()) +2 Σ
Tkk m) Wm
m=l-d
-TT((((()2] +2
Σ Tik mlwm
(246)
Thus study the learning dynamics of convolutional linear models, the object that needs to be focused on is Tik m). the spectrum of it. Thus, the learning dynamics of w will be first along higher eigenvalues of T, and then lower eigen ones.
Interpretation of T matrix Note that the Tik m] matrix can regarded as the spatially averaged version of the covariance matrix, especially by averaging local cross covariances of pixels at distance k, m.
Optimal solution
The fixed point of the equation is the following linear equation
Σ
T k, mmwm Wm = Tr(S(k)2]
m
Let the vector Rs = Tr(S(k)E), k= {1 d,... - 1}
w² =
(247)
Here we successfully derived the matrix equation for gradient flow for weights w.
Cyclic case (circular convolution)
Now consider the eyclic shift matrix
0
1
0
0
0
0
P
0
0
1
1
0
0
83

---

then pT =P-1, so we have
T T - (o [)S(H)]
=Tr[P" (0²1+ E
111pm=+(021- Σ]]
= $ + Tr| Tmppm-
Note that
Tr(Pm-k] = N mm
So
Tik, m] N Tr(pp--k))
= =N 8mnk + Tr Tr[pm->E]
N
Note that the second term the shift-averaged version of the covariance matrix
1
Tr
N
N
km mod
N
Let's define a new averaged cross-covariance vector
x[k]
1
1
TT[P
N
Diitk mod
N
Note x[k] = x[-A]
R = Nx(=]
Tik m] = N N(02mmk 8mk + x(m - ki)
We want to solve the linear system
Tik mlwm Ru
Note that Tik m| a Toeplitz matrix.
In the circulant case the optimal solution satisfies,
K
(No2dem + Rs Rk-m)w", R
K
where R& is the averaged version of the covariance matrix. So again it's circulant regression with Ridge-like penalty.
Rs = Tr(Pk
let Σxm 1o2 Ra-m 222+1x2K+1 be the Toeplitz patch covariance matrix, then we can
write Rg as the center column of it
Σ
(No"dam + NEK. NÊ..O
m=-K
The solution can be written in vector form as
(No"I WoOI+NEE NÊJ=" NÊ.0
wº (No21 + N -I (021++8)-11
84

---

Thus the denoiser equals the central column from the Gaussian solution to the patch covariance. The gradient flow dynamics is
0C
d-1
2 +
Σ
8km
(248)
Dwk
m=1-d
We can write the flow dynamics in vector form
dw
OL
1I Ow
dr
2n (NE..0 Nio21 + Σ)W)
-nnNN (Σ:.0- (021 E)w)
Solution will be
W + eexp ( 2n)N E)) (w(0) w*)
(249)
with w = -0-2++)))11))
G.4.2 Sampling dynamics of patch linear convolutional net
Consider kernel of width 2K + 1,
Σ
2ri N Wk
n
kE K:K
K
2xi #t + w- ke 2 #
wi
+
Whe
k=1
K
2
2mTkl
w0
Wkk w- k) cos
i(w-k wx) sin
N
k=1
For symmetric filter weights wA = W-kk we have
K
2Tkkl
71 = wi +
2wh cos
(250)
N
During sampling have
1k(0)
do Ck
Ch. (x(
(251)
The key integral governing the variance
(f" (x()) - dx)
$ exp
(252)
2rrk exp wo(d) + Σ (w)(A) cos
N
dd)
Thus it integrates these sinusoidal modulations in the frequency domain.
G.5
Appendix: Useful math
Discrete Fourier Transformation (DFT) matrix
1
Fjs-
2xi jk/N
(253)
N
Property
85

---

Symmetry
F=
Conjugacy
FF =1
N
Fjk: Fim
e
2mi jk jk/N e22rkkm/N
N
k
2nd Mimm
jm
Properties of Circulant and DFT
Consider matrices of the following form M = FAF*
(254)
Explicitly, the matrix reads
M,k
Σ
Am mik
jm
1
mk
exp
exp
2
VN
N
VN
N
1
jm N
.ak
Am exp
2nT
2mi
N
N
1
Σ Am exp
m(j
k))
2Ti
N
N
Thus we see Myk depends only on j - k, hence is circulant. Let's define the circulant coefficients G-x MjAA-
Then we have for A 0, 1,
1 Σ Am exp
mA N
CA
2
(255)
N
The special case is the "DC" non-oscillating term, which are the diagonal values in M
CD
Σ
Am
N
the adjacent sub-diagonal values in M are
-2mmi m
C1 =
Am exp
N
N
-2 m(N 1)
CN-1
Am exp
N
N
1
(+2Ti) m
Am exp
N
N
More generally C = c are complex conjugates, or equal in the real case.
Properties of the A spectrum Since E is real symmetric, the eigenvalues exhibit mirror, i.e. even symmetry d Ax for k 1,2, ..., -
There are one or two standalone eigenvalues: when N odd, Ag is the zero-th frequency, DC component, which is unpaired.
When N even, An (DC component) and AN//2 (Nyquist frequency) are both standing alone, which are both unpaired.
86

---

H
Detailed derivation of Flow Matching model
Consider the objective of flow matching [28], at certain t
L = (xx.) x1rp1 [1((x;;1) - (xi - xn)
(256)
(x. x1 [[[([1 t)xo + txiit) -(x1 xx))l
(257)
x (1 - 1)x0 + lx1
(258)
Given linear function approximator of the velocity field,
u(x:t) t) = W + b:
(259)
[ = Exxr~((0.I), x1mm [Wo(1 - t)xo + txx) - bt - (x1 xo)ll
Exx,~((O.1), x1mp1 (Wa((1- - t)xo + (x1) + b, - (x1 xo)) " (Wr((1 - (t)o + +xxx + b, - (x1 - xo))
-Exo~N(0.I), x1mpi Tr ((1 t)xo + tx1) T W W:((1 t)xo + tx1) + -b + (x1 - xx)) Xo)
- -2(x1 - xx) - 2(x1 - xx)² Wr((1 - x)x + (x)) 2b ((1-(t))x)))) + +txi)
-Oxo=~(O.TT), Ximpi Tr
W W1(11 t)xo + tx ((1 - t)xo + tx1 )T + bb²
b + (x1 - xo) Xo)
- -2(x1 - xx) 3 bj 2w t)xo + (x1 xx0 + 2b W ((1-))x + tx1)
Similar to the diffusion case, it will also depend only on the mean and covariance of pr.
Exo.x [(1 t)xo [x1] ty
Exoxx [x1 xo] = /
Exp,x; [[x] x0)²"xx xo)] = =Ex,,x1 [x{x] x1 - 2x xo + x xol
= =Tr(= + u111 +1
Ex, x1 1((1 t)xo + (xx) xo) t(E u4i²) (1 - +t]I
Ex [((1 t)xu + tx1)((1 t)xo =x,)" 1² ( + ) +11+1+11111111)1) (1-t)² II
Taking full expectation, the average loss reads.
L= W. (P (10=)++p+))) + (1- + bbbb+ (Σ + +1p4 +1)
2p by 2w. ( (1²) (1 t)I) W
(260)
The gradients with respect to parameters are
VnL 2lb [WW m]
(261)
VwL 2 (f + ppr + +1 ()-)) (t mp²² (1-t)1) + tbby
(262)
Simplifying case = 0
Note the special case == 0
VbC 2b
(263)
VwC = 2 w (
(+Σ (1 -1
(264)
Optimal solution The optimal solution to the full case is b" μ p--t- tW* μ
(265)
W (IE (1 - f)11 + (1 (1-))21)
(266)
We can represent it on the eigenbasis of Σ, [u]]
tbp (1 -t)
W Σ
t2 + (1 - +) u
(267)
87

---

Asymptotics
Consider the limit - O,
Wi-0 1-0 = -
b =
Consider the limit ++ 1,
Wi-1 =1
b = - -- =0
H.1
Solution to the flow matching sampling ODE with optimal linear solution
Solving the sampling ODE of flow matching integrating from 0 to with the linear vector field
dx u(x:t)
(268)
dt
Under the linear solution case
dx Wi + bj
(269)
di
Simplified zero mean case 1
=
dx
tAg (1- 12 + (1 - 2)2
Un uuf x
dt
Solving the flow-matching sampling ODE mode by mode
u dx
Dd (1 1)
11
dt
td + (1 t)2
dcc(t)
the (1 11-t) t)
12 + (1 - 1)2 (t)
dl
(1) = = In| 12 + (1 - -1) + C.
(t) = CVP + +(--t)²
(270)
(1(1) = vpAx + t12
(0)
Cn (1() Ax (0)
This is the correct scaling xf x. The sampling trajectory of x; reads
x, Σ (t)u:
(271)
Σ Vetee + + u) xo
(272)
Thus, time the covariance of the sampled points is
E[xxx]]
Σ (t=xe de + u
(273)
with variance X = A2Ax + (1 -t)² along cigenmode uş
88

---

General case + 0
Wix + -W dx dt
th (1 t) t)
=M+ 13 (1 1)2 Uu u. (x tp)
(274)
can also be solved mode by mode. Redefine variable yt = Xt - 1p Then each mode can be solved accordingly, (t)) = y Yr
dy: dt
dx: 1
dt
Σ A2 + (1 - t)² un u
td -(1 -(1-t)
dcx(t)
tds (1
dt 12 + (1 ()2 (1)
Using the same solution above, we get the full solution for the sampling equation with any μ
X = 1u + Σ ex (()us
A
ty E
VP + (1 t)3 X
(275)
H.2 Learning dynamics of flow matching objective (single layer)
Simplifying case = 0, single layer network Note the special case μ = 0
VbL 2b
(276)
Vw -00---00000-00---000 (-1]
-
(277)
dW
dr WC
(278)
JW
dr 2n (111))=(11000000---11000
(279)
Using the eigenbasis projection
dWWu;
2n W + (1-t))1) -IE - (1 t]]] uk
(280)
dr
-2n Wu: (t Ax + (1-t) - (A) (1 - t))
(281)
= -2n > + +11 (1-1)3) (Wi: 12 (1 f)2 u1
Ld - -(1 - t)
(282)
W(T)uE
th - -(1 (1- - 1)
= Aexp (-207 (P>x d + (1 (1-)))))
(283)
ehx + (1 t)2
W(T)us = th t1
th - (1-t)
BA + (1 - 1) U + (W(O)u:= - 12x + (1 - 1)2 Us} exp (-217 (12 A + (1-
(284)
The full solutions of the weight and bias are
W(7) W + Σ ((O)):
th - (1 - t)
u T exp ( (- 1n227((2x A + +(1 (1--)))) (285)
e\ + t)2 U1
R
b(T) = -b0)= cx)()21p
(286)
is easy to see this solution has a similar structure that of the denoising score matching objective for diffusion models.
89

---

Remarks
The learning dynamics of weight eigenmodes at different times t were visualized in Fig. SA.
Note that the convergence speed of each mode is exp(-2n7 + (1- 1)))). Similarly, at the same time t, the higher the data variance Ag, the faster the convergence speed.
At smaller t, all eigenmodes converge at similar speed.
At larger t ~ 1, the eigenmodes resolved at distinct speeds depending the eigenvalues.
Note that for each eigenvalue λv. there is special time point where the convergence speed is maximized: t" = 1/(AA + 1).
H.2.1 Interaction of weight learning and flow sampling
Consider the sampling dynamics of flow matching model with learned weights:
dx W(T, t)x++)+++bl b(T,
(287)
Assume the weight initialization is aligned and the same across
w(O,t) Euuu
(288)
then
W(r,t)= = W + Σ (W()) - W uuuu exp
2t X1 + (1 (1-t1) (1-t)²))
(289)
thx (1 - 1)
u Σ (o
DA - (1-
un exp (-2 217 (1 A2 + (1- (11-19))
a2 + (1 2)2
a2 + - -t)2
k
(290)
Ignoring the bias part, consider the weight integration along C₁ (t) 1z,x(1)
d
1)= - (1 12 + (1- + (ox f2 + (1 - 1)2
DA (1 (1-1) exp ( - (t>x + (1 (1-)))) cx(t)
di (t)
(1-
(291)
We have the integration of the coefficient:
L' dt
LA= (1 -
(Qx:
Dh (1 11)
- -1177(12x (1 -1)1]]
P2 + 11-
+
fb (1 t)2
exp
20rdy Ac-1
(erf
V2
yr
+ erf
Aot1
Auti
4 4/11
1/T (Ax 1)
(Bi(-2nTT) - Bi (-21TTAR) log (Ax:)
1
AurAg Que
log (Ax) + (Bi(-2
Bi 21yTT
Ax+1
2
21/T (Ax + 1)
erf (V
2hpT
+ erf 1tp/
21TT
(292)
dv +1
do +
90

---

(1)
- exp((1)
exp S dt 1² (1 t)² + (Qx lu - (1 -
Lxx (1-1)
exp (-2n7 (te X + +1 ([--t)))]]
Aa + (1 1)2
(5 log (Ax) + 2 (Ei( (5(()117) Ei (- 211T 2x)) + 1
π
Que 2ur
cxp
2
2n (Ax + 1)
21/T
21/T
crf
erf A
A +
Xx +
= -eeyy A /
exp
(E((-21TT=- - Ei Ei(- (-2 (x)
exp
2
20yT (Au 1)
π
Qxe
Znr Ap+
erf
A +T+ + erf dey
21/T
201T
(293)
d +
A
(1)2
A
X
exp
(Ei( (Eis-2n5) - --ii-- ( xxt)
(294)
exp
2np (Au + 1) (V
Que Ax+1
2nrdz
erff
(v
2nf /2 erf
Xx +1+
(Sa// du + ))
2n7
Remarks
The learning dynamics of generated variance were visualized in Fig. 8B.
The power law relationship between the convergence time T of generated variance and the target variance A was shown in Fig. 9. For the harmonic mean criterion, the power law coefficient was also close to -1
H.3 Learning dynamics of flow matching objective (two layers)
Let PPT, Given the general loss,
(295)
= (Vw (VwL) = (Vw)T P
V,L T= (PE++1 (P (IE - (1 - -1)1]] P
- (IE (1 (1-1)
= -1-2(1--1111))
pp PP- + (1 13 p + (1 (1-+(111))))pp))))
91

---

Similarly, let let o/:
n 2 -2? (1 (1-
PP (1 + (1 T) + u (1-))))
=2 2(1) - - -1))) P + 0k Σ
pT Um Am2 + (1- (1-1)?)u??)
m
+ (P)x + (1- - (1-t)))ut t)² Σ pT Un u
= 2 2(11 - (1 - -1-1)) + ok Σ
4m Am + (1 (1-1))) + (f + (1 (1-1)))a Σ gn
m
V. C 1(t2; (1 f))gg= +2 " T9 9m) (f=) + (1 (1-t))) + =2(P>x
+(1 +1--t)9 > (K am dn
-44 (1 4)) +2
(t Am + &
+ t
+ (1 (1-1)
(1-) (g 9m))m
= -44(1>+ (1 () +2 Σ Am + +²xx + +[11 - 113) (a% qm) qm
(296)
m
Simplifying assumption: aligned initialization V C - - (1- g) + 4(2 + (1 - (t))))(( (
Assume at initialization of 9m 0, Vk: + m:
(297)
The learning dynamics follow:
dr Qx -1IV ax C
=4n (1As - (1 (111)) - -PPA + +1-- (1))))) (w)) 5x
1
(9% 9x) = 4n An (t>s - (1 (1-th) -)) - (t2xs Ao + (1 +1-t)) (n0) (ol 9b)
(298)
2 dr
Using abbreviations:
A = 8n((1x - (1 (1-11 f))
B = 871(1 A + (1 - ()?)
r(T) = |1 1
we can see the core structure:
dr(r) Ar Br
(299)
dr
With initialization r(0) 0% (T 0) Qk. the solution reads:
A
1
11qu || (T)
B + (Qu
1
A)e-Ar
t)x - -(1 -11-t)
Qk
e>x + -12
Qx + 1de (1-))
Qu)e SeT (1-))(
(Puu(1I-))
Q
=Q
(300)
Q2 +Qi Qx)e 8ny (11)-(1-1))
where
1Ax (1-1)
Q
(301)
12 + (1 1)2
92

---

So the overall weight dynamics read:
W(r;()= = Σ |lull (r)))(())uuu,
Qs
Σ
2i uu
Qx)e 807 (che (1-x-(1-0)
Qk (Q
Qx
typ - (1
(302)
Qx)e 88y (1-2-(1-0)
t² + (1- t)2
Q+ + (Q
Note that the optimal weight (301) is not positive definite, but the two-layer symmetric weight with denoiser loss, there are scenarios:
network is. So, different from the diffusion
When 1>
1 AN+I
A> 0 and Q" > 0. The solution converges to Q": lim 11 [1² =
Qi
When &
XEFI A < 0 and Q < 0. In this case, the optimal solution is "non-achievable"
by a two-layer symmetric network. lim,- e cc, so lim, || ( 0. In other words, 0 becomes stable fixed point instead of Q", and the solution ||gx || (T) will be attracted to and stuck at
Thus,
thp t1 12 + (1 -)- 1)2
lim W(r;t)= 10000
Σ
Us u
(303)
R, where A A
Because of this, asymptotically speaking, the symmetric network architecture pd P will not approxi- mate this vector field very well.
Thus, for the purpose of studying the learning dynamics of flow matching models, some extension beyond the symmetric two-layer linear network required for a thorough analysis.
93

---

I Detailed Experimental Procedure
I.1 Computational Resources
All experiments were conducted on research cluster. Model training was performed on single A100 H100 GPU. MLP training experiments took 20mins-2hrs while CNN based UNet training experiments took 5-8 hours, using around 20GB RAM.
Evaluations were also done single A100 H100 GPU, with heavy covariance computation done with CUDA and trajectory plotting and fitting on CPU. Covariance computation for generated samples generally took few minutes.
1.2 MLP architecture inspired by UNet
We used the following custom architecture inspired by UNet in [26] and [62] paper. The basic block the following
class (EEMMPPBoocktorrch...o :onoulleee. def __init__
in_features, out_features, meb_feeaures, dropout=8, skip_scale=1, eps=1e-5, adaptive_scalesTrue, init=dict(), init_zero=dict(), ):
super()...int__4 self.in in_features
self. huffoout_feature = out fea_ures
self.emb_feature emb_features emb_features
self. delf.dropout: dropout
self.skip_scale = skip_scale
self aelf.adaptive_scal. = adaptive_scale
self.normè mm.LayerNormiin__etture. eps=eps)
eps=eps)
self. fce = Linear ftt__ffaturesii_feetturee self.affine Linear (inferrurtetttesseen_eeee. if adaptive_scale else 1), *xinit)
features=out_featurre **init)
out uur_Feestrrssouutffaureesss
self.norml mm._ayerNorrmmuut_feetttte featuueu eps=eps)
eps=eps)
self.fc1 = Linear (in_fraauresssottfeetrre,,, out_featuressott_eeatrrs
**init_zero)
self.skip None
if out features != in_features:
selfskin Linear *winit)
out_.fatturesoott_eettrees
def forward(self, x, emb):
orig = self.fco(F...
if h.FFfroff.sllussaf.nomeee//
params = self (em)) to(x. .unsqueeze(1)
scale, shift params. chunk (chunkss, dim=1)
waw.l(ttoocc.addcul@seeee.... self.normi (x), scale 1))
else:
x = F.silu(self.non normi (parammm))))))))
fel @el..fcc(F.dropouu(x asself.dropout saadinggself..raininnn
x.add_(sell skip(orig) if self.skip None else orig)
x skip.skipsccall
return
94

---

and the full architecture backbone
class UNetBlockStyleMLP harpdyyyeeee___aabbbobbbeee Laass (nn. Modulee) "* time-dependent score-based model.
def ..initt__self,, ndim=2, nlayers=5, nhidden=64, ltr___eeeeedmmmeeee :
super (aaerrri...i____
self. embed =
scale=1)
layers = nn ModuleeeLsst layers
PBlock (ndim, nhidden, Sime_embed__dee _dmm)
for in range range (nlayers-2)-
ayvers.appenddUNettie P0 nhidden, (ime_embed_dim)) Lin. Linear (nhidden, ndim)) self.net = layers
def forward(self, x, t_enc, cond=None):
# t_enc ireconditioned version of sigma, usually # In_std_vec torch. log(ssd_vec) / if cond is not None:
raise NotImplementedError ("Conditional training not implemented t_embed self. wewf.embed embed(t_enc)
for layer in self.net[:-1]
x layer(x, t_embed) pred self.net(-1](x)
return pred
class EDMPrecondKrapper (odulle:
def __init___seell model, 1ggma_data=0.5, 1igma_min=0.cm 802, sigma_max=80,
rho=7.0):
super()..iiit__...
self.model model
self.sigma_dar sigma_data
self.sigma_min sigma_min
self.sigma_nd max sigma_max
self.rho = rho
def forward(self, X, sigma, cond=None, ):
si sigma[sigma == 0] self.sigma_mt
## edm preconditioning for input and output
## https://github
py3
unsqueze sigma have same dimension
(which may have dim)
sigma_vee vec signssiiaaddaat _attaaaaaaaaaa ** / (sigma_vec ** self.sigms selffsiima_daar data ** 2)
view[[I,, ] + [1, 1 (X.ndim 1))
c_skip = self.sigma c_out sigma_vec * self. sigma_aaaa data / (sigma_vec 2 self dataadata **
2). sq)
c_in 1 (self.sigma_dee _ata
sigma_vec sg.sqrt()
c_noise Aigma.log() /
model_out self w.t...ooelll..... (c_in X, c_noise, cond=cond)
return c_skip * X + c_out model_out
This architecture can efficiently learn point cloud distributions. More details about the architecture and training can found code supplementary.
1.3 EDM Loss Function
We employ the loss function LEDM introduced the Elucidated Diffusion Model (EDM) paper [26], which is specific weighting scheme for training diffusion models.
For each data point x E TRd, the loss computed follows. The noise level for each data point sampled from log-normal distribution with hyperparameters Paaza and Pad (c.g.. P Paaza -1.2 and
95

---

Inbal sample Traiiin samples
Figure 29: Example of learning to generate low-dimensional manifold with Song UNet-inspired MLP denoiser.
Pau 1.2). Specifically, the noise level is sampled via cxp (Pmean + Pau c); e~ 10.
The weighting function per noise scale is defined as:
a² data )
w(o)
(a ddma
with hyperparameter daaa (e.g., data 0.5). The noisy input y is created by the following.
y = x + on, n n~~ No,, La);
Let Do(y, labels) denote the output of the denoising network when given the noisy input y, the noise level a, and optional conditioning labels. The EDM loss per data point can be computed as:
[(x) = w(a) DD + on, a, labels) - x||
Taking expectation the data points and noise scales, the overall loss reads
LEOM = Exxpaane Bn E.r w( (Dx( on, a, labels) - x1]
(304)
class EDMLoss:
def sinit__selll P_mean=-1.2, P_std=1.2, sigma_data=0.9
self.P_meat P_mean
self.P_std P_std
self.sigma_date sigma_data
def (self, net, X, labels=None, ):
rnd_normal rorch.ra shapee( [1, ] * (X.ndim 1), devicesX.device)
# unsqueeze to match the ndim of sigma (rnd_normal * self.P_std self. P__mean)........ .epp))
weight = (sigma ** 2 + self. sigmadigmaddta ** 2) / (sigma * self sigmasignaddat, ** 2
# maybe augment
sigma D_yn = net(X n, sigma, cond=labels,
torch.randn_like(C
loss = weight * ((D_yn x) ** 2)
return loss
96

---

1.4
Experiment 1: Diffusion Learning of High-dimensional Gaussian Data
I.411
Data Generation and Covariance Specification
We consider learning a score-based generative model on synthetic data drawn from a high-dimensional Gaussian distribution of dimension d 128, 256, 512. Specifically, we first sample vector of variances
o² = (2...
where each o7 is drawn from torch. (tarct randn(. .))). We then them in descending order and normalize these variances to have mean equals to fix the overall scale. Denoting
log-normal distribution (implemented via
D diag(o)......
we generate a random rotation matrix R E Rdxd by performing a QR decomposition of matrix of i..d. Gaussian entries. This allows us to construct the covariance
Σ = R.DRR.
This rotation matrix R is the eigenbasis of the true covariance matrix. To obtain training samples {x:} C IR, we draw x, from // (0, Σ). In practice, we generate total of 10,000 samples and stack them as pnts. We compute the empirical covariance of the training set, Zemp Cov(pnts), and verify that it is close to the prescribed true covariance E.
14.2 Network Architecture and Training Setup
We train a multi-layer perceptron (MLP) approximate the noise conditional score function. The base network, implemented as
(ndim=d, nlayers=5, nhidden=256, time_embed__dim225
maps a data vector x E Fd and time embedding to a vector of the same dimension Rd. This backbone then wrapped in EDM-style preconditioner via:
model_precd EDMP (DDdereccdd data = 0.5, Fmin 0.002, maax 80, P 7.0).
which standardizes and scales the imput according the EDM framework [26].
We use EDM loss with hyperparameters P_mean- = 11.2, P. P__sdd = 1.2, and data = 0.5. We train the model for 5000 steps using mini-batches of size 1024. The Adam optimizer is used with learning rate 1r = 10-4 Each training step processes a batch of data from pnts, adds noise with randomized noise scales, and backpropagates through the EDM loss. The loss values each training steps are recorded.
L.433 Sampling and Trajectory Visualization
To visualize the sampling evolution, sample from the diffusion model using the Heun's 's 2nd order deterministic sampler, starting from z ~ Lz)
sassampler(modee z, num_steps = 20, Omin = 0.002, Fmax = 80, = 7).
We store these samples in sample_s letrllsssseeeeeeeeeeeeeeeeeeeeee to track how sampled distribution evolves over training.
1.4.4 Covariance Evaluation in the True Eigenbasis
To measure how well the trained model captures the true covariance structure, compute the sample covariance from the final generated samples, denoted sample- We then project Esaple and the true Σ into the eigenbasis of Σ. Specifically, letting R be the rotation used above, we compute
R
S sample R and R ER.
Since Σ R T is diagonal in that basis, then compare the diagonal elements of R. Sample R with dDa As training proceeds, we track the ratio diag(R' R ESample ())/diag(D) D) to observe convergence toward across the spectrum.
intermediate results, including loss values and sampled trajectories, are stored disk for later analysis.
97

---

1,5
Experiment 2: Diffusion Learning of MNIST MLLP
15.11 Data Preprocessing
For our second experiment, we apply the same EDM architecture to several natural image datasets: MNIST, CIFAR, AFHQ32, FFHQ32, FFHQS2-fixword, UFHQ32-randomword. All dataset except for MNIST are RGB images with 32 resolution, while MNIST BW images with 28 resolution. These images were flattened vectors, (784d for MNIST, 3072 for others) and stacked pnts matrix. We normalize these intensities from [0, I] to [-1, by x H x-0.5 The resulting data tensor pnts is then transferred to GPU memory for training, and we estimate its empirical covariance Eemp = Cov(pnts) for reference.
0,5
1.5.2 Network Architecture and Training Setup
Since the natural dataset is higher dimensional than the synthetic data in the previous experiment, use a deeper MI.P network: For MNIST:
model
UNeeBBlockstyleeM__baakeeee (atp:/eeee___aacaabceeennnnni 784, nlayers 8, nhidden = 1024, time_embed___ee dim 128).
For others
model UNetBlockStyleM.P hetrelolckktyyleee.p..bcoooonond 3072, nlayers 8, nhidden 3072, time_embed dim 128).
We again wrap this MLP in an EDM preconditioner:
model_precd EDMP EDMPreconddraaperr((ddell daata 0.5, min 0.002, m 80, 7.0).
The model is trained using the EDMLoss described the previous section, with parameters P P_mean- -1.2, P_std- 1.2, and data 0.5. We set the training hyperparameters to Ir 10-4, n_stpps 100000, and batch_size 2048.
1.5.3 Sampling and Analysis
As before, we define a callback function Gaussian noise ~ (O, I784) and applies the EDM sampler to produce generated samples. These intermediate samples are stored in sample_store for later analysis.
that periodically draws i.i.d.
In addition, we assess convergence the mean the generated samples by computing
||[[_out]] - Eipnts]1
and we track how this mean mrea-squared error evolves over training steps. We also examine the sample covariance Esample of the final outputs, comparing its diagonal a given eigenbasis to a target spectrum (e.g. the diagonal variances of the training data reference covariance).
All trajectories and intermediate statistics are saved disk for further inspection. In particular, plot the difference between Example and Σ an eigenbasis to illustrate whether the learned samples capture the underlying covariance structure of the training data.
I.6 Experiment 3: Diffusion learning of Image Datasets with EDM-style CNN UNet
We used model configuration similar to /ttps:///iihhe https://gittubb.cm//1aaaaa www..ooo.//leeeeeeeeeeee but with simplified training code more similar previous experiments.
For the MNIST dataset, trained UNet-based CNN (with four blocks, each containing one layer, no attention, and channel multipliers of 1, 2, 3, and 4) on MNIST for 50,000 steps using a batch size of 2,048. learning rate 10 16 base model channels, and an evaluation sample size 5,000.
For the CIFAR-10 10 dataset, we trained a UNet model (with three blocks, each containing one layer, wide channels of size 128, and attention at resolution 16) for 50,000 steps using batch size of 512, learning rate 10-4, and an evaluation sample size of 2,000 (evaluated in batches of 1,024) with 20 sampling steps.
For the AFHQ, FFHQ (32 pixels) dataset, we used the same UNet architecture and training setup, with four blocks, wide channels of size 128, and attention at resolution 8, trained for 50,000 steps
98

---

with batch size of 256 and a learning rate of x 10 10-4. Evaluation was conducted on 2,000 samples batches of 512.
For the AFHQ, FFHQ (64 pixels) dataset, we trained UNet model with four blocks (each containing one layer, wide channels of size 128, and attention resolution 8) for 250,000 steps using batch size of 256, learning rate of 1 x 10-A, and an evaluation sample size of 2,000 (evaluated in batches of 512).
1.7 Architectural Ablation: Diffusion Learning of Image Datasets with E.M-style CNN Res.Net
To systematically examine the effects of network depth and width diffusion learning dynamics, conducted a controlled set of experiments using EDM-style CNN ResNet architectures trained on the FFHQ-32x32 dataset.
We designed a simplified single-resolution ResNet denoiser "SongUNNe without skip connec- tions or attention. It consists of a stack of residual convolutional blocks conditioned on positional timestep embeddings, followed by normalization and output convolution. This minimal EDM-style architecture isolates the effects of depth and width diffusion learning dynamics.
All models followed the EDM training configuration [26] with simplified code consistent with previous experiments, and were trained for 50,000 steps using batch size of 256, Adam optimizer, learning rate of 10-4 2,000 samples are generated given training steps for evaluation. No attention layers were used.
We systematically varied residual layers per block (L E {1, 2, 3, 5)); and (2) network width, by varying the base channel dimension (C [4,6, 6,8, 8, 12, 16, 32, 128, 256}). Each configuration
architectural factors: (1) network depth, by adjusting the number of
trained independently with identical optimizer settings and no data augmentations to isolate the contribution of architecture to convergence and spectral scaling behavior.
99


Reviews:
[{"review_id": "5nvi2C0kO1", "reviewer": "Reviewer_d8oW", "summary": "This paper analyzes the spectral bias in denoising score matching. They analyze the learning dynamics of diffusion models, focusing on how the generated distribution evolves throughout training under some strong assumptions to ensure analytical traceability. Additionally, they also empirically test their insights extracted from the theoretical analysis in a more realistic setting. They provide both theoretical and empirical evidence for spectral bias i.e., weights converge faster along the eigenmodes or Fourier modes with high variance, and the learned distribution recovers the true variance first along the top eigenmodes.", "questions": "See weaknesses", "limitations": "Limitations are discussed in Section 7.", "rating": 5, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 3, "significance": 3, "originality": 3, "strengths_and_weaknesses": "# Strenghts\n1. This paper studies an interesting problem with a very sensible methodology. They first start by considering the simple toy setting with linear denoisers where everything is analytically tractable, and then empirically validate the results in a more realistic setting.\n2. The presentation is clear. And the paper is well structured.\n# Weaknesses\n1. It seems that the experimental results deviate more and more from the predicted theory as we move towards the more realistic settings\n -  Even before introduction of Convolutions in Experiment 2 with MLPs on natural images the observed power-law exponent is significantly attenuated(from the theoretical prediction of -1 to -0.48)\n -  In the most realistic experiment(experiment 3) the observed power-law exponent is approximately 0, meaning that the spectral bias is absent.\n\nIf in realistic settings the purported phenomenon(spectral bias) is no longer observable then I am not sure what is the upshot of the analysis in this paper. Is it just an idiosycracy of simplistic setting studied in section 4? Or is there reason to believe that it is a more general phenomena for all diffusion models? I think a more detailed discussion connecting the analysis connecting the analytically tractable setting and the realistic setting would be highly beneficial, especially since the empirical analysis is inconclusive.", "was_revised": true, "final_justification": "I concur with the other reviewers that this is a decent contribution so I suggest acceptance.", "rebuttal": "We appreciate the reviewer’s positive assessment of our work and the insightful comment regarding the attenuation of spectral bias in more realistic settings.\n\nWe 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.\n\n**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. \n\n**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.\n\nOur 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. \nTo 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.)\nSurprisingly, 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.\n\nThese 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.\n\nIn 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.\n\n----\n**Numerical simulation of the patch convolution learning dynamics**\nWe 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. \n\nDuring 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**. \n\nIf 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. \n\nIf 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. \n\nIn 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. \n\nTab.S1 Spectral convergence time computed for patch linear convolutional networks | FFHQ32\n\n|  | Increase scaling | Decease scaling |\n| --- | --- | --- |\n| P=3 | $0.06\\lambda^{0.48}$ | $11.88\\lambda^{-0.21}$ |\n| P=5 | $0.05\\lambda^{0.24}$ | $3.66\\lambda^{-0.36}$ |\n| P=7 | $0.05\\lambda^{0.01}$ | $3.24\\lambda^{-0.37}$ |\n| P=11 | $0.07 \\lambda^{-0.22}$ | $2.61 \\lambda^{-0.40}$ |\n| P=15 | $0.08\\lambda^{-0.28}$ | $2.57\\lambda^{-0.40}$ |\n| P=25 | $0.09 \\lambda^{-0.30}$ | $2.54 \\lambda^{-0.40}$ |\n\nWe 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."}, {"review_id": "E0mTzNSEqn", "reviewer": "Reviewer_gTYo", "summary": "The paper presents an analytical study of how diffusion models learn the data distribution over time during training. Focusing on linear denoisers trained with full-batch gradient flow, the authors derive closed-form expressions for both the model parameters and the generated distribution. They find that modes with higher variance are learned faster, leading to what they describe as a spectral bias. This inverse-variance effect is extended to deeper linear networks and linear convolutional architectures. The authors also run experiments with MLP and convolutional UNet models on synthetic and image datasets to compare against the theoretical predictions.", "questions": "**1.** Is there any difference between the EDM and DSM formulations, or are they equivalent in the context considered here?\n\n\n**2.** Do the authors consider feasible and useful to carry a similar analysis with a random features model parametrising the score? While still linear in the parameters, the model is non linear in the data.\n\n\n**3.** I believe there is a typo in Line 130: it should refer to Appendix C.4, not C.3.\n\n**4.** In the discussion, the authors mention local convolutional structure. Could they elaborate on how one might, in principle, build a theoretical framework to handle this?", "limitations": "Yes.", "rating": 5, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 4, "significance": 3, "originality": 3, "strengths_and_weaknesses": "**Quality**\n\nThe analytical results are clearly stated under well-specified assumptions. The paper links training dynamics to a mode-wise description of learning speed, which is insightful. The experimental validation is well done, spanning both synthetic and natural image data. \n\nA limitation is that the theory relies on idealized conditions (linear networks, full-batch gradient flow, and aligned initialization) which restricts its direct applicability to practical settings. The experiments with CNNs show deviations from the theory, which are acknowledged.\n\n\n**Clarity**\n\nThe paper is generally well-written and well-structured. Assumptions are stated clearly, and the figures effectively illustrate the main points. The appendix is very detailed, and it's clear that the authors have put significant effort into it.\n\n**Significance**\n\nUnderstanding how diffusion models learn the data distribution over time is an important problem, and the paper provides theoretical insights on this front. The proposed framework gives a nice picture of how different parts of the data are learned at different speeds, which could be relevant for understanding failure modes or improving training strategies.\n\nThat said, the practical impact is somewhat limited by the idealized assumptions behind the theory, though qualitative features of the analysis are noted in more practical settings.\n\n\n**Originality**\n\nThe work offers an interesting perspective on spectral bias in diffusion models, focusing on the training phase rather than the sampling process. The mode-wise view of convergence time and its link to inverse-variance scaling seems novel. Some elements of the setup, like the linear score approximation and Gaussian assumption, are standard, but the way they’re combined here leads to new observations about the learning dynamics in this context.", "was_revised": true, "final_justification": "I thank the authors for their careful and detailed rebuttal addressing all questions and remarks. I believe that simplified solvable settings are highly valuable for providing insights into complex systems. As shown by the authors, some trends of more complicated models are captured. Therefore, I am raising my evaluation.", "rebuttal": "We appreciate the reviewer’s positive assessment and detailed comment on our paper! Here are the response to the questions. \n\n**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.). \n\nThe 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. \n\nWe 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. \n\n**Q2**:  Random feature model is indeed an interesting and tractable next step, and we were planning to carry similar analysis for them. \n\nAs 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. \n\nIn 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. \n\nOne 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. \n\n**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. \n\n**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.  \n\n**Numerical simulation of the patch convolution learning dynamics**\nSpecifically, 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. \n\nDuring 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**. \n\nIf 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. \n\nIf 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. \n\nIn 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. \n\nTab.S1 Spectral convergence time computed for patch linear convolutional networks | FFHQ32\n\n|  | Increase scaling | Decease scaling |\n| --- | --- | --- |\n| P=3 | $0.06\\lambda^{0.48}$ | $11.88\\lambda^{-0.21}$ |\n| P=5 | $0.05\\lambda^{0.24}$ | $3.66\\lambda^{-0.36}$ |\n| P=7 | $0.05\\lambda^{0.01}$ | $3.24\\lambda^{-0.37}$ |\n| P=11 | $0.07 \\lambda^{-0.22}$ | $2.61 \\lambda^{-0.40}$ |\n| P=15 | $0.08\\lambda^{-0.28}$ | $2.57\\lambda^{-0.40}$ |\n| P=25 | $0.09 \\lambda^{-0.30}$ | $2.54 \\lambda^{-0.40}$ |\n\nWe 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."}, {"review_id": "yRG9uZs80z", "reviewer": "Reviewer_PL6U", "summary": "This paper presents a theoretical analysis of the learning dynamics in diffusion models, focusing on the simplified setting of linear denoisers trained on data equivalent to a single anisotropic Gaussian distribution. The authors derive a closed-form solution for the weight dynamics under gradient flow (Equation 7) and, assuming commutativity of the weight matrices across noise scales, provide an analytical expression for the generated data distribution at any point during training (Proposition 4.2).\n\nThe central implication of such is the existence of a \"spectral bias\", where the time required to learn an eigenmode is inversely proportional to its variance ($\\tau_k\\propto \\lambda_k^{-1}$). This results in a coarse-to-fine learning trajectory, with high-variance modes being learned much faster than low-variance ones. The analysis is extended to deeper models, such as a two-layer symmetric linear network (Proposition 5.1), and to linear convolutional networks, where a similar spectral bias is predicted in the Fourier domain (Propositions 5.3 and 5.4).\n\nTo validate these theoretical predictions, the authors conduct empirical studies using MLP and CNN-based denoisers. The results show that the learning dynamics of MLPs align well with the theory. In contrast, CNNs show a different behavior where many modes are learned near-simultaneously, suggesting that the bias shifts from the global image covariance to the patch-level covariance, a phenomenon the authors attribute to the inductive bias of local convolutions.", "questions": "Can you elaborate on Figure 17~20? While I agree they resemble the increase-decrease separation in Figure 2, the tailing eigenmodes seem not to show strong incremental learning, also the scatter points seem not to agree with the fitted line (both increase fit and decrease fit).", "limitations": "While I commend the authors for their effort in building an analytical theory, the paper offers limited practical insight due to its restrictive core assumptions. Contemporary work is beginning to offer precise insights into phenomena like the effect of early stopping (https://arxiv.org/abs/2505.17638), but this paper's impact is constrained by its idealized setting.", "rating": 4, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 2, "significance": 2, "originality": 2, "strengths_and_weaknesses": "Strengths:\n* The paper provides a thorough and rigorous gradient-flow analysis for the learning dynamics of linear denoisers, with extension to two-layer and convolutional cases.\n* The theoretical claims are supported by empirical validation on both synthetic and real-world image datasets (Figures 3 and 9).\n\n\nWeaknesses:\n* The primary limitation is the reliance on strong simplifying assumptions, a linear denoiser and data drawn from a single Gaussian distribution which may not fully capture the complexity of practical diffusion models.\n* The extensions to symmetric linear and convolutional networks feel less developed. For example, while the authors attribute the faster convergence of CNNs to the theoretical \"N-fold speed-up\" from weight sharing, this connection is not rigorously established. The empirical speed-up could be influenced by many confounding factors in a modern U-Net architecture and optimization.\n* The presentation could be streamlined to better emphasize the central theoretical contributions. For instance, the key result for weight dynamics (Equation 7) is presented as an intermediate step rather than a standalone proposition or theorem. The frequent use of subsections for brief interpretive points also fragments the narrative flow.", "was_revised": true, "final_justification": "While I still have reservations about the practical utility of linear denoisers, I recognize that the paper is of high quality regarding extensive experiment results and rigorous derivations. Therefore, I am raising my score to 4.", "rebuttal": "We thank the reviewer for their thorough summary and thoughtful feedback!\n\n**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.\n\n**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. \n\nAdditionally, 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. \n\n**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.\n\n**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.\n\n**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.\n\nIn 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.\n\n----\n**Numerical simulation of the patch convolution learning dynamics**\n\nWe 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. \n\nDuring 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**. \n\nIf 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. \n\nIf 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. \n\nIn 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. \n\nTab.S1 Spectral convergence time computed for patch linear convolutional networks | FFHQ32\n\n|  | Increase scaling | Decease scaling |\n| --- | --- | --- |\n| P=3 | $0.06\\lambda^{0.48}$ | $11.88\\lambda^{-0.21}$ |\n| P=5 | $0.05\\lambda^{0.24}$ | $3.66\\lambda^{-0.36}$ |\n| P=7 | $0.05\\lambda^{0.01}$ | $3.24\\lambda^{-0.37}$ |\n| P=11 | $0.07 \\lambda^{-0.22}$ | $2.61 \\lambda^{-0.40}$ |\n| P=15 | $0.08\\lambda^{-0.28}$ | $2.57\\lambda^{-0.40}$ |\n| P=25 | $0.09 \\lambda^{-0.30}$ | $2.54 \\lambda^{-0.40}$ |\n\nWe 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."}, {"review_id": "6eDx0D5bit", "reviewer": "Reviewer_d6G6", "summary": "This manuscript delivers a rigorous analysis of diffusion model training dynamics under a linear denoiser assumption. Assuming the data follow a Gaussian distribution, the optimal denoiser becomes a linear operator whose convergence rate is governed by the eigenvalue spectrum of the data covariance matrix. While the theoretical results hold exactly for linear denoisers, the authors demonstrate through experiments that nonlinear neural networks exhibit qualitatively similar convergence behaviour.", "questions": "1. In line 94, does $s_\\theta$ denote the score function and $D_\\theta$ the denoiser? Please clarify to avoid confusion.\n2. Equation (6) defines $W_\\sigma=\\Sigma(\\Sigma + \\sigma^2 I)^{-1}$, but Equation (3) suggests $W_\\sigma=(\\Sigma + \\sigma^2 I)^{-1}\\Sigma$. Which form is correct?\n3. In Section 5.2, Proposition 5.3 uses a Fourier basis, whereas Proposition 5.4 relies on principal components of local patches. What justifies this change given that the filter width is the only differing parameter?", "limitations": "yes", "rating": 5, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 4, "clarity": 3, "significance": 4, "originality": 4, "strengths_and_weaknesses": "**Strengths**\n\n- Provides a sound theoretical framework linking covariance eigenvalues to denoiser convergence speed.\n- Empirical validation supports the theoretical predictions and reveals analogous trends in nonlinear denoisers.\n- Section 6.2’s experiments highlighting distinct CNN behaviors enrich the manuscript.\n- This work will be important in helping researchers understanding the generalization power of diffusion models.\n\n**Weaknesses**\n\n- The analysis is limited to linear denoisers. However, the authors do acknowledge this limitation and I wouldn't consider it a major drawback, as the empirical results suggest that the theory is still relevant for nonlinear denoisers.", "was_revised": true, "final_justification": "The paper is already solid and the authors provide further explanations to my questions especially regarding Fourier basis vs local patches. I will hence maintain my score.", "rebuttal": "We thank the reviewer for their careful reading and insightful feedback. Below, we address each of the questions in turn:\n\n**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.\n\n**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$. \n\nThus 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. \n\n**Q3**: Thanks for noticing this differences! \n\nAt 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. \n\nBy 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*. \n\nWe will expand our discussion of this technical choice in the camera-ready version to make the motivation and technical benefits clearer.\n\nWe appreciate the reviewer’s suggestions and will incorporate these clarifications into the final manuscript."}]

Author rebuttal:

Meta-review: This paper develops an analytical framework for analyzing the evolution of the learned distribution during diffusion‐model training. Via Gaussian equivalence,  the authors obtain closed-form gradient-flow dynamics of weights for one- and two-layer linear and linear-convolutional denoisers, showing convergence along principal components (linear) and Fourier modes (convolutional).
It proves a pronounced spectral bias: the convergence time of each mode follows an inverse power law in that mode’s variance. The reviewers appreciate the sound theoretical framework, empirical validation mirroring the theory, and the clear exposition. The authors should incorporate the clarifications and numerical simulations provided during the discussion and update the final manuscript to reflect the thoughtful suggestions made by the reviewers.
Decision: Accept (spotlight)