# The Law of Parsimony in Gradient Descent for Learning Deep Linear Networks

Can Yaras111The first two authors contributed to this work equally. Department
of Electrical Engineering & Computer Science, University of Michigan  Peng
Wang  Department of Electrical Engineering & Computer Science, University of
Michigan  Wei Hu  Department of Electrical Engineering & Computer Science,
University of Michigan  Zhihui Zhu  Department of Computer Science &
Engineering, Ohio State University  Laura Balzano  Department of Electrical
Engineering & Computer Science, University of Michigan  Qing Qu  Department of
Electrical Engineering & Computer Science, University of Michigan

###### Abstract

Over the past few years, an extensively studied phenomenon in training deep
networks is the implicit bias of gradient descent towards parsimonious
solutions. In this work, we investigate this phenomenon by narrowing our focus
to deep linear networks. Through our analysis, we reveal a surprising “law of
parsimony” in the learning dynamics when the data possesses low-dimensional
structures. Specifically, we show that the evolution of gradient descent
starting from orthogonal initialization only affects a minimal portion of
singular vector spaces across all weight matrices. In other words, the
learning process happens only within a small invariant subspace of each weight
matrix, despite the fact that all weight parameters are updated throughout
training. This simplicity in learning dynamics could have significant
implications for both efficient training and a better understanding of deep
networks. First, the analysis enables us to considerably improve training
efficiency by taking advantage of the low-dimensional structure in learning
dynamics. We can construct smaller, equivalent deep linear networks without
sacrificing the benefits associated with the wider counterparts. Second, it
allows us to better understand deep representation learning by elucidating the
linear progressive separation and concentration of representations from
shallow to deep layers. We also conduct numerical experiments to support our
theoretical results. The code for our experiments can be found at
https://github.com/cjyaras/lawofparsimony.

##  1 Introduction

In recent years, deep learning has demonstrated remarkable success across a
wide range of applications in engineering and science [1]. Numerous studies
have shown that the effectiveness of deep learning is partially due to the
implicit bias of its learning dynamics, which favors some particular solutions
that generalize exceptionally well without overfitting in the over-
parameterized setting [2, 3, 4, 5]. To gain insight into the implicit bias of
deep networks, a line of recent work has shown that gradient descent (GD)
tends to learn simple functions [6, 7, 8, 9, 10]. For instance, some studies
have shown that gradient descent is biased towards max-margin solutions in
linear networks trained for binary classification via separable data [11, 10].
In addition to the simplicity bias, another line of work showed that deep
networks trained by GD exhibit a bias towards low-rank solutions [5, 12]. The
works [13, 14] demonstrated that adding depth to a matrix factorization
enhances an implicit tendency towards low-rank solutions, leading to more
accurate recovery.

Despite the abundant empirical evidence in practical nonlinear networks, most
theoretical results are developed based on over-parameterized linear models
[8, 11, 5, 14]. Notably, deep linear networks (DLNs), defined by multiple
hidden layers with identity activations between layers, have been widely used
as prototypes of practical deep networks for studying their nonlinear learning
dynamics [15, 16, 17, 13]. Moreover, despite its simplicity, some properties
of DLNs resemble those of their nonlinear counterparts. For example, for both
linear and nonlinear networks, the work [5] empirically showed that the low-
rank bias exists at both initialization and after training and is resilient to
the choice of hyper-parameters and learning methods. The work [16] showed that
a DLN exhibits a striking, hierarchical progressive differentiation of
structures in its internal hidden representations, which is similar to its
nonlinear counterparts. Additionally, [18, 19] demonstrated the practical
applications of DLNs.

Figure 1: Efficient training of deep linear networks. Left: Principal
components of end-to-end GD trajectories for overparameterized DLN and its
equivalent compressed network. Right: Training loss vs. wall-time comparison.

Figure 2: Progressive collapse with linear decay on deep _linear_ and
_nonlinear_ networks. The $x$-axis denotes the layer index and the $y$-axis
denotes the separation measure in (13).

#### Contributions.

In this work, we investigate a particular implicit bias of GD during the
training of DLNs. When the cross-correlation matrix of the training data
exhibits a low-dimensional structure, which manifests as either a low-rank or
wide matrix in our context, we show that the dynamics of GD tend towards
parsimonious solutions. By examining the evolution of GD started from
orthogonal initialization, we unveil a surprising, yet common phenomenon:

_The learning process takes place only within a minimal invariant subspace of
each weight matrix, while the remaining singular subspaces stay unaffected
throughout training._

Here, “unaffected” means that the remaining singular subspaces are not
influenced or altered by the GD updates throughout the entire training
process. Notably, this phenomenon, which we term the “law of parsimony”,
persists despite the fact that GD updates all parameters of the weight
matrices. When we use initialization of small scale, our work sheds new light
on the implicit low-rank bias of the GD trajectory, explicitly explaining why
deeper networks favor low-rank solutions throughout the entire training
process [5]. Moreover, we demonstrate that such a phenomenon could have far-
reaching implications for both understanding and improving the training
efficiency of DLNs, which we highlight below.

  * •

Dramatically more efficient training by constructing “equivalent” smaller
networks. As illustrated in Figure 2, since learning only happens within a
small _invariant_ subspace of the weights, we can construct and train
significantly smaller DLNs that share the same learning dynamics as their
wider counterparts. As such, we can significantly reduce the computational
complexity of training deep networks without sacrificing the benefits of the
associated wider networks, achieving the best of both worlds. We
experimentally demonstrate such benefits on applications like deep matrix
completion [14], showing that we can obtain the same sampling complexity as
the original network while significantly improving the optimization
efficiency. We believe such a finding could also have broad applications in
training and fine-tuning practical deep networks [20].

  * •

New theoretical insights into the progressive data separation in deep
representation learning. We show that the law of parsimony plays an important
role in deciphering representation learning in the context of multi-class
classification problems. Recent works [21, 22, 23, 24, 25] have experimentally
shown that the features across layers of a deep network exhibit a certain law
of data separation that manifests in modern deep architectures during the
terminal phases of training. Specifically, each layer of a trained network
roughly improves a certain measure of data separation by an equal
multiplicative factor [22], which is shown in Figure 2. In this work, we
theoretically investigate this phenomenon based on DLNs with orthogonal data.
Thanks to the “law of parsimony” in the GD trajectory, we can precisely
characterize the _linear_ progressive separation and concentration of
representations from shallow to deep layers, potentially shedding new light on
understanding the phenomenon in deep nonlinear networks.

#### Notations and Organization.

Let $\mathbb{R}^{n}$ be the $n$-dimensional Euclidean space and $\|\cdot\|$ be
the Euclidean norm. Given any $n\in\mathbb{N}$, We use
$\bm{I}_{n}\in\mathbb{R}^{n\times n}$ to denote an identity matrix of size
$n$, and $\bm{1}_{n}\in\mathbb{R}^{n}$ denote an all one vector of length $n$.
Given any $L\in\mathbb{N}$, we use $[L]$ to denote the index set
$\\{1,\cdots,L\\}$. Let $\mathcal{O}^{m\times
n}=\\{\bm{X}\in\mathbb{R}^{m\times n}:\bm{X}^{\top}\bm{X}=\bm{I}_{n}\\}$
denote the set of all $m\times n$ orthogonal matrices and $\mathcal{O}^{n}$
the set of all $n\times n$ orthogonal matrices.

The rest of the paper is organized as follows. In Section 2, we introduce the
basic problem setup. We present our main theoretical result in Section 3. The
applications in deep matrix factorization and progressive collapse are
demonstrated in Section 4.1 and Section 4.2, respectively.

##  2 Problem Formulation

#### Basic Setup of Deep Linear Networks.

Suppose that we have $N$ training samples
$\\{(\bm{x}_{i},\bm{y}_{i})\\}_{i=1}^{N}\subset\mathbb{R}^{d_{x}}\times\mathbb{R}^{d_{y}}$.
Let $\bm{X}=\left[\bm{x}_{1}\ \bm{x}_{2}\ \dots\
\bm{x}_{N}\right]\in\mathbb{R}^{d_{x}\times N}$ and $\bm{Y}=\left[\bm{y}_{1}\
\bm{y}_{2}\ \dots\ \bm{y}_{N}\right]\in\mathbb{R}^{d_{y}\times N}$ and define
$\bm{Y}\bm{X}^{\top}$ to be the cross-correlation matrix. The goal of training
a deep network is to learn a parameterized, hierarchical function
$f:\mathbb{R}^{d_{x}}\mapsto\mathbb{R}^{d_{y}}$ that maps an input
$\bm{x}_{i}\in\mathbb{R}^{d_{x}}$ to its corresponding label
$\bm{y}_{i}\in\mathbb{R}^{d_{y}}$ for all $1\leq i\leq N$. In this work, we
consider an $L$-layer ($L\geq 2$) linear network
$f_{\bm{\Theta}}(\cdot):\mathbb{R}^{d_{x}}\rightarrow\mathbb{R}^{d_{y}}$,
parameterized by $\bm{\Theta}=\\{\bm{W}_{l}\\}_{l=1}^{L}$ with input
$\bm{x}\in\mathbb{R}^{d_{x}}$, i.e.,

$f_{\bm{\Theta}}(\bm{x})\;:=\;\bm{W}_{L}\cdots\bm{W}_{1}\bm{x}\;=\;\bm{W}_{L:1}\bm{x},$
(1)

where $\bm{W}_{1}\in\mathbb{R}^{d_{1}\times d_{x}}$,
$\bm{W}_{l}\in\mathbb{R}^{d_{l}\times d_{l-1}}\;(l=2,\dots,L-1)$, and
$\bm{W}_{L}\in\mathbb{R}^{d_{y}\times d_{L-1}}$ are weight matrices. For
convenience, throughout the paper we adopt the abbreviations
$\bm{W}_{j:i}=\bm{W}_{j}\cdots\bm{W}_{i}$ and
$\bm{W}_{j:i}^{\top}=\bm{W}_{i}^{\top}\cdots\bm{W}_{j}^{\top}$ for $j\geq i$,
where both are identity if $j<i$.

To learn the network parameters $\bm{\Theta}$, we consider minimizing the
$\ell_{2}$ loss on the training data
$\\{(\bm{x}_{i},\bm{y}_{i})\\}_{i=1}^{N}\subseteq\mathbb{R}^{d_{x}}\times\mathbb{R}^{d_{y}}$
as

$\displaystyle\min_{\bm{\Theta}}\
\ell(\bm{\Theta})=\frac{1}{2}\sum_{i=1}^{N}\left\|f_{\bm{\Theta}}(\bm{x}_{i})-\bm{y}_{i}\right\|_{F}^{2}=\frac{1}{2}\left\|{\bm{W}}_{L:1}\bm{X}-\bm{Y}\right\|_{F}^{2}.$
(2)

#### Training DLNs via GD.

As the network is often over-parameterized, a common approach to enforce
implicit regularization in solving the problem is to utilize GD starting from
small initialization [14, 13]. Here, we use $\bm{W}_{l}(t)$ to denote the
weight of the $l$-th layer at the $t$-th iteration.

  * •

Orthogonal initialization. We initialize the weight matrices $\bm{W}_{l}(0)$
for all $l\in[L]$ using $\varepsilon$-scaled orthogonal matrices for some
$\varepsilon>0$, i.e.,

$\displaystyle\bm{W}_{l}^{\top}(0)\bm{W}_{l}(0)=\varepsilon^{2}\bm{I}\quad\mbox{or}\quad\bm{W}_{l}(0)\bm{W}_{l}^{\top}(0)=\varepsilon^{2}\bm{I},\quad\forall
l\in[L],$ (3)

which depends on the size of $\bm{W}_{l}$. It is worth noting that orthogonal
weight initialization is a commonly employed technique in neural network
training due to its ability to speed up the convergence of GD [26, 27, 15, 28,
29].

  * •

Learning dynamics of GD. We update all weights via GD for $t=1,2,\dots$ as

$\displaystyle\bm{W}_{l}(t)=(1-\eta\lambda)\bm{W}_{l}(t-1)-\eta\nabla_{\bm{W}_{l}}\ell(\bm{\Theta}(t-1)),\;\;\forall\;l\in[L],$

where $\lambda\geq 0$ is an optional weight decay parameter and $\eta>0$ is
the learning rate. Substituting the explicit form of the gradient
$\nabla_{\bm{W}_{l}}\ell(\bm{\Theta})$ of (2) into the above equation, we
obtain the analytical form of GD as

$\displaystyle\bm{W}_{l}(t)=(1-\eta\lambda)\bm{W}_{l}(t-1)-\eta\bm{W}_{L:l+1}^{\top}(t-1){\bm{\Gamma}}(t-1)\bm{W}_{l-1:1}^{\top}(t-1),$
(4)

where we denote
${\bm{\Gamma}}(t)=\left(\bm{W}_{L:1}(t)\bm{X}-\bm{Y}\right)\bm{X}^{\top}$ for
simplicity.

##  3 A Law of Parsimony in Gradient Descent of DLNs

Before stating our main result, we make the following assumption for ease of
exposition.

######  Assumption 1.

The weight matrices are square except the last layer, i.e.,
$d_{1}=d_{2}=\cdots=d_{L-1}=d$ for some $d\in\mathbb{N}_{+}$. In particular,
we have $d_{x}=d$. Also, the input data is whitened in the sense that
$\bm{X}\bm{X}^{\top}=\bm{I}_{d_{x}}$.222For any full rank
$\bm{X}\in\mathbb{R}^{d_{x}\times N}$ with $N\geq d_{x}$, whitened data can
always be obtained with a data pre-processing step such as preconditioning.

Note that the above assumptions can be relaxed to give very similar results to
the ones we prove here. For example, Assumption 1 can be potentially relaxed
to $d_{L-1}\geq d_{L-2}\geq\cdots\geq d_{1}$. Moreover, empirical evidence
suggests that our main results should approximately hold for any well-
conditioned $\bm{X}$, not necessarily restricted to whitened $\bm{X}$ – we
leave this study for future work. Based on the assumption above, we show that
all iterates $\bm{W}_{l}(t)$ along the GD trajectory exhibit parsimonious
structures when the cross-correlation matrix possesses low-dimensional
structures.

Figure 3: Evolution of SVD of weight matrices. We visualize the SVD of the
first layer weight matrix of an $L=3$ layer deep linear network with
$d_{x}=d_{y}=30$ and $r=3$ (Case 1) throughout GD without weight decay. Left:
Magnitude of the $i$-th singular value $\sigma_{i}(t)$ at iteration $t$.
Middle: Angle $\angle(\bm{v}_{i}(t),\bm{v}_{i}(0))$ between the $i$-th right
singular vector at iteration $t$ and initialization. Right: Angle
$\angle(\bm{u}_{i}(t),\bm{u}_{i}(0))$ between the $i$-th left singular vector
at iteration $t$ and initialization.

######  Theorem 1.

Suppose that an $L$-layer $f_{\bm{\Theta}}(\cdot)$ and the training data
$(\bm{X},\bm{Y})$ satisfy Assumption 1. We run GD (4) with weight decay
parameter $\lambda$ and learning rate $\eta$ to train $f_{\bm{\Theta}}(\cdot)$
starting from $\varepsilon$-scaled orthogonal initialization (3). Then, the
iterates $\left\\{\bm{W}_{l}(t)\right\\}_{l=1}^{L}$ for all $t\geq 0$ possess
parsimonious structures in the following sense:

  * Case 1.

Suppose the cross-correlation $\bm{Y}\bm{X}^{\top}\in\mathbb{R}^{d_{y}\times
d_{x}}$ is of rank $r\in\mathbb{N}_{+}$ with $d_{y}=d_{x}$, and
$m:=d_{x}-2r>0$. Then there exist orthogonal matrices
$\\{\bm{U}_{l}\\}_{l=1}^{L}\subseteq\mathcal{O}^{d}$ and
$\\{\bm{V}_{l}\\}_{l=1}^{L}\subseteq\mathcal{O}^{d}$ satisfying
$\bm{V}_{l+1}=\bm{U}_{l}$ for all $l\in[L-1]$, such that $\bm{W}_{l}(t)$
admits the following decomposition

$\displaystyle\bm{W}_{l}(t)=\bm{U}_{l}\begin{bmatrix}\widetilde{\bm{W}}_{l}(t)&\bm{0}\\\
\bm{0}&\rho(t)\bm{I}_{m}\end{bmatrix}\bm{V}_{l}^{\top}$ (5)

for all $l\in[L]$ and $t\geq 0$, where
$\widetilde{\bm{W}}_{l}(t)\in\mathbb{R}^{2r\times 2r}$ for all $l\in[L]$ with
$\widetilde{\bm{W}}_{l}(0)=\varepsilon\bm{I}_{2r}$, and

$\displaystyle\rho(t)=\rho(t-1)\left(1-\eta\lambda-\eta\cdot\rho(t-1)^{2(L-1)}\right)$
(6)

for all $t\geq 1$ with $\rho(0)=\varepsilon$.

  * Case 2.

Suppose the cross-correlation $\bm{Y}\bm{X}^{\top}\in\mathbb{R}^{d_{y}\times
d_{x}}$ with $d_{y}=r$ satisfies $m:=d_{x}-2d_{y}>0$. Then, $\bm{W}_{l}(t)$
admits the same decomposition as in (5) for all $l\in[L-1]$ and $t\geq 0$
except

$\displaystyle\rho(t)=\varepsilon\left(1-\eta\lambda\right)^{t},\ \forall
t\geq 0.$ (7)

We defer comprehensive comparison with prior arts and detailed proofs to
Section 5 and Appendix B, respectively. To help the reader gain more insights
from our results, we make several remarks in the following.

#### Dynamics of singular values and vectors of weight matrices.

It is worth noting that the decomposition (5) is closely related to the
singular value decomposition (SVD) of $\bm{W}_{l}(t)$. Specifically, let
$\bm{U}_{l}=[\bm{U}_{l,1}\ \bm{U}_{l,2}]$, $\bm{V}_{l}=[\bm{V}_{l,1}\
\bm{V}_{l,2}]$, where $\bm{U}_{l,1},\bm{V}_{l,1}\in\mathcal{O}^{d\times 2r}$,
$\bm{U}_{l,2},\bm{V}_{l,2}\in\mathcal{O}^{d\times(d-2r)}$. Let
$\widetilde{\bm{W}}_{l}(t)=\widetilde{\bm{U}}_{l}(t)\widetilde{\bm{\Sigma}}_{l}(t)\widetilde{\bm{V}}_{l}^{\top}(t)$
be an SVD of $\widetilde{\bm{W}}_{l}(t)$, where
$\widetilde{\bm{U}}_{l}(t),\widetilde{\bm{V}}_{l}(t)\in\mathcal{O}^{2r}$ and
$\widetilde{\bm{\Sigma}}_{l}(t)\in\mathbb{R}^{2r\times 2r}$ is a diagonal
matrix. Then, we can rewrite (5) into

