figuret

# Principal Components Bias in Over-parameterized Linear Models, and its
Manifestation in Deep Neural Networks

Guy Hacohen<EMAIL_ADDRESS>
The School of Computer Science and Engineering
Edmond and Lily Safra Center for Brain Sciences
The Hebrew University of Jerusalem
Jerusalem, Israel
Daphna Weinshall<EMAIL_ADDRESS>
The School of Computer Science and Engineering
The Hebrew University of Jerusalem
Jerusalem, Israel

###### Abstract

Recent work suggests that convolutional neural networks of different
architectures learn to classify images in the same order. To understand this
phenomenon, we revisit the over-parametrized deep linear network model. Our
asymptotic analysis, assuming that the hidden layers are wide enough, reveals
that the convergence rate of this model’s parameters is exponentially faster
along the directions of the larger principal components of the data, at a rate
governed by the corresponding singular values. We term this convergence
pattern the _Principal Components bias (PC-bias)_. Empirically, we show how
the _PC-bias_ streamlines the order of learning of both linear and non-linear
networks, more prominently at earlier stages of learning. We then compare our
results to the simplicity bias, showing that both biases can be seen
independently, and affect the order of learning in different ways. Finally, we
discuss how the _PC-bias_ may explain some benefits of early stopping and its
connection to PCA, and why deep networks converge more slowly with random
labels.

Keywords: deep linear networks, learning dynamics, PC-bias, simplicity bias,
learning order

##  1 Introduction

The dynamics of learning in deep neural networks is an intriguing subject, not
yet sufficiently understood. Recent empirical studies of learning dynamics
(Hacohen et al., 2020; Pliushch et al., 2021) showed that neural networks
learn the training examples of natural datasets in a consistent order, and
further impose a consistent order on the test set. Below we call this effect
_Learning Order Constancy_ (LOC). Currently, the characteristics of visual
data, which may explain this phenomenon, remain unclear. Surprisingly, this
universal order persists despite the variability introduced into the training
of different models and architectures.

To understand this phenomenon, we start by analyzing the deep linear network
model (Saxe et al., 2014, 2019), defined by the concatenation of linear
operators in a multi-class classification setting. Accordingly, in Section 3
we prove that the convergence of the weights of deep linear networks is
governed by the eigendecomposition of the raw data, which is blind to the
labels of the data, in a phenomenon we term _PC-bias_. These asymptotic
results are valid when the hidden layers are wide enough and can be seen as an
extension of the known behavior of the single-layer convex linear model (Le
Cun et al., 1991). In Section 4, we empirically show that this pattern of
convergence is indeed observed in deep linear networks of moderate width,
validating the plausibility of our assumptions. We continue by showing that
the _LOC-effect_ in deep linear networks is determined solely by their _PC-
bias_. We prove a similar (weaker) result for the non-linear two-layer ReLU
model trained on a binary classification problem, introduced by Allen-Zhu et
al. (2019), where this model is presented as a certain extension of NTK (Jacot
et al., 2018).

In Section 5, we extend the study empirically to non-linear networks and
investigate the relationship between the _PC-bias_ and the _LOC-effect_ in
general deep networks. We first show that the order by which examples are
learned by linear networks is highly correlated with the order induced by
prevalent deep CNN models. We then show directly that the learning order of
non-linear CNN models is affected by the principal decomposition of the data.
Moreover, the _LOC-effect_ diminishes when the data is whitened, indicating a
tight connection between the _PC-bias_ and the _LOC-effect_.

Our results are reminiscent of another phenomenon, termed _Spectral bias_ (see
Section 2.2), which associates the learning dynamics of neural networks with
the Fourier decomposition of functions in the hypothesis space. In Section 5.3
we investigate the relation between the _PC-bias_ , _spectral bias_ , and the
_LOC-effect_. We find that the _LOC-effect_ is very robust: (i) when we
neutralize the _spectral bias_ by using low complexity models such as deep
linear networks, the effect is still observed; (ii) when we neutralize the
_PC-bias_ by using whitened data, the _LOC-effect_ persists. We hypothesize
that at the beginning of learning, the learning dynamics of neural models is
governed by the eigendecomposition of the raw data. As learning proceeds,
control of the dynamics slowly shifts to other factors.

The PC-bias has implications beyond the _LOC-effect_ , as expanded in Section
6:

##### 1\. Early stopping.

It is often observed that when training deep networks with real data, the
highest generalization accuracy is obtained before convergence. Consequently,
early stopping is often prescribed to improve generalization. Following the
commonly used assumption that in natural images the lowest principal
components correspond to noise (Torralba and Oliva, 2003), our results predict
the benefits of early stopping, and relate it to PCA. In Section 6 we
investigate the relevance of this conclusion to real non-linear networks (see
Basri et al., 2019; Li et al., 2020, for complementary accounts).

##### 2\. Slower convergence with random labels.

Zhang et al. (2017) showed that neural networks can learn any label
assignment. However, training with random label assignments is known to
converge slower as compared to training with the original labels (Krueger et
al., 2017). We report a similar phenomenon when training deep linear networks.
Our analysis shows that when the principal eigenvectors are correlated with
class identity, as is often the case in natural images, the loss decreases
faster when given true label assignments as against random label assignments.
In Section 6 we investigate this hypothesis empirically in linear and non-
linear networks.

##### 3\. Weight initialization.

Different weight initialization schemes have been proposed to stabilize the
learning and minimize the hazard of "exploding gradients" (e.g., Glorot and
Bengio, 2010; He et al., 2015). Our analysis (see Appendix A) identifies a
related variant, which eliminates the hazard when all the hidden layers are
roughly of equal width. In the deep linear model, it can be proven that the
proposed normalization variant in a sense minimizes repeated gradient
amplification.

##  2 Scientific Background and Previous Work

In Section 2.1 we put our theoretical results in the context of extensive
related work. In Section 2.2 we briefly discuss work on the _simplicity bias_
, one expression of which being the _spectral_ bias alluded to above.

###  2.1 Theoretical background

A large body of work concerns the analysis of deep over-parameterized linear
networks. While not a universal approximator, this model is nevertheless
trained by minimizing a non-convex objective function with a multitude of
minima. The investigation of such networks is often employed to shed light on
the learning dynamics when complex geometric landscapes are explored by GD
(Fukumizu, 1998; Arora et al., 2018; Wu et al., 2019a; Du and Hu, 2019; Du et
al., 2019; Hu et al., 2020b; Yun et al., 2021). Note that while such networks
provably achieve a simple low-rank solution (Ji and Telgarsky, 2019; Du et
al., 2018) when given linearly separable data, in our work nothing like this
is assumed about the data.

Early on Baldi and Hornik (1989) characterized the optimization landscape of
the over-parameterized linear model and its relation to PCA. More recent work
suggests that with sufficient over-parameterization, this landscape is well
behaved (Kawaguchi, 2016; Zhou and Liang, 2018) and all its local minima are
global (Laurent and von Brecht, 2017). Deep linear networks are also used to
study biases induced by architecture or optimization (Ji and Telgarsky, 2019;
Wu et al., 2019b).

Several related studies investigated the dynamics of learning in over-
parameterized models by way of spectral analysis. However, while we employ the
eigenvectors of the data covariance matrix111$X$ and $Y$ denote the matrices
whose columns are the data points and one-hot label vectors respectively, see
notations in Section 3. $\Sigma_{\scriptscriptstyle XX}=XX^{\top}$, most of
the behavior identified in these studies is driven by the eigenvectors of the
cross-covariance matrix $\Sigma_{\scriptscriptstyle YX}=YX^{\top}$. This
difference is crucial, as $\Sigma_{\scriptscriptstyle XX}$ is blind to the
class labels, unlike $\Sigma_{\scriptscriptstyle YX}$. In addition, quite a
few of the studies reviewed below assumed a shallow 2-layer network, while our
analysis accommodates over-parametrized networks of any depth, further showing
that convergence depends on the eigenvectors of $\Sigma_{\scriptscriptstyle
XX}$ rather than $\Sigma_{\scriptscriptstyle YX}$ without any additional
assumption.

More specifically, Saxe et al. (2019) analyzed the evolution of learning
dynamics in a shallow $2$-layers linear model, showing that convergence is
guided by the eigenvectors of $\Sigma_{\scriptscriptstyle YX}$. This analysis
assumed $\Sigma_{\scriptscriptstyle XX}=I_{d}$, thus obscuring the dependence
of convergence on the raw data eigenvectors of $\Sigma_{\scriptscriptstyle
XX}$, as all the eigenvalues are now identical by assumption. Likewise, Arora
et al. (2018) also assume that $\Sigma_{\scriptscriptstyle XX}=I_{d}$ while
investigating continuous gradients of the deep linear model and allowing for a
more general loss. Although there is some superficial resemblance between
their pre-conditioning and our _gradient scale matrices_ as defined below,
they are defined quite differently and have different properties.

Similarly to Saxe et al. (2019), Gidel et al. (2019) also investigated a
shallow 2-layers linear model, but no longer assumed that
$\Sigma_{\scriptscriptstyle XX}=I_{d}$. Instead, they assumed that
$\Sigma_{\scriptscriptstyle XX}$ and $\Sigma_{\scriptscriptstyle YX}$ can be
jointly decomposed, having the same eigenvectors. Under these conditions, they
were able to show that convergence is guided by the eigenvectors of
$\Sigma_{\scriptscriptstyle YX}$, which correspond to the eigenvectors of
$\Sigma_{\scriptscriptstyle XX}$ by the joint decomposition assumption.
However, $\Sigma_{\scriptscriptstyle XX}$ in this case still depends on the
class labels by assumption, which is not realistic for many datasets.
Recently, Nguyen (2021) showed that the learning dynamics of $2$ layered non-
linear auto-encoders converge faster along the eigenvectors of
$\Sigma_{\scriptscriptstyle XX}$. While their results resemble ours, unlike us
they focus on auto-encoders and assume that the data is sampled from a
Gaussian distribution.

Relations between shallow and deep over-parameterized linear networks are also
often studied through the lens of gradient flow (Arora et al., 2018, 2019a).
Arora et al. (2019a) derived an equation for the gradient-flow of over-
parameterized deep linear networks, which was later expanded to shallow ReLU
networks by Williams et al. (2019). Bah et al. (2019) showed that this
gradient flow can be re-interpreted as a Riemannian gradient flow of matrices
of some fixed ranked, hinting at a strong connection between the optimization
process of shallow and deep linear networks, in accordance with our results.

Another line of related work was pioneered by Jacot et al. (2018), who showed
how neural networks can be analyzed using kernel theory and introduced the
Neural Tangent Kernel (NTK). In this framework, it was shown that convergence
is fastest along the largest kernel’s principal components. In the one-layer
linear model, this result implies that convergence depends on the eigenvectors
of $\Sigma_{\scriptscriptstyle XX}$, but to the best of our knowledge, it was
not extended to the over-parameterized deep model analyzed here. While similar
results to ours may be achieved in the future by using kernel theory, we
provide here direct proof, thus bypassing possible limitations of this kernel
theory (Chizat et al., 2019; Yehudai and Shamir, 2019). This direct proof
further enables us to examine the behavior outside the limit of infinite
width, and carefully examine our assumptions using empirical methods.

###  2.2 Simplicity Bias

Diverse empirical observations seem to support the hypothesis that neural
networks start by learning a simple model, which then gains complexity as
learning proceeds (Gunasekar et al., 2018; Soudry et al., 2018; Hu et al.,
2020a; Kalimeris et al., 2019; Gissin et al., 2020; Heckel and Soltanolkotabi,
2020; Ulyanov et al., 2018; Pérez et al., 2019; Cao et al., 2021; Jin and
Montúfar, 2020). This phenomenon is sometimes called _simplicity bias_ (Dingle
et al., 2018; Shah et al., 2020).

In a related line of work, Rahaman et al. (2019) empirically demonstrated that
the complexity of classifiers learned by ReLU networks increases with time.
Basri et al. (2019, 2020) showed theoretically, by way of analyzing elementary
neural network models, that these models first fit the data with low-frequency
functions, and gradually add higher frequencies to improve the fit.
Nevertheless, the _spectral bias_ and _PC-bias_ are inherently different.
Indeed, the eigendecomposition of raw images is closely related to the Fourier
analysis of images as long as the statistical properties of images are
(approximately) translation-invariant (Simoncelli and Olshausen, 2001;
Torralba and Oliva, 2003). Still, the _PC-bias_ is guided by spectral
properties of the raw data and is therefore blind to class labels. In
contrast, the _spectral bias_ , as well as the related _frequency bias_ that
has been shown to characterize NTK models (Basri et al., 2020), are all guided
by spectral properties of the learned hypothesis, which crucially depend on
label assignment.

##  3 Theoretical Analysis

In this section we analyze the over-parameterized deep linear networks model,
in which the principal component bias fully governs the learning dynamics.

###  3.1 Notations and Definitions

Let $\mathbb{X}=\\{(\bm{x}_{i},\bm{y}_{i})\\}_{i=1}^{n}$ denote the training
data, where $\bm{x}_{i}\in\mathbb{R}^{q}$ denotes the i-th data point and hot
vector $\bm{y}_{i}\in\\{0,1\\}^{K}$ denotes its corresponding label. Let
$\frac{1}{n_{i}}\bm{m}_{i}$ denote the centroid (mean) of class $i$ with
$n_{i}$ points
$\frac{1}{n_{i}}m_{i}=\frac{1}{n_{i}}\sum_{j=1}^{n}{x_{j}\mathbbm{1}_{\left[y_{j}=i\right]}}$,
and $M=[\bm{m}_{1}\ldots\bm{m}_{K}]^{\top}$ – the matrix whose rows are
vectors $\bm{m}_{i}^{\top}$. Finally, let $X$ and $Y$ denote the matrices
whose $i^{th}$ column is $\bm{x}_{i}$ and $\bm{y}_{i}$ respectively.
$\Sigma_{\scriptscriptstyle XX}=XX^{\top}$ and $\Sigma_{\scriptscriptstyle
YX}=YX^{\top}$ denote the covariance matrix of $X$ and cross-covariance of $X$
and $Y$ respectively. We note that $\Sigma_{\scriptscriptstyle XX}$ captures
the structure of the data irrespective of class identity, and that
$\Sigma_{\scriptscriptstyle YX}=M$.

####  3.1.1 Definitions

######  Definition 1 (Principal coordinate system)

The coordinate system obtained by rotating the data in $\mathbb{R}^{q}$ by an
orthonormal matrix $U^{\top}$, where $SVD(\Sigma_{\scriptscriptstyle
XX})\\!=\\!UDU^{\top}$.

In this system $\Sigma_{\scriptscriptstyle XX}\\!=\\!D$, a diagnoal matrix
whose elements are the singular values of $XX^{\top}$, arranged in decreasing
order $d_{1}\geq d_{2}\geq\ldots\geq d_{q}\geq 0$.

We analyze a deep linear network with $L$ layers, whose parameters are
computed by minimizing (using gradient descent) the following loss
$L(\mathbb{X})$

$\displaystyle
L(\mathbb{X})=\frac{1}{2}\|{\bm{W}}X-Y\|_{F}^{2}\quad\quad{\bm{W}}\coloneqq\prod_{l=L}^{1}W_{l},\quad
W_{l}\in\mathbb{R}^{m_{l}\times m_{l-1}}$ (1)

Above $m_{l}$ denotes the number of neurons in layer $l$, where $m_{0}=q$ and
$m_{L}=K$.

######  Definition 2 (Compact representation)

Given a deep linear network $\bm{y}=W_{L}\cdot\ldots\cdot W_{1}\bm{x}$, its
compact representation is denoted ${\bm{W}}\in\mathbb{R}^{K\times q}$ where
${\bm{W}}=\prod_{l=L}^{1}W_{l}$.

######  Definition 3 (Gradient matrix)

For a deep linear network whose compact representation is ${\bm{W}}$, the
gradient matrix is $G_{r}={\bm{W}}\Sigma_{\scriptscriptstyle
XX}-\Sigma_{\scriptscriptstyle YX}$.

In the principal coordinate system, $G_{r}={\bm{W}}D-M$.

######  Definition 4 (STD initialization)

Given a sequence of random matrices $\\{W_{l}\\}_{{l}=1}^{L}$ where
$W_{l}\in\mathbb{R}^{m_{l}\times m_{{l}-1}}$, which are chosen i.i.d from a
distribution with mean 0 and variance $\sigma^{2}_{l}$, define

$\sigma^{2}_{l}=\frac{2}{m_{{l}-1}+m_{l}}~{}\forall{l}\in[2\ldots{L}-1],\qquad\sigma^{2}_{1}=\frac{1}{m_{1}},\qquad\sigma^{2}_{L}=\frac{1}{m_{{L}-1}}$

Def. 4 is a variant of the Glorot initialization (Glorot and Bengio, 2010).
For more analytical and empirical results on _STD initialization_ , refer to
App. B.2 and App. C.1 respectively.

####  3.1.2 Asymptotic Analysis

Our analysis involves two asymptotic quantities: the learning rate (or
gradient step size) $\mu$, and the network’s width $\mathbbm{m}$
($\mathbbm{m}$ denotes the width of the smallest hidden layer). Accordingly,
our analysis neglects terms of magnitude $O(\mu^{2})$ and
$O\left(\frac{1}{\mathbbm{m}}\right)$, and employs convergence in probability
when $\mathbbm{m}\to\infty$. In Fig. 3 we show the plausibility of these
assumptions, as the predicted dynamics is seen to hold even when the width
$\mathbbm{m}$ is smaller than the input size, and the step size $\mu$ permits
convergence in a relatively short time.

To simplify the presentation and to allow for modular analysis, we use the
following conventions:  (i) $O(\mu^{2})$and
$O\left(\frac{1}{\mathbbm{m}}\right)$ denote terms which are upper bounded by
$c\mu^{2}$ when $\mu$ is small and $c\frac{1}{\mathbbm{m}}$ when $\mathbbm{m}$
is large respectively, where $c\in\mathbb{R}$ denotes a fixed constant that
does _not_ depend on either $\mu$ or $\mathbbm{m}$.  (ii) Depending on the
context, if matrices are involved then $O(\mu^{2})$ and
$O(\frac{1}{\mathbbm{m}})$ denote full rank matrices with element-wise
asymptotic quantities defined as above.  (iii) The notation $\xrightarrow{p}$
denotes convergence in probability when $\mathbbm{m}\to\infty$, where the
probabilistic lower bound on width $\mathbbm{m}$ does _not_ depend on $\mu$.
When matrices are involved, convergence in probability occurs element-wise.

The treatment of the two asymptotic quantities $\mu$ and $\mathbbm{m}$ is
consolidated in the formal analysis of the network’s dynamics, as detailed in
the proofs of Thms 15-16 in Appendix A.2, and the proofs of Thms 7-8 in
Appendix B.2. Importantly, the analysis is valid at each time step $t$, where
the only thing that changes with time is the magnitude of the constants
governing the asymptotic behavior.

###  3.2 The Dynamics of Deep Over-Parametrized Linear Networks

Extending the analysis described in Du and Hu (2019), we derive here the
temporal dynamics of ${\bm{W}}$, denoted ${\bm{W}}^{(t)}$, as it changes with
each GD step.

###### Proposition 5

Let $G_{r}^{(t)}$ (Def. 3) denote the gradient matrix at time ${t}$. Let
${B}_{{l}}^{(t)}$ and ${A}_{{l}}^{(t)}$ denote the gradient scale matrices,
which are defined as follows

$\begin{split}{B}_{{l}}^{(t)}\coloneqq\Big{(}\prod_{j=l}^{1}W_{j}^{(t)}\Big{)}^{\top}\Big{(}\prod_{j=l}^{1}W_{j}^{(t)}\Big{)}\in\mathbb{R}^{q\times
q}\qquad{A}_{{l}}^{(t)}\coloneqq\Big{(}\prod_{j={L}}^{{l}+1}W_{j}^{(t)}\Big{)}\Big{(}\prod_{j={L}}^{{l}+1}W_{j}^{(t)}\Big{)}^{\top}\in\mathbb{R}^{K\times
K}\end{split}$ (2)

