# Average gradient outer product as a mechanism for deep neural collapse

Daniel Beaglehole  Peter Súkeník  Marco Mondelli  Mikhail Belkin

###### Abstract

Deep Neural Collapse (DNC) refers to the surprisingly rigid structure of the
data representations in the final layers of Deep Neural Networks (DNNs).
Though the phenomenon has been measured in a wide variety of settings, its
emergence is only partially understood. In this work, we provide substantial
evidence that DNC formation occurs primarily through deep feature learning
with the average gradient outer product (AGOP). This takes a step further
compared to efforts that explain neural collapse via feature-agnostic
approaches, such as the unconstrained features model. We proceed by providing
evidence that the right singular vectors and values of the weights are
responsible for the majority of within-class variability collapse in DNNs. As
shown in recent work, this singular structure is highly correlated with that
of the AGOP. We then establish experimentally and theoretically that AGOP
induces neural collapse in a randomly initialized neural network. In
particular, we demonstrate that Deep Recursive Feature Machines, a method
originally introduced as an abstraction for AGOP feature learning in
convolutional neural networks, exhibits DNC.

neural collapse, recursive feature machines, neural networks

##  1 Introduction

How Deep Neural Networks (DNNs) learn a transformation of the data to form a
prediction and what are the properties of this transformation constitute
fundamental questions in the theory of deep learning. A promising avenue to
understand the hidden mechanisms of DNNs is Neural Collapse (NC) (Papyan et
al., 2020). NC is a widely observed structural property of overparametrized
DNNs, occurring in the terminal phase of gradient descent training on
classification tasks. This property is linked with the performance of DNNs,
such as generalization and robustness (Papyan et al., 2020; Su et al., 2023).
NC is defined by four properties that describe the geometry of feature
representations of the training data in the last layer. Of these, in this
paper we focus on the following two, because they are relevant in the context
of deeper layers. The within-class variability collapse (NC1) states that
feature vectors of the training samples within a single class collapse to the
common class-mean. The orthogonality or simplex equiangular tight frame
property (NC2) states that these class-means form an orthogonal or simplex
equiangular tight frame. While initially observed in just the final hidden
layer, these properties were measured for the intermediate layers of DNNs as
well (Rangamani et al., 2023; He & Su, 2022), indicating that NC progresses
depth-wise. This has lead to the definition and study of the Deep Neural
Collapse (DNC), a direct depth-wise extension of the standard Neural Collapse
(Súkeník et al., 2023).

The theoretical explanations of the formation of (D)NC have mostly relied on a
data-agnostic model, known as unconstrained features model (UFM) (Mixon et
al., 2020; Lu & Steinerberger, 2022). This assumes that DNNs are infinitely
expressive and, thus, optimizes for the feature vectors in the last layer
directly. The deep UFM (DUFM) was then introduced to account for more layers
(Súkeník et al., 2023; Tirer & Bruna, 2022). The (D)UFM has been identified as
a useful analogy for neural network collapse, as (D)NC emerges naturally in
this setting for a variety of loss functions and training regimes (see Section
2). The (D)NC has been proven to emerge or be optimal under (D)UFM in various
settings – under different loss functions, training regimes and in different
senses, see Section 2. However, (D)UFM discards the training data and most of
the network completely, thus ignoring also the role of feature learning in
DNC. This is a serious gap to fully understand the formation of DNC in the
context of the full DNN pipepline.

In this paper, we argue that feature learning is responsible for inducing DNC
in trained neural networks. To argue this, we first give substantial evidence
that within-class variability collapse in DNNs is primarily caused by the
application of the right singular vectors of the weight matrix and the
subsequent multiplication with its singular values (Section 4). This is in
opposition to the belief that neural collapse improves through layers as a
consequence of the application of the non-linearity (Tirer & Bruna, 2022).
Indeed, if $W$ is a full rank matrix, then the common metric of within-class
variability for $X$ (the neural representations that are multiplied by $W$)
and $WX$ are equal. Thus, the only way in which the application of a linear
layer reduces the within-class variability is by erasing the within-class
differences via the small singular directions of $W$. This operation can be
viewed as discarding irrelevant information, or denoising, a form of feature
learning.

The small singular directions of $W$ are defined by its right singular vectors
and singular values. Therefore, it is fully deducible from the Gram matrix of
the weights at each layer, $W^{\top}W.$ Radhakrishnan et al. (2022);
Beaglehole et al. (2023) have identified that, in many settings and across all
layers of the network, $W^{\top}W$ is approximately proportional to the
average gradient outer product (AGOP) with respect to the inputs to that
layer, in a statement termed the Neural Feature Ansatz (NFA). The NFA suggests
that neural networks extract features from the data representations at every
layer by projection onto the AGOP with respect to that layer.

Recent work then incorporated layerwise linear transformation with the AGOP
into kernel machines, such as the neural tangent kernel (Jacot et al., 2018),
to model the deep feature learning mechanism of neural networks (Beaglehole et
al., 2023). The output of their backpropagation-free method, Deep Recursive
Feature Machines (Deep RFM), is a standard, randomly initialized neural
network, where each random weight matrix is additionally right-multiplied by
the AGOP of a kernel machine. Their method was shown to improve performance of
convolutional kernels on vision datasets and recreate the edge detection
ability of convolutional neural networks.

Strikingly, we show that the neural network generated by this very same
method, Deep RFM, exhibits deep neural collapse (Section 5). We link the
mechanism of DNC formation between DNNs and Deep RFM by demonstrating that
projection onto the AGOP is responsible for most of the DNC metric improvement
in Deep RFM both empirically and theoretically. Our results demonstrate that
_(i)_ AGOP is the mechanism of DNC in Deep RFM, and _(ii)_ when the NFA holds,
the right singular structure of $W$, and therefore the AGOP, induces within-
class variability collapse in DNNs. We thus establish feature learning through
the AGOP as a plausible mechanism for DNC in DNNs.

We verify these claims with an extensive experimental evaluation, which
demonstrates a consistent formation of DNC in Deep RFMs across several vision
datasets. Additionally, we provide a theoretical analysis that explains this
mechanism in the asymptotic high-dimensional regime and the smaller data
regime under additional assumptions. The first analysis is based on the
approximate linearity of kernel matrices in the regime of large number of data
samples and dimension. We argue in our second analysis that the formation of
DNC in Deep RFMs is optimal in the kernel ridge regression sense. Further, if
the data is already partially grouped into clusters (as assumed in (Seleznova
et al., 2023)), we show how the AGOP can strengthen the DNC from one layer to
the next.

##  2 Related work

#### Neural collapse (NC).

Since the seminal paper by Papyan et al. (2020), the phenomenon of neural
collapse has attracted significant attention. Hui et al. (2022) argue that NC
does not generalize to the test set. Galanti et al. (2022) use NC to improve
generalization bounds for transfer learning, while Wang et al. (2023) discuss
transferability capabilities of pre-trained models based on their NC on the
target distribution. Haas et al. (2022) argue that NC can lead to improved OOD
detection. Connections between robustness and NC are discussed in (Su et al.,
2023). For a survey, we refer the interested reader to (Kothapalli, 2023).

The main theoretical framework to study the NC is the Unconstrained Features
Model (UFM) (Mixon et al., 2020; Lu & Steinerberger, 2022). Under this
framework, many works have shown the emergence of the NC, either via
optimality and/or loss landscape analysis (Weinan & Wojtowytsch, 2022; Lu &
Steinerberger, 2022; Zhou et al., 2022), or via the study of gradient-based
optimization (Ji et al., 2022; Han et al., 2022; Wang et al., 2022). Some
papers also focus on the generalized class-imbalanced setting, where the NC
does not emerge in its original formulation (Thrampoulidis et al., 2022; Fang
et al., 2021; Hong & Ling, 2023). Deviating from the strong assumptions of the
UFM model, a line of work focuses on gradient descent in homogeneous networks
(Poggio & Liao, 2020; Xu et al., 2023; Kunin et al., 2022), while others
introduce perturbations (Tirer et al., 2022).

Crucially, it was soon noticed that NC propagates to earlier layers of DNNs,
as opposed to just being a phenomenon of the last layer. He & Su (2022);
Rangamani et al. (2023); Parker et al. (2023) show the presence of the DNC in
DNNs in various training regimes. Theoretical works have then tried to extend
the theory to reach beyond the last layer. A deep linear UFM is considered by
(Dang et al., 2023) and the optimality of DNC is proved under various regimes.
The first work to study deep non-linear ReLU-activated UFM is by Súkeník et
al. (2023), where the authors show the optimality of DNC across the whole
depth for the deep UFM in binary classification.

#### AGOP feature learning.

The NFA was shown to capture many aspects of feature learning in neural
networks including improved performance and interesting structural properties.
The initial work on the NFA connects it to the emergence of spurious features
and simplicity bias, why pruning DNNs may increase performance, the “lottery
ticket hypothesis”, and a mechanism for grokking in vision settings
(Radhakrishnan et al., 2022). Zhu et al. (2023) connect large learning rates
and catapults in neural network training to better alignment of the AGOP with
the true features. Radhakrishnan et al. (2024) demonstrate that the AGOP
recovers the implicit of bias of deep linear networks toward low-rank
solutions in matrix completion. Beaglehole et al. (2024) identify that the NFA
is due to alignment of the weights to the pre-activation tangent kernel and
introduce an intervention to improve NFA feature learning.

##  3 Preliminaries

###  3.1 Notation

For simplicity, we assume a class-balanced setting, where $N=Kn$ is the total
number of training samples, with $K$ being the number of classes and $n$ the
number of samples per class. Let us order the training samples such that the
labels $Y\in\mathbb{R}^{K\times N}$ can be written as
$I_{K}\otimes\mathbf{1}_{n}^{\top},$ where $I_{K}$ denotes a $K\times K$
identity matrix, $\otimes$ denotes the Kronecker product and $\mathbf{1}_{n}$
a row vector of all-ones of length $n$. For a matrix $A$ and a column vector
$v$, the operation $A\odiv v$, divides each row of $A$ by the corresponding
element of $v$.