$\displaystyle\bm{W}_{l}(t)=\begin{bmatrix}\bm{U}_{l,1}\widetilde{\bm{U}}_{l}(t)&\bm{U}_{l,2}\end{bmatrix}\begin{bmatrix}\widetilde{\bm{\Sigma}}_{l}(t)&\bm{0}\\\
\bm{0}&\rho(t)\bm{I}_{m}\end{bmatrix}\begin{bmatrix}\bm{V}_{l,1}\widetilde{\bm{V}}_{l}(t)&\bm{V}_{l,2}\end{bmatrix}^{\top},$
(8)

which is essentially an SVD of $\bm{W}_{l}(t)$ (besides the ordering of
singular values). According to this, we can verify that the (repeated)
singular value $\rho(t)$ undergoes minimal changes across iterations when
$\varepsilon$ is small according to (6) and (7) – this is illustrated in
Figure 3 (left).

#### Low-rank implicit bias.

We emphasize that our result sheds new light on the implicit low-rank bias of
GD. Specifically, it follows from (6) and (7) that
$\lim_{\varepsilon\rightarrow 0}\rho(t)=0$ for all $t\geq 0$. This, together
with (8), implies that the dynamics of GD are inherently biased towards
finding low-rank solutions with a rank of at most $2r$. In contrast to
existing work that demonstrates the tendency of GD to find low nuclear-norm
solutions [12, 14], we directly show that GD tends to find low-rank solutions.

#### Invariance of subspaces.

According to (8), it is evident that the subspace of dimension $d-2r$ formed
by left (resp. right) singular vectors in $\bm{U}_{l,2}$ (resp.
$\bm{V}_{l,2}$) corresponding to the singular values in $\rho(t)$ remains
unchanged during iterations; see Figure 3. This indicates that the learning
process occurs only within an invariant subspace of dimension $2r$. This
result allows us to gain insights into and improve the training efficiency of
DLNs for deep matrix completion (Section 4.1), as well as elucidate the linear
progressive separation of representations of DLNs (Section 4.2).

#### Comparison to prior arts.

In our analysis, we specifically investigate the impact of weight decay on the
implicit bias of GD. Unlike previous work on implicit bias [30, 31, 14, 32],
which did not explicitly consider weight decay, we carefully examine the
effect of this regularization technique. In particular, when weight decay
regularizer $\lambda>0$ is applied, we observe that the singular value
$\rho(t)$ tends to zero asymptotically as $t$ goes to infinity. This finding
indicates that gradient descent with weight decay is biased towards finding
low-rank solutions. Moreover, we believe our result can be generalized to
other optimization methods beyond GD, such as Adam [33], AdaGrad [34], and
RMSprop [35].

##  4 Applications and Experiments

In this section, we focus on showcasing two specific applications of Theorem
1. In Section 4.1, we present the utilization of Case 1 in Theorem 1 to
improve the speed of training in low-rank deep matrix completion. In Section
4.2, we apply the findings from Case 2 in Theorem 1 to gain insights into the
phenomenon of progressive feature collapse in representation learning.

###  4.1 Application I: Accelerating Deep Low-Rank Matrix Completion

First, we demonstrate how the parsimonious structures of GD in DLNs can be
applied to dramatically improve the optimization efficiency of solving deep
matrix completion [14].

#### Problem Setup.

We consider the low-rank matrix completion problem [36, 37, 38] with ground-
truth $\bm{\Phi}\in\mathbb{R}^{d\times d}$ with $r:=\mbox{rank}(\bm{\Phi})\ll
d$. Our goal is to recover $\bm{\Phi}$ from as few number of observations as
possible, where the observed entries are encoded by an index matrix
$\bm{\Omega}\in\\{0,1\\}^{d\times d}$. To solve the problem, we consider the
recent deep matrix factorization approach [14], by optimizing a variant of
Problem (2) with $\bm{Y}=\bm{\Phi}$ and identity input $\bm{X}=\bm{I}_{d}$,
where the objective333The recovery error is defined by flipping the entries in
$\bm{\Omega}$ above, i.e., the error in the unobserved entries.  is defined as

$\ell_{\mathrm{mc}}(\bm{\Theta}):=\frac{1}{2}\|\bm{\Omega}\odot(\bm{W}_{L:1}-\bm{\Phi})\|_{F}^{2}.$
(9)

When the complete observation $\bm{\Omega}=\bm{1}_{d}\bm{1}_{d}^{\top}$ is
available, the above problem simplifies to deep matrix factorization, as
depicted in Problem (2). Moreover, if the network depth is $L=2$, Problem (9)
reduces to a Burer-Monteiro factorization [39]. Despite its nonconvexity,
significant advances have been made in understanding its global optimality and
GD convergence under various settings in the past few years [40, 41, 42, 43,
44, 45, 12, 46, 47, 48, 49].

#### Benefits of Deep Networks.

In practice, the true rank $r$ is often not known exactly. Instead, we may
have a rough estimate of its upper bound $\widehat{r}$, i.e.,
$r\leq\widehat{r}\ll d$. When we consider Problem (9) in the over-
parameterized regime by overestimating the rank as $\widehat{r}\geq r$, as
recent work has shown that GD starting from small initialization has an
implicit bias towards the true rank solution [12], we do not need an explicit
regularization and set the weight decay $\lambda=0$ when we run GD for
optimizing Problem (9). Furthermore, in the over-parameterized regime, more
recent work [14] demonstrated that using a deeper network (i.e., $L\geq 3$) in
solving Problem (9) enjoys several substantial benefits over the shallow
counterpart $L=2$.

  * •

Benefits of depth. As shown in Figure 5 (left), when we increase over-
parameterization, training deeper networks ($L=3$) with GD is less prone to
overfitting. Additionally, the work [14] has shown that deeper networks
improve sample complexity over their shallow counterparts. It has also been
shown that GD for deeper networks has a stronger implicit bias towards low-
rank solutions [5].

  * •

Benefits of width. On the other hand, increasing the width of the network
results in accelerated convergence of GD in terms of iterations. As shown in
Figure 5 (right), increasing the network width for a 3-layer network reduces
the number of GD iterations needed to converge.

Nonetheless, the advantages of deeper and wider networks are accompanied by
computational trade-offs: Increasing both the depth and width of a network
significantly increases the number of parameters that need to be optimized,
thereby causing the per-iteration cost of minimizing Problem (9) via GD to be
much greater than that of shallower and narrower networks. In the following,
we show that this computational challenge associated with training deeper and
wider networks can be mitigated by invoking the law of parsimony in Theorem 1,
where we can construct an approximately equivalent but considerably smaller
network to speed up the optimization process.

Figure 4: Benefits of depth and width in overparameterized matrix completion.
Left: Recovery error using $L=2$ vs. $L=3$. Right: Number of GD iterations to
converge.

Figure 5: Comparison of efficiency of compressed networks vs. narrow network
with different overestimated $\widehat{r}$. Left: Number of iterations to
converge. Right: Time to converge.

#### Network Compression via Law of Parsimony.

To build up intuition, we first describe our approach in a simplified setting
of deep matrix factorization, where we have full observation of $\bm{\Phi}$.
Then we extend the idea to over-parameterized low-rank matrix completion.

_$\blacktriangleright$ The Deep Matrix Factorization Setting:_ With
$\bm{\Omega}=\bm{1}_{d}\bm{1}_{d}^{\top}$, (9) now reduces to the vanilla
problem (2) where we can apply Theorem 1 with $m=d-2\widehat{r}>0$.444We note
that we can replace $r$ in Theorem 1 with $\widehat{r}$ and maintain the same
dynamics provided that $d-2\widehat{r}>0$.  Based on (5) as well as the fact
that $\bm{V}_{l+1}=\bm{U}_{l}$ for all $l\in[L-1]$ from Theorem 1, we can
always write the end-to-end matrix $\bm{W}_{L:1}(t)$ as

$\displaystyle\bm{W}_{L:1}(t)=\begin{bmatrix}\bm{U}_{L,1}&\bm{U}_{L,2}\end{bmatrix}\begin{bmatrix}\widetilde{\bm{W}}_{L:1}(t)&\bm{0}\\\
\bm{0}&\rho^{L}(t)\bm{I}_{m}\end{bmatrix}\begin{bmatrix}\bm{V}_{1,1}^{\top}\\\
\bm{V}_{1,2}^{\top}\end{bmatrix}=\bm{U}_{L,1}\widetilde{\bm{W}}_{L:1}(t)\bm{V}_{1,1}^{\top}+\rho^{L}(t)\bm{U}_{L,2}\bm{V}_{1,2}^{\top}$

for all $t\geq 0$, where
$\widetilde{\bm{W}}_{L:1}(t)\in\mathbb{R}^{2\widehat{r}\times 2\widehat{r}}$
is the end-to-end matrix for all the compressed weights,
$\bm{U}_{L,1},\bm{V}_{1,1}\in\mathcal{O}^{d\times 2\widehat{r}}$, and
$\bm{U}_{L,2},\bm{V}_{1,2}\in\mathcal{O}^{d\times m}$. Then, our claim for
deep matrix factorization is that:

_For optimizing ( 9) with small initialization, running GD on the original
weights $\left\\{\bm{W}_{l}\right\\}_{l=1}^{L}\subseteq\mathbb{R}^{d\times d}$
is almost equivalent to running GD on the compressed weights
$\\{\widetilde{\bm{W}}_{l}\\}_{l=1}^{L}\subseteq\mathbb{R}^{2\widehat{r}\times
2\widehat{r}}$._

As a result, since the compressed weights consist of only $4L\widehat{r}^{2}$
parameters as opposed to the $Ld^{2}$ parameters of the original weights,
their optimization can be significantly more efficient when $\widehat{r}\ll
d$. To numerically verify our claim, we train the compressed network for
target $\bm{\Phi}$ with $d=1000$, $\widehat{r}=r=5$ in Figure 2. It can be
observed that the end-to-end GD trajectory of optimizing the compressed
network closely follows the trajectory of the original network, while
converging to the optimal solution an entire order of magnitude earlier. While
our claim may not possess complete rigor, let us briefly outline two
fundamental components to elucidate why our claim above holds true:

  * •

The effects of small initialization $\varepsilon$ and depth $L$. From (6), we
know that $\rho(t)$ depends on $\varepsilon$, with
$\lim_{\varepsilon\rightarrow 0}\rho(t)=0$. As such, the term
$\rho^{L}(t)\cdot\bm{U}_{L,2}\bm{V}_{1,2}^{\top}\approx\bm{0}$ when we use
small initialization $\varepsilon\approx 0$, and its size decreases with
larger depth $L$. Therefore, we have

$\displaystyle\bm{W}_{L:1}(t)\;=\;\bm{U}_{L,1}\widetilde{\bm{W}}_{L:1}(t)\bm{V}_{1,1}^{\top}+\rho^{L}(t)\bm{U}_{L,2}\bm{V}_{1,2}^{\top}\;\approx\;\bm{U}_{L,1}\widetilde{\bm{W}}_{L:1}(t)\bm{V}_{1,1}^{\top},\quad\forall
t\geq 0.$ (10)

  * •

Invariance of weight subspaces and GD dynamics throughout training. Moreover,
we know from Theorem 1 that all $\left\\{\bm{U}_{l,1}\right\\}_{l=1}^{L}$ and
$\left\\{\bm{V}_{l,1}\right\\}_{l=1}^{L}$ remain unchanged throughout the GD
dynamics, and gradients commute under orthogonal transformations across
different layers. Therefore, running GD on the original weights
$\bm{W}_{l}(t)$ is essentially equivalent to running GD on the compressed
weights $\widetilde{\bm{W}}_{l}(t)$ for all $t\geq 0$ and $l\in[L]$.

Figure 6: Ablation of $\gamma$ in compressed network for deep matrix
completion. Comparison of original overparameterized 3-layer network and its
approximate compressed networks $\bm{W}_{\mathrm{comp}}^{(\gamma)}$
($\gamma>0$) and $\bm{W}_{\mathrm{comp}}^{(0)}$ ($\gamma=0$). Left: Principal
components of end-to-end trajectories of each network. Middle: Recovery error
vs. iteration comparison. Right: Recovery error vs wall-time comparison.

_$\blacktriangleright$ Extension to Deep Matrix Completion:_ We should mention
that our result for deep matrix factorization cannot be directly applied to
deep matrix completion – this is due to the fact that the observation matrix
$\bm{\Omega}\odot\bm{\Phi}$ for a generic $\bm{\Omega}$ does not necessarily
have the low-rank structure of $\bm{\Phi}$, and hence the $\bm{V}_{1,1}$ and
$\bm{U}_{L,1}$ factors in (10) do _not_ remain unchanged throughout the GD
iterations for all $t\geq 0$. To deal with this issue, we propose to update
both $\bm{V}_{1,1}(t)$ and $\bm{U}_{L,1}(t)$ factors via GD with learning rate
$\gamma\eta$ in

$\bm{W}_{\mathrm{comp}}^{(\gamma)}(t):=\bm{U}_{L,1}(t)\widetilde{\bm{W}}_{L:1}(t)\bm{V}_{1,1}^{\top}(t).$
(11)

This is done simultaneously with the GD updates on the subnetwork
$\widetilde{\bm{W}}_{L:1}(t)$, which uses the original learning rate $\eta$.
We call $\bm{W}_{\mathrm{comp}}^{(\gamma)}(t)$ the _compressed network_ ,
where $\gamma\in(0,1)$ denotes the discrepancy in the learning rate. More
specifically, we make the following modifications to our method:

  * •

Initialize $(\bm{V}_{1,1}(0),\bm{U}_{L,1}(0))$. We initialize the factor
$(\bm{V}_{1,1}(0),\bm{U}_{L,1}(0))$ using the ones described in (10), but they
are calculated based upon $\bm{\Omega}\odot\bm{\Phi}$ instead of $\bm{\Phi}$.

  * •

Update $(\bm{V}_{1,1}(t),\bm{U}_{L,1}(t))$ with small learning rates. Although
the subspaces $(\bm{V}_{1,1}(t),\bm{U}_{L,1}(t))$ are changing, they are
changing much more slowly compared to the compressed weights
$\widetilde{\bm{W}}_{l}(t)$. Therefore, we update $\bm{V}_{1,1}(t)$ and
$\bm{U}_{L,1}(t)$ using a discrepant and smaller learning rate $\gamma\eta$
with ($\gamma\ll 1$) for all $t\geq 0$, compared to the learning rate $\eta$
used for the compressed weights $\widetilde{\bm{W}}_{l}(t)$.

As illustrated in Figure 6 (left), monitoring the end-to-end GD trajectory on
the compressed network $\bm{W}_{\mathrm{comp}}^{(\gamma)}(t)$ in (11) reveals
that maintaining $(\bm{U}_{L,1}(t),\bm{V}_{1,1}(t))$ unchanged (i.e.,
$\gamma=0$) results in the GD trajectory of $\bm{W}_{\mathrm{comp}}^{(0)}(t)$
deviating from that of the original network during the later stages of
training, leaving the final test error fairly high. This is due to the
accumulation of approximation error throughout the final training phase. On
the other hand, updating $(\bm{V}_{1,1}(t),\bm{U}_{L,1}(t))$ using a small
learning rate $\gamma\eta$ with $\gamma=0.01$ ensures that the GD trajectory
of $\bm{W}_{\mathrm{comp}}^{(\gamma)}(t)$ closely mirrors that of the original
network, resulting in a significantly reduced recovery error upon convergence
of training.

Regarding computation, it is worth noting that the modified approach requires
optimizing an additional $4\widehat{r}d$ parameters compared to the deep
matrix factorization scenario – this is due to the additional updates on
$(\bm{V}_{1,1}(t),\bm{U}_{L,1}(t))$. Nevertheless, compared to optimizing the
original network, our modified approach remains considerably more efficient,
optimizing only $O(d)$ parameters as opposed to the $O(d^{2})$ parameters
required for the original network. Our experimental result in Figure 6 (right)
also supports this, demonstrating that the compressed network converges
significantly faster in terms of time compared to the original network.

#### Compressed Networks vs. Narrow Networks.

Despite the intriguing observation that the effective rank of the changing
subspace is capped at $2\widehat{r}$, it still prompts the following question:
_Does this imply that optimizing a narrow network of the same width
$2\widehat{r}$ would perform just as efficiently as the compressed network
with a true width of $d\gg\widehat{r}$?_

Our experimental results suggest that the answer is _no_ in general – we
compare the training efficiency of a $2\widehat{r}$-compressed network (within
a wide network of width $d\gg\widehat{r}$) versus a narrow network with width
2$\widehat{r}$ under different over-parameterized estimates $\widehat{r}$. As
depicted in Figure 5 (left), the compressed network requires fewer iterations
to reach convergence, and the number of iterations necessary is almost
unaffected by $\widehat{r}$. Consequently, _training compressed networks is
considerably more time-efficient than training narrow networks of the same
size_ , provided that $\widehat{r}$ is not significantly larger than $r$. The
distinction between the compressed and narrow networks underscores the
benefits of wide networks, as previously demonstrated and discussed in Figure
5 (right), where increasing the network width results in faster convergence.
However, increasing the network width alone also increases the number of
parameters. By employing our network compression methodology, we can achieve
the best of both worlds.

#### Experimental Setups.

Regarding our experimental setups, for Figure 5 we consider matrix completion
problem with $d=100$, $r=10$, and $30$% of entries observed, and for Figure 5
and Figure 6, we consider a matrix completion problem with $d=1000$, $r=5$,
and $20$% of entries observed. We optimize deep networks with depth $L=3$ via
GD starting from a small orthogonal initialization of scale
$\varepsilon=10^{-3}$, until the objective (9) achieves a value less than
$10^{-10}$.

###  4.2 Application II: Understanding Progressive Feature Collapse in DLN

For the multi-class classification problem, we employ our general result in
Theorem 1 to demystify the intriguing phenomenon illustrated in Figure 2. We
refer the reader to Appendix C for the proofs of the results presented in this
section.

#### Problem Setup for Multi-Class Classification.