The compact representation ${\bm{W}}^{(t)}$ obeys the following dynamics:

${\bm{W}}^{(t+1)}={\bm{W}}^{(t)}-\mu\sum_{{l}=1}^{L}{A}_{{l}}^{(t)}\cdot
G_{r}^{(t)}\cdot{B}_{{l}-1}^{(t)}+O(\mu^{2})$ (3)

The proof can be found in Appendix B.1. Note that
${B}_{0}^{(t)}={A}_{{L}}^{(t)}=I$,
${B}_{{L}}^{(t)}={{\bm{W}}^{(t)}}^{\top}{\bm{W}}^{(t)}$, and
${A}_{0}^{(t)}={\bm{W}}^{(t)}{{\bm{W}}^{(t)}}^{\top}$.

##### Gradient scale matrices.

When the number of hidden layers is $0$ ($L=1$), both gradient scale matrices
reduce to the identity matrix and the dynamics in (3) is reduced to the known
result that ${\bm{W}}^{(t+1)}={\bm{W}}^{(t)}-\mu G_{r}^{(t)}$ (e.g., Le Cun et
al., 1991). Recall, however, that the focus of this paper is the over-
parameterized linear model with $L>1$, in which the loss is not convex. Since
the difference between the convex linear model and the over-parametrized deep
model boils down to these matrices, our convergence analysis henceforth
focuses on the dynamics of the _gradient scale matrices_. In accordance, we
analyze the evolution of the _gradient scale matrices_ as learning proceeds.

###### Theorem 6

Let ${\bm{W}}$ denote the compact representation of a deep linear network,
where $\mathbbm{m}=\min\left({m_{1},...,m_{L-1}}\right)$ denotes the size of
its smallest hidden layer and $\mathbbm{m}\geq\\{m_{0},m_{L}\\}$. At each
layer $l$, assume weight initialization $W^{(0)}_{l}$ obtained by sampling
from a distribution with mean $0$ and variance $\sigma_{l}^{2}$, normalized as
specified in Def. 4. Let
$\left({B}^{(t)}_{l}(\mathbbm{m})\right)_{\mathbbm{m}=1}^{\infty}$ and
$\left({A}^{(t)}_{l}(\mathbbm{m})\right)_{\mathbbm{m}=1}^{\infty}$ denote two
sequences of gradient scale matrices as defined in (2), where the
$\mathbbm{m}^{th}$ element of each series corresponds to a network whose
smallest hidden layer has $\mathbbm{m}$ neurons. Let $\xrightarrow{p}$ denote
element-wise convergence in probability as $\mathbbm{m}\to\ \infty$. Then
$\forall{t},l$:

$\begin{split}&{B}^{(t)}_{l}(m)\xrightarrow{p}I+O(\mu^{2})~{}~{}\forall{l}\in[1\ldots{L}-1],\qquad{B}^{(t)}_{{L}}(m)\xrightarrow{p}\frac{m_{L}}{\mathbbm{m}}I+O(\mu^{2})\xrightarrow{p}O(\mu^{2})\\\
&\mathrm{var}[{B}^{(t)}_{l}(m)]=O(\mu^{2})+O\left(\frac{1}{\mathbbm{m}}\right)\qquad\forall{l}\end{split}$
(4)

and

$\begin{split}&{A}^{(t)}_{l}(m)\xrightarrow{p}I+O(\mu^{2})~{}~{}\forall{l}\in[1\ldots{L}-1],\qquad{A}^{(t)}_{0}(\mathbbm{m})\xrightarrow{p}\frac{m_{0}}{\mathbbm{m}}I+O(\mu^{2})\xrightarrow{p}O(\mu^{2})\\\
&\mathrm{var}[{A}^{(t)}_{l}(m)]=O(\mu^{2})+O\left(\frac{1}{\mathbbm{m}}\right)\qquad\forall{l}\end{split}$
(5)

##### Proof sketch.

The proof goes by way of induction on time $t$. To begin with, we recall that
all the weight matrices $\\{W^{(0)}_{l}\\}_{{l}=1}^{L}$ are initialized by
sampling from a distribution with mean $0$ and variance
$\sigma_{l}^{2}=O(\frac{1}{\mathbbm{m}})$. In Appendix A.1 we prove some
statistical properties of such matrices, and their corresponding _gradient
scale matrices_ $B_{l}^{(0)}$ and $A_{l}^{(0)}$. This analysis shows that (4)
and (5) are true at $t=0$. To prove the induction step, we resort to results
stated in Appendix A.2, where we analyze random matrices that are defined
similarly to the gradient scale matrices, and whose dynamics correspond to the
update rule defined in (3). The main result can be used to show directly that
if the theorem’s assertion holds for $B_{{l}}^{(t)}$, it also holds for
$B_{{l}}^{(t+1)}$. The full proof and details are given in Appendix B.2,

##### Relation to the one-layer linear model.

From Thm 6, if the hidden layers are sufficiently wide (as measured by
$\mathbbm{m}$), we can assume that $B_{{l}}^{(t)}(\mathbbm{m})\approx I$ and
$A_{{l}}^{(t)}(\mathbbm{m})\approx I$ $\forall{l}$. In this case, the dynamics
of the over-parameterized model is identical to the dynamics of the convex
linear model, ${\bm{W}}^{(t+1)}={\bm{W}}^{(t)}-\mu G_{r}^{(t)}$. This is shown
more formally below in Thm 7.

##### Note about convergence, and convergence rate.

In Thm 6 we prove two asymptotic results: ${B}^{(t)}_{l}(m)\xrightarrow{p}I$
and ${A}^{(t)}_{l}(m)\xrightarrow{p}I$ up to small deviations of magnitude
$O(\mu^{2})$. This is true at each gradient step and $\forall t$, irrespective
of the convergence of the network. The only thing that changes with $t$ is the
magnitude of the constants governing this asymptotic behavior, hidden in the
notations $O\left(\frac{1}{\mathbbm{m}}\right)$ and $O(\mu^{2})$ \- this
magnitude increases with time. For a fixed network, there will be a time when
the asymptotic behavior is no longer applicable, and this may happen before or
after the network has converged to its final solution.

Importantly, these constants are significantly different in the two sequences.
Specifically, in (4) and (5) we see that the convergence of ${B}^{(t)}_{L}(m)$
is governed by $O\left(\frac{K}{\mathbbm{m}}\right)$, while the convergence of
${A}^{(0)}_{l}(m)$ is governed by $O\left(\frac{q}{\mathbbm{m}}\right)$.
Typically $q\gg K$, as $q$ is the dimension of the data space which can be
fairly large, while $K$ is the typically small number of classes. Considered
in the context of the proof of Thm 6, the induction implies that these
constants are amplified in each iteration $t$. As a result, when $\mathbbm{m}$
is large but fixed, ${A}^{(t)}_{l}(m)$ is expected to deviate from $I$ sooner
than ${B}^{(t)}_{l}(m)$ as the GD steps are reiterated.

##### Empirical validation.

To understand the practical significance of the observations stated above,
about the difference in convergence rate between the sequences
${B}^{(t)}_{l}(m)$ and ${A}^{(t)}_{l}(m)$, we resort to simulations whose
results are shown in Fig. 3. These empirical results, recounting linear
networks with 4 hidden layers of width 1024, clearly show that during a
significant part of the training both _gradient scale matrices_ remain
approximately $I$. The difference between the convergence rate of
$B_{{l}}^{(t)}$ and $A_{{l}}^{(t)}$ is seen later on, when $A_{{l}}^{(t)}$
starts to deviate from $I$ shortly before convergence, while $B_{{l}}^{(t)}$
remains essentially the same throughout.

(a) Diagonal layer 3

(b) Off-diagonal layer 3

Figure 1: The dynamics of $B_{3}^{(t)}$ and $A^{(t)}_{3}$ when training $10$
$5$-layered linear networks on the small mammals dataset (see App. D.4). (a)
The empirical $L_{2}$-distance of the diagonal elements of $B_{3}^{(t)}$ and
$A^{(t)}_{3}$ from333More precisely, the distance is computed respectively to
their analytical values of $\alpha_{3}^{(t)}$ and $\beta_{3}^{(t)}$, computed
as $diag\left(A^{(t)}_{3}-\alpha_{3}^{(t)}I\right)$ and
$diag\left(B_{3}^{(t)}-\beta_{3}^{(t)}I\right)$, where
$\alpha_{i}^{(t)},\beta_{i}^{(t)}$ are defined in Thm 10, §A.1.$~{}1$. (b) The
empirical $L_{2}$-distance of the off-diagonal elements of $B_{3}^{(t)}$ and
$A^{(t)}_{3}$ from 0. The networks reach maximal test accuracy in epoch
$t=100$, before the divergence of $A^{(t)}_{3}$. In these experiments
$B_{i}^{(t)}$ never reaches the point of divergence. These result are typical
of all the layers in the networks, see Fig. 15 in Appendix C.2.

###  3.3 Weight Evolution

Next, we investigate the implications of Thm 6 in regard to the evolution of
the compact representation ${\bm{W}}$. In Thm 7 we show that at the beginning
of training, when the assumption that both $B^{(t)}_{l}\approx I$ and
$A^{(t)}_{l}\approx I$ is applicable, this evolution resembles the single-
layer linear model and is governed by the eigendecomposition of the data. In
Thm 8 we show that later on, when only $B^{(t)}_{l}\approx I$ is applicable
and the evolution no longer resembles the single-layer linear model, it is
still governed by the eigendecomposition of the data.

More specifically, let $\bm{w}_{j}^{opt}$ denote the $j^{\mathrm{th}}$ column
of the optimal solution of (1), and $\bm{w}_{j}^{(0)}$ the $j^{\mathrm{th}}$
column of the initial weight matrix. Thm 7 states that if the hidden layers of
the network are all wider than a certain fixed width $\hat{\mathbbm{m}}$ and
if the learning rate is slower than $\hat{\mu}$, then at the beginning of
learning (time $t\in[0\ldots\hat{t}]$), the $j^{\mathrm{th}}$ column of the
compact representation ${\bm{W}}^{(t)}$ is roughly the following linear
combination of $\bm{w}_{j}^{(0)}$ and $\bm{w}_{j}^{opt}$:

$\bm{w}_{j}^{(t+1)}\approx\lambda_{j}^{t}\bm{w}_{j}^{(0)}+[1-\lambda_{j}^{t}]\bm{w}_{j}^{opt}\quad\quad\lambda_{j}=1-\mu
d_{j}L$ (6)

with high probability (at least $(1-\delta)$) and low error (at most
$\varepsilon$). Constant $\lambda_{j}$ depends on the $j^{\mathrm{th}}$
singular value of the data $d_{j}$ and the number of layers $L$.

Result (6) is reminiscent of the well-understood dynamics of training the
convex one-layer linear model. It is composed of two additive terms, revealing
two parallel and independent processes:  (i) The dependence on random
initialization tends to 0 exponentially fast as a function of time ${t}$, when
$\lambda_{j}^{t}$ goes to $0$.  (ii) The optimal solution is reached
exponentially fast as a function of time ${t}$, when $1-\lambda_{j}^{t}$ goes
to $1$.  In either case, convergence is fastest for the largest singular
eigenvalue, or the first column of ${\bm{W}}$, and slowest for the smallest
singular value. This behavior is visualized in Fig. 2(a).

(a) 2-layer linear network

(b) 2-layer ReLU network

Figure 2: Empirical evaluation of the dependence of convergence rate on the
eigendecomposition in a _binary_ classification problem. Here the classifier
is a vector, denoted $\bm{w}$. Each line corresponds to a specific principal
eigenvector, plotting in log-log scale the std of element $w_{i}$
($i=1,10,20,30,100,3072$, see legend) over 10 independently trained networks.

Figure 3: Mean accuracy of 10 st-VGG networks evaluated on test data projected
to dimensionality $\\{1,10,100,1000\\}$.

Formally, the theorem can be stated as follows (see proof in Appendix B.2):

###### Theorem 7

Let $\bm{w}_{j}^{(t)}$ denote the $j^{\mathrm{th}}$ column of the compact
representation matrix ${\bm{W}}^{(t)}$, and $\bm{w}_{j}^{opt}$ the
$j^{\mathrm{th}}$ column of the optimal solution of (1). Assume that the data
is rotated to its principal coordinate system, and let $d_{j}$ denote the
$j^{\mathrm{th}}$ singular value of the data. Then there exists $\hat{t}$ such
that $\forall\delta,\varepsilon$ and $\forall{t}\leq\hat{t}$,
$\exists\hat{\mathbbm{m}},\hat{\mu}$ such that
$\forall\mu<\hat{\mu},\mathbbm{m}\geq\hat{\mathbbm{m}}$

$Prob\bigg{(}\Big{\|}\bm{w}_{j}^{(t+1)}-\big{[}\lambda_{j}^{t}\bm{w}_{j}^{(0)}+[1-\lambda_{j}^{t}]\bm{w}_{j}^{opt}\big{]}\Big{\|}<\varepsilon\bigg{)}>(1-\delta)\quad\quad\lambda_{j}=1-\mu
d_{j}L$ (7)

##### Proof sketch.

We first prove by induction on time $t$ that $\forall t$,
${\bm{W}}^{(t+1)}={\bm{W}}^{(t)}-\mu LG_{r}^{(t)}$ with high probability and
small error. We then shift to the principal coordinate system, and note that
this is true separately for each column of matrix ${\bm{W}}$. The evolution of
each column is then shown to be a telescopic series, whose solution is
$\bm{w}_{j}^{(t)}=\lambda_{j}^{t}\bm{w}_{j}^{(0)}+[1-\lambda_{j}^{t}]\bm{w}_{j}^{opt}$
for $\lambda_{j}=1-\mu Ld_{j}$, with high probability and small error.

Next we state Thm 8, which is relevant at later stages of learning as it
imposes no asymptotic constraints on the evolution of ${A}^{(t)}_{l}$. This is
the case discussed in the proof of Thm 6 and illustrated in Fig. 3, when
${B}^{(t)}_{l}\approx I$ remains a good approximation while ${A}^{(t)}_{l}\neq
I$. More specifically, let $A^{(t)}=\sum_{{l}=1}^{L}A_{{l}}^{(t)}$. Thm 8
asserts that if the hidden layers of the network are all wider than
$\breve{\mathbbm{m}}$ and if the learning rate is slower than $\breve{\mu}$,
then at time $t\in[0\ldots\breve{t}]$,444Presumably, in comparison to Thm 7
and if we fix the network’s width $\breve{m}=\hat{m}$, then
$\breve{t}\gg\hat{t}$. the $j^{\mathrm{th}}$ column of the compact
representation ${\bm{W}}^{(t)}$ is approximately the following

$\bm{w}_{j}^{(t+1)}\approx\prod_{{t^{\prime}}=1}^{{t}}(I-\mu
d_{j}A^{({t^{\prime}})})\bm{w}_{j}^{(0)}+\mu
d_{j}\left[\sum_{{t^{\prime}}=1}^{{t}}\prod_{{t^{\prime\prime}}={t^{\prime}}+1}^{{t}}(I-\mu
d_{j}A^{({t^{\prime\prime}})})A^{({t^{\prime}})}\right]\bm{w}_{j}^{opt}$

with high probability (at least $(1-\delta)$) and low error (at most
$\varepsilon$). Note that the $j^{\mathrm{th}}$ column of the compact
representation ${\bm{W}}^{(t)}$ is still a linear combination of
$\bm{w}_{j}^{(0)}$ and $\bm{w}_{j}^{opt}$, but now the coefficients are much
more complex and depend on $A_{{l}}^{(t)}$.

Formally, the theorem can be stated as follows (see proof in Appendix B.2):

###### Theorem 8

Let $\bm{w}_{j}^{(t)}$ denote the $j^{\mathrm{th}}$ column of the compact
representation matrix ${\bm{W}}^{(t)}$, and $\bm{w}_{j}^{opt}$ the
$j^{\mathrm{th}}$ column of the optimal solution of (1). Assume that the data
is rotated to its principal coordinate system, and let $d_{j}$ denote the
$j^{\mathrm{th}}$ singular value of the data. Then there exists $\breve{t}$
such that $\forall\delta,\varepsilon$ and $\forall{t}\leq\breve{t}$,
$\exists\breve{\mathbbm{m}},\breve{\mu}$ such that
$\forall\mu<\breve{\mu},\mathbbm{m}\geq\breve{\mathbbm{m}}$

$\begin{split}Prob\Bigg{(}\bigg{\|}\bm{w}_{j}^{(t+1)}-&\bigg{[}\prod_{{t^{\prime}}=1}^{{t}}(I-\mu
d_{j}A^{({t^{\prime}})})\bm{w}_{j}^{(0)}+\\\ &\mu
d_{j}\Big{(}\sum_{{t^{\prime}}=1}^{{t}}\prod_{{t^{\prime\prime}}={t^{\prime}}+1}^{{t}}(I-\mu
d_{j}A^{({t^{\prime\prime}})})A^{({t^{\prime}})}\Big{)}\bm{w}_{j}^{opt}\bigg{]}\bigg{\|}<\varepsilon\Bigg{)}>(1-\delta)\end{split}$
(8)

Although the dynamics now depend on matrices $A^{(t)}_{l}$ as well, it is
still the case that the convergence of each column is governed by its singular
value $d_{j}$. This suggests that while the _PC-bias_ is more pronounced at
earlier stages of learning, its effect persists throughout.

###  3.4 Adding ReLU Activation

We extend the results above to a relatively simple non-linear model suggested
and analyzed by Allen-Zhu et al. (2019); Arora et al. (2019b); Basri et al.
(2019), here trained on a classification task rather than a regression task.
Specifically, it is a two-layer model with ReLU activation, where only the
weights of the first layer are being learned. Similarly to (1), the loss is
defined as

$L(\mathbb{X})=\frac{1}{2}\sum_{i=1}^{n}\|f(\bm{x}_{i})-\bm{y}_{i}\|^{2}\qquad
f(\bm{x}_{i})=\bm{a}^{\intercal}\cdot\sigma(W\bm{x}_{i}),~{}~{}\bm{a}\in\mathbb{R}^{\mathbbm{m}},~{}W\in\mathbb{R}^{\mathbbm{m}\times
d}$

$\mathbbm{m}$ denotes the number of neurons in the hidden layer. We consider a
binary classification problem with 2 classes, where $y_{i}=1$ for
$\bm{x}_{i}\in C_{1}$, and $y_{i}=-1$ for $\bm{x}_{i}\in C_{2}$.
$\sigma(\cdot)$ denotes the ReLU activation function $\sigma(u)=\max(u,0)$.
Unlike deep linear networks, the analysis in this case requires two additional
symmetry assumptions:

  1. 1.

Each point $\bm{x}_{i}$ is drawn from a symmetric distribution $\mathcal{D}$
with density $f_{\mathcal{D}}({\mathbf{X}})$, such that:
$f_{\mathcal{D}}(\bm{x}_{i})=f_{\mathcal{D}}(-\bm{x}_{i})$.

  2. 2.

$W$ and $\bm{a}$ are initialized symmetrically so that
$\bm{w}^{(0)}_{{\textstyle\mathstrut}2i}=-\bm{w}^{(0)}_{{\textstyle\mathstrut}2i-1}$
and $a_{{\textstyle\mathstrut}2i}=-a_{{\textstyle\mathstrut}2i-1}$ $\forall
i\in[\frac{\mathbbm{m}}{2}]$.

###### Theorem 9

Retaining the assumptions stated above, at the beginning of the learning, the
temporal dynamics of the model can be shown to obey the following update rule:

$W^{(t+1)}\approx
W^{(t)}-\mu\frac{1}{2}\Big{[}(\bm{a}\bm{a}^{\intercal})W^{(t)}\Sigma_{\scriptscriptstyle
XX}-\tilde{M}^{(t)}\Big{]}$