Both a DNN and Deep RFM of depth $L$ can be written as:

$\displaystyle
f(x)=W_{L+1}\sigma(m_{L}\sigma(m_{L-1}\dots\sigma(m_{1}(x))\dots)),$

where $m_{l}$ is a matrix and $\sigma$ is the ReLU activation function. If the
matrix $m_{l}$ is defined by a single matrix, as in fully-connected networks,
and not a product, we denote it by a weight matrix $W_{l}$. Each weight matrix
is decomposed via its SVD as $W_{l}=U_{l}S_{l}V_{l}^{\top}.$ The training data
$X_{1}\in\mathbb{R}^{d_{1}\times N}$ is stacked into columns and we let
$X_{l},\tilde{X}_{l},l\geq 2$ be the feature representations of the training
data after $l$ layers of the DNN after and before the non-linearity,
respectively.

###  3.2 Deep Neural Collapse

We define the feature vectors $x^{l}_{ci}$ of the $i$-th sample of the $c$-th
class as the input to $m_{l}$. Then, $\tilde{x}^{l}_{ci}$ are the feature
vectors produced from $m_{l-1}$ before the application of $\sigma$, and are
defined only for $l>1$. Let
$\mu_{c}^{l}:=\frac{1}{n}\sum_{i=1}^{n}x_{ci}^{l}$,
$\tilde{\mu}_{c}^{l}:=\frac{1}{n}\sum_{i=1}^{n}\tilde{x}_{ci}^{l}$ be the
corresponding class means, and $H_{l},\tilde{H}_{l}\in\mathbb{R}^{d_{l}\times
K}$ the matrices of class means. We will also consider the centered and
normalized class means
$\bar{H}_{l}=\left(H_{l}-\mu^{l}_{G}\right)\odiv\left(H_{l}-\mu^{l}_{G}\right)$,
where $\mu^{l}_{G}$ is the global feature mean at layer $l.$ NC1 is defined in
terms of two metrics:

$\displaystyle\Sigma_{W}^{l}$
$\displaystyle=\frac{1}{N}\sum_{c=1}^{K}\sum_{i=1}^{n}(x_{ci}^{l}-\mu^{l}_{c})(x_{ci}^{l}-\mu^{l}_{c})^{\top},$
$\displaystyle\Sigma_{B}^{l}$
$\displaystyle=\frac{1}{K}\sum_{c=1}^{K}(\mu^{l}_{c}-\mu^{l}_{G})(\mu^{l}_{c}-\mu^{l}_{G})^{\top}~{},$

i.e., the _within-class_ and _between-class_ variability at layer $l$.

######  Definition 3.1.

In the context of this work, the neural collapse at layer $l$ is achieved if
the following two properties are satisfied:

  * •

DNC1: The within-class variability at layer $l$ is $0$. This property can be
stated for the features either after or before the application of the
activation function $\sigma$. In the former case, the condition requires
$x^{l}_{ci}=x^{l}_{cj}$ for all $i,j\in\\{1,\ldots,n\\}\triangleq[n]$ (or, in
matrix notation, $X_{l}=H_{l}\otimes\mathbf{1}_{n}^{\top}$); in the latter,
$\tilde{x}^{l}_{ci}=\tilde{x}^{l}_{cj}$ for all $i,j\in[n]$ (or, in matrix
notation, $\tilde{X}_{l}=\tilde{H}_{l}\otimes\mathbf{1}_{n}^{\top}$).

  * •

DNC2: The class-means at the $l$-th layer form either an orthogonal basis or a
simplex equiangular tight frame (ETF). As for DNC1, this property can be
stated for features after or before $\sigma$. We write the ETF class-mean
covariance
$\Sigma_{\mathrm{ETF}}\triangleq(1+\frac{1}{K-1})I_{K}-\frac{1}{K-1}\mathbf{1}_{K}\mathbf{1}_{K}^{\top}.$
In the former case, the condition requires either $H_{l}^{\top}H_{l}\propto
I_{K}$ or $H_{l}^{\top}H_{l}\propto\Sigma_{\mathrm{ETF}}$; in the latter
$\tilde{H}_{l}^{\top}\tilde{H}_{l}\propto I_{K}$ or
$\tilde{H}_{l}^{\top}\tilde{H}_{l}\propto\Sigma_{\mathrm{ETF}}.$

Figure 1: Feature variability collapse from different singular value
decomposition components on CIFAR-10. We measure the reduction in the NC1
metric throughout training at each layer of an MLP with $L=5$. Each layer is
decomposed into its input, $\Phi(X)$, projection onto the right singular space
of $W$, $SV^{\top}\Phi(X)$, and then $U$, the left singular vectors of $W$,
and application of ReLU.

###  3.3 Average gradient outer product

The AGOP is a statistic defined with respect to a dataset
$X\in\mathbb{R}^{d_{0}\times N}$ and any model, that accepts inputs from the
data domain $f:\mathbb{R}^{d_{0}\times 1}\rightarrow\mathbb{R}^{K\times 1}$,
where $K$ is the number of outputs. By writing the Jacobian of the model
outputs with respect to its inputs at layer $l$ as $\frac{\partial
f(x)}{\partial x^{l}}\in\mathbb{R}^{d_{0}\times K}$, we define the AGOP as an
operator on $f$ and $X$:

$\displaystyle\operatorname{AGOP}(f,X)\triangleq\frac{1}{N}\sum_{c=1}^{K}\sum_{i=1}^{N}\frac{\partial
f(x_{ci})}{\partial x}\frac{\partial f(x_{ci})}{\partial x}^{\top}~{}.$

This object has important implications for learning, in part, because the AGOP
of a predictor that has somewhat learned the target function will resemble the
expected gradient outer product (EGOP) of the target function (Yuan et al.,
2023). While the AGOP is specific to a model and training data, the EGOP is an
operator on the population distribution and the function to be estimated. The
EGOP contains specific useful information about the target and data
distribution, such as low-rank structure, that improves prediction.

Remarkably, it was identified in (Radhakrishnan et al., 2022; Beaglehole et
al., 2023) that neural networks will automatically estimate the AGOP in the
Gram matrix of the weights in all layers of the neural network. Stated
formally, the authors observe the Neural Feature Ansatz, here posed for fully-
connected architectures.

######  Ansatz 3.2 (Neural Feature Ansatz (Radhakrishnan et al., 2022)).

Let $f$ be a depth-$L$ neural network trained on data $X$ using stochastic
gradient descent. Then, for all layers $l\in[L]$,

$\displaystyle
W_{l}^{\top}W_{l}\,\propto\,\frac{1}{N}\sum_{c=1}^{K}\sum_{i=1}^{n}\frac{\partial
f(x_{ci})}{\partial x^{l}}\frac{\partial f(x_{ci})}{\partial x^{l}}^{\top}.$

Crucially, the AGOP can be defined for any differentiable estimator, not
necessarily a neural network. For example, prior work considers the AGOP of
general non-parametric methods (Trivedi et al., 2014). In (Radhakrishnan et
al., 2022), the authors introduce a kernel learning algorithm, the Recursive
Feature Machine (RFM), that performs kernel regression and estimates the AGOP
in an alternating fashion, modeling the NFA feature learning mechanism of
fully-connected networks.

Subsequent work introduced the Deep Recursive Feature Machine (Deep RFM) to
model deep feature learning in neural networks (Beaglehole et al., 2023). Deep
RFM recursively generates representations by transforming the input to that
layer with the AGOP of the model w.r.t. this input, itself estimated from
shallow RFM, and then applying a random feature map. To define the Deep RFM,
let $\\{k_{l}\\}_{l=1}^{L+1}:\mathbb{R}^{d_{l}\times
d_{l}}\xrightarrow[]{}\mathbb{R}$ be a set of kernels. For the ease of
notation, we will write $k_{l}(X_{l},X_{l})$ for a matrix of kernel
evaluations between columns of $X_{l}.$ Then, the Deep RFM is described in
Algorithm 1.

Algorithm 1 Deep Recursive Feature Machine (Deep RFM)

0: $X_{1},Y,\\{k_{l}\\}_{l=1}^{L+1},L,\\{\Phi_{l}\\}_{l=1}^{L}$ {kernels:
$\\{k_{l}\\}_{l}$, depth: $L$, feature maps: $\\{\Phi_{l}\\}_{l}$}

0: $\alpha_{L+1},\\{M_{l}\\}_{l=1}^{L}$

for $l\in 1,\dots,L$ do

Normalize the data, $X_{l}\xleftarrow[]{}X_{l}\odiv\left\|X_{l}\right\|$

Learn coefficients, $\alpha_{l}=Y(k_{l}(X_{l},X_{l})+\mu I)^{-1}$.

Construct predictor, $f_{l}(\cdot)=\alpha_{l}k_{l}(\cdot,X_{l})$.

Compute AGOP: $M_{l}=\sum_{c,i=1}^{K,n}\nabla f_{l}(x^{l}_{ci})\nabla
f_{l}(x^{l}_{ci})^{\top}$.

Transform the data $X_{l+1}\xleftarrow[]{}\Phi_{l}(M_{l}^{1/2}X_{l})$.

end for

Normalize the data,
$X_{L+1}\xleftarrow[]{}X_{L+1}\odiv{\left\|X_{L+1}\right\|}$

Learn coefficients, $\alpha_{L+1}=Y(k_{L+1}(X_{L+1},X_{L+1})+\mu I)^{-1}$.

Figure 2: Neural collapse with Deep RFM on 2000 points from CIFAR-10. The
$2000\times 2000$ matrix of inner products of all pairs of points in $X_{l}$
extracted from layers $l\in\\{1,\ldots,19\\}$ of Deep RFM following centering
by the global mean $\mu^{l}_{G}$ and normalizing each point to the unit
sphere. The data are ordered in blocks so that points of the same class are
adjacent to one another, arranged from class $1$ to class $10$. All points of
the same class have inner products approximately equal to $1$ at the final
layer $l=19$ (are equal to each other), while points of different classes have
inner products approximately $-\frac{1}{K-1}$, for each of the $K=10$ classes.