We consider a $K$-class classification problem with training data samples
$\\{(\bm{x}_{k,i},\bm{y}_{k})\\}_{i\in[n_{k}],k\in[K]}$, where
$\bm{x}_{k,i}\in\mathbb{R}^{d}$ is the $i$-th sample in the $k$-th class,
$\bm{y}_{k}\in\mathbb{R}^{K}$ is an one-hot label vector555A one-hot label,
e.g., $\bm{y}_{k}$, has only the $k$-th entry equal to $1$ with the remaining
entries equal to $0$., and the number of samples $n_{k}$ in each class is
balanced with $n_{1}=\cdots=n_{K}=n$. We denote the total number of samples by
$N=nK$. Based upon the above, we train a $L$-layer DLN to learn weights
$\bm{\Theta}=\\{\bm{W}_{l}\\}_{l=1}^{L}$ via minimizing the $\ell_{2}$ loss in
(2), where $\bm{W}_{l}\in\mathbb{R}^{d_{l}\times d_{l-1}}$, $d_{0}=d$, and
$d_{L}=K$. We write the feature $\bm{z}_{k,i}^{l}$ of an input sample
$\bm{x}_{k,i}$ in the $l$-th layer as

$\displaystyle\bm{z}_{k,i}^{l}:=\bm{W}_{l}\dots\bm{W}_{1}\bm{x}_{k,i}=\bm{W}_{l:1}\bm{x}_{k,i},\
\forall l\in[L],$ (12)

and we denote $\bm{z}_{k,i}^{0}=\bm{x}_{k,i}$. To characterize the network’s
capability to separate data separation across layers on the whole training
dataset, we use a metric introduced in [50] as

$\displaystyle
D_{l}\;:=\;\mathrm{Tr}(\bm{\Sigma}_{W}^{l})/\mathrm{Tr}(\bm{\Sigma}_{B}^{l}),$
(13)

where $\bm{\Sigma}_{W}^{l}$ and $\bm{\Sigma}_{B}^{l}$ characterize the
between-class and with-class variabilities for the $l$-th layer respectively
as

$\displaystyle\bm{\Sigma}_{W}^{l}=\frac{1}{N}\sum_{k=1}^{K}\sum_{i=1}^{n}\left(\bm{z}_{k,i}^{l}-\bar{\bm{z}}_{k}^{l}\right)\left(\bm{z}_{k,i}^{l}-\bar{\bm{z}}_{k}^{l}\right)^{\top},\
\bm{\Sigma}_{B}^{l}=\frac{1}{K}\sum_{k=1}^{K}\left(\bar{\bm{z}}_{k}^{l}-\bar{\bm{z}}^{l}\right)\left(\bar{\bm{z}}_{k}^{l}-\bar{\bm{z}}^{l}\right)^{\top}$
(14)

for all $l\in[L-1]$, where
$\bar{\bm{z}}_{k}^{l}=\frac{1}{n}\sum_{i=1}^{n}\bm{z}_{i,k}^{l}$ denotes the
sample of the $k$-th class for the $l$-th layer’s feature, and
$\bar{\bm{z}}^{l}=\frac{1}{K}\sum_{k=1}^{K}\bar{\bm{z}}_{k}^{l}$ denotes the
corresponding global sample mean. Intuitively,
$\mathrm{Tr}(\bm{\Sigma}_{W}^{l})$ measures how well the features
$\bm{z}_{k,i}^{l}$ in the $l$-th layer collapse to their means in each class,
and $\mathrm{Tr}(\bm{\Sigma}_{B})$ measures the discrimination between
classes.

Therefore, by using $\mathrm{Tr}(\bm{\Sigma}_{B})$ as a normalization factor,
the metric $D_{l}$ measures how the data are separated and concentrated to the
class means in the $l$-th layer. The smaller the value of $D_{l}$ is, the more
collapsed the features are. Additionally, it should be noted that the metric
$D_{l}$ can be viewed as a simplification of the original metric
$\mathrm{Tr}(\bm{\Sigma}_{W}^{l}(\bm{\Sigma}_{B}^{l})^{\dagger})$ extensively
studied in [51, 52, 22], by replacing the pseudoinverse with a trace division.

#### Theoretical Result: Progressive Data Separation with Linear Decay.

Based upon the above setup, we are ready to theoretically justify the
progressive data separation phenomenon in Figure 2, where the metric $D_{l}$
decays at least with a linear rate across layers. Our result is based on DLNs
under some mild conditions. We defer all proofs to Appendix C.

Figure 7: Linear decay of feature separation in trained deep networks. Varying
the depth $L$ and activation type (Linear vs ReLU), we plot the separation
measure $D_{l}$ in (13) at each layer $l$ up to the penultimate layer ($L-1$),
along with the best log-linear fit line.

######  Theorem 2.

For a $K$-class classification problem on a balanced dataset, suppose that the
input dataset $\bm{X}\in\mathbb{R}^{d\times N}$ is square and orthogonal. For
an $L$-layer DLN with parameters $\bm{\Theta}=\\{\bm{W}_{l}\\}_{l=1}^{L}$
described in (1) with $d_{l}=d>2K$ for all $l\in[L-1]$, suppose that
$\bm{\Theta}$ satisfies
(i) Global Optimality: $\bm{W}_{L:1}\bm{X}=\bm{Y}.$
(ii) Balancedness: There exists a positive constant
$\varepsilon\leq\left\\{\frac{n^{1/2L}}{\sqrt{30}L\sqrt[4]{d-K}},\frac{(n/2)^{1/4L}}{\sqrt[4]{d-K}},\frac{1}{\sqrt{2(\sqrt{K}+1)}}\right\\}$
such that $\bm{W}_{l+1}^{\top}\bm{W}_{l+1}=\bm{W}_{l}\bm{W}_{l}^{\top},\forall
l\in[L-2],\
\|\bm{W}_{L}^{\top}\bm{W}_{L}-\bm{W}_{L-1}\bm{W}_{L-1}^{\top}\|_{F}\leq\varepsilon^{2}\sqrt{d-K}.$
(iii) Unchanged Spectrum: There exists a positive constant $\varepsilon>0$ and
an index set $\mathcal{A}\subseteq[d]$ with $|\mathcal{A}|=d-2K$ such that for
all $l\in[L-1]$ that $\sigma_{i}(\bm{W}_{l})=\varepsilon,\ \forall
i\in\mathcal{A}.$
Then, it holds for all $l=0,1,\dots,L-2$ that

$\displaystyle{D_{l+1}}/D_{l}\leq 2(\sqrt{K}+1)\varepsilon^{2}.$ (15)

Note that (iii) depends on Theorem 1. Based on the above theorem, we can
directly verify that DLNs trained by gradient flow (i.e., $\eta\to 0_{+}$ in
(4)) with proper initialization exhibit linear progressive collapse.

######  Corollary 1.

Consider the setting in Theorem 2. Suppose that we employ gradient flow to
train Problem (2), and initialize $\bm{W}_{l}(0)$ for all $l\in[L-1]$ to be
arbitrary $\varepsilon$-scaled orthogonal matrices satisfying (3) and
$\bm{W}_{L}(0)=\varepsilon[\bm{U}_{L}\ \bm{0}]$, where
$\bm{U}_{L}\in\mathcal{O}^{K}$ and
$\varepsilon\leq\left\\{\frac{n^{1/2L}}{\sqrt{30}L\sqrt[4]{d-K}},\frac{(n/2)^{1/4L}}{\sqrt[4]{d-K}},\frac{1}{\sqrt{2(\sqrt{K}+1)}}\right\\}$.
If an optimal solution $\bm{\Theta}^{*}$ is found, then (15) holds.

#### Experimental Results.

We conduct numerical experiments to verify our theoretical results. In the
multi-class classification problem, we consider $K=5$ classes and $n=10$
samples per class, and randomly generate the input data
$\bm{X}\in\mathbb{R}^{d\times N}$ with $d=N=50$ such that
$\bm{X}\bm{X}^{\top}=\bm{I}$. We draw the network weights from random
$\varepsilon$-scale orthogonal matrices with $\varepsilon=0.5$, and vary the
depth $L\in\\{6,8,10\\}$ as well as activation
$\sigma\in\\{\mbox{Linear},\mbox{ReLU}\\}$. The results are shown in Figure 7.
First, we can see that the linear progressive decay phenomenon persists across
networks of different depths and even in nonlinear networks. Second, as
suggested by (13), the decay rate shown in Figure 7 is independent of the
depth $L$, suggesting that deeper networks yield greater feature separation at
the penultimate layer. Additional experimental results can be found in
Appendix A.

##  5 Conclusion

This paper offers a comprehensive analysis of the law of parsimony in gradient
descent for learning DLNs, contributing to the ongoing pursuit of more
efficient and effective deep-learning models. By uncovering the mechanisms
that drive parsimonious solutions, we hope to inspire future research and the
development of advanced techniques that harness the power of simplicity in
deep learning. Ultimately, our goal is to bridge the gap between theory and
practice, enabling practitioners to design and train deep learning models with
improved efficiency and effectiveness. Finally, we conclude by providing a
survey of related works and subsequently discuss the connections and
distinctions between our results and the existing literature.

#### Linear networks.

Due to their relative simplicity, deep linear networks are widely used as an
alternative approach for investigating the optimization, generalization, and
representation characteristics of non-linear networks. Regarding optimization
properties, previous works such as [53, 54] have studied the optimization
landscape of deep linear networks. Some recent works [55, 17, 56] have
established convergence guarantees of gradient descent for training deep
linear networks. Many seminal works are devoted to explaining the
generalization ability of deep networks via different approaches, such as
sharpness [57], neural tangent kernel [58, 59], and implicit bias [8, 60, 61].
In regards to understanding deep representations, recent seminal works [51,
62, 52] studied an intriguing phenomenon termed neural collapse, which is
prevalent across different network architectures, datasets, and training
losses during the terminal phase of training. Notably, many researchers have
provided theoretical explanations for this phenomenon by assuming the
unconstrained feature model [63, 52, 64, 65, 66, 67, 68]. This assumption
simplifies over-parameterized nonlinear networks into two-layer linear
networks, enabling a deeper understanding of the underlying mechanisms.

#### Implicit bias.

In recent years, implicit bias (a.k.a. implicit regularization) has played an
important role in understanding the phenomenon that deep neural networks in
the over-parameterized setting often generalize well even when trained without
any explicit regularization. Numerous studies have been dedicated to
unraveling the mysteries of implicit bias, investigating it from both
theoretical and empirical perspectives. In particular, simplicity bias, low-
rank bias, and spectral bias have been extensively explored in the literature.
Simplicity bias refers to the tendency of gradient descent (GD) for training
deep networks to learn only the simplest features over other useful but more
complex features [6, 31, 69], e.g., it has been shown [11, 70, 8] that GD
favors max-margin solutions in linear models for classification problems. Low-
rank bias refers to the notion that deep networks trained by GD are biased
toward low-rank solutions [12, 14, 60]. For instance, [12] showed that GD for
solving matrix factorization is biased towards minimum nuclear norm solutions.
The line of works studied the robustness with overparameterization with
implicit bias [71, 72, 73]. The works [13, 14] demonstrated that adding depth
to matrix factorization enhances an implicit tendency towards minimum nuclear
norm solutions, leading to more accurate recovery. Recently, [60] provided
substantial empirical observations and concluded that deep networks exhibit an
inductive bias towards solutions with lower effective ranks. As for spectral
bias, it describes the phenomenon that the learning dynamics of deep networks
tends to find low-frequency functions. [74, 75] studied the spectral bias of
deep networks using tools from Fourier analysis.

#### Comparison to existing works.

Here, we would like to highlight the differences and connections between
existing work and our own. First, most of the existing works analyze the
implicit bias by investigating the dynamics of gradient flow [14, 31, 76],
whereas our work directly focuses on the dynamics of GD by utilizing the low-
dimensional structure of the cross-correlation matrix $\bm{Y}\bm{X}^{\top}$.
Our main proof idea is to show that the gradient updates across all iterates
are low rank and share a common nullspace throughout the entire training
process. Since the weights are initialized orthogonally, the right singular
vectors can be chosen arbitrarily to align with the nullspace of gradient
updates so that they remain fixed. Second, it is worth noting that we can
incorporate weight decay regularization into our analysis. Unlike previous
work on implicit bias [30, 31, 14, 32] which did not explicitly consider
weight decay, we carefully examine the effect of this regularization
technique. In particular, when weight decay regularizer $\lambda>0$ is
applied, we observe that the singular value $\rho(t)$ tends to zero
asymptotically as $t$ goes to infinity. This finding indicates that gradient
descent with weight decay is also biased towards finding low-rank solutions.
Third, compared to the dynamics analysis of SVD of the product matrix (i.e.,
$\bm{W}_{L:1}$) in [14], which characterizes its singular values and vectors
through partial differential equations, Theorem 1 directly characterizes the
dynamics of SVD of each weight matrix (i.e., $\bm{W}_{l}$) via revealing the
structures of singular values and vectors. More precisely, our results not
only accurately predict the singular values of each weight matrix during the
training process, as shown in (6) and (7), but also demonstrate that the
subspace corresponding to the singular values in $\rho(t)$ remains unchanged
throughout the iterations.

## Acknowledgment

CY and QQ acknowledge support from U-M START & PODS grants, NSF CAREER
CCF-2143904, NSF CCF-2212066, and NSF CCF-2212326. QQ and PW also acknowledge
support from ONR N00014-22-1-2529, an AWS AI Award, and a gift grant from KLA.
PW and LB acknowledge support from DoE award DE-SC0022186, ARO YIP
W911NF1910027, and NSF CAREER CCF-1845076. ZZ acknowledges support from NSF
grant CCF-2240708. WH acknowledges support from the Google Research Scholar
Program.