Above $\tilde{M}^{(t)}$ denotes the difference between the centroids of the 2
classes, computed in the half-space defined by
$\bm{w}^{(t)}_{r}\cdot\bm{x}\geq 0$.

The complete proof can be found in App. B.3. Thm 9 is reminiscent of the
single-layer linear model dynamics ${\bm{W}}^{(t+1)}={\bm{W}}^{(t)}-\mu
G_{r}^{(t)}$, and we may conclude that when it holds, using the principal
coordinate system, the rate of convergence of the $j$-th column of $W^{(t)}$
is governed by the singular value $d_{j}$. Fig. 2(b) empirically demonstrated
this result, showing the weight convergence of $10$ ReLU models trained on the
small mammals dataset (5 classes out of CIFAR-100, see App. D.4 for details),
along different principal directions of the data.

##  4 PC-Bias: Empirical Study

In this section, we first analyze deep linear networks, showing that the
convergence rate is indeed governed by the principal singular values of the
data, which demonstrates the plausibility of the assumptions made in Section
3. We continue by extending the scope of the investigation to non-linear
neural networks, finding evidence for the _PC-bias_ mostly in the earlier
stages of learning.

###  4.1 Methodology

We say that a linear network is $L$-layered when it has $L-1$ hidden fully
connected (FC) layers (without convolutional layers). In our empirical study,
we relaxed some assumptions of the theoretical study, to increase the
resemblance of the trained networks to networks in common use. Specifically,
we changed the initialization to the commonly used Glorot initialization,
replaced the $L_{2}$ loss with the cross-entropy loss, and employed SGD
instead of the deterministic GD. As to be expected, the original assumptions
yielded similar results (see §C.3.1). The results presented summarize
experiments with networks of equal width across all hidden layers, fixing the
moderate555Note that for most image datasets, $q$ (the input dimension) is
much larger than $1024$, and therefore $\frac{q}{\mathbbm{m}}>>1$. In this
sense $\mathbbm{m}=1024$ is moderate, as one can no longer assume
$A_{l}^{(t)}(\mathbbm{m})\approx I$. value of $\mathbbm{m}=1024$ in order to
assess the relevance of our asymptotic results obtained with
$\mathbbm{m}\to\infty$. Using a different width for each layer yielded similar
qualitative results. Details regarding the hyper-parameters, architectures,
and datasets can be found in §D.1, §D.3 and §D.4 respectively.

###  4.2 PC-bias in Deep Linear Networks

In this section, we train $L$-layered linear networks, then compute their
compact representations ${\bm{W}}$ rotated to align with the canonical
coordinate system (Def. 1). Note that each row $\bm{w}(r)$ in ${\bm{W}}$
essentially defines the one-vs-all separating hyper-plane corresponding to
class $r$.

To examine both the variability between models and their convergence rate, we
inspect $\bm{w}(r)$ at different time points during learning. The rate of
convergence can be measured directly, by observing the changes in the weights
of each element in $\bm{w}(r)$. These weight values666We note that the weights
tend to start larger for smaller principal components, see Fig. 5(a) left. are
compared with the corresponding optimal weights $w_{opt}$. The variability
between models is measured by calculating the standard deviation (std) of each
row vector $\bm{w}(r)$ across $N$ models, obtained with different random
initializations.

We begin with linear networks. We trained $10$ $5$-layered FC linear networks,
and $10$ linear st-VGG convolutional networks. When analyzing the compact
representation of such networks we observe a similar behavior – weights
corresponding to larger principal components converge faster to the optimal
value, and their variability across models converges faster to 0 (Figs.
5(a),5(b)). Thus, while the theoretical results are asymptotic, the _PC-bias_
is empirically seen throughout the entire learning process of deep linear
networks.

We note that in the principal coordinate system, the absolute values of the
optimal solution are often higher in directions corresponding to lower
principal components. In App. C.3.2, we discuss this phenomenon in detail. We
report a faster convergence in the directions corresponding to the higher
principal components when normalizing this effect.

(a) 5-Layered linear network

(b) Linear convolutional network

(c) 5-Layered linear network, whitened data

(d) Non-linear convolutional network

Figure 5: Convergence of the compact representation along the principal
directions in different epochs. The value of the $X$-axis corresponds to the
index of a principal eigenvalue, from the most significant to the least
significant. (a) $10$ $5$-layered linear networks trained on the cats and dogs
dataset. 3 plots are provided, corresponding to snapshots taken at different
stages of learning: the beginning (epoch $1$, left), intermediate stage (epoch
$100$, middle), and close to convergence (epoch $1500$ right). Bottom panel:
average distance of the weights in $\bm{w}(1)$ from the optimal linear
classifier; top panel: respective std. (b) Similarly, for $10$ linear st-VGG
convolutional networks, trained on CIFAR-10 (epochs plotted: $1,100,500$). (c)
Similarly, for $10$ $5$-layered linear networks, trained on the cats and dogs
dataset, with ZCA-whitening (epochs plotted: $1,2,1000$). (d) Similarly, for
$10$ non-linear st-VGG networks trained on the cats and dogs dataset. Here the
distance to the optimal solution is not well defined and we therefore only
show the std (epochs plotted: $1,10,100$). In all plots, the data dimension is
$3072$, while only the largest $2000$ principal components are shown for
clarity. A visualization of the values of $\lambda_{i}$ can be found in Fig.
18.

##### Whitened data.

The _PC-bias_ is neutralized when the data is whitened, at which point
$\Sigma_{\scriptscriptstyle XX}$ is the scaled identity matrix. In Fig. 5(c),
we plot the results of the same experimental protocol while using a ZCA-
whitened dataset. As predicted, the networks no longer show any bias towards
any principal direction. Weights in all directions are scaled similarly, and
the std over all models is the same in each epoch, irrespective of the
principal direction. (Additional experiments show that this is _not_ an
artifact, due to the lack of uniqueness when deriving the principal
eigenvectors of a white signal).

###  4.3 PC-Bias in General CNNs

In this section, we investigate the manifestation of the _PC-bias_ in non-
linear deep convolutional networks. As we cannot directly track the learning
dynamics separately in each principal direction of non-linear networks, we
adopt two different evaluation mechanisms:

##### (I) Linear approximation.

We considered several linear approximations, but since all of them showed the
same qualitative behavior, we report results with the simplest one.
Specifically, to obtain a linear approximation of a non-linear network,
without max-pooling or batch-normalization layers, we follow the definition of
the compact representation from Section 3 while ignoring any non-linear
activation. We then align this matrix with the canonical coordinate system
(Def. 1), and observe the evolution of the weights and their std across models
along the principal directions during learning. Note that now the networks do
not converge to the same compact representation, which is not unique.
Nevertheless, we see that the _PC-bias_ governs the weight dynamics to a
noticeable extent.

More specifically, in these networks, a large fraction of the lowest principal
components hardly ever changes during learning. Nevertheless, the _PC-bias_
affects the higher principal components, most notably at the beginning of
training (see Fig. 5(d)). Thus weights corresponding to higher principal
components converge faster, and the std across models of such weights
decreases faster for higher principal components.

##### (II) Projection to higher PC’s.

We created a modified _test-set_ , by projecting each test example on the span
of the first $P$ principal components. This is equivalent to reducing the
dimensionality of the test set to $P$ using PCA. We trained an ensemble of
$N$=$100$ st-VGG networks on the original training set of the small mammals
dataset (see App. D.4), then evaluated these networks during training on 4
versions of the test-set, reduced to $P$=$1$,$10$,$100$,$1000$ dimensions
respectively. Mean accuracy is plotted in Fig. 3. Similar results are obtained
when training VGG-19 networks on CIFAR-10, see §C.4.

Taking a closer look at Fig. 3, we see that when evaluated on lower
dimensionality test-data ($P$=$1$,$10$), the networks’ accuracy peaks after a
few epochs, at which point performance starts to decrease. This result
suggests that the networks rely more heavily on these dimensions in the
earlier phases of learning, and then continue to learn other things. In
contrast, when evaluated on higher dimensionality test-data
($P$=$100$,$1000$), accuracy continues to rise, longer so for larger $P$. This
suggests that significant learning of the additional dimensions continues in
later stages of the learning.

##  5 PC-Bias: Learning Order Constancy

In this section, we show that the _PC-bias_ is significantly correlated with
the learning order of deep neural networks, and can therefore partially
account for the _LOC-effect_ described in Section 1. Following Hacohen et al.
(2020), we measure the "speed of learning" of each example by computing its
_accessibility_ score. This score is computed per example, and characterizes
how fast an ensemble of $N$ networks learns it. Formally,
$accessibility(\bm{x})=\mathbb{E}\left[\mathbbm{1}(f_{i}^{e}(\bm{x})=y(\bm{x})\right)]$,
where $f_{i}^{e}(\bm{x})$ denotes the outcome of the $i$-th network trained
over $e$ epochs, and the mean is taken over networks and epochs. For the set
of datapoints $\\{(\bm{x}_{j},\bm{y}_{j})\\}_{j=1}^{n}$, _Learning Order
Constancy_ is manifested by the high correlation between 2 instances of
$accessibility(\bm{x})$, each computed from a different ensemble.

_PC-bias_ is shown to pertain to _LOC_ in two ways: First, in Section 5.1 we
show high correlation between the learning order in deep linear and non-linear
networks. Since the _PC-bias_ fully accounts for _LOC_ in deep linear
networks, this suggests that it also accounts (at least partially) for the
observed _LOC_ in non-linear networks. Comparison with the _critical principal
component_ verifies this assertion. Second, we show in Section 5.2 that when
the _PC-bias_ is neutralized, _LOC_ diminishes as well. In Section 5.3 we
discuss the relationship between the _spectral bias_ , _PC-bias_ and the _LOC-
effect_.

###  5.1 PC-Bias is Correlated with LOC

We first compare the order of learning of non-linear models and deep linear
networks by computing the correlation between the _accessibility_ scores of
both models. This comparison reveals high correlation ($r=0.85$,
$p<10^{-45}$), as seen in Fig. 6(a). To investigate directly the connection
between the _PC-bias_ and _LOC_ , we define the _critical principal component_
of an example to be the first principal component $P$, such that a linear
classifier trained on the original data can classify the example correctly
when projected to $P$ principal components. We trained $N$=$100$ st-VGG
networks on the cats and dogs dataset (2 classes out of CIFAR-10, see §D.4 for
details), and computed for each example its _accessibility_ score and its
_critical principal component_. In Fig. 6(b) we see a strong negative
correlation between the two scores ($p$=$-0.93$, $r\\!<\\!10^{-4}$),
suggesting that the _PC-bias_ affects the order of learning as measured by
_accessibility_.

(a)

(b)

Figure 6: (a) Correlation between the _accessibility_ score of $N$=$100$ st-
VGG networks trained with a low learning rate88footnotemark: 8, and $N$=$100$
linear st-VGG networks, trained on small mammals (see §D.3, §D.4 for details).
(b) Correlation between the accessibility score of $N$=$100$ st-VGG networks
trained on cats and dogs, and the _critical principal component_ score. The
_accessibility_ plot is smoothed by moving average of width 10. Error bars
indicate standard error.

(c) Linear networks

(d) Non-linear networks

Figure 7: _LOC_ as measured with and without _PC-bias_. Each bar represents
the correlation between the learning order of 2 collections of $10$ networks
trained on CIFAR-10. Orange bars represent natural images, in which the _PC-
bias_ is present, while blue bars represent whitened data, in which the _PC-
bias_ is eliminated. As _PC-bias_ is more prominent earlier on, we compare
these correlations for the entire data (right 2 bars), and for the subset of
$20\%$ "fastest learned" examples (left 2 bars).

###  5.2 Neutralizing the PC-bias Leads to Diminishing LOC

Whitening the data eliminates the _PC-bias_ as shown in Fig. 5(c), since all
the singular values are now identical. Here we use this observation to further
probe into the dependency of the _Learning Order Constancy_ on the _PC-bias_.
Starting with the linear case, we trained 4 ensembles of $N$=$10$ $2$-layered
linear networks (following the same methodology as Section 4.1) on the cats
and dogs dataset, $2$ with and $2$ without ZCA-whitening. We computed the
_accessibility_ score for each ensemble separately, and correlated the scores
of the $2$ ensembles in each test case. Each correlation captures the
consistency of the _LOC-effect_ for the respective condition. This correlation
is expected to be very high for natural images. Low correlation implies that
the _LOC-effect_ is weak, as training the same network multiple times yields a
different learning order.

Fig. 6(c) shows the results for deep linear networks. As expected, the
correlation when using natural images is very high. However, when using
whitened images, correlation plummets, indicating that the _LOC-effect_ is
highly dependent on the _PC-bias_. We note that the drop in the correlation is
much higher when considering only the $20\%$ "fastest learned" examples,
suggesting that the _PC-bias_ affects learning order more strongly at earlier
stages of learning.

Fig. 6(d) shows the results when repeating this experiment with non-linear
networks, training $2$ collections of $N$=$10$ VGG-19 networks on CIFAR-10. We
report that the elimination of the _PC-bias_ in this case affects _LOC_ much
less, suggesting that the _PC-bias_ can only partially account for the _LOC-
effect_ in the non-linear case. However, we note that at the beginning of
learning, when the _PC-bias_ is most pronounced, once again the drop is much
larger and very significant (half).

###  5.3 Spectral Bias, PC-bias and LOC

The _spectral bias_ (Rahaman et al., 2019) characterizes the dynamics of
learning in neural networks differently, asserting that initially neural
models can be described by low frequencies only. This may provide an
alternative explanation to LOC. Recall that LOC is manifested in the
consistency of the _accessibility_ score across networks. To compare between
the _spectral bias_ and _accessibility_ score, we first need to estimate for
each example whether it can be correctly classified by a low frequency model.
Accordingly, we define for each example a _discriminability_ measure – the
percentage out of its $k$ neighbors that share with it class identity.
Intuitively, an example has a low _discriminability_ score when it is
surrounded by examples from other classes, which would presumably force the
learned boundary to incorporate high frequencies. In §C.5 we show that in the
2D case analyzed by Rahaman et al. (2019), this measure strongly correlates
($r$=$-0.8$, $p<10^{-2}$) with the spectral bias.

We trained several networks (VGG-19 and st-VGG) on several real datasets,
including small mammals, STL-10, CIFAR-10/100, and a subset of ImageNet-20
(see App. D.4). For each network and dataset, we computed the _accessibility_
score as well as the _discriminability_ of each example. The vector space, in
which discriminability is evaluated, is either the raw data or the network’s
perceptual space (penultimate layer activation). The correlation between these
scores is shown in Table 5.3.

Table 1: Correlation between _accessibility_ and _discriminability_.

Dataset | Raw data | Penultimate
---|---|---
Small mammals | 0.46 | 0.85
ImageNet 20 | 0.01 | 0.54
CIFAR-100 | 0.51 | 0.85
STL10 | 0.44 | 0.7

Using raw data, low correlation is still seen between the _accessibility_ and
_discriminability_ scores when inspecting the smaller datasets (small mammals,
CIFAR-100, and STL10). This correlation vanishes when considering the larger
ImageNet-20 dataset. It would appear that on its own, the _spectral bias_
cannot adequately explain the _LOC-effect_. On the other hand, in the
perceptual space, the correlation between _discriminability_ and
_accessibility_ is quite significant for all datasets. Differently from the
supposition of the spectral bias, it seems that networks learn a
representation where the _spectral bias_ is evident, but this bias does not
necessarily govern its learning before the representation has been obtained.

##### Discussion.

In the limit of infinite width, the spectral bias phenomenon in deep linear
networks follows from the PC-bias. In function space, the training process of
neural networks can be decomposed along different directions defined by the
eigenfunctions of the neural tangent kernel, where each direction has its
convergence rate and the rate is determined by the corresponding eigenvalue
(Cao et al., 2021). Since the Neural Tangent Kernel for a deep linear network
is proportional to $\Sigma_{XX}$, the spectral bias, in this case, corresponds
to the statement that the right singular vectors of $X$ are learned at rates
corresponding to their (squared) singular values. Looking at the spectral
decomposition of ${\bm{W}}X$, we see that the PC Bias implies in function
space that the right singular vectors of $X$ are learned at rates
corresponding to their singular values. Thus the PC Bias implies the Spectral
Bias in this model.

Figure 8: Visualization of the small mammals dataset, with amplification of
$1.5\%$ of its principal components by a factor of $10$. Top: original data;
middle: data amplified along the highest principal components; bottom: data
amplified along the lowest principal components

##  6 PC-Bias: Further Implications

##### Early Stopping and the Generalization Gap.

Considering natural images, it is often assumed that the least significant
principal components of the data represent noise (Torralba and Oliva, 2003).
In such cases, our analysis predicts that as noise dominates the components
learned later in learning, early stopping is likely to be beneficial. To test
this hypothesis directly, we manipulated CIFAR-10 to amplify the signal in
either the $1.5\%$ most significant (higher) or $1.5\%$ least significant
(lower) principal components (see examples in Fig. 8). Accuracy over the
original test set, after training $10$ st-VGG and linear st-VGG networks on
these manipulated images, can be seen in Fig. 10. Both in linear and non-
linear networks, early stopping is more beneficial when lower principal
components are amplified, and significantly less so when higher components are
amplified, as predicted by the _PC-bias_.

##### Slower Convergence with Random Labels.

Deep neural models can learn any random label assignment to a given training
set (Zhang et al., 2017). However, when trained on randomly labeled data,
convergence appears to be much slower (Krueger et al., 2017). Assume, as
before, that in natural images the lower principal components are dominated by
noise. We argue that the _PC-bias_ now predicts this empirical result, since
learning randomly labeled examples requires a signal present in lower
principal components. To test this hypothesis directly, we trained $10$
$2$-layered linear networks (following the same methodology as in Section 4.1)
using datasets of natural images. Indeed, these networks converge slower with
random labels (see Fig. 9(c)). In Fig. 9(d) we repeat this experiment after
having whitened the images, to neutralize the _PC-bias_. Now convergence rate
is identical, whether the labels are original or shuffled. Clearly, in deep
linear networks, the _PC-bias_ gives a full account for this phenomenon.

(a) linear network

(b) non-linear network

Figure 9: Comparing the accuracy trajectory when amplifying the highest (blue
line) and lowest (orange line) principal components.

(c) original data

(d) whitened data

Figure 10: Learning curves when using real and shuffled labels: $10$
$2$-layered linear networks, (a) before and (b) after whitening.

To further check the relevance of this account to non-linear networks, we
artificially generate datasets where only the first $P$ principal components
are discriminative, while the remaining components become noise by design. We
constructed two such datasets: in one the labels are correlated with the
original labels, while in the other they are not. Specifically, PCA is used to
reduce the dimensionality of a two-class dataset to $P$, and the optimal
linear separator in the reduced representation is computed. Next, all the
labels of points that are incorrectly classified by the optimal linear
separator are switched, so that the train and test sets are linearly separable
by this separator. Note that the modified labels are still highly correlated
with the original labels (for $P=500$: $p=0.82$, $r<10^{-10}$). The second
dataset is generated by repeating the process while starting from randomly
shuffled labels. This dataset is likewise fully separable when projected to
the first $P$ components, but its labels are uncorrelated with the original
labels (for $P=500$: $p=0.06$, $r<10^{-10}$).

The mean training accuracy of $10$ non-linear networks with
$P$=$10$,$50$,$500$ is plotted in Fig. 11(a) (first dataset) and Fig. 11(b)
(second dataset). In both cases, the lower $P$ is (namely, only the first few
principal components are discriminative), the faster the data is learned by
the non-linear network. Whether the labels are real or shuffled makes little
qualitative difference, as predicted by the _PC-bias_.

(a) Original labels

(b) Shuffled labels

Figure 11: Learning curves of st-VGG networks trained on 3 datasets, which are
linearly separable after projection to the highest $P$ principal components
(see legend).