Note that the Deep RFM as defined here is less general than the version in
(Beaglehole et al., 2023), as we are considering only one AGOP estimation step
in the inner-optimization of RFM. The high-dimensional feature maps
$\Phi_{l}(\cdot)$ are usually realized as $\sigma(W_{l}\cdot),$ where $W_{l}$
is a matrix with standard Gaussian entries, thus $\Phi_{l}$ serves as a random
features generator.

The single loop in the Deep RFM represents a reduced RFM learning process. The
RFM itself is based on kernel ridge regression, therefore we introduce the
ridge coefficient $\mu.$ The training procedure of RFM is meant to perform
well in the corresponding kernel ridge regression which will be particularly
important in Section 5.4.

##  4 Deep Neural Collapse as an effect of linear denoising

There are precisely two ways that DNNs can discard unnecessary information in
the data, such as noise or spurious features. One of them is through non-
linearities, i.e., ReLU, where negative entries are erased. The other one is
through the small eigendirections of $W_{l}^{\top}W_{l}$ or, approximately,
the non-trivial null space of $W_{l}$, where components of vectors that belong
to the null space are discarded. If we write again the SVD
$W_{l}=U_{l}S_{l}V_{l}^{\top}$, then $S_{l}V_{l}^{\top}$ is the part with
possibly non-trivial null space. This right singular component can be viewed
as the feature extraction part of the linear layer. In contrast, the
multiplication with $U_{l}$ simply rotates or reflects the prior
representation.

All parts of the process, together with the ReLU, are important for neural
network predictions. However, the DNC1 can only improve through the
$S_{l}V_{l}^{\top}$ multiplication and the ReLU application. In particular,
considering the standard DNC1 metric $tr(\Sigma_{W})/tr(\Sigma_{B})$
(Rangamani et al., 2023), a simple computation gives

$\frac{tr(\hat{\Sigma}_{W})}{tr(\hat{\Sigma}_{B})}=\frac{tr(U\Sigma_{W}U^{\top})}{tr(U\Sigma_{B}U^{\top})}=\frac{tr(\Sigma_{W})}{tr(\Sigma_{B})}.$

Here, the $\Sigma$ matrices refer to the within-class variability after the
application of $SV^{\top},$ while the $\hat{\Sigma}$ matrices correspond to
the output of the full linear layer. The first equality holds due to the
linearity and the second due to the cyclic property of the trace.

While both ReLU and $S_{l}V_{l}^{\top}$ can influence the DNC1 metric, our
experiments show that $S_{l}V_{l}^{\top}$ induces the majority of within-class
variability collapse. We verify this claim by plotting the DNC1 metrics of all
layers of a fully-connected network trained on CIFAR-10 (Figure 1, see
Appendix B for similar results on MNIST and SVHN), where each layer is
decomposed into its different parts – (1) the input to that layer, (2) the
embedding after multiplication with $S_{l}V_{l}^{\top}$, and (3) the embedding
after multiplication with the left singular vectors and application of the
ReLU. We see here the the ratio of within-class to between-class variability
decreases mostly between steps (1) and (2), on application of the right
singular structure. This has profound implications for both neural collapse
and feature learning. It shows that, counter-intuitively, DNNs cancel most of
the intra-class sample differences (which, assuming neural collapse preserves
task-relevant information, is the irrelevant variation in the data) by the
linear weight matrix, and not the non-linear ReLU.

Since these results are in a regime where the NFA holds with high correlation
(see Appendix B), we know that $S_{l}V_{l}^{\top}$ will project into a very
similar space as the AGOP. Therefore, AGOP performs linear denoising that
significantly decreases the within-class variability, further supporting its
importance in feature learning.

Unlike DNC1, we have observed a strong DNC2 progression only in the very last
layer. This is due to the fact that DNC2 was observed in the related work
(Rangamani et al., 2023; Parker et al., 2023) in a regime of large
initialization. In contrast, the NFA holds in a regime of small initialization
(Radhakrishnan et al., 2022; Beaglehole et al., 2023), and in this small
initialization (or feature learning) regime, the low-rank bias prevents strong
collapse from occurring, see Li et al. (2020) for a good intuition on the low-
rank bias.

##  5 Deep Neural Collapse in Deep Recursive Feature Machines

We have given evidence that the improvement in neural collapse for DNNs is
induced by the AGOP. We now demonstrate that AGOP is a mechanism for DNC in
the Deep RFM model, which was proposed in prior work as an abstraction for
deep feature learning in DNNs via the layerwise introduction of AGOP. We first
demonstrate empirically that Deep RFM exhibits DNC - NC1 and NC2 - and that
the induction of DNC is due to the application of AGOP. We additionally
formally establish that the reduction in within-class variability cannot be
due to the application of the random feature map. Finally, we prove
theoretically in an asymptotic and non-asymptotic setting that DNC will occur
as a consequence of Deep RFM.

###  5.1 Empirical evidence

We first visualize the formation of DNC in the input vectors of across each
layer of Deep RFM. In Figure 2, we show the Gram matrix,
$\bar{X}_{l}\bar{X}_{l}^{\top}$, of centered and normalized feature vectors
$\bar{X}_{l}=\left(X_{l}-\mu_{G}^{l}\right)\odiv\|X_{l}-\mu_{G}^{l}\|$
extracted from several layers of Deep RFM, exactly as presented in Algorithm 1
and trained on $2000$ points from CIFAR-10. After 18 layers of DeepRFM, the
final Gram matrix is approximately equal to that of perfectly collapsed data -
all points are at exactly their class means, and the class means form an ETF.

We then apply Deep RFM to full CIFAR-10, and show collapse occurring as a
consequence of projection onto the square root of the AGOP, $M^{1/2}_{l}$
(Figure 3). We observe that the improvement in NC1 is entirely due to
$M^{1/2}_{l}$, and the random feature map actually worsens the NC1 value.
Similarly, the majority of improvement in NC2, measured by
$\frac{1}{K^{2}}\|\bar{H}_{l}\bar{H}_{l}^{\top}-\Sigma_{\mathrm{ETF}}\|_{F}$,
is due to multiplication with this matrix. These results are consistent with
the crucial point that Deep RFM deviates from a simple deep random feature
model by additionally projecting onto the AGOP, and we do not expect a random
feature model to induce neural collapse (as we will derive in Section 5.2).
Note that these observations are consistent across datasets (see Appendix B).

###  5.2 Random features cannot significantly reduce within-class variability

In the following proposition, we prove that the random features mapping with
ReLU influences the within- and between-class variabilities only marginally
and if anything, it makes the DNC1 metrics slightly worse.

(a)

(b)

Figure 3: Neural collapse with Deep RFM on full CIFAR-10 across 20 layers. We
plot the NC1 and NC2 metrics across feature vectors $X_{l}$,
$M_{l}^{1/2}X_{l}$, and $\sigma(W_{l}M_{l}^{1/2}X_{l})$. Plot (a) is on a log
scale, while plot (b) is with a linear scale.

######  Proposition 5.1.

Let $x,y\in\mathbb{R}^{d}$ be two fixed vectors of unit length. Let
$W\in\mathbb{R}^{D\times d}$ be a weight matrix whose entries are initialized
i.i.d. from $\mathcal{N}(0,2/D).$ Assume $x^{\top}y=r.$ Then,

$\displaystyle\mathbb{E}(\sigma(Wx)^{\top}\sigma(Wy))=$
$\displaystyle\frac{(1-r^{2})^{3/2}}{\pi}+\frac{r^{2}\sqrt{1-r^{2}}}{\pi}$ (1)
$\displaystyle+$
$\displaystyle\frac{r}{2}+\frac{r}{\pi}\arctan\left(\frac{r}{\sqrt{1-r^{2}}}\right).$

Moreover, the variance of this dot product scales with $1/D,$ and the product
is sub-exponential, making it well-concentrated.

The proof is a direct computation, which is deferred to Appendix A.

This proposition gives a relationship between the dot product $x^{\top}y$ of
two vectors and the dot product of the outputs of their corresponding random
feature maps $\sigma(Wx)^{\top}\sigma(Wy).$ This is relevant because for unit
vectors, the dot product is a direct indication of the distance between the
two vectors. In Figure 4 (left), we plot a function which takes the distance
between two unit vectors $x,y$ and outputs the expected distance between their
feature maps $\sigma(Wx),\sigma(Wy).$ While the distance in the input space
ranges between $0$ (aligned vectors) and $2$ (counter-aligned vectors), in the
output space the maximum is $\sqrt{2}$ due to the ReLU activation. However,
importantly, we see that the expected distance between $\sigma(Wx)$ and
$\sigma(Wy)$ is bigger than the corresponding linear interpolation, and the
ratio

$\frac{\mathbb{E}\left[\left\|\sigma(Wx)-\sigma(Wy)\right\|\right]}{\left\|x-y\right\|},$

is monotonically decreasing and lies in an interval $[1/\sqrt{2},1],$ as
indicated in the Figure 4 (right), see also (Liu & Hui, 2023) for similar
considerations.

This means that the distances between data points generally increase, but not
drastically. They increase irrespective of the angle between data points,
however they tend to expand relatively more, if the data points are already
aligned. Therefore, the DNC1 metric should on average marginally increase or
stay constant. This is in agreement with our Deep RFM measurements on all
vision datasets (Figure 3 and Appendix B). Importantly, this also supports the
connection between neural networks and Deep RFM, since in DNNs the DNC1 metric
also doesn’t significantly change after applying the ReLU (Figure 1).

The key phenomenon in the DNC formation is thus the feature learning through
the AGOP (explicitly in Deep RFM), which is aligned with the right singular
space and singular values in DNNs due to the NFA.

Figure 4: The relationship between $\left\|x-y\right\|$ and
$\mathbb{E}(\left\|\sigma(Wx)-\sigma(Wy)\right\|)$. On the left, we see the
direct relationship between these two quantities compared with the linear
interpolation (green). The theoretical computation (orange) and the simulation
result (blue) match almost perfectly.

###  5.3 Asymptotic analysis