## References

  * [1] Yann LeCun, Yoshua Bengio, and Geoffrey Hinton.  Deep learning.  nature, 521(7553):436–444, 2015.
  * [2] Behnam Neyshabur.  Implicit regularization in deep learning.  arXiv preprint arXiv:1709.01953, 2017.
  * [3] Preetum Nakkiran, Gal Kaplun, Yamini Bansal, Tristan Yang, Boaz Barak, and Ilya Sutskever.  Deep double descent: Where bigger models and more data hurt.  Journal of Statistical Mechanics: Theory and Experiment, 2021(12):124003, 2021.
  * [4] Mikhail Belkin, Daniel Hsu, Siyuan Ma, and Soumik Mandal.  Reconciling modern machine-learning practice and the classical bias–variance trade-off.  Proceedings of the National Academy of Sciences, 116(32):15849–15854, 2019.
  * [5] Minyoung Huh, Hossein Mobahi, Richard Zhang, Brian Cheung, Pulkit Agrawal, and Phillip Isola.  The low-rank simplicity bias in deep networks.  Transactions on Machine Learning Research, 2023.
  * [6] Harshay Shah, Kaustav Tamuly, Aditi Raghunathan, Prateek Jain, and Praneeth Netrapalli.  The pitfalls of simplicity bias in neural networks.  Advances in Neural Information Processing Systems, 33:9573–9585, 2020.
  * [7] Guillermo Valle-Perez, Chico Q. Camargo, and Ard A. Louis.  Deep learning generalizes because the parameter-function map is biased towards simple functions.  In International Conference on Learning Representations, 2019.
  * [8] Suriya Gunasekar, Jason D Lee, Daniel Soudry, and Nati Srebro.  Implicit bias of gradient descent on linear convolutional networks.  Advances in neural information processing systems, 31, 2018.
  * [9] Ziwei Ji and Matus Telgarsky.  Gradient descent aligns the layers of deep linear networks.  arXiv preprint arXiv:1810.02032, 2018.
  * [10] Daniel Kunin, Atsushi Yamamura, Chao Ma, and Surya Ganguli.  The asymmetric maximum margin bias of quasi-homogeneous neural networks.  arXiv preprint arXiv:2210.03820, 2022.
  * [11] 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\.
  * [12] Suriya Gunasekar, Blake E Woodworth, Srinadh Bhojanapalli, Behnam Neyshabur, and Nati Srebro.  Implicit regularization in matrix factorization.  Advances in Neural Information Processing Systems, 30, 2017.
  * [13] 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, 2019.
  * [14] Sanjeev Arora, Nadav Cohen, Wei Hu, and Yuping Luo.  Implicit regularization in deep matrix factorization.  Advances in Neural Information Processing Systems, 32, 2019.
  * [15] Andrew M Saxe, James L McClelland, and Surya Ganguli.  Exact solutions to the nonlinear dynamics of learning in deep linear neural networks.  arXiv preprint arXiv:1312.6120, 2013.
  * [16] 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.
  * [17] 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. PMLR, 2018.
  * [18] Sefi Bell-Kligler, Assaf Shocher, and Michal Irani.  Blind super-resolution kernel estimation using an internal-gan.  Advances in Neural Information Processing Systems, 32, 2019.
  * [19] Pierre De Handschutter, Nicolas Gillis, and Xavier Siebert.  Deep matrix factorizations.  arXiv preprint arXiv:2010.00380, 2020.
  * [20] Edward J Hu, yelong shen, Phillip Wallis, Zeyuan Allen-Zhu, Yuanzhi Li, Shean Wang, Lu Wang, and Weizhu Chen.  LoRA: Low-rank adaptation of large language models.  In International Conference on Learning Representations, 2022.
  * [21] Vardan Papyan.  Traces of class/cross-class structure pervade deep learning spectra.  The Journal of Machine Learning Research, 21(1):10197–10260, 2020\.
  * [22] Hangfeng He and Weijie J Su.  A law of data separation in deep learning.  arXiv preprint arXiv:2210.17020, 2022.
  * [23] Xiao Li, Sheng Liu, Jinxin Zhou, Xinyu Lu, Carlos Fernandez-Granda, Zhihui Zhu, and Qing Qu.  Principled and efficient transfer learning of deep models via neural collapse.  arXiv preprint arXiv:2212.12206, 2022.
  * [24] Akshay Rangamani, Marius Lindegaard, Tomer Galanti, and Tomaso Poggio.  Feature learning in deep classifiers through intermediate neural collapse.  Technical report, Center for Brains, Minds and Machines (CBMM), 2023.
  * [25] Shuo Xie, Jiahao Qiu, Ankita Pasad, Li Du, Qing Qu, and Hongyuan Mei.  Hidden state variability of pretrained language models can guide computation reduction for transfer learning.  arXiv preprint arXiv:2210.10041, 2022.
  * [26] Minmin Chen, Jeffrey Pennington, and Samuel Schoenholz.  Dynamical isometry and a mean field theory of rnns: Gating enables signal propagation in recurrent neural networks.  In International Conference on Machine Learning, pages 873–882. PMLR, 2018.
  * [27] Jeffrey Pennington, Samuel Schoenholz, and Surya Ganguli.  The emergence of spectral universality in deep networks.  In International Conference on Artificial Intelligence and Statistics, pages 1924–1932. PMLR, 2018.
  * [28] Lechao Xiao, Yasaman Bahri, Jascha Sohl-Dickstein, Samuel Schoenholz, and Jeffrey Pennington.  Dynamical isometry and a mean field theory of cnns: How to train 10,000-layer vanilla convolutional neural networks.  In International Conference on Machine Learning, pages 5393–5402. PMLR, 2018.
  * [29] Wei Hu, Lechao Xiao, and Jeffrey Pennington.  Provable benefit of orthogonal initialization in optimizing deep linear networks.  arXiv preprint arXiv:2001.05992, 2020.
  * [30] Hancheng Min, Salma Tarmoun, René Vidal, and Enrique Mallada.  Convergence and implicit bias of gradient flow on overparametrized linear networks.  arXiv preprint arXiv:2105.06351, 2022.
  * [31] Daniel Gissin, Shai Shalev-Shwartz, and Amit Daniely.  The implicit bias of depth: How incremental learning drives generalization.  arXiv preprint arXiv:1909.12051, 2019.
  * [32] Gal Vardi and Ohad Shamir.  Implicit regularization in relu networks with the square loss.  In Conference on Learning Theory, pages 4224–4258. PMLR, 2021.
  * [33] Diederik P Kingma and Jimmy Ba.  Adam: A method for stochastic optimization.  arXiv preprint arXiv:1412.6980, 2014.
  * [34] John Duchi, Elad Hazan, and Yoram Singer.  Adaptive subgradient methods for online learning and stochastic optimization.  Journal of machine learning research, 12(7), 2011.
  * [35] Tijmen Tieleman and Geoffrey Hinton.  Rmsprop: Divide the gradient by a running average of its recent magnitude. coursera: Neural networks for machine learning.  COURSERA Neural Networks Mach. Learn, 17, 2012.
  * [36] Emmanuel Candes and Benjamin Recht.  Exact matrix completion via convex optimization.  Communications of the ACM, 55(6):111–119, 2012.
  * [37] Emmanuel J Candès and Terence Tao.  The power of convex relaxation: Near-optimal matrix completion.  IEEE Transactions on Information Theory, 56(5):2053–2080, 2010\.
  * [38] Mark A Davenport and Justin Romberg.  An overview of low-rank matrix recovery from incomplete observations.  IEEE Journal of Selected Topics in Signal Processing, 10(4):608–622, 2016.
  * [39] Samuel Burer and Renato DC Monteiro.  A nonlinear programming algorithm for solving semidefinite programs via low-rank factorization.  Mathematical Programming, 95(2):329–357, 2003.
  * [40] Prateek Jain, Praneeth Netrapalli, and Sujay Sanghavi.  Low-rank matrix completion using alternating minimization.  In Proceedings of the forty-fifth annual ACM symposium on Theory of computing, pages 665–674, 2013.
  * [41] Qinqing Zheng and John Lafferty.  Convergence analysis for rectangular matrix completion using burer-monteiro factorization and gradient descent.  arXiv preprint arXiv:1605.07051, 2016.
  * [42] Ruoyu Sun and Zhi-Quan Luo.  Guaranteed matrix completion via non-convex factorization.  IEEE Transactions on Information Theory, 62(11):6535–6579, 2016\.
  * [43] Rong Ge, Jason D Lee, and Tengyu Ma.  Matrix completion has no spurious local minimum.  Advances in neural information processing systems, 29, 2016.
  * [44] Srinadh Bhojanapalli, Behnam Neyshabur, and Nati Srebro.  Global optimality of local search for low rank matrix recovery.  Advances in Neural Information Processing Systems, 29, 2016.
  * [45] Rong Ge, Chi Jin, and Yi Zheng.  No spurious local minima in nonconvex low rank problems: A unified geometric analysis.  In International Conference on Machine Learning, pages 1233–1242. PMLR, 2017.
  * [46] Qiuwei Li, Zhihui Zhu, and Gongguo Tang.  The non-convex geometry of low-rank matrix optimization.  Information and Inference: A Journal of the IMA, 8(1):51–96, 2019\.
  * [47] Yuejie Chi, Yue M Lu, and Yuxin Chen.  Nonconvex optimization meets low-rank matrix factorization: An overview.  IEEE Transactions on Signal Processing, 67(20):5239–5269, 2019\.
  * [48] Yuanzhi Li, Tengyu Ma, and Hongyang Zhang.  Algorithmic regularization in over-parameterized matrix sensing and neural networks with quadratic activations.  In Conference On Learning Theory, pages 2–47. PMLR, 2018.
  * [49] Mahdi Soltanolkotabi, Dominik Stöger, and Changzhi Xie.  Implicit balancing and regularization: Generalization and convergence guarantees for overparameterized asymmetric matrix sensing.  arXiv preprint arXiv:2303.14244, 2023.
  * [50] Tom Tirer, Haoxiang Huang, and Jonathan Niles-Weed.  Perturbation analysis of neural collapse.  arXiv preprint arXiv:2210.16658, 2022.
  * [51] Vardan Papyan, XY Han, and David L Donoho.  Prevalence of neural collapse during the terminal phase of deep learning training.  Proceedings of the National Academy of Sciences, 117(40):24652–24663, 2020.
  * [52] Zhihui Zhu, Tianyu Ding, Jinxin Zhou, Xiao Li, Chong You, Jeremias Sulam, and Qing Qu.  A geometric analysis of neural collapse with unconstrained features.  Advances in Neural Information Processing Systems, 34, 2021.
  * [53] Kenji Kawaguchi.  Deep learning without poor local minima.  Advances in neural information processing systems, 29, 2016.
  * [54] Haihao Lu and Kenji Kawaguchi.  Depth creates no bad local minima.  arXiv preprint arXiv:1702.08580, 2017.
  * [55] Sanjeev Arora, Nadav Cohen, Noah Golowich, and Wei Hu.  A convergence analysis of gradient descent for deep linear neural networks.  arXiv preprint arXiv:1810.02281, 2018.
  * [56] Armin Eftekhari.  Training linear neural networks: Non-local convergence and complexity results.  In International Conference on Machine Learning, pages 2836–2847. PMLR, 2020.
  * [57] Behnam Neyshabur, Srinadh Bhojanapalli, David McAllester, and Nati Srebro.  Exploring generalization in deep learning.  Advances in neural information processing systems, 30, 2017.
  * [58] Arthur Jacot, Franck Gabriel, and Clément Hongler.  Neural tangent kernel: Convergence and generalization in neural networks.  Advances in neural information processing systems, 31, 2018.
  * [59] Nikhil Vyas, Yamini Bansal, and Preetum Nakkiran.  Limitations of the ntk for understanding generalization in deep learning.  arXiv preprint arXiv:2206.10012, 2022.
  * [60] Minyoung Huh, Hossein Mobahi, Richard Zhang, Brian Cheung, Pulkit Agrawal, and Phillip Isola.  The low-rank simplicity bias in deep networks.  arXiv preprint arXiv:2103.10427, 2021.
  * [61] Guillermo Valle-Perez, Chico Q Camargo, and Ard A Louis.  Deep learning generalizes because the parameter-function map is biased towards simple functions.  arXiv preprint arXiv:1805.08522, 2018.
  * [62] Cong Fang, Hangfeng He, Qi Long, and Weijie J Su.  Exploring deep neural networks via layer-peeled model: Minority collapse in imbalanced training.  Proceedings of the National Academy of Sciences, 118(43), 2021.
  * [63] Can Yaras, Peng Wang, Zhihui Zhu, Laura Balzano, and Qing Qu.  Neural collapse with normalized features: A geometric analysis over the riemannian manifold.  arXiv preprint arXiv:2209.09211, 2022.
  * [64] XY Han, Vardan Papyan, and David L Donoho.  Neural collapse under mse loss: Proximity to and dynamics on the central path.  arXiv preprint arXiv:2106.02073, 2021.
  * [65] Jinxin Zhou, Xiao Li, Tianyu Ding, Chong You, Qing Qu, and Zhihui Zhu.  On the optimization landscape of neural collapse under mse loss: Global optimality with unconstrained features.  arXiv preprint arXiv:2203.01238, 2022.
  * [66] Jinxin Zhou, Chong You, Xiao Li, Kangning Liu, Sheng Liu, Qing Qu, and Zhihui Zhu.  Are all losses created equal: A neural collapse perspective.  In Alice H. Oh, Alekh Agarwal, Danielle Belgrave, and Kyunghyun Cho, editors, Advances in Neural Information Processing Systems, 2022.
  * [67] Peng Wang, Huikang Liu, Can Yaras, Laura Balzano, and Qing Qu.  Linear convergence analysis of neural collapse with unconstrained features.  In OPT 2022: Optimization for Machine Learning (NeurIPS 2022 Workshop), 2022.
  * [68] Vignesh Kothapalli.  Neural collapse: A review on modelling principles and generalization.  Transactions on Machine Learning Research, 2023.
  * [69] Depen Morwani, Jatin Batra, Prateek Jain, and Praneeth Netrapalli.  Simplicity bias in 1-hidden layer neural networks.  arXiv preprint arXiv:2302.00457, 2023.
  * [70] Mor Shpigel Nacson, Jason Lee, Suriya Gunasekar, Pedro Henrique Pamplona Savarese, Nathan Srebro, and Daniel Soudry.  Convergence of gradient descent on separable data.  In The 22nd International Conference on Artificial Intelligence and Statistics, pages 3420–3428. PMLR, 2019.
  * [71] Wei Hu, Zhiyuan Li, and Dingli Yu.  Simple and effective regularization methods for training on noisily labeled data with generalization guarantee.  In International Conference on Learning Representations, 2020.
  * [72] Chong You, Zhihui Zhu, Qing Qu, and Yi Ma.  Robust recovery via implicit bias of discrepant learning rates for double over-parameterization.  Advances in Neural Information Processing Systems, 33:17733–17744, 2020.
  * [73] Sheng Liu, Zhihui Zhu, Qing Qu, and Chong You.  Robust training under label noise by over-parameterization.  In International Conference on Machine Learning, pages 14153–14172. PMLR, 2022.
  * [74] Yuan Cao, Zhiying Fang, Yue Wu, Ding-Xuan Zhou, and Quanquan Gu.  Towards understanding the spectral bias of deep learning.  arXiv preprint arXiv:1912.01198, 2019.
  * [75] 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.
  * [76] Hancheng Min, Salma Tarmoun, René Vidal, and Enrique Mallada.  On the explicit role of initialization on the convergence and implicit bias of overparametrized linear networks.  In International Conference on Machine Learning, pages 7760–7768. PMLR, 2021.
  * [77] Gilbert W Stewart and Ji-guang Sun.  Matrix perturbation theory.  Academic press, 1990.

Appendix

In the appendix, we present additional experiments in Appendix A, and provide
complete proofs for the technical results of Sections 3 and 4 in Appendix B
and Appendix C respectively. Before we proceed, we introduce some further
notation. For simplicity, let $\ell(t)=\ell(\bm{\Theta}(t))$. We use
$\sigma_{\max}(\bm{A})$ (or $\|\bm{A}\|$), $\sigma_{i}(\bm{A})$, and
$\sigma_{\min}(\bm{A})$ to denote the largest, the $i$-th largest, and the
smallest singular values of a matrix $\bm{A}$, respectively. Given weight
matrices $\bm{W}_{1},\dots,\bm{W}_{L}$, we denote
$\bm{W}_{j:i}:=\bm{W}_{j}\bm{W}_{j-1}\dots\bm{W}_{i}$ if $1\leq i\leq j\leq
L$, and $\bm{W}_{0:1}:=\bm{I}$ and $\bm{W}_{L+1:L}=\bm{I}$. We denote the
Kronecker product by $\otimes$.

##  Appendix A Additional Experiments

In this section, we present additional experimental results to supplement
those presented in the main paper. We note that all experiments in this work
were carried out on a single NVIDIA Tesla V100 GPU.

Figure 8: Evolution of SVD of weight matrices for all layers. We visualize the
SVD of the weight matrix for each layer of an $L=3$ layer deep linear network
with $d_{x}=d_{y}=30$ and $r=3$ (Case 1) throughout GD without weight decay.
Left: Magnitude of the $i$-th singular value $\sigma_{i}(t)$ at iteration $t$.
Middle: Angle $\angle(\bm{v}_{i}(t),\bm{v}_{i}(0))$ between the $i$-th right
singular vector at iteration $t$ and initialization. Right: Angle
$\angle(\bm{u}_{i}(t),\bm{u}_{i}(0))$ between the $i$-th left singular vector
at iteration $t$ and initialization.  Figure 9: Evolution of SVD of weight
matrices with momentum. We visualize the SVD of the first layer weight matrix
of an $L=3$ layer deep linear network with $d_{x}=d_{y}=30$ and $r=3$ (Case 1)
throughout GD with momentum (16). Left: Magnitude of the $i$-th singular value
$\sigma_{i}(t)$ at iteration $t$. Middle: Angle
$\angle(\bm{v}_{i}(t),\bm{v}_{i}(0))$ between the $i$-th right singular vector
at iteration $t$ and initialization. Right: Angle
$\angle(\bm{u}_{i}(t),\bm{u}_{i}(0))$ between the $i$-th left singular vector
at iteration $t$ and initialization.

First, we extend the visualization in Figure 3 to all layers, rather than
simply the first layer – the results are shown in Figure 8. We see that the
same parsimonious structures persist across all layers throughout GD, as
implied by Theorem 1.

Next, we demonstrate that the law of parsimony generalizes to other variants
of GD beyond weight decay, such as GD with momentum, i.e., we update all
weights for $t=2,3,\dots$ as

$\displaystyle\bm{W}_{l}(t)=\bm{W}_{l}(t-1)-\eta\nabla_{\bm{W}_{l}}\ell(\bm{\Theta}(t-1))+\mu(\bm{W}_{l}(t-1)-\bm{W}_{l}(t-2))$
(16)

for all $l\in[L]$, where $0\leq\mu<1$ is the momentum parameter and $\eta>0$
is the learning rate. We see in Figure 9 that introducing the momentum term
still results in a low-dimensional trajectory along minimal singular
subspaces. Using a similar approach to the analysis presented in this work, it
should be fairly straightforward to extend our results to GD with momentum,
but this is left as future work.

Figure 10: Evolution of SVD of weight matrices without whitened data. We
visualize the SVD of the first layer weight matrix of an $L=3$ layer deep
linear network with $d_{x}=30$, $d_{y}=3$, and $N=300$ throughout GD without
weight decay. Left: Magnitude of the $i$-th singular value $\sigma_{i}(t)$ at
iteration $t$. Middle: Angle $\angle(\bm{v}_{i}(t),\bm{v}_{i}(0))$ between the
$i$-th right singular vector at iteration $t$ and initialization. Right: Angle
$\angle(\bm{u}_{i}(t),\bm{u}_{i}(0))$ between the $i$-th left singular vector
at iteration $t$ and initialization.

Besides the low-dimensional structure of $\bm{Y}\bm{X}^{\top}$, the other data
assumption made in this work is that the input data is whitened, i.e.,
$\bm{X}\bm{X}^{\top}=\bm{I}$. Here, we demonstrate that the law of parsimony
continues to hold approximately when $\bm{X}$ is not whitened or preprocessed
in any way – in particular, we draw the entries of
$\bm{X}\in\mathbb{R}^{d_{x}\times N}$ i.i.d. from the standard normal
distribution. The results are shown in Figure 10. We can see that most of the
singular value remain unchanged from initialization as with before, and the
corresponding singular subspaces evolve very little throughout optimization.

To conclude this section, we present several extensions to the experimental
results in Section 4.2. First, to investigate the effect of initialization
scale $\varepsilon$ on the rate of linear progressive collapse, we fix the
depth $L=4$ and vary $\varepsilon\in\\{0.5,0.25,0.125\\}$, keeping all other
aspects of the setup the same as in Section 4.2. The results are shown in
Figure 11. As implied by (15), we see that decreasing the initialization scale
$\varepsilon$ leads to a steeper decay in the measure $D_{l}$ (13) in both
linear and nonlinear networks, demonstrating that there is indeed a necessary
dependence on $\varepsilon$ in the upper bound (15). Therefore, the
initialization scale should be chosen carefully when training deep networks in
practice, particularly due to the real implications of progressive collapse in
transfer learning [23].

Figure 11: Linear decay of feature separation in trained deep networks with
varying initialization scale $\varepsilon$. Varying the initialization scale
$\varepsilon$ for networks with both linear and ReLU activations for an $L=4$
deep network, we plot the separation measure $D_{l}$ in (13) at each layer $l$
up to the penultimate layer ($L-1$), along with the best log-linear fit line.

Throughout this work, we have assumed that the network weights are initialized
orthogonally for the sake of analysis – we now empirically verify that the
linear progressive collapse phenomenon is independent of the type of
initialization, i.e., the distribution from which the weight matrices are
drawn. With the same experimental set-up as in Section 4.2, we fix the depth
$L=8$ and vary the initialization type among orthogonal, normal, and uniform
initializations. The results are shown in Figure 12. To keep the
initialization scale comparable between different initialization schemes, for
the normally distributed weights, we draw the entries independently from
$\mathcal{N}\left(0,\frac{\varepsilon^{2}}{d}\right)$, and for the uniformly
distributed weights, we draw the entries independently from
$\mathcal{U}\left(-\frac{3\varepsilon}{d},\frac{3\varepsilon}{d}\right)$. This
guarantees that drawing $\bm{W}\in\mathbb{R}^{d\times d}$ from any
distribution gives $\mathop{\rm
E\,\\!}\nolimits[\|\bm{W}\|_{F}^{2}]=d\cdot\varepsilon^{2}$. As a result, we
see that progressive collapse occurs in both linear and nonlinear networks,
and in fact decays at the same rate, regardless of the kind of initialization
used.

Figure 12: Linear decay of feature separation in trained deep networks with
different initialization types. Varying the initialization type (orthogonal vs
normal vs uniform) and activation type (Linear vs ReLU) for an $L=8$ deep
network, we plot the separation measure $D_{l}$ in (13) at each layer $l$ up
to the penultimate layer ($L-1$), along with the best log-linear fit line.

##  Appendix B Proofs in Section 3

In this section, we work towards proving Theorem 1 by establishing Lemmas 1
and 2, which imply the conclusions of Case 1 and Case 2 respectively. Before
proceeding, we note that by the assumption that the data
$\bm{X}\in\mathbb{R}^{d_{x}\times n}$ is whitened, i.e.,
$\bm{X}\bm{X}^{\top}=\bm{I}_{d_{x}}$, we have that ${\bm{\Gamma}}(t)$ in (4)
takes the form