##  7 Summary and Discussion

When trained with gradient descent, the convergence rate of the over-
parameterized deep linear network model is provably governed by the
eigendecomposition of the data. Specifically, we show that parameters
corresponding to the most significant principal components converge faster
than the least significant components. Empirical evidence is provided for the
relevance of this result to more realistic non-linear networks. We term this
effect _PC-bias_. This result provides a complementary account for some
prevalent empirical observations, including the benefit of early stopping and
the slower convergence rate with shuffled labels.

We use the _PC-bias_ to explain the _Learning Order Constancy (LOC)_.
Different empirical schemes are used to show that examples learned at earlier
stages are more distinguishable by the data’s higher principal components,
which may indicate that networks’ training relies more heavily on higher
principal components early on. A causal link between the _PC-bias_ and the
_LOC-effect_ is demonstrated, as the _LOC-effect_ diminishes when the _PC-
bias_ is eliminated by whitening the images. Finally, we analyze these
findings given a related phenomenon termed _spectral bias_. While the _PC-
bias_ may be more prominent early on, the _spectral bias_ may be more
important in later stages of learning.

## References

  * Allen-Zhu et al. (2019) Zeyuan Allen-Zhu, Yuanzhi Li, and Yingyu Liang.  Learning and generalization in overparameterized neural networks, going beyond two layers.  In Hanna M. Wallach, Hugo Larochelle, Alina Beygelzimer, Florence d’Alché-Buc, Emily B. Fox, and Roman Garnett, editors, _Advances in Neural Information Processing Systems 32: Annual Conference on Neural Information Processing Systems 2019, NeurIPS 2019, December 8-14, 2019, Vancouver, BC, Canada_ , pages 6155–6166, 2019.
  * Arora et al. (2018) Sanjeev Arora, Nadav Cohen, and Elad Hazan.  On the optimization of deep networks: Implicit acceleration by overparameterization.  In _International Conference on Machine Learning_ , pages 244–253, 2018.
  * Arora et al. (2019a) Sanjeev Arora, Nadav Cohen, Noah Golowich, and Wei Hu.  A convergence analysis of gradient descent for deep linear neural networks.  In _7th International Conference on Learning Representations, ICLR 2019, New Orleans, LA, USA, May 6-9, 2019_. OpenReview.net, 2019a.
  * Arora et al. (2019b) Sanjeev Arora, Simon Du, Wei Hu, Zhiyuan Li, and Ruosong Wang.  Fine-grained analysis of optimization and generalization for overparameterized two-layer neural networks.  In _International Conference on Machine Learning_ , pages 322–332, 2019b.
  * Bah et al. (2019) Bubacarr Bah, Holger Rauhut, Ulrich Terstiege, and Michael Westdickenberg.  Learning deep linear neural networks: Riemannian gradient flows and convergence to global minimizers.  _CoRR_ , abs/1910.05505, 2019.
  * Baldi and Hornik (1989) Pierre Baldi and Kurt Hornik.  Neural networks and principal component analysis: Learning from examples without local minima.  _Neural networks_ , 2(1):53–58, 1989.
  * Basri et al. (2019) Ronen Basri, David W. Jacobs, Yoni Kasten, and Shira Kritchman.  The convergence rate of neural networks for learned functions of different frequencies.  In _Advances in Neural Information Processing Systems_ , pages 4761–4771, 2019.
  * Basri et al. (2020) Ronen Basri, Meirav Galun, Amnon Geifman, David Jacobs, Yoni Kasten, and Shira Kritchman.  Frequency bias in neural networks for input of non-uniform density.  In _International Conference on Machine Learning_ , pages 685–694. PMLR, 2020.
  * Cao et al. (2021) Yuan Cao, Zhiying Fang, Yue Wu, Ding-Xuan Zhou, and Quanquan Gu.  Towards understanding the spectral bias of deep learning.  In Zhi-Hua Zhou, editor, _Proceedings of the Thirtieth International Joint Conference on Artificial Intelligence, IJCAI 2021, Virtual Event / Montreal, Canada, 19-27 August 2021_ , pages 2205–2211. ijcai.org, 2021.  doi: 10.24963/ijcai.2021/304.
  * Chizat et al. (2019) Lénaïc Chizat, Edouard Oyallon, and Francis R. Bach.  On lazy training in differentiable programming.  In Hanna M. Wallach, Hugo Larochelle, Alina Beygelzimer, Florence d’Alché-Buc, Emily B. Fox, and Roman Garnett, editors, _Advances in Neural Information Processing Systems 32: Annual Conference on Neural Information Processing Systems 2019, NeurIPS 2019, December 8-14, 2019, Vancouver, BC, Canada_ , pages 2933–2943, 2019.
  * Dingle et al. (2018) Kamaludin Dingle, Chico Q Camargo, and Ard A Louis.  Input–output maps are strongly biased towards simple outputs.  _Nature communications_ , 9(1):1–7, 2018.
  * Du and Hu (2019) Simon Du and Wei Hu.  Width provably matters in optimization for deep linear neural networks.  In _International Conference on Machine Learning_ , pages 1655–1664. PMLR, 2019.
  * Du et al. (2019) Simon Du, Jason Lee, Haochuan Li, Liwei Wang, and Xiyu Zhai.  Gradient descent finds global minima of deep neural networks.  In _International Conference on Machine Learning_ , pages 1675–1685. PMLR, 2019.
  * Du et al. (2018) Simon S. Du, Wei Hu, and Jason D. Lee.  Algorithmic regularization in learning deep homogeneous models: Layers are automatically balanced.  In Samy Bengio, Hanna M. Wallach, Hugo Larochelle, Kristen Grauman, Nicolò Cesa-Bianchi, and Roman Garnett, editors, _Advances in Neural Information Processing Systems 31: Annual Conference on Neural Information Processing Systems 2018, NeurIPS 2018, December 3-8, 2018, Montréal, Canada_ , pages 382–393, 2018.
  * Fukumizu (1998) Kenji Fukumizu.  Effect of batch learning in multilayer neural networks.  _Gen_ , 1(04):1E–03, 1998.
  * Gidel et al. (2019) Gauthier Gidel, Francis Bach, and Simon Lacoste-Julien.  Implicit regularization of discrete gradient dynamics in linear neural networks.  _Advances in Neural Information Processing Systems_ , 32:3202–3211, 2019.
  * Gissin et al. (2020) Daniel Gissin, Shai Shalev-Shwartz, and Amit Daniely.  The implicit bias of depth: How incremental learning drives generalization.  In _8th International Conference on Learning Representations, ICLR 2020, Addis Ababa, Ethiopia, April 26-30, 2020_. OpenReview.net, 2020.
  * Glorot and Bengio (2010) Xavier Glorot and Yoshua Bengio.  Understanding the difficulty of training deep feedforward neural networks.  In _Proceedings of the thirteenth international conference on artificial intelligence and statistics_ , pages 249–256, 2010.
  * Gunasekar et al. (2018) Suriya Gunasekar, Jason D. Lee, Daniel Soudry, and Nati Srebro.  Implicit bias of gradient descent on linear convolutional networks.  In Samy Bengio, Hanna M. Wallach, Hugo Larochelle, Kristen Grauman, Nicolò Cesa-Bianchi, and Roman Garnett, editors, _Advances in Neural Information Processing Systems 31: Annual Conference on Neural Information Processing Systems 2018, NeurIPS 2018, December 3-8, 2018, Montréal, Canada_ , pages 9482–9491, 2018.
  * Hacohen et al. (2020) Guy Hacohen, Leshem Choshen, and Daphna Weinshall.  Let’s agree to agree: Neural networks share classification order on real datasets.  In _International Conference on Machine Learning_ , pages 3950–3960. PMLR, 2020.
  * He et al. (2015) Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun.  Delving deep into rectifiers: Surpassing human-level performance on imagenet classification.  In _Proceedings of the IEEE International Conference on Computer Vision (ICCV)_ , December 2015.
  * Heckel and Soltanolkotabi (2020) Reinhard Heckel and Mahdi Soltanolkotabi.  Denoising and regularization via exploiting the structural bias of convolutional generators.  In _8th International Conference on Learning Representations, ICLR 2020, Addis Ababa, Ethiopia, April 26-30, 2020_. OpenReview.net, 2020.
  * Hu et al. (2020a) Wei Hu, Lechao Xiao, Ben Adlam, and Jeffrey Pennington.  The surprising simplicity of the early-time learning dynamics of neural networks.  In Hugo Larochelle, Marc’Aurelio Ranzato, Raia Hadsell, Maria-Florina Balcan, and Hsuan-Tien Lin, editors, _Advances in Neural Information Processing Systems 33: Annual Conference on Neural Information Processing Systems 2020, NeurIPS 2020, December 6-12, 2020, virtual_ , 2020a.
  * Hu et al. (2020b) Wei Hu, Lechao Xiao, and Jeffrey Pennington.  Provable benefit of orthogonal initialization in optimizing deep linear networks.  In _8th International Conference on Learning Representations, ICLR 2020, Addis Ababa, Ethiopia, April 26-30, 2020_. OpenReview.net, 2020b.
  * Jacot et al. (2018) Arthur Jacot, Clément Hongler, and Franck Gabriel.  Neural tangent kernel: Convergence and generalization in neural networks.  In Samy Bengio, Hanna M. Wallach, Hugo Larochelle, Kristen Grauman, Nicolò Cesa-Bianchi, and Roman Garnett, editors, _Advances in Neural Information Processing Systems 31: Annual Conference on Neural Information Processing Systems 2018, NeurIPS 2018, December 3-8, 2018, Montréal, Canada_ , pages 8580–8589, 2018.
  * Ji and Telgarsky (2019) Ziwei Ji and Matus Telgarsky.  Gradient descent aligns the layers of deep linear networks.  In _7th International Conference on Learning Representations, ICLR 2019, New Orleans, LA, USA, May 6-9, 2019_. OpenReview.net, 2019.
  * Jin and Montúfar (2020) Hui Jin and Guido Montúfar.  Implicit bias of gradient descent for mean squared error regression with wide neural networks.  _CoRR_ , abs/2006.07356, 2020.
  * Kalimeris et al. (2019) Dimitris Kalimeris, Gal Kaplun, Preetum Nakkiran, Benjamin L. Edelman, Tristan Yang, Boaz Barak, and Haofeng Zhang.  SGD on neural networks learns functions of increasing complexity.  In Hanna M. Wallach, Hugo Larochelle, Alina Beygelzimer, Florence d’Alché-Buc, Emily B. Fox, and Roman Garnett, editors, _Advances in Neural Information Processing Systems 32: Annual Conference on Neural Information Processing Systems 2019, NeurIPS 2019, December 8-14, 2019, Vancouver, BC, Canada_ , pages 3491–3501, 2019.
  * Kawaguchi (2016) Kenji Kawaguchi.  Deep learning without poor local minima.  In Daniel D. Lee, Masashi Sugiyama, Ulrike von Luxburg, Isabelle Guyon, and Roman Garnett, editors, _Advances in Neural Information Processing Systems 29: Annual Conference on Neural Information Processing Systems 2016, December 5-10, 2016, Barcelona, Spain_ , pages 586–594, 2016.
  * Krueger et al. (2017) David Krueger, Nicolas Ballas, Stanislaw Jastrzebski, Devansh Arpit, Maxinder S. Kanwal, Tegan Maharaj, Emmanuel Bengio, Asja Fischer, and Aaron C. Courville.  Deep nets don’t learn via memorization.  In _5th International Conference on Learning Representations, ICLR 2017, Toulon, France, April 24-26, 2017, Workshop Track Proceedings_. OpenReview.net, 2017.
  * Laurent and von Brecht (2017) Thomas Laurent and James von Brecht.  Deep linear neural networks with arbitrary loss: All local minima are global.  _CoRR_ , abs/1712.01473, 2017.
  * Le Cun et al. (1991) Yann Le Cun, Ido Kanter, and A Sara Solla.  Second order properties of error surfaces learning time and generalization.  _Advances in neural information processing systems_ , 3:918–924, 1991.
  * Li et al. (2020) Mingchen Li, Mahdi Soltanolkotabi, and Samet Oymak.  Gradient descent with early stopping is provably robust to label noise for overparameterized neural networks.  In _International Conference on Artificial Intelligence and Statistics_ , pages 4313–4324. PMLR, 2020.
  * Nguyen (2021) Phan-Minh Nguyen.  Analysis of feature learning in weight-tied autoencoders via the mean field lens.  _arXiv preprint arXiv:2102.08373_ , 2021.
  * Pérez et al. (2019) Guillermo Valle Pérez, Chico Q. Camargo, and Ard A. Louis.  Deep learning generalizes because the parameter-function map is biased towards simple functions.  In _7th International Conference on Learning Representations, ICLR 2019, New Orleans, LA, USA, May 6-9, 2019_. OpenReview.net, 2019.
  * Pliushch et al. (2021) Iuliia Pliushch, Martin Mundt, Nicolas Lupp, and Visvanathan Ramesh.  When deep classifiers agree: Analyzing correlations between learning order and image statistics.  _arXiv preprint arXiv:2105.08997_ , 2021.
  * Rahaman et al. (2019) Nasim Rahaman, Aristide Baratin, Devansh Arpit, Felix Draxler, Min Lin, Fred Hamprecht, Yoshua Bengio, and Aaron Courville.  On the spectral bias of neural networks.  In _International Conference on Machine Learning_ , pages 5301–5310. PMLR, 2019.
  * Russakovsky et al. (2015) Olga Russakovsky, Jia Deng, Hao Su, Jonathan Krause, Sanjeev Satheesh, Sean Ma, Zhiheng Huang, Andrej Karpathy, Aditya Khosla, Michael Bernstein, et al.  Imagenet large scale visual recognition challenge.  _International journal of computer vision_ , 115(3):211–252, 2015.
  * Saxe et al. (2014) Andrew M. Saxe, James L. McClelland, and Surya Ganguli.  Exact solutions to the nonlinear dynamics of learning in deep linear neural networks.  In Yoshua Bengio and Yann LeCun, editors, _2nd International Conference on Learning Representations, ICLR 2014, Banff, AB, Canada, April 14-16, 2014, Conference Track Proceedings_ , 2014.
  * Saxe et al. (2019) Andrew M Saxe, James L McClelland, and Surya Ganguli.  A mathematical theory of semantic development in deep neural networks.  _Proceedings of the National Academy of Sciences_ , 116(23):11537–11546, 2019.
  * Shah et al. (2020) Harshay Shah, Kaustav Tamuly, Aditi Raghunathan, Prateek Jain, and Praneeth Netrapalli.  The pitfalls of simplicity bias in neural networks.  In Hugo Larochelle, Marc’Aurelio Ranzato, Raia Hadsell, Maria-Florina Balcan, and Hsuan-Tien Lin, editors, _Advances in Neural Information Processing Systems 33: Annual Conference on Neural Information Processing Systems 2020, NeurIPS 2020, December 6-12, 2020, virtual_ , 2020.
  * Simoncelli and Olshausen (2001) Eero P Simoncelli and Bruno A Olshausen.  Natural image statistics and neural representation.  _Annual review of neuroscience_ , 24(1):1193–1216, 2001.
  * Soudry et al. (2018) Daniel Soudry, Elad Hoffer, Mor Shpigel Nacson, Suriya Gunasekar, and Nathan Srebro.  The implicit bias of gradient descent on separable data.  _The Journal of Machine Learning Research_ , 19(1):2822–2878, 2018.
  * Torralba and Oliva (2003) Antonio Torralba and Aude Oliva.  Statistics of natural image categories.  _Network: computation in neural systems_ , 14(3):391–412, 2003.
  * Ulyanov et al. (2018) Dmitry Ulyanov, Andrea Vedaldi, and Victor Lempitsky.  Deep image prior.  In _Proceedings of the IEEE conference on computer vision and pattern recognition_ , pages 9446–9454, 2018.
  * Williams et al. (2019) Francis Williams, Matthew Trager, Daniele Panozzo, Cláudio T. Silva, Denis Zorin, and Joan Bruna.  Gradient dynamics of shallow univariate relu networks.  In Hanna M. Wallach, Hugo Larochelle, Alina Beygelzimer, Florence d’Alché-Buc, Emily B. Fox, and Roman Garnett, editors, _Advances in Neural Information Processing Systems 32: Annual Conference on Neural Information Processing Systems 2019, NeurIPS 2019, December 8-14, 2019, Vancouver, BC, Canada_ , pages 8376–8385, 2019.
  * Wu et al. (2019a) Lei Wu, Qingcan Wang, and Chao Ma.  Global convergence of gradient descent for deep linear residual networks.  In Hanna M. Wallach, Hugo Larochelle, Alina Beygelzimer, Florence d’Alché-Buc, Emily B. Fox, and Roman Garnett, editors, _Advances in Neural Information Processing Systems 32: Annual Conference on Neural Information Processing Systems 2019, NeurIPS 2019, December 8-14, 2019, Vancouver, BC, Canada_ , pages 13368–13377, 2019a.
  * Wu et al. (2019b) Yifan Wu, Barnabas Poczos, and Aarti Singh.  Towards understanding the generalization bias of two layer convolutional linear classifiers with gradient descent.  In _The 22nd International Conference on Artificial Intelligence and Statistics_ , pages 1070–1078. PMLR, 2019b.
  * Yehudai and Shamir (2019) Gilad Yehudai and Ohad Shamir.  On the power and limitations of random features for understanding neural networks.  _Advances in Neural Information Processing Systems_ , 32:6598–6608, 2019.
  * Yun et al. (2021) Chulhee Yun, Shankar Krishnan, and Hossein Mobahi.  A unifying view on implicit bias in training linear neural networks.  In _9th International Conference on Learning Representations, ICLR 2021, Virtual Event, Austria, May 3-7, 2021_. OpenReview.net, 2021.
  * Zhang et al. (2017) Chiyuan Zhang, Samy Bengio, Moritz Hardt, Benjamin Recht, and Oriol Vinyals.  Understanding deep learning requires rethinking generalization.  In _5th International Conference on Learning Representations, ICLR 2017, Toulon, France, April 24-26, 2017, Conference Track Proceedings_. OpenReview.net, 2017.
  * Zhou and Liang (2018) Yi Zhou and Yingbin Liang.  Critical points of linear neural networks: Analytical forms and landscape properties.  In _Proc. Sixth International Conference on Learning Representations (ICLR)_ , 2018.

##  A Random Matrices

###  A.1 Multiplication of Random Matrices

In this section, we present and prove some statistical properties of general
random matrices and their multiplications. Let
$\\{Q_{n}\in\mathbb{R}^{m_{{n}}\times m_{{n}-1}}\\}_{n=1}^{N}$ denote a set of
random matrix whose elements are sampled iid from a distribution with mean $0$
and variance $\sigma_{n}^{2}$, whose kurtosis is bounded by $c$, and whose
support is compact where the norm of each element is bounded by
$c^{\prime}\sigma_{n}^{2}$ for fixed constants $c,c^{\prime}$. Let

$\displaystyle\bm{Q}^{l}$
$\displaystyle=\prod_{{n}={l}}^{1}Q_{n}=Q_{{l}}\cdot\ldots\cdot Q_{1},$
$\displaystyle{B}^{l}={\bm{Q}^{l}}^{\top}\bm{Q}^{{l}}\in\mathbb{R}^{m_{0}\times
m_{0}}$ (9) $\displaystyle\textsf{Q}^{l}$
$\displaystyle=\prod_{{n}={N}}^{{l}+1}Q_{n}=Q_{N}\cdot\ldots\cdot Q_{{l}+1},$
$\displaystyle{A}^{l}=\textsf{Q}^{{l}}{\textsf{Q}^{{l}}}^{\top}\in\mathbb{R}^{m_{N}\times
m_{N}}$ (10)

###### Theorem 10

For random matrices ${A}^{l}$ and ${B}^{l}$ as defined in (9)-(10)