Many works have derived that, under mild conditions on the data and kernel,
non-linear kernel matrices of high dimensional data are well approximated by
linear kernels with a non-linear perturbation (Karoui, 2010; Adlam et al.,
2019; Hu & Lu, 2022). In this section, we provide a proof that Deep RFM will
induce DNC when the non-linear kernels $\\{k_{l}\\}_{l}$ and feature maps
$\\{\Phi_{l}\\}_{l}$ satisfy this property. In particular, in the high-
dimensional setting in these works, non-linear kernel matrices written
$k(X):=k(X,X)$ for simplicity, are well approximated by

$\displaystyle
k(X)\approx\gamma_{k}\mathbf{1}\mathbf{1}^{\top}+X^{\top}X+\lambda_{k}I,$

where $\lambda_{k},\gamma_{k}>0$. Following (Adlam & Pennington, 2020), for
additional simplicity we consider centered kernels,

$\displaystyle k_{\mathrm{lin}}(X)\equiv X^{\top}X+\lambda_{k}I~{}.$

We analyze the Deep RFM with two slight modifications: _(i)_ we don’t
normalize the data at each iteration, and _(ii)_ we scale the data by a value
$\kappa^{-1}$ at each iteration. The recursive procedure begins with a dataset
$X_{l}$ at layer $l$ in Deep RFM, and constructs $\tilde{X}_{l}$ from the AGOP
$M_{l}$ of the kernel ridgeless regression solution,

$\displaystyle\tilde{X}_{l+1}^{\top}\tilde{X}_{l+1}=\kappa^{-1}X_{l}^{\top}M_{l}X_{l},$

for $\kappa=1-2\lambda_{k}(1+\lambda_{\Phi}^{-1}).$ Then, we will apply a non-
linear feature map $\Phi_{\mathrm{lin}}$ corresponding to a kernel $k_{\Phi}$
that exhibits a similar linearization on the data. For $k_{\Phi}$, we
introduce a new parameter $\lambda_{\Phi}$ capturing the deviation of this
feature map from linearity, and write the corresponding Gram matrix

$\displaystyle X_{l+1}^{\top}X_{l+1}$
$\displaystyle=\Phi_{\mathrm{lin}}(\tilde{X}_{l+1})^{\top}\Phi_{\mathrm{lin}}(\tilde{X}_{l+1})$
$\displaystyle\equiv k_{\Phi}(\tilde{X}_{l+1})$
$\displaystyle=\tilde{X}_{l+1}^{\top}\tilde{X}_{l+1}+\lambda_{\Phi}I,$

where $\lambda_{k},\lambda_{\Phi}$ are finite and
$\lambda_{\Phi}\gg\lambda_{k}>0$.

Intuitively, the ratio $\lambda_{k}/\lambda_{\Phi}$ determines how close to
linear $k_{\mathrm{lin}}$ is when applied to the dataset, and therefore the
rate of NC. By Assumption 5.3, the scaled gram matrix

$\bar{G}_{l}\equiv\lambda_{\Phi}^{-1}X^{\top}_{l}X_{l}$

has minimum eigenvalue at least $1$. Meanwhile, the scaled kernel matrix has
form,

$\lambda_{\Phi}^{-1}\left(X^{\top}_{l}X_{l}+\lambda_{k}I\right)=\bar{G}_{l}+\lambda_{k}/\lambda_{\Phi}\cdot
I~{},$

If $\lambda_{k}/\lambda_{\Phi}$ were exactly $0$, then
$\left(\bar{G}_{l}+\lambda_{k}/\lambda_{\Phi}\cdot I\right)^{-1}$ is equal to
$\bar{G}^{-1}_{l}=\lambda_{\Phi}\left(X^{\top}_{l}X_{l}\right)^{-1}$, the
scaled linear regression solution. In this case, collapse will occur in just
one layer of Deep RFM, as the AGOP of a linear function is exactly the outer
product of the solution with itself. Otherwise, if this ratio were large, i.e.
$\lambda_{k}/\lambda_{\Phi}=\Omega(1)$,
$\left(\bar{G}_{l}+\lambda_{k}/\lambda_{\Phi}\cdot I\right)^{-1}$ can be far
from the linear regression solution. In this case, the kernel regression
solution contains non-linear variation that increases the rank of the AGOP
(preventing collapse). Between these limits, a larger $\lambda_{k}$ delays
collapse as a function of the number of layers of Deep RFM, while larger
$\lambda_{\Phi}$ accelerates collapse. Note that this theory offers an
explanation why a non-linear activation is needed for DNC: $\lambda_{\Phi}=0$
with a linear feature map.

We provide a simple interpretation for why non-linear feature maps introduce
an identity term $\lambda_{\Phi}I$ to the data Gram matrix. As in our
experiments on Deep RFM where $\Phi(x)=\sigma(Wx)$, the random feature map
separates all pairs of data points, especially nearby ones, lowering the off-
diagonal entries in the empirical kernel and increasing the effective rank of
this matrix (Figure 4).

Having established our simplified Deep RFM, we prove that DNC occurs in the
following full-rank, high-dimensional setting.

######  Assumption 5.2 (Data is high dimensional).

We assume that the data has dimension $d\geq n$.

######  Assumption 5.3 (Data is full rank).

We assume that the initial data Gram matrix $X_{1}^{\top}X_{1}$ has minimum
eigenvalue at least $\lambda_{\phi}>0$. Otherwise the theorem holds for $l>2$,
after the non-linear feature map is applied.

Our analysis requires these assumptions to guarantee
$(X_{l}^{\top}X_{l})^{-1}$ is stable (i.e. it exists and has maximum
eigenvalue $\lambda_{\Phi}^{-1}$ at every iteration) for all $l$. Our setup
can be viewed as a “beyond-worst-case” theoretical analysis (Roughgarden,
2019), where inputs are taken to be arbitrary followed by perturbation with a
“smoothing” parameter $\lambda_{\Phi}$ (Spielman & Teng, 2004). Our empirical
evidence for Deep RFM on real datasets, where $n>d$, suggests that these
assumptions can be relaxed using a more refined analysis.

We present our main theorem.

######  Theorem 5.4 (Deep RFM exhibits neural collapse).

Suppose we apply Deep RFM on any dataset with $\\{\Phi_{l}\\}_{l}$ and
$\\{k_{l}\\}_{l}$ all chosen as the feature map $\Phi_{\mathrm{lin}}$ and
kernel $k_{\mathrm{lin}}$ above, with no ridge parameter ($\mu=0$). Then,
there exists a universal constant $C>0$, such that for any $0<\epsilon\leq 1$,
provided
$\lambda_{k}\leq\frac{C\lambda_{\Phi}}{2(1+\lambda_{\Phi}^{-1})n}(1-\epsilon)$,
and for all layers $l\in\\{2,\ldots,L\\}$ of Deep RFM,

$\displaystyle\|\tilde{X}_{l+1}^{\top}\tilde{X}_{l+1}-yy^{\top}\|\leq(1-\epsilon)\|\tilde{X}_{l}^{\top}\tilde{X}_{l}-yy^{\top}\|+O(\lambda_{k}^{2}\lambda_{\Phi}^{-2}).$

Note that this theorem immediately implies exponential convergence to NC1 and
NC2 up to error $O(L\lambda_{k}^{2}\lambda_{\Phi}^{-2})$, a small value
determined by the parameters of the problem. As a validation of our theory, we
see that this exponential improvement in the DNC metric occurs in all layers
(Figure 3).

The proof for this theorem is roughly as follows. Recall that, by our argument
earlier in this section, a linear kernel will cause collapse within just one
layer of Deep RFM. However, in the more general case we consider, a small non-
linear deviation from a linear kernel, $\lambda_{k}$, is introduced. Beginning
with the first layer, partial collapse occurs if the ratio
$\lambda_{k}\lambda_{\Phi}^{-1}$ is sufficiently small. Following the partial
collapse, in subsequent layers, the data Gram matrix will sufficiently
resemble the collapsed matrix, so that the non-linear kernel solution on the
collapsed data will behave like the linear solution, leading to further
convergence to the collapsed data Gram matrix, a fixed point of the Deep RFM
iteration (Lemma A.2).

###  5.4 Non-asymptotic analysis

Next, we show that the formation of the neural collapse is not only implicitly
given by the specific optimization procedure of the Deep RFM, but is also
implicitly regularized for in the parametrized kernel ridge regression, a
model class that includes RFM (i.e., a single layer of Deep RFM). Consider a
parametrized kernel on $\mathbb{R}\times\mathbb{R}$ of the form
$k_{M}(x,y)=\phi(\left\|x-y\right\|^{2}_{M}),$ where
$\left\|x-y\right\|^{2}_{M}=(x-y)^{\top}M(x-y)$, $\phi$ is a strictly
decreasing, strictly positive univariate function s.t. $\phi(0)=1$, and
$k_{M}$ is a valid kernel for any symmetric positive semi-definite $M$. Let
$\mathcal{H}_{M}$ be the corresponding reproducing kernel Hilbert space
(RKHS). Then, by the parametrized kernel ridge regression, we mean the
following optimization problem from Chen et al. (2023) over all possible $M$
and $\alpha\in\mathbb{R}^{n\times K}$:

$\displaystyle\underset{f,M}{\min}\frac{1}{2}\left\|f(X)-Y\right\|^{2}+\frac{\mu}{2}\left\|f\right\|^{2}_{\mathcal{H}_{M}}.$
(2)

Using the representer theorem (Schölkopf et al., 2001), the optimal solution
of Problem 2 can be written as

$f^{*}(z)=\sum_{c,i=1,1}^{K,n}\alpha_{ci}^{*}k_{M^{*}}(x_{ci},z).$

We can re-formulate Problem 2 as the following finite-dimensional optimization
problem:

$\displaystyle\underset{\alpha,M}{\min}\,\alpha^{\top}(k_{M}^{2}+\mu
k_{M})\alpha-2\alpha^{\top}k_{M}y,$ (3)

where $k_{M}:=k_{M}(X,X)$ is the matrix of pair-wise kernel evaluations on the
data and $y$ are labels organized into a column. Since we do not explicitly
regularize $M,$ we will drop the dependence on it, treat $k$ as a free
optimization variable and compute the optimal value of the following relaxed
problem:

$\displaystyle\underset{\alpha,k}{\min}\,\alpha^{\top}(k^{2}+\mu
k)\alpha-2\alpha^{\top}ky.$ (4)

We make two final simple assumptions. First, we assume that the kernel $k$ is
positive definite. This assumption is only to allow for simpler mathematical
analysis and could be dropped using continuity arguments. Note that even a
low-rank $M$ does not naturally produce singular $k$, since $k$ always has
$1$’s on the diagonal. Second, we assume that the task is binary
classification with one-dimensional labels equal to $\pm 1.$ This assumption
is only didactic and can be dropped. Now, we are ready to state the optimality
result.

######  Theorem 5.5.

The optimal limiting solution of the Problem 4 is
$k^{*}=I_{2}\otimes(\mathbf{1}_{n}\mathbf{1}_{n}^{\top})=\frac{1}{2}(yy^{\top}+\mathbf{1}_{2n}\mathbf{1}_{2n}^{\top}).$
Any positive definite $k$ sufficiently close to this (positive semi-definite)
solution is therefore nearly optimal.

###### Proof.

We provide a sketch here and defer the detailed proof to Appendix A. Assuming
positive definiteness of $k,$ the conditionally optimal $\alpha^{*}|k$ is
explicitly computable. Plugging back into the objective and simplifying, we
get:

$\underset{k}{\max}\hskip 2.84526pty^{\top}(k+\mu I)^{-1}ky.$

Let $VSV^{\top}$ be the eigenvalue decomposition of $k.$ Then, we can
equivalently rewrite this objective to:

$\sum_{i=1}^{2n}(y^{\top}v_{i})^{2}\frac{\lambda_{i}}{\lambda_{i}+\mu}.$ (5)

By continuity of the Problem 4 in both variables, we can without loss of
generality plug $k^{*}$ into this, despite it being low-rank (for details see
Appendix A). The function value in this case would equal $2n\frac{n}{n+\mu}.$
We show that any symmetric p.s.d. kernel $k$ cannot achieve better objective
value.

Note that $\sum_{i=1}^{2n}(y^{\top}v_{i})^{2}=2n$, since $(v_{i})_{i=1}^{2n}$
forms an orthogonal basis and, therefore, a tight frame with constant $1.$
Moreover, $\sum_{i=1}^{2n}\lambda_{i}=2n,$ since the diagonal entries of $k$
are by definition all equal to 1. After normalizing $y$ to unit norm, it can
be shown that we have a weighted average of $\lambda_{i}/(\lambda_{i}+\mu),$
where the weights correspond to the the alignment between $y$ and the
corresponding eigenvector $v_{i}$. We want $y$ to be most aligned with the
eigenvectors corresponding to the biggest eigenvalues. Note that this is the
case for $k^{*},$ where $y^{\top}v_{i}=0$ for all $i\geq 3,$ thus all
eigenvectors of zero eigenvalue and $(y^{\top}v_{i})^{2}=n$ for the leading
two eigenvalues that are both equal $n.$

Increasing the value of (5) with a different kernel matrix $k$ would require
having at least one of its eigenvalues be bigger than $n$ (the leading one)
and the corresponding alignment $(y^{\top}v_{1})^{2}\geq n$ (see Appendix A).
However, this is impossible due to the Perron-Frobenius theorem, which states
that a positive matrix has a simple leading eigenvalue with strictly positive
eigenvector (see Appendix A for computations). ∎

This result itself, though very general, requires further reasoning to explain
the DNC formation through Deep RFM. First, it is not hard to see that there
exist matrices $M$ that cause $k_{M}$ to resemble $k^{*}$, if the data is
already sufficiently clustered (e.g. after a few layers of Deep RFM or a pre-
trained embedding). For example, if $\left\|x_{c,i}-x_{c,j}\right\|^{2}\leq
A^{2},\hskip 2.84526pt\forall c\in[K];i,j\in[n]$ and
$\left\|x_{c,i}-x_{c^{\prime},j}\right\|^{2}\geq R^{3}A^{2},\hskip
2.84526pt\forall c\neq c^{\prime}\in[K];i,j\in[n],$ with $A>0$, $R>1$, then
just taking $M=\frac{1}{AR}I$ yields data for which the distances within
cluster do not exceed $R^{-2}$ while distances between clusters exceed $R.$
Then, $k_{M}(x_{c,i},x_{c,j})\geq\exp(-R^{-2})\overset{\sim}{=}1-R^{-2}$ and
$k_{M}(x_{c,i},x_{c^{\prime},j})\leq\exp(-R).$ For sufficiently large $R$ this
would yield the ratio of between- and within-class square distances rising
from at least $R^{3}$ to at least $R^{4}.$

Second, we argue that RFM naturally minimizes both Problem 3 and Problem 4.
This claim is natural since the RFM is a parametrized kernel regression model,
where $M$ is set to be the AGOP, which should minimize (3) as AGOP captures
task-specific low-rank structure (Chen et al., 2023). Then, given any setting
of $M$, RFM is conditionally optimal w.r.t. the parameter $\alpha$, as this
method simply solves the original kernel regression problem for each fixed
$M$. Further, we note that choosing $k$ to be a dot product kernel and $M$ to
be the AGOP of an interpolating linear classifier implies that $k_{M}=k^{*}$,
the optimal solution outlined in Theorem 5.5, by a similar argument to the
proof of Theorem 5.4.

##  6 Conclusion

This work is the first to explain the apparent yet never properly understood
connection between feature learning and (deep) neural collapse. First, we
experimentally show that the DNC1 metric progression through the layers is
mostly due to the linear denoising via the application of the center-right
singular structure, as opposed to the application of the ReLU as suggested by
the feature learning free UFM (Tirer & Bruna, 2022). Under the neural feature
ansatz, this highlights the importance of AGOP in the formation of the DNC.

We demonstrate the importance of AGOP in DNC by demonstrating this object
induces neural collapse in randomly initialized neural networks through Deep
RFM. We validate that AGOP induces NC in Deep RFM both empirically and via
asymptotic and non-asymptotic theoretical analyses. The connection is
theoretically reinforced by showing that AGOP (as opposed to the random
features) is indeed the only part of the Deep RFMs that can cause collapse.
Our findings bridge the unsatisfactory gap between the collapse and the data –
with previous work mostly only focusing on the very end of the network and
ignoring the training data and its structures.

## Acknowledgements

We acknowledge support from the National Science Foundation (NSF) and the
Simons Foundation for the Collaboration on the Theoretical Foundations of Deep
Learning through awards DMS-2031883 and #814639 as well as the TILOS institute
(NSF CCF-2112665). This work used the programs (1) XSEDE (Extreme science and
engineering discovery environment) which is supported by NSF grant numbers
ACI-1548562, and (2) ACCESS (Advanced cyberinfrastructure coordination
ecosystem: services & support) which is supported by NSF grants numbers
#2138259, #2138286, #2138307, #2137603, and #2138296. Specifically, we used
the resources from SDSC Expanse GPU compute nodes, and NCSA Delta system, via
allocations TG-CIS220009. Marco Mondelli is supported by the 2019 Lopez-Loreta
prize.