${\bm{\Gamma}}(t)=\bm{W}_{L:1}(t)-\bm{\Phi},\ \text{where}\
\bm{\Phi}=\bm{Y}\bm{X}^{\top}\in\mathbb{R}^{d_{y}\times d_{x}}$ (17)

In addition, we note that all statements quantified by $i$ in this section
implicity hold for all $i\in[m]$ (as defined in Theorem 1) for the sake of
notational brevity.

######  Lemma 1.

Under the setting of Theorem 1 (Case 1), there exist orthonormal sets
$\\{\bm{u}_{i}^{(l)}\\}_{i=1}^{m}\subset\mathbb{R}^{d}$ and
$\\{\bm{v}_{i}^{(l)}\\}_{i=1}^{m}\subset\mathbb{R}^{d}$ for $l\in[L]$
satisfying $\bm{v}_{i}^{(l+1)}=\bm{u}_{i}^{(l)}$ for all $l\in[L-1]$ such that
the following hold for all $t\geq 0$:

$\displaystyle\mathcal{A}(t)$
$\displaystyle:\bm{W}_{l}(t)\bm{v}_{i}^{(l)}=\rho(t)\bm{u}_{i}^{(l)}\quad$
$\displaystyle\forall l\in[L],$ $\displaystyle\mathcal{B}(t)$
$\displaystyle:\bm{W}_{l}^{\top}(t)\bm{u}_{i}^{(l)}=\rho(t)\bm{v}_{i}^{(l)}\quad$
$\displaystyle\forall l\in[L],$ $\displaystyle\mathcal{C}(t)$
$\displaystyle:\bm{\Phi}^{\top}\bm{W}_{L:l+1}(t)\bm{u}_{i}^{(l)}=\bm{0}\quad$
$\displaystyle\forall l\in[L],$ $\displaystyle\mathcal{D}(t)$
$\displaystyle:\bm{\Phi}\bm{W}^{\top}_{l-1:1}(t)\bm{v}_{i}^{(l)}=\bm{0}\quad$
$\displaystyle\forall l\in[L],$

where
$\rho(t)=\rho(t-1)\left(1-\eta\lambda-\eta\cdot\rho(t-1)^{2(L-1)}\right)$ for
all $t\geq 1$ with $\rho(0)=\varepsilon$.

###### Proof.

Define $\bm{\Psi}:=\bm{W}^{\top}_{L:2}(0)\bm{\Phi}$. Since the rank of
$\bm{\Phi}$ is at most $r$, we have that the rank of
$\bm{\Psi}\in\mathbb{R}^{d\times d}$ is at most $r$, which implies that
$\mathrm{dim}\,\mathcal{N}\left(\bm{\Psi}\right)=\mathrm{dim}\,\mathcal{N}\left(\bm{\Psi}^{\top}\right)\geq
d-r$. We define the subspace

$\displaystyle\mathcal{S}:=\mathcal{N}\left(\bm{\Psi}\right)\cap\mathcal{N}\left(\bm{\Psi}^{\top}\bm{W}_{1}(0)\right)\subset\mathbb{R}^{d}.$

Since $\bm{W}_{1}(0)\in\mathbb{R}^{d\times d}$ is nonsingular, we have

$\displaystyle\mathrm{dim}(\mathcal{S})\geq 2(d-r)-d=m.$

Let $\\{\bm{v}_{i}^{(1)}\\}_{i=1}^{m}$ denote an orthonormal basis for
$\mathcal{S}$ and set
$\bm{u}_{i}^{(1)}:=\bm{W}_{1}(0)\bm{v}_{i}^{(1)}/\varepsilon$, where
$\varepsilon>0$ is the constant in (3) – since $\bm{W}_{1}(0)$ is orthogonal,
$\\{\bm{u}_{i}^{(1)}\\}_{i=1}^{m}$ is also an orthonormal set. Then we
trivially have $\bm{W}_{1}(0)\bm{v}_{i}^{(1)}=\varepsilon\bm{u}_{i}^{(1)}$,
which, together with (3), implies
$\bm{W}_{1}^{\top}(0)\bm{u}_{i}^{(1)}=\varepsilon\bm{v}_{i}^{(1)}$. It follows
from $\bm{v}_{i}^{(1)}\in\mathcal{S}$ that $\bm{\Psi}\bm{v}_{i}^{(1)}=\bm{0}$
and $\bm{\Psi}^{\top}\bm{W}_{1}(0)\bm{v}_{i}^{(1)}=\bm{0}$, which is
equivalent to $\bm{W}^{\top}_{L:2}(0)\bm{\Phi}\bm{v}_{i}^{(1)}=\bm{0}$ and
$\bm{\Phi}^{\top}\bm{W}_{L:2}(0)\bm{W}_{1}(0)\bm{v}_{i}^{(1)}=\varepsilon\bm{\Phi}^{\top}\bm{W}_{L:2}(0)\bm{u}_{i}^{(1)}=\bm{0}$
respectively. Since $\bm{W}_{L:2}^{\top}(0)$ is full column rank, we further
have that $\bm{\Phi}\bm{v}_{i}^{(1)}=\bm{0}$.

Now let $\mathcal{E}(l)$ be the event that we have orthonormal sets
$\\{\bm{u}_{i}^{(l)}\\}_{i=1}^{m}$ and $\\{\bm{v}_{i}^{(l)}\\}_{i=1}^{m}$
satisfying $\bm{W}_{l}(0)\bm{v}_{i}^{(l)}=\varepsilon\bm{u}_{i}^{(l)}$,
$\bm{W}_{l}^{\top}(0)\bm{u}_{i}^{(l)}=\varepsilon\bm{v}_{i}^{(l)}$,
$\bm{\Phi}^{\top}\bm{W}_{L:l+1}(0)\bm{u}_{i}^{(l)}=\bm{0}$, and
$\bm{\Phi}\bm{W}_{l-1:1}^{\top}(0)\bm{v}_{i}^{(l)}=\bm{0}$. From the above
arguments, we have that $\mathcal{E}(1)$ holds – now suppose $\mathcal{E}(k)$
holds for some $1\leq k<L$. Set $\bm{v}_{i}^{(k+1)}:=\bm{u}_{i}^{(k)}$ and
$\bm{u}_{i}^{(k+1)}:=\bm{W}_{k+1}(0)\bm{v}_{i}^{(k+1)}/\varepsilon$. This,
together with (3), implies that
$\bm{W}_{k+1}(0)\bm{v}_{i}^{(k+1)}=\varepsilon\bm{u}_{i}^{(k+1)}$ and
$\bm{W}_{k+1}^{\top}(0)\bm{u}_{i}^{(k+1)}=\varepsilon\bm{v}_{i}^{(k+1)}$.
Moreover, we have

$\displaystyle\bm{\Phi}^{\top}\bm{W}_{L:(k+1)+1}(0)\bm{u}_{i}^{(k+1)}$
$\displaystyle=\bm{\Phi}^{\top}\bm{W}_{L:k+1}(0)\bm{W}_{k+1}^{\top}(0)\bm{u}_{i}^{(k+1)}/\varepsilon^{2}$
$\displaystyle=\bm{\Phi}^{\top}\bm{W}_{L:k+1}(0)\bm{v}_{i}^{(k+1)}/\varepsilon$
$\displaystyle=\bm{\Phi}^{\top}\bm{W}_{L:k+1}(0)\bm{u}_{i}^{(k)}/\varepsilon=\bm{0},$

where the first two equalities follow from (3) and
$\bm{u}_{i}^{(k+1)}=\bm{W}_{k+1}(0)\bm{v}_{i}^{(k+1)}/\varepsilon$, and the
last equality is due to $\bm{v}_{i}^{(k+1)}=\bm{u}_{i}^{(k)}$. Similarly, we
have

$\displaystyle\bm{\Phi}\bm{W}_{(k+1)-1:1}^{\top}(0)\bm{v}_{i}^{(k+1)}$
$\displaystyle=\bm{\Phi}\bm{W}_{k-1:1}^{\top}(0)\bm{W}_{k}^{\top}(0)\bm{v}_{i}^{(k+1)}$
$\displaystyle=\bm{\Phi}\bm{W}_{k-1:1}^{\top}(0)\bm{W}_{k}^{\top}(0)\bm{u}_{i}^{(k)}$
$\displaystyle=\varepsilon\bm{\Phi}\bm{W}_{k-1:1}^{\top}(0)\bm{v}_{i}^{(k)}=\bm{0},$

where the second equality follows from $\bm{v}_{i}^{(k+1)}=\bm{u}_{i}^{(k)}$
and the third equality is due to
$\bm{W}_{k}^{\top}(0)\bm{u}_{i}^{(k)}=\varepsilon\bm{v}_{i}^{(k)}$. Therefore
$\mathcal{E}(k+1)$ holds, so we have $\mathcal{E}(l)$ for all $l\in[L]$. As a
result, we have shown the base cases $\mathcal{A}(0)$, $\mathcal{B}(0)$,
$\mathcal{C}(0)$, and $\mathcal{D}(0)$.

Now we proceed by induction on $t\geq 0$. Suppose that $\mathcal{A}(t)$,
$\mathcal{B}(t)$, $\mathcal{C}(t)$, and $\mathcal{D}(t)$ hold for some $t\geq
0$. First, we show $\mathcal{A}(t+1)$ and $\mathcal{B}(t+1)$. We have

$\displaystyle\bm{W}_{l}(t+1)\bm{v}_{i}^{(l)}$
$\displaystyle=\left[(1-\eta\lambda)\bm{W}_{l}(t)-\eta\bm{W}_{L:l+1}^{\top}(t)\bm{\Gamma}(t)\bm{W}_{l-1:1}^{\top}(t)\right]\bm{v}_{i}^{(l)}$
$\displaystyle=\left[(1-\eta\lambda)\bm{W}_{l}(t)-\eta\bm{W}_{L:l+1}^{\top}(t)\left(\bm{W}_{L:1}(t)-\bm{\Phi}\right)\bm{W}_{l-1:1}^{\top}(t)\right]\bm{v}_{i}^{(l)}$
$\displaystyle=(1-\eta\lambda)\bm{W}_{l}(t)\bm{v}_{i}^{(l)}-\eta\bm{W}_{L:l+1}^{\top}(t)\bm{W}_{L:1}(t)\bm{W}_{l-1:1}^{\top}(t)\bm{v}_{i}^{(l)}$
$\displaystyle=(1-\eta\lambda)\bm{W}_{l}(t)\bm{v}_{i}^{(l)}-\eta\cdot\rho(t)^{2(L-1)}\bm{W}_{l}(t)\bm{v}_{i}^{(l)}$
$\displaystyle=\rho(t)\left(1-\eta\lambda-\eta\cdot\rho(t)^{2(L-1)}\right)\bm{u}_{i}^{(l)}=\rho(t+1)\bm{u}_{i}^{(l)}$

for all $l\in[L]$, where the first equality follows from (4), the second
equality follows from (17), the third equality follows from $\mathcal{D}(t)$,
and the fourth equality follows from $\mathcal{A}(t)$ and $\mathcal{B}(t)$
applied repeatedly along with $\bm{v}_{i}^{(l+1)}=\bm{u}_{i}^{(l)}$ for all
$l\in[L-1]$, proving $\mathcal{A}(t+1)$. Similarly, we have

$\displaystyle\bm{W}_{l}^{\top}(t+1)\bm{u}_{i}^{(l)}$
$\displaystyle=\left[(1-\eta\lambda)\bm{W}_{l}^{\top}(t)-\eta\bm{W}_{l-1:1}(t)\bm{\Gamma}^{\top}(t)\bm{W}_{L:l+1}(t)\right]\bm{u}_{i}^{(l)}$
$\displaystyle=\left[(1-\eta\lambda)\bm{W}_{l}^{\top}(t)-\eta\bm{W}_{l-1:1}(t)\bm{\left
missing}(\bm{W}_{L:1}^{\top}(t)-\bm{\Phi}^{\top}\right)\bm{W}_{L:l+1}(t)\right]\bm{u}_{i}^{(l)}$
$\displaystyle=(1-\eta\lambda)\bm{W}_{l}^{\top}(t)\bm{u}_{i}^{(l)}-\eta\bm{W}_{l-1:1}(t)\bm{W}_{L:1}^{\top}(t)\bm{W}_{L:l+1}(t)\bm{u}_{i}^{(l)}$
$\displaystyle=(1-\eta\lambda)\bm{W}_{l}^{\top}(t)\bm{u}_{i}^{(l)}-\eta\cdot\rho(t)^{2(L-1)}\bm{W}_{l}^{\top}(t)\bm{u}_{i}^{(l)}$
$\displaystyle=\rho(t)\left(1-\eta\lambda-\eta\cdot\rho(t)^{2(L-1)}\right)\bm{v}_{i}^{(l)}=\rho(t+1)\bm{v}_{i}^{(l)}$

for all $l\in[L]$, where the third equality follows from $\mathcal{C}(t)$, and
the fourth equality follows from $\mathcal{A}(t)$ and $\mathcal{B}(t)$ applied
repeatedly along with $\bm{v}_{i}^{(l+1)}=\bm{u}_{i}^{(l)}$ for all
$l\in[L-1]$, proving $\mathcal{B}(t+1)$. Now, we show $\mathcal{C}(t+1)$. For
any $k\in[L-1]$, it follows from $\bm{v}_{i}^{(k+1)}=\bm{u}_{i}^{(k)}$ and
$\mathcal{A}(t+1)$ that

$\bm{W}_{k+1}(t+1)\bm{u}_{i}^{(k)}=\bm{W}_{k+1}(t+1)\bm{v}_{i}^{(k+1)}=\rho(t+1)\bm{u}_{i}^{(k+1)}.$

Repeatedly applying the above equality for $k=l,l+1,\dots,L-1$, we obtain

$\bm{\Phi}^{\top}\bm{W}_{L:l+1}(t)\bm{u}_{i}^{(l)}=\rho(t+1)^{L-l}\bm{\Phi}^{\top}\bm{u}_{i}^{(L)}=\bm{0}$

which follows from $\mathcal{C}(t)$, proving $\mathcal{C}(t+1)$. Finally, we
show $\mathcal{D}(t+1)$. For any $k\in\\{2,\dots,L\\}$, it follows from
$\bm{v}_{i}^{(k)}=\bm{u}_{i}^{(k-1)}$ and $\mathcal{B}(t+1)$ that

$\bm{W}_{k-1}^{\top}(t+1)\bm{v}_{i}^{(k)}=\bm{W}_{k-1}^{\top}(t+1)\bm{u}_{i}^{(k-1)}=\rho(t+1)\bm{v}_{i}^{(k-1)}.$

Repeatedly applying the above equality for $k=l,l-1,\dots,2$, we obtain

$\bm{\Phi}\bm{W}^{\top}_{l-1:1}(t)\bm{v}_{i}^{(l)}=\bm{\Phi}\bm{v}_{i}^{(1)}=\bm{0}$

which follows from $\mathcal{D}(t)$. Thus we have proven $\mathcal{D}(t+1)$,
concluding the proof. ∎

######  Lemma 2.

Under the setting of Theorem 1 (Case 2), there exist orthonormal sets
$\\{\bm{u}_{i}^{(l)}\\}_{i=1}^{m}\subset\mathbb{R}^{d}$ for $l\in[L-1]$ and
$\\{\bm{v}_{i}^{(l)}\\}_{i=1}^{m}\subset\mathbb{R}^{d}$ for $l\in[L]$
satisfying $\bm{v}_{i}^{(l+1)}=\bm{u}_{i}^{(l)}$ for all $l\in[L-1]$ such that
the following hold for all $t\geq 0$:

$\displaystyle\mathcal{A}(t)$
$\displaystyle:\bm{W}_{l}(t)\bm{v}_{i}^{(l)}=\rho(t)\bm{u}_{i}^{(l)}\quad$
$\displaystyle\forall l\in[L-1],$ $\displaystyle\mathcal{B}(t)$
$\displaystyle:\bm{W}_{l}^{\top}(t)\bm{u}_{i}^{(l)}=\rho(t)\bm{v}_{i}^{(l)}\quad$
$\displaystyle\forall l\in[L-1],$ $\displaystyle\mathcal{C}(t)$
$\displaystyle:\bm{W}_{L:l+1}(t)\bm{u}_{i}^{(l)}=\bm{0}\quad$
$\displaystyle\forall l\in[L-1],$ $\displaystyle\mathcal{D}(t)$
$\displaystyle:\bm{\Gamma}(t)\bm{W}_{l-1:1}^{\top}(t)\bm{v}_{i}^{(l)}=\bm{0}\quad$
$\displaystyle\forall l\in[L],$

where $\rho(t)=\rho(t-1)(1-\eta\lambda)$ for all $t\geq 1$ with
$\rho(0)=\varepsilon$.

###### Proof.

Since $\bm{\Gamma}(0)\in\mathbb{R}^{d_{y}\times d_{x}}$, the rank of

$\displaystyle\nabla_{\bm{W}_{1}}\ell(0)=\bm{W}_{L:2}^{\top}(0)\bm{\Gamma}(0)\in\mathbb{R}^{d\times
d}$ (18)

is at most $d_{y}$, which implies that
$\mathrm{dim}\,\mathcal{N}\left(\nabla_{\bm{W}_{1}}\ell(0)\right)=\mathrm{dim}\,\mathcal{N}\left(\nabla_{\bm{W}_{1}}^{\top}\ell(0)\right)\geq
d-d_{y}$. We define the subspace

$\displaystyle\mathcal{S}:=\mathcal{N}\left(\nabla_{\bm{W}_{1}}\ell(0)\right)\cap\mathcal{N}\left(\nabla_{\bm{W}_{1}}^{\top}\ell(0)\bm{W}_{1}(0)\right)\subset\mathbb{R}^{d}.$

Since $\bm{W}_{1}(0)\in\mathbb{R}^{d\times d}$ is nonsingular, we have

$\displaystyle\mathrm{dim}(\mathcal{S})\geq 2(d-d_{y})-d=m.$