$\displaystyle\mathbb{E}({B}^{l})$
$\displaystyle=\beta_{l}I\quad\quad\beta_{l}=\prod\limits_{{n}=1}^{{l}}m_{n}\sigma_{n}^{2}$
(11) $\displaystyle\mathbb{E}({A}^{l})$
$\displaystyle=\alpha_{l}I\quad\quad\alpha_{l}=\prod\limits_{{n}={l}+1}^{{N}}m_{{n}-1}\sigma_{n}^{2}$
(12)

Proof  We only prove (11), as the proof of (12) is similar. To simplify the
presentation, we use the following auxiliary notations: $V=Q_{1}$,
$U=\prod_{{n}={l}}^{2}Q_{n}\implies\bm{Q}^{l}=U{V}$.

The proof proceeds by induction on ${l}$.

  * •

${l}=1$:

$\begin{split}\mathbb{E}[{B}^{1}_{ij}]&=\mathbb{E}[\sum_{k=1}^{m_{1}}{V}_{ki}{V}_{kj}]=\sum_{k=1}^{m_{1}}\mathbb{E}[{V}_{ki}]\mathbb{E}[{V}_{kj}]=0\qquad{i\neq
j}\\\
\mathbb{E}[{B}^{1}_{ii}]&=\mathbb{E}[\sum_{k=1}^{m_{1}}{V}_{ki}{V}_{ki}]=\sum_{k=1}^{m_{1}}\mathbb{E}[{V}_{ki}^{2}]=m_{1}\sigma_{1}^{2}\\\
\end{split}$

Thus $\mathbb{E}({B}^{1})=\beta_{1}I$.

  * •

Assume that (11) holds for ${l}-1$.

${B}^{l}_{ij}=\sum_{k}\bm{Q}^{l}_{ki}\bm{Q}^{l}_{kj}=\sum_{k}\sum_{\nu}U_{k\nu}{V}_{\nu
i}\sum_{\rho}U_{k\rho}{V}_{\rho j}\\\ $

and therefore

$\mathbb{E}[{B}^{l}_{ij}]=\sum_{k}\sum_{\nu}\sum_{\rho}\mathbb{E}[U_{k\nu}{V}_{\nu
i}U_{k\rho}{V}_{\rho j}]=\sum_{\nu}\sum_{\rho}\mathbb{E}[{V}_{\nu i}{V}_{\rho
j}]\sum_{k}\mathbb{E}[U_{k\nu}U_{k\rho}]$ (13)

The last transition follows from the independence of $U$ and ${V}$. From (13),
where we denote ${B}^{\prime}=U^{\top}U$

$\begin{split}\mathbb{E}[{B}^{l}_{ij}]&=\sum_{\nu}\mathbb{E}[{V}_{\nu
i}]\sum_{\rho}\mathbb{E}[{V}_{\rho
j}]\mathbb{E}[(U^{\top}U)_{\nu\rho}]=0\qquad{i\neq j}\\\
\mathbb{E}[{B}^{l}_{ii}]&=\sum_{\nu}\sum_{\rho}\mathbb{E}[{V}_{\nu i}{V}_{\rho
i}]\mathbb{E}[(U^{\top}U)_{\nu\rho}]=\sum_{\nu=1}^{m_{1}}\mathbb{E}[{V}_{\nu
i}^{2}]\mathbb{E}[{B}^{\prime}_{\nu\nu}]=m_{1}\sigma_{1}^{2}\prod\limits_{{n}=2}^{{l}}m_{{n}}\sigma_{n}^{2}=\beta_{l}\end{split}$

where the last transition above follows from the induction assumption applied
to ${B}^{\prime}=U^{\top}U$. Thus $\mathbb{E}({B}^{l})=\beta_{l}I$ and (11)
follows.

Let $\mathbbm{m}$ denote the width of the smallest hidden layer,
$\mathbbm{m}=\min\left({m_{1},\ldots,m_{{L}-1}}\right)$, and assume that
$\max\left({m_{1},\ldots,m_{{L}-1}}\right)\leq\mathbbm{m}+\Delta m$ for a
fixed constant $\Delta m$. $m_{0},m_{L}$ are likewise fixed constants ($m_{0}$
corresponds to the input dimension $q$ and $m_{L}$ corresponds to the number
of classes $K$), while $\mathbbm{m}$ is not bounded. Assume that the
distribution of the random matrices $\\{Q_{n}\\}_{{n}=1}^{L}$ is normalized as
specified in Def. 4, where specifically

$\sigma^{2}_{n}=\frac{2}{m_{{n}-1}+m_{n}}~{}\forall{n}\in[2\ldots{L}-1],\qquad\sigma^{2}_{1}=\frac{1}{m_{1}},\qquad\sigma^{2}_{L}=\frac{1}{m_{{L}-1}}$
(14)

It follows that asymptotically, when $\mathbbm{m}\to\infty$, we can write

$\displaystyle m_{n}\sigma_{n}^{2}$
$\displaystyle=1+O\left(\frac{1}{\mathbbm{m}}\right)~{}~{}{n}\in$
$\displaystyle[1$ $\displaystyle\ldots{L}-1]$ $\displaystyle\qquad
m_{{L}}\sigma_{{L}}^{2}=\frac{m_{{L}}}{\mathbbm{m}}+O\left(\frac{1}{\mathbbm{m}}\right)$
$\displaystyle m_{{n}-1}\sigma_{n}^{2}$
$\displaystyle=1+O\left(\frac{1}{\mathbbm{m}}\right)~{}~{}{n}\in$
$\displaystyle[2$ $\displaystyle\ldots{L}]$ $\displaystyle\qquad
m_{0}\sigma_{1}^{2}=\frac{m_{0}}{\mathbbm{m}}+O\left(\frac{1}{\mathbbm{m}}\right)$

###### Corollary 11

From Thm 10, using the initialization scheme specified in Def. 4 and (14)

$\displaystyle\mathbb{E}({B}^{l})$
$\displaystyle=I+O\left(\frac{1}{\mathbbm{m}}\right)~{}~{}\forall{l}\in$
$\displaystyle[1$ $\displaystyle\ldots{L}-1]$
$\displaystyle\qquad\mathbb{E}({B}^{{L}})=\frac{m_{L}}{\mathbbm{m}}I+O\left(\frac{1}{\mathbbm{m}}\right)=O\left(\frac{1}{\mathbbm{m}}\right)$
$\displaystyle\mathbb{E}({A}^{l})$
$\displaystyle=I+O\left(\frac{1}{\mathbbm{m}}\right)~{}~{}\forall{l}\in$
$\displaystyle[1$ $\displaystyle\ldots{L}-1]$
$\displaystyle\qquad\mathbb{E}({A}^{0})=\frac{m_{0}}{\mathbbm{m}}I+O\left(\frac{1}{\mathbbm{m}}\right)=O\left(\frac{1}{\mathbbm{m}}\right)$

Recall that in our asymptotic matrix notations,
$O\left(\frac{1}{\mathbbm{m}}\right)$ is short-hand for a matrix which, for
large enough $\mathbbm{m}$, is upper bounded element-wise by
$C\frac{1}{\mathbbm{m}}$ where $C$ denotes a fixed full rank matrix, and
similarly $O(\mu^{2})$ for small enough $\mu$. When ${B}^{l}$ is concerned
$C\in\mathbb{R}^{m_{0}\times m_{0}}$, and when ${A}^{l}$ is concerned
$C\in\mathbb{R}^{m_{L}\times m_{L}}$. More generally, for clarity and
simplicity of notations, when we discuss asymptotic properties of functions of
a certain matrix $M$, including its expected value $\mathbb{E}(M)$ or variance
$\mathrm{var}(M)$, it is to be understood that the properties are considered
element-wise.

###### Theorem 12

For random matrices ${A}^{l}$ and ${B}^{l}$ as defined in (9)-(10), using the
initialization scheme specified in Def. 4 and (14)

$\mathrm{var}({B}^{l})=O\left(\frac{1}{\mathbbm{m}}\right),\qquad\mathrm{var}({A}^{l})=O\left(\frac{1}{\mathbbm{m}}\right)\quad\forall
l$ (15)

Proof  Once again, we only provide a detailed proof for
$\mathrm{var}({B}^{l})$, as the proof for $\mathrm{var}({A}^{l})$ is similar.
From Corr 11, and since element-wise
$Var(Z)=\mathbb{E}(Z^{2})-[\mathbb{E}(Z)]^{2}$, it is sufficient to show that
the following (somewhat stronger) assertion is valid:

$\mathbb{E}[({B}^{l}_{ij})^{2}]=\begin{cases}O\left(\frac{1}{\mathbbm{m}}\right)&i\neq
j\\\ 1+O\left(\frac{1}{\mathbbm{m}}\right)&i=j,{l}<{L}\\\
O\left(\frac{1}{\mathbbm{m}}\right)&i=j,{l}={L}\end{cases},\qquad\mathbb{E}[{B}^{l}_{ii}{B}^{l}_{jj}]=1+O\left(\frac{1}{\mathbbm{m}}\right)~{}~{}i\neq
j$ (16)

The proof proceeds by induction on ${l}$.

##### Base case ${l}=1$.

Below $Q$ stands for $Q_{1}$, to simplify index notations.

$\mathbb{E}[({B}^{1}_{ij})^{2}]=\mathbb{E}\bigg{[}\sum_{\nu=1}^{m_{1}}Q_{\nu
i}Q_{\nu j}\sum_{\rho=1}^{m_{1}}Q_{\rho i}Q_{\rho
j}\bigg{]}=\begin{cases}\sum\limits_{\nu=1}^{m_{1}}\mathbb{E}[Q_{\nu
i}^{2}]\mathbb{E}[Q_{\nu j}^{2}]=\frac{1}{m_{1}}&i\neq j\\\
\sum\limits_{\nu=1}^{m_{1}}\sum\limits_{{\begin{subarray}{c}\rho=1\\\
\rho\neq\nu\end{subarray}}}^{m_{1}}\mathbb{E}[Q_{\nu i}^{2}]\mathbb{E}[Q_{\rho
i}^{2}]+\sum\limits_{\nu=1}^{m_{1}}\mathbb{E}[Q_{\nu
i}^{4}]=1+O\left(\frac{1}{\mathbbm{m}}\right)&i=j\end{cases}$

Above we use the assumed fixed bound on the kurtosis of the distribution of
$Q$.

$\mathbb{E}[{B}^{l}_{ii}{B}^{l}_{jj}]=\mathbb{E}\bigg{[}\sum_{\nu=1}^{m_{1}}Q_{\nu
i}Q_{\nu i}\sum_{\rho=1}^{m_{1}}Q_{\rho j}Q_{\rho
j}\bigg{]}=\sum\limits_{\nu=1}^{m_{1}}\sum\limits_{\rho=1}^{m_{1}}\frac{1}{m_{1}}\frac{1}{m_{1}}=1$

##### Induction step.

Assume that (16) holds for ${l}-1$, and similarly to the above, let $Q$ stand
for $Q_{l}$ to simplify index notations. Let $\frac{1}{\tilde{m}_{l}}$ denote
the variance of $Q_{l}$ as defined in (14).

$\begin{split}\mathbb{E}[({B}^{{l}}_{ij})^{2}]&=\mathbb{E}\bigg{[}\sum_{\nu=1}^{m_{l}}\sum_{\rho=1}^{m_{l}}Q_{\nu
i}{B}^{{l}-1}_{\nu\rho}Q_{\rho
j}\sum_{\alpha=1}^{m_{l}}\sum_{\beta=1}^{m_{l}}Q_{\alpha
i}{B}^{{l}-1}_{\alpha\beta}Q_{\beta
j}\bigg{]}\quad\quad\mathrm{{\color[rgb]{0.36,0.54,0.66}\definecolor[named]{pgfstrokecolor}{rgb}{0.36,0.54,0.66}(by~{}independence)}}\\\
&=\sum_{\nu=1}^{m_{l}}\sum_{\rho=1}^{m_{l}}\sum_{\alpha=1}^{m_{l}}\sum_{\beta=1}^{m_{l}}~{}\mathbb{E}[Q_{\nu
i}Q_{\rho j}Q_{\alpha i}Q_{\beta
j}]\mathbb{E}[{B}^{{l}-1}_{\nu\rho}{B}^{{l}-1}_{\alpha\beta}]\end{split}$ (17)

When ${i\neq j}$, using the independence of the elements of $Q_{l}$ and the
induction assumption

$\begin{split}\mathbb{E}[({B}^{{l}}_{ij})^{2}]&=\sum_{\begin{subarray}{c}\nu=1\\\
\alpha=\nu\end{subarray}}^{m_{l}}\sum_{\begin{subarray}{c}\rho=1\\\
\beta=\rho\end{subarray}}^{m_{l}}\mathbb{E}[Q_{\nu i}^{2}]\mathbb{E}[Q_{\rho
j}^{2}]\mathbb{E}[({B}^{{l}-1}_{\nu\rho})^{2}]\\\
&=\sum_{\nu=1}^{m_{l}}\sum_{\begin{subarray}{c}\rho=1\\\
\rho\neq\nu\end{subarray}}^{m_{l}}\frac{1}{\tilde{m}_{l}}\frac{1}{\tilde{m}_{l}}\mathbb{E}[({B}^{{l}-1}_{\nu\rho})^{2}]+\sum_{\begin{subarray}{c}\nu=1\\\
\rho=\nu\end{subarray}}^{m_{l}}\frac{1}{\tilde{m}_{l}}\frac{1}{\tilde{m}_{l}}\mathbb{E}[({B}^{{l}-1}_{\nu\nu})^{2}]=O\left(\frac{1}{\mathbbm{m}}\right)\quad\forall
l\end{split}$

When $i=j$, we collect below all the terms in (17) which are not 0:

$\begin{split}E[({B}^{{l}}_{ii})^{2}]&=\sum\limits_{\begin{subarray}{c}\nu=1\\\
\rho=\nu\end{subarray}}^{m_{l}}\sum\limits_{{\begin{subarray}{c}\alpha=1\\\
\beta=\alpha\\\ \alpha\neq\nu\end{subarray}}}^{m_{l}}\mathbb{E}[Q_{\nu
i}^{2}]\mathbb{E}[Q_{\alpha
i}^{2}]\mathbb{E}[{B}^{{l}-1}_{\nu\nu}{B}^{{l}-1}_{\alpha\alpha}]+\sum\limits_{\begin{subarray}{c}\nu=1\\\
\alpha=\nu\end{subarray}}^{m_{l}}\sum\limits_{{\begin{subarray}{c}\rho=1\\\
\beta=\rho\\\ \rho\neq\nu\end{subarray}}}^{m_{l}}\mathbb{E}[Q_{\nu
i}^{2}]\mathbb{E}[Q_{\rho i}^{2}]\mathbb{E}[({B}^{{l}-1}_{\nu\rho})^{2}]\\\
&+\sum\limits_{\begin{subarray}{c}\nu=1\\\
\beta=\nu\end{subarray}}^{m_{l}}\sum\limits_{{\begin{subarray}{c}\rho=1\\\
\alpha=\rho\\\ \rho\neq\nu\end{subarray}}}^{m_{l}}\mathbb{E}[Q_{\nu
i}^{2}]\mathbb{E}[Q_{\rho
i}^{2}]\mathbb{E}[({B}^{{l}-1}_{\nu\rho})^{2}]+\sum\limits_{\begin{subarray}{c}\nu=1\\\
\rho=\alpha=\beta=\nu\end{subarray}}^{m_{l}}\mathbb{E}[Q_{\nu
i}^{4}]\mathbb{E}[({B}^{{l}-1}_{\nu\nu})^{2}]\end{split}$

From the induction assumption and since the kurtosis of $Q$ is bounded by
assumption

$\begin{split}E[({B}^{{l}}_{ii})^{2}]&=\sum\limits_{\nu=1}^{m_{l}}\sum\limits_{{\begin{subarray}{c}\alpha=1\\\
\alpha\neq\nu\end{subarray}}}^{m_{l}}\mathbb{E}[Q_{\nu
i}^{2}]\mathbb{E}[Q_{\alpha
i}^{2}]\mathbb{E}[{B}^{{l}-1}_{\nu\nu}{B}^{{l}-1}_{\alpha\alpha}]+O\left(\frac{1}{\mathbbm{m}}\right)\\\
&=\sum\limits_{\nu=1}^{m_{l}}\sum\limits_{{\begin{subarray}{c}\alpha=1\\\
\alpha\neq\nu\end{subarray}}}^{m_{l}}\frac{1}{\tilde{m}_{l}}\frac{1}{\tilde{m}_{l}}\mathbb{E}[{B}^{{l}-1}_{\nu\nu}{B}^{{l}-1}_{\alpha\alpha}]+O\left(\frac{1}{\mathbbm{m}}\right)=\begin{cases}1+O\left(\frac{1}{\mathbbm{m}}\right)&{l}<{L\vspace{0.2cm}}\\\
O\left(\frac{1}{\mathbbm{m}}\right)&{l}={L}\end{cases}\end{split}$ (18)

To justify the last transition, recall that the initialization scheme
specified in Def. 4 implies that
$m_{l}\frac{1}{\tilde{m}_{l}}=1+O\left(\frac{1}{\mathbbm{m}}\right)~{}\forall{l}\in[1\ldots{L}-1]$.
Making use once again of the induction assumption, it follows from (18) that
now $E[({B}^{{l}}_{ii})^{2}]=1+O\left(\frac{1}{\mathbbm{m}}\right)$. However,
since $m_{L}$ is fixed whereas
$\frac{1}{\tilde{m}_{L}}=O\left(\frac{1}{\mathbbm{m}}\right)$, it follows that
$E[({B}^{{L}}_{ii})^{2}]=O\left(\frac{1}{\mathbbm{m}}\right)$.

A similar argument will show that
$\mathbb{E}[{B}^{l}_{ii}{B}^{l}_{jj}]=1+O\left(\frac{1}{\mathbbm{m}}\right)$
when $i\neq j$.

###### Theorem 13

Let $\\{X(\mathbbm{m})\\}$ denote a sequence of random matrices where
$\mathbb{E}[X(\mathbbm{m})]=F+O\left(\frac{1}{\mathbbm{m}}\right)$ and
$\mathrm{var}[X(\mathbbm{m})]=O\left(\frac{1}{\mathbbm{m}}\right)$. Then
$X(\mathbbm{m})\xrightarrow{p}F$, where $\xrightarrow{p}$ denotes element-wise
convergence in probability as $\mathbbm{m}\to\infty$.

Proof  For every element $(i,j)$ of matrix $X(\mathbbm{m})$, we need to show
that $\forall\varepsilon,\delta>0~{}\exists\mathbbm{m}^{\prime}\in\mathbb{N}$,
such that $\forall\mathbbm{m}>\mathbbm{m}^{\prime}$

$Prob\left(|X_{ij}(m)-F_{ij}|>\varepsilon\right)<\delta$

Henceforth we use $x,f$ as shorthand for $X_{ij}(\mathbbm{m}),F_{ij}$
respectively. Since
$\mathbb{E}[X(\mathbbm{m})]=F+O\left(\frac{1}{\mathbbm{m}}\right)$ where by
definition the asymptotic behavior occurs element-wise, it follows that
$\forall\varepsilon\\!>\\!0$ $\exists\mathbbm{m}_{1}\in\mathbb{N}$ such that
$\forall\mathbbm{m}>\mathbbm{m}_{1}$ we have

$|\mathbb{E}(x)-f|<\frac{\varepsilon}{2}$

In which case

$Prob\left(|x-f|>\varepsilon\right)\leq
Prob\left(|x-\mathbb{E}(x)|>\frac{\varepsilon}{2}\right)$

Since $\mathrm{var}[X(\mathbbm{m})]=O\left(\frac{1}{\mathbbm{m}}\right)$, it
follows that
$\forall\varepsilon,\delta>0,~{}\exists\mathbbm{m}_{2}\in\mathbb{N}~{}~{}\ni~{}\forall\mathbbm{m}>\mathbbm{m}_{2}$