## References

  * Adlam & Pennington (2020) Adlam, B. and Pennington, J.  The neural tangent kernel in high dimensions: Triple descent and a multi-scale theory of generalization.  In _International Conference on Machine Learning_ , pp. 74–84. PMLR, 2020.
  * Adlam et al. (2019) Adlam, B., Levinson, J., and Pennington, J.  A random matrix perspective on mixtures of nonlinearities for deep learning.  _arXiv preprint arXiv:1912.00827_ , 2019.
  * Beaglehole et al. (2023) Beaglehole, D., Radhakrishnan, A., Pandit, P., and Belkin, M.  Mechanism of feature learning in convolutional neural networks.  _arXiv preprint arXiv:2309.00570_ , 2023.
  * Beaglehole et al. (2024) Beaglehole, D., Mitliagkas, I., and Agarwala, A.  Gradient descent induces alignment between weights and the empirical ntk for deep non-linear networks.  _arXiv preprint arXiv:2402.05271_ , 2024.
  * Chen et al. (2023) Chen, Y., Li, Y., Liu, K., and Ruan, F.  Kernel learning in ridge regression” automatically” yields exact low rank solution.  _arXiv preprint arXiv:2310.11736_ , 2023.
  * Dang et al. (2023) Dang, H., Nguyen, T., Tran, T., Tran, H., and Ho, N.  Neural collapse in deep linear network: From balanced to imbalanced data.  _arXiv preprint arXiv:2301.00437_ , 2023.
  * Fang et al. (2021) Fang, C., He, H., Long, Q., and Su, W. J.  Exploring deep neural networks via layer-peeled model: Minority collapse in imbalanced training.  In _Proceedings of the National Academy of Sciences (PNAS)_ , volume 118, 2021.
  * Galanti et al. (2022) Galanti, T., György, A., and Hutter, M.  Improved generalization bounds for transfer learning via neural collapse.  In _First Workshop on Pre-training: Perspectives, Pitfalls, and Paths Forward at ICML_ , 2022.
  * Haas et al. (2022) Haas, J., Yolland, W., and Rabus, B. T.  Linking neural collapse and l2 normalization with improved out-of-distribution detection in deep neural networks.  _Transactions on Machine Learning Research (TMLR)_ , 2022.
  * Han et al. (2022) Han, X. Y., Papyan, V., and Donoho, D. L.  Neural collapse under mse loss: Proximity to and dynamics on the central path.  In _International Conference on Learning Representations (ICLR)_ , 2022.
  * He & Su (2022) He, H. and Su, W. J.  A law of data separation in deep learning.  _arXiv preprint arXiv:2210.17020_ , 2022.
  * Hong & Ling (2023) Hong, W. and Ling, S.  Neural collapse for unconstrained feature model under cross-entropy loss with imbalanced data.  _arXiv preprint arXiv:2309.09725_ , 2023.
  * Hu & Lu (2022) Hu, H. and Lu, Y. M.  Universality laws for high-dimensional learning with random features.  _IEEE Transactions on Information Theory_ , 69(3):1932–1964, 2022.
  * Hui et al. (2022) Hui, L., Belkin, M., and Nakkiran, P.  Limitations of neural collapse for understanding generalization in deep learning.  _arXiv preprint arXiv:2202.08384_ , 2022.
  * Jacot et al. (2018) Jacot, A., Gabriel, F., and Hongler, C.  Neural tangent kernel: Convergence and generalization in neural networks.  In _Conference on Neural Information Processing Systems (NeurIPS)_ , 2018.
  * Ji et al. (2022) Ji, W., Lu, Y., Zhang, Y., Deng, Z., and Su, W. J.  An unconstrained layer-peeled perspective on neural collapse.  In _International Conference on Learning Representations (ICLR)_ , 2022.
  * Karoui (2010) Karoui, N. E.  The spectrum of kernel random matrices.  _The Annals of Statistics_ , pp. 1–50, 2010.
  * Kothapalli (2023) Kothapalli, V.  Neural collapse: A review on modelling principles and generalization.  In _Transactions on Machine Learning Research (TMLR)_ , 2023.
  * Kunin et al. (2022) Kunin, D., Yamamura, A., Ma, C., and Ganguli, S.  The asymmetric maximum margin bias of quasi-homogeneous neural networks.  _arXiv preprint arXiv:2210.03820_ , 2022.
  * Li et al. (2020) Li, Z., Luo, Y., and Lyu, K.  Towards resolving the implicit bias of gradient descent for matrix factorization: Greedy low-rank learning.  _arXiv preprint arXiv:2012.09839_ , 2020.
  * Liu & Hui (2023) Liu, C. and Hui, L.  Relu soothes the ntk condition number and accelerates optimization for wide neural networks.  _arXiv preprint arXiv:2305.08813_ , 2023.
  * Lu & Steinerberger (2022) Lu, J. and Steinerberger, S.  Neural collapse under cross-entropy loss.  In _Applied and Computational Harmonic Analysis_ , volume 59, 2022.
  * Mixon et al. (2020) Mixon, D. G., Parshall, H., and Pi, J.  Neural collapse with unconstrained features.  _arXiv preprint arXiv:2011.11619_ , 2020.
  * Papyan et al. (2020) Papyan, V., Han, X. Y., and Donoho, D. L.  Prevalence of neural collapse during the terminal phase of deep learning training.  In _Proceedings of the National Academy of Sciences (PNAS)_ , volume 117, 2020.
  * Parker et al. (2023) Parker, L., Onal, E., Stengel, A., and Intrater, J.  Neural collapse in the intermediate hidden layers of classification neural networks.  _arXiv preprint arXiv:2308.02760_ , 2023.
  * Poggio & Liao (2020) Poggio, T. and Liao, Q.  Explicit regularization and implicit bias in deep network classifiers trained with the square loss.  _arXiv preprint arXiv:2101.00072_ , 2020.
  * Radhakrishnan et al. (2022) Radhakrishnan, A., Beaglehole, D., Pandit, P., and Belkin, M.  Feature learning in neural networks and kernel machines that recursively learn features.  _arXiv preprint arXiv:2212.13881_ , 2022.
  * Radhakrishnan et al. (2024) Radhakrishnan, A., Belkin, M., and Drusvyatskiy, D.  Linear recursive feature machines provably recover low-rank matrices.  _arXiv preprint arXiv:2401.04553_ , 2024.
  * Rangamani et al. (2023) Rangamani, A., Lindegaard, M., Galanti, T., and Poggio, T.  Feature learning in deep classifiers through intermediate neural collapse.  _Technical Report_ , 2023.
  * Roughgarden (2019) Roughgarden, T.  Beyond worst-case analysis.  _Communications of the ACM_ , 62(3):88–96, 2019.
  * Schölkopf et al. (2001) Schölkopf, B., Herbrich, R., and Smola, A. J.  A generalized representer theorem.  In _International conference on computational learning theory_ , pp. 416–426. Springer, 2001.
  * Seleznova et al. (2023) Seleznova, M., Weitzner, D., Giryes, R., Kutyniok, G., and Chou, H.-H.  Neural (tangent kernel) collapse.  _arXiv preprint arXiv:2305.16427_ , 2023.
  * Spielman & Teng (2004) Spielman, D. A. and Teng, S.-H.  Smoothed analysis of algorithms: Why the simplex algorithm usually takes polynomial time.  _Journal of the ACM (JACM)_ , 51(3):385–463, 2004.
  * Su et al. (2023) Su, J., Zhang, Y. S., Tsilivis, N., and Kempe, J.  On the robustness of neural collapse and the neural collapse of robustness.  _arXiv preprint arXiv:2311.07444_ , 2023.
  * Súkeník et al. (2023) Súkeník, P., Mondelli, M., and Lampert, C.  Deep neural collapse is provably optimal for the deep unconstrained features model.  _arXiv preprint arXiv:2305.13165_ , 2023.
  * Thrampoulidis et al. (2022) Thrampoulidis, C., Kini, G. R., Vakilian, V., and Behnia, T.  Imbalance trouble: Revisiting neural-collapse geometry.  In _Conference on Neural Information Processing Systems (NeurIPS)_ , 2022.
  * Tirer & Bruna (2022) Tirer, T. and Bruna, J.  Extended unconstrained features model for exploring deep neural collapse.  In _International Conference on Machine Learning (ICML)_ , 2022.
  * Tirer et al. (2022) Tirer, T., Huang, H., and Niles-Weed, J.  Perturbation analysis of neural collapse.  _arXiv preprint arXiv:2210.16658_ , 2022.
  * Trivedi et al. (2014) Trivedi, S., Wang, J., Kpotufe, S., and Shakhnarovich, G.  A consistent estimator of the expected gradient outerproduct.  In _UAI_ , pp. 819–828, 2014.
  * Wang et al. (2022) Wang, P., Liu, H., Yaras, C., Balzano, L., and Qu, Q.  Linear convergence analysis of neural collapse with unconstrained features.  In _OPT 2022: Optimization for Machine Learning (NeurIPS 2022 Workshop)_ , 2022.
  * Wang et al. (2023) Wang, Z., Luo, Y., Zheng, L., Huang, Z., and Baktashmotlagh, M.  How far pre-trained models are from neural collapse on the target dataset informs their transferability.  In _Proceedings of the IEEE/CVF International Conference on Computer Vision_ , pp. 5549–5558, 2023.
  * Weinan & Wojtowytsch (2022) Weinan, E. and Wojtowytsch, S.  On the emergence of simplex symmetry in the final and penultimate layers of neural network classifiers.  In _Mathematical and Scientific Machine Learning_ , 2022.
  * Woodbury (1950) Woodbury, M. A.  _Inverting modified matrices_.  Department of Statistics, Princeton University, 1950.
  * Xu et al. (2023) Xu, M., Rangamani, A., Liao, Q., Galanti, T., and Poggio, T.  Dynamics in deep classifiers trained with the square loss: Normalization, low rank, neural collapse, and generalization bounds.  In _Research_ , volume 6, 2023.
  * Yuan et al. (2023) Yuan, G., Xu, M., Kpotufe, S., and Hsu, D.  Efficient estimation of the central mean subspace via smoothed gradient outer products.  _arXiv preprint arXiv:2312.15469_ , 2023.
  * Zhou et al. (2022) Zhou, J., Li, X., Ding, T., You, C., Qu, Q., and Zhu, Z.  On the optimization landscape of neural collapse under mse loss: Global optimality with unconstrained features.  In _International Conference on Machine Learning (ICML)_ , 2022.
  * Zhu et al. (2023) Zhu, L., Liu, C., Radhakrishnan, A., and Belkin, M.  Catapults in sgd: spikes in the training loss and their impact on generalization through feature learning.  _arXiv preprint arXiv:2306.04815_ , 2023.

##  Appendix A Proofs

See 5.1

###### Proof.

It can easily be seen that the separate indices $(Wx,Wy)_{i=1}^{D}$ are
i.i.d., therefore it suffices to analyze the
$\mathbb{E}(\sigma(w^{T}x)\sigma(w^{T}y))$ for a single row of $W.$ For
simplicity of computation, we will re-scale the $w$ to have variance 1 entries
and later re-scale back again. From the spherical symmetry of the multivariate
standard normal distribution, we can without loss of generality assume that
$x=e_{1}=(1,0,\dots,0)$ and
$y=re_{1}+\sqrt{1-r^{2}}e_{2}=(r,\sqrt{1-r^{2}},0,\dots,0).$ With this, the
expectation reduces to
$\mathbb{E}(\sigma(w_{1})\sigma(rw_{1}+\sqrt{1-r^{2}}w_{2})).$ We can further
compute (while renaming $x\equiv w_{1};y\equiv w_{2}$):

$\displaystyle\mathbb{E}(\sigma(w_{1})\sigma(rw_{1}+\sqrt{1-r^{2}}w_{2}))=\int_{\mathbb{R}}\int_{\mathbb{R}}\sigma(x)\sigma(rx+\sqrt{1-r^{2}}y)\frac{1}{2\pi}e^{-\frac{x^{2}}{2}}e^{-\frac{y^{2}}{2}}dxdy=$
$\displaystyle\sqrt{1-r^{2}}\int_{\mathbb{R}}\sigma(x)\frac{1}{\sqrt{2\pi}}e^{-\frac{x^{2}}{2}}\left(\int_{\mathbb{R}}\sigma\left(\frac{r}{\sqrt{1-r^{2}}}x+y\right)\frac{1}{\sqrt{2\pi}}e^{-\frac{y^{2}}{2}}dy\right)dx$

For the inner-integral, denote $a=\frac{r}{\sqrt{1-r^{2}}}x$ for now. Then the
inner integral can be computed as:

$\int_{-a}^{\infty}(y+a)\frac{1}{\sqrt{2\pi}}e^{-\frac{y^{2}}{2}}dy=\int_{-a}^{\infty}y\frac{1}{\sqrt{2\pi}}e^{-\frac{y^{2}}{2}}dy+a\int_{-a}^{\infty}\frac{1}{\sqrt{2\pi}}e^{-\frac{y^{2}}{2}}dy=\frac{1}{\sqrt{2\pi}}e^{-\frac{a^{2}}{2}}+a\Phi(a),$