Let $\\{\bm{v}_{i}^{(1)}\\}_{i=1}^{m}$ denote an orthonormal basis for
$\mathcal{S}$ and set
$\bm{u}_{i}^{(1)}:=\bm{W}_{1}(0)\bm{v}_{i}^{(1)}/\varepsilon$, where
$\varepsilon>0$ is the constant in (3) – since $\bm{W}_{1}(0)$ is orthogonal,
$\\{\bm{u}_{i}^{(1)}\\}_{i=1}^{m}$ is also an orthonormal set. Then we
trivially have $\bm{W}_{1}(0)\bm{v}_{i}^{(1)}=\varepsilon\bm{u}_{i}^{(1)}$,
which, together with (3), implies
$\bm{W}_{1}^{\top}(0)\bm{u}_{i}^{(1)}=\varepsilon\bm{v}_{i}^{(1)}$. It follows
from $\bm{v}_{i}^{(1)}\in\mathcal{S}$ that
$\nabla_{\bm{W}_{1}}\ell(0)\bm{v}_{i}^{(1)}=\bm{0}$ and
$\nabla_{\bm{W}_{1}}^{\top}\ell(0)\bm{W}_{1}(0)\bm{v}_{i}^{(1)}=\bm{0}$, which
is equivalent to $\bm{W}_{L:2}^{\top}(0)\bm{\Gamma}(0)\bm{v}_{i}^{(1)}=\bm{0}$
and
$\bm{\Gamma}^{\top}(0)\bm{W}_{L:2}(0)\bm{W}_{1}(0)\bm{v}_{i}^{(1)}=\varepsilon\bm{\Gamma}^{\top}(0)\bm{W}_{L:2}(0)\bm{u}_{i}^{(1)}=\bm{0}$
respectively by (18). Since $\bm{W}^{\top}_{L:2}(0)$ and
$\bm{\Gamma}^{\top}(0)$ are full column rank, we have
$\bm{\Gamma}(0)\bm{v}_{i}^{(1)}=\bm{0}$ and
$\bm{W}_{L:2}(0)\bm{u}_{i}^{(1)}=\bm{0}$. It then follows that

$\displaystyle\bm{0}=\bm{\Gamma}(0)\bm{v}_{i}^{(1)}=\left(\bm{W}_{L:1}(0)-\bm{\Phi}\right)\bm{v}_{i}^{(1)}=\varepsilon\bm{W}_{L:2}(0)\bm{u}_{i}^{(1)}-\bm{\Phi}\bm{v}_{i}^{(1)}=-\bm{\Phi}\bm{v}_{i}^{(1)}.$

Therefore, we have $\bm{\Phi}\bm{v}_{i}^{(1)}=\bm{0}$.

Now let $\mathcal{E}(l)$ be the event that we have orthonormal sets
$\\{\bm{u}_{i}^{(l)}\\}_{i=1}^{m}$ and $\\{\bm{v}_{i}^{(l)}\\}_{i=1}^{m}$
satisfying $\bm{W}_{l}(0)\bm{v}_{i}^{(l)}=\varepsilon\bm{u}_{i}^{(l)}$,
$\bm{W}_{l}^{\top}(0)\bm{u}_{i}^{(l)}=\varepsilon\bm{v}_{i}^{(l)}$,
$\bm{W}_{L:l+1}(0)\bm{u}_{i}^{(l)}=\bm{0}$, and
$\bm{\Gamma}(0)\bm{W}_{l-1:1}^{\top}(0)\bm{v}_{i}^{(l)}=\bm{0}$. From the
above arguments, we have that $\mathcal{E}(1)$ holds – now suppose
$\mathcal{E}(k)$ holds for some $1\leq k<L-1$. Set
$\bm{v}_{i}^{(k+1)}:=\bm{u}_{i}^{(k)}$ and
$\bm{u}_{i}^{(k+1)}:=\bm{W}_{k+1}(0)\bm{v}_{i}^{(k+1)}/\varepsilon$. This,
together with (3), implies that
$\bm{W}_{k+1}(0)\bm{v}_{i}^{(k+1)}=\varepsilon\bm{u}_{i}^{(k+1)}$ and
$\bm{W}_{k+1}^{\top}(0)\bm{u}_{i}^{(k+1)}=\varepsilon\bm{v}_{i}^{(k+1)}$.
Moreover, we have

$\displaystyle\bm{W}_{L:(k+1)+1}(0)\bm{u}_{i}^{(k+1)}$
$\displaystyle=\bm{W}_{L:k+1}(0)\bm{W}_{k+1}^{\top}(0)\bm{u}_{i}^{(k+1)}/\varepsilon^{2}$
$\displaystyle=\bm{W}_{L:k+1}(0)\bm{v}_{i}^{(k+1)}/\varepsilon$
$\displaystyle=\bm{W}_{L:k+1}(0)\bm{u}_{i}^{(k)}/\varepsilon=\bm{0},$

where the first two equalities follow from (3) and
$\bm{u}_{i}^{(k+1)}=\bm{W}_{k+1}(0)\bm{v}_{i}^{(k+1)}/\varepsilon$, and the
last equality is due to $\bm{v}_{i}^{(k+1)}=\bm{u}_{i}^{(k)}$. Similarly, we
have

$\displaystyle\bm{\Gamma}(0)\bm{W}_{(k+1)-1:1}^{\top}(0)\bm{v}_{i}^{(k+1)}$
$\displaystyle=\bm{\Gamma}(0)\bm{W}_{k-1:1}^{\top}(0)\bm{W}_{k}^{\top}(0)\bm{v}_{i}^{(k+1)}$
$\displaystyle=\bm{\Gamma}(0)\bm{W}_{k-1:1}^{\top}(0)\bm{W}_{k}^{\top}(0)\bm{u}_{i}^{(k)}$
$\displaystyle=\varepsilon\bm{\Gamma}(0)\bm{W}_{k-1:1}^{\top}(0)\bm{v}_{i}^{(k)}=\bm{0},$

where the second equality follows from $\bm{v}_{i}^{(k+1)}=\bm{u}_{i}^{(k)}$
and the third equality is due to
$\bm{W}_{k}^{\top}(0)\bm{u}_{i}^{(k)}=\varepsilon\bm{v}_{i}^{(k)}$. Therefore
$\mathcal{E}(k+1)$ holds, so we have $\mathcal{E}(l)$ for all $l\in[L-1]$.
Finally, setting $\bm{v}_{i}^{(L)}=\bm{u}_{i}^{(L-1)}$, we have

$\displaystyle\bm{\Gamma}(0)\bm{W}_{L-1:1}^{\top}(0)\bm{v}_{i}^{(L)}$
$\displaystyle=\bm{\Gamma}(0)\bm{W}_{L-2:1}^{\top}(0)\bm{W}_{L-1}^{\top}(0)\bm{v}_{i}^{(L)}$
$\displaystyle=\bm{\Gamma}(0)\bm{W}_{L-1:1}^{\top}(0)\bm{W}_{L-1}^{\top}(0)\bm{u}_{i}^{(L-1)}$
$\displaystyle=\varepsilon\bm{\Gamma}(0)\bm{W}_{L-1:1}^{\top}(0)\bm{v}_{i}^{(L-1)}=\bm{0}.$

As a result, we have shown the base cases $\mathcal{A}(0)$, $\mathcal{B}(0)$,
$\mathcal{C}(0)$, and $\mathcal{D}(0)$.

Now we proceed by induction on $t\geq 0$. Suppose that $\mathcal{A}(t)$,
$\mathcal{B}(t)$, $\mathcal{C}(t)$, and $\mathcal{D}(t)$ hold for some $t\geq
0$. First, we show $\mathcal{A}(t+1)$ and $\mathcal{B}(t+1)$. We have

$\displaystyle\bm{W}_{l}(t+1)\bm{v}_{i}^{(l)}$
$\displaystyle=\left[(1-\eta\lambda)\bm{W}_{l}(t)-\eta\bm{W}_{L:l+1}^{\top}(t){\bm{\Gamma}}(t)\bm{W}_{l-1:1}^{\top}(t)\right]\bm{v}_{i}^{(l)}$
$\displaystyle=(1-\eta\lambda)\bm{W}_{l}(t)\bm{v}_{i}^{(l)}=\rho(t)(1-\eta\lambda)\bm{u}_{i}^{(l)}=\rho(t+1)\bm{u}_{i}^{(l)}$

for all $l\in[L-1]$, where the first equality follows from (4), the second
equality uses $\mathcal{D}(t)$, and the last equality is due to
$\mathcal{A}(t)$, proving $\mathcal{A}(t+1)$. Similarly, we have

$\displaystyle\bm{W}_{l}^{\top}(t+1)\bm{u}_{i}^{(l)}$
$\displaystyle=\left[(1-\eta\lambda)\bm{W}_{l}^{\top}(t)-\eta\bm{W}_{l-1:1}(t)\bm{\Gamma}^{\top}(t)\bm{W}_{L:l+1}(t)\right]\bm{u}_{i}^{(l)}$
$\displaystyle=(1-\eta\lambda)\bm{W}_{l}^{\top}(t)\bm{u}_{i}^{(l)}=\rho(t)(1-\eta\lambda)\bm{v}_{i}^{(l)}=\rho(t+1)\bm{v}_{i}^{(l)}$

for all $l\in[L-1]$, where the second equality uses $\mathcal{C}(t)$, and the
last equality is due to $\mathcal{B}(t)$, proving $\mathcal{B}(t+1)$. Now, we
show $\mathcal{C}(t+1)$. For any $k\in[L-2]$, it follows from
$\bm{v}_{i}^{(k+1)}=\bm{u}_{i}^{(k)}$ and $\mathcal{A}(t+1)$ that

$\bm{W}_{k+1}(t+1)\bm{u}_{i}^{(k)}=\bm{W}_{k+1}(t+1)\bm{v}_{i}^{(k+1)}=\rho(t+1)\bm{u}_{i}^{(k+1)}.$

Repeatedly applying the above equality for $k=l,l+1,\dots,L-2$, we obtain

$\displaystyle\bm{W}_{L:l+1}(t+1)\bm{u}_{i}^{(l)}$
$\displaystyle=\rho(t+1)^{L-l-1}\bm{W}_{L}(t+1)\bm{u}_{i}^{(L-1)}$
$\displaystyle=\rho(t+1)^{L-l-1}\left[(1-\eta\lambda)\bm{W}_{L}(t)-\eta\bm{\Gamma}(t)\bm{W}_{L-1:1}^{\top}(t)\right]\bm{u}_{i}^{(L-1)}=\bm{0}$

where $\bm{W}_{L}(t)\bm{u}_{i}^{(L-1)}=\bm{0}$ follows from $\mathcal{C}(t)$
and $\bm{\Gamma}(t)\bm{W}_{L-1:1}^{\top}(t)\bm{u}_{i}^{(L-1)}=\bm{0}$ follows
from $\bm{u}_{i}^{(L-1)}=\bm{v}_{i}^{(L)}$ and $\mathcal{D}(t)$, proving
$\mathcal{C}(t+1)$. Finally, we show $\mathcal{D}(t+1)$. For any
$k\in\\{2,\dots,L\\}$, it follows from $\bm{v}_{i}^{(k)}=\bm{u}_{i}^{(k-1)}$
and $\mathcal{B}(t+1)$ that

$\displaystyle\bm{W}_{k-1}^{\top}(t+1)\bm{v}_{i}^{(k)}=\bm{W}_{k-1}^{\top}(t+1)\bm{u}_{i}^{(k-1)}=\rho(t+1)\bm{v}_{i}^{(k-1)}.$

Repeatedly applying the above equality for $k=l,l-1,\dots,2$, we obtain

$\displaystyle\bm{\Gamma}(t+1)\bm{W}_{l-1:1}^{\top}(t+1)\bm{v}_{i}^{(l)}$
$\displaystyle=\rho(t+1)^{l-1}\bm{\Gamma}(t+1)\bm{v}_{i}^{(1)}$
$\displaystyle=\rho(t+1)^{l-1}(\bm{W}_{L:1}(t+1)-\bm{\Phi})\bm{v}_{i}^{(1)}=\bm{0}$

where the last line follows from $\bm{\Phi}\bm{v}_{i}^{(1)}=\bm{0}$ as well as

$\bm{W}_{L:1}(t+1)\bm{v}_{i}^{(1)}=\bm{W}_{L:2}(t+1)\bm{W}_{1}(t+1)\bm{v}_{i}^{(1)}=\rho(t+1)\bm{W}_{L:2}(t+1)\bm{u}_{i}^{(1)}=\bm{0}$

by $\mathcal{A}(t+1)$ and $\mathcal{C}(t+1)$. Thus we have proven
$\mathcal{D}(t+1)$, concluding the proof. ∎

###### Proof of Theorem 1.

We show the result for Case 1 – the proof is nearly identical for Case 2.

By $\mathcal{A}(t)$ and $\mathcal{B}(t)$ of Lemma 1, there exists orthonormal
matrices $\\{\bm{U}_{l,2}\\}\subset\mathbb{R}^{d\times m}$ and
$\\{\bm{V}_{l,2}\\}\subset\mathbb{R}^{d\times m}$ for $l\in[L]$ satisfying
$\bm{U}_{l+1,2}=\bm{V}_{l,2}$ for all $l\in[L-1]$ as well as

$\bm{W}_{l}(t)\bm{V}_{l,2}=\rho(t)\bm{U}_{l,2}\quad\mbox{and}\quad\bm{W}_{l}(t)^{\top}\bm{U}_{l,2}=\rho(t)\bm{V}_{l,2}$
(19)

for all $l\in[L]$ and $t\geq 0$, where $\rho(t)$ satisfies (6) for $t\geq 1$
with $\rho(0)=\varepsilon$. First, complete $\bm{V}_{1,2}$ to an orthonormal
basis for $\mathbb{R}^{d}$ as
$\bm{V}_{1}=\begin{bmatrix}\bm{V}_{1,1}&\bm{V}_{1,2}\end{bmatrix}\in\mathcal{O}^{d}$.
Then for each $l\in[L-1]$, set
$\bm{U}_{l}=\begin{bmatrix}\bm{U}_{l,1}&\bm{U}_{l,2}\end{bmatrix}\in\mathcal{O}^{d}$
where $\bm{U}_{l,1}=\bm{W}_{l}(0)\bm{V}_{l,1}/\varepsilon$ and
$\bm{V}_{l+1}=\begin{bmatrix}\bm{V}_{l+1,1}&\bm{V}_{l+1,2}\end{bmatrix}\in\mathcal{O}^{d}$
where $\bm{V}_{l+1,1}=\bm{U}_{l,1}$, and finally set
$\bm{U}_{L}=\begin{bmatrix}\bm{U}_{L,1}&\bm{U}_{L,2}\end{bmatrix}\in\mathcal{O}^{d}$
where $\bm{U}_{L,1}=\bm{W}_{L}(0)\bm{V}_{L,1}/\varepsilon$. We note that
$\bm{V}_{l+1}=\bm{U}_{l}$ for each $l\in[L-1]$. Then we have

$\bm{U}_{l,1}^{\top}\bm{W}_{l}(t)\bm{V}_{l,2}=\rho(t)\bm{U}_{l,1}^{\top}\bm{U}_{l,2}=\bm{0}$
(20)

for all $l\in[L]$, where the first equality follows from (19). Similarly, we
also have

$\bm{U}_{l,2}^{\top}\bm{W}_{l}(t)\bm{V}_{l,1}=\rho(t)\bm{V}_{l,2}^{\top}\bm{V}_{l,1}=\bm{0}$
(21)

for all $l\in[L]$, where the first equality also follows from (19). Therefore,
combining (19), (20), and (21) yields

$\bm{U}_{l}^{\top}\bm{W}_{l}(t)\bm{V}_{l}=\begin{bmatrix}\bm{U}_{l,1}&\bm{U}_{l,2}\end{bmatrix}^{\top}\bm{W}_{l}(t)\begin{bmatrix}\bm{V}_{l+1,1}&\bm{V}_{l+1,2}\end{bmatrix}=\begin{bmatrix}\widetilde{\bm{W}}_{l}(t)&\bm{0}\\\
\bm{0}&\rho(t)\bm{I}_{m}\end{bmatrix}$

for all $l\in[L]$, where $\widetilde{\bm{W}}_{l}(0)=\varepsilon\bm{I}_{2r}$ by
construction of $\bm{U}_{l,1}$. This directly implies (5), completing the
proof. ∎

##  Appendix C Proofs in Section 4

Suppose that $\bm{\Theta}=\\{\bm{W}_{l}\\}_{l=1}^{L}$ satisfies
(i) Global Optimality:

$\displaystyle\bm{W}_{L:1}\bm{X}=\bm{Y}.$ (22)

(ii) Balancedness: There exists a positive constant
$\varepsilon\leq\left\\{\frac{n^{1/2L}}{\sqrt{30}L\sqrt[4]{d-K}},\frac{(n/2)^{1/4L}}{\sqrt[4]{d-K}},\frac{1}{\sqrt{2(\sqrt{K}+1)}}\right\\}$
such that

$\displaystyle\bm{W}_{l+1}^{\top}\bm{W}_{l+1}=\bm{W}_{l}\bm{W}_{l}^{\top},\forall
l\in[L-2],\
\|\bm{W}_{L}^{\top}\bm{W}_{L}-\bm{W}_{L-1}\bm{W}_{L-1}^{\top}\|_{F}\leq\varepsilon^{2}\sqrt{d-K}.$
(23)

(iii) Unchanged Spectrum: There exists a positive constant $\varepsilon>0$ and
an index set $\mathcal{A}\subseteq[d]$ with $|\mathcal{A}|=d-2K$ such that for
all $l\in[L-1]$ that

$\displaystyle\sigma_{i}(\bm{W}_{l})=\varepsilon,\ \forall i\in\mathcal{A}.$
(24)

For ease of exposition, we introduce some additional notation. In our
analysis, we can assume $\bm{Y}=\bm{I}_{K}\otimes\bm{1}_{n}^{\top}$ without
loss of generality. This, together with (22) and $\bm{X}$ is full rank, yields
that the rank of $\bm{W}_{l}$ is at least $K$ for all $l\in[L]$. Let

$\displaystyle\bm{W}_{l}=\bm{U}_{l}\bm{\Sigma}_{l}\bm{V}_{l}^{\top}=\begin{bmatrix}\bm{U}_{l,1}&\bm{U}_{l,2}\end{bmatrix}\begin{bmatrix}\bm{\Sigma}_{l,1}&\bm{0}\\\
\bm{0}&\bm{\Sigma}_{l,2}\end{bmatrix}\begin{bmatrix}\bm{V}_{l,1}^{\top}\\\
\bm{V}_{l,2}^{\top}\end{bmatrix}=\bm{U}_{l,1}\bm{\Sigma}_{l,1}\bm{V}_{l,1}^{\top}+\bm{U}_{l,2}\bm{\Sigma}_{l,2}\bm{V}_{l,2}^{\top},$
(25)