$\mathrm{var}(x)<\frac{\varepsilon^{2}}{4}\delta$

From the above, and using Chebyshev inequality

$Prob\left(|x-f|>\varepsilon\right)<\frac{4\mathrm{var}(x)}{\varepsilon^{2}}<\delta$

$\forall\mathbbm{m}>\mathbbm{m}^{\prime}$, where
$\mathbbm{m}^{\prime}=\max\\{\mathbbm{m}_{1},\mathbbm{m}_{2}\\}$.

###### Corollary 14

Let ${A}^{l}(m)$ and ${B}^{l}(\mathbbm{m})$ denote a sequence of random
matrices as defined in (9)-(10), corresponding to multi-layer linear models
for which $m=\min\left({m_{1},...,m_{L-1}}\right)$. Then

$\displaystyle{B}^{l}(m)\xrightarrow{p}I~{}~{}\forall{l}\in$ $\displaystyle[1$
$\displaystyle\ldots{L}-1]$ $\displaystyle{B}^{{L}}(m)\xrightarrow{p}0$
$\displaystyle{A}^{l}(m)\xrightarrow{p}I~{}~{}\forall{l}\in$ $\displaystyle[1$
$\displaystyle\ldots{L}-1]$ $\displaystyle{A}^{(0)}(m)\xrightarrow{p}0$

Proof  This result follows from Corr 11, Thm 12, and Thm 13.

###  A.2 The Dynamics of Random Matrices

We now formalize a dynamical system, which captures the evolution of the
weight matrices $\\{W_{l}\\}$ by gradient descent as seen in (33)-(34).
Specifically, given random matrices as defined in (9)-(10), consider a
dynamical system whereby $Q_{j}\rightarrow Q_{j}-\Delta Q_{j}~{}\forall j$,
and

$\begin{split}\Delta
Q_{j}&=\mu\Big{(}\prod_{{n}={L}}^{j+1}Q_{n}\Big{)}^{\top}G_{r}\Big{(}\prod_{{n}=j-1}^{1}Q_{n}\Big{)}^{\top}=\mu[\textsf{Q}^{j}]^{\top}G_{r}[\bm{Q}^{j-1}]^{\top}\\\
G_{r}&=\bm{Q}^{{L}}\Sigma_{\scriptscriptstyle XX}-\Sigma_{\scriptscriptstyle
YX}\end{split}$ (19)

Denoting
$\bm{Q}^{l}\rightarrow\bm{Q}^{l}-\Delta\bm{Q}^{l},~{}{B}^{l}\rightarrow{B}^{l}-\Delta{B}^{l}$
and applying the product rule

$\displaystyle\Delta\bm{Q}^{l}$
$\displaystyle=\sum_{j=1}^{{l}}\Big{(}\prod_{{n}={l}}^{j+1}Q_{n}\Big{)}\Delta
Q_{j}\Big{(}\prod_{{n}=j-1}^{1}Q_{n}\Big{)}=\sum_{j=1}^{{l}}\Big{(}\prod_{{n}={l}}^{j+1}Q_{n}\Big{)}\Delta
Q_{j}\bm{Q}^{j-1}$ (20) $\displaystyle\Delta{B}^{l}$
$\displaystyle=[{\Delta\bm{Q}^{l}}^{\top}\bm{Q}^{l}+{\bm{Q}^{l}}^{\top}\Delta\bm{Q}^{l}]$
(21)

Recall that in our notations, $O(\mu^{2})$ is short-hand for a matrix which,
for small enough $\mu$, is upper bounded element-wise by $C\mu^{2}$ where $C$
denotes a _fixed_ full rank matrix. Similarly,
$O\left(\frac{1}{\mathbbm{m}}\right)$ is short-hand for a matrix which, for
large enough $\mathbbm{m}$, is upper bounded element-wise by
$C\frac{1}{\mathbbm{m}}$. Since $m_{0}=m_{q}$ and $m_{L}=m_{K}$, when
${B}^{l}$ is concerned $C\in\mathbb{R}^{q\times q}$, and when ${A}^{l}$ is
concerned $C\in\mathbb{R}^{K\times K}$. In addition, we will use the notation
$O(\varepsilon)$ as short-hand for a matrix that is upper bounded element-wise
by $C\varepsilon$, where $C$ denotes a _fixed_ full rank matrix.

###### Theorem 15

Let
${B}^{l}(\mathbbm{m})={\bm{Q}^{l}(\mathbbm{m})}^{\top}\bm{Q}^{{l}}(\mathbbm{m})$
denote a sequence of random matrices as defined in (9) for
$\\{Q_{n}\\}_{{n}=1}^{L}$, whose dynamics is captured by (19)-(21). Assume
that ${B}^{l}(m)$ is full rank $\forall l,\mathbbm{m}$. If

$\begin{split}&{B}^{l}(m)\xrightarrow{p}I+O(\mu^{2})\qquad\forall{l}\in[1\ldots{L}-1],\qquad\qquad{B}^{{L}}(m)\xrightarrow{p}O(\mu^{2})\\\
&\mathrm{var}[{B}^{l}(m)]=O(\mu^{2})+O\left(\frac{1}{\mathbbm{m}}\right)\qquad\forall{l}\\\
&\mathbb{E}[\bm{Q}^{L}(m)]=O(\mu)+O\left(\frac{1}{\mathbbm{m}}\right)\end{split}$

then

$\displaystyle\Delta{B}^{l}(m)\xrightarrow{p}O(\mu^{2})$
$\displaystyle\qquad\mathrm{var}[\Delta{B}^{l}(m)]=O(\mu^{2})+O\left(\frac{1}{\mathbbm{m}}\right)\qquad\forall{l}$
$\displaystyle\mathbb{E}[\Delta\bm{Q}^{L}(m)]=O(\mu)+O\left(\frac{1}{\mathbbm{m}}\right)$

Proof  ${B}^{l}(m)\xrightarrow{p}I+O(\mu^{2})$ and
${B}^{{L}}(m)\xrightarrow{p}O(\mu^{2})$ implies that
$\forall\varepsilon,\delta>0~{}\exists\hat{\mathbbm{m}}\in\mathbb{N}$, such
that $\forall\mathbbm{m}>\hat{\mathbbm{m}}$ and with probability larger than
$1-\delta$:

${B}^{l}(m)=I+O(\mu^{2})+O(\varepsilon)\quad\forall{l}\leq
L-1,\qquad\qquad{B}^{{L}}(m)=O(\mu^{2})+O(\varepsilon)$ (22)

To evaluate $\Delta{B}^{l}$ from (21), we start from (20) to obtain

${\bm{Q}^{l}}^{\top}\Delta\bm{Q}^{l}=\sum_{j=1}^{{l}}{\bm{Q}^{l}}^{\top}\Big{(}\prod_{{l}}^{j+1}Q_{n}\Big{)}\Delta
Q_{j}\bm{Q}^{j-1}$ (23)

We make use of the Moore–Penrose pseudo-inverse of matrices
$\\{\bm{Q}^{l}\\}$, defined as

$[\bm{Q}^{l}]^{+}=({\bm{Q}^{l}}^{\top}\bm{Q}^{{l}})^{-1}{\bm{Q}^{l}}^{\top}$

This definition is valid since by assumption ${B}^{l}$ is full rank
$\forall{l}$, and is therefore invertible. Additionally

$[\bm{Q}^{l}]^{+}\cdot[(\bm{Q}^{l})^{\top}]^{+}=({\bm{Q}^{l}}^{\top}\bm{Q}^{{l}})^{-1}{\bm{Q}^{l}}^{\top}\cdot{\bm{Q}^{l}}({\bm{Q}^{l}}^{\top}\bm{Q}^{{l}})^{-1}=({\bm{Q}^{l}}^{\top}\bm{Q}^{{l}})^{-1}=({B}^{l})^{-1}$
(24)

We next use (9) and (24) to simplify $T_{j}$ – the $j^{\mathrm{th}}$ term in
the sum (23)

$\displaystyle T_{j}$
$\displaystyle=\mu{\bm{Q}^{l}}^{\top}\left(\prod_{{l}}^{j+1}Q_{n}\right)[\textsf{Q}^{j}]^{\top}G_{r}[\bm{Q}^{j-1}]^{\top}{\bm{Q}^{j-1}}$
$\displaystyle=\mu{\bm{Q}^{l}}^{\top}\left(\prod_{{l}}^{j+1}Q_{n}\right)\cdot\bm{Q}^{j}[\bm{Q}^{j}]^{+}\cdot[(\bm{Q}^{j})^{\top}]^{+}[\bm{Q}^{j}]^{\top}\cdot[\textsf{Q}^{j}]^{\top}G_{r}{B}^{j-1}$
$\displaystyle\mathrm{{\color[rgb]{0.36,0.54,0.66}\definecolor[named]{pgfstrokecolor}{rgb}{0.36,0.54,0.66}({Lemma~{}}\ref{lem:mat-
flip})}}$
$\displaystyle=\mu{\bm{Q}^{l}}^{\top}\bm{Q}^{l}[\bm{Q}^{j}]^{+}[(\bm{Q}^{j})^{\top}]^{+}\bm{Q}^{L}G_{r}{B}^{j-1}$
(25)
$\displaystyle=\mu{B}^{l}[{B}^{j}]^{-1}{\bm{Q}^{{L}}}^{\top}[\bm{Q}^{{L}}\Sigma_{\scriptscriptstyle
XX}-\Sigma_{\scriptscriptstyle YX}]{B}^{j-1}$
$\displaystyle\mathrm{{\color[rgb]{0.36,0.54,0.66}\definecolor[named]{pgfstrokecolor}{rgb}{0.36,0.54,0.66}(using~{}(\ref{eq:deltaQ}),(\ref{eq:pseudo}))}}$
$\displaystyle=\mu{B}^{l}[{B}^{j}]^{-1}[{B}^{{L}}\Sigma_{\scriptscriptstyle
XX}-{\bm{Q}^{{L}}}^{\top}\Sigma_{\scriptscriptstyle YX}]{B}^{j-1}$
$\displaystyle=-\mu{\bm{Q}^{{L}}}^{\top}\Sigma_{\scriptscriptstyle
YX}+O(\mu^{2})+O(\varepsilon)$
$\displaystyle\mathrm{{\color[rgb]{0.36,0.54,0.66}\definecolor[named]{pgfstrokecolor}{rgb}{0.36,0.54,0.66}(using~{}(\ref{eq:Beps}))}}$

By assumption $\mathbb{E}[\bm{Q}^{L}]=O(\mu)+O(\varepsilon)$, and therefore

$\mathbb{E}(T_{j})=O(\mu^{2})+O(\varepsilon)\quad\implies\quad\mathbb{E}[{\bm{Q}^{l}}^{\top}\Delta\bm{Q}^{l}]=\sum_{j=1}^{{l}}\mathbb{E}(T_{j})=O(\mu^{2})+O(\varepsilon)$
(26)

Finally, since
${\Delta\bm{Q}^{l}}^{\top}\bm{Q}^{l}=[{\bm{Q}^{l}}^{\top}\Delta\bm{Q}^{l}]^{\top}$,
from (21) and (26)

$\mathbb{E}[\Delta{B}^{l}]=\mathbb{E}[{\Delta\bm{Q}^{l}}^{\top}\bm{Q}^{l}]+\mathbb{E}[{\Delta\bm{Q}^{l}}^{\top}\bm{Q}^{l}]^{\top}=O(\mu^{2})+O(\varepsilon)$
(27)

To conclude the proof, we will show that
$\forall\varepsilon^{\prime},\delta^{\prime}>0~{}\exists\hat{\mathbbm{m}}^{\prime}\in\mathbb{N}$,
such that $\forall\mathbbm{m}>\hat{\mathbbm{m}}^{\prime}$

$Prob\left(|\Delta{B}^{l}-O(\mu^{2})|>\varepsilon^{\prime}\right)<\delta^{\prime}\qquad\mathrm{element\mbox{-}wise}$
(28)

Let $b$ denote an element of matrix $\Delta{B}^{l}-O(\mu^{2})$. Recall that
(27) is true element-wise with probability $(1-\delta)$
$\forall\varepsilon,\delta$ and $\forall\mathbbm{m}>\hat{\mathbbm{m}}$.
Therefore, $\forall\varepsilon^{\prime},\delta^{\prime}>0$, we can choose
$\varepsilon,\delta$ and the corresponding $\hat{\mathbbm{m}}^{\prime}$ such
that

$Prob\left(|\mathbb{E}(b)|<\frac{\varepsilon^{\prime}}{2}\right)>(1-\frac{\delta^{\prime}}{2})\quad\forall\mathbbm{m}>\hat{\mathbbm{m}}^{\prime}$
(29)
$\begin{split}Prob\left(|b|>\varepsilon^{\prime}\right)&=Prob\left(|b|>\varepsilon^{\prime},~{}|\mathbb{E}(b)|<\frac{\varepsilon^{\prime}}{2}\right)+Prob\left(|b|>\varepsilon^{\prime},~{}|\mathbb{E}(b)|\geq\frac{\varepsilon^{\prime}}{2}\right)\\\
&\leq
Prob\left(|b-\mathbb{E}(b)|>\frac{\varepsilon^{\prime}}{2}\right)+Prob\left(|\mathbb{E}(b)|\geq\frac{\varepsilon^{\prime}}{2}\right)\end{split}$

$\mathrm{var}(\Delta{B}^{l})=O(\mu^{2})+O\left(\frac{1}{\mathbbm{m}}\right)$
implies that
$\forall\varepsilon^{\prime},\delta^{\prime}>0~{}\exists\hat{\mathbbm{m}}^{\prime\prime}\in\mathbb{N}$,
such that $\forall\mathbbm{m}>\hat{\mathbbm{m}}^{\prime\prime}$

$\mathrm{var}(b)<\frac{\delta^{\prime}\varepsilon^{\prime 2}}{8}$ (30)

Using Chebychev inequality, (29) and (30), we finally have

$Prob\left(|b|>\varepsilon^{\prime}\right)\leq\frac{4\mathrm{var}(b)}{\varepsilon^{\prime
2}}+\frac{\delta^{\prime}}{2}<\delta^{\prime}\qquad\forall\mathbbm{m}>\max\\{\hat{\mathbbm{m}}^{\prime},\hat{\mathbbm{m}}^{\prime\prime}\\}$

We can now conclude that (28) is true element-wise.

We will only provide a sketch of the proof for
$\mathrm{var}[\Delta{B}^{l}(m)]$ and $\mathbb{E}[\Delta\bm{Q}^{L}(m)]$, as a
detailed proof follows very similar steps, and principles, to the proof
presented above. To analyze the variance of $\Delta{B}^{l}(m)$, we start from
(25) and observe that
$\Delta{B}^{l}(m)=-{l}\mu{\bm{Q}^{{L}}}^{\top}\Sigma_{\scriptscriptstyle
YX}+O(\mu^{2})+O(\varepsilon)$, from which (and the theorem’s assumptions) it
can be shown that
$\mathrm{var}[\Delta{B}^{l}(m)]=O(\mu^{2})+O\left(\frac{1}{\mathbbm{m}}\right)$.
Similarly to (25), we can derive that $\mathbb{E}[\Delta\bm{Q}^{l}(m)]=-\mu
LG_{r}+O(\mu^{2})+O\left(\frac{1}{\mathbbm{m}}\right)$, from which it follows
that
$\mathbb{E}[\Delta\bm{Q}^{L}(m)]=O(\mu)+O\left(\frac{1}{\mathbbm{m}}\right)$.

###### Theorem 16

Let
${A}^{l}(\mathbbm{m})=\textsf{Q}^{{l}}(\mathbbm{m}){\textsf{Q}^{l}(\mathbbm{m})}^{\top}$
denote a sequence of random matrices as defined in (10) but where
${A}^{(0)}(m)=\textsf{Q}^{0}\Sigma_{\scriptscriptstyle
XX}{\textsf{Q}^{0}}^{\top}$, whose dynamics is captured by (19). Assume that
${A}^{l}(m)$ is full rank $\forall l,\mathbbm{m}$. If

$\begin{split}&{A}^{l}(m)\xrightarrow{p}I+O(\mu^{2})~{}~{}\forall{l}\in[1\ldots{L}-1],\qquad{A}^{0}(\mathbbm{m})\xrightarrow{p}O(\mu^{2})\\\
&\mathrm{var}[{A}^{l}(m)]=O(\mu^{2})+O\left(\frac{1}{\mathbbm{m}}\right)\qquad\forall{l}\\\
&\textsf{Q}^{0}(m)\xrightarrow{p}O(\mu)+O\left(\frac{1}{\mathbbm{m}}\right)\qquad\mathrm{var}[\textsf{Q}^{0}(m)]=O(\mu^{2})+O\left(\frac{1}{\mathbbm{m}}\right)\end{split}$

then

$\displaystyle\Delta{A}^{l}(m)\xrightarrow{p}O(\mu^{2})$
$\displaystyle\qquad\mathrm{var}[\Delta{A}^{l}(m)]=O(\mu^{2})+O\left(\frac{1}{\mathbbm{m}}\right)\qquad\forall{l}$
$\displaystyle\Delta\textsf{Q}^{0}(m)\xrightarrow{p}O(\mu)+O\left(\frac{1}{\mathbbm{m}}\right)$
$\displaystyle\qquad\mathrm{var}[\Delta\textsf{Q}^{0}(m)]=O(\mu^{2})+O\left(\frac{1}{\mathbbm{m}}\right)$

The proof is mostly similar to Thm 15.

###  A.3 Some Useful Lemmas

###### Lemma 17

Given function $G(W)=\frac{1}{2}\|{U}W{V}X-Y\|_{F}^{2}$, its derivative is the
following

$\frac{dG(W)}{dW}={U}^{\top}{U}W{V}X({V}X)^{\top}-{U}^{\top}Y({V}X)^{\top}={U}^{\top}[{U}W{V}\Sigma_{\scriptscriptstyle
XX}-\Sigma_{\scriptscriptstyle YX}]{V}^{\top}$ (31)

###### Lemma 18

Assume $\bm{Q}=\prod_{{n}={N}}^{1}Q_{n}$, where
$Q_{n}\in\mathbb{R}^{m_{{n}}\times m_{{n}-1}}$ denotes a random matrix whose
elements are sampled iid from a distribution with mean $0$ and variance
$\sigma_{n}^{2}$ $\forall i,j$, initialized as defined in (14). Then

$\mathbb{E}[\bm{Q}_{ij}]=0\qquad\mathrm{var}[\bm{Q}_{ij}]=O\left(\frac{1}{\mathbbm{m}}\right)$
(32)

Proof  By induction on ${N}$. Clearly for ${N}=1$:

$\mathbb{E}[\bm{Q}_{ij}]=\mathbb{E}[(Q_{1})_{ij}]=0\qquad\mathrm{var}[\bm{Q}_{ij}]=\mathrm{var}[(Q_{1})_{ij}]=\sigma_{1}^{2}$

Assume that (32) holds for ${N}-1$. Let $V=\prod_{{n}={N}-1}^{1}Q_{n}$,
$U=Q_{N}$. It follows that

$\mathbb{E}[\bm{Q}_{ij}]=\mathbb{E}[(UV)_{ij}]=\sum_{k}\mathbb{E}[U_{ik}V_{kj}]=\sum_{k}\mathbb{E}[U_{ik}]\mathbb{E}[V_{kj}]=0$

where the last transition follows from the independence of $U$ and ${V}$. In a
similar manner

$\begin{split}\mathrm{var}[\bm{Q}_{ij}]&=\mathbb{E}[\bm{Q}_{ij}^{2}]=\mathbb{E}[(\sum_{k}U_{ik}V_{kj})^{2}]=\mathbb{E}[\sum_{k}U_{ik}V_{kj}\sum_{l}U_{il}V_{lj}]=\sum_{k}\mathbb{E}[U_{ik}^{2}]\mathbb{E}[V_{kj}]^{2}\
\\\
&=m_{{N}-1}\sigma_{{N}}^{2}\frac{1}{m_{{N}-1}}\prod_{{n}=1}^{{N}-1}m_{n}\cdot\sigma^{2}_{n}=\frac{1}{m_{{N}}}\prod_{{n}=1}^{{N}}m_{n}\cdot\sigma^{2}_{n}\end{split}$