where $\phi(\cdot),\Phi(\cdot)$ will denote the pdf and cdf of the standard
normal, respectively. Plugging back into the outer integral we get:

$\displaystyle(\cdots)=\sqrt{1-r^{2}}\int_{0}^{\infty}\frac{x}{2\pi}e^{-\frac{x^{2}}{2(1-r^{2})}}dx+r\int_{0}^{\infty}\frac{x^{2}}{\sqrt{2\pi}}e^{-\frac{x^{2}}{2}}\Phi\left(\frac{rx}{\sqrt{1-r^{2}}}\right)dx.$

We will compute the two summands separately:

$\int_{0}^{\infty}\frac{x}{2\pi}e^{-\frac{x^{2}}{2(1-r^{2})}}dx=\frac{\sqrt{1-r^{2}}}{\sqrt{2\pi}}\int_{0}^{\infty}\frac{x}{\sqrt{1-r^{2}}\sqrt{2\pi}}e^{-\frac{x^{2}}{2(1-r^{2})}}dx=\frac{1-r^{2}}{2\pi},$

where the last equality is a well-known formula for the expectation of ReLU of
a Gaussian. For the second integral, denote $a=\frac{r}{\sqrt{1-r}}$ and,
using the integration by parts for $x\phi(x)$ and $x\Phi(ax)$ we can write

$\displaystyle\int_{0}^{\infty}x^{2}\phi(x)\Phi(ax)dx$
$\displaystyle=\left[-x\phi(x)\Phi(ax)\right]_{0}^{\infty}+\int_{0}^{\infty}\phi(x)\Phi(ax)dx+a\int_{0}^{\infty}x\phi(x)\phi(ax)dx$
$\displaystyle=\int_{0}^{\infty}\phi(x)\Phi(ax)dx+a\int_{0}^{\infty}x\phi(x)\phi(ax)dx.$

The second integal evaluates to the same expression as the integral in the
first summand above, but gets multiplied by $a$ and thus equals
$\frac{r\sqrt{1-r^{2}}}{2\pi}.$ For the second integral, first note that this
parametric integral is differentiable w.r.t. $a$ and satisfies the condition
on exchanging the derivative with the integral, so we will compute it using
differentiation w.r.t. $a$. The value of the integral for $a=0$ is trivially
$1/4.$ After taking the derivative, it is easy to see that the integral
evaluates to the same one we already computed twice, so denoting
$f(a)=\int_{0}^{\infty}\phi(x)\Phi(ax)$ we get

$f^{\prime}(a)=\frac{1}{2\pi(1+a^{2})}.$

We see that the $f^{\prime}$ is just half of the pdf of standard Cauchy
distribution, while $f(0)$ is half of its cdf, therefore $f$ is half of its
cdf and we can write:

$f(a)=\frac{1}{4}+\frac{1}{2\pi}\arctan(a).$

Plugging back the $a$ and putting everything together we get the final:

$\mathbb{E}(\sigma(w_{1})\sigma(rw_{1}+\sqrt{1-r^{2}}w_{2}))=\frac{(1-r^{2})^{\frac{3}{2}}}{2\pi}+\frac{r^{2}\sqrt{1-r^{2}}}{2\pi}+\frac{r}{4}+\frac{r}{2\pi}\arctan\left(\frac{r}{\sqrt{1-r^{2}}}\right).$

We see that after scaling the weights with the scaling indicated in the
statement, we get the final result. Since $\sigma(w^{T}x)^{T}\sigma(w^{T}y)$
is clearly sub-exponential, the sum of $D$ such i.i.d. variables will be as
well. ∎

###  A.1 Non-asymptotic results

See 5.5

###### Proof.

Denote $\mathcal{L}(\alpha,k)=\alpha^{T}(k^{2}+\mu k)\alpha-2\alpha^{T}ky.$ If
$k$ is positive definite, we can compute the conditionally optimal
$\alpha^{*}|k$ is explicitly by solving for the gradient equal to 0:

$\nabla_{\alpha}\mathcal{L}(\alpha,k)=2(k^{2}+\mu
k)\alpha-2ky=0\iff\alpha=(k+\mu I)^{-1}y.$

Plugging this back into the objective and simplifying, we get:

$\underset{k}{\max}\,y^{T}(k+\mu I)^{-1}ky.$

Let $k=VSV^{T}$ be the eigenvalue decomposition of $k.$ Then, using $(k+\mu
I)^{-1}k=V(S+\mu I)^{-1}V^{T}VSV^{T},$ we can equivalently rewrite this
objective to:

$\sum_{i=1}^{2n}(y^{T}v_{i})^{2}\frac{\lambda_{i}}{\lambda_{i}+\mu}.$

By continuity of the Problem 3 in both variables, we can without loss of
generality plug in $k^{*}$ into this despite being low-rank. This can be seen
by contradiction – if the $\mathcal{L}(\alpha^{*},k^{*})$ would not be equal
to the $y^{T}(k+\mu I)^{-1}ky$, we can find a converging sequence $k_{i}$ of
symmetric PD kernel matrices that converge to $k^{*}$ for which the two
objectives are equal and by continuity they converge to the objectives
evaluated at $\alpha^{*},k^{*},$ thus they must be equal as well.

Therefore, $\mathcal{L}(\alpha^{*},k^{*})$ evaluates to $2n\frac{n}{n+\mu}.$
We show that any symmetric PD kernel $k$ cannot achieve better objective
value. Note that $\sum_{i=1}^{2n}(y^{T}v_{i})^{2}=2n$ since
$(v_{i})_{i=1}^{2n}$ forms an orthogonal basis and therefore a tight frame
with constant $1.$ Moreover, $\sum_{i=1}^{2n}\lambda_{i}=2n,$ since the
diagonal entries of $k$ are by definition all equal 1. After normalizing $y$
to unit norm, looking back at Equation 5 we see that we have a weighted
average of $\lambda_{i}/(\lambda_{i}+\mu),$ where the weights correspond to
the the alignment between $y$ and $v_{i}$ \- the corresponding eigenvector. We
want $y$ to be most aligned with the eigenvectors corresponding to the biggest
eigenvalues. Note that this is the case for the $k^{*},$ where $y^{T}v_{i}=0$
for all $i\geq 3,$ thus all eigenvectors of the zero eigenvalue and
$(y^{T}v_{i})^{2}=n$ for the leading two eigenvalues that are both equal $n.$

Beating the value of 5 with another $k$ would require having at least one
eigenvalue bigger than $n$ (the leading one) and the corresponding alignment
$(y^{T}v_{1})^{2}\geq n$. This can be seen by concavity of the function
$g(x)=\frac{x}{x+\mu},$ since

$\sum_{i=1}^{2n}\omega_{i}g(\lambda_{i})\leq
g\left(\sum_{i=1}^{2n}\omega_{i}\lambda_{i}\right)\leq
g\left(\omega_{1}\lambda_{1}+(1-\omega_{1})(2n-\lambda_{1})\right).$

From the rightmost expression it is clear that $\omega_{1}$ must be at least a
half. However, this is impossible due to the Perron-Frobenius theorem, which
states that a positive matrix has a simple leading eigenvalue with strictly
positive eigenvector, therefore $(y^{T}v_{1})<n$ since $v_{1}$ is entry-wise
positive. ∎

###  A.2 Asymptotic results

######  Lemma A.1 (Woodbury Inverse Formula (Woodbury, 1950)).

$\displaystyle\left(P+UV^{\top}\right)^{-1}=P^{-1}-P^{-1}U\left(I+V^{\top}P^{-1}U\right)^{-1}V^{\top}P^{-1}.$

######  Lemma A.2 (Fixed point of collapse).

Let $A^{*}=yy^{\top}+\lambda_{\Phi}I$, the collapsed data gram matrix
following an application of the non-linear random feature map $\Phi$. Then,

$\displaystyle\left(A^{*}\right)^{-1}=\lambda_{\Phi}^{-1}I-\lambda_{\Phi}^{-1}\left(\lambda_{\Phi}+n\right)^{-1}yy^{\top}.$

###### Proof of Lemma A.2.

We write $yy^{\top}=nUU^{\top}$, where $U\in\mathbb{R}^{n\times K}$ is the
matrix of normalized eigenvectors of $yy^{\top}$. By Lemma A.1 and that
$U^{\top}U=I$,

$\displaystyle\left(A^{*}\right)^{-1}$
$\displaystyle=\left(\lambda_{\Phi}I+nUU^{\top}\right)^{-1}$
$\displaystyle=\lambda_{\Phi}^{-1}I-\frac{n}{\lambda_{\Phi}^{2}}U\left(I+n\lambda_{\Phi}^{-1}U^{\top}U\right)^{-1}U^{\top}$
$\displaystyle=\lambda_{\Phi}^{-1}I-\frac{n}{\lambda_{\Phi}^{2}}\cdot\frac{1}{1+n\lambda_{\Phi}^{-1}}UU^{\top}$
$\displaystyle=\lambda_{\Phi}^{-1}I-\frac{n}{\lambda_{\Phi}(\lambda_{\Phi}+n)}UU^{\top}$
$\displaystyle=\lambda_{\Phi}^{-1}I-\frac{\lambda_{\Phi}^{-1}}{\lambda_{\Phi}+n}yy^{\top}.$

∎

###### Proof of Proposition 5.4.

We have solved the kernel ridgeless regression problem to get coefficients
$\alpha=\left(X_{l}X_{l}^{\top}+\lambda_{k}I\right)^{-1}y$. Then, the
predictor at layer $l$ evaluated on the training data is,

$\displaystyle f_{l}(X_{l})=(X_{l}^{\top}X_{l}+\lambda_{k}I)\alpha~{}.$

As in Deep RFM, let $M$ be the $\operatorname{AGOP}$. Then,