be a singular value decomposition (SVD) of $\bm{W}_{l}$, where
$\bm{\Sigma}_{l}\in\mathbb{R}^{d\times d}$ is diagonal,
$\bm{\Sigma}_{l,1}=\operatorname{diag}\left(\sigma_{l,1},\dots,\sigma_{l,K}\right)$
with $\sigma_{l,1}\geq\cdots\geq\sigma_{l,K}>0$ being the singular values, and
$\bm{\Sigma}_{l,2}=\operatorname{diag}\left(\sigma_{l,K+1},\dots,\sigma_{l,d}\right)$
with $\sigma_{l,K+1}\geq\cdots\geq\sigma_{l,d}\geq 0$ being the remaining
singular values of $\bm{W}_{l}$; $\bm{U}_{l}\in\mathcal{O}^{d}$ with
$\bm{U}_{l,1}\in\mathbb{R}^{d\times K}$,
$\bm{U}_{l,2}\in\mathbb{R}^{d\times(d-K)}$; $\bm{V}_{l}\in\mathcal{O}^{d}$
with $\bm{V}_{l,1}\in\mathbb{R}^{d\times K}$,
$\bm{V}_{l,2}\in\mathbb{R}^{d\times(d-K)}$. Noting that
$\bm{W}_{L}\in\mathbb{R}^{K\times d}$, let

$\displaystyle\bm{W}_{L}=\bm{U}_{L}\bm{\Sigma}_{L}\bm{V}_{L}^{\top}=\bm{U}_{L}\begin{bmatrix}\bm{\Sigma}_{L,1}&\bm{0}\end{bmatrix}\begin{bmatrix}\bm{V}_{L,1}^{\top}\\\
\bm{V}_{L,2}^{\top}\end{bmatrix}=\bm{U}_{L}\bm{\Sigma}_{L,1}\bm{V}_{L,1}^{\top}$
(26)

be its singular value decomposition (SVD), where
$\bm{\Sigma}_{L}\in\mathbb{R}^{K\times d}$, and
$\bm{\Sigma}_{L,1}=\operatorname{diag}\left(\sigma_{L,1},\dots,\sigma_{L,K}\right)$
with $\sigma_{L,1}\geq\cdots\geq\sigma_{L,K}>0$ being the singular values;
$\bm{U}_{L}\in\mathcal{O}^{K}$, and $\bm{V}_{L}\in\mathcal{O}^{d}$ with
$\bm{V}_{L,1}\in\mathbb{R}^{d\times K}$,
$\bm{V}_{L,2}\in\mathbb{R}^{d\times(d-K)}$.

Before we prove Theorem 2, we need some preliminary results. We first show
that when $\\{\bm{W}_{l}\\}_{l=1}^{L}$ satisfy (22) and (23), we can bound the
leading $K$ singular values of $\bm{W}_{l}$ for all $l\in[L]$. Throughout this
section, let

$\displaystyle\delta:=\varepsilon^{2}\sqrt{d-K}\leq\frac{n^{1/L}}{30L^{2}},$
(27)

where the inequality follows from

$\displaystyle\varepsilon\leq\min\left\\{\frac{n^{1/2L}}{\sqrt{30}L\sqrt[4]{d-K}},\frac{(n/2)^{1/4L}}{\sqrt[4]{d-K}},\frac{1}{\sqrt{2(\sqrt{K}+1)}}\right\\}.$
(28)

######  Lemma 3.

Suppose that the weights $\\{\bm{W}_{l}\\}_{l=1}^{L}$ satisfy (22) and (23).
Then, it holds that

$\displaystyle\left(\sqrt{\frac{n}{2}}\right)^{1/L}\leq\sigma_{K}\left(\bm{W}_{l}\right)\leq\sigma_{1}\left(\bm{W}_{l}\right)\leq\left(\sqrt{2n}\right)^{1/L},\
\forall l\in[L].$ (29)

###### Proof of Lemma 3.

It follows from (22) and the fact that $\bm{X}\in\mathbb{R}^{d\times d}$ is
orthogonal that $\bm{W}_{L:1}=\bm{Y}\bm{X}^{\top}$. This, together with
$\bm{Y}=\bm{I}_{K}\otimes\bm{1}^{\top}_{n}$ and $\bm{X}\in\mathcal{O}^{d}$,
yields that $\bm{W}_{L:1}$ is of rank $K$ and

$\displaystyle\sigma_{i}(\bm{W}_{L:1})=\sqrt{n},\ \forall i\in[K].$ (30)

This, together with [55, Lemma 6] and (27), implies
$\sigma_{1}(\bm{W}_{l})\leq\left(\sqrt{2n}\right)^{1/L}$. Using
$\bm{W}_{l+1}^{T}\bm{W}_{l+1}=\bm{W}_{l}\bm{W}_{l}^{T}$ for all $l\in[L-2]$ in
(23), we obtain

$\displaystyle\|\bm{W}_{L:1}^{\top}\bm{W}_{L:1}-\left(\bm{W}_{1}^{\top}\bm{W}_{1}\right)^{L}\|_{F}$
$\displaystyle=\|\bm{W}_{L:1}^{\top}\bm{W}_{L:1}-\bm{W}_{L-1:1}^{\top}\bm{W}_{L-1}\bm{W}_{L-1}^{\top}\bm{W}_{L-1:1}\|_{F}$
$\displaystyle=\|\bm{W}_{L-1:1}^{\top}\left(\bm{W}_{L}^{\top}\bm{W}_{L}-\bm{W}_{L-1}\bm{W}_{L-1}^{\top}\right)\bm{W}_{L-1:1}\|_{F}$
$\displaystyle\leq\|\bm{W}_{L}^{\top}\bm{W}_{L}-\bm{W}_{L-1}\bm{W}_{L-1}^{\top}\|_{F}\prod_{l=1}^{L-1}\|\bm{W}_{l}\|^{2}$
$\displaystyle\leq\delta\left(2n\right)^{(L-1)/L}\leq\frac{n}{15L^{2}},$

where the second inequality uses (23) and
$\|\bm{W}_{l}\|\leq(\sqrt{2n})^{1/L}$ for all $l\in[L-1]$, and the last
inequality follows from (27). Using (30) and Weyl’s inequality, we obtain

$\displaystyle\sigma_{K}\left(\left(\bm{W}_{1}^{\top}\bm{W}_{1}\right)^{L}\right)$
$\displaystyle\geq\sigma_{K}\left(\bm{W}_{L:1}^{\top}\bm{W}_{L:1}\right)-\|\bm{W}_{L:1}^{\top}\bm{W}_{L:1}-\left(\bm{W}_{1}^{\top}\bm{W}_{1}\right)^{L}\|$
$\displaystyle\geq
n-\|\bm{W}_{L:1}^{\top}\bm{W}_{L:1}-\left(\bm{W}_{1}^{\top}\bm{W}_{1}\right)^{L}\|_{F}\geq\left(1-\frac{1}{15L^{2}}\right)n.$

Therefore, we have
$\sigma_{K}(\bm{W}_{1})\geq\left(\left(1-\frac{1}{15L^{2}}\right)n\right)^{1/2L}$.
Using this and (23), we obtain
$\sigma_{K}(\bm{W}_{l})\geq\left(\left(1-\frac{1}{15L^{2}}\right)n\right)^{1/2L}$
for all $l\in[L-1]$. This, together with Weyl’s inequality and
$\|\bm{W}_{L}^{T}\bm{W}_{L}-\bm{W}_{L-1}\bm{W}_{L-1}^{T}\|_{F}\leq\delta$ in
(23), yields

$\displaystyle\sigma_{K}\left(\bm{W}_{L}^{T}\bm{W}_{L}\right)$
$\displaystyle\geq\sigma_{K}\left(\bm{W}_{L-1}^{T}\bm{W}_{L-1}\right)-\|\bm{W}_{L}^{T}\bm{W}_{L}-\bm{W}_{L-1}\bm{W}_{L-1}^{T}\|$
$\displaystyle\geq\left(1-\frac{1}{15L^{2}}\right)^{1/L}n^{1/L}-\delta\geq\left(\frac{n}{2}\right)^{1/L},$

where the last inequality follows from

$\displaystyle
1-\frac{1}{15L^{2}}\geq\left(\frac{1}{30L^{2}}+\left(\frac{1}{2}\right)^{1/L}\right)^{L},\
\forall L\geq 1.$

Since $1-1/15L^{2}\geq 1/2$ for all $L\geq 1$, we also have
$\sigma_{K}(\bm{W}_{l})\geq\left(\frac{n}{2}\right)^{1/2L}$. Combining all
this together yields (29). ∎

Then, we show that if the weight matrices $\\{\bm{W}_{l}\\}_{l=1}^{L}$ satisfy
(22) and (23), the right singular vectors of $\bm{W}_{l}$ are equal or close
to the left singular vectors of $\bm{W}_{l+1}$ for all $l\in[L-1]$.

######  Lemma 4.

Suppose that $\\{\bm{W}_{l}\\}_{l=1}^{L}$ satisfy (22) and (23) with $\delta$
satisfying (27), and admit the SVD in (25) and (26). Then, it holds that

$\displaystyle\bm{\Sigma}_{l+1}=\bm{\Sigma}_{l},\
\|\bm{\Sigma}_{l,2}^{2}\|_{F}\leq\delta,\ \forall l\in[L-1],\
\bm{V}_{l+1}^{\top}\bm{U}_{l}=\bm{I}_{d},\ \forall l\in[L-2],$ (31)
$\displaystyle\|\bm{V}_{L,1}^{T}\bm{U}_{L-1,2}\|_{F}\leq\frac{2\sqrt{\delta}\sqrt[4]{K}}{n^{1/2L}},\
\sigma_{\min}(\bm{V}_{L,1}^{\top}\bm{U}_{L-1,1})\geq
1-\frac{2\sqrt{\delta}\sqrt[4]{K}}{n^{1/2L}}.$ (32)

###### Proof.

It follows from $\bm{W}_{l+1}^{\top}\bm{W}_{l+1}=\bm{W}_{l}\bm{W}_{l}^{\top}$
for all $l\in[L-2]$ and (25) that

$\displaystyle\bm{V}_{l+1}\bm{\Sigma}_{l+1}^{\top}\bm{\Sigma}_{l+1}\bm{V}_{l+1}^{\top}=\bm{U}_{l}\bm{\Sigma}_{l}\bm{\Sigma}_{l}^{\top}\bm{U}_{l}^{\top}.$

This, together with the argument in [17, Proof of Theorem 1], implies (31).
Using
$\|\bm{W}_{L}^{\top}\bm{W}_{L}-\bm{W}_{L-1}\bm{W}_{L-1}^{\top}\|_{F}\leq\delta$,
(25), and (26), we have

$\displaystyle\delta$
$\displaystyle\geq\|\bm{V}_{L}\bm{\Sigma}_{L}^{\top}\bm{\Sigma}_{L}\bm{V}_{L}^{\top}-\bm{U}_{L-1}\bm{\Sigma}_{L-1}\bm{\Sigma}_{L-1}^{\top}\bm{U}_{L-1}^{\top}\|_{F}$
$\displaystyle=\|\bm{V}_{L,1}\bm{\Sigma}^{2}_{L,1}\bm{V}_{L,1}^{T}-\bm{U}_{L-1,1}\bm{\Sigma}^{2}_{L-1,1}\bm{U}_{L-1,1}^{T}-\bm{U}_{L-1,2}\bm{\Sigma}_{L-1,2}^{2}\bm{U}_{L-1,2}^{T}\|_{F}.$
(33)

Obviously, we have
$\|\bm{\Sigma}_{L}^{T}\bm{\Sigma}_{L}-\bm{V}_{L}^{T}\bm{U}_{L-1}\bm{\Sigma}_{L-1}^{2}\bm{U}_{L-1}^{T}\bm{V}_{L}\|_{F}\leq\delta$
due to $\bm{V}_{L}\in\mathcal{O}^{d}$. This, together with [55, Lemma 4],
implies

$\displaystyle\|\bm{\Sigma}_{L}^{T}\bm{\Sigma}_{L}-\bm{\Sigma}_{L-1}^{2}\|_{F}\leq\delta.$

This, together with (31), implies $\|\bm{\Sigma}_{l,2}^{2}\|_{F}\leq\delta$
for all $l\in[L-1]$. Using this and the structures of $\bm{\Sigma}_{L}$ and
$\bm{\Sigma}_{L-1}$, we obtain

$\displaystyle\|\bm{\Sigma}_{L,1}^{2}-\bm{\Sigma}_{L-1,1}^{2}\|_{F}\leq\delta,\
\|\bm{\Sigma}_{L-1,2}^{2}\|_{F}\leq\delta.$ (34)

Using (33), $\|\bm{U}^{\top}\bm{A}\bm{U}\|_{F}\leq\|\bm{A}\|_{F}$ for any
$\bm{U}\in\mathcal{O}^{d\times(d-K)}$, we further obtain

$\displaystyle\|\bm{U}_{L-1,2}^{T}\bm{V}_{L,1}\bm{\Sigma}^{2}_{L,1}\bm{V}_{L,1}^{T}\bm{U}_{L-1,2}-\bm{\Sigma}_{L-1,2}^{2}\|_{F}\leq\delta.$

This, together with (34), yields

$\displaystyle\|\bm{U}_{L-1,2}^{T}\bm{V}_{L,1}\bm{\Sigma}^{2}_{L,1}\bm{V}_{L,1}^{T}\bm{U}_{L-1,2}\|_{F}\leq
2\delta.$

It follows from this and Lemma 5 that

$\displaystyle\|\bm{\Sigma}_{L,1}\bm{V}_{L,1}^{T}\bm{U}_{L-1,2}\|_{F}\leq\sqrt{2\delta}\sqrt[4]{K}.$
(35)

Noting that
$\|\bm{\Sigma}_{L,1}\bm{V}_{L,1}^{T}\bm{U}_{L-1,2}\|_{F}\geq\sigma_{\min}(\bm{\Sigma}_{L,1})\|\bm{V}_{L,1}^{T}\bm{U}_{L-1,2}\|_{F}$,
we have

$\displaystyle\|\bm{V}_{L,1}^{T}\bm{U}_{L-1,2}\|_{F}\leq\frac{\|\bm{\Sigma}_{L,1}\bm{V}_{L,1}^{T}\bm{U}_{L-1,2}\|_{F}}{\sigma_{\min}(\bm{\Sigma}_{L,1})}\leq\frac{\sqrt{2\delta}\sqrt[4]{K}}{\left(\frac{n}{2}\right)^{1/2L}}\leq\frac{2\sqrt{\delta}\sqrt[4]{K}}{n^{1/2L}},$
(36)

where the second inequality follows from (35) and Lemma 3. Then, we compute

$\displaystyle\sigma_{\min}^{2}(\bm{V}_{L,1}^{T}\bm{U}_{L-1,1})$
$\displaystyle=\min_{\|\bm{x}\|=1}\|\bm{U}_{L-1,1}^{T}\bm{V}_{L,1}\bm{x}\|^{2}=\min_{\|\bm{x}\|=1}\bm{x}^{T}\bm{V}_{L,1}^{T}\bm{U}_{L-1,1}\bm{U}_{L-1,1}^{T}\bm{V}_{L,1}\bm{x}$
$\displaystyle=\min_{\|\bm{x}\|=1}\bm{x}^{T}\bm{V}_{L,1}^{T}\left(\bm{I}-\bm{U}_{L-1,2}\bm{U}_{L-1,2}^{T}\right)\bm{V}_{L,1}\bm{x}=1-\max_{\|\bm{x}\|=1}\|\bm{U}_{L-1,2}^{T}\bm{V}_{L,1}\bm{x}\|^{2}$
$\displaystyle\geq 1-\sigma_{\max}^{2}(\bm{U}_{L-1,2}^{T}\bm{V}_{L,1})\geq
1-\|\bm{U}_{L-1,2}^{T}\bm{V}_{L,1}\|_{F}^{2},$

where the third equality follows from $\bm{U}_{L-1}\in\mathcal{O}^{d}$. This,
together with (36), implies

$\displaystyle\sigma_{\min}(\bm{V}_{L,1}^{T}\bm{U}_{L-1,1})\geq\sqrt{1-\|\bm{U}_{L-1,2}^{T}\bm{V}_{L,1}\|_{F}^{2}}\geq
1-\|\bm{U}_{L-1,2}^{T}\bm{V}_{L,1}\|_{F}\geq
1-\frac{4\delta\sqrt{K}}{n^{1/2L}}.$

∎

Recall that $\bar{\bm{x}}_{k}=\sum_{i=1}^{n}\bm{x}_{k,i}/n$ and
$\bar{\bm{x}}=\sum_{k=1}^{K}\bar{\bm{x}}_{k}/K$. Let

$\displaystyle\bm{\Delta}_{W}=\left[\bm{\delta}_{1,1},\dots,\bm{\delta}_{1,n},\dots,\bm{\delta}_{K,1},\dots,\bm{\delta}_{K,n}\right]\in\mathbb{R}^{d\times
N},\ \text{where}\ \bm{\delta}_{k,i}=\bm{x}_{k,i}-\bar{\bm{x}}_{k},\ \forall
k,i,$ (37)
$\displaystyle\bm{\Delta}_{B}=\left[\bar{\bm{\delta}}_{1},\dots,\bar{\bm{\delta}}_{K}\right]\in\mathbb{R}^{d\times
K},\ \text{where}\ \bar{\bm{\delta}}_{k}=\bar{\bm{x}}_{k}-\bar{\bm{x}},\
\forall k\in[K].$ (38)

This, together with (12), (14), and (13), yields

$\displaystyle
D_{0}=\frac{K}{N}\frac{\|\bm{\Delta}_{W}\|_{F}^{2}}{\|\bm{\Delta}_{B}\|_{F}^{2}},\quad
D_{l}=\frac{K}{N}\frac{\|\bm{W}_{l:1}\bm{\Delta}_{W}\|_{F}^{2}}{\|\bm{W}_{l:1}\bm{\Delta}_{B}\|_{F}^{2}}.$
(39)

###### Proof of Theorem 2.

Using (31) in Lemma 4, there exists a diagonal matrix
$\tilde{\bm{\Sigma}}=\mathrm{diag}(\sigma_{1},\dots,\sigma_{d})$ with
$\sigma_{1}\geq\dots\geq\sigma_{d}$ such that

$\displaystyle\bm{\Sigma}_{l}=\tilde{\bm{\Sigma}},\ \forall l\in[L-1].$ (40)