With the initialization scheme defined in (14),
$\mathrm{var}(\bm{Q}_{ij})=O\left(\frac{1}{\mathbbm{m}}\right)$.

###### Lemma 19

Let $C\in\mathbb{R}^{k\times m}$ and $V\in\mathbb{R}^{m\times k}$ both of rank
$k$, where $k<m$. Then

$CV=I~{}\implies~{}C=CVV^{+}$

Proof  By definition $C=V^{+}$ , hence

$C=V^{+}VC=CVV^{+}$

##  B Supplementary Proofs

###  B.1 Deep Linear networks

Here we prove Prop. 5 as defined in Section 3.2.

Proposition 5. Let $G_{r}^{(t)}$ (Def. 3) denote the gradient matrix at time
${t}$. Let $A_{{l}}^{(t)}$ and $B_{{l}}^{(t)}$ denote the gradient scale
matrices, which are defined in (2). The compact representation
${\bm{W}}^{(t)}$ obeys the following dynamics:

${\bm{W}}^{(t+1)}={\bm{W}}^{(t)}-\mu\sum_{{l}=1}^{L}{A}_{{l}}^{(t)}\cdot
G_{r}^{(t)}\cdot{B}_{{l}-1}^{(t)}+O(\mu^{2})$

Proof  At time ${t}$, the gradient step $\Delta W_{l}^{(t)}$ of layer $l$ is
defined by differentiating $L(\mathbb{X})$ with respect to $W_{l}^{(t)}$.
Henceforth we omit index ${t}$ for clarity. First, we rewrite $L(\mathbb{X})$
as follows:

$L(\mathbb{X};W_{l})=\frac{1}{2}\|\left(\prod_{j=L}^{l+1}W_{j}\right)W_{l}\left(\prod_{j=l-1}^{1}W_{j}\right)X-Y\|_{F}^{2}$

Differentiating $L(\mathbb{X};W_{l})$ to obtain the gradient $\Delta
W_{l}=\frac{\partial L(\mathbb{X};W_{l})}{\partial W_{l}}$, using Lemma 17
above, we get

$\Delta
W_{l}=\left(\prod_{j=L}^{l+1}W_{j}\right)^{\top}[{\bm{W}}\Sigma_{\scriptscriptstyle
XX}-\Sigma_{\scriptscriptstyle YX}]\left(\prod_{l-1}^{1}W_{j}\right)^{\top}$
(33)

Finally

$\Delta{\bm{W}}=\prod_{l=L}^{1}(W_{l}-\mu\Delta
W_{l})-\prod_{l=L}^{1}W_{l}=-\mu\sum_{{l}=1}^{L}\left(\prod_{{n}=L}^{{l}+1}W_{n}\right)\Delta
W_{l}\left(\prod_{{n}={l}-1}^{1}W_{n}\right)+O(\mu^{2})$ (34)

Substituting $\Delta W_{l}$ and $G_{r}$ (as specified in Def. 3) into the
above completes the proof.

###  B.2 Weight Evolution

The proofs in this section assume that the norm of the data matrix $X$ is
bounded, that the norms of the initial random weight matrices $W^{(0)}_{l}$
are bounded $\forall l,\mathbbm{m}$, and that the norms of the gradient scale
matrices are also bounded $\forall l,\mathbbm{m}$ in the relevant time range
${t}\leq\bar{t}$. The last assertion follows from our initial assumption, that
the support of the distribution of the weight matrices is compact, and
additionally that the norm of each element is bounded by $c\sigma_{n}^{2}$ for
fixed a constant $c$.999These assumptions can be relaxed to having bounds with
high probability, which follow from the analysis of random matrices in
Appendix A, but for simplicity we assume fix bounds.

We start by proving Thm 6, which is stated in Section 3.2.
Theorem 6. Let ${\bm{W}}$ denote the compact representation of a deep linear
network, where $\mathbbm{m}=\min\left({m_{1},...,m_{L-1}}\right)$ denotes the
size of its smallest hidden layer and $\mathbbm{m}\geq\\{m_{0},m_{L}\\}$. At
each layer $l$, assume weight initialization $W^{(0)}_{l}$ obtained by
sampling from a distribution with mean $0$ and variance $\sigma_{l}^{2}$,
normalized as specified in Def. 4. Let
$\left({B}^{(t)}_{l}(\mathbbm{m})\right)_{\mathbbm{m}=1}^{\infty}$ and
$\left({A}^{(t)}_{l}(\mathbbm{m})\right)_{\mathbbm{m}=1}^{\infty}$ denote two
sequences of gradient scale matrices as defined in (2), where the
$\mathbbm{m}^{th}$ element of each series corresponds to a network whose
smallest hidden layer has $\mathbbm{m}$ neurons. Let $\xrightarrow{p}$ denote
element-wise convergence in probability as $\mathbbm{m}\to\ \infty$. Then
$\forall{t},l$:

$\begin{split}&{B}^{(t)}_{l}(m)\xrightarrow{p}I+O(\mu^{2})~{}~{}\forall{l}\in[1\ldots{L}-1],\qquad{B}^{(t)}_{{L}}(m)\xrightarrow{p}\frac{m_{L}}{\mathbbm{m}}I+O(\mu^{2})\xrightarrow{p}O(\mu^{2})\\\
&\mathrm{var}[{B}^{(t)}_{l}(m)]=O(\mu^{2})+O\left(\frac{1}{\mathbbm{m}}\right)\qquad\forall{l}\end{split}$

and

$\begin{split}&{A}^{(t)}_{l}(m)\xrightarrow{p}I+O(\mu^{2})~{}~{}\forall{l}\in[1\ldots{L}-1],\qquad{A}^{(t)}_{0}(\mathbbm{m})\xrightarrow{p}\frac{m_{0}}{\mathbbm{m}}I+O(\mu^{2})\xrightarrow{p}O(\mu^{2})\\\
&\mathrm{var}[{A}^{(t)}_{l}(m)]=O(\mu^{2})+O\left(\frac{1}{\mathbbm{m}}\right)\qquad\forall{l}\end{split}$

Proof  We will only work out the detailed proof of the first assertion
concerning ${B}^{(t)}_{l}(m)$, as the proof of the second assertion is
similar. Specifically, we prove by induction on time ${t}$ a stronger claim:

$\begin{split}&{B}^{(t)}_{l}(m)\xrightarrow{p}I+O(\mu^{2})~{}~{}\forall{l}\in[1\ldots{L}-1],\qquad{B}^{(t)}_{{L}}(m)\xrightarrow{p}O(\mu^{2})\\\
&\mathrm{var}[{B}^{(t)}_{l}(m)]=O(\mu^{2})+O\left(\frac{1}{\mathbbm{m}}\right)\qquad\forall{l}\\\
&\mathbb{E}[{{\bm{W}}}^{(t)}(m)]=O(\mu)+O\left(\frac{1}{\mathbbm{m}}\right)\end{split}$
(35)

For the corresponding series $\\{{\bm{W}}(m)\\}$.

##### Base case $t=0$.

From (2)

${B}_{{l}}^{(0)}\coloneqq\Big{(}\prod_{j=l}^{1}W_{j}^{(0)}\Big{)}^{\top}\Big{(}\prod_{j=l}^{1}W_{j}^{(0)}\Big{)}$

Recall that all weight matrices $\\{W^{(0)}_{l}\\}_{{l}=1}^{L}$ are
initialized by sampling from a distribution with mean $0$ and variance
$\sigma_{l}^{2}=O(\frac{1}{\mathbbm{m}})$. In Appendix A.1 we prove some
statistical properties of such matrices, and their corresponding _gradient
scale matrices_ $A_{l}^{(0)}$ and $B_{l}^{(0)}$. This analysis culminates in
Corr 11 and Thm 12, which can now be used directly to infer that

$\begin{split}&\mathbb{E}({B}^{(0)}_{l})=I+O\left(\frac{1}{\mathbbm{m}}\right)~{}\forall{l}\in[1\ldots{L}-1],\qquad\mathbb{E}({B}^{(0)}_{{L}})=O\left(\frac{1}{\mathbbm{m}}\right)\\\
&\mathrm{var}[{B}_{{l}}^{(0)}(\mathbbm{m})]=O\left(\frac{1}{\mathbbm{m}}\right)\forall
l\end{split}$

Finally, Thm 13 further proves that such matrices also satisfy

${B}^{(0)}_{l}(m)\xrightarrow{p}I~{}~{}\forall{l}\in[1\ldots{L}-1],\qquad{B}^{(0)}_{{L}}(m)\xrightarrow{p}0$

Similarly, from Lemma 18 we have that
$\mathbb{E}[{{\bm{W}}}^{(0)}(m)]=0~{}\forall m$, and we therefore conclude
that the assertion in (35) is true at $t=0$.

##### Induction step.

In Appendix A.2 we analyze random matrices which are defined similarly to the
gradient scale matrices, and whose dynamics correspond to the update rule
defined in (3). The main result is stated in Thm 15, which can be used
directly now to show that if assertion (35) holds for $B_{{l}}^{(t)}$ and
${{\bm{W}}}^{(t)}(m)$, it also holds for $B_{{l}}^{(t+1)}$ and
${{\bm{W}}}^{(t+1)}(m)$.

First, let us verify that the conditions of Thm 15 are met. By construction,
the dynamics captured by (19)-(21) describes the dynamics of $\\{W_{l}\\}$, a
result that follows from the proof of Prop. 5, where specifically (33)-(34)
imply (19)-(20). In addition, as by assumption $\mathbbm{m}\geq m_{0}$ and
$m_{0}=q$, ${B}^{(t)}_{l}(m)\in\mathbb{R}^{q\times q}$ is full rank $\forall
l,\mathbbm{m}$ with probability 1. We may therefore use Thm 15 to conclude,
based on the induction assumption, that

$\begin{split}&\Delta{B}^{(t)}_{l}(m)\xrightarrow{p}O(\mu^{2}),\quad\mathrm{var}[\Delta{B}^{(t)}_{l}(m)]=O(\mu^{2})+O\left(\frac{1}{\mathbbm{m}}\right)\qquad\forall{l}\\\
&\mathbb{E}[\Delta{\bm{W}}^{(t)}(m)]=O(\mu)+O\left(\frac{1}{\mathbbm{m}}\right)\end{split}$

Noting that
${B}^{(t+1)}_{l}(m)={B}^{(t)}_{l}(m)+\Delta{B}^{(t)}(m)~{}\forall{l}$ and
${\bm{W}}^{(t+1)}(m)={\bm{W}}^{(t)}(m)+\Delta{\bm{W}}^{(t)}_{l}(m)$, the
assertion in (35) follows for all $t$.

We proceed to prove Thm 7, which is stated in Section 3.3.
Theorem 7. Let $\bm{w}_{j}^{(t)}$ denote the $j^{\mathrm{th}}$ column of the
compact representation matrix ${\bm{W}}^{(t)}$, and $\bm{w}_{j}^{opt}$ the
$j^{\mathrm{th}}$ column of the optimal solution of (1). Assume that the data
is rotated to its principal coordinate system, and let $d_{j}$ denote the
$j^{\mathrm{th}}$ singular value of the data. Then there exists $\hat{t}$ such
that $\forall\delta,\varepsilon$ and $\forall{t}\leq\hat{t}$,
$\exists\hat{\mathbbm{m}},\hat{\mu}$ such that
$\forall\mu<\hat{\mu},\mathbbm{m}\geq\hat{\mathbbm{m}}$

$Prob\bigg{(}\Big{\|}\bm{w}_{j}^{(t+1)}-\big{[}\lambda_{j}^{t}\bm{w}_{j}^{(0)}+[1-\lambda_{j}^{t}]\bm{w}_{j}^{opt}\big{]}\Big{\|}<\varepsilon\bigg{)}>(1-\delta)\quad\quad\lambda_{j}=1-\mu
d_{j}L$

Proof  First let us derive a bound on the total gradient magnitude
$\|G_{r}^{(t)}\|$. Since by assumption the norm of the data $\|X\|$ and the
norms of $W^{(0)}_{l}~{}\forall l$ are bounded, it follows that the loss at
time $t=0$ is also bounded. Let $U_{1}$ denote this bound, and $U_{2}$ denote
the data bound $\|X\|^{2}\leq U_{2}$. Let $L^{(t)}(\mathbb{X})$ denote the
loss at time ${t}$. Since the loss is decreasing, it follows that
$L^{(t)}(\mathbb{X})\leq U_{1}~{}\forall{t}$, and therefore

$\|G_{r}^{(t)}\|^{2}=\|[{\bm{W}}^{(t)}X-Y]X^{\top}\|^{2}\leq
2L^{(t)}(\mathbb{X})\|X\|^{2}\leq 2U_{1}U_{2}$

We now introduce the notation ${B}^{(t)}_{l}=I+O(\mu^{2})+\Delta^{(t)}_{B{l}}$
and ${A}^{(t)}_{l}=I+O(\mu^{2})+\Delta^{(t)}_{A{l}}$. Let $U_{3}$ denote a
tight bound so that
$\\{\|\Delta^{(t)}_{A{l}}\|^{2},\|\Delta^{(t)}_{B{l}}\|^{2}\\}\leq
U_{3}~{}\forall{l},{t}\leq\bar{t}$, where we further assume that $U_{3}\leq 1$
(this last assumption will be justified later). Starting from (3)

$\displaystyle{\bm{W}}^{(t+1)}=$
$\displaystyle{\bm{W}}^{(t)}-\mu\sum_{{l}=1}^{L}A_{{l}}^{(t)}\cdot
G_{r}^{(t)}\cdot B_{{l}-1}^{(t)}+O(\mu^{2})$ $\displaystyle=$
$\displaystyle{\bm{W}}^{(t)}-\mu\sum_{{l}=1}^{L}[I+O(\mu^{2})+\Delta^{(t)}_{A{l}}]~{}G_{r}^{(t)}~{}[I+O(\mu^{2})+\Delta^{(t)}_{B{l}}]+O(\mu^{2})$
(36) $\displaystyle=$
$\displaystyle{\bm{W}}^{(t)}-\mu\left[LG_{r}^{(t)}+\bigg{(}\sum_{{l}=1}^{L}\Delta^{(t)}_{A{l}}\bigg{)}G_{r}^{(t)}+G_{r}^{(t)}\sum_{{l}=1}^{L}\Delta^{(t)}_{B{l}}+\sum_{{l}=1}^{L}\Delta^{(t)}_{A{l}}~{}G_{r}^{(t)}~{}\Delta^{(t)}_{B{l}}\right]+O(\mu^{2})$

The last transition is valid because the norms of
$G_{r}^{(t)},\Delta^{(t)}_{A{l}},\Delta^{(t)}_{B{l}}$ are bounded. It follows
that

$\|{\bm{W}}^{(t+1)}-[{\bm{W}}^{(t)}-\mu LG_{r}^{(t)}]\|^{2}\leq\mu
L2U_{1}U_{2}[U_{3}+U_{3}+U_{3}^{2}]+O(\mu^{2})\leq\mu
6LU_{1}U_{2}U_{3}+O(\mu^{2})$

Thus $\exists\hat{\mu}$ such that $\forall\mu<\hat{\mu}$

$\|{\bm{W}}^{(t+1)}-[{\bm{W}}^{(t)}-\mu
LG_{r}^{(t)}]\|^{2}\leq\hat{\mu}12LU_{1}U_{2}U_{3}$ (37)

From Thm 6, $\forall{t}$ and $\forall l\leq{L}-1$,
${B}^{(t)}_{l}(m)\xrightarrow{p}I+O(\mu^{2})$ and
${A}^{(t)}_{l}(m)\xrightarrow{p}I+O(\mu^{2})$ . In other words, if we fix
$\hat{t}$ and consider all the iterations leading to $\hat{t}$, then
$\forall\varepsilon,\delta>0~{}\exists\hat{\mathbbm{m}}\in\mathbb{N}$, such
that $\forall l,\mathbbm{m}>\hat{\mathbbm{m}},{t}\leq\hat{t}$

$\begin{split}&Prob\left(\Big{|}{B}^{(t)}_{l}(m)-[I+O(\mu^{2})]\Big{|}^{2}\leq\min\bigg{\\{}\frac{\varepsilon}{\hat{\mu}12LU_{1}U_{2}},1\bigg{\\}}\right)>(1-\delta)\qquad\mathrm{element\mbox{-}wise}\\\
&Prob\left(\Big{|}{A}^{(t)}_{l}(m)-[I+O(\mu^{2})]\Big{|}^{2}\leq\min\bigg{\\{}\frac{\varepsilon}{\hat{\mu}12LU_{1}U_{2}},1\bigg{\\}}\right)>(1-\delta)\qquad\mathrm{element\mbox{-}wise}\\\
&\implies\qquad
Prob\left(U_{3}\leq\min\bigg{\\{}\frac{\varepsilon}{\hat{\mu}12LU_{1}U_{2}},1\bigg{\\}}\right)>(1-\delta)\quad\mathrm{{\color[rgb]{0.36,0.54,0.66}\definecolor[named]{pgfstrokecolor}{rgb}{0.36,0.54,0.66}(\mathit{U}_{3}~{}is~{}tight)}}\end{split}$

Finally, from (37) and $\forall\mu<\hat{\mu}$

$Prob\left(\|{\bm{W}}^{(t+1)}-[{\bm{W}}^{(t)}-\mu
LG_{r}^{(t)}]\|^{2}\leq\varepsilon\right)>(1-\delta)$ (38)

Next, we evaluate

$\tilde{\bm{W}}^{(t+1)}={\bm{W}}^{(t)}-\mu LG_{r}^{(t)}$ (39)

We first shift to the principal coordinate system defined in Def 1. In this
representation $G_{r}^{(t)}=W^{(t)}D-M$, where $D=diag(\\{d_{j}\\}_{j=1}^{q})$
is a diagonal matrix whose elements are the principal eigenvalues of the data
$\\{d_{j}\\}_{j=1}^{q}$, arranged in decreasing order. Since $D$ is diagonal,
(39) can be written separately for each column of $\tilde{\bm{W}}^{(t)}$, and
we get for each column $\bm{w}_{j}^{(t)}$

$\tilde{\bm{w}}_{j}^{(t+1)}=\bm{w}_{j}^{(t)}-\mu
L[\bm{w}_{j}^{(t)}d_{j}-\bm{m}_{j}]\qquad j\in K$

This is a telescoping series, whose solution is

$\begin{split}\tilde{\bm{w}}_{j}^{(t+1)}&=(1-\mu Ld_{j})\bm{w}_{j}^{(t)}+\mu
L\bm{m}_{j}=\ldots\\\ &=(1-\mu Ld_{j})^{t}\bm{w}_{j}^{(0)}+\mu
Ld_{j}\left[\sum_{\nu=1}^{{t}}(1-\mu
Ld_{j})^{\nu-1}\right]\frac{\bm{m}_{j}}{d_{j}}\\\
&=\lambda_{j}^{t}\bm{w}_{j}^{(0)}+[1-\lambda_{j}^{t}]\frac{\bm{m}_{j}}{d_{j}},\qquad\lambda_{j}=1-\mu
Ld_{j}\end{split}$

As $\bm{w}_{j}^{opt}=\frac{\bm{m}_{j}}{d_{j}}$ and using (38), the assertion
in the theorem follows.

