Year: 2022

Abstract:
The recently discovered Neural Collapse (NC) phenomenon occurs pervasively in today’s deep net training paradigm of driving cross-entropy (CE) loss towards zero. During NC, last-layer features collapse to their class-means, both classifiers and class-means collapse to the same Simplex Equiangular Tight Frame, and classifier behavior collapses to the nearest-class-mean decision rule. Recent works demonstrated that deep nets trained with mean squared error (MSE) loss perform comparably to those trained with CE. As a preliminary, we empirically establish that NC emerges in such MSE-trained deep nets as well through experiments on three canonical networks and five benchmark datasets. We provide, in a Google Colab notebook, PyTorch code for reproducing MSE-NC and CE-NC: here. The analytically-tractable MSE loss offers more mathematical opportunities than the hard-to-analyze CE loss, inspiring us to leverage MSE loss towards the theoretical investigation of NC. We develop three main contributions: (I) We show a new decomposition of the MSE loss into (A) terms directly interpretable through the lens of NC and which assume the last-layer classifier is exactly the least-squares classifier; and (B) a term capturing the deviation from this least-squares classifier. (II) We exhibit experiments on canonical datasets and networks demonstrating that term-(B) is negligible during training. This motivates us to introduce a new theoretical construct: the central path, where the linear classifier stays MSE-optimal for feature activations throughout the dynamics. (III) By studying renormalized gradient flow along the central path, we derive exact dynamics that predict NC.

heading: 1 INTRODUCTION
text: Modern deep learning includes paradigmatic procedures that are commonly adopted, but not entirely understood. Some examples are multi-layered architectures, stochastic gradient descent, batch normalization, cross-entropy (CE) loss, and training past zero error towards zero loss. Analyzing the properties of these practices is an important research task.
In this work, we theoretically investigate behavior of last-layer features in classification deep nets. In particular, we consider the training of deep neural networks on datasets containing images from C different classes with N examples in each class. After passing the i-th example in the c-th class through all layers except the last-layer of the network, the network outputs some last-layer features hi,c ∈ RP . The last-layer of the network—which, for each class c, possesses a classifier wc ∈ RP and bias bc ∈ R—then predicts a label for the example using the rule argmaxc′ (⟨wc′ ,hi,c⟩+ bc′). The network’s performance is evaluated by calculating the error defined by
Error = Ave i,c 1{c ̸= argmax c′ (⟨wc′ ,hi,c⟩+ bc′)},
while the weights, biases, and other parameters of the network (that determine the behavior of the layers before the last layer) in the network are updated by minimizing the CE loss defined by
CE = −Ave i,c log exp{⟨wc,hi,c⟩+ bc}∑C
c′=1 exp{⟨wc′ ,hi,c⟩+ bc′} ,
∗Equal contribution. Listed alphabetically.
where Ave is the operator that averages over its subscript indices.
Prior works such as Zhang et al. (2016); Belkin et al. (2019) have shown that overparameterized classifiers (such as deep nets) can “memorize” their training set without harming performance on unseen test data. Moreover, works such as Soudry et al. (2018) have further shown that continuing to train networks past memorization can still lead to performance improvements1,2. Papyan, Han, and Donoho (2020) recently examined this setting, referring to the phase during which one trains past zero-error towards zero-CE-loss as the Terminal Phase of Training (TPT). During TPT, they exposed a phenomenon called Neural Collapse (NC).