$\displaystyle M=\sum_{i=1}^{n}\nabla f(x^{(i)})(\nabla
f(x^{(i)}))^{\top}=X_{l}\alpha\alpha^{\top}X_{l}^{\top}~{}.$

We drop the subscript $l$ for simplicity. Therefore,

$\displaystyle
X^{\top}MX=X^{\top}X\left(X^{\top}X+\lambda_{k}I\right)^{-1}yy^{\top}\left(X^{\top}X+\lambda_{k}I\right)^{-1}X^{\top}X~{}.$

Let $A=X^{\top}X$. Applying Lemma A.1,

$\displaystyle\left(A+\lambda_{k}I\right)^{-1}=A^{-1}-\lambda_{k}A^{-1}\left(I+\lambda_{k}A^{-1}\right)^{-1}A^{-1}~{}.$

Therefore,

$\displaystyle
X^{\top}MX=\left(I-\lambda_{k}\left(I+\lambda_{k}A^{-1}\right)^{-1}A^{-1}\right)yy^{\top}\left(I-\lambda_{k}\left(I+\lambda_{k}A^{-1}\right)^{-1}A^{-1}\right)~{}.$

Applying that, by assumption, $\lambda_{k}\lambda_{\Phi}^{-1}$ is a small
value $\epsilon_{\mathrm{lin}}$, and as $\lambda_{\Phi}$ is the minimum
eigenvalue of $A$, we have that
$\lambda_{k}\|A^{-1}\|\leq\lambda_{k}\lambda_{\Phi}^{-1}\triangleq\epsilon_{\mathrm{lin}}$.
Therefore,
$\left(I+\lambda_{k}A^{-1}\right)^{-1}\lambda_{k}A^{-1}=\lambda_{k}A^{-1}-\left(\lambda_{k}A^{-1}\right)^{2}+\cdots=\lambda_{k}A^{-1}+O(\epsilon_{\mathrm{lin}}^{2})$,
where $O(\epsilon_{\mathrm{lin}}^{2})$ refers to some matrix of spectral norm
at most $C\epsilon_{\mathrm{lin}}^{2}$ for some constant $C>0$. Then,

$\displaystyle
X^{\top}MX=\left(I-\lambda_{k}A^{-1}+O(\epsilon_{\mathrm{lin}}^{2})\right)yy^{\top}\left(I-\lambda_{k}A^{-1}+O(\epsilon_{\mathrm{lin}}^{2})\right)=yy^{\top}-\lambda_{k}A^{-1}yy^{\top}-\lambda_{k}yy^{\top}A^{-1}+O(\epsilon_{\mathrm{lin}}^{2})~{}.$

Let $\tilde{A}=\frac{A-A^{*}}{\|A-A^{*}\|}$, and $\epsilon_{A}=\|A-A^{*}\|$.
Then, applying Lemma A.1,

$\displaystyle
A^{-1}=\left(A^{*}+\epsilon_{A}\tilde{A}\right)^{-1}=\left(I+\epsilon_{A}\left(A^{*}\right)^{-1}\tilde{A}\right)^{-1}\left(A^{*}\right)^{-1}=\left(A^{*}\right)^{-1}-\epsilon_{A}\left(A^{*}\right)^{-1}\left(I+\epsilon_{A}\tilde{A}\left(A^{*}\right)^{-1}\right)^{-1}\tilde{A}\left(A^{*}\right)^{-1}~{}.$

Let
$\tilde{\Psi}=\left(I+\epsilon_{A}\tilde{A}\left(A^{*}\right)^{-1}\right)^{-1}$.
Assuming partial collapse has already occurred, i.e.,
$\epsilon_{A}\lambda_{\Phi}^{-1}<1/2$,

$\displaystyle\|\tilde{\Psi}\|\leq\left(1-\epsilon_{A}\lambda_{\Phi}^{-1}\right)^{-1}\leq
2~{}.$

In this case, using that
$yy^{\top}\left(A^{*}\right)^{-1}=\left(A^{*}\right)^{-1}yy^{\top}=(1+\lambda_{\Phi}^{-1})yy^{\top}$
(Lemma A.2),

$\displaystyle X^{\top}MX$
$\displaystyle=yy^{\top}-\lambda_{k}A^{-1}yy^{\top}-\lambda_{k}yy^{\top}A^{-1}+O(\epsilon_{\mathrm{lin}}^{2})$
$\displaystyle=yy^{\top}-\lambda_{k}\left(A^{*}\right)^{-1}yy^{\top}-\lambda_{k}yy^{\top}\left(A^{*}\right)^{-1}+O(\epsilon_{\mathrm{lin}}^{2})$
$\displaystyle-\lambda_{k}\epsilon_{A}\left(A^{*}\right)^{-1}\tilde{\Psi}\tilde{A}\left(A^{*}\right)^{-1}yy^{\top}-\lambda_{k}\epsilon_{A}yy^{\top}\left(A^{*}\right)^{-1}\tilde{A}\tilde{\Psi}\left(A^{*}\right)^{-1}$
$\displaystyle=yy^{\top}-2\lambda_{k}(1+\lambda_{\Phi}^{-1})yy^{\top}+O(\epsilon_{\mathrm{lin}}^{2})$
$\displaystyle-\lambda_{k}\epsilon_{A}(1+\lambda_{\Phi}^{-1})\left(A^{*}\right)^{-1}\tilde{\Psi}\tilde{A}yy^{\top}-\lambda_{k}\epsilon_{A}(1+\lambda_{\Phi}^{-1})yy^{\top}\tilde{A}\tilde{\Psi}\left(A^{*}\right)^{-1}~{}.$

Therefore, where $\kappa=1-2\lambda_{k}(1+\lambda_{\Phi}^{-1})>1/2$ by choice
of $\lambda_{k}\cdot 2\lambda_{\Phi}^{-1}(1+\lambda_{\Phi}^{-1})n<1-\epsilon$,

$\displaystyle\|X_{l+1}^{\top}X_{l+1}-A^{*}\|$
$\displaystyle=\|\kappa^{-1}X^{\top}MX-yy^{\top}\|$
$\displaystyle=\left\|O(\epsilon_{\mathrm{lin}}^{2})+\epsilon_{A}\cdot\lambda_{k}\cdot
2\kappa^{-1}(1+\lambda_{\Phi}^{-1})\left(A^{*}\right)^{-1}\tilde{\Psi}\tilde{A}yy^{\top}\right\|$
$\displaystyle\leq
O(\epsilon_{\mathrm{lin}}^{2})+\epsilon_{A}\cdot\lambda_{k}\cdot
2\lambda_{\Phi}^{-1}(1+\lambda_{\Phi}^{-1})n$
$\displaystyle<O(\epsilon_{\mathrm{lin}}^{2})+\epsilon_{A}(1-\epsilon)$
$\displaystyle=O(\epsilon_{\mathrm{lin}}^{2})+(1-\epsilon)\|X_{l}X_{l}^{\top}-A^{*}\|~{}.$

It remains to show that partial collapse,
$\epsilon_{A}\lambda_{\Phi}^{-1}<1/2$, happens in the first iteration. To
ensure this condition, recall,

$\displaystyle
X^{\top}MX=yy^{\top}-\lambda_{k}A^{-1}yy^{\top}-\lambda_{k}yy^{\top}A^{-1}+O(\epsilon_{\mathrm{lin}}^{2})~{}.$

Therefore, because $\|A^{-1}\|\leq\lambda_{\Phi}^{-1}$, where $\Psi$ is some
error matrix with norm $1$,

$\displaystyle
X^{\top}MX=yy^{\top}+\lambda_{k}\lambda_{\Phi}^{-1}n\Psi+O(\lambda_{k}^{2}\lambda_{\Phi}^{-2})~{}.$

Therefore,

$\displaystyle\|X^{\top}MX-yy^{\top}\|<1/2~{},$

provided $\lambda_{k}\lambda_{\Phi}^{-1}n<C_{1}/2$, for some universal
constant $C_{1}$. ∎

##  Appendix B Additional plots

###  B.1 Deep RFM

(a)

(b)

Figure 5: Neural collapse with Deep RFM on SVHN.

(a)

(b)

Figure 6: Neural collapse with Deep RFM on MNIST.

###  B.2 Neural networks

(a)

(b)

Figure 7: Feature variability collapse (NC1) from different singular value
decomposition components.

###  B.3 NFA metrics for neural networks

(a)

(b)

(c)

Figure 8: Train loss and NFA correlations. In the first column, we plot the
evolution of $\operatorname{tr}{\Sigma_{W}}/\operatorname{tr}{\Sigma_{T}}$ and
$\operatorname{tr}{\Sigma_{B}}/\operatorname{tr}{\Sigma_{T}}$, where
$\Sigma_{T}\triangleq\Sigma_{W}+\Sigma_{B}$. In the second column we plot the
development of the NFA, where correlation is measured between
$W_{l}^{\top}W_{l}$ and the square root of the AGOP with respect to the inputs
at layer $l$, $X_{l}$. The third column is the train loss of the neural
network.

##  Appendix C Experimental details

In both Deep RFM and neural networks, we one hot encode labels and use $\pm 1$
for within/outside of each class.

For the neural network experiments, we use 5 hidden layer MLP networks with no
biases and ReLU activation function, initialized using Xavier Uniform with
initialization scale $0.01$ in all hidden layers. We use default
initialization for the linear readout layer. The layers all use width $512$.
We train with Adam using no weight decay and learning rate $2\times 10^{-5}$
for 250 epochs and measure the AGOP and NC metrics every 10 epochs.

For the Deep RFM experiments, we use the Laplace kernel $k$, which evaluates
two datapoints $x,z\in\mathbb{R}^{d}$ as
$k(x,z)=\exp\left(-\|x-z\|_{2}/L\right)$, where $L$ is a bandwidth parameter.
In these experiments, we set $L=2.0$. We perform ridgeless regression to
interpolate the data, i.e. $\mu=0.$ We use width $512$ for the random feature
map in Deep RFM with ReLU activation function.

Unless otherwise specified, we perform experiments using the first 50,000
points from each of MNIST, SVHN, and CIFAR-10 loaded with Torchvision
datasets.