We conclude by proving Thm 8. To this end we introduce another notation -
$A^{(t)}=\sum_{{l}=1}^{L}A_{{l}}^{(t)}$, and let $U_{4}$ denote the bound on
${A}_{l}$ where $\|A_{{l}}^{(t)}\|^{2}\leq
U_{4}~{}~{}\forall{l},{t}\leq\bar{t}$.
Theorem 8. Let $\bm{w}_{j}^{(t)}$ denote the $j^{\mathrm{th}}$ column of the
compact representation matrix ${\bm{W}}^{(t)}$, and $\bm{w}_{j}^{opt}$ the
$j^{\mathrm{th}}$ column of the optimal solution of (1). Assume that the data
is rotated to its principal coordinate system, and let $d_{j}$ denote the
$j^{\mathrm{th}}$ singular value of the data. Then there exists $\breve{t}$
such that $\forall\delta,\varepsilon$ and $\forall{t}\leq\breve{t}$,
$\exists\breve{\mathbbm{m}},\breve{\mu}$ such that
$\forall\mu<\breve{\mu},\mathbbm{m}\geq\breve{\mathbbm{m}}$

$\begin{split}Prob\Bigg{(}\bigg{\|}\bm{w}_{j}^{(t+1)}-&\bigg{[}\prod_{{t^{\prime}}=1}^{{t}}(I-\mu
d_{j}A^{({t^{\prime}})})\bm{w}_{j}^{(0)}+\\\ &\mu
d_{j}\Big{(}\sum_{{t^{\prime}}=1}^{{t}}\prod_{{t^{\prime\prime}}={t^{\prime}}+1}^{{t}}(I-\mu
d_{j}A^{({t^{\prime\prime}})})A^{({t^{\prime}})}\Big{)}\bm{w}_{j}^{opt}\bigg{]}\bigg{\|}<\varepsilon\Bigg{)}>(1-\delta)\end{split}$
(40)

Proof  We use the same bound notations as in the proof of Thm 7, but where
only $\|\Delta^{(t)}_{B{l}}\|^{2}\leq U_{3}~{}\forall{l},{t}\leq\bar{t}$ a
tight bound, $U_{3}\leq 1$. Similarly to (36),

$\begin{split}&{\bm{W}}^{(t+1)}={\bm{W}}^{(t)}-\mu\left[\bigg{(}\sum_{{l}=1}^{L}{A}^{(t)}_{l}\bigg{)}G_{r}^{(t)}+\sum_{{l}=1}^{L}{A}^{(t)}_{l}~{}G_{r}^{(t)}~{}\Delta^{(t)}_{B{l}}\right]+O(\mu^{2})\\\
\implies\qquad&\|{\bm{W}}^{(t+1)}-[{\bm{W}}^{(t)}-\mu
A^{(t)}G_{r}^{(t)}]\|^{2}\leq\mu 2LU_{1}U_{2}U_{4}+O(\mu^{2})\end{split}$

From Thm 6,
$\forall\varepsilon,\delta>0~{}\exists\breve{\mathbbm{m}},\breve{t},\breve{\mu}$,
such that
$\forall\mathbbm{m}\geq\breve{\mathbbm{m}},{t}\leq\breve{t},\mu\leq\breve{\mu}$

$\displaystyle
Prob\left(U_{3}\leq\min\bigg{\\{}\frac{\varepsilon}{\breve{\mu}4LU_{1}U_{2}U_{4}},1\bigg{\\}}\right)>(1-\delta)$
$\displaystyle\implies\quad$ $\displaystyle
Prob\left(\|{\bm{W}}^{(t+1)}-[{\bm{W}}^{(t)}-\mu
A^{(t)}G_{r}^{(t)}]\|^{2}\leq\varepsilon\right)>(1-\delta)$ (41)

As before, we evaluate

$\tilde{\bm{W}}^{(t+1)}={\bm{W}}^{(t)}-\mu A^{(t)}G_{r}^{(t)}$

In the principal coordinate system this expression is separable by columns,
thus

$\tilde{\bm{w}}_{j}^{(t+1)}=\bm{w}_{j}^{(t)}-\mu\sum_{{l}=1}^{L}A_{{l}}^{(t)}(d_{j}\bm{w}_{j}^{(t)}-\bm{m}_{j})\qquad
j\in[K]$

Once again we have a telescoping series, whose solution is

$\tilde{\bm{w}}_{j}^{(t+1)}=\prod_{{t^{\prime}}=1}^{{t}}\big{(}I-\mu
d_{j}A^{({t^{\prime}})}\big{)}\bm{w}_{j}^{(0)}+\mu
d_{j}\Big{[}\sum_{{t^{\prime}}=1}^{{t}}\prod_{{t^{\prime\prime}}={t^{\prime}}+1}^{{t}}\big{(}I-\mu
d_{j}A^{({t^{\prime\prime}})}\big{)}A^{({t^{\prime}})}\Big{]}\frac{\bm{m}_{j}}{d_{j}}$
(42)

As $\bm{w}_{j}^{opt}=\frac{\bm{m}_{j}}{d_{j}}$ and using (41), the assertion
in the theorem follows.

###  B.3 Adding Non-Linear ReLU Activation

In this section, we analyze the two-layer model with ReLU activation, where
only the weights of the first layer are being learned (Arora et al., 2019b).
Similarly to (1), the loss is defined as

$L(\mathbb{X})=\frac{1}{2}\sum_{i=1}^{n}\|f(\bm{x}_{i})-\bm{y}_{i}\|^{2}\qquad
f(\bm{x}_{i})=\bm{a}^{\intercal}\cdot\sigma(W\bm{x}_{i}),~{}~{}\bm{a}\in\mathbb{R}^{\mathbbm{m}},~{}W\in\mathbb{R}^{\mathbbm{m}\times
d}$

$\mathbbm{m}$ denotes the number of neurons in the hidden layer. We consider a
binary classification problem with 2 classes, where $y_{i}=1$ for
$\bm{x}_{i}\in C_{1}$, and $y_{i}=-1$ for $\bm{x}_{i}\in C_{2}$.
$\sigma(\cdot)$ denotes the ReLU activation function applied element-wise to
vectors, where $\sigma(u)=u$ if $u\geq 0$, and $0$ otherwise.

At time ${t}$, each gradient step is defined by differentiating
$L(\mathbb{X})$ with respect to $W$. Due to the non-linear nature of the
activation function $\sigma(\cdot)$, we separately101010Since the ReLU
function is not everywhere differentiable, the following may be considered the
definition of the update rule. differentiate each row of $W$, denoted
$\bm{w}_{r}$ where $r\in[\mathbbm{m}]$, as follows:

$\begin{split}\bm{w}^{(t+1)}_{r}&-\bm{w}^{(t)}_{r}=-\mu\frac{\partial
L(\mathbb{X})}{\partial\bm{w}_{r}}{\bigg{|}}_{{\textstyle\mathstrut}{\bm{w}_{r}=\bm{w}^{(t)}_{r}}}=-\mu\sum_{i=1}^{n}\Big{[}\bm{a}^{\intercal}\cdot\sigma(W^{(t)}\bm{x}_{i})-y_{i}\Big{]}\frac{\partial
f(\bm{x}_{i})}{\partial\bm{w}_{r}}{\bigg{|}}_{{\textstyle\mathstrut}\bm{w}_{r}=\bm{w}^{(t)}_{r}}\\\
&=-\mu\sum_{i=1}^{n}\Big{[}\sum_{j=1}^{\mathbbm{m}}a_{j}\sigma(\bm{w}^{(t)}_{j}\cdot\bm{x}_{i})-y_{i}\Big{]}a_{r}\bm{x}_{i}^{\intercal}\mathbbm{1}_{{\textstyle\mathstrut}{\bm{w}_{r}^{(t)}}}(\bm{x}_{i})\\\
&=-\mu
a_{r}\sum_{i=1}^{n}\mathbbm{1}_{{\textstyle\mathstrut}{\bm{w}_{r}^{(t)}}}(\bm{x}_{i})\Big{[}\Psi^{(t)}(\bm{x}_{i})\cdot\bm{x}_{i}-y_{i}\Big{]}\bm{x}_{i}^{\intercal}\qquad\qquad\textrm{where}~{}\Psi^{(t)}(\bm{x}_{i})=\sum_{j=1}^{\mathbbm{m}}a_{j}\bm{w}^{(t)}_{j}\mathbbm{1}_{{\textstyle\mathstrut}{\bm{w}_{j}^{(t)}}}(\bm{x}_{i})\end{split}$

Above $\mathbbm{1}_{{\textstyle\mathstrut}{\bm{w}_{r}^{(t)}}}(\bm{x}_{i})$
denotes the indicator function that equals 1 when
$\bm{w}^{(t)}_{r}\cdot\bm{x}_{i}\geq 0$, and $0$ otherwise.

To proceed, we make two assumptions:

  1. 1.

Each point $\bm{x}_{i}$ is drawn from a symmetric distribution $\mathcal{D}$
with density $f_{\mathcal{D}}({\mathbf{X}})$, such that:
$f_{\mathcal{D}}(\bm{x}_{i})=f_{\mathcal{D}}(-\bm{x}_{i})$.

  2. 2.

$W$ and $\bm{a}$ are initialized so that
$\bm{w}^{(0)}_{{\textstyle\mathstrut}2i}=-\bm{w}^{(0)}_{{\textstyle\mathstrut}2i-1}$
and $a_{{\textstyle\mathstrut}2i}=-a_{{\textstyle\mathstrut}2i-1}$ $\forall
i\in[\frac{\mathbbm{m}}{2}]$.

Theorem 9. Retaining the assumptions stated above, at the beginning of the
learning, the temporal dynamics of the model can be shown to obey the
following update rule:

$W^{(t+1)}\approx
W^{(t)}-\mu\frac{1}{2}\Big{[}(\bm{a}\bm{a}^{\intercal})W^{(t)}\Sigma_{\scriptscriptstyle
XX}-\tilde{M}^{(t)}\Big{]}$

Above $\tilde{M}^{(t)}$ denotes the difference between the centroids of the 2
classes, computed in the half-space defined by
$\bm{w}^{(t)}_{r}\cdot\bm{x}\geq 0$.

Proof  It follows from Assumption 2 that at the beginning of training
$\mathbbm{1}_{{\textstyle\mathstrut}{\bm{w}_{2j}^{(0)}}}(\bm{x}_{i})+\mathbbm{1}_{{\textstyle\mathstrut}{\bm{w}_{2j-1}^{(0)}}}(\bm{x}_{i})=1$,
$\forall\bm{x}_{i}$ such that
$\bm{w}_{{\textstyle\mathstrut}2j-1}\bm{x}_{i}\neq\bm{w}_{{\textstyle\mathstrut}2j}\bm{x}_{i}\neq
0$, and $\forall j\in[\frac{\mathbbm{m}}{2}]$. Consequently

$\Psi^{(0)}(\bm{x}_{i})=\sum_{j=1}^{\mathbbm{m}}a_{j}\bm{w}^{(0)}_{j}\mathbbm{1}_{{\textstyle\mathstrut}{\bm{w}_{j}^{(t)}}}(\bm{x}_{i})=\frac{1}{2}\sum_{j=1}^{\mathbbm{m}}a_{j}\bm{w}^{(0)}_{j}=\frac{1}{2}\bm{a}^{\intercal}W^{(0)}$

$\forall\bm{x}_{i}$ such that
$\bm{w}_{{\textstyle\mathstrut}2j-1}\bm{x}_{i}\neq\bm{w}_{{\textstyle\mathstrut}2j}\bm{x}_{i}\neq
0$. Finally

$\bm{w}^{(1)}_{r}-\bm{w}^{(0)}_{r}=-\mu
a_{r}\big{[}\frac{1}{2}\bm{a}^{\intercal}W^{(0)}\sum_{\begin{subarray}{c}i=1\\\
\bm{w}^{(0)}_{r}\bm{x}_{i}\geq
0\end{subarray}}^{n}\bm{x}_{i}\bm{x}_{i}^{\intercal}-\sum_{\begin{subarray}{c}i=1\\\
\bm{w}^{(0)}_{r}\bm{x}_{i}\geq
0\end{subarray}}^{n}y_{i}\bm{x}_{i}^{\intercal}\big{]}$

Next, we note that Assumption 1 implies

$\mathbb{E}[\sum_{\begin{subarray}{c}i=1\\\ {\bm{w}\cdot\bm{x}_{i}\geq
0}\end{subarray}}^{n}\bm{x}_{i}\bm{x}_{i}^{\intercal}]=\frac{1}{2}\mathbb{E}[\sum_{i=1}^{n}\bm{x}_{i}\bm{x}_{i}^{\intercal}]=\frac{1}{2}\mathbb{E}[\Sigma_{\scriptscriptstyle
XX}]$

for any vector $\bm{w}$. Thus, if the sample size $n$ is large enough, at the
beginning of training we expect to see

$\bm{w}^{(t+1)}_{r}-\bm{w}^{(t)}_{r}\approx-\mu\frac{a_{r}}{2}[\bm{a}^{\intercal}W^{(t)}\Sigma_{\scriptscriptstyle
XX}-\tilde{\bm{m}}^{(t)}_{r}]\quad\forall r$

where row vector $\tilde{\bm{m}}^{(t)}_{r}$ denotes the vector difference
between the centroids of classes $C_{1}$ and $C_{2}$, computed in the half-
space defined by $\bm{w}^{(t)}_{r}\cdot\bm{x}\geq 0$. Finally (for small
${t}$)

$W^{(t+1)}-W^{(t)}\approx-\mu\frac{1}{2}\Big{[}(\bm{a}\bm{a}^{\intercal})W^{(t)}\Sigma_{\scriptscriptstyle
XX}-\tilde{M}^{(t)}\Big{]}$

where $\tilde{M}^{(t)}$ denotes the matrix whose $r$-th row is
$a_{r}\tilde{\bm{m}}^{(t)}_{r}$. This equation is reminiscent of the single-
layer linear model dynamics ${\bm{W}}^{(t+1)}={\bm{W}}^{(t)}-\mu G_{r}^{(t)}$,
and we may conclude that when it holds and using the principal coordinate
system, the rate of convergence of the $j$-th column of $W^{(t)}$ is governed
by the singular value $d_{j}$ .

##  C Additional Empirical Results

###  C.1 Weight Initialization

We evaluate empirically the weight initialization scheme from Def. 4 and (14).
When compared to Glorot uniform initialization, the only difference between
the two schemes lies in how the first and last layers are scaled. Thus, to
highlight the difference between the methods, we analyze a fully connected
linear network with a single hidden layer, whose dimension (the number of
hidden neurons) is much larger than the input and output dimensions. We
trained $N$=$10$ such networks on a binary classification problem, once with
the initialization suggested in (14), and again with Glorot uniform
initialization. While both initialization schemes achieve the same final
accuracy upon convergence, our proposed initialization variant converges
faster on both train and test datasets (see Fig. 14).

(a) train

(b) test

Figure 13: Learning curves of a fully connected linear network with one hidden
layer, trained on the dogs and cats dataset, and initialized by either Glorot
uniform initialization (orange), or the initialization proposed in (14)
(blue).

Figure 14: Evaluations on test-sets projected to the first $P$ principal
components, for different values of $P$ (see legend) of $10$ VGG-19 models
trained on CIFAR-10

###  C.2 Divergence of Gradient Scale Matrices

We now discuss some additional results, complementary to Section 3.3. In Fig.
3, we visualize how the gradient scale matrices $B_{l}^{(t)}(\mathbbm{m})$
remain approximately equal to $I$ much longer than $A_{l}^{(t)}(\mathbbm{m})$,
for a specific layer in a 5-layered model. In Fig. 15 we show a similar trend
in all 5 layers of the network. Note that in the input and output layers,
$B_{l}^{(t)}(\mathbbm{m})$ and $A_{l}^{(t)}(\mathbbm{m})$ are $I$ by
definition, and hence their convergence is self-evident.

(a) Diagonal layer 1

(b) Diagonal layer 2

(c) Diagonal layer 4

(d) Diagonal layer 5

(e) Off-diagonal layer 1

(f) Off-diagonal layer 2

(g) Off-diagonal layer 4

(h) Off-diagonal layer 5

Figure 15: Similar to Fig. 3, showing the dynamics of $A^{(t)}_{l}$ and
$B_{l}^{(t)}$ in additional layers, when training $10$ $5$-layered linear
networks on the small mammals dataset. (a-d) The empirical $L_{2}$-distance of
the diagonal elements of $A^{(t)}_{3}$ and $B_{3}^{(t)}$ from their analytical
value of $\alpha_{l}^{(t)}$ and $\beta_{l}^{(t)}$ –
$diag\left(A^{(t)}_{l}-\alpha_{l}^{(t)}I\right)$ and
$diag\left(B_{l}^{(t)}-\beta_{l}^{(t)}I\right)$ – respectively. (e-h) The
empirical $L_{2}$-distance of the off-diagonal elements of $A^{(t)}_{l}$ and
$B_{l}^{(t)}$ from 0. The networks reach their maximal test accuracy in epoch
$t=100$, before the divergence of $A^{(t)}_{l}$.

###  C.3 Empirical Validation of the PC-bias on Original Model

####  C.3.1 Models with $L_{2}$ Loss

In Section 4, we relaxed some of the theoretical assumptions while pursuing
the empirical investigation, in order to match more commonly used models.
Specifically, we changed the initialization to the commonly used Glorot
initialization, replaced the $L_{2}$ loss with the cross-entropy loss, and
employed SGD instead of the deterministic GD. As can be expected from the
theoretical results, without this relaxation the PC-bias becomes even more
pronounced. To support this claim, we repeated the experiments of Section 4.2
and Fig. 5 with the assumptions of the theoretical analysis, showing the
results in Fig. 16.

To examine both the variability between models and their convergence rate, we
inspect $\bm{w}(r)$ at different time points during learning. The rate of
convergence can be measured directly, by observing the changes in the weights
of each element in $\bm{w}(r)$. These weight values are compared with the
optimal values of the $L_{2}$ solution in each row. The variability between
models is measured by calculating the standard deviation (std) of each
$\bm{w}(r)$ across $N$ models which differ in the random initialization.

Specifically, we trained $10$ repetitions of a $5$-layered linear network with
$L_{2}$ loss and _STD initialization_ (see Def.4). We report a more pronounced
phenomenon – convergence in the directions corresponding to larger principal
components converge faster.

Figure 16: Same as Fig. 5, using the $L_{2}$ loss and the initialization
scheme proposed in Def. 4. Epochs plotted: $1,5,10,50,100,500$.

####  C.3.2 Distance of Convergence in Each Principal Direction

In Figs. 5,16, we see that weights in directions corresponding to larger
principal components converge faster, in the sense that the std across
different repetitions drops to zero faster along these directions. However,
the optimal solution along each direction is drastically different. Directions
corresponding to larger principal components tend to have lower values in
their optimal solution. This could serve as a possible explanation to the
convergence phenomenon – it might be possible that in all directions the
distance of the current solution drops at the same speed, but as directions
corresponding to higher principal values start nearer to the optimal solution,
they seem to converge faster.

To test this hypothesis, we repeated the experiment from the previous section
(plotted in Fig. 16), and computed the distance in each direction from the
optimal solution, normalized by the maximal difference achieved in each
direction. We plot these results in Fig. 17. As suggested by our theory, we
report that the normalized distance in each direction also drops faster in
directions corresponding to larger principal components.

Figure 17: Rate of convergence along the principal directions in different
epochs. The value of the $X$-axis corresponds to the index of a principal
eigenvalue, from the most significant to the least significant. The value of
the $Y$-axis is the $L_{2}$ distance between the mean value of the weights and
the optimal solution, normalized by the maximum distance across all epochs.
The results are plotted for epochs: $1,2,10,50,200$, for $10$ $5$-layered
linear networks trained on the cats and dogs dataset with the $L_{2}$ loss.

(a) $\lambda_{j}$ for Fig. 5(a)

(b) $\lambda_{j}$ for Fig. 5(c)

Figure 18: In the same experimental setup as in Fig. 5, we show the evolution
of $\lambda_{j}$ with epochs.

###  C.4 Projection to Higher PC’s

In Section 4.3 we described an evaluation methodology, based on the creation
of a modified _test-set_ by projecting each test example on the span of the