heading: 1.1 NEURAL COLLAPSE
text: NC is defined relative to the feature global mean,
µG = Ave i,c hi,c,
the feature class-means,
µc = Ave i hi,c, c = 1, . . . , C,
the feature within-class covariance,
ΣW = Ave i,c
(hi,c − µc)(hi,c − µc)⊤, (1)
and the feature between-class covariance,
ΣB = Ave c
(µc − µG)(µc − µG)⊤.
It is characterized by the following four limiting behaviors where limits take place with increasing training epoch t: (NC1) Within-class variability collapse3:
Σ†BΣW → 0, where † denotes the Moore-Penrose pseudoinverse. (NC2) Convergence to Simplex ETF:
⟨µc − µG,µc′ − µG⟩ ∥µc − µG∥2∥µc′ − µG∥2
→ {
1, c = c′ −1 C−1 , c ̸= c ′
∥µc − µG∥2 − ∥µc′ − µG∥2 → 0 ∀c ̸= c′
(NC3) Convergence to self-duality:
wc ∥wc∥2 − µc − µG ∥µc − µG∥2 → 0
(NC4): Simplification to nearest class center:
argmax c′ ⟨wc′ ,h⟩+ bc′ → argmin c′ ∥h− µc′∥2
1A similar phenomenon had been previously observed in boosting (Bartlett et al., 1998). 2An alternative line of work on “early stopping” (Prechelt, 1998; Li et al., 2020; Rice et al., 2020) advocates for terminating the training process early and shows that it could be beneficial when training on noisy data or small datasets. Such datasets are out of the scope of this paper.
3Our characterization of (NC1) is more precise than that given by Papyan, Han, and Donoho (2020), which only states ΣW→0 for rhetorical simplicity. The convergence of the trace-norm of Σ†BΣW to zero is the actual quantity measured by Figure 6 of Papyan, Han, and Donoho (2020) when demonstrating (NC1). Σ†BΣW
The (NC2) property captures convergence to a simple geometric structure called an Equiangular Tight Frame (ETF). An ETF is a collection of vectors {vc}Cc=1 having equal lengths and equal, maximally separated pair-wise angles. In classification deep nets, last-layer features are of higher dimension than the number of classes i.e. P>C. In this setting, the maximal angles are given by
⟨vc,vc′⟩ ∥vc∥2∥vc′∥2 =
{ 1, for c = c′
− 1C−1 , for c ̸= c ′ ,
and the ETF is called a Simplex ETF4. Observe that as C increases, the ETF approaches a (partial) orthogonal matrix. Thus, when C is large, this translates to the intuitive notion that classifiers and class-means tend with training to (near) orthogonality.

heading: 1.2 DEEP NET CLASSIFICATION WITH MSE LOSS
text: While classification deep nets are typically trained with CE loss, Demirkaya et al. (2020) and Hui & Belkin (2020) recently reported deep nets trained with squared error (MSE) loss,
L(W , b,H) =1 2 Ave i,c ∥Whi,c + b− yi,c∥22 + λ 2 (∥W ∥2F + ∥b∥22) (2)
= 1
2CN ∥WH + b1⊤CN − Y ∥2F +
λ 2 (∥W ∥2F + ∥b∥22),
achieve comparable test performance as those trained with CE loss. Above, H ∈ RP×CN and Y ∈ RC×CN are matrices resulting from stacking5 together the feature vectors hi,c and one-hot vectors yi,c as their respective columns; W ∈ RC×P is the matrix resulting from the stacking of classifiers wc as rows; b ∈ RC is the vector resulting from concatenating the scalars {bc}Cc=1; and 1CN is the length-CN vector of ones. Table 1 in Appendix A shows supplementary measurements affirming the findings of Hui & Belkin (2020); Demirkaya et al. (2020), i.e. the measurements affirm that MSE-trained networks indeed achieve accuracies on testing data comparable to those of CE-trained networks (cf. the analogous Table 1 in Papyan, Han, and Donoho (2020)). The analytically-tractable MSE loss offers more mathematical opportunities than the hard-to-analyze CE loss and inspires this paper’s main theoretical contributions explicitly characterizing MSE-NC.

heading: 1.3 CONTRIBUTIONS
text: Our main contributions are as follows:
• We propose a new decomposition of the MSE loss: L = LNC1 + LNC2/3 + L⊥LS. The terms LNC1 and LNC2/3 possess interpretations attributable to NC phenomena and assume the classifier W is exactly the least-squares optimal classifier WLS (relative to the given H); and L⊥LS captures the deviation of W from WLS. (Section 2)
• We provide empirical measurements of our decomposition (Figure 2) on realistic datasetnetwork combinations showing that L⊥LS becomes negligible during training—leading us to define the central path where L⊥LS=zero. (Section 2)
• We reveal a key invariance property on the central path. The invariance motivates the examination of a representative set of features, X=Σ− 1 2
W H , that we call renormalized features. (Sections 3.1-3.2)
• We study the gradient flow of those renormalized features along the central path and derive exact, closed-form dynamics that imply NC. The dynamics are explicit in terms of the
captures the subtlety that the size of the within-class “noise” (captured by ΣW ) should be viewed relative to the size of the class-means (captured by ΣB).
4Traditional research on ETFs (cf. Strohmer & Heath Jr (2003); Fickus & Mixon (2015)) examines the P≤C setting where C ETF-vectors span their ambient RP space. However, in the P>C setting of classification deep nets, the vectors can not span RP . Following the precedent of Papyan, Han, and Donoho (2020), we interpret ETFs as equal-length and maximally-equiangular vectors that are not necessarily spanning.
5Assume the stacking is performed in i-then-c order: the first column is (i, c) = (1, 1), the second is (i, c) = (2, 1),..., the N -th is (i, c) = (N, 1), the (N + 1)-st is (i, c) = (1, 2), and so on. This matters for formalizing Equation 11 in Corollary 2.
singular value decomposition of the renormalized feature class-means at initialization. (Section 3.3)
Additionally, we complement this paper with new, extensive measurements on five benchmark datasets—in particular, the MNIST, FashionMNIST, CIFAR10, SVHN, and STL10 datasets6 (Deng et al., 2009; Krizhevsky & Hinton, 2009; LeCun et al., 2010; Xiao et al., 2017)—and three canonical deep nets—in particular, the VGG, ResNet, and DenseNet networks (He et al., 2016; Huang et al., 2017; Simonyan & Zisserman, 2014)—that verify the empirical reality of MSE-NC i.e. they show (NC1)-(NC4) indeed occur for networks trained with MSE loss. These experiments establish that theoretical modeling of MSE-NC is empirically well-motivated. They are lengthy—together spanning four pages with seven figures and a table—so we collectively defer them to Appendix A.

heading: 2 DECOMPOSITION OF MSE LOSS
text: Inspired by the community’s interest in the role of the MSE loss in deep net training (Demirkaya et al., 2020; Hui & Belkin, 2020; Mixon et al., 2020; Poggio & Liao, 2020a;b), we derive a new decomposition of the MSE loss that gives insights into the NC phenomenon. First, absorb the bias vector into the weight matrix—by defining the extended weight matrix W̃ = [W , b] ∈ RC×(P+1)
and the extended feature vector h̃i,c = [hi,c; 1] ∈ RP+1—so that Equation 2 can be rewritten as
L(W̃ , H̃) = 1 2 Ave i,c ∥W̃ h̃i,c − yi,c∥22 + λ 2 ∥W̃ ∥2F . (3)
Using { h̃i,c } , define the entities µ̃c, µ̃G, H̃ , and Σ̃W analogously to those in Section 1.1. We
further define the extended total covariance and extended class-means matrices, respectively, as
Σ̃T =Ave i,c
(h̃i,c − µ̃G)(h̃i,c − µ̃G)⊤ ∈ R(P+1)×(P+1)
M̃ = [µ̃1, . . . , µ̃C ] ∈ R(P+1)×C .
Next, we reformulate, with weight decay incorporated, a classic result of Webb & Lowe (1990):
Proposition 1 (Webb & Lowe (1990) with Weight Decay). For fixed extended features H̃ , the optimal classifier minimizing the MSE loss L(W̃ , H̃) is
W̃LS = 1
C M̃⊤(Σ̃T + λI) −1,
where I is the identity matrix. Note that W̃LS depends on H̃ only.
Note that we can interpret W̃LS as a function of H̃ , leading us to identify the following decomposition of L where one term depends on H̃ only.
Theorem 1. (Decomposition of MSE Loss; Proof in Appendix B) The MSE loss, L(W̃ , H̃), can be decomposed into two terms, L(W̃ , H̃) = LLS(H̃) + L⊥LS(W̃ , H̃), where
LLS(H̃) = 1
2 Ave i,c
∥W̃LSh̃i,c − yi,c∥22 + λ
2 ∥W̃LS∥2F ,
and L⊥LS(W̃ , H̃) = 1 2 tr { (W̃ − W̃LS) ( Σ̃T + µ̃Gµ̃ ⊤ G + λI ) (W̃ − W̃LS)⊤ } .
In the above, LLS(H̃) is independent of W̃ . Intuitively, it is the MSE-performance of the optimal classifiers W̃LS (rather than the “real classifiers” W̃ ) on input H̃ .
6CIFAR100 and ImageNet are omitted because Demirkaya et al. (2020); Hui & Belkin (2020) showed that they require an additional scaling-heuristic on top of the traditional MSE loss. Rigorous investigation into this scaling-heuristic would have demanded mathematical analysis outside the scope of this paper as well as experimentation that would have consumed expensive amounts of computational resources.
The component L⊥LS(W̃ , H̃) is non-negative7 and is zero only when W̃ = W̃LS. Therefore, L⊥LS(H̃) quantifies the distance of W̃ from W̃LS. In short, the least-squares component, LLS(H̃), captures the behavior of the network when the classifier possesses optimal, least squares behavior. Meanwhile, the deviation component, L⊥LS(W̃ , H̃), captures the divergence from that behavior.
We can further decompose LLS(H̃) into two terms, one capturing activation collapse (NC1) and the other capturing convergence to Simplex ETF of both features and classifiers ((NC2) and (NC3)). Theorem 2. (Decomposition of Least-Squares Component; Proof in Appendix C) The least-squares component, LLS(H̃), of the MSE decomposition in Theorem 1 can be further decomposed into LLS(H̃) = LNC1(H̃) + LNC2/3(H̃), where
LNC1(H̃) = 1 2 tr { W̃LS [ Σ̃W + λI ] W̃⊤LS } ,
LNC2/3(H̃) = 1
2C ∥W̃LSM̃ − I∥2F .
Inspection of these terms is revealing. First, observe that LNC2/3 is a function of the class-means and MSE-optimal classifiers. Minimizing LNC2/3 will push the (unextended) class-means and classifiers towards the same Simpex ETF matrix8 i.e. (NC2)-(NC3). Next, note that the within-class variation is independent of the means. Thus, despite the fact that classifiers are converging towards some (potentially large) ETF matrix, we can always reduce LNC1 by pushing ΣW towards zero, which corresponds to (NC1) .
Figure 2 measures the empirical values of the above decomposition terms on five canonical datasets and three prototypical networks. It shows that L⊥LS(W̃ , H̃) becomes negligible compared to LLS(H̃) when training canonical deep nets on benchmark datasets. In other words, L(W̃ , H̃)≈LLS(H̃) starting from an early epoch in training and persists into TPT.
This motivates us to formulate the central path: P = { (W̃LS(H̃), H̃) | H̃ ∈ R(P+1)×CN } , (4)
where the notation W̃LS(·) makes explicit the fact that W̃LS is a function of H̃ only. Intuitively, for a classifier-features pair to lie on the central path, i.e. (W̃ , H̃) ∈ P , means that the “real classifier” W̃ will exactly equal W̃LS(H̃)—the optimal classifier that would result from fixing H̃ and minimizing L w.r.t. just the classifier. Combined with Theorem 1, we see (W̃ , H̃) ∈ P if and only if L(W̃ , H̃)=LLS(H̃). Figure 2 shows that classifier-features pairs lie approximately on the central path during TPT, allowing us to shift focus from L to LLS.

heading: 3 EXACT CLOSED-FORM ANALYSIS ON CENTRAL PATH
text: Consider λ = 0 in this section. For the subsequent theory, we adopt the unconstrained features (Mixon et al., 2020) or layer-peeled (Fang et al., 2021) modeling perspective9. We consider continuous training dynamics with gradient flow where time-of-training is denoted by the variable10 t≥0. Within this model, we analyze the dynamics of the renormalized features X=Σ− 1 2
W H on the central path d dt X = −ΠTXX (∇XLLS(X)) , (5)
7The middle term is the sum of positive-semidefinite matrices—one of which is positive-definite. 8Appendices C.1-C.2 elaborates upon this interpretation of LNC2/3 in more detail. 9Both the unconstrained features model (Mixon et al., 2020) and the (1-)layer-peeled model (Fang et al., 2021) advocate for the mathematical analysis of deep nets through variants of gradient flow directly on the last-layer features and classifiers. The layer-peeled model is so-named because it proposes to analyze deep nets by "peeling away" one layer at a time. The unconstrained features model is so-name because features are not constrained to be the output of a deep net forward pass but are rather free variables that can be optimized directly: This captures the nearly unlimited flexibility afforded to feature engineering transformations in modern deep nets by the many nonlinear, overparameterized layers.
10Intuitively, one should interpret t=0 as the time-of-training when feature-classifier pairs effectively enter the central path.
where the operator ΠTXX projects the gradient onto the tangent space of the manifold 11, X , of all identity-covariance features. For Σ− 1 2
W H to be well-defined, we assume that ΣW remains full-rank 12
during training. We call Equation 5 the continually renormalized gradient flow and will motivate it in more detail in later subsections. We now state our main theorem:
Theorem 3. (NC Under Continually Renormalized Gradient Flow) Consider the continually renormalized gradient flow (Equation 5) in which the dynamics are restricted to the central path, and the features are renormalized to identity within-class covariance. Assume the features are initialized with zero global mean. Then, Neural Collapse emerges on the renormalized features with explicit dynamics of the flow given in Proposition 2 as well as Corollaries 1 and 2 later in this paper.
Three assumptions standout which we discuss below:
(A1) Restriction to central path: Figure 2 shows that—in practice, during TPT—L⊥LS(W̃ , H̃) is observed to be negligible compared to the dominant term, LLS(H̃). This indicates (W̃ , H̃) is close to the central path (Equation 4) where L(W̃ , H̃) = LLS(H̃). One the central path, one is effectively training with LLS as the loss. (A2) Zero global mean: Proximity to the central path implies that the last-layer biases become close to bLS. Having biases b = bLS is effectively performing classification with no bias term (i.e. using only a linear classifier W ) on zero-mean data. To see this, define the globally-centered features,
11X is a variant of the well-known generalized Stiefel manifold (Absil et al., 2009, Chapter 3.4). 12Given the definition in Equation 1, the positive-definiteness of ΣW is unsurprising in deep net contexts because CN ≫ P i.e. the number of examples is larger than the width of the last-layer. Note this does not contradict (NC1). (NC1) is characterized by ΣW approaching the zero-matrix with training (relative to ΣB). In practice (for example, in our supplementary Appendix A experiments), we observed that the zero-matrix is never actually achieved during the training dynamics and ΣW indeed remains full-rank.
globally-centered means, and total-covariance matrix in unextended coordinates:
ĎH = H − µG1⊤CN ; ĎM = [µ1 − µG, . . . ,µC − µG]; ΣT = 1
CN ĎHĎH⊤.
Proposition 1 then reduces to the following forms for the unextended weights and biases: W̃LS = [WLS, bLS] = [ C−1 ĎM⊤Σ−1T , C −11C − C−1 ĎM⊤Σ−1T µG ] . (6)
For an arbitrary activation-target pair (h,y), the prediction error obeys
WLSh+ bLS − y = C−1 ĎM⊤Σ−1T h+ C −11C − C−1 ĎM⊤Σ−1T µG − y
= WLS(h− µG)− (y − C−11C). The last line exhibits two terms in parentheses, both traceable to the action of the bias bLS. The term h−µG demonstrates that the bias induces a global-mean subtraction on h, while the term y−C−11 demonstrates that the bias also induces a global-mean subtraction on the one-hot targets.
(A3) Renormalization to identity covariance: Renormalization is ubiquitous in the empirical deep learning literature through works such batch normalization and its variants (Ioffe & Szegedy, 2015; Salimans & Kingma, 2016; Wu & He, 2018; Ulyanov et al., 2016; Ba et al., 2016; Krizhevsky et al., 2012). It is also ideologically precedented in the theoretical ML literature through works such as Douglas et al. (1998); Banburski et al. (2019); Poggio & Liao (2020a;b). Within this paper, the covariance-renormalization is inspired by an invariance property of the loss as well as the intuitive concept of the signal-to-noise ratio from mathematical statistics (discussed later in Sections 3.1-3.2).
In the subsequent sections, we motivate and introduce the theoretical constructs necessary for proving Theorem 3 and for presenting the explicit dynamics of the flow.

heading: 3.1 INVARIANCE PROPERTY
text: By Assumption (A2), classification on the central path is equivalent to classification on globallycentered features using only the linear classifier WLS. Then, using Equation 6, we can equivalently represent the central path (Equation 4) as follows: Definition 1 (Zero Global Mean Central Path).
sP = {( WLS, ĎH ) ∣∣∣WLS = C−1 ĎM⊤Σ−1T } .
An invariance property holds on the central path. Let A denote a symmetric, full-rank matrix. Then, the explicit form in Equation 6 for WLS implies
WLS(AĎH)AĎH = C −1 (A ĎM)⊤ [(AĎH) (AĎH)⊤]−1 AĎH
= C−1 ĎM⊤A [ A−1 ( ĎHĎH⊤ )−1 A−1 ] AĎH
= C−1 ĎM⊤ ( ĎHĎH⊤ )−1
ĎH = C−1 ĎM⊤Σ−1T ĎH = WLS
( ĎH ) ĎH,
(7)
where the notation WLS(·) makes explicit that WLS is a function of a set of input features. Thus, the actual predictions made by the least squares classifier are invariant to choice of the coordinates in which we express ĎH , i.e. all features AĎH are equivalently performing. Among those coordinate systems, we will prefer the one in which the “noise” is “whitened” or “sphered.” Recall that we assume ΣW (ĎH)—where the notation ΣW (·) expresses the within-class covariance as a function of the features—is positive-definite. Consider the coordinate transformation ĎH 7→AĎH with A = Σ− 1 2
W ( ĎH).
In these coordinates, the features are “renormalized” to spherical covariance,
ΣW (AĎH) = AΣW (ĎH)A = I,
so we will call AĎH the renormalized features. The class-means of AĎH are
ĎM(AĎH) = A ĎM(ĎH) = Σ − 12 W ( ĎH) ĎM(ĎH),
where the notation ĎM(·) expresses that the means are a function of features. Thus, Equation 5 describes continuous training dynamics where we continually “renormalize” features to their equivalently-performing representer.

heading: 3.2 SIGNAL-TO-NOISE RATIO MATRIX ON THE CENTRAL PATH
text: The Σ− 1 2
W ĎM term evokes the canonical idea of a signal-to-noise ratio from statistics. To see this, note
that the classifier in Equation 6 can be rewritten as
WLS = C −1 ĎM⊤ ( C−1 ĎM ĎM⊤ +ΣW )−1 .
The combined presence of the terms ĎM ĎM⊤ and ΣW in the denominator, along with ĎM in the numerator, signals to us the presence of the informative signal-to-noise (SNR) ratio matrix:
SNR ≡ Σ− 1 2
W ĎM . (8)
Why this terminology? The globally-centered class-means represent the overall “signal” indicating that the class-means are separated from each other: if they were non-separated—e.g. all equal— ĎM would be zero. However, a classifier decision must keep in view noise—captured by ΣW —which might confuse the classifier. So, one wishes the signal be large compared to the noise. It may help to consider the mental model—which is not at all necessary to the correctness of our analysis—under which the activations are normally distributed with hi,c ∼ N (µc,ΣW ). Under that model, a linear classifier will indeed get confused if the norm ∥µc − µG∥2 is “small compared” to ΣW . Replacing this somewhat vague statement by a discussion of quantitative properties of the SNR matrix is well understood by mathematical statisticians to be decisive for understanding classification performance in the normal case. Additionally, observe that SNR = ĎM(Σ− 1 2
W ĎH). Thus, connecting back to
Section 3.1, the SNR matrix is simply the class-means matrix of the renormalized features.
The SNR can be further understood through its singular value decomposition (SVD). Definition 2 (SVD of SNR Matrix). Denote the SVD of the SNR matrix (Equation 8) as follows:
SNR = UΩV ⊤ = C−1∑ j=1 ωjujv ⊤ j .
Here, {ωc}C−1c=1 are the non-zero singular values of SNR; Ω = diag ( {ωc}C−1c=1 , 0 ) ∈ RC×C ; and the left and right singular-vectors, U ∈ RP×C and V ∈ RC×C , are partial orthogonal and orthogonal matrices, respectively, with columns {uj}Cj=1 and {vj}Cj=1. The SNR matrix is rank C − 1 since ΣW is assumed full-rank and ĎM is rank C − 1 due to global-mean subtraction.
The non-zero singular values, {ωj}C−1j=1 , are decisive for understanding the separation performance of the least-squares linear classifier. Good performance demands that the singular values be large. Works such as Fukunaga (1972, Chapter 10) and Hastie & Tibshirani (1996) show that the magnitude of each singular value of the SNR matrix is the size of the class separation in the direction of the corresponding singular vector. In this sense, the smallest singular value captures the smallest separation between classes. Consequently, driving-larger the smallest non-zero singular value makes the task of linear separation more immune to (spherical) noise and to Euclidean-norm-constrained adversarial noise. As we will see in the next section, examining the dynamics of these singular-values when the features undergo gradient flow in renormalized coordinates leads to precise characterizations of Neural Collapse13.

heading: 3.3 DYNAMICS OF SNR
text: For any matrix Z, we use Zt in this section to denote the state of that matrix at time t. Similarly, we denote with {ωj(t)}C−1j=1 the state of the non-zero SNR singular values at t. We now present Proposition 2 and Corollaries 1-2 that provide the explicit dynamics referenced by Theorem 3. Proposition 2 (Dynamics of Singular Values of SNR Matrix; Proof in Appendix D.4). Continually renormalized gradient flow on the central path (Equation 5) induces the following closed-form dynamics on the SNR singular values (Definition 2):
c1 log(ωj(t)) + c2ω 2 j (t) + c3ω 4 j (t) = aj + t, t ≥ 0, for all j = 1, . . . , C − 1. (9)
c1, c2, and c3 are positive constants independent of j, and aj is a constant depending on ωj(0). 13For additional, supplementary intuitions about invariance, renormalization, and continually renormalized gradient flow, see the exposition in Appendices D.1-D.3.
Corollary 1 (Properties of SNR Singular Values; Proof in Appendix D.5). SNR singular values (Definition 2) following the Equation 9 dynamics satisfy the following limiting behaviors:
1. limt→∞ ωj(t) = ∞ and limt→∞ ωj(t)4√t/c3 = 1, for all j = 1, . . . , C − 1.
2. limt→∞ maxj ωj(t) minj ωj(t) = 1.
The first fact shows that all non-zero singular values of the SNR matrix, Σ− 1 2
W ĎM , grow to infinity at an asymptotic rate of 4 √ t/c3. Intuitively, this shows that the “signal” becomes infinitely large compared to the “noise”, i.e. (NC1). The second fact shows that the singular values of Σ
− 12 W
ĎM (which has zero-mean columns) approach equality—implying that the limiting matrix is a Simplex ETF (Lemma 7, Appendix D.6) i.e. (NC2). Finally, Papyan, Han, and Donoho (2020, Theorem 1) show that (NC1-2) imply (NC3-4) on the central path. Corollary 2 formalizes these intuitions.
Corollary 2 (Neural Collapse Under MSE Loss; Proof in Appendix D.6). Under continually renormalized gradient flow (Equation 5), the SNR matrix (Equation 8) converges to
lim t→∞
1
ωmax(t) SNRt = Û0V̂
⊤ 0 , (10)
where Û0 ∈ RP×(C−1) and V̂0 ∈ RC×(C−1) are the left and right singular vectors of the SNR matrix (Definition 2) at t=0 corresponding to the non-zero singular values; and ωmax(t) is the largest singular value at time t. Furthermore, Corollary 1 implies the occurrence of (NC1)-(NC4) i.e. renormalized gradient flow on the central path leads to Neural Collapse.
Moreover, denoting the Kronecker product with ⊗, the renormalized features matrix converges to
lim t→∞
1
ωmax(t) Σ
− 12 W,t ĎHt = (Û0V̂ ⊤ 0 )⊗ 1⊤N . (11)

heading: 4 RELATED WORKS
text: Since the publication of Papyan, Han, and Donoho (2020), several works (Mixon et al., 2020; Lu & Steinerberger, 2020; E & Wojtowytsch, 2020; Poggio & Liao, 2020a;b; Fang et al., 2021; Ergen & Pilanci, 2020; Zhu et al., 2021) proposed mathematical frameworks ratifying Neural Collapse.
Among them, Mixon et al. (2020) and Poggio & Liao (2020a;b) also examine MSE-NC. Our contributions are distinct from the MSE-NC analyses of Mixon et al. (2020) and Poggio & Liao (2020a;b): The work of Mixon et al. (2020) relies on a linearization of the unconstrained features model ODE while Poggio & Liao (2020a;b)—in a special section co-authored with Andrzej Banburski—examine homogeneous, weight-normalized networks. In contrast, this paper considers the dynamics of the renormalized gradient flow on the central path (motivated by the discussions in Sections 2-3). Neither Mixon et al. (2020) nor Poggio & Liao (2020a;b) provide exact, closed-form dynamics as we do.
Outside the MSE setting, Lu & Steinerberger (2020); E & Wojtowytsch (2020); Ergen & Pilanci (2020); Fang et al. (2021); Zhu et al. (2021) examine the emergence of the NC properties under variants of the unconstrained features/layer-peeled model trained with CE loss. These works focus on characterizing the loss landscape or the global minima without describing the dynamics.
We provide a detailed survey of all above-mentioned papers in Appendix E.

heading: 5 CONCLUSION
text: In this paper, after verifying that NC occurs when prototypical classification deep nets are trained with MSE loss on canonical datasets, we then derive and measure a novel decomposition of the MSE loss. We observed that the last-layer classifier tends to the least-squares classifier—motivating us to define the central path on which the classifier behaves exactly as optimal least-squares classifier. On the central path, we showed invariance properties that inspired us to examine the renormalized features and their corresponding continually renormalized gradient flow. This flow induces closed-form dynamics that imply the occurrence of Neural Collapse.

heading: REPRODUCIBILITY STATEMENT
text: This paper is reproducible. Experimental details about all empirical results described in this paper are provided in Appendix A. Additionally, we provide PyTorch (Paszke et al., 2019) code for reproducing NC under both MSE and CE—as well as the Figure 2 MSE decomposition—in the following Google Colaboratory notebook: here. Experimental measurements used to generated the plots in this paper have been deposited in the Stanford Digital Repository, here. All datasets and networks used in our experiments originate from the PyTorch Model Zoo with pre-processing details described in Appendix A. Formal statements and proofs of all our theoretical results are provided in Appendices B-D.

heading: ACKNOWLEDGEMENTS
text: Some of the computing for this project was performed on the Sherlock cluster. We would like to thank Stanford University and the Stanford Research Computing Center for providing computational resources and support that contributed to these research results.
We acknowledge the support of the Natural Sciences and Engineering Research Council of Canada (NSERC), [funding reference number 512236 and 512265]. This research was funded by a Connaught New Researcher Award at the University of Toronto. This research was enabled in part by support provided by Compute Ontario (www.computeontario.ca) and Compute Canada (www.computecanada.ca). This work was also partially supported by NSF Division of Mathematical Sciences Grants 1407813, 1418362, and 1811614 and by private donors.

heading: APPENDIX

heading: A MSE Neural Collapse experiments 14
text: A.1 Datasets . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
A.2 Networks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
A.3 Optimization methodology . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
A.4 Computational resources . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
A.5 Formatting . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
A.6 Experimental results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
A.7 Exceptions: STL10-ResNet and STL10-DenseNet . . . . . . . . . . . . . . . . . . 21
A.8 Open Questions: Weight-decay, batch normalization, SGD, generalization, test data 21

heading: B Theorem 1 22

heading: C Theorem 2 22
text: C.1 Intuitions for Theorem 2 in unextended coordinates . . . . . . . . . . . . . . . . . 23
C.2 Additional intuitions for Theorem 2 terms . . . . . . . . . . . . . . . . . . . . . . 23

heading: D Theorem 3 23
text: D.1 Implications of invariance . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
D.2 Aligned SNR coordinates . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
D.3 Continually renormalized gradient flow . . . . . . . . . . . . . . . . . . . . . . . 27
D.4 Proof of Proposition 2 (gradient flow in aligned SNR coordinates) . . . . . . . . . 28
D.5 Proof of Corollary 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
D.6 Proof of Corollary 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31

heading: E Related works examining Neural Collapse 33
text: E.1 Mixon, Parshall, and Pi (2020) . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34
E.2 Lu and Steinerberger (2020) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34
E.3 E and Wojtowytsch (2020) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34
E.4 Poggio & Liao (2020a;b) (with Banburski) . . . . . . . . . . . . . . . . . . . . . . 35
E.5 Ergen & Pilanci (2020) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35
E.6 Fang, He, Long, and Su (2021) . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35
E.7 Zhu, Ding, Zhou, Li, You, Sulam, and Qu (2021) . . . . . . . . . . . . . . . . . . 36

heading: A MSE NEURAL COLLAPSE EXPERIMENTS
text: The experiments in this section examine the properties of Neural Collapse (NC) on deep nets trained using MSE loss. The direct MSE-analogues to the cross-entropy (CE) loss table and figures in Papyan, Han, and Donoho (2020) are in Table 1 and Figures 3-9 here. Furthermore, Figures 10-11 compares the MSE-NC experiments observed in this paper with the CE-NC behaviors in Papyan, Han, and Donoho (2020). Experimental descriptions and discussions are within the captions. Subsections A.1-A.5 describe formatting and experimental details. Then, Subsection A.6 collectively presents all experimental figures for this section. Subsection A.7 addresses outlier behaviors observed in the plots. Finally, Subsection A.8 discusses open questions inspired by NC.

heading: A.1 DATASETS
text: We consider the MNIST, FashionMNIST, CIFAR10, SVHN, and STL10 datasets with pre-processing the same as in Papyan, Han, and Donoho (2020).

heading: A.1.1 CHOICE OF DATASETS
text: Demirkaya et al. (2020) and Hui & Belkin (2020) showed multi-class classification deep nets trained with MSE loss—sometimes with a heuristical-scaling—demonstrate comparable test-error performance to those trained with CE loss. Even without additional scaling, i.e. with vanilla MSE loss, we found this true for ten class classification datasets—such as MNIST, FashionMNIST, SVHN, CIFAR10, and STL10—so we focus on these five datasets. For more than ten classes (i.e. datasets such as CIFAR100 and ImageNet), both Demirkaya et al. (2020) and Hui & Belkin (2020) showed that the additional scaling-heuristic does need to be applied to the loss before comparable test-performance can be achieved. In informal, exploratory experiments not reported here, we were able to reproduce their results on datasets with more than ten classes. But, we feel these scaling-heuristics merit further scientific investigation of their own and are beyond the scope of this article14, so we do not include datasets that require scaling-heuristic tuning.

heading: A.2 NETWORKS
text: We train VGG, ResNet, and DenseNet models with the same depth-selection procedures and architecture-specification choices as in Papyan, Han, and Donoho (2020). In the MSE setting, the final chosen depths were as follows:
DATASET VGG RESNET DENSENET
MNIST VGG11 RESNET18 DENSENET40 FASHIONMNIST VGG11 RESNET18 DENSENET250 SVHN VGG13 RESNET34 DENSENET40 CIFAR10 VGG11 RESNET50 DENSENET100 STL10 VGG11 RESNET18 DENSENET201

heading: A.3 OPTIMIZATION METHODOLOGY
text: The optimization algorithm, parameters, and hyperparameter tuning are the same as in Papyan, Han, and Donoho (2020).

heading: A.4 COMPUTATIONAL RESOURCES
text: The experiments were run on Stanford University’s Sherlock computing cluster. Each dataset-network combination was trained on a single GPU attached to a CPU with at least 32GB of RAM—the specific types of the CPUs/GPUs vary according to whichever was first assigned to us by the HPC cluster scheduler.
14Thorough experimentation with this scaling-heuristic would also consume prohibitively expensive amounts of computational resources.

heading: A.5 FORMATTING
text: The coloring and formatting of the plots are the same as in Papyan, Han, and Donoho (2020).

heading: DATASET NET ACC. 0-ERROR ACC. FINAL

heading: A.6 EXPERIMENTAL RESULTS

heading: A.7 EXCEPTIONS: STL10-RESNET AND STL10-DENSENET
text: The STL10-ResNet and STL10-DenseNet dataset-network pairs stand out among the experiments described in this section—sometimes displaying outlier behavior either by converging more slowly to Neural Collapse or by exhibiting trends inconsistent with those found in other dataset-network pairs.
STL10 stands apart from the other canonical datasets for multiple reasons15 and is a “less-typical” benchmark compared to more “compulsory” datasets like CIFAR10. In this particular case, we hypothesize that its outlier behavior might be caused by the size of the STL10 images (96 × 96 compared to the 32 × 32, after padding, of the other datasets)—leading to a higher-dimensional problem, which, in turn, induces a harder non-convex optimization problem, which might make SGD less likely to converge to a useful optimum or might make the convergence much slower and harder to observe under a fixed computational budget.
It was previously noted in Hui & Belkin (2020) and Demirkaya et al. (2020) that more challenging classification problems—in those projects, problems with more classes to be labeled—may require modifications to the MSE loss. We decided that such modifications are unnecessary in the ten class problems examined in this paper’s experiments. But perhaps the outlier nature of the STL10-ResNet and STL10-DenseNet combinations signal that the MSE loss modifications proposed by Hui & Belkin (2020) and Demirkaya et al. (2020) ought to have been used.

heading: A.8 OPEN QUESTIONS: WEIGHT-DECAY, BATCH NORMALIZATION, SGD, GENERALIZATION, TEST DATA
text: The experiments in this paper as well as those in Papyan, Han, and Donoho (2020) were conducted under the canonical setting where deep nets were trained using SGD with weight-decay and batch normalization—leading one to wonder about the roles that these particular ingredients play on the NC phenomena. Additionally, most observables focus on the train data—raising the question of how the NC phenomena behave on test data as well as their relationship to generalization. We find these questions intriguing, but feel each deserves careful experimentation outside the scope of this paper.
Nonetheless, we note here that several existing papers have already begun insightful investigations in these directions. For example, weight-decay, weight-normalization (as a proxy for batchnormalization), and SGD play key roles in the analyses of Banburski et al. (2019); Poggio & Liao (2020a;b) that, along other things, lead to the prediction of NC in homogeneous deep nets. Banburski et al. (2021) also explores the connection between NC and margin distributions with generalization.
Another example is Zhu et al. (2021) in which the authors conduct several experiments on ResNets trained on CIFAR10. On their models, the authors examine NC-related properties relative to the train data, test data, and randomly labeled data; they also conducted ablation studies varying control parameters (weight-decay, width, etc.) and the optimization algorithm (SGD, ADAM, L-BFGS). Based on their results, Zhu et al. (2021) propose thought-provoking conjectures on the role of each of these components (see Section E.7 for a brief survey of their findings).
As a preliminary exploration, we include Figure 12 showing the variability collapse (NC1) behavior on test data for the networks trained with MSE loss (used in this paper) as well as those trained with CE loss (released16 by Papyan, Han, and Donoho (2020)). From Figure 12, we see that variability collapse occurs much slower on the test data than on the train data. Although not shown here, the other NC phenomena behave similarly.
In this Appendix, we show Figure 12 since it is concise, within reach, and may be of interest to readers who share our curiosity. We refrain from presenting the entire series of multi-figure NC measurements on test data because it is outside the scope of this paper and would significantly lengthen this appendix. We also refrain from speculating on any underlying explanations until more meticulous experiments are conducted.
15For example, one sense in which STL-10 is clearly different from the other datasets is the comparatively small number of labeled training examples per class.
16The CE test data results were released by Papyan, Han, and Donoho (2020) on the paper’s corresponding Stanford Data Repository entry.

heading: B THEOREM 1
text: Theorem 1. (Decomposition of MSE Loss) The MSE loss, L(W̃ , H̃), can be decomposed into two terms, L(W̃ , H̃) = LLS(H̃) + L⊥LS(W̃ , H̃), where
LLS(H̃) = 1
2 Ave i,c
∥W̃LSh̃i,c − yi,c∥22 + λ
2 ∥W̃LS∥2F ,
and L⊥LS(W̃ , H̃) = 1 2 tr { (W̃ − W̃LS) ( Σ̃T + µ̃Gµ̃ ⊤ G + λI ) (W̃ − W̃LS)⊤ } .
Proof. Recall the MSE loss:
L(W̃ , H̃) = 1 2 Ave i,c ∥W̃ h̃i,c − yi,c∥22 + λ 2 ∥W̃ ∥2F . (12)
Consider a specific extended activations matrix H̃; the least-squares solution W̃ minimizing L(W̃ , H̃) with H̃ kept fixed must obey the following first-order optimality condition:
Ave i,c
(W̃LSh̃i,c − yi,c)h̃⊤i,c + λW̃LS = 0. (13)
From Proposition 1 in the main manuscript, the solution to the above is given by
W̃LS = C −1M̃⊤(Σ̃T + λI) −1.
The loss can be rewritten as 1
2 Ave i,c
∥W̃LSh̃i,c − yi,c + (W̃ − W̃LS)h̃i,c∥22 + λ
2 ∥W̃ ∥2F .
Combining the above with Equation 13 gives, after rearranging: 1
2 Ave i,c
∥W̃LSh̃i,c − yi,c∥22 + 1
2 Ave i,c
∥(W̃ − W̃LS)h̃i,c∥22 − λ tr { W̃LS(W̃ − W̃LS)⊤ } + λ
2 ∥W̃ ∥2F ,
which is equivalent to 1
2 Ave i,c
∥W̃LSh̃i,c − yi,c∥22 + 1
2 Ave i,c
∥(W̃ − W̃LS)h̃i,c∥22 + λ
2 ∥W̃LS∥2F +
λ 2 ∥W̃LS − W̃ ∥2F .
Using the above, the loss can indeed be decomposed as
L(W̃ , H̃) = LLS(H̃) + L⊥LS(W̃ , H̃), where
LLS(H̃) ≡ 1
2 Ave i,c
∥W̃LSh̃i,c − yi,c∥22 + λ
2 ∥W̃LS∥2F (14)
and
L⊥LS(W̃ , H̃) ≡ 1
2 Ave i,c
∥(W̃ − W̃LS)h̃i,c∥22 + λ
2 ∥W̃LS − W̃ ∥2F
= 1 2 tr { (W̃ − W̃LS) ( Σ̃T + µ̃Gµ̃ ⊤ G + λI ) (W̃ − W̃LS)⊤ } .
This completes the proof.

heading: C THEOREM 2
text: Theorem 2. (Decomposition of Least-Squares Component) The least-squares component, LLS(H̃), of the MSE decomposition in Theorem 1 can be further decomposed into LLS(H̃) = LNC1(H̃) + LNC2/3(H̃), where
LNC1(H̃) = 1 2 tr { W̃LS [ Σ̃W + λI ] W̃⊤LS } ,
LNC2/3(H̃) = 1
2C ∥W̃LSM̃ − I∥2F .
Proof. Under Euclidean distance, perturbations to the extended activations matrix H̃ which affect only the class-means or global-mean are orthogonal to perturbations which affect only the within-class covariance. Using this, LLS(H̃) can be further decomposed via the Pythagorean theorem:
LLS(H̃) = 1
2 Ave i,c
∥W̃LS(h̃i,c − µ̃c)∥22 + λ
2 ∥W̃LS∥2F +
1 2 Ave c ∥W̃LSµ̃c − yi,c∥22.
The first and second terms above merge into LNC1(H̃), while the third term becomes LNC2/3(H̃):
LNC1(H̃) = 1 2 tr { W̃LS ( Σ̃W + λI ) W̃⊤LS } (15)
LNC2/3(H̃) = 1
2C ∥W̃LSM̃ − I∥2F .
This completes the proof.
C.1 INTUITIONS FOR THEOREM 2 IN UNEXTENDED COORDINATES
For intuition, consider when λ = 0, i.e. the no weight-decay case. Proposition 1 in the main text gives
WLS = C −1 ĎM⊤Σ−1T , (16)
where ĎM ∈ RP×C is the matrix with columns µc − µG. Returning to the unextended coordinates, LNC1(ĎH) simplifies to
LNC1(ĎH) = 1 2 tr { WLSΣWW ⊤ LS } , (17)
where ĎH ∈ RP×CN has columns hi,c − µG. The term LNC2/3(ĎH) also simplifies instructively to
LNC2/3(ĎH) = 1
2C ∥WLS ĎM −Φ∥2F , (18)
where Φ ∈ RC×C is the standard Simplex ETF:
Φ = I − 1 C 11⊤.

heading: C.2 ADDITIONAL INTUITIONS FOR THEOREM 2 TERMS
text: The expressions in Equation 15 and Equation 17 further evoke the intuition that
LNC1(ĎH) is a variance term that goes to zero only under activation collapse (NC1) .
More specifically, we see that, with the class-means held constant, the only way to have LNC1(ĎH) → 0 is for ΣW → 0. Similarly, by examining Equation 18, we see that
LNC2/3(ĎH) quantifies deviations of WLS ĎM—i.e. the matrix of class-mean predictions—from the standard Simplex ETF.
Under (NC2) and (NC3) , both WLS and ĎM tend to jointly aligned ETF’s, possibly in an alternate pose; so WLS → ΦUT and ĎM → UΦ for a partial orthogonal matrix U satisfying UTU = IC . Since Φ2 = Φ, (NC2) and (NC3) together demand WLS ĎM → Φ. Thus, LNC2/3(ĎH) can be interpreted as a semi-metric reflecting the distance from achieving (NC2) and (NC3) .

heading: D THEOREM 3
text: D.1 IMPLICATIONS OF INVARIANCE
In this section, we will discuss in more detail how the invariance observed Section 3.1 leads the the continually renormalized gradient flow (Equation 5). In particular, we have seen in Section 3.1 of the main text that, on the central path, both the predictions WLS(ĎH)ĎH and the MSE loss
L(WLS(ĎH), ĎM(ĎH)) are invariant under the transformation ĎH 7→ AĎH when A is a invertible matrix. Since our interest mainly lies with A = Σ− 1 2
W , we will restrict the intuitive discussion in this section to positive-definite A.
Now, consider the set H of features ĎH with with non-singular within-class covariances. We view H as a fiber bundle17, where—on each fiber—live (apparently) different activations, ĎH , that (in fact) generate the exact same class predictions and the exact same MSE loss (see Figure 13). The base space—we will denote it X—can be taken to be the collection of X ∈ H where ΣW (X) = I . Every ĎH ∈ H is representable as ĎH = AX where A = Σ− 1 2
W ( ĎH) and ΣW (X) = I . The activation
matrices X in the base space might be variously called “pre-whitened,” “normalized,” “standardized,” or “sphered.” We shall call X a normalized features manifold18. Invoking invariance, we can make the following observations about our optimization task:
1. If, at a specific t, Ht happens to be a normalized set of activations (i.e. Ht ∈ X is currently located in the base space X ), there is no performance benefit to leaving the base space. There is also no performance benefit to moving along a fiber; only by moving from fiber-to-fiber can we improve performance.
2. In some sense, we waste time except when jumping from fiber-to-fiber; we might prefer to stay in the base space all the time by forcing the dynamics to jump from fiber-to-fiber.
On the other hand, we started with the viewpoint of studying gradient flow of the original H . So, we consider the following natural model for the application of gradient descent in X:
1. For a given initial feature activations, ĎH0 ∈ H, we renormalize—obtaining a starting point X0 ∈ X ; here X0 = Σ − 12 W (
ĎH0)ĎH0. Class predictions and MSE loss performance for the MSE-optimal classifier do not change from this renormalization.
2. We compute the usual gradient of MSE loss, at X0, obtaining a step ∆X0 = −η∇XLLS(X0), where η is a step size.
3. We take a step, going to ĎH1 = X0 +∆X0.
17An even larger space could be considered, where ΣW (ĎH) is not necessarily full rank, and then other components of this space having covariances of rank C, C − 1,...,1 could be considered. We just discuss the full-rank case here.
18This is sometimes called the Generalized Stiefel Manifold and will be defined more formally in Section D.2.
4. ĎH1 will not necessarily be normalized. We map back along whatever fiber we have landed upon, obtaining a corresponding point in the base space X , call this X1. Here,
X1 = Σ − 12 W ( ĎH1)ĎH1.
Class predictions and MSE loss performance do not change from this renormalization. 5. Repeat steps 2, 3, 4 at X1, obtaining thereby X2; and so on.
This process might be described as gradient descent with continual renormalization:
1. Renormalize the initial features activation matrix; 2. Compute the ordinary gradient of MSE loss and take the gradient descent step; 3. Renormalize again after each such step; and 4. Repeat steps 2 and 3.
While the continual renormalization process differs from usual gradient flow, it is both intuitively understandable and sensible. See discussion of assumption (A3) in Section 3

heading: D.2 ALIGNED SNR COORDINATES
text: Analysis of continually renormalized gradient flow can be simplified, without loss of generality, by applying a change of basis in both row and column space: Definition 3 (Transformation into Aligned SNR Coordinates). Consider the SVD decomposition of the SNR from Definition 2 outputting left singular vectors U ∈ RP×C and right singular vectors V ∈ RC×C . For any matrix Z ∈ RP×C , define the transformation Z → U⊤ZV as the transformation into aligned-SNR coordinates.
The aligned-SNR coordinates is so-named because it diagonalizes i.e. aligns the SNR matrix: Definition 4 (Aligned SNR Matrix). Consider the SVD decomposition of the SNR from Definition 2. Combined with Equation 8, observe that
Ω = diag ( {ωc}C−1c=1 , 0 ) = U⊤Σ − 12 W ĎMV
We call Ω the aligned SNR matrix.
The following facts about the aligned SNR matrix becomes useful later in our derivations: Observation 1 (SVD of Aligned SNR Matrix). The SVD of the aligned SNR matrix itself is simply
Ω = C−1∑ j=1 ωjeje ⊤ j ,
with the canonical basis {ej}Cj=1 as its singular vectors.
As discussed in Section 3.2, the {ωc}C−1c=1 that comprise Ω are decisive for understanding separation performance. In fact, when λ=0, the MSE loss can be entirely characterized by the SNR singular values on the central path: Lemma 1 (Spectral Representation of MSE Loss). When λ = 0, MSE loss obeys
L(WLS(ĎH), ĎH) = LLS(WLS(ĎH), ĎH) = 1
2 C−1∑ j=1
1
ω2j + C = L
( {ωj}C−1j=1 ) ,
on the central path, and its decomposition obeys
LNC1(ĎH) = 1
2 C−1∑ j=1 ω2j (C + ω2j ) 2 = LNC1 ( {ωj}C−1j=1 ) LNC2/3(ĎH) = 1
2 C−1∑ j=1 1 C
( ω2j
ω2j + C − 1
)2 = LNC2/3 ( {ωj}C−1j=1 ) .
Proof. Using Equation 16, Equation 17, and Equation 18, the loss on the central path equals
L(WLS, ĎM) = 1
2C ∥C−1 ĎM⊤Σ−1T ĎM −Φ∥ 2 F +
1 2 tr { C−1 ĎM⊤Σ−1T ΣWΣ −1 T ĎMC−1 } . (19)
Observe ĎM⊤Σ−1T ĎM and Φ are simultaneously diagonalizable since they share the same [1, . . . , 1] ⊤ null space and Φ has C − 1 equal eigenvalues. Also, observe that
ΣT = ΣW + 1
C ĎM ĎMT .
Let Ω̂ be submatrix of the aligned-SNR matrix Ω corresponding to the non-zero singular values i.e. Ω̂ = diag ( {ωj}C−1c=1 ) . After taking simultaneous-diagonalizations and canceling partial orthogonal
matrices, Equation 19 can be written solely in terms of Ω̂:
L(Ω̂) = 1 2C ∥∥∥∥Ω̂⊤ (Ω̂Ω̂⊤ + CIC−1)−1 Ω̂− IC−1∥∥∥∥2 F + 1 2 tr { Ω̂⊤(Ω̂Ω̂⊤ + CIC−1) −2Ω̂ } .
This can be further simplified into
L({ωj}C−1c=1 }) = 1
2 C−1∑ j=1 1 C
( ω2j
ω2j + C − 1
)2 +
ω2j (C + ω2j ) 2
= 1
2 C−1∑ j=1
1
ω2j + C . (20)
This completes the proof for Lemma 1.
Thus, the aligned SNR matrix, U⊤Σ− 1 2
W ĎMV , maximally simplifies analysis by reducing the vanilla
SNR matrix to a diagonal matrix (Definition 4) whose entries determine the MSE loss (Lemma 1).
We introduce the alignment of features into their own notion of aligned SNR coordinates: Definition 5 (Renormalization to aligned SNR Coordinates). Consider ĎH ∈ RP×CN with corresponding SNR left and right singular vectors U and V (Definition 2). Then, define the SNRaligned renormalization of the features as
X = U⊤ (ΣW ) − 12 ĎH(V ⊗ IN ) = U⊤ ( ĎHCĎH⊤ )− 12 ĎH(V ⊗ IN ) ∈ RC×NC ,
where ⊗ is the Kronecker product and C is a class-centering matrix defined as
C = 1
CN
( ICN − 1
N Y ⊤Y
) . (21)
As discussed in Section 3.1, the transformation ĎH → Σ− 1 2
W ĎH does not change the predictions of the
least squares classifier. Hence, it does not change the MSE loss. It is easy to check that this still holds for X i.e. invariance is preserved by the change into SNR-aligned coordinates.
It is also easy to check that the sphericity of the within-class covariance is preserved as well: Observation 2 (Sphericity of Within-Class Covariance of SNR Coordinates).
XCX⊤ = IC
Finally, routine applications of Definitions and canceling out of (partial) orthogonal matrix multiplications in aligned SNR coordinates lead to the following simplified representation for Ω and its differential: Observation 3 (Relation Between SNR Coordinates and SNR Matrix). The aligned SNR matrix, Ω (Definition 4) can be expressed as a function of the features in aligned SNR-coordinates X (Definition 5) as simply
Ω(X) = 1
N XY ⊤.
Moreover, let dΩij (·) : TXX → R be the differential 1-form19 associated with the ij-th entry of Ω. Define dΩ (·) : TXX → RC×C such that, for some Z ∈ TXX , each entry matrix of the matrix dΩ (Z) is dΩi,j (Z). Then, it follows that
dΩ (Z) = 1
N ZY ⊤, ∀Z ∈ TXX .
From a geometric perspective, the transformation in Definition 5 maps the features, H , to X belonging to a Normalized Features Manifold: Definition 6 (Normalized Features Manifold).
X = {X ∈ RC×CN |XCX⊤ = IC}
To understand the dynamics of the singular values of Ω—which is connected to the dynamics of the features through Lemma 1—we will analyze the gradient flow on this manifold. To this end, we first identify its tangent space at a particular X as well as the projection operator onto that tangent space. Proposition 3 (Tangent Space of Normalized Features Manifold).
TXX = {Z ∈ RC×CN |XCZ⊤ +ZCX = 0} Proposition 4 (Projection Onto Tangent Space of Normalized Features Manifold).
ΠTXX(Z) = Z − 1 2 (XCZ⊤ +ZCX⊤)X

heading: D.3 CONTINUALLY RENORMALIZED GRADIENT FLOW
text: In fact, gradient flow on the normalized features manifold (Definition 6) is the continuous analogue of the intuitive, discrete algorithm in Subsection D.1. In particular, the intuitive algorithm consists of a gradient step, ĎH1 = X0 +∆X0 followed by a mapping to the base space of the fiber bundle, X1 = Σ(ĎH1)
− 12 ĎH1. For small ∆X0, these two steps are in fact equivalent (up to a negligible term) to taking a TXX -projected gradient step on the manifold X as proven below: Lemma 2. Assuming X0 ∈ X , a renormalized gradient step,
X1 = Σ − 12 W ( ĎH1)ĎH1 = Σ − 12 W (X0 +∆X0) · (X0 +∆X0), is equivalent, up to an O ( ∥∆X0∥2 ) term, to a TXX -projected gradient step, i.e.
X1 =X0 +ΠTXX(∆X0) +O ( ∥∆X0∥2 ) .
Proof. Notice that
X1 =Σ − 12 W ( ĎH1)ĎH1
= ( (X0 +∆X0)C(X0 +∆X0) ⊤)− 12 (X0 +∆X0) = ( I +∆X0CX ⊤ 0 +X0C∆X ⊤ 0 +∆X0C∆X T 0
)− 12 (X0 +∆X0), where in the last step we used our assumption that X0 ∈ X , i.e., X0CX0 = I . By Taylor’s Theorem,
(I +A)− 1 2 = I − 1
2 A+O(∥A∥2).
Therefore, we get
X1 =
( I − 1
2
( ∆X0CX ⊤ 0 +X0C∆X ⊤ 0 +∆X0C∆X T 0 ) +O ( ∥∆X0∥2 )) (X0 +∆X0)
=X0 +∆X0 − 1
2
( ∆X0CX ⊤ 0 −X0C∆X⊤0 ) X0 +O ( ∥∆X0∥2 ) =X0 +ΠTXX(∆X0) +O ( ∥∆X0∥2 ) ,
where the last step follows from Proposition 4.
As we transition from a discrete gradient descent to a continuous gradient flow, the step size η → 0, and the residual O ( ∥∆X0∥2 ) in the above lemma becomes negligible, since ∆X0 =
−η∇XLLS(X0). This motivates us to study, in the next section, the continually renormalized gradient flow (as defined in Equation 5 of the main text) of X on X :
d dt X = −ΠTXX (∇XLLS(X)) .

heading: D.4 PROOF OF PROPOSITION 2 (GRADIENT FLOW IN ALIGNED SNR COORDINATES)
text: To prove Proposition 2, we first set up some auxiliary notation and lemmas. Denote the j-th row (transposed) of X and Y , respectively, by
xj = X ⊤ej (22)
and yj = Y ⊤ej .
Recall that—when differentiating SVDs—the derivative of the j-th singular value, ωj , only depends on the j-th singular subspace (see Equation 17 of Townsend (2016)). As a consequence, we will see in the following lemma that, for any differential dX , the corresponding dωj only depends on dxj .
Lemma 3 (Derivative of Singular Values With Respect to SNR Coordinates). Given the differential 1-form dωj (·) : TXX → R, the following holds for all19 dX ∈ RC×CN :
dωj (ΠTXX(dX)) =
( 1
N y⊤j − ωjx⊤j C
) dxj , ∀j = 1, . . . , C.
where (consistent with Definition 2) we adopt the convention that ωC=0 is the C-th singular value of the SNR matrix.
Proof. Without loss of generality, assume that X is represented in aligned SNR coordinates (Section D.2). Using Observation 1 and the differential of the singular value decomposition (see Equation 17 of Townsend (2016)):
dωj (ΠTXX(dX)) = e ⊤ j dΩ (ΠTXX(dX)) ej , (23)
where dΩ (·) is defined as in Observation 3. Next, Observation 3 implies
dΩ (ΠTXX(dX)) = 1
N ΠTXX(dX)Y
⊤. (24)
Using the projection onto the tangent space, given in Proposition 4 above,
ΠTXX(dX) = dX − 1 2 (dXCX⊤ +XC dX⊤)X. (25)
Combining Equation 23, Equation 24, and Equation 25, we obtain:
dωj (ΠTXX(dX)) =e ⊤ j dΩ (ΠTXX(dX)) ej (26)
= 1
N e⊤j ΠTXX(dX)Y ⊤ej
= 1
N e⊤j
( dX − 1
2
( dXCX⊤ +XC dX⊤ ) X ) Y ⊤ej .
19Notation: For real-valued functions f , we use df (·) : TXX → R to denote the associated differential 1-form i.e. the function outputting the directional derivative of f in the direction of the argument; For matrices Z ∈ Rm×n, we use dZ ∈ Rm×n—without succeeding parenthesis—to denote the matrix differential (cf. Townsend (2016)). These notions are equivalent: df (·) can be represented as a vector df in the dual-space of TXX ; When Z is a function of some vector v, dZ can be represented as a collection of 1-forms by defining dZij (·) = 〈 dZij dv , · 〉 for each entry Zij .
Moreover, using Observation 1 and 3, we can simplify these expressions into
1
N e⊤j dX Y
⊤ej = 1
N y⊤j dxj
1
N e⊤j dXCX ⊤XY ⊤ej =ωje ⊤ j dXCX ⊤ej = ωjx ⊤ j C dxj
1
N e⊤j XC dX ⊤ XY ⊤ej =ωje ⊤ j XC dX ⊤ ej = ωjx ⊤ j C dxj .
Finally, substituting the above three expressions into Equation 26, we get:
dωj (ΠTXX(dX)) =
( 1
N y⊤j − ωjx⊤j C
) dxj ,
which proves the claim.
Using Lemma 3, we can now obtain the dynamics of xj under gradient flow: Lemma 4. Under continually renormalized gradient flow (Equation 5),
d dt xj = ωj (C + ω2j ) 2
( 1
N yj − ωjCxj
) , ∀j = 1, . . . , C,
where (consistent with Definition 2) we adopt the convention that ωC=0 is the C-th singular value of the SNR matrix.
Proof. Recall xj ∈ RCN is the j-th row of X ∈ RC×CN , which we assume to be in aligned SNR coordinates without loss of generality (Section D.2). Applying the flow definition in Equation 5 to the i-th element of xj gives
d dt xji = −e⊤j ΠTXX ( dLLS dX ) ei
= − 〈 ΠTXX ( dLLS dX ) , eje ⊤ i 〉 F
= − 〈 dLLS dX ,ΠTXX ( eje ⊤ i )〉 F ,
where dLLSdX ∈ R C×CN is the standard derivative within the ambient RC×CN -space, ei ∈ RCN and ej ∈ RC are the canonical basis vectors, and ⟨·, ·⟩F is the Frobenius matrix inner-product. Next, observe that the differential 1-form dLLS (·) : TXX → R satisfies the following:
dLLS (Z) = 〈 dLLS dX ,Z 〉 F , ∀Z ∈ TXX .
Taking Z = ΠTXX ( eje ⊤ i ) then leads to
d dt xji = −dLLS
( ΠTXX ( eje ⊤ i )) = −
C−1∑ k=1 dLLS dωk dωk ( ΠTXX ( eje ⊤ i )) = −dLLS
dωj dωj
( ΠTXX ( eje ⊤ i )) where the second step follows by the chain rule, and the last step follows from Lemma 3 in that dωk ( ΠTXX ( eje ⊤ i )) = 0 if k ̸= j. Moreover, Lemma 3 also gives
dωj ( ΠTXX ( eje ⊤ i )) =
( 1
N y⊤j − ωjx⊤j C
) ei.
It then follows that
d dt xj =− dLLS dωj
( 1
N yj − ωjCxj ) =
ωj (C + ω2j ) 2
( 1
N yj − ωjCxj
) ,
where the last step follows from differentiating the expression in Lemma 1.
The gradient flows of xj induce the dynamics of ωj , which are described below.
Lemma 5 (Induced Dynamics on SVD of SNR). Under continually renormalized gradient flow (Equation 5), the dynamics for the j-th non-zero SNR singular value (Definition 2) is
d dt ωj = 1 N ωj (C + ω2j ) 2 , ∀j = 1, . . . , C − 1. (27)
Proof. Without loss of generality, assume X is represented in aligned SNR coordinates (Section D.2). Observation 1, the differential of the singular value decomposition (see Equation 17 of Townsend (2016)), and Observation 3 collectively imply that for any20 dX ∈ TXX ,
dωj (dX) = e ⊤ j dΩ (dX) ej =
1
N e⊤j dX Y
⊤ej = 1
N y⊤j dxj . (28)
Recall that the differential of ωj only depends on the differential of xj (see Lemma 3). Then, we can apply the chain rule:
d dt ωj = dωj dxj dxj dt = 1 N ωj (C + ω2j ) 2 y⊤j
( 1
N yj − ωjCxj
) , (29)
where the last equality results substituting-in Lemma 4 and Equation 28. Since Y Y ⊤ = NI ,
1
N2 y⊤j yj =
1
N .
Applying the same relation and Equation 21:
x⊤j Cyj = e ⊤ j XCY ⊤ej = 1
CN e⊤j X
( I − 1
N Y ⊤Y
) Y ⊤ej = 0.
Combining all the above equations, we obtain
d dt ωj = 1 N ωj (C + ω2j ) 2 .
The closed-form given in Proposition 2 now directly follows.
Proposition 2 (Dynamics of Singular Values of SNR Matrix). Continually renormalized gradient flow on the central path (Equation 5) induces the following closed-form dynamics on the SNR singular values (Definition 2):
c1 log(ωj(t)) + c2ω 2 j (t) + c3ω 4 j (t) = aj + t, t ≥ 0, for all j = 1, . . . , C − 1. (9)
c1, c2, and c3 are positive constants independent of j, and aj is a constant depending on ωj(0).
Proof. Follows from symbolically solving the ODE in Lemma 5 with routine methods. The constants are c1 = C2N , c2 = CN , and c3 = N4 .
20For our goal of applying the chain rule in Equation 29, we can restrict our attention to just dX ∈ TXX : The definition of renormalized gradient flow (Equation 5) ensures that dX
dt will always be in TXX .

heading: D.5 PROOF OF COROLLARY 1
text: Corollary 1 (Properties of SNR Singular Values). SNR singular values (Definition 2) following the Equation 9 dynamics satisfy the following limiting behaviors:
1. limt→∞ ωj(t) = ∞ and limt→∞ ωj(t)4√t/c3 = 1, for all j = 1, . . . , C − 1.
2. limt→∞ maxj ωj(t) minj ωj(t) = 1.
Proof. As t tends to infinity, the right-hand side of Equation 9 diverges to infinity and therefore so does the left-hand side (LHS). As the LHS approaches infinity, the logarithmic terms become negligible compared to the dominant quartic term, implying ω4j (t) → ∞. Since ωj(t) are singular values, they must be non-negative—implying ωj(t) → ∞. Based on the same argument, observe that limt→∞ ωj(t) 4 √
t c3
= 1 for all j. Since the constant c3 is independent of j, it follows that
limt→∞ maxj ωj(t) minj ωj(t) = 1.

heading: D.6 PROOF OF COROLLARY 2
text: Lemma 6. Under continually renormalized gradient flow (Equation 5), the left and right singular vectors of the SNR matrix remain constant i.e. they are independent of t.
Proof. Without loss of generality, assume X is in represented in aligned SNR coordinates (Section D.2). Recall that, in this coordinate system, the corresponding SNR matrix Ω = N−1XY ⊤ is diagonal, and the left and right singular vectors are simply partial identity matrices.
To show that the singular values of the SNR remain constant, it then suffices to show that dΩdt = dX dt Y is a diagonal matrix as well21. Lemma 4 gives d
dt X =
1
N D1Y −D2XC where D1 and D2 are C×C diagonal matrices with the diagonal entries {
ωj (C+ω2j ) 2 }C c=1
and{ ω2j
(C+ω2j ) 2 }C c=1 , respectively. Then,
d dt Ω =
( d
dt X
) Y ⊤ = 1
N D1Y Y
⊤ −D2XCY ⊤ = D1,
where, in the last equality, we used the easy-to-check identities that Y Y ⊤=NIC and CY =0. In more detail, the first identity follows from the fact that Y is the one-hot label vectors stacked as columns; the second identity follows from the facts that C is the class-centering matrix, and the one-hot label is the same for examples from the same class. Thus, we have shown dΩdt is diagonal, which concludes our proof.
Lemma 7. A matrix E ∈ RP×C is a Simplex ETF if and only if
1. E has exactly C−1 non-zero singular values, which are all equal.
2. E has a rank-1 nullspace spanned by the ones vector, i.e., E1C=0. In other words, E has zero-mean columns.
Proof. First, recall from Papyan, Han, and Donoho (2020, Definition 1) that a P×C matrix E is called a Simplex ETF if it satisfies
E⊤E = α
( C
C − 1 I − 1 C − 1 1C1 ⊤ C ) 21Intuitively, changing the diagonal SNR matrix by a diagonal matrix increment implies the updated SNR matrix is still diagonal. This, in turn, implies the left and right singular vectors of the updated SNR matrix are still the left and right partial identity matrices, respectively.
for some scaling α>0. We now prove the equivalence.
Simplex ETF implies 1-2: Consider the SVD decomposition E=UESEV ⊤E , where UE is a P×(C−1) partial orthogonal matrix satisfying U⊤EUE=I , SE is the (C−1)×(C−1) diagonal matrix of singular values, and V ⊤E is a (C−1)×C partial orthogonal matrix satisfying V ⊤E VE=I . Since E is a Simplex ETF, by definition,
E⊤E = α
( C
C − 1 I − 1 C − 1 1C1 ⊤ C
) ,
and according to the SVD decomposition
E⊤E = VESEU ⊤ EUESEV ⊤ E = VES 2 EV ⊤ E .
Therefore,
VES 2 EV ⊤ E = α
( C
C − 1 I − 1 C − 1 1C1 ⊤ C
) .
Notice that the right-hand-side has C−1 equal singular values. This implies SE—and, thus, E as well—has C−1 equal singular values. In other words, SE is a (C−1)×(C−1) diagonal matrix equal to diag(s, . . . , s) for some scalar s>0.
Next, notice that on the one hand
E⊤E1C = α
( C
C − 1 I − 1 C − 1 1C1 ⊤ C
) 1C = 0.
On the other hand,
E⊤E1C = VESEU ⊤ EUESEV ⊤ E 1C = VES 2 EV ⊤ E 1C = s 2VEV ⊤ E 1C .
Combining the above, we get VEV ⊤E 1C=0. Since VEV ⊤ E is a C×C matrix of rank C−1, we deduce that the rank-1 nullspace of V ⊤E is spanned by the C-dimensional ones vector, i.e. V ⊤ E 1C=0. Consequently, E has zero-mean columns since E1C=UESEV ⊤E 1C=0.
1-2 implies Simplex ETF: Assume E has exactly C−1 non-zero, equal singular values as well as a rank-1 nullspace spanned by the ones-vector, i.e., E1C=0. Then, E=UESEV ⊤E where
• UE is a P×(C−1) partial orthogonal matrix satisfying U⊤EUE=I;
• SE = diag(s, ..., s) is a (C−1)×(C−1) diagonal matrix for some scalar s>0; and
• V ⊤E is a (C−1)×C partial orthogonal matrix satisfying V ⊤E VE=I and also satisfying V ⊤E 1C=0.
Therefore, E⊤E = VESEU ⊤ EUESEV ⊤ E = VES 2 EV ⊤ E = s 2VEV ⊤ E .
Using our assumptions on VE ,
VEV ⊤ E +
1 C 1C1 ⊤ C = I,
where we divide each of the ones-vectors by √ C to create a unit vector. Thus, we conclude
E⊤E = s2 ( I − 1
C 1C1
⊤ C
) = α ( C
C − 1 I − 1 C − 1 1C1 ⊤ C
) ,
where α = s2C−1C . Hence, by definition, E is a Simplex ETF.
Corollary 2 (Neural Collapse Under MSE Loss). Under continually renormalized gradient flow (Equation 5), the SNR matrix (Equation 8) converges to
lim t→∞
1
ωmax(t) SNRt = Û0V̂
⊤ 0 , (10)
where Û0 ∈ RP×(C−1) and V̂0 ∈ RC×(C−1) are the left and right singular vectors of the SNR matrix (Definition 2) at t=0 corresponding to the non-zero singular values; and ωmax(t) is the
largest singular value at time t. Furthermore, Corollary 1 implies the occurrence of (NC1)-(NC4) i.e. renormalized gradient flow on the central path leads to Neural Collapse.
Moreover, denoting the Kronecker product with ⊗, the renormalized features matrix converges to
lim t→∞
1
ωmax(t) Σ
− 12 W,t ĎHt = (Û0V̂ ⊤ 0 )⊗ 1⊤N . (11)
Proof. Derivation of Equation 10: Corollary 1 proves the singular values j = 1, . . . , C − 1 of SNRt diverge to infinity and that their ratio tends to one. By Lemma 6, the renormalized gradient flow will not change the singular vectors of the SNR matrix, and—combined with the fact that the singular values converge to equality (second fact of Corollary 1)—we get the limit in Equation 10.
Derivation of (NC1): Let sj(·) denote the j-th singular value of its argument. Then, observe that tr ( Σ†B,tΣW,t ) = C tr (( ĎMt ĎM ⊤ t )† ΣW,t ) (Def. of ΣB,t)
= C tr ( Σ 1 2
W,t
( ĎMt ĎM ⊤ t )† Σ 1 2 W,t ) (Cyclic property of trace)
= C C−1∑ j=1 sj ( Σ 1 2 W,t ( ĎMt ĎM ⊤ t )† Σ 1 2 W,t ) (Trace is sum of singular values)
= C C−1∑ j=1 s−1j ( Σ − 12 W,t ( ĎMt ĎM ⊤ t ) Σ − 12 W,t ) (Def. of pseudoinverse)
= C C−1∑ j=1 s−2j ( Σ − 12 W,t ĎMt )
= C−1∑ j=1 1 ω2j (t) t→∞−−−→ 0. (Corollary 1)
By Horn & Johnson (2012, Theorem 1.3.22),
λj(Σ † B,tΣW,t) = λj(Σ 0.5 W,tΣ † B,tΣ 0.5 W,t) ≥ 0,
where λj(·) denotes the j-th eigenvalue of its argument, and the last inequality follows from the positive-semidefiniteness of Σ0.5W,tΣ † B,tΣ 0.5 W,t. Since the trace is the sum of eigenvalues, the only way for the trace of Σ†B,tΣW,t to tend to zero is if all eigenvalues also tend to zero. All eigenvalues tending to zero implies the matrix itself tends to zero, i.e.
lim t→∞
Σ†B,tΣW,t = 0,
which is the definition of (NC1) in Section 1.1.
Derivation of (NC2)-(NC4): Recall that the SNR matrix has zero-mean columns and rank C−1 (see Definition 2). This, combined with the fact that the singular values converge to equality (second fact of Corollary 1) imply, by Lemma 7, that the renormalized class-means converge to a Simplex ETF i.e. (NC2) .
From Theorem 1 of Papyan, Han, and Donoho (2020), we then know that (NC3) and (NC4) follow from (NC1) and (NC2) on the central path.
Derivation of Equation 11: Combining the limit in Equation 10 with (NC1) proves the limit in Equation 11.

heading: E RELATED WORKS EXAMINING NEURAL COLLAPSE
text: In this section, we discuss the contributions and limitations of seven recent works that propose and analyze theoretical abstractions of Neural Collapse. These works are only available in preprint, and may not yet be peer-reviewed. Thus, they might ultimately appear with very different claims or results. Additionally, works such as Poggio & Liao (2020a;b); Ergen & Pilanci (2020) also analyze behaviors other than NC; we will only discuss the parts relevant to Neural Collapse here.

heading: E.1 MIXON, PARSHALL, AND PI (2020)
text: Mixon et al. (2020) considered the unconstrained features model in Equation 2 (without weight decay) where, under gradient flow, (W ,H) evolve according to a nonlinear ordinary differential equation (ODE). They followed a two-step strategy for studying Neural Collapse. First, they linearized the ODE—claiming nonlinear terms are negligible for models initialized near the origin—and proved the simplified ODE converges to a subspace of (W ,H) satisfying (NC1) and (NC3) . Second, they proved that gradient flow, restricted to that subspace, converges to (NC2) .
The assumption of small weights and classifiers leading to the linearized ODE is not aligned with today’s paradigm. Specifically, the most commonly used He initialization (He et al., 2015) is designed: (i) to create weights with non-negligible magnitude; and (ii) to preserve the magnitude of features, as they propagate throughout the layers of the network, exactly so that last-layer features would have non-negligible magnitude. Moreover, the analysis of Mixon et al. essentially assumes that (NC1) and (NC3) occur much sooner than (NC2) . However, from the experiments in both Papyan, Han, and Donoho (2020) and this paper, there is no empirical evidence that (NC2) happens slower than (NC1) and (NC3) in practice.

heading: E.2 LU AND STEINERBERGER (2020)
text: While the MSE loss provides a mathematically natural setting for analysis, the modern paradigm in multi-class classification with deep learning involves training with CE loss, which is more challenging to analyze than MSE.
Lu & Steinerberger (2020) studied the (one-example-per-class) unconstrained22 features model with CE loss:
min W ,M
CE(W ,M) s.t. ∥wc∥2 = ∥µc∥2 = 1.
Since under linear separability the CE loss can be driven arbitrarily close to zero, just by re-scaling the norms of W and M , the authors further imposed a norm constraint on wc and µc. Lu & Steinerberger observe that the global minimizer of this optimization problem is only achieved once W and M are the same Simplex ETF. This derivation is suggestive, but it does not identify closed-form dynamics which would get gradient flow to such a global minimizer, nor does it address the rate of convergence to Neural Collapse. Additionally, the constraint on µc possesses no immediate or direct analogy to standard deep net training—where procedures often control the norm of the weights W , but not features H—nor class-means M .

heading: E.3 E AND WOJTOWYTSCH (2020)
text: E & Wojtowytsch (2020) also consider the unconstrained features model22 with CE loss,
min W ,H
CrossEntropy(WH) s.t. ∥W ∥2 ≤ 1, ∥hi,c∥2 ≤ 1,
where they adopt a more technical, spectral norm constraint on W to specify their model.
Building on the results of Chizat & Bach (2018; 2020), E & Wojtowytsch also construct a simple counter-example showing that Neural Collapse need not occur in two-layer, infinite-width networks— which have been the focus of intense recent study in the theoretical deep learning community (Mei et al., 2018; Rotskoff & Vanden-Eijnden, 2018; Arora et al., 2019). Thus, E & Wojtowytsch’s counterexample suggests the alternative perspective that, despite the expressiveness of infinite-width, two-layer networks, such abstractions do not capture key aspects of trained deep nets.
As with Lu & Steinerberger (2020), standard deep net training does not possess any direct analogies for constraining the norm of features (as opposed to weights)—nor are there any paradigmatic regularizations that correspond to controlling the spectral norm on W . Moreover, the work does not characterize any closed-form dynamics or the rate of convergence to Neural Collapse.
22The model is unconstrained in the sense that the features are allowed to move directly with gradient flow and are not constrained to be the output of a forward pass—not in the sense that there are no constraints on the optimization problem.

heading: E.4 POGGIO & LIAO (2020A;B) (WITH BANBURSKI)
text: Distinguished from the simplified unconstrained features models in the previously mentioned works is the theoretical analysis of Poggio & Liao (2020a;b) (in a special section, co-authored with Andrzej Banburski).
The authors study deep homogeneous classification networks, with weight normalization layers, trained with stochastic gradient descent and weight decay. This is much closer to today’s training paradigm, but the setting still differs from the one in which Neural Collapse has been empirically observed in Papyan, Han, and Donoho (2020) and in Section A of this paper. In particular, they replace batch normalization with weight normalization and consider deep homogeneous networks; homogeneous networks can not have bias vectors nor skip connections, which are present both in ResNet and DenseNet. Moreover, the work gives explicit descriptions of neither the dynamics nor the rate of convergence to Neural Collapse.
E.5 ERGEN & PILANCI (2020)
While the above-described works tend to focus on either the used-in-practice CE loss or the theoretically-insightful MSE loss, Ergen & Pilanci (2020) observed that these are both instances of the general class of convex loss functions and, thus, one could derive insights from the classical convex analysis literature. Moreover, compared to Mixon et al. (2020); Lu & Steinerberger (2020); E & Wojtowytsch (2020), this work studies the optimization starting from the second-to-last layer features rather than the last-layer features. In particular, the authors use a strong-duality argument to show that NC emerges in the optimal solution of an equivalent proxy-model to the following optimization:
min HL−1,WL−1,WL,γ,α
L ( WL (BNγ,α (WL−1HL−1))+ ,Y ) + λ
2
( ∥γ∥22 + ∥α∥ 2 2 + ∥WL∥ 2 F ) ,
where L(·) is a general convex loss, HL−1 are the second-to-last layer activations, WL−1 are the second-to-last layer weights, WL are the network classifiers, Y are the training targets, λ is a weight-decay parameter, BNγ,α(·) is a batch-norm operator parameterized by α and γ, and (·)+ is a ReLU. The incorporation of batch-normalization and weight-decay ensures the existence of bounded, well-defined optimal solutions—serving a similar role to that of weight-normalization and weight decay in Poggio & Liao (2020a;b) as well as the norm constraints in the other aforementioned related works.
Since strong-duality only characterizes properties of the converged optimal solution of an optimization model, Ergen & Pilanci (2020) does not provide insights into the dynamics with which that solution is achieved which training.
E.6 FANG, HE, LONG, AND SU (2021)
In Fang et al. (2021), the authors introduce the (N -)layer-peeled model in which one considers only the direct optimization of the N -th-to-last layer features of a deep net along with the weights that come after the N -th-to-last layer. The motivating philosophy is that, after raw inputs are passed through some initial number of layers, the overparameterization of those layers would allow us to effectively model the N -th-to-last layer features as freely-moving in some subset of Euclidean space. In this terminology, the concurrent works of Mixon et al. (2020); Lu & Steinerberger (2020); E & Wojtowytsch (2020) on the unconstrained (last-layer) features model could be considered instances of a 1-layer-peeled model; while the model of Ergen & Pilanci (2020) could be considered as a 2-layer-peeled model. This perspective is attractive because it gives a name and organization to a common modeling philosophy behind the above-described body of independent works. For comparison, the work of Poggio & Liao (2020a;b) is a non-example of layer-peeled modeling as it considers optimization on the weights of homogeneous deep nets and not the input features.
Fang et al. (2021) then analyzes a convex relaxation of the 1-layer-peeled model—with norm constraints on the weights and features—into a semidefinite program. Not only do the authors show that this model exhibits Neural Collapse in the canonical setting of balanced examples-per-class, but they also analyze the behavior of this model under imbalanced classes. While, in the imbalanced case, one would intuitively expect the Simplex ETF to “skew” to have bigger angles around overrepresented classes and smaller angles around under-represented ones; Fang et al. (2021) identifies the
surprising phenomenon—named minority collapse—in their model where, when the imbalances pass a certain threshold, the last-layer features and classifiers of the under-represented classes collapse to be exactly the same. However, their work does not provide closed-form dynamics or rates at which collapse—in either the balanced or imbalance case—occurs.

heading: E.7 ZHU, DING, ZHOU, LI, YOU, SULAM, AND QU (2021)
text: In Zhu et al. (2021), the authors examine the following unconstrained features model:
min W ,H,b CrossEntropy (WH + b,Y ) + λW 2 ∥W ∥2F + λH 2 ∥H∥2F + λb 2 ∥b∥2F ,
where (W , b) are the classifier weights and biases, H are the last layer features, and (λW , λH , λb) are weight-decay parameters. On this model, the authors not only prove that all minima exhibit Neural Collapse but also that all local minima are global minima. In comparison to our current paper, Zhu et al. (2021) focus on characterizing the landscape of the loss and, thus, do not explore the dynamics and rate at which such minimizers of the loss are achieved.
Zhu et al. (2021) also make notable empirical contributions by conducting a series of experiments on the MNIST and CIFAR10 datasets trained on MLPs, ResNet18, and ResNet50. Their measurements give evidence for the following novel NC-related phenomena in deep classification networks:
1. NC is algorithm independent: NC emerges in realistic classification deep net training regardless of whether the algorithm is SGD, ADAM, or L-BFGS.
2. NC occurs on random labels: NC emerges even when the one-hot target vectors are completely shuffled.
3. Width improves NC: Increasing network width expedites NC when training with random labels.
The authors of Zhu et al. (2021) moreover conducted ablation experiments suggesting that the following substitutions can be made to deep neural net architectures without affecting performance:
1. Weight-decay substitution: Replacing (A) weight-decay on the norm of all network parameters with (B) weight-decay just on the norm of the last-layer features and classifiers.
2. Classifier substitution: Replacing (A) the last-layer classifiers that are trained with SGD with (B) Simplex ETF classifiers that are fixed throughout training.
While the authors only demonstrated these new behaviors on limited network-dataset combinations, these experiments indeed inspire interesting conjectures about the generalization behavior of deep nets as well as potential architecture design improvements; Zhu et al. (2021) discuss many of these conjectures and related open-questions in detail.

Review summaries:
This paper extends the recent work on Neural Collapse, using Mean Squared Error (MSE) instead of CE, as MSE is easier for analysis. With this, the paper shows that the least square loss can be decomposed into one that corresponds to a so called 'central' path (namely a set of optimal tuples (W, b, H) given H for which there is an optimal loss), and a perpendicular loss. The paper shows that the perpendicular loss is much smaller than the optimal least square loss and thus Neural Collapse appears due to the fact that the optimizer focuses on the central path. Empirical studies confirm the findings.

Recently, Papyan, Han and Donoho (2020) have observed that training neural networks beyond zero training error leads to simplex arrangements of the features. The submission studies this phenomenon, called neural collapse, from a theoretical perspective, when training is done via minimizing the square loss.  In particular, it is shown that the square loss can be split, such that one summand corresponds to the quality of the features (quantified by the loss of the MSE classifier on them) and that the other summand corresponds to the quality of the classification layer (quantified by the its deviationfrom the least squares classifier). Furthermore, the first summand quantifies the closeness of the features to neural collapse. In a second step, the authors investigate the gradient flow induced from the first summand and derive closed form solutions.

This paper studies the phenomenon of Neural Collapse (NC) and empirically shows that it occurs during the training of deep networks with the MSE loss. Then it theoretically analyzes NC with MSE loss by decomposing it and introducing the notion of central path. It shows a closed-form dynamics predicts NC in this setting.