According to Lemma 3, we have
$\sigma_{K}(\bm{W}_{l})\geq\left(\sqrt{n/2}\right)^{1/L}\geq\delta$, where the
last inequality follows from (27) and (28). It follows from this and (24) that
$\mathcal{A}\subseteq\\{K+1,\dots,d\\}$. For simplicity, we write
$\tilde{\bm{\Sigma}}=\begin{bmatrix}\tilde{\bm{\Sigma}}_{1}&\bm{0}\\\
\bm{0}&\tilde{\bm{\Sigma}}_{2}\end{bmatrix}$ satisfying
$\tilde{\bm{\Sigma}}_{1}=\operatorname{diag}\left(\sigma_{1},\dots,\sigma_{K}\right)$
and
$\tilde{\bm{\Sigma}}_{2}=\operatorname{diag}\left(\sigma_{K+1},\dots,\sigma_{d}\right)$.
Using (32) and (24), we have

$\displaystyle\sum_{i=K+1}^{d}\sigma_{i}^{4}\leq\delta^{2}=\varepsilon^{4}(d-K).$

Using (24) and letting
$\mathcal{A}^{c}=\\{K+1,\dots,d\\}\setminus\mathcal{A}$, we have

$\displaystyle\sum_{i\in\mathcal{A}^{c}}\sigma_{i}^{4}+(d-2K)\varepsilon^{4}\leq\varepsilon^{4}(d-K),$

which implies $\sigma_{i}\leq\sqrt[4]{K}\varepsilon$ for all
$i\in\mathcal{A}^{c}$. Therefore, we have

$\displaystyle\sigma_{i}\leq\sqrt[4]{K}\varepsilon,\ i=K+1,\dots,d.$ (41)

It follows from (31) in Lemma 4 that

$\displaystyle\bm{W}_{l}\dots\bm{W}_{1}=\bm{U}_{l}\tilde{\bm{\Sigma}}^{l}\bm{V}_{1}^{T},\
\forall l\in[L-1].$ (42)

It follows from (22), (37), and $\bm{Y}=\bm{I}_{K}\otimes\bm{1}^{T}_{n}$ that

$\displaystyle\bm{W}_{L:1}\bm{\Delta}_{W}$
$\displaystyle=\bm{W}_{L:1}\begin{bmatrix}\bm{x}_{1,1}&\dots&\bm{x}_{1,n}&\dots&\bm{x}_{K,1}&\dots&\bm{x}_{K,n}\end{bmatrix}$
$\displaystyle\quad-\bm{W}_{L:1}\begin{bmatrix}\bar{\bm{x}}_{1}&\dots&\bar{\bm{x}}_{1}&\dots&\bar{\bm{x}}_{K}&\dots&\bar{\bm{x}}_{K}\end{bmatrix}=\bm{0}.$
(43)

Substituting (25), (31), and (40) into (C) yields

$\displaystyle\bm{0}=\bm{W}_{L}\bm{U}_{L-1}\tilde{\bm{\Sigma}}^{L-1}\bm{V}_{1}^{T}\bm{\Delta}_{W}=\bm{W}_{L}\left(\bm{U}_{L-1,1}\tilde{\bm{\Sigma}}_{1}^{L-1}\bm{V}_{1,1}^{T}+\bm{U}_{L-1,2}\tilde{\bm{\Sigma}}_{2}^{L-1}\bm{V}_{1,2}^{T}\right)\bm{\Delta}_{W}.$

As a result, we obtain for all $l\in[L-1]$,

$\displaystyle\|\bm{W}_{L}\bm{U}_{L-1,1}\tilde{\bm{\Sigma}}_{1}^{L-1}\bm{V}_{1,1}^{T}\bm{\Delta}_{W}\|_{F}$
$\displaystyle=\|\bm{U}_{L}\bm{\Sigma}_{L,1}\bm{V}_{L,1}^{T}\bm{U}_{L-1,2}\tilde{\bm{\Sigma}}_{2}^{L-1}\bm{V}_{1,2}^{T}\bm{\Delta}_{W}\|_{F}$
$\displaystyle\leq\|\bm{\Sigma}_{L,1}\|\|\bm{V}_{L,1}^{T}\bm{U}_{L-1,2}\|\|\tilde{\bm{\Sigma}}_{2}^{L-1-l}\|\|\tilde{\bm{\Sigma}}_{2}^{l}\bm{V}_{1,2}^{T}\bm{\Delta}_{W}\|_{F}$
$\displaystyle\leq(2n)^{1/2L}\frac{2\sqrt{\delta}\sqrt[4]{K}}{n^{1/2L}}\left(\sqrt[4]{K}\varepsilon\right)^{L-1-l}\|\tilde{\bm{\Sigma}}_{2}^{l}\bm{V}_{1,2}^{T}\bm{\Delta}_{W}\|_{F}$
$\displaystyle\leq
2\sqrt{2}K^{(L-l)/4}\sqrt[4]{d-K}\varepsilon^{L-l}\|\tilde{\bm{\Sigma}}_{2}^{l}\bm{V}_{1,2}^{T}\bm{\Delta}_{W}\|_{F},$
(44)

where the second inequality follows from (29), (32), and (41), and the last
inequality uses (27). Next, we compute

$\displaystyle\sigma_{\min}(\bm{W}_{L}\bm{U}_{L-1,1})$
$\displaystyle=\sigma_{\min}(\bm{U}_{L}\bm{\Sigma}_{L,1}\bm{V}_{L,1}^{T}\bm{U}_{L-1,1})=\sigma_{\min}(\bm{\Sigma}_{L,1}\bm{V}_{L,1}^{T}\bm{U}_{L-1,1})$
$\displaystyle\geq\sigma_{\min}(\bm{V}_{L,1}^{T}\bm{U}_{L-1,1})\sigma_{\min}(\bm{\Sigma}_{L,1})\geq\left(1-\frac{2\sqrt{\delta}\sqrt[4]{K}}{n^{1/2L}}\right)\left(\frac{n}{2}\right)^{1/2L}$
$\displaystyle\geq\frac{1}{2}\left(\frac{n}{2}\right)^{1/2L},$ (45)

where the second inequality follows from (29) and (32), and the last
inequality uses (27) and (28). Then, we compute

$\displaystyle\|\bm{W}_{L}\bm{U}_{L-1,1}\tilde{\bm{\Sigma}}_{1}^{L-1}\bm{V}_{1,1}^{T}\bm{\Delta}_{W}\|_{F}$
$\displaystyle\geq\sigma_{\min}(\bm{W}_{L}\bm{U}_{L-1,1}\tilde{\bm{\Sigma}}_{1}^{L-l-2})\|\tilde{\bm{\Sigma}}_{1}^{l+1}\bm{V}_{1,1}^{T}\bm{\Delta}_{W}\|_{F}$
$\displaystyle\geq\sigma_{\min}(\bm{W}_{L}\bm{U}_{L-1,1})\sigma_{\min}(\tilde{\bm{\Sigma}}_{1}^{L-l-2})\|\tilde{\bm{\Sigma}}_{1}^{l+1}\bm{V}_{1,1}^{T}\bm{\Delta}_{W}\|_{F}$
$\displaystyle\geq\frac{1}{2}\left(\frac{n}{2}\right)^{(L-l-1)/2L}\|\tilde{\bm{\Sigma}}_{1}^{l+1}\bm{V}_{1,1}^{T}\bm{\Delta}_{W}\|_{F},$
(46)

where the last inequality follows from (29) and (C). This, together with (C),
implies for $l\in[L-2]$,

$\displaystyle\frac{\|\tilde{\bm{\Sigma}}_{1}^{l+1}\bm{V}_{1,1}^{T}\bm{\Delta}_{W}\|_{F}}{\|\tilde{\bm{\Sigma}}_{2}^{l}\bm{V}_{1,2}^{T}\bm{\Delta}_{W}\|_{F}}\leq\kappa_{1}\varepsilon^{L-l},\
\text{where}\
\kappa_{1}:=\frac{4\sqrt{2}K^{(L-l)/4}\sqrt[4]{d-K}}{(n/2)^{(L-l-1)/2L}}.$
(47)

It follows from (25), (40), and (42) that for all $l\in[L-1]$,

$\displaystyle\|\bm{W}_{l:1}\bm{\Delta}_{W}\|_{F}^{2}=\|\bm{U}_{l}\tilde{\bm{\Sigma}}^{l}\bm{V}_{1}^{T}\bm{\Delta}_{W}\|_{F}^{2}=\|\tilde{\bm{\Sigma}}^{l}_{1}\bm{V}_{1,1}^{T}\bm{\Delta}_{W}\|_{F}^{2}+\|\tilde{\bm{\Sigma}}^{l}_{2}\bm{V}_{1,2}^{T}\bm{\Delta}_{W}\|_{F}^{2}.$

This, together with (41), yields

$\displaystyle\frac{\|\bm{W}_{l+1:1}\bm{\Delta}_{W}\|_{F}^{2}}{\|\bm{W}_{l:1}\bm{\Delta}_{W}\|_{F}^{2}}$
$\displaystyle=\frac{\|\tilde{\bm{\Sigma}}^{l+1}_{1}\bm{V}_{1,1}^{T}\bm{\Delta}_{W}\|_{F}^{2}+\|\tilde{\bm{\Sigma}}^{l+1}_{2}\bm{V}_{1,2}^{T}\bm{\Delta}_{W}\|_{F}^{2}}{\|\tilde{\bm{\Sigma}}^{l}_{1}\bm{V}_{1,1}^{T}\bm{\Delta}_{W}\|_{F}^{2}+\|\tilde{\bm{\Sigma}}^{l}_{2}\bm{V}_{1,2}^{T}\bm{\Delta}_{W}\|_{F}^{2}}$
$\displaystyle\leq\frac{\kappa_{1}\varepsilon^{L-l}\|\tilde{\bm{\Sigma}}^{l}_{2}\bm{V}_{1,2}^{T}\bm{\Delta}_{W}\|_{F}+\sqrt{K}\varepsilon^{2}\|\tilde{\bm{\Sigma}}^{l}_{2}\bm{V}_{1,2}^{T}\bm{\Delta}_{W}\|_{F}^{2}}{\|\tilde{\bm{\Sigma}}^{l}_{2}\bm{V}_{1,2}^{T}\bm{\Delta}_{W}\|_{F}^{2}}$
$\displaystyle\leq\kappa_{1}\varepsilon^{L-l}+\sqrt{K}\varepsilon^{2}\leq(\sqrt{K}+1)\varepsilon^{2},$
(48)

where the first inequality follows from (41) and (47). Suppose that we have

$\displaystyle\|\bm{V}_{1,2}^{T}\bm{\Delta}_{B}\|_{F}\leq\|\bm{V}_{1,1}^{T}\bm{\Delta}_{B}\|_{F}.$
(49)

Then, we have

$\displaystyle\|\tilde{\bm{\Sigma}}^{l+1}_{2}\bm{V}_{1,2}^{T}\bm{\Delta}_{B}\|_{F}\leq\|\bm{V}_{1,2}^{T}\bm{\Delta}_{B}\|_{F}\leq\|\bm{V}_{1,1}^{T}\bm{\Delta}_{B}\|_{F}\leq\|\tilde{\bm{\Sigma}}^{l+1}_{1}\bm{V}_{1,2}^{T}\bm{\Delta}_{B}\|_{F},$

where the first inequality follows from $\sigma_{i}\leq 1$ for all
$i=K+1,\dots,d$ due to (41) and (28), the second inequality uses (49), and the
last inequality holds by (29). Using (42) and a similar argument as above, we
compute

$\displaystyle\frac{\|\bm{W}_{l:1}\bm{\Delta}_{B}\|_{F}^{2}}{\|\bm{W}_{l+1:1}\bm{\Delta}_{B}\|_{F}^{2}}$
$\displaystyle=\frac{\|\tilde{\bm{\Sigma}}^{l}_{1}\bm{V}_{1,1}^{T}\bm{\Delta}_{B}\|_{F}^{2}+\|\tilde{\bm{\Sigma}}^{l}_{2}\bm{V}_{1,2}^{T}\bm{\Delta}_{B}\|_{F}^{2}}{\|\tilde{\bm{\Sigma}}^{l+1}_{1}\bm{V}_{1,1}^{T}\bm{\Delta}_{B}\|_{F}^{2}+\|\tilde{\bm{\Sigma}}^{l+1}_{2}\bm{V}_{1,2}^{T}\bm{\Delta}_{B}\|_{F}^{2}}$
$\displaystyle\leq\frac{\|\tilde{\bm{\Sigma}}^{l}_{1}\bm{V}_{1,1}^{T}\bm{\Delta}_{B}\|_{F}^{2}+\|\tilde{\bm{\Sigma}}^{l}_{2}\bm{V}_{1,2}^{T}\bm{\Delta}_{B}\|_{F}^{2}}{\|\tilde{\bm{\Sigma}}^{l+1}_{1}\bm{V}_{1,1}^{T}\bm{\Delta}_{B}\|_{F}^{2}}\leq\frac{2\|\tilde{\bm{\Sigma}}^{l}_{1}\bm{V}_{1,1}^{T}\bm{\Delta}_{B}\|_{F}^{2}}{\|\tilde{\bm{\Sigma}}^{l+1}_{1}\bm{V}_{1,1}^{T}\bm{\Delta}_{B}\|_{F}^{2}}$
$\displaystyle\leq\frac{2\|\tilde{\bm{\Sigma}}^{l}_{1}\bm{V}_{1,1}^{T}\bm{\Delta}_{B}\|_{F}^{2}}{\sigma_{\min}(\tilde{\bm{\Sigma}})\|\tilde{\bm{\Sigma}}^{l}_{1}\bm{V}_{1,1}^{T}\bm{\Delta}_{B}\|_{F}^{2}}\leq
2,$

where the last inequality uses $\sigma_{\min}(\tilde{\bm{\Sigma}})\geq 1$ due
to (29). This, together with (39) and (C), implies (15).

The rest of the proof is devoted to showing (49). For ease of exposition, let
$\bm{V}:=\bm{X}\bm{Y}^{T}/\sqrt{n}\in\mathbb{R}^{d\times K}$. One can verify
$\bm{V}^{T}\bm{V}=\bm{I}_{K}$. In addition, we can compute

$\displaystyle\bm{V}\bm{V}^{T}\bm{\Delta}_{B}=\bm{\Delta}_{B},$ (50)

where the equality follows from (38) and $\bm{X}\in\mathcal{O}^{d}$. Using
(22) and $\bm{X}\in\mathcal{O}^{d}$, we have

$\displaystyle\bm{W}_{L:1}=\bm{Y}\bm{X}^{T}=\sqrt{n}\bm{V}^{T}.$ (51)

This, together with (26) and (42), yields

$\displaystyle\bm{U}_{L}\bm{\Sigma}_{L,1}\bm{V}_{L,1}^{T}\bm{U}_{L-1}\tilde{\bm{\Sigma}}^{L-1}\bm{V}_{1}^{T}=\sqrt{n}\bm{V}^{T}.$

Therefore, we obtain

$\displaystyle n\bm{V}\bm{V}^{T}$
$\displaystyle=\bm{V}_{1}\tilde{\bm{\Sigma}}^{L-1}\bm{U}_{L-1}^{T}\bm{V}_{L,1}\bm{\Sigma}_{L,1}^{2}\bm{V}_{L,1}^{T}\bm{U}_{L-1}\tilde{\bm{\Sigma}}^{L-1}\bm{V}_{1}^{T}$
$\displaystyle=\left(\bm{V}_{1,1}\tilde{\bm{\Sigma}}_{1}^{L-1}\bm{U}_{L-1,1}^{T}+\bm{V}_{1,2}\tilde{\bm{\Sigma}}_{2}^{L-1}\bm{U}_{L-1,2}^{T}\right)\bm{V}_{L,1}\bm{\Sigma}_{L,1}^{2}\bm{V}_{L,1}^{T}\left(\bm{U}_{L-1,1}\tilde{\bm{\Sigma}}_{1}^{L-1}\bm{V}_{1,1}^{T}+\bm{U}_{L-1,2}\tilde{\bm{\Sigma}}_{2}^{L-1}\bm{V}_{1,2}^{T}\right).$

Using (32) and Davis-Kahan Theorem [77, Theorem V.3.6], we have

$\displaystyle\|\bm{V}\bm{V}^{T}-\bm{V}_{1,1}\bm{V}_{1,1}^{T}\|_{F}\leq\frac{1}{2}.$
(52)

Then, we compute

$\displaystyle\|\bm{V}_{1,2}^{T}\bm{\Delta}_{B}\|_{F}$
$\displaystyle=\|(\bm{I}-\bm{V}_{1,1}\bm{V}_{1,1}^{T})\bm{\Delta}_{B}\|_{F}=\|(\bm{I}-\bm{V}\bm{V}^{T})\bm{\Delta}_{B}+(\bm{V}\bm{V}^{T}-\bm{V}_{1,1}\bm{V}_{1,1}^{T})\bm{\Delta}_{B}\|_{F}$
$\displaystyle=\|(\bm{V}\bm{V}^{T}-\bm{V}_{1,1}\bm{V}_{1,1}^{T})\bm{\Delta}_{B}\|_{F}\leq\frac{1}{2}\|\bm{\Delta}_{B}\|_{F}.$

where the third equality uses (50), and the inequality follows from (52).
This, together with $\bm{V}_{1}\in{\cal O}^{d}$, directly implies (49). ∎

##  Appendix D Auxiliary Results

######  Lemma 5.

Given a matrix $\bm{A}\in\mathbb{R}^{m\times n}$ of rank $r\geq 0$, we have

$\displaystyle\|\bm{A}^{T}\bm{A}\|_{F}\leq\|\bm{A}\|_{F}^{2}\leq\sqrt{r}\|\bm{A}^{T}\bm{A}\|_{F}.$
(53)

###### Proof.

Let $\bm{A}=\bm{U}\bm{\Sigma}\bm{V}^{T}$ be a singular value decomposition of
$\bm{A}$, where

$\displaystyle\bm{\Sigma}=\begin{bmatrix}\tilde{\bm{\Sigma}}&\bm{0}\\\
\bm{0}&\bm{0}\end{bmatrix},\
\tilde{\bm{\Sigma}}=\operatorname{diag}(\sigma_{1},\dots,\sigma_{r}),$

$\bm{U}\in\bm{O}^{m}$, and $\bm{V}\in\bm{O}^{n}$. Then we compute

$\displaystyle\|\bm{A}^{T}\bm{A}\|_{F}^{2}=\|\bm{\Sigma}^{T}\bm{\Sigma}\|_{F}^{2}=\sum_{i=1}^{r}\sigma_{i}^{4},\
\|\bm{A}\|_{F}^{2}=\|\bm{\Sigma}\|_{F}^{2}=\sum_{i=1}^{r}\sigma_{i}^{2},$

which, together with the AM-QM inequality, directly implies (53). ∎