# Towards Resolving the Implicit Bias of
Gradient Descent for Matrix Factorization: Greedy Low-Rank Learning

Zhiyuan Li, Yuping Luo
Department of Computer Science
Princeton University
Princeton, NJ 08544
<EMAIL_ADDRESS>
&Kaifeng Lyu11footnotemark: 1
Institute for Interdisciplinary Information Sciences
Tsinghua University
Beijing, China
<EMAIL_ADDRESS>Alphabet ordering.

###### Abstract

Matrix factorization is a simple and natural test-bed to investigate the
implicit regularization of gradient descent. Gunasekar et al. [2018]
conjectured that Gradient Flow with infinitesimal initialization converges to
the solution that minimizes the nuclear norm, but a series of recent papers
argued that the language of norm minimization is not sufficient to give a full
characterization for the implicit regularization. In this work, we provide
theoretical and empirical evidence that for depth-2 matrix factorization,
gradient flow with infinitesimal initialization is mathematically equivalent
to a simple heuristic rank minimization algorithm, Greedy Low-Rank Learning,
under some reasonable assumptions. This generalizes the rank minimization view
from previous works to a much broader setting and enables us to construct
counter-examples to refute the conjecture from [Gunasekar et al., 2018]. We
also extend the results to the case where depth $\geq 3$, and we show that the
benefit of being deeper is that the above convergence has a much weaker
dependence over initialization magnitude so that this rank minimization is
more likely to take effect for initialization with practical scale.

##  1 Introduction

There are usually far more learnable parameters in deep neural nets than the
number of training data, but still deep learning works well on real-world
tasks. Even with explicit regularization, the model complexity of state-of-
the-art neural nets is so large that they can fit randomly labeled data easily
[Zhang et al., 2017]. Towards explaining the mystery of generalization, we
must understand what kind of implicit regularization does Gradient Descent
(GD) impose during training. Ideally, we are hoping for a nice mathematical
characterization of how GD constrains the set of functions that can be
expressed by a trained neural net.

As a direct analysis for deep neural nets could be quite hard, a line of works
turned to study the implicit regularization on simpler problems to get
inspirations, for example, low-rank matrix factorization, a fundamental
problem in machine learning and information process. Given a set of
observations about an unknown matrix ${\bm{W}}^{*}\in\mathbb{R}^{d\times d}$
of rank $r^{*}\ll d$, one needs to find a low-rank solution ${\bm{W}}$ that is
compatible with the given observations. Examples include matrix sensing,
matrix completion, phase retrieval, robust principal component analysis, just
to name a few (see Chi et al. 2019 for a survey). When ${\bm{W}}^{*}$ is
symmetric and positive semidefinite, one way to solve all these problems is to
parameterize ${\bm{W}}$ as ${\bm{W}}={\bm{U}}{\bm{U}}^{\top}$ for
${\bm{U}}\in\mathbb{R}^{d\times r}$ and optimize
$\mathcal{L}({\bm{U}}):=\frac{1}{2}f({\bm{U}}{\bm{U}}^{\top})$, where
$f(\,\cdot\,)$ is some empirical risk function depending on the observations,
and $r$ is the rank constraint. In theory, if the rank constraint is too
loose, the solutions do not have to be low-rank and we may fail to recover
${\bm{W}}^{*}$. However, even in the case where the rank is unconstrained
(i.e., $r=d$), GD with small initialization can still get good performance in
practice. This empirical observation reveals that the implicit regularization
of GD exists even in this simple matrix factorization problem, but its
mechanism is still on debate. Gunasekar et al. [2017] proved that Gradient
Flow (GD with infinitesimal step size, a.k.a., GF) with infinitesimal
initialization finds the minimum nuclear norm solution in a special case of
matrix sensing, and further conjectured this holds in general.

######  Conjecture 1.1 (Gunasekar et al. 2017, informal).

With sufficiently small initialization, GF converges to the minimum nuclear
norm solution of matrix sensing.

Subsequently, Arora et al. [2019a] challenged this view by arguing that a
simple mathematical norm may not be a sufficient language for characterizing
implicit regularization. One example illustrated in Arora et al. [2019a] is
regarding matrix sensing with a single observation. They showed that GD with
small initialization enhances the growth of large singular values of the
solution and attenuates that of smaller ones. This enhancement/attenuation
effect encourages low-rank, and it is further intensified with depth in deep
matrix factorization (i.e., GD optimizes $f({\bm{U}}_{1}\cdots{\bm{U}}_{L})$
for $L\geq 2$). However, these are not captured by the nuclear norm alone.
Gidel et al. [2019], Gissin et al. [2020] further exploited this idea and
showed in the special case of full-observation matrix sensing that GF learns
solutions with gradually increasing rank. Razin and Cohen [2020] showed in a
simple class of matrix completion problems that GF decreases the rank along
the trajectory while any norm grows towards infinity. More aggressively, they
conjectured that the implicit regularization can be explained by rank
minimization rather than norm minimization.

##### Our Contributions.

In this paper, we move one further step towards resolving the implicit
regularization in the matrix factorization problem. Our theoretical results
show that GD performs rank minimization via a greedy process in a broader
setting. Specifically, we provide theoretical evidence that GF with
infinitesimal initialization is in general mathematically equivalent to
another algorithm called Greedy Low-Rank Learning (GLRL). At a high level,
GLRL is a greedy algorithm that performs rank-constrained optimization and
relaxes the rank constraint by $1$ whenever it fails to reach a global
minimizer of $f(\,\cdot\,)$ with the current rank constraint. As a by-product,
we refute 1.1 by demonstrating an counterexample (Example 5.9).

We also extend our results to deep matrix factorization Section 6, where we
prove that the trajectory of GF with infinitesimal identity initialization
converges to a deep version of GLRL, at least in the early stage of the
optimization. We also use this result to confirm the intuition achieved on toy
models [Gissin et al., 2020], that benefits of depth in matrix factorization
is to encourage rank minimization even for initialization with a relatively
larger scale, and thus it is more likely to happen in practice. This shows
that describing the implicit regularization using GLRL is more expressive than
using the language of norm minimization. We validate all our results with
experiments in Appendix C.

##  2 Related Works

##### Norm Minimization.

The view of norm minimization, or the closely related view of margin
maximization, has been explored in different settings. Besides the nuclear
norm minimization for matrix factorization [Gunasekar et al., 2017] discussed
in the introduction, previous works have also studied the norm
minimization/margin maximization for linear regression [Wilson et al., 2017,
Soudry et al., 2018a, b, Nacson et al., 2019b, c, Ji and Telgarsky, 2019b],
deep linear neural nets [Ji and Telgarsky, 2019a, Gunasekar et al., 2018],
homogeneous neural nets [Nacson et al., 2019a, Lyu and Li, 2020], ultra-wide
neural nets [Jacot et al., 2018, Arora et al., 2019b, Chizat and Bach, 2020].

##### Small Initialization and Rank Minimization.

The initialization scale can greatly influence the implicit regularization. A
sufficiently large initialization can make the training dynamics fall into the
lazy training regime defined by Chizat et al. [2019] and diminish test
accuracy. Using small initialization is particularly important to bias
gradient descent to low-rank solutions for matrix factorization, as
empirically observed by Gunasekar et al. [2017]. Arora et al. [2019a], Gidel
et al. [2019], Gissin et al. [2020], Razin and Cohen [2020] studied how
gradient flow with small initialization encourages low-rank in simple
settings, as discussed in the introduction. Li et al. [2018] proved recovery
guarantees for gradient flow solving matrix sensing under Restricted Isometry
Property (RIP), but the proof cannot be generalized easily to the case without
RIP. Belabbas [2020] made attempts to prove that gradient flow is
approximately rank-1 in the very early phase of training, but it does not
exclude the possibility that the approximation error explodes later and
gradient flow is not converging to low-rank solutions. Compared to these
works, the current paper studies how GF encourages low-rank in a much broader
setting.

##  3 Background

##### Notations.

For two matrices ${\bm{A}},{\bm{B}}$, we define
$\left<{\bm{A}},{\bm{B}}\right>:=\operatorname{Tr}({\bm{A}}{\bm{B}}^{\top})$
as their inner product. We use
$\left\|{\bm{A}}\right\|_{\mathrm{F}},\left\|{\bm{A}}\right\|_{\ast}$ and
$\left\|{\bm{A}}\right\|_{2}$ to denote the Frobenius norm, nuclear norm and
the largest singular value of ${\bm{A}}$ respectively. For a matrix
${\bm{A}}\in\mathbb{R}^{d\times d}$, we use
$\lambda_{1}({\bm{A}}),\dots,\lambda_{d}({\bm{A}})$ to denote the eigenvalues
of ${\bm{A}}$ in decreasing order (if they are all reals). We define
${\mathbb{S}}_{d}$ as the set of symmetric $d\times d$ matrices and
${\mathbb{S}}^{+}_{d}\subseteq{\mathbb{S}}_{d}$ as the set of positive
semidefinite (PSD) matrices. We write ${\bm{A}}\succeq{\bm{B}}$ or
${\bm{B}}\preceq{\bm{A}}$ if ${\bm{A}}-{\bm{B}}$ is PSD. We use
${\mathbb{S}}^{+}_{d,r}$, ${\mathbb{S}}^{+}_{d,\leq r}$ to denote the set of
$d\times d$ PSD matrices with rank $=r,\leq r$ respectively.

##### Matrix Factorization.

Matrix factorization problem asks one to optimize
$\mathcal{L}({\bm{U}},{\bm{V}}):=\frac{1}{2}f({\bm{U}}{\bm{V}}^{\top})$ among
${\bm{U}},{\bm{V}}\in\mathbb{R}^{d\times r}$, where $f:\mathbb{R}^{d\times
d}\to\mathbb{R}$ is a convex function. A notable example is matrix sensing.
There is an unknown rank-$r^{*}$ matrix ${\bm{W}}^{*}\in\mathbb{R}^{d\times
d}$ with $r^{*}\ll d$. Given $m$ measurements
${\bm{X}}_{1},\dots,{\bm{X}}_{m}\in\mathbb{R}^{d\times d}$, one can observe
$y_{i}:=\left<{\bm{X}}_{i},{\bm{W}}^{*}\right>$ through each measurement. The
goal of matrix sensing is to reconstruct ${\bm{W}}^{*}$ via minimizing
$f({\bm{W}}):=\frac{1}{2}\sum_{i=1}^{m}\left(\left<{\bm{W}},{\bm{X}}_{i}\right>-y_{i}\right)^{2}$.
Matrix completion is a notable special case of matrix sensing in which every
measurement has the form
${\bm{X}}_{i}={\bm{e}}_{p_{i}}{\bm{e}}_{q_{i}}^{\top}$, where
$\\{{\bm{e}}_{1},\cdots,{\bm{e}}_{d}\\}$ stands for the standard basis (i.e.,
exactly one entry is observed through each measurement).

For technical simplicity, in this paper we focus on the symmetric case as in
previous works [Gunasekar et al., 2017]. Given a $\mathcal{C}^{3}$-smooth
convex function $f:\mathbb{R}^{d\times d}\to\mathbb{R}$, we aim to find a low-
rank solution for the convex optimization problem (P):

$\min f({\bm{W}})\quad\text{s.t.}\quad{\bm{W}}\succeq{\bm{0}}$ (P)

For this, we parameterize ${\bm{W}}$ as ${\bm{W}}={\bm{U}}{\bm{U}}^{\top}$ for
${\bm{U}}\in\mathbb{R}^{d\times r}$ and optimize
$\mathcal{L}({\bm{U}}):=\frac{1}{2}f({\bm{U}}{\bm{U}}^{\top})$. We assume WLOG
throughout this paper that $f({\bm{W}})=f({\bm{W}}^{\top})$; otherwise, we can
set
$f^{\prime}({\bm{W}})=\frac{1}{2}\left(f({\bm{W}})+f({\bm{W}}^{\top})\right)$
so that $f^{\prime}({\bm{W}})=f^{\prime}({\bm{W}}^{\top})$ while
$\mathcal{L}({\bm{U}})=\frac{1}{2}f^{\prime}({\bm{U}}{\bm{U}}^{\top})$ is
unaffected. This assumption makes $\nabla f({\bm{W}})$ symmetric for every
symmetric ${\bm{W}}$.

Note that matrix factorization in the general case can be reduced to this
symmetric case: let ${\bm{U}}^{\prime}=\left[\begin{smallmatrix}{\bm{U}}\\\
{\bm{V}}\end{smallmatrix}\right]\in\mathbb{R}^{2d\times r}$,
$f^{\prime}\left(\left[\begin{smallmatrix}{\bm{A}}&{\bm{B}}\\\
{\bm{C}}&{\bm{D}}\end{smallmatrix}\right]\right)=\frac{1}{2}f({\bm{B}})+\frac{1}{2}f({\bm{C}})$,
then
$f({\bm{U}}{\bm{V}}^{\top})=f^{\prime}({\bm{U}}^{\prime}{\bm{U}}^{\prime\top})$.
So focusing on the symmetric case does not lose generality.

##### Gradient Flow.

In this paper, we analyze Gradient Flow (GF) on symmetric matrix
factorization, which is defined by the following ODE for
${\bm{U}}(t)\in\mathbb{R}^{d\times r}$:

$\frac{\textup{{d}}{\bm{U}}}{\textup{{d}}t}=-\nabla\mathcal{L}({\bm{U}})=-\nabla
f({\bm{U}}{\bm{U}}^{\top}){\bm{U}}.$ (1)

Let ${\bm{W}}(t)={\bm{U}}(t){\bm{U}}(t)^{\top}\in\mathbb{R}^{d\times d}$. Then
the following end-to-end dynamics holds for ${\bm{W}}(t)$:

$\frac{\textup{{d}}{\bm{W}}}{\textup{{d}}t}=-{\bm{W}}\nabla f({\bm{W}})-\nabla
f({\bm{W}}){\bm{W}}=:{\bm{g}}({\bm{W}}).$ (2)

We use $\phi({\bm{W}}_{0},t)$ to denote the matrix ${\bm{W}}(t)$ in (2) when
${\bm{W}}(0)={\bm{W}}_{0}\succeq{\bm{0}}$. Throughout this paper, we assume
$\phi({\bm{W}}_{0},t)$ exists for all
$t\in{\mathbb{R}},{\bm{W}}_{0}\succeq{\bm{0}}$. It is easy to prove that
${\bm{U}}$ is a stationary point of $\mathcal{L}(\,\cdot\,)$ (i.e.,
$\nabla\mathcal{L}({\bm{U}})={\bm{0}}$) iff ${\bm{W}}={\bm{U}}{\bm{U}}^{\top}$
is a critical point of (2) (i.e., ${\bm{g}}({\bm{W}})={\bm{0}}$); see Lemma
A.1 for a proof. If ${\bm{W}}$ is a minimizer of $f(\,\cdot\,)$ in
${\mathbb{S}}^{+}_{d}$ (i.e., ${\bm{W}}$ is a minimizer of (P)), then
${\bm{W}}$ is a critical point of (2), but the reverse may not be true, e.g.,
${\bm{g}}({\bm{0}})={\bm{0}}$, but ${\bm{0}}$ is not necessarily a minimizer.

In this paper, we particularly focus on the overparameterized case, where
$r=d$, to understand the implicit regularization of GF when there is no rank
constraint for the matrix ${\bm{W}}$.

##  4 Warmup Examples

Before introducing our main results, we illustrate how GD performs greedy
learning using two warmup examples.

##### Linearization Around the Origin.

In general, for a loss function
$\mathcal{L}({\bm{U}})=\frac{1}{2}f({\bm{U}}{\bm{U}}^{\top})$, we can always
apply Taylor expansion $f({\bm{W}})\approx f({\bm{0}})+\left<{\bm{W}},\nabla
f({\bm{0}})\right>$ around the origin to approximate it with a linear
function. This motivates us to study the linear case:
$f({\bm{W}}):=f_{0}-\left<{\bm{W}},{\bm{Q}}\right>$ for some symmetric matrix
${\bm{Q}}$. In this case, the matrix ${\bm{U}}$ follows the ODE,
$\frac{\textup{{d}}{\bm{U}}}{\textup{{d}}t}={\bm{Q}}{\bm{U}}$, which can be
understood as a continuous version of the classical power iteration method for
solving the top eigenvector. Let
${\bm{Q}}:=\sum_{i=1}^{d}\mu_{i}{\bm{v}}_{i}{\bm{v}}_{i}^{\top}$ be the
eigendecomposition of ${\bm{Q}}$, where
$\mu_{1}\geq\mu_{2}\geq\cdots\geq\mu_{d}$ and
${\bm{v}}_{1},\dots,{\bm{v}}_{d}$ are orthogonal to each other. Then we can
write the solution as:

${\bm{U}}(t)=e^{t{\bm{Q}}}{\bm{U}}(0)=\left(\sum\nolimits_{i=1}^{d}e^{\mu_{i}t}{\bm{v}}_{i}{\bm{v}}_{i}^{\top}\right){\bm{U}}(0).$
(3)

When $\mu_{1}>\mu_{2}$, the ratio between $e^{\mu_{1}t}$ and $e^{\mu_{i}t}$
for $i\neq 1$ increases exponentially fast. As $t\to+\infty$, ${\bm{U}}(t)$
and ${\bm{W}}(t)$ become approximately rank-1 as long as
${\bm{v}}_{i}^{\top}{\bm{U}}(0)\neq{\bm{0}}$, i.e.,

$\lim_{t\to\infty}e^{-\mu_{1}t}{\bm{U}}(t)={\bm{v}}_{1}{\bm{v}}_{1}^{\top}{\bm{U}}(0),\qquad\lim_{t\to\infty}e^{-2\mu_{1}t}{\bm{W}}(t)=({\bm{v}}_{1}^{\top}{\bm{W}}(0){\bm{v}}_{1}){\bm{v}}_{1}{\bm{v}}_{1}^{\top}.$
(4)

The analysis for the simple linear case reveals that GD encourages low-rank
through a process similar to power iteration. However, $f({\bm{W}})$ is non-
linear in general, and the linear approximation is close to $f({\bm{W}})$ only
if ${\bm{W}}$ is very small. With sufficiently small initialization, we can
imagine that GD still resembles the above power iteration in the early phase
of the optimization. But what if ${\bm{W}}(t)$ grows to be so large that the
linear approximation is far from the actual $f({\bm{W}})$?

##### Full-observation Matrix Sensing.

To understand the dynamics of GD when the linearization fails, we now consider
a well-studied special case [Gissin et al., 2020]:
$\mathcal{L}({\bm{U}})=\frac{1}{2}f({\bm{U}}{\bm{U}}^{\top}),f({\bm{W}})=\frac{1}{2}\|{\bm{W}}-{\bm{W}}^{*}\|_{\mathrm{F}}^{2}$
for some unknown PSD matrix ${\bm{W}}^{*}$. GF in this case can be written as:

$\frac{\textup{{d}}{\bm{U}}}{\textup{{d}}t}=({\bm{W}}^{*}-{\bm{U}}{\bm{U}}^{\top}){\bm{U}},\qquad\frac{\textup{{d}}{\bm{W}}}{\textup{{d}}t}=({\bm{W}}^{*}-{\bm{W}}){\bm{W}}+{\bm{W}}({\bm{W}}^{*}-{\bm{W}}).$
(5)

Let ${\bm{W}}^{*}:=\sum_{i=1}^{d}\mu_{i}{\bm{v}}_{i}{\bm{v}}_{i}^{\top}$ be
the eigendecomposition of ${\bm{W}}^{*}$. Our previous analysis shows that the
dynamics is approximately
$\frac{\textup{{d}}{\bm{U}}}{\textup{{d}}t}={\bm{W}}^{*}{\bm{U}}$ in the early
phase and thus encourages low-rank.

To get a sense for the later phases, we simplify the setting by specifying
${\bm{U}}(0)=\sqrt{\alpha}{\bm{I}}$ for a small number $\alpha$. We can write
${\bm{W}}(0)$ and ${\bm{W}}^{*}$ as diagonal matrices
${\bm{W}}(0)=\mathrm{diag}(\alpha,\alpha,\cdots,\alpha),{\bm{W}}^{*}=\mathrm{diag}(\mu_{1},\mu_{2},\cdots,\mu_{d})$
with respect to the basis ${\bm{v}}_{1},\dots,{\bm{v}}_{d}$. It is easy to see
that ${\bm{W}}(t)$ is always a diagonal matrix, since the time derivatives of
non-diagonal coordinates stay $0$ during training. Let
${\bm{W}}(t)=\mathrm{diag}(\sigma_{1}(t),\sigma_{2}(t),\cdots,\sigma_{d}(t))$,
then $\sigma_{i}(t)$ satisfies the dynamical equation
$\frac{\textup{{d}}}{\textup{{d}}t}\sigma_{i}(t)=2\sigma_{i}(t)(\mu_{i}-\sigma_{i}(t))$,
and thus
$\sigma_{i}(t)=\frac{\alpha\mu_{i}}{\alpha+(\mu_{i}-\alpha)e^{-2\mu_{i}t}}$.
This shows that every $\sigma_{i}(t)$ increases from $\alpha$ to $\mu_{i}$
over time. As $\alpha\to 0$, every $\sigma_{i}(t)$ has a sharp transition from
near $0$ to near $\mu_{i}$ at time roughly
$(\frac{1}{2\mu_{i}}+o(1))\log\frac{1}{\alpha}$, which can be seen from the
following limit:

$\lim_{\alpha\to
0}\sigma_{i}\left((\tfrac{1}{2\mu_{i}}+c)\log(1/\alpha)\right)=\lim_{\alpha\to
0}\frac{\alpha\mu_{i}}{\alpha+(\mu_{i}-\alpha)\alpha^{1+2c\mu_{i}}}=\begin{cases}0&\quad
c\in(-\frac{1}{2\mu_{i}},0),\\\ \mu_{i}&\quad c\in(0,+\infty).\end{cases}$

This means for every $q\in(\frac{1}{2\mu_{i}},\frac{1}{2\mu_{i+1}})$ for
$i=1,\dots,d-1$ (or $q\in(\frac{1}{2\mu_{i}},+\infty)$ for $i=d$),
$\lim_{\alpha\to
0}{\bm{W}}(q\log(1/\alpha))=\mathrm{diag}(\mu_{1},\mu_{2},\dots,\mu_{i},0,0,\cdots,0)$.
Therefore, when the initialization is sufficiently small, GF learns each
component of ${\bm{W}}^{*}$ one by one, according to the relative order of
eigenvalues. At a high level, this shows a greedy nature of GD: GD starts
learning with simple models; whenever it underfits, it increases the model
complexity (which is rank in our case). This is also called sequential
learning or incremental learning in the literature [Gidel et al., 2019, Gissin
et al., 2020].

However, it is unclear how and why this sequential learning/incremental
learning can occur in general. Through the first warmup example, we may
understand why GD learns a rank-1 matrix in the early phase, but does GD
always learn solutions with rank $2,3,4,\dots$ sequentially? If true, what is
the mechanism behind this? The current paper answers the questions by
providing both theoretical and empirical evidence that the greedy learning
behavior does occur in general with a similar reason as for the first warmup
example.

##  5 Greedy Low-Rank Learning (GLRL)

In this section, we present a trajectory-based analysis for the implicit bias
of GF on matrix factorization. Our main result is that GF with infinitesimal
initialization is generically the same as that of a simple greedy algorithm,
_Greedy Low-Rank Learning_ (GLRL, Algorithm 1).

The GLRL algorithm consists of several phases, numbered from $1$. In phase
$r$, GLRL increases the rank constraint to $r$ and optimizes
$\mathcal{L}({\bm{U}}_{r}):=\frac{1}{2}f({\bm{U}}_{r}{\bm{U}}_{r}^{\top})$
among ${\bm{U}}_{r}\in\mathbb{R}^{d\times r}$ via GD until it reaches a
stationary point ${\bm{U}}_{r}(\infty)$, i.e.,
$\nabla\mathcal{L}({\bm{U}}_{r}(\infty))={\bm{0}}$. At convergence,
${\bm{W}}_{r}:={\bm{U}}_{r}(\infty){\bm{U}}^{\top}_{r}(\infty)$ is a critical
point of (2), and we call it the $r$-th critical point of GLRL. If
${\bm{W}}_{r}$ is further a minimizer of $f(\,\cdot\,)$ in
${\mathbb{S}}^{+}_{d}$, or equivalently, $\lambda_{1}(-\nabla
f({\bm{W}}_{r}))\leq 0$ (see Lemma A.2), then GLRL returns ${\bm{W}}_{r}$;
otherwise GLRL enters phase $r+1$.

parameter : step size $\eta>0$; small $\epsilon>0$

$r\leftarrow 0,{\bm{W}}_{0}\leftarrow{\bm{0}}\in\mathbb{R}^{d\times d}$, and
${\bm{U}}_{0}(\infty)\in\mathbb{R}^{d\times 0}$ is an empty matrix

while  _$\lambda_{1}(-\nabla f({\bm{W}}_{r})) >0$_ do

$r\leftarrow r+1$

${\bm{u}}_{r}\leftarrow$ unit top eigenvector of $-\nabla f({\bm{W}}_{r-1})$

${\bm{U}}_{r}(0)\leftarrow[{\bm{U}}_{r-1}(\infty)\ \
\sqrt{\epsilon}{\bm{u}}_{r}]\in\mathbb{R}^{d\times r}$

for  _$t=0,1,\ldots$_ do

${\bm{U}}_{r}(t+1)\leftarrow{\bm{U}}_{r}(t)-\eta\nabla\mathcal{L}({\bm{U}}_{r}(t))$

${\bm{W}}_{r}\leftarrow{\bm{U}}_{r}(\infty){\bm{U}}^{\top}_{r}(\infty)$ 111In
practice, we approximate the infinite time limit by running sufficiently many
steps.

return ${\bm{W}}_{r}$

Algorithm 1 Greedy Low-Rank Learning (GLRL)

To set the initial point of GD in phase $r$, GLRL appends a small column
vector ${\bm{\delta}}_{r}\in\mathbb{R}^{d}$ to the resulting stationary point
${\bm{U}}_{r-1}(\infty)$ from the last phase, i.e.,
${\bm{U}}_{r}(0)\leftarrow\left[{\bm{U}}_{r-1}(\infty)~{}~{}{\bm{\delta}}_{r}\right]\in\mathbb{R}^{d\times
r}$ (in the case of $r=1$,
${\bm{U}}_{1}(0)\leftarrow\left[{\bm{\delta}}_{1}\right]\in\mathbb{R}^{d\times
1}$). In this way,
${\bm{U}}_{r}(0){\bm{U}}^{\top}_{r}(0)={\bm{W}}_{r-1}+{\bm{\delta}}_{r}{\bm{\delta}}_{r}^{\top}$
is perturbed away from the $(r-1)$-th critical point. In GLRL, we set
${\bm{\delta}}_{r}=\sqrt{\epsilon}{\bm{u}}_{r}$, where ${\bm{u}}_{r}$ is the
top eigenvector of $-\nabla f({\bm{W}}_{r})$ with unit norm
$\|{\bm{u}}_{r}\|_{2}=1$, and $\epsilon>0$ is a parameter controlling the
magnitude of perturbation (preferably very small). Note that it is guaranteed
that $\lambda_{1}(-\nabla f({\bm{W}}_{r-1}))>0$; otherwise ${\bm{W}}_{r-1}$ is
a minimizer of the convex function $f(\,\cdot\,)$ in ${\mathbb{S}}^{+}_{d}$
and GLRL exits before phase $r$. Expanding $f(\,\cdot\,)$ around
${\bm{W}}_{r-1}$ shows that the loss is decreasing in this choice of
${\bm{\delta}}_{r}$.

$\displaystyle\mathcal{L}({\bm{U}}_{r}(0))=\frac{1}{2}f({\bm{W}}_{r-1}+{\bm{\delta}}_{r}{\bm{\delta}}_{r}^{\top})$
$\displaystyle=\mathcal{L}({\bm{U}}_{r-1}(\infty))+\frac{1}{2}{\bm{\delta}}_{r}^{\top}\nabla
f({\bm{W}}_{r-1}){\bm{\delta}}_{r}+O(\|{\bm{\delta}}_{r}\|_{2}^{4})$
$\displaystyle=\mathcal{L}({\bm{U}}_{r-1}(\infty))-\frac{\epsilon}{2}\lambda_{1}(-\nabla
f({\bm{W}}_{r-1}))+O(\epsilon^{2}).$

##### Trajectory of GLRL.

We define the (limiting) trajectory of GLRL by taking the learning rate
$\eta\to 0$. The goal is to show that the trajectory of GLRL is close to that
of GF with infinitesimal initialization. Recall that $\phi({\bm{W}}_{0},t)$
stands for the solution ${\bm{W}}(t)$ in (2) when ${\bm{W}}(0)={\bm{W}}_{0}$.

######  Definition 5.1 (Trajectory of GLRL).

Let $\overline{{}{\bm{W}}}_{0,\epsilon}:={\bm{0}}$ be the 0th critical point
of GLRL. For every $r\geq 1$, if the $(r-1)$-th critical point
$\overline{{}{\bm{W}}}_{r-1,\epsilon}$ exists and is not a minimizer of
$f(\,\cdot\,)$ in ${\mathbb{S}}^{+}_{d}$, we define
${\bm{W}}^{\mathrm{G}}_{r,\epsilon}(t):=\phi(\overline{{}{\bm{W}}}_{r-1,\epsilon}+\epsilon{\bm{u}}_{r,\epsilon}{\bm{u}}_{r,\epsilon}^{\top},t)$,
where ${\bm{u}}_{r,\epsilon}$ is a top eigenvector of $\nabla
f(\overline{{}{\bm{W}}}_{r-1,\epsilon})$ with unit norm,
$\|{\bm{u}}_{r,\epsilon}\|_{2}=1$. We define
$\overline{{}{\bm{W}}}_{r,\epsilon}:=\lim_{t\to+\infty}{\bm{W}}^{\mathrm{G}}_{r,\epsilon}(t)$
to be the $r$-th critical point of GLRL if the limit exists.

Throughout this paper, we always focus on the case where the top eigenvalue of
every $\nabla f(\overline{{}{\bm{W}}}_{r-1,\epsilon})$ is unique. In this
case, the trajectory of GLRL is unique for every $\epsilon>0$, since the
normalized top eigenvectors can only be $\pm{\bm{u}}_{r,\epsilon}$, and both
of them lead to the same ${\bm{W}}^{\mathrm{G}}_{r,\epsilon}(t)$.

##### Comparison to existing greedy algorithms for rank-constrained
optimization.

The most related one to GLRL (Algorithm 1) is probably _Rank-1 Matrix Pursuit_
(R1MP) proposed by Wang et al. [2014] for matrix completion, which was later
generalized to general convex loss in [Yao and Kwok, 2016]. R1MP maintains a
set of rank-1 matrices as the basis, and in phase $r$, R1MP adds the same
${\bm{u}}_{r}{\bm{u}}_{r}^{\top}$ as defined in Algorithm 1 into its basis and
solve
$\min_{{\bm{\alpha}}}f(\sum_{i=1}^{r}\alpha_{i}{\bm{u}}_{i}{\bm{u}}_{i}^{\top})$
for rank-$r$ estimation. The main difference between R1MP and GLRL is that the
optimization in each phase of R1MP is performed on the coefficients
${\bm{\alpha}}$, while the entire ${\bm{U}}_{r}$ evolves with GD in each phase
of GLRL. In LABEL:fig:group_fight, we provide empirical evidence that GLRL
generalizes better than R1MP when ground truth is low-rank, although GLRL may
have a higher computational cost depending on $\eta,\epsilon$.

Similar to R1MP, Greedy Efficient Component Optimization (GECO, Shalev-Shwartz
and Singer 2010) also chooses the $r$-th component of its basis as the top
eigenvector of $-\nabla f({\bm{W}}_{r})$, while it solves
$\min_{{\bm{\beta}}}f(\sum_{1\leq i,j\leq
r}\beta_{ij}{\bm{u}}_{i}{\bm{u}}_{j}^{\top})$ for the rank-$r$ estimation.
Khanna et al. [2017] provided convergence guarantee for GECO assuming strong
convexity. Haeffele and Vidal [2019] proposed a local-descent meta algorithm,
of which GLRL can be viewed as a specific realization.

###  5.1 The Limiting Trajectory: A General Theorem for Dynamical System

To prove the equivalence between GF and GLRL, we first introduce our high-
level idea by analyzing the behavior of a more general dynamical system around
its critical point, say ${\bm{0}}$. A specific example is (2) if we set
${\bm{\theta}}$ to be the vectorization of ${\bm{W}}$.

$\frac{\textup{{d}}{\bm{\theta}}}{\textup{{d}}t}={\bm{g}}({\bm{\theta}}),\quad\text{where}\quad{\bm{g}}({\bm{0}})={\bm{0}}.$
(6)

We use $\phi({\bm{\theta}}_{0},t)$ to denote the value of ${\bm{\theta}}(t)$
in the case of ${\bm{\theta}}(0)={\bm{\theta}}_{0}$. We assume that
${\bm{g}}({\bm{\theta}})$ is $\mathcal{C}^{2}$-smooth with
${\bm{J}}({\bm{\theta}})$ being the Jacobian matrix and
$\phi({\bm{\theta}}_{0},t)$ exists for all ${\bm{\theta}}_{0}$ and $t$. For
ease of presentation, in the main text we assume ${\bm{J}}({\bm{0}})$ is
diagonalizable over $\mathbb{R}$ and defer the same result for the general
case into Section E.3. Let
${\bm{J}}({\bm{0}})=\tilde{{\bm{V}}}\tilde{{\bm{D}}}\tilde{{\bm{V}}}^{-1}$ be
the eigendecomposition, where $\tilde{{\bm{V}}}$ is an invertible matrix and
$\tilde{{\bm{D}}}=\mathrm{diag}(\tilde{\mu}_{1},\dots,\tilde{\mu}_{d})$ is the
diagonal matrix consisting of the eigenvalues
$\tilde{\mu}_{1}\geq\tilde{\mu}_{2}\geq\cdots\geq\tilde{\mu}_{d}$. Let
$\tilde{{\bm{V}}}=(\tilde{{\bm{v}}}_{1},\dots,\tilde{{\bm{v}}}_{d})$ and
$\tilde{{\bm{V}}}^{-1}=(\tilde{{\bm{u}}}_{1},\dots,\tilde{{\bm{u}}}_{d})^{\top}$,
then $\tilde{{\bm{u}}}_{i}$, $\tilde{{\bm{v}}}_{i}$ are left and right
eigenvectors associated with $\tilde{\mu}_{i}$ and
$\tilde{{\bm{u}}}_{i}^{\top}\tilde{{\bm{v}}}_{j}=\delta_{ij}$. We can rewrite
the eigendecomposition as
${\bm{J}}({\bm{0}})=\sum_{i=1}^{d}\tilde{\mu}_{i}\tilde{{\bm{v}}}_{i}\tilde{{\bm{u}}}_{i}^{\top}$.

We also assume the top eigenvalue $\tilde{\mu}_{1}$ is positive and unique.
Note $\tilde{\mu}_{1}>0$ means the critical point ${\bm{\theta}}={\bm{0}}$ is
unstable, and in matrix factorization it means ${\bm{0}}$ is a strict saddle
point of $\mathcal{L}(\,\cdot\,)$.

The key observation is that if the initialization is infinitesimal, the
trajectory is almost uniquely determined. To be more precise, we need the
following definition:

######  Definition 5.2.

For any ${\bm{\theta}}_{0}\in\mathbb{R}^{d}$ and ${\bm{u}}\in\mathbb{R}^{d}$,
we say that $\\{{\bm{\theta}}_{\alpha}\\}_{\alpha\in(0,1)}$ converges to
${\bm{\theta}}_{0}$ with positive alignment with ${\bm{u}}$ if
$\lim\limits_{\alpha\to 0}{\bm{\theta}}_{\alpha}={\bm{\theta}}_{0}$ and
$\liminf\limits_{\alpha\to
0}\left<\frac{{\bm{\theta}}_{\alpha}-{\bm{\theta}}_{0}}{\|{\bm{\theta}}_{\alpha}-{\bm{\theta}}_{0}\|_{2}},{\bm{u}}\right>>0$.

A special case is that the direction of
${\bm{\theta}}_{\alpha}-{\bm{\theta}}_{0}$ converges, i.e.,
$\bar{{\bm{\theta}}}:=\lim_{\alpha\to
0}\frac{{\bm{\theta}}_{\alpha}-{\bm{\theta}}_{0}}{\|{\bm{\theta}}_{\alpha}-{\bm{\theta}}_{0}\|_{2}}$
exists. In this case, $\\{{\bm{\theta}}_{\alpha}\\}$ has positive alignment
with either ${\bm{u}}$ or $-{\bm{u}}$ except for a zero-measure subset of
$\bar{{\bm{\theta}}}$. This means any convergent sequence generically falls
into either of these two categories.

The following theorem shows that if the initial point ${\bm{\theta}}_{\alpha}$
converges to ${\bm{0}}$ with positive alignment with $\tilde{{\bm{u}}}_{1}$ as
$\alpha\to 0$, the trajectory starting with ${\bm{\theta}}_{\alpha}$ converges
to a unique trajectory
${\bm{z}}(t):=\phi(\alpha\tilde{{\bm{v}}}_{1},t+\frac{1}{\tilde{\mu}_{1}}\log\frac{1}{\alpha})$.
By symmetry, there is another unique trajectory for sequences
$\\{{\bm{\theta}}_{\alpha}\\}$ with positive alignment to
$-\tilde{{\bm{u}}}_{1}$, which is
${\bm{z}}^{\prime}(t):=\phi(-\alpha\tilde{{\bm{v}}}_{1},t+\frac{1}{\tilde{\mu}_{1}}\log\frac{1}{\alpha})$.
This is somewhat surprising: different initial points should lead to very
different trajectories, but our analysis shows that generically there are only
two limiting trajectories for infinitesimal initialization. We will soon see
how this theorem helps in our analysis for matrix factorization in Sections
5.2 and 5.3.

######  Theorem 5.3.

Let
${\bm{z}}_{\alpha}(t):=\phi(\alpha\tilde{{\bm{v}}}_{1},t+\frac{1}{\tilde{\mu}_{1}}\log\frac{1}{\alpha})$
for every $\alpha>0$, then ${\bm{z}}(t):=\lim_{\alpha\to
0}{\bm{z}}_{\alpha}(t)$ exists and is also a solution of
$\eqref{eq:dynamical}$, i.e., ${\bm{z}}(t)=\phi({\bm{z}}(0),t)$. If
${\bm{\delta}}_{\alpha}$ converges to ${\bm{0}}$ with positive alignment with
$\tilde{{\bm{u}}}_{1}$ as $\alpha\to 0$, then $\forall t\in\mathbb{R}$, there
is a constant $C>0$ such that

$\left\|\phi\left({\bm{\delta}}_{\alpha},t+\tfrac{1}{\tilde{\mu}_{1}}\log\tfrac{1}{\left<{\bm{\delta}}_{\alpha},\tilde{{\bm{u}}}_{1}\right>}\right)-{\bm{z}}(t)\right\|_{2}\leq
C\cdot\|{\bm{\delta}}_{\alpha}\|_{2}^{\frac{\tilde{\gamma}}{\tilde{\mu}_{1}+\tilde{\gamma}}},$
(7)

for every sufficiently small $\alpha$, where
$\tilde{\gamma}:=\tilde{\mu}_{1}-\tilde{\mu}_{2}>0$ is the eigenvalue gap.

###### Proof sketch.

The main idea is to linearize the dynamics near origin as we have done for the
first warmup example. For sufficiently small ${\bm{\theta}}$, by Taylor
expansion of ${\bm{g}}({\bm{\theta}})$, the dynamics is approximately
$\frac{\textup{{d}}{\bm{\theta}}}{\textup{{d}}t}\approx{\bm{J}}({\bm{0}}){\bm{\theta}}$,
which can be understood as a continuous version of power iteration. If the
linear approximation is exact, then
${\bm{\theta}}(t)=e^{t{\bm{J}}(0)}{\bm{\theta}}(0)$. For large enough $t_{0}$,
$e^{t_{0}{\bm{J}}(0)}=\sum_{i=1}^{d}e^{\tilde{\mu}_{i}t_{0}}\tilde{{\bm{v}}}_{i}\tilde{{\bm{u}}}_{i}^{\top}=e^{\tilde{\mu}_{1}t_{0}}\tilde{{\bm{v}}}_{1}\tilde{{\bm{u}}}_{1}^{\top}+O(e^{\tilde{\mu}_{2}t_{0}})$.
Therefore, as long as the initial point ${\bm{\theta}}(0)$ has a positive
inner product with $\tilde{{\bm{u}}}_{1}$, ${\bm{\theta}}(t_{0})$ should be
very close to $\epsilon\tilde{{\bm{v}}}_{1}$ for some $\epsilon>0$, and the
rest of the trajectory after $t_{0}$ should be close to the trajectory
starting from $\epsilon\tilde{{\bm{v}}}_{1}$. However, here is a tradeoff: we
should choose $t_{0}$ to be large enough so that the power iteration takes
effect; but if $t_{0}$ is so large that the norm of ${\bm{\theta}}(t_{0})$
reaches a constant scale, then the linearization fails unavoidably.
Nevertheless, if the initialization scale is sufficiently small, we show via a
careful error analysis that there is always a suitable choice of $t_{0}$ such
that ${\bm{\theta}}(t_{0})$ is well approximated by
$\epsilon\tilde{{\bm{v}}}_{1}$ and the difference between
${\bm{\theta}}(t_{0}+t)$ and $\phi(\epsilon\tilde{{\bm{v}}}_{1},t)$ is bounded
as well. We defer the details to Appendix E. ∎

###  5.2 Equivalence Between GD and GLRL: Rank-One Case

Now we establish the equivalence between GF and GLRL in the first phase. The
main idea is to apply Theorem 5.3 on (2). For this, we need the following
lemma on the eigenvalues and eigenvectors.

######  Lemma 5.4.

Let ${\bm{g}}({\bm{W}}):=-{\bm{W}}\nabla f({\bm{W}})-\nabla
f({\bm{W}}){\bm{W}}$ and ${\bm{J}}({\bm{W}})$ be its Jacobian. Then
${\bm{J}}({\bm{0}})$ is symmetric and thus diagonalizable. Let $-\nabla
f({\bm{0}})=\sum_{i=1}^{d}\mu_{i}{\bm{u}}_{1[i]}{\bm{u}}_{1[i]}^{\top}$ be the
eigendecomposition of the symmetric matrix $-\nabla f({\bm{0}})$, where
$\mu_{1}\geq\mu_{2}\geq\cdots\geq\mu_{d}$. Then ${\bm{J}}({\bm{0}})$ has the
form:

${\bm{J}}({\bm{0}})[{\bm{\Delta}}]=\sum_{i=1}^{d}\sum_{j=1}^{d}(\mu_{i}+\mu_{j})\left<{\bm{\Delta}},{\bm{u}}_{1[i]}{\bm{u}}_{1[j]}^{\top}\right>{\bm{u}}_{1[i]}{\bm{u}}_{1[j]}^{\top},$
(8)

where ${\bm{J}}({\bm{0}})[{\bm{\Delta}}]$ stands for the resulting matrix
produced by left-multiplying ${\bm{J}}({\bm{0}})$ to the vectorization of
${\bm{\Delta}}$. For every pair of $1\leq i\leq j\leq d$, $\mu_{i}+\mu_{j}$ is
an eigenvalue of ${\bm{J}}({\bm{0}})$ and
${\bm{u}}_{1[i]}{\bm{u}}_{1[j]}^{\top}+{\bm{u}}_{1[j]}{\bm{u}}_{1[i]}^{\top}$
is a corresponding eigenvector. All the other eigenvalues are $0$.

We simplify the notation by letting ${\bm{u}}_{1}:={\bm{u}}_{1[1]}$. A direct
corollary of Lemma 5.4 is that ${\bm{u}}_{1}{\bm{u}}_{1}^{\top}$ is the top
eigenvector of ${\bm{J}}({\bm{0}})$. According to Theorem 5.3, now there are
only two types of trajectories, which correspond to infinitesimal
initialization ${\bm{W}}_{\alpha}\to{\bm{0}}$ with positive alignment with
${\bm{u}}_{1}{\bm{u}}_{1}^{\top}$ or $-{\bm{u}}_{1}{\bm{u}}_{1}^{\top}$. As
the initialization must be PSD, ${\bm{W}}_{\alpha}\to{\bm{0}}$ cannot have
positive alignment with $-{\bm{u}}_{1}{\bm{u}}_{1}^{\top}$. For the former
case, Theorem 5.6 below states that, for every fixed time $t$, the GF solution
$\phi({\bm{W}}_{\alpha},T({\bm{W}}_{\alpha})+t)$ after shifting by a time
offset
$T({\bm{W}}_{\alpha}):=\frac{1}{2\mu_{1}}\log(\langle{\bm{W}}_{\alpha},{\bm{u}}_{1}{\bm{u}}_{1}^{\top}\rangle^{-1})$
converges to the GLRL solution ${\bm{W}}^{\mathrm{G}}_{1}(t)$ as
${\bm{W}}_{\alpha}\to{\bm{0}}$. The only assumption for this result is that
${\bm{0}}$ is not a minimizer of $f(\,\cdot\,)$ in ${\mathbb{S}}^{+}_{d}$
(which is equivalent to $\lambda_{1}(-\nabla f({\bm{0}}))>0$) and $-\nabla
f({\bm{0}})$ has an eigenvalue gap. In the full observation case, this
assumption is satisfied easily if the ground-truth matrix has a unique top
eigenvalue. The proof for Theorem 5.6 is deferred to Section G.1.

######  Assumption 5.5.

$\mu_{1}>\max\\{\mu_{2},0\\}$, where $\mu_{1}:=\lambda_{1}(-\nabla
f({\bm{0}})),\mu_{2}:=\lambda_{2}(-\nabla f({\bm{0}}))$.

######  Theorem 5.6.

Under 5.5, the following limit ${\bm{W}}^{\mathrm{G}}_{1}(t)$ exists and is a
solution of (2).

${\bm{W}}^{\mathrm{G}}_{1}(t):=\lim_{\epsilon\to
0}{\bm{W}}^{\mathrm{G}}_{1,\epsilon}\left(\tfrac{1}{2\mu_{1}}\log\tfrac{1}{\epsilon}+t\right)=\lim_{\epsilon\to
0}\phi\left(\epsilon{\bm{u}}_{1}{\bm{u}}_{1}^{\top},\tfrac{1}{2\mu_{1}}\log\tfrac{1}{\epsilon}+t\right).$
(9)

Let $\\{{\bm{W}}_{\alpha}\\}\subseteq{\mathbb{S}}^{+}_{d}$ be PSD matrices
converging to ${\bm{0}}$ with positive alignment with
${\bm{u}}_{1}{\bm{u}}_{1}^{\top}$ as $\alpha\to 0$, that is, $\lim_{\alpha\to
0}{\bm{W}}_{\alpha}={\bm{0}}$ and $\exists\alpha_{0},q>0$ such that
$\left<{\bm{W}}_{\alpha},{\bm{u}}_{1}{\bm{u}}_{1}^{\top}\right>\geq
q\left\|{\bm{W}}_{\alpha}\right\|_{\mathrm{F}}$ for all $\alpha<\alpha_{0}$.
Then $\forall t\in\mathbb{R}$, there is a constant $C>0$ such that

$\left\|\phi\left({\bm{W}}_{\alpha},\tfrac{1}{2\mu_{1}}\log\tfrac{1}{\left<{\bm{W}}_{\alpha},{\bm{u}}_{1}{\bm{u}}_{1}^{\top}\right>}+t\right)-{\bm{W}}^{\mathrm{G}}_{1}(t)\right\|_{\mathrm{F}}\leq
C\left\|{\bm{W}}_{\alpha}\right\|_{\mathrm{F}}^{\frac{\tilde{\gamma}}{2\mu_{1}+\tilde{\gamma}}}$
(10)

for every sufficiently small $\alpha$, where
$\tilde{\gamma}:=2\mu_{1}-(\mu_{1}+\mu_{2})=\mu_{1}-\mu_{2}$.

It is worth to note that ${\bm{W}}^{\mathrm{G}}_{1}(t)$ has rank $\leq 1$ for
any $t\in\mathbb{R}$, since every ${\bm{W}}^{\mathrm{G}}_{1,\epsilon}(t)$ has
rank $\leq 1$ and the set ${\mathbb{S}}^{+}_{d,\leq 1}$ is closed. This
matches with the first warmup example: GD does start learning with rank-1
solutions. Interestingly, in the case where the limit
$\overline{{}{\bm{W}}}_{1}:=\lim_{t\to+\infty}{\bm{W}}^{\mathrm{G}}_{1}(t)$
happens to be a minimizer of $f(\,\cdot\,)$ in ${\mathbb{S}}^{+}_{d}$, GLRL
should exit with the rank-1 solution $\overline{{}{\bm{W}}}_{1}$ after the
first phase, and the following theorem shows that this is also the solution
found by GF.

######  Assumption 5.7.

$f({\bm{W}})$ is locally analytic at each point.

######  Theorem 5.8.

Under Assumptions 5.5 and 5.7, if
$\|{\bm{W}}^{\mathrm{G}}_{1}(t)\|_{\mathrm{F}}$ is bounded for all $t\geq 0$,
then the limit
$\overline{{}{\bm{W}}}_{1}:=\lim_{t\to+\infty}{\bm{W}}^{\mathrm{G}}_{1}(t)$
exists. Further, if $\overline{{}{\bm{W}}}_{1}$ is a minimizer of
$f(\,\cdot\,)$ in ${\mathbb{S}}^{+}_{d}$, then for PSD matrices
$\\{{\bm{W}}_{\alpha}\\}\subseteq{\mathbb{S}}^{+}_{d}$ converging to
${\bm{0}}$ with positive alignment with ${\bm{u}}_{1}{\bm{u}}_{1}^{\top}$ as
$\alpha\to 0$, it holds that $\lim_{\alpha\to
0}\lim_{t\to+\infty}\phi({\bm{W}}_{\alpha},t)=\overline{{}{\bm{W}}}_{1}$.

5.7 is a natural assumption, since $f(\,\cdot\,)$ in most cases of matrix
factorization is a quadratic or polynomial function (e.g., matrix sensing,
matrix completion). In general, it is unlikely for a gradient-based
optimization process to get stuck at saddle points [Lee et al., 2017, Panageas
et al., 2019]. Thus, we should expect to see in general that GLRL finds the
rank-1 solution if the problem is feasible with rank-1 matrices. This means at
least for this subclass of problems, the implicit regularization of GD is
rather unrelated to norm minimization. Below is a concrete example:

######  Example 5.9 (Counter example of 1.1, Gunasekar et al. 2017).

Theorem 5.8 enables us to construct counterexamples of the _implicit nuclear
norm regularization_ conjecture in [Gunasekar et al., 2017]. The idea is to
construct a problem where every rank-1 stationary point of
$\mathcal{L}({\bm{U}})$ (i.e., $\nabla\mathcal{L}({\bm{U}})={\bm{0}}$ and
${\bm{U}}\in\mathbb{R}^{d\times d}$ is rank-1) attains the global minimum but
none of them is minimizing the nuclear norm. Below we give a concrete matrix
completion problem that meets the above requirement. Let ${\bm{M}}$ be a
partially observed matrix to be recovered, where the entries in
$\Omega=\\{(1,3),(1,4),(2,3),(3,1),(3,2),(4,1)\\}$ are observed and the others
(marked with “?”) are unobserved. The optimization problem is defined formally
by
$\mathcal{L}({\bm{U}})=\frac{1}{2}f({\bm{U}}{\bm{U}}^{\top}),f({\bm{W}})=\frac{1}{2}\sum_{(i,j)\in\Omega}(W_{ij}-M_{ij})^{2}$.

${\bm{M}}=\begin{bmatrix}?&?&1&R\\\ ?&?&R&?\\\ 1&R&?&?\\\
R&?&?&?\end{bmatrix},{\bm{M}}_{\mathrm{norm}}=\begin{bmatrix}R&1&1&R\\\
1&R&R&1\\\ 1&R&R&1\\\
R&1&1&R\end{bmatrix},{\bm{M}}_{\mathrm{rank}}=\begin{bmatrix}1&R&1&R\\\
R&R^{2}&R&R^{2}\\\ 1&R&1&R\\\ R&R^{2}&R&R^{2}\end{bmatrix}.$

Here $R>1$ is a large constant, e.g., $R=100$. The minimum nuclear norm
solution is the rank-2 matrix ${\bm{M}}_{\mathrm{norm}}$, which has
$\|{\bm{M}}_{\mathrm{norm}}\|_{\ast}=4R$ (which is $400$ when $R=100$).
${\bm{M}}_{\mathrm{rank}}$ is a rank-1 solution with much larger nuclear norm,
$\|{\bm{M}}_{\mathrm{norm}}\|_{\ast}=2R^{2}+2$ (which is $20002$ when
$R=100$). We can verify that $f(\,\cdot\,)$ satisfies Assumptions 5.5 and 5.7
and ${\bm{W}}^{\mathrm{G}}_{1}(t)$ converges to the rank-1 solution
${\bm{M}}_{\mathrm{rank}}$. Therefore, GF with infinitesimal initialization
converges to ${\bm{M}}_{\mathrm{rank}}$ rather than
${\bm{M}}_{\mathrm{norm}}$, which refutes the conjecture in [Gunasekar et al.,
2017]. See Appendix B for a formal statement.

###  5.3 Equivalence between GD and GLRL: General Case

Theorem 5.6 shows that for any _fixed_ time $t$, the trajectory of GLRL in the
first phase approximates GF with infinitesimal initialization, i.e.,
${\bm{W}}^{\mathrm{G}}_{1}(t)=\lim_{\alpha\to
0}\widehat{{\bm{W}}}_{\alpha}(t)$, where
$\widehat{{\bm{W}}}_{\alpha}(t):=\phi({\bm{W}}_{\alpha},\frac{1}{2\mu_{1}}\log(\langle{\bm{W}}_{\alpha},{\bm{u}}_{1}{\bm{u}}_{1}^{\top}\rangle^{-1})+t)$.
However, ${\bm{W}}^{\mathrm{G}}_{1}(\infty)\neq\lim_{\alpha\to
0}\widehat{{\bm{W}}}_{\alpha}(\infty)$ does not hold in general, unless the
prerequisite in Theorem 5.8 is satisfied, i.e., unless
$\overline{{}{\bm{W}}}_{1}={\bm{W}}^{\mathrm{G}}_{1}(\infty)$ is a minimizer
of $f(\,\cdot\,)$ in ${\mathbb{S}}^{+}_{d}$. This is because of the well-known
result that GD converges to local minimizers [Lee et al., 2016, 2017]. We
adapt Theorem 2 of Lee et al. [2017] to the setting of GF (Theorem G.5) and
obtain the following result (Theorem 5.10); see Section G.4 for the proof.

######  Theorem 5.10.

Let $f:\mathbb{R}^{d\times d}\to\mathbb{R}$ be a convex
$\mathcal{C}^{2}$-smooth function. (1). All stationary points of
$\mathcal{L}:\mathbb{R}^{d\times
d}\to\mathbb{R},\mathcal{L}({\bm{U}})=\frac{1}{2}f({\bm{U}}{\bm{U}}^{\top})$
are either strict saddles or global minimizers; (2). For any random
initialization, GF (1) converges to strict saddles of $\mathcal{L}({\bm{U}})$
with probability $0$.

Therefore, for convex $f(\,\cdot\,)$ such as matrix sensing and completion,
suppose $f(\,\cdot\,)$ has no rank-$1$ PSD minimizer, then no matter how small
$\alpha$ is, $\widehat{{\bm{W}}}_{\alpha}(\infty)$ (if exists) is a minimizer
of $f(\,\cdot\,)$ with a higher rank and thus away from the rank-1 matrix
$\overline{{}{\bm{W}}}_{1}$. In other words, ${\bm{W}}^{\mathrm{G}}_{1}(t)$
only describes the limiting trajectory of GF in the first phase, i.e., when GF
goes from near ${\bm{0}}$ to near $\overline{{}{\bm{W}}}_{1}$. After a
sufficiently long time (which depends on $\alpha$), GF escapes the critical
point $\overline{{}{\bm{W}}}_{1}$, but this part is not described by
${\bm{W}}^{\mathrm{G}}_{1}(t)$.

To understand how GF escapes $\overline{{}{\bm{W}}}_{1}$, a priori, we need to
know how GF approaches $\overline{{}{\bm{W}}}_{1}$. Using a similar argument
for Theorem 5.3, Theorem 5.11 shows that generically GF only escapes in the
direction of ${\bm{v}}_{1}{\bm{v}}_{1}^{\top}$, where ${\bm{v}}_{1}$ is the
(unique) top eigenvector of $-\nabla f(\overline{{}{\bm{W}}}_{1})$, and thus
the limiting trajectory exactly matches with that of GLRL in the second phase
until GF gets close to another critical point
$\overline{{}{\bm{W}}}_{2}\in{\mathbb{S}}^{+}_{d,\leq 2}$. If
$\overline{{}{\bm{W}}}_{2}$ is still not a minimizer of $f(\,\cdot\,)$ in
${\mathbb{S}}^{+}_{d}$ (but it is a local minimizer in
${\mathbb{S}}^{+}_{d,\leq 2}$ generically), then GF escapes
$\overline{{}{\bm{W}}}_{2}$ and the above process repeats until
$\overline{{}{\bm{W}}}_{K}$ is a minimizer in ${\mathbb{S}}^{+}_{d}$ for some
$K$. Here by “generically” we hide some technical assumptions and we elaborate
on them in Appendix H. See Figure 1 and Figure 2 for experimental verification
of the equivalence between GD and GLRL. We end this section with the following
characterization of GF:

######  Theorem 5.11 (Theorem G.2, informal).

Let $\overline{{}{\bm{W}}}$ be a critical point of (2) satisfying that
$\overline{{}{\bm{W}}}$ is a local minimizer of $f(\,\cdot\,)$ in
${\mathbb{S}}^{+}_{d,\leq r}$ for some $r\geq 1$ but not a minimizer in
${\mathbb{S}}^{+}_{d}$. Let $-\nabla
f(\overline{{}{\bm{W}}})=\sum_{i=1}^{d}\mu_{i}{\bm{v}}_{i}{\bm{v}}_{i}^{\top}$
be the eigendecomposition of $-\nabla f(\overline{{}{\bm{W}}})$. If
$\mu_{1}>\mu_{2}$ and if there exists time $T_{\alpha}\in\mathbb{R}$ for every
$\alpha$ so that $\phi({\bm{W}}_{\alpha},T_{\alpha})$ converges to
$\overline{{}{\bm{W}}}$ with positive alignment with the top principal
component ${\bm{v}}_{1}{\bm{v}}_{1}^{\top}$ as $\alpha\to 0$, then for every
fixed $t$, $\lim\limits_{\alpha\to
0}\phi({\bm{W}}_{\alpha},T_{\alpha}+\frac{1}{2\mu_{1}}\log\frac{1}{\left<\phi({\bm{W}}_{\alpha},T_{\alpha}),{\bm{v}}_{1}{\bm{v}}_{1}^{\top}\right>}+t)$
exists and is equal to ${\bm{W}}^{\mathrm{G}}(t):=\lim_{\epsilon\to
0}\phi(\overline{{}{\bm{W}}}+\epsilon{\bm{v}}_{1}{\bm{v}}_{1}^{\top},\frac{1}{2\mu_{1}}\log\frac{1}{\epsilon}+t)$.

##### Characterization of the trajectory of GF.

Generically, the trajectory of GF with small initialization can be split into
$K$ phases by $K+1$ critical points of (2),
$\\{\overline{{}{\bm{W}}}_{r}\\}_{r=0}^{K}$
($\overline{{}{\bm{W}}}_{0}={\bm{0}}$), where in phase $r$ GF escapes from
$\overline{{}{\bm{W}}}_{r-1}$ in the direction of the top principal component
of $-\nabla f(\overline{{}{\bm{W}}}_{r-1})$ and gets close to
$\overline{{}{\bm{W}}}_{r}$. Each $\overline{{}{\bm{W}}}_{r}$ is a local
minimizer of $f(\,\cdot\,)$ in ${\mathbb{S}}^{+}_{d,\leq r}$, but none of them
is a minimizer of $f(\,\cdot\,)$ in ${\mathbb{S}}^{+}_{d}$ except
$\overline{{}{\bm{W}}}_{K}$. The smaller the initialization is, the longer GF
stays around each $\overline{{}{\bm{W}}}_{r}$. Moreover,
$\\{\overline{{}{\bm{W}}}_{r}\\}_{r=0}^{K}$ corresponds to
$\\{\overline{{}{\bm{W}}}_{r,{\epsilon}}\\}_{r=0}^{K}$ in Definition 5.1 with
infinitesimal ${\epsilon}>0$.

Figure 1: The trajectory of depth-$2$ GD, ${\bm{W}}_{\text{GD}}(t)$, converges
to the trajectory of GLRL, ${\bm{W}}_{\text{GLRL}}(t)$, as the initialization
scale goes to $0$. We plot
$\mathrm{dist}(t)=\min_{t^{\prime}\in\mathcal{T}}\|{\bm{W}}_{\text{GD}}(t)-{\bm{W}}_{\text{GLRL}}(t^{\prime})\|_{\mathrm{F}}$
for different initialization scale $\|{\bm{W}}(0)\|_{\mathrm{F}}$, where
$\mathcal{T}$ is a discrete subset of $\mathbb{R}$ that $\delta$-covers the
entire trajectory of GLRL:
$\max_{t}\min_{t^{\prime}\in\mathcal{T}}\left\|{\bm{W}}_{\text{GLRL}}(t)-{\bm{W}}_{\text{GLRL}}(t^{\prime})\right\|_{\mathrm{F}}\leq\delta$
for $\delta\approx 0.00042$. For each $\|{\bm{W}}(0)\|_{\mathrm{F}}$, we run
$20$ random seeds and plot them separately. The ground truth
${\bm{W}}^{*}\in\mathbb{R}^{20\times 20}$ is a randomly generated rank-$3$
matrix with $\|{\bm{W}}^{*}\|_{\mathrm{F}}=20$. $30\%$ entries are observed.
See more in Section C.1.

##  6 Benefits of Depth: A View from GLRL

In this section, we consider matrix factorization problems with depth $L\geq
3$. Our goal is to understand the effect of the depth-$L$ parametrization
${\bm{W}}={\bm{U}}_{1}{\bm{U}}_{2}\cdots{\bm{U}}_{L}$ on the implicit bias —
how does depth encourage GF to find low rank solutions? We take the standard
assumption in existing analysis for the end-to-end dynamics that the weight
matrices have a balanced initialization, i.e.
${\bm{U}}_{i}^{\top}(0){\bm{U}}_{i}(0)={\bm{U}}_{i+1}(0){\bm{U}}_{i+1}^{\top}(0),\
\forall 1\leq i\leq L-1$. Arora et al. [2018] showed that if
$\\{{\bm{U}}_{i}\\}_{i=1}^{L}$ is balanced at initialization, then we have the
following end-to-end dynamics. Similar to the depth-2 case, we use
$\phi({\bm{W}}(0),t)$ to denote ${\bm{W}}(t)$, where

$\frac{\textup{{d}}{\bm{W}}}{\textup{{d}}t}=-\sum\nolimits_{i=0}^{L-1}({\bm{W}}{\bm{W}}^{\top})^{\frac{i}{L}}\nabla
f({\bm{W}})({\bm{W}}^{\top}{\bm{W}})^{1-\frac{i+1}{L}}.$ (11)

The lemma below is the foundation of our analysis for the deep case, which
greatly simplifies (11). Due to the space limit, we defer its derivations and
applications into Appendix I.

######  Lemma 6.1.

If ${\bm{W}}(t)$ is a symmetric solution of (11), then for
${\bm{M}}(t):={\bm{W}}(t)^{2/L}$, we have

$\frac{\textup{{d}}{\bm{M}}}{\textup{{d}}t}=-\nabla
f({\bm{M}}^{L/2}){\bm{M}}^{L/2}-{\bm{M}}^{L/2}\nabla f({\bm{M}}^{L/2}).$ (12)

parameter : step size $\eta>0$; small $\epsilon>0$

$\epsilon^{\prime}\leftarrow\epsilon^{1/L}$,
$\mathcal{L}\left({\bm{U}}_{1},\cdots,{\bm{U}}_{L}\right):=f({\bm{W}}_{1}\cdots{\bm{W}}_{L})$.

${\bm{W}}_{0}\leftarrow{\bm{0}}\in\mathbb{R}^{d\times d}$, and
${\bm{U}}_{0,1}(\infty),\dots,{\bm{U}}_{0,L}(\infty)\in\mathbb{R}^{d\times 0}$
are empty matrices

while  _$\lambda_{1}(-\nabla f\left({\bm{W}}_{r}\right)) >0$_ do

$r\leftarrow r+1$

let ${\bm{u}}_{r}$ be a top (unit) eigenvector of $-\nabla
f\left({\bm{W}}_{r-1}\right)$

${\bm{U}}_{r,1}(0)\leftarrow\begin{bmatrix}{\bm{U}}_{r-1,1}(\infty)&\epsilon^{\prime}{\bm{u}}_{r}\end{bmatrix}\in\mathbb{R}^{d\times
r}$

${\bm{U}}_{r,k}(0)\leftarrow\begin{bmatrix}{\bm{U}}_{r-1,k}(\infty)&{\bm{0}}\\\
{\bm{0}}&\epsilon^{\prime}\end{bmatrix}\in\mathbb{R}^{r\times r}$ for all
$2\leq k\leq L-1$

${\bm{U}}_{r,L}(0)\leftarrow\begin{bmatrix}{\bm{U}}_{r-1,L}(\infty)\\\
\epsilon^{\prime}{\bm{u}}_{r}^{\top}\end{bmatrix}\in\mathbb{R}^{r\times d}$

for  _$t=0,1,\ldots$_ do

${\bm{U}}_{r,i}(t+1)\leftarrow{\bm{U}}_{r,i}(t)-\eta\nabla_{{\bm{U}}_{i}}\mathcal{L}\left({\bm{U}}_{r,1}(t),\cdots,{\bm{U}}_{r,L}(t)\right)$,
$\forall 1\leq i\leq L$.

${\bm{W}}_{r}\leftarrow{\bm{U}}_{r,1}(\infty)\cdots{\bm{U}}_{r,L}(\infty)$

return ${\bm{W}}_{r}$

Algorithm 2 Deep Greedy Low-Rank Learning (Deep GLRL)

Our main result, Theorem 6.2, gives a characterization of the limiting
trajectory for deep matrix factorization with infinitesimal identity
initialization. Here $\overline{{}{\bm{W}}}(t):=\lim_{\alpha\to
0}{\bm{W}}^{\mathrm{G}}_{\alpha}(t)$ is the trajectory of deep GLRL, where
${\bm{W}}^{\mathrm{G}}_{\alpha}(t):=\phi(\alpha{\bm{e}}_{1}{\bm{e}}_{1}^{\top},\frac{{\alpha}^{-(1-\frac{1}{P})}}{2\mu_{1}(P-1)}+t)$
(see Algorithm 2). The dynamics for general initialization is more
complicated. Please see discussions in Appendix J.

######  Theorem 6.2.

Let $P=\frac{L}{2}$, $L\geq 3$. Suppose $\left\|\nabla
f({\bm{0}})\right\|_{2}=\lambda_{1}(-\nabla
f({\bm{0}}))>\max\\{\lambda_{2}(-\nabla f({\bm{0}})),0\\}$,222$\left\|\nabla
f({\bm{0}})\right\|_{2}=\lambda_{1}(-\nabla f({\bm{0}}))$ is a technical
assumption which we believe could be removed with a more refined analysis.

$\text{for every fixed
}t\in\mathbb{R},\quad\left\|\phi\left({\alpha}{\bm{I}},\tfrac{{\alpha}^{-(1-1/P)}}{2\mu_{1}(P-1)}+t\right)-\overline{{}{\bm{W}}}(t)\right\|_{\mathrm{F}}=O(\alpha^{\frac{1}{P(P+1)}}),$
(13)

and for any $2\leq k\leq d$,

$\text{for every fixed
}t\in\mathbb{R},\quad\lambda_{k}\left(\phi\left({\alpha}{\bm{I}},\tfrac{{\alpha}^{-(1-1/P)}}{2\mu_{1}(P-1)}+t\right)\right)=O(\alpha).$
(14)

##### So how does depth encourage GF to find low-rank solutions?

When the ground truth is low-rank, say rank-$k$, our experiments (Figure 2)
suggest that GF with small initialization deep matrix sensing finds solutions
with smaller $k$-low-rankness compared to the depth-2 case, thus achieving
better generalization. At first glance, this is contradictory to what Theorem
6.2 suggests, i.e., the convergence rate of deep GLRL at a constant time gets
slower as the depth increases. However, it turns out the uniform upper bound
for the distance between GF and GLRL is not the ideal metric for the eventual
$k$-low-rankness of learned solution. Below we will illustrate why _the
$r$-low-rankness of GF within each phase $r$_ is a better metric and how they
are different.

######  Definition 6.3 ($r$-low-rankness).

For matrix ${\bm{M}}\in\mathbb{R}^{d\times d}$, we define the _$r$ -low-
rankness_ of ${\bm{M}}$ as $\sqrt{\sum_{i=r+1}^{d}\sigma^{2}_{i}({\bm{M}})}$,
where $\sigma_{i}({\bm{M}})$ is the $i$-th largest singular value of
${\bm{M}}$.

Suppose $f(\,\cdot\,)$ admits a unique minimizer ${\bm{W}}_{0}$ in
${\mathbb{S}}^{+}_{d,1}$, and we run GF from $\alpha{\bm{I}}$ for both
depth-$2$ and depth-$L$ cases. Intuitively, the $1$-low-rankness of the
depth-2 solution is $\Omega(\alpha^{1-\mu_{2}/\mu_{1}})$, which can be seen
from the second warmup example in Section 4. For the depth-$L$ solution,
though it may diverge from the trajectory of deep GLRL more than the depth-2
solution does, its $1$-low-rankness is only $O(\alpha)$, as shown in Theorem
6.4. The key idea is to show that there is a basin in the manifold of rank-$1$
matrices around ${\bm{W}}_{0}$ such that any GF starting within the basin
converges to ${\bm{W}}_{0}$. Based on this, we can prove that starting from
any matrix $O(\alpha)$-close to the basin, GF converges to a solution
$O(\alpha)$-close to ${\bm{W}}_{0}$. See Appendix K for more details.

######  Theorem 6.4.

In the same settings as Theorem 6.2, if $\overline{{}{\bm{W}}}(\infty)$ exists
and is a minimizer of $f(\,\cdot\,)$ in ${\mathbb{S}}^{+}_{d,\leq 1}$, under
the additional regularity assumption K.1, we have

$\inf_{t\in\mathbb{R}}\left\|\phi\left({\alpha}{\bm{I}},t\right)-\overline{{}{\bm{W}}}(\infty)\right\|_{F}=O(\alpha).$
(15)

##### Interpretation for the advantage of depth with multiple phases.

For depth-2 GLRL, the low-rankness is raised to some power less then $1$ per
phase (depending on the eigengap). For deep GLRL, we show the low-rankness is
only multiplied by some constant for the first phase and speculate it to be
true for later phases. This conjecture is supported by our experiments; see
Figure 2. Interestingly, our theory and experiments (Figure 5) suggest that
while being deep is good for generalization, being much deeper may not be much
better: once $L\geq 3$, increasing the depth does not improve the order of
low-rankness significantly. While this theoretic result is only for identity
initialization, Theorem D.1 and Corollary D.2 further show that the dynamics
of GF (11) with any initialization pointwise converges as $L\to\infty$, under
a suitable time rescaling. See Figure 6 for experimental verification.

Figure 2: GD passes by the same set of critical points as GLRL when the
initialization scale is small, and gets much closer to the critical points
when $L\geq 3$. Depth-$2$ GD requires a much smaller initialization scale to
maintain small low-rankness. Here the ground truth matrix
${\bm{W}}^{*}\in\mathbb{R}^{20\times 20}$ is of rank $3$ as stated in Section
C.1. In this case, GLRL has $3$ phases and $4$ critical points
$\\{\overline{{}{\bm{W}}}_{r}\\}_{r=0}^{3}$, where
$\overline{{}{\bm{W}}}_{0}={\bm{0}}$ and
$\overline{{}{\bm{W}}}_{3}={\bm{W}}^{*}$. For each depth $L$ and
initialization scale $\|{\bm{W}}(0)\|_{\mathrm{F}}$, we plot the distance
between the current step of GD and the closest critical point of GLRL,
$\min_{0\leq r\leq
3}\|{\bm{W}}_{\text{GD}}(t)-\overline{{}{\bm{W}}}_{r}\|_{\mathrm{F}}$, the
norm of full gradient,
$\|\nabla_{{\bm{U}}_{1:L}}\mathcal{L}({\bm{U}}_{1:L})\|_{\mathrm{F}}$ and the
$(r+1)$-low-rankness of ${\bm{W}}_{\text{GD}}(t)$ with
$r:=\operatorname*{arg\,min}_{0\leq r\leq
3}\|{\bm{W}}_{\text{GD}}(t)-\overline{{}{\bm{W}}}_{r}\|_{\mathrm{F}}$.

##  7 Conclusion and Future Directions

In this work, we connect gradient descent to Greedy Low-Rank Learning (GLRL)
to explain the success of using gradient descent to find low-rank solutions in
the matrix factorization problem. This enables us to construct counterexamples
to the implicit nuclear norm conjecture in Gunasekar et al. [2017]. Taking the
view of GLRL can also help us understand the benefits of depth.

Our result on the equivalence between gradient flow with infinitesimal
initialization and GLRL is based on some regularity conditions that we expect
to hold generically. We leave it a future work to justify these condition,
possibly through a smoothed analysis on the objective $f(\,\cdot\,)$. Another
interesting future direction is to find the counterpart of GLRL in training
deep neural nets. This could be one way to go beyond the view of norm
minimization in the study of the implicit regularization of gradient descent.

#### Acknowledgments

The authors thank Sanjeev Arora and Jason D. Lee for helpful discussions. The
authors also thank Runzhe Wang for useful suggestions on writing. ZL and YL
acknowledge support from NSF, ONR, Simons Foundation, Schmidt Foundation,
Mozilla Research, Amazon Research, DARPA and SRC. ZL is also supported by
Microsoft PhD Fellowship.

## References

  * Agrawal et al. [2018] Akshay Agrawal, Robin Verschueren, Steven Diamond, and Stephen Boyd.  A rewriting system for convex optimization problems.  _Journal of Control and Decision_ , 5(1):42–60, 2018.
  * Arora et al. [2018] Sanjeev Arora, N Cohen, and Elad Hazan.  On the optimization of deep networks: Implicit acceleration by overparameterization.  In _35th International Conference on Machine Learning_ , 2018.
  * Arora et al. [2019a] Sanjeev Arora, Nadav Cohen, Wei Hu, and Yuping Luo.  Implicit regularization in deep matrix factorization.  In H. Wallach, H. Larochelle, A. Beygelzimer, F. d’ Alché-Buc, E. Fox, and R. Garnett, editors, _Advances in Neural Information Processing Systems 32_ , pages 7411–7422. Curran Associates, Inc., 2019a.
  * Arora et al. [2019b] Sanjeev Arora, Simon S Du, Wei Hu, Zhiyuan Li, Russ R Salakhutdinov, and Ruosong Wang.  On exact computation with an infinitely wide neural net.  In H. Wallach, H. Larochelle, A. Beygelzimer, F. d’ Alché-Buc, E. Fox, and R. Garnett, editors, _Advances in Neural Information Processing Systems 32_ , pages 8139–8148. Curran Associates, Inc., 2019b.
  * Belabbas [2020] Mohamed Ali Belabbas.  On implicit regularization: Morse functions and applications to matrix factorization.  _arXiv preprint arXiv:2001.04264_ , 2020.
  * Bezanson et al. [2012] Jeff Bezanson, Stefan Karpinski, Viral B Shah, and Alan Edelman.  Julia: A fast dynamic language for technical computing.  _arXiv preprint arXiv:1209.5145_ , 2012.
  * Chi et al. [2019] 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.
  * Chizat and Bach [2020] Lénaïc Chizat and Francis Bach.  Implicit bias of gradient descent for wide two-layer neural networks trained with the logistic loss.  volume 125 of _Proceedings of Machine Learning Research_ , pages 1305–1338. PMLR, 09–12 Jul 2020.
  * Chizat et al. [2019] Lénaïc Chizat, Edouard Oyallon, and Francis Bach.  On lazy training in differentiable programming.  In H. Wallach, H. Larochelle, A. Beygelzimer, F. d'Alché-Buc, E. Fox, and R. Garnett, editors, _Advances in Neural Information Processing Systems 32_ , pages 2937–2947. Curran Associates, Inc., 2019.
  * Clarke et al. [2008] Francis H. Clarke, Yuri S. Ledyaev, Ronald J. Stern, and Peter R. Wolenski.  _Nonsmooth analysis and control theory_ , volume 178.  Springer Science & Business Media, 2008.
  * Clarke [1975] Frank H. Clarke.  Generalized gradients and applications.  _Transactions of the American Mathematical Society_ , 205:247–262, 1975.
  * Clarke [1990] Frank H Clarke.  _Optimization and Nonsmooth Analysis_.  Society for Industrial and Applied Mathematics, 1990.  doi: 10.1137/1.9781611971309.
  * Diamond and Boyd [2016] Steven Diamond and Stephen Boyd.  CVXPY: A Python-embedded modeling language for convex optimization.  _Journal of Machine Learning Research_ , 17(83):1–5, 2016.
  * Du and Lee [2018] Simon Du and Jason Lee.  On the power of over-parametrization in neural networks with quadratic activation.  In Jennifer Dy and Andreas Krause, editors, _Proceedings of the 35th International Conference on Machine Learning_ , volume 80 of _Proceedings of Machine Learning Research_ , pages 1329–1338, Stockholmsmässan, Stockholm Sweden, 10–15 Jul 2018. PMLR.
  * Gidel et al. [2019] Gauthier Gidel, Francis Bach, and Simon Lacoste-Julien.  Implicit regularization of discrete gradient dynamics in linear neural networks.  In H. Wallach, H. Larochelle, A. Beygelzimer, F. d’ Alché-Buc, E. Fox, and R. Garnett, editors, _Advances in Neural Information Processing Systems 32_ , pages 3196–3206. Curran Associates, Inc., 2019.
  * Gissin et al. [2020] Daniel Gissin, Shai Shalev-Shwartz, and Amit Daniely.  The implicit bias of depth: How incremental learning drives generalization.  In _International Conference on Learning Representations_ , 2020.
  * Gunasekar et al. [2017] Suriya Gunasekar, Blake E Woodworth, Srinadh Bhojanapalli, Behnam Neyshabur, and Nati Srebro.  Implicit regularization in matrix factorization.  In I. Guyon, U. V. Luxburg, S. Bengio, H. Wallach, R. Fergus, S. Vishwanathan, and R. Garnett, editors, _Advances in Neural Information Processing Systems 30_ , pages 6151–6159. Curran Associates, Inc., 2017.
  * Gunasekar et al. [2018] Suriya Gunasekar, Jason D Lee, Daniel Soudry, and Nati Srebro.  Implicit bias of gradient descent on linear convolutional networks.  In S. Bengio, H. Wallach, H. Larochelle, K. Grauman, N. Cesa-Bianchi, and R. Garnett, editors, _Advances in Neural Information Processing Systems 31_ , pages 9482–9491. Curran Associates, Inc., 2018.
  * Haeffele and Vidal [2019] Benjamin D. Haeffele and René Vidal.  Structured Low-Rank Matrix Factorization: Global Optimality, Algorithms, and Applications.  _IEEE Transactions on Pattern Analysis and Machine Intelligence (PAMI)_ , 42(6):1468–1482, 2019.
  * Hiriart-Urruty and Lewis [1999] Jean-Baptiste Hiriart-Urruty and A. S. Lewis.  The clarke and michel-penot subdifferentials of the eigenvalues of a symmetric matrix.  _Comput. Optim. Appl._ , 13(1-3):13–23, 1999\.  doi: 10.1023/A:1008644520093.
  * Jacot et al. [2018] Arthur Jacot, Franck Gabriel, and Clement Hongler.  Neural tangent kernel: Convergence and generalization in neural networks.  In S. Bengio, H. Wallach, H. Larochelle, K. Grauman, N. Cesa-Bianchi, and R. Garnett, editors, _Advances in Neural Information Processing Systems 31_ , pages 8571–8580. Curran Associates, Inc., 2018.
  * Ji and Telgarsky [2019a] Ziwei Ji and Matus Telgarsky.  Gradient descent aligns the layers of deep linear networks.  In _International Conference on Learning Representations_ , 2019a.
  * Ji and Telgarsky [2019b] Ziwei Ji and Matus Telgarsky.  A refined primal-dual analysis of the implicit bias.  _arXiv preprint arXiv:1906.04540_ , 2019b.
  * Khanna et al. [2017] Rajiv Khanna, Ethan Elenberg, Alexandros G Dimakis, and Sahand Negahban.  On approximation guarantees for greedy low rank optimization.  _arXiv preprint arXiv:1703.02721_ , 2017.
  * Lee et al. [2016] Jason D Lee, Max Simchowitz, Michael I Jordan, and Benjamin Recht.  Gradient descent only converges to minimizers.  In _Conference on learning theory_ , pages 1246–1257, 2016.
  * Lee et al. [2017] Jason D Lee, Ioannis Panageas, Georgios Piliouras, Max Simchowitz, Michael I Jordan, and Benjamin Recht.  First-order methods almost always avoid saddle points.  _arXiv preprint arXiv:1710.07406_ , 2017.
  * Li et al. [2018] Yuanzhi Li, Tengyu Ma, and Hongyang Zhang.  Algorithmic regularization in over-parameterized matrix sensing and neural networks with quadratic activations.  In Sébastien Bubeck, Vianney Perchet, and Philippe Rigollet, editors, _Proceedings of the 31st Conference On Learning Theory_ , volume 75 of _Proceedings of Machine Learning Research_ , pages 2–47. PMLR, 06–09 Jul 2018.
  * Lyu and Li [2020] Kaifeng Lyu and Jian Li.  Gradient descent maximizes the margin of homogeneous neural networks.  In _International Conference on Learning Representations_ , 2020.
  * Nacson et al. [2019a] Mor Shpigel Nacson, Suriya Gunasekar, Jason Lee, Nathan Srebro, and Daniel Soudry.  Lexicographic and depth-sensitive margins in homogeneous and non-homogeneous deep models.  In Kamalika Chaudhuri and Ruslan Salakhutdinov, editors, _Proceedings of the 36th International Conference on Machine Learning_ , volume 97 of _Proceedings of Machine Learning Research_ , pages 4683–4692, Long Beach, California, USA, 09–15 Jun 2019a. PMLR.
  * Nacson et al. [2019b] Mor Shpigel Nacson, Jason Lee, Suriya Gunasekar, Pedro Henrique Pamplona Savarese, Nathan Srebro, and Daniel Soudry.  Convergence of gradient descent on separable data.  In Kamalika Chaudhuri and Masashi Sugiyama, editors, _Proceedings of Machine Learning Research_ , volume 89 of _Proceedings of Machine Learning Research_ , pages 3420–3428. PMLR, 16–18 Apr 2019b.
  * Nacson et al. [2019c] Mor Shpigel Nacson, Nathan Srebro, and Daniel Soudry.  Stochastic gradient descent on separable data: Exact convergence with a fixed learning rate.  In Kamalika Chaudhuri and Masashi Sugiyama, editors, _Proceedings of Machine Learning Research_ , volume 89 of _Proceedings of Machine Learning Research_ , pages 3051–3059. PMLR, 16–18 Apr 2019c.
  * Panageas et al. [2019] Ioannis Panageas, Georgios Piliouras, and Xiao Wang.  First-order methods almost always avoid saddle points: The case of vanishing step-sizes.  In H. Wallach, H. Larochelle, A. Beygelzimer, F. d'Alché-Buc, E. Fox, and R. Garnett, editors, _Advances in Neural Information Processing Systems 32_ , pages 6474–6483. Curran Associates, Inc., 2019.
  * Paszke et al. [2019] Adam Paszke, Sam Gross, Francisco Massa, Adam Lerer, James Bradbury, Gregory Chanan, Trevor Killeen, Zeming Lin, Natalia Gimelshein, Luca Antiga, Alban Desmaison, Andreas Kopf, Edward Yang, Zachary DeVito, Martin Raison, Alykhan Tejani, Sasank Chilamkurthy, Benoit Steiner, Lu Fang, Junjie Bai, and Soumith Chintala.  Pytorch: An imperative style, high-performance deep learning library.  In H. Wallach, H. Larochelle, A. Beygelzimer, F. d'Alché-Buc, E. Fox, and R. Garnett, editors, _Advances in Neural Information Processing Systems 32_ , pages 8024–8035. Curran Associates, Inc., 2019.
  * Perko [2013] Lawrence Perko.  _Differential equations and dynamical systems_ , volume 7.  Springer Science & Business Media, 2013.
  * Razin and Cohen [2020] Noam Razin and Nadav Cohen.  Implicit regularization in deep learning may not be explainable by norms.  _arXiv preprint arXiv:2005.06398_ , 2020.
  * Shalev-Shwartz and Singer [2010] Shai Shalev-Shwartz and Yoram Singer.  On the equivalence of weak learnability and linear separability: New relaxations and efficient boosting algorithms.  _Machine learning_ , 80(2-3):141–163, 2010.
  * Soudry et al. [2018a] Daniel Soudry, Elad Hoffer, Mor Shpigel Nacson, Suriya Gunasekar, and Nathan Srebro.  The implicit bias of gradient descent on separable data.  _Journal of Machine Learning Research_ , 19(70):1–57, 2018a.
  * Soudry et al. [2018b] Daniel Soudry, Elad Hoffer, and Nathan Srebro.  The implicit bias of gradient descent on separable data.  In _International Conference on Learning Representations_ , 2018b.
  * Tieleman and Hinton [2012] Tijmen Tieleman and Geoffrey Hinton.  Lecture 6.5-rmsprop: Divide the gradient by a running average of its recent magnitude.  _COURSERA: Neural networks for machine learning_ , 4(2):26–31, 2012.
  * Wang et al. [2014] Zheng Wang, Ming-Jun Lai, Zhaosong Lu, Wei Fan, Hasan Davulcu, and Jieping Ye.  Rank-one matrix pursuit for matrix completion.  In _International Conference on Machine Learning_ , pages 91–99, 2014\.
  * Wilson et al. [2017] Ashia C Wilson, Rebecca Roelofs, Mitchell Stern, Nati Srebro, and Benjamin Recht.  The marginal value of adaptive gradient methods in machine learning.  In I. Guyon, U. V. Luxburg, S. Bengio, H. Wallach, R. Fergus, S. Vishwanathan, and R. Garnett, editors, _Advances in Neural Information Processing Systems 30_ , pages 4148–4158. Curran Associates, Inc., 2017.
  * Yao and Kwok [2016] Quanming Yao and James Tin Yau Kwok.  Greedy learning of generalized low-rank models.  In _IJCAI International Joint Conference on Artificial Intelligence_ , 2016.
  * Zhang et al. [2017] Chiyuan Zhang, Samy Bengio, Moritz Hardt, Benjamin Recht, and Oriol Vinyals.  Understanding deep learning requires rethinking generalization.  In _International Conference on Learning Representations_ , 2017.
  * Łojasiewicz [1965] Stanisław Łojasiewicz.  Ensembles semi-analytiques.  _IHES notes_ , 1965.

##  Appendix A Preliminary Lemmas

######  Lemma A.1.

For ${\bm{U}}_{0}\in\mathbb{R}^{d\times r}$ and
${\bm{W}}_{0}:={\bm{U}}_{0}{\bm{U}}_{0}^{\top}$, the following statements are
equivalent:

  1. (1).

${\bm{U}}_{0}$ is a stationary point of
$\mathcal{L}({\bm{U}})=\frac{1}{2}f({\bm{U}}{\bm{U}}^{\top})$;

  2. (2).

$\nabla f({\bm{W}}_{0}){\bm{W}}_{0}={\bm{0}}$;

  3. (3).

${\bm{W}}_{0}:={\bm{U}}_{0}{\bm{U}}_{0}^{\top}$ is a critical point of (2).

###### Proof.

(2) $\Rightarrow$ (3) is trivial. We only prove (1) $\Rightarrow$ (2), (3)
$\Rightarrow$ (1).

##### Proof for (1) ${\bm{\Rightarrow}}$ (2).

If ${\bm{U}}_{0}$ is a stationary point, then
${\bm{0}}=\nabla\mathcal{L}({\bm{U}}_{0})=\nabla f({\bm{W}}_{0}){\bm{U}}_{0}$.
So

$\nabla f({\bm{W}}_{0}){\bm{W}}_{0}=\left(\nabla
f({\bm{W}}_{0}){\bm{U}}_{0}\right){\bm{U}}_{0}^{\top}={\bm{0}}.$

##### Proof for (3) ${\bm{\Rightarrow}}$ (1).

If ${\bm{W}}_{0}$ is a critical point, then

$0=\left<{\bm{g}}({\bm{W}}_{0}),\nabla
f({\bm{W}}_{0})\right>=-2\operatorname{Tr}(\nabla
f({\bm{W}}_{0}){\bm{W}}_{0}\nabla f({\bm{W}}_{0}))=-2\|\nabla
f({\bm{W}}_{0}){\bm{U}}_{0}\|_{\mathrm{F}}^{2},$

which implies $\nabla\mathcal{L}({\bm{U}}_{0})={\bm{0}}$. ∎

######  Lemma A.2.

For a stationary point ${\bm{U}}_{0}\in\mathbb{R}^{d\times r}$ of
$\mathcal{L}({\bm{U}})=\frac{1}{2}f({\bm{U}}{\bm{U}}^{\top})$ where
$f(\,\cdot\,)$ is convex, ${\bm{W}}_{0}:={\bm{U}}_{0}{\bm{U}}_{0}^{\top}$
attains the global minimum of $f(\,\cdot\,)$ in
${\mathbb{S}}^{+}_{d}:=\\{{\bm{W}}:{\bm{W}}\succeq{\bm{0}}\\}$ iff $\nabla
f({\bm{W}}_{0})\succeq{\bm{0}}$.

###### Proof.

Since $f({\bm{W}})$ is a convex function and ${\mathbb{S}}^{+}_{d}$ is convex,
we know that ${\bm{W}}_{0}$ is a global minimizer of $f({\bm{W}})$ in
${\mathbb{S}}^{+}_{d}$ iff

$\left<\nabla f({\bm{W}}_{0}),{\bm{W}}-{\bm{W}}_{0}\right>\geq
0,\qquad\forall{\bm{W}}\succeq{\bm{0}}.$ (16)

Note that $\left<\nabla
f({\bm{W}}_{0}),{\bm{W}}_{0}\right>=\operatorname{Tr}(\nabla
f({\bm{W}}_{0}){\bm{W}}_{0})$. By Lemma A.1, $\left<\nabla
f({\bm{W}}_{0}),{\bm{W}}_{0}\right>=0$. Combining this with (16), we know that
${\bm{W}}_{0}$ is a global minimizer iff

$\left<\nabla f({\bm{W}}_{0}),{\bm{W}}\right>\geq
0,\qquad\forall{\bm{W}}\succeq{\bm{0}}.$ (17)

It is easy to check that this condition is equivalent to $\nabla
f({\bm{W}}_{0})\succeq{\bm{0}}$. ∎

##  Appendix B Proofs for Counter-example

######  Conjecture B.1 (Formal Statement, Gunasekar et al. [2017]).

Suppose $f:\mathbb{R}^{d\times d}\to\mathbb{R}$ is a quadratic function and
$\min\limits_{{\bm{W}}\succeq\bm{0}}f({\bm{W}})=0$. Then for any
${\bm{W}}_{\mathrm{init}}\succ\bm{0}$ if
$\overline{{}{\bm{W}}}_{1}=\lim\limits_{\alpha\to
0}\lim\limits_{t\to+\infty}\phi(\alpha{\bm{W}}_{\mathrm{init}},t)$ exists and
$f(\overline{{}{\bm{W}}}_{1})=0$, then
$\|\overline{{}{\bm{W}}}_{1}\|_{\ast}=\min\limits_{{\bm{W}}\succeq{\bm{0}}}\|{\bm{W}}\|_{\ast}~{}~{}\text{s.t.}~{}~{}f({\bm{W}})=0$.

######  Propsition B.2 (Formal Statement for Example 5.9).

For constant $R>1$, let

${\bm{M}}=\begin{bmatrix}?&?&1&R\\\ ?&?&R&?\\\ 1&R&?&?\\\
R&?&?&?\end{bmatrix},{\bm{M}}_{\mathrm{norm}}=\begin{bmatrix}R&1&1&R\\\
1&R&R&1\\\ 1&R&R&1\\\ R&1&1&R\end{bmatrix},\textrm{ and
}{\bm{M}}_{\mathrm{rank}}=\begin{bmatrix}1&R&1&R\\\ R&R^{2}&R&R^{2}\\\
1&R&1&R\\\ R&R^{2}&R&R^{2}\end{bmatrix}.$

and

$\mathcal{L}({\bm{U}})=\frac{1}{2}f({\bm{U}}{\bm{U}}^{\top}),\quad
f({\bm{W}})=\frac{1}{2}\sum_{(i,j)\in\Omega}(W_{ij}-M_{ij})^{2}$

where $\Omega=\\{(1,3),(1,4),(2,3),(3,1),(3,2),(4,1)\\}$.

Then for any
${\bm{W}}_{\mathrm{init}}\succeq\bm{0},~{}~{}\text{s.t.}~{}~{}{\bm{u}}_{1}^{\top}{\bm{W}}_{\mathrm{init}}{\bm{u}}_{1}>0$,

$\lim\limits_{\alpha\to
0}\lim\limits_{t\to+\infty}\phi(\alpha{\bm{W}}_{\mathrm{init}},t)={\bm{M}}_{\mathrm{rank}}.$

Moreover, we have

$\left\|{\bm{M}}_{\mathrm{rank}}\right\|_{\ast}=2R^{2}+2>4R=\left\|{\bm{M}}_{\mathrm{norm}}\right\|_{\ast}=\min\limits_{{\bm{W}}\geq\bm{0}}\left\|{\bm{W}}\right\|_{\ast}~{}~{}\text{s.t.}~{}~{}f({\bm{W}})=0.$

###### Proof.

We define ${\bm{W}}^{\mathrm{G}}_{1,\epsilon}(t),{\bm{W}}^{\mathrm{G}}_{1}(t)$
in the same way as in Definition 5.1, Theorem 5.6.

$\displaystyle{\bm{W}}^{\mathrm{G}}_{1,\epsilon}(t)$
$\displaystyle:=\phi\left(\epsilon{\bm{u}}_{1}{\bm{u}}_{1}^{\top},t\right),$
$\displaystyle{\bm{W}}^{\mathrm{G}}_{1}(t)$ $\displaystyle:=\lim_{\epsilon\to
0}{\bm{W}}^{\mathrm{G}}_{1,\epsilon}(\tfrac{1}{2\mu_{1}}\log\tfrac{1}{\epsilon}+t).$

Below we will show

  1. 1.

5.7 and 5.5 are satisfied.

  2. 2.

$\left\|{\bm{W}}^{\mathrm{G}}_{1}(t)\right\|_{\mathrm{F}}$ bounded for $t\geq
0$;

  3. 3.

$\lim_{t\to+\infty}{\bm{W}}^{\mathrm{G}}_{1}(t)={\bm{M}}_{\mathrm{rank}}$;

  4. 4.

${\bm{M}}_{\mathrm{norm}}=\operatorname*{arg\,min}\limits_{{\bm{W}}\geq\bm{0}}\left\|{\bm{W}}\right\|_{\ast}~{}~{}\text{s.t.}~{}~{}f({\bm{W}})=0.$

Thus Since ${\bm{M}}_{\mathrm{rank}}$ is a global minimizer of $f(\,\cdot\,)$,
applying Theorem 5.8 finishes the proof.

##### Proof for Item 1.

Let ${\bm{M}}_{0}:=\nabla f({\bm{0}})$, then

${\bm{M}}_{0}=\begin{bmatrix}0&0&1&R\\\ 0&0&R&0\\\ 1&R&0&0\\\
R&0&0&0\end{bmatrix}.$

Let ${\bm{A}}:=\left[\begin{smallmatrix}1&R\\\ R&0\\\
\end{smallmatrix}\right]$, then we have
$\lambda_{1}({\bm{A}})=\frac{1+\sqrt{1+R^{2}}}{2},\lambda_{2}({\bm{A}})=\frac{1-\sqrt{1+R^{2}}}{2}$,
thus $\lambda_{1}({\bm{A}})>|\lambda_{2}({\bm{A}})|>0>\lambda_{2}({\bm{A}})$.
As a result, $\lambda_{1}({\bm{A}})=\left\|{\bm{A}}\right\|_{2}$. Let
${\bm{v}}_{1}\in\mathbb{R}^{2}$ be the top eigenvector of ${\bm{A}}$. We claim
that ${\bm{u}}_{1}=\left[\begin{smallmatrix}{\bm{v}}_{1}\\\ {\bm{v}}_{1}\\\
\end{smallmatrix}\right]\in\mathbb{R}^{4}$ is the top eigenvector of $\nabla
f({\bm{0}})$. First by definition it is easy to check that
${\bm{M}}_{0}{\bm{u}}_{1}=\lambda_{1}({\bm{A}}){\bm{u}}_{1}$. Further noticing
that ${\bm{M}}_{0}^{2}=\left[\begin{smallmatrix}{\bm{A}}^{2}&{\bm{0}}\\\
{\bm{0}}&{\bm{A}}^{2}\\\ \end{smallmatrix}\right]$, we know
$\lambda^{2}_{i}({\bm{M}}_{0})\in\\{\lambda_{1}^{2}({\bm{A}}),\lambda_{2}^{2}({\bm{A}})\\}$,
$\forall i\in\\{1,2,3,4\\}$. That is,
$\lambda_{1}({\bm{M}}_{0})=\lambda_{1}({\bm{A}})$,
$\lambda_{2}({\bm{M}}_{0})=-\lambda_{2}({\bm{A}})$,
$\lambda_{3}({\bm{M}}_{0})=\lambda_{2}({\bm{A}})$, and
$\lambda_{4}({\bm{M}}_{0})=-\lambda_{1}({\bm{A}})$. Thus 5.5 is satisfied.
Also note that $f$ is quadratic, thus analytic, i,e., 5.7 is also satisfied.

##### Proof for Item 2.

Let $(x_{\epsilon}(t),y_{\epsilon}(t))\in\mathbb{R}^{2}$ be the gradient flow
of $g(x,y)=\frac{1}{2}(x^{2}-1)^{2}+(xy-R)^{2}$ starting from
$(x_{\epsilon}(0),y_{\epsilon}(0))=\sqrt{\epsilon}{\bm{v}}_{1}$.

$\begin{split}\frac{\textup{{d}}x(t)}{\textup{{d}}t}&=(1-x(t)^{2})x(t)-2y(t)(x(t)y(t)-R)\\\
\frac{\textup{{d}}y(t)}{\textup{{d}}t}&=-2x(t)(x(t)y(t)-R)\end{split}$ (18)

Let ${\bm{W}}_{\epsilon}(t)$ be the following matrix:

${\bm{W}}_{\epsilon}(t):=\begin{bmatrix}x_{\epsilon}(t)\\\ y_{\epsilon}(t)\\\
x_{\epsilon}(t)\\\
y_{\epsilon}(t)\end{bmatrix}\begin{bmatrix}x_{\epsilon}(t)&y_{\epsilon}(t)&x_{\epsilon}(t)&y_{\epsilon}(t)\end{bmatrix}.$

Then it is easy to verify that
${\bm{W}}_{\epsilon}(0)={\bm{W}}^{\mathrm{G}}_{1,\epsilon}(0)$ and
${\bm{W}}_{\epsilon}(t)$ satisfies (2). Thus by the existence and uniqueness
theorem, we have
${\bm{W}}_{\epsilon}(t)={\bm{W}}^{\mathrm{G}}_{1,\epsilon}(t)$ for all $t$.
Taking the limit $\epsilon\to 0$, we know that ${\bm{W}}^{\mathrm{G}}_{1}(t)$
can also be written in the following form:

${\bm{W}}^{\mathrm{G}}_{1}(t)=\begin{bmatrix}x(t)\\\ y(t)\\\ x(t)\\\
y(t)\end{bmatrix}\begin{bmatrix}x(t)&y(t)&x(t)&y(t)\end{bmatrix},$

and $(x_{\epsilon}(t),y_{\epsilon}(t))\in\mathbb{R}^{2}$ is a gradient flow of
$g(x,y)=\frac{1}{2}(x^{2}-1)^{2}+(xy-R)^{2}$.

Since $g(x(t),y(t))$ is non-increasing overtime, and
$\lim\limits_{t\to-\infty}g(x(-t),y(-t))=g(x(-\infty),y(-\infty))=g(0,0)=R^{2}+0.5$,
we know $\lvert x(t)y(t)\rvert\leq 3R$ for all $t$. So whenever
$y^{2}(t)-x^{2}(t)\geq 9R^{2}$, we have
$x^{2}(t)\leq\frac{9R^{2}}{y^{2}(t)}\leq\frac{9R^{2}}{y^{2}(t)-x^{2}(t)}\leq
1$. In this case,
$\frac{\textup{{d}}(y^{2}(t)-x^{2}(t))}{\textup{{d}}t}=2x^{2}(t)(x^{2}(t)-1)\leq
0$. Combining this with $y(-\infty)^{2}-x(-\infty)^{2}=0\leq 9R^{2}$, we have
$y^{2}(t)-x^{2}(t)\leq 9R^{2}$ for all $t$, which also implies that $y(t)$ is
bounded. Noticing that $9R^{2}\geq g(x(t),y(t))\geq(x^{2}(t)-1)^{2}$, we know
$x^{2}(t)$ is also bounded. Therefore, ${\bm{W}}_{1}^{G}(t)$ is bounded.

##### Proof for Item 3.

Note that $(x(\infty),y(\infty))$ is a stationary point of $g(x,y)$. It is
clear $g(x,y)$ only has $3$ stationary points — $(0,0)$, $(1,R)$ and
$(-1,-R)$. Thus $\overline{{}{\bm{W}}}_{1}$ can only be ${\bm{0}}$ or
${\bm{M}}_{\mathrm{rank}}$. However, since for all $t$,
$f({\bm{W}}^{\mathrm{G}}_{1}(t))<f({\bm{0}})$,
$\overline{{}{\bm{W}}}_{1}=\lim_{t\to\infty}{\bm{W}}^{\mathrm{G}}_{1}(t)$
cannot be ${\bm{0}}$. So $\overline{{}{\bm{W}}}_{1}$ must be
${\bm{M}}_{\mathrm{rank}}$.

##### Proof for Item 4.

Let $m_{ij}$ be $(i,j)$th element of ${\bm{M}}$. Suppose
${\bm{M}}\succeq{\bm{0}}$, we have

$\displaystyle({\bm{e}}_{1}-{\bm{e}}_{4})^{\top}{\bm{M}}({\bm{e}}_{1}-{\bm{e}}_{4})\geq
0$ $\displaystyle\Longrightarrow m_{11}+m_{44}\geq m_{14}+m_{41}=2R$
$\displaystyle({\bm{e}}_{2}-{\bm{e}}_{3})^{\top}{\bm{M}}({\bm{e}}_{2}-{\bm{e}}_{3})\geq
0$ $\displaystyle\Longrightarrow m_{22}+m_{33}\geq m_{23}+m_{32}=2R$

Thus
$4R=\min\limits_{{\bm{W}}\geq\bm{0}}\left\|{\bm{W}}\right\|_{*}~{}~{}\text{s.t.}~{}~{}f({\bm{W}})=0$,
where the equality is only attained at $m_{ii}=R,i=1,2,3,4$. Otherwise, either
$\begin{bmatrix}m_{11}&m_{14}\\\ m_{41}&m_{44}\end{bmatrix}$ or
$\begin{bmatrix}m_{22}&m_{23}\\\ m_{32}&m_{33}\end{bmatrix}$ will have
negative eigenvalues. Contradiction to that ${\bm{M}}\succeq{\bm{0}}$.

Below we will show the rest unknown off-diagonal entries must be $1$. Let
${\bm{V}}=\begin{bmatrix}1&-1&0&0\\\ 0&0&1&0\\\ 0&0&0&1\end{bmatrix}$, we have

${\bm{M}}\succeq{\bm{0}}\Longrightarrow{\bm{V}}{\bm{M}}{\bm{V}}^{\top}\succeq{\bm{0}}\Longrightarrow\begin{bmatrix}0&m_{13}-m_{23}&m_{14}-m_{24}\\\
m_{31}-m_{32}&R&R\\\ m_{41}-m_{42}&R&R\end{bmatrix}\succeq{\bm{0}},$

which implies $m_{13}=m_{23}$, $m_{14}=m_{24}$.

With the same argument for ${\bm{V}}=\begin{bmatrix}1&0&0&0\\\ 0&1&0&0\\\
0&0&1&-1\end{bmatrix}$, we have $m_{13}=m_{14}$, $m_{23}=m_{24}$. Also note
${\bm{M}}$ is symmetric and $m_{13}=1$, thus $m_{ij}=m_{ji}=1,\ \forall
i=1,2,j=3,4$. Thus
${\bm{M}}_{\mathrm{norm}}=\operatorname*{arg\,min}_{{\bm{W}}\succeq\bm{0}}\left\|{\bm{W}}\right\|_{\ast}~{}~{}\text{s.t.}~{}~{}f({\bm{W}})=0$,
which is unique. ∎

##  Appendix C Experiments

###  C.1 General Setup

The code is written in Julia [Bezanson et al., 2012] and PyTorch [Paszke et
al., 2019].

The ground-truth matrix ${\bm{W}}_{*}$ is low-rank by construction: we sample
a random orthogonal matrix ${\bm{U}}$, a diagonal matrix ${\bm{S}}$ with
Frobenius norm $\left\|{\bm{S}}\right\|_{\mathrm{F}}=1$ and set
${\bm{W}}_{*}={\bm{U}}{\bm{S}}{\bm{U}}^{\top}$. Each measurement ${\bm{X}}$ in
${\bm{X}}_{1},\dots,{\bm{X}}_{m}$ is generated by sampling two one-hot vectors
${\bm{u}}$ and ${\bm{v}}$ uniformly and setting
${\bm{X}}=\frac{1}{2}{\bm{u}}{\bm{v}}^{\top}+\frac{1}{2}{\bm{v}}{\bm{u}}^{\top}$.

In Figures 1, 2, 5, 4, 3 and 7, the ground truth matrix ${\bm{W}}_{*}$ has
shape $20\times 20$ and rank $3$, where
$\left\|{\bm{W}}^{*}\right\|_{\mathrm{F}}=20$,
$\lambda_{1}({\bm{W}}^{*})=17.41$, $\lambda_{2}({\bm{W}}^{*})=8.85$,
$\lambda_{3}({\bm{W}}^{*})=4.31$ and $\lambda_{1}(-\nabla
f({\bm{0}}))=6.23,\lambda_{2}(-\nabla f({\bm{0}}))=5.41$. $p=0.3$ is used for
generating measurements, except $p=0.25$ in Figure 3, i.e., each pair of
entries of ${\bm{W}}^{*}_{ij}$ and ${\bm{W}}^{*}_{ji}$ is observed with
probability $p$.

##### Gradient Descent.

Let $\tilde{\epsilon}>0$ be the Frobenius norm of the target random
initialization. For the depth-2 case, we sample $2$ orthogonal matrices
${\bm{V}}_{1},{\bm{V}}_{2}$ and a diagonal matrix ${\bm{D}}$ with Frobenius
norm $\tilde{\epsilon}$, and we set
${\bm{U}}={\bm{V}}_{1}{\bm{D}}^{1/2}{\bm{V}}_{2}^{\top}$; for the depth-$L$
case with $L\geq 3$, we sample $L$ orthogonal matrices
${\bm{V}}_{1},\dots,{\bm{V}}_{L}$ and a diagonal matrix ${\bm{D}}$ with
Frobenius norm $\tilde{\epsilon}$, and we set
${\bm{U}}_{i}:={\bm{V}}_{i}{\bm{D}}^{1/L}{\bm{V}}_{i+1}^{\top}$
(${\bm{V}}_{L+1}={\bm{V}}_{1}$). In this way, we can guarantee that the end-
to-end matrix ${\bm{W}}={\bm{U}}_{1}\cdots{\bm{U}}_{L}$ is symmetric and the
initialization is balanced for $L\geq 3$.

Depth ($L$) | Simulation method
---|---
2 | Constant LR, $\eta=10^{-3}$ for $10^{6}$ iterations
3 | Adaptive LR, $\eta=2\times 10^{-5}$ and $\varepsilon=10^{-4}$ for $10^{6}$ iterations
4 | Adaptive LR, $\eta=3\times 10^{-4}$ and $\varepsilon=10^{-3}$ for $10^{6}$ iterations
Table 1: Choice of hyperparameters for simulating gradient flow. For $L=2$,
gradient descent escapes saddles in $O(\log\frac{1}{{\epsilon}})$ time, where
${\epsilon}$ is the distance between the initialization and the saddle.

We discretize the time to simulate gradient flow. When $L>2$, gradient flow
stays around saddle points for most of the time, therefore we use full-batch
GD with adaptive learning rate $\tilde{\eta_{t}}$, inspired by RMSprop
[Tieleman and Hinton, 2012], for faster convergence:

$\displaystyle v_{t+1}$ $\displaystyle=\alpha
v_{t}+(1-\alpha)\left\|\nabla\mathcal{L}({\bm{\theta}}_{t})\right\|_{2}^{2},$
$\displaystyle\tilde{\eta}_{t}$
$\displaystyle=\frac{\eta}{\sqrt{\frac{v_{t+1}}{1-\alpha^{t+1}}}+\varepsilon},$
$\displaystyle{\bm{\theta}}_{t+1}$
$\displaystyle={\bm{\theta}}_{t}-\tilde{\eta}_{t}\nabla\mathcal{L}({\bm{\theta}}_{t}),$

where $\alpha=0.99$, $\eta$ is the (unadjusted) learning rate. The choices of
hyperparameters are summarized in Table 1. The continuous time for
${\bm{\theta}}_{t}$ is measured as $\sum_{i=0}^{t-1}\tilde{\eta}_{i}$.

##### GLRL.

In Figures 1, 2, 3 and 4, the GLRL’s trajectory is obtained by running
Algorithm 1 with $\epsilon=10^{-7}$ and $\eta=10^{-3}$. The stopping criterion
is that if the loop has been iterated for $10^{7}$ times.

###  C.2 Experimental Equivalence between GLRL and Gradient Descent

Here we provide experimental evidence supporting our theoretical claims about
the equivalence between GLRL and GF for both cases, $L=2$ and $L\geq 3$.

In Figure 1, we show the distance from every point on GF (simulated by GD)
from random initialization is close to the trajectory of GLRL. In Figure 2, we
first run GLRL and obtain the critical points
$\\{\overline{{}{\bm{W}}}_{r}\\}_{r=0}^{3}$ passed by GLRL. We also define the
distance of a matrix ${\bm{W}}$ to the critical points to be $\min_{0\leq
r\leq 3}\|{\bm{W}}-\overline{{}{\bm{W}}}_{r}\|_{\mathrm{F}}$.

###  C.3 How well does GLRL work?

We compare GLRL with gradient descent (with not-so-small initialization),
nuclear norm minimization and R1MP [Wang et al., 2014]. We use CVXPY [Diamond
and Boyd, 2016, Agrawal et al., 2018] for finding the nuclear norm solution.
The results are shown in Figure 3. GLRL can fully recover the ground truth,
while others have difficulty doing so.

Figure 3: GD with small initialization outperforms R1MP and minimal nuclear
norm solution on synthetic data with low-rank ground truth. Solid (dotted)
curves correspond to test (training) loss. Here the loss
$f({\bm{W}}):=\frac{1}{d^{2}}\|{\bm{W}}-{\bm{W}}^{*}\|_{\mathrm{F}}^{2}$ and
$f({\bm{0}})=1$. We run 10 random seeds for GD and plot them separately (most
of them overlap).

###  C.4 How does initialization affect the convergence rate to the rank-1
GLRL trajectory?

We use the general setting in Section C.1. In these experiments, we use the
constant learning rate $10^{-5}$ for $4\times 10^{7}$ iterations. The
reference matrix ${\bm{W}}_{\text{ref}}$ is obtained by running the first
stage of GLRL with $\left\|{\bm{W}}(0)\right\|_{\mathrm{F}}=10^{-48}$ and we
pick one matrix in the trajectory with
$\left\|{\bm{W}}_{\text{ref}}\right\|_{\mathrm{F}}$ about 0.6.

For every ${\epsilon}=10^{i},i\in\\{-1,-2,-3,-4,-5\\}$, we run both gradient
descent and the first phase of GLRL with
$\left\|{\bm{W}}(0)\right\|_{\mathrm{F}}={\epsilon}$. For gradient descent, we
use random initialization so $\left\|{\bm{W}}(0)\right\|_{\mathrm{F}}$ is full
rank w.p. 1. The distance of a trajectory to ${\bm{W}}_{\text{ref}}$ is
defined as $\min_{t\geq
0}\left\|{\bm{W}}(t)-{\bm{W}}_{\text{ref}}\right\|_{\mathrm{F}}$. In practice,
as we discretized time to simulate gradient flow, we check every $t$ during
simulation to compute the distance. As a result, the estimation might be
inaccurate when a trajectory is really close to ${\bm{W}}_{\text{ref}}$.

The result is shown at Figure 4. We observe that GLRL trajectories are closer
to the reference matrix ${\bm{W}}_{\text{ref}}$ by magnitudes. Thus the take
home message here is that GLRL is in general a more computational efficient
method to simulate the trajectory of GF (GD) with infinitesimal
initialization, as one can start GLRL with a much larger initialization, while
still maintaining high precision.

Figure 4:  Using ${\epsilon}{\bm{v}}_{1}{\bm{v}}_{1}^{\top}$ (denoted by “rank
$1$”) as initialization makes GD much closer to GLRL compared to using random
initialization (denoted by “rank $d$”), where ${\bm{v}}_{1}$ is the top
eigenvector of $-\nabla f({\bm{0}})$. We take a fixed reference matrix on the
trajectory of GLRL with constant norm and plot the distance of GD with each
initialization to it respectively..  Figure 5: Deep matrix factorization
encourages GF to find low rank solutions at a much practical initialization
scale, e.g. $10^{-3}$. Here the ground truth is rank-$3$. For each setting, we
run 5 different random seeds. The solid curves are the mean and the shaded
area indicates one standard deviation. We observe that performance of GD is
quite robust to its initialization. Note that for $L>2$, the shaded area with
initialization scale $10^{-7}$ is large, as the sudden decrement of loss
occurs at quite different continuous times for different random seeds in this
case.

###  C.5 Benefit of Depth: polynomial vs exponential dependence on
initialization

To verify the our theory in Section 6, we run gradient descent with different
depth and initialization. The results are shown in Figure 5. We can see that
as the initialization becomes smaller, the final solution gets closer to the
ground truth. However, a depth-2 model requires exponentially small
initialization, while deeper models require polynomial small initialization,
though it takes much longer to converge.

##  Appendix D The marginal value of being deeper

Theorem D.1 shows that the end-to-end dynamics (19) converges point-wise while
$L\to\infty$ if the product of learning rate and depth, $\eta L$, is fixed as
constant. Interestingly, (19) also allows us to simulate the dynamics of
${\bm{W}}(t)$ for all depths $L$ while the computation time is independent of
$L$. In Figure 6, we compare the effect of depth while fixing the
initialization and $\eta L$. We can see that deeper models converge faster.
The difference between $L=1,2$, and $4$ is large, while difference among
$L\geq 16$ is marginal.

Figure 6:  The marginal value of being deeper. The trajectory of GD converges
when depth goes to infinity. Solid (dotted) curves correspond to test (train)
loss. The $x$-axis stands for the normalized continuous time $t$ (multiplied
by $L$).

######  Theorem D.1.

Suppose ${\bm{W}}=\tilde{\bm{U}}\tilde{\bm{\Sigma}}\tilde{\bm{V}}^{\top}$ is
the SVD decomposition of ${\bm{W}}$ where
$\tilde{\bm{\Sigma}}=\mathrm{diag}(\sigma_{1},\dots,\sigma_{d})$. The dynamics
of $L$-layer linear net is the following, $\circ$ denotes the entry-wise
multiplication:

$\frac{\textup{{d}}{\bm{W}}}{\textup{{d}}t}=-L\tilde{\bm{U}}\left(\left(\tilde{\bm{U}}^{\top}\nabla
f({\bm{W}})\tilde{\bm{V}}\right)\circ{\bm{K}}^{(L)}\right)\tilde{\bm{V}}^{\top},$
(19)

where $K^{(L)}_{i,i}=\sigma_{i}^{2-2/L}$,
$K^{(L)}_{i,j}=\frac{\sigma_{i}^{2}-\sigma_{j}^{2}}{L\sigma_{i}^{2/L}-L\sigma_{j}^{2/L}}$
for $i\neq j$.

###### Proof.

We start from (11):

$\displaystyle\frac{\textup{{d}}{\bm{W}}}{\textup{{d}}t}$
$\displaystyle=-\sum_{l=0}^{L-1}({\bm{W}}{\bm{W}}^{\top})^{\frac{l}{L}}\nabla
f({\bm{W}})({\bm{W}}^{\top}{\bm{W}})^{\frac{L-1-l}{L}}$
$\displaystyle=-\sum_{l=0}^{L-1}\tilde{\bm{U}}\tilde{\bm{\Sigma}}^{\frac{2l}{L}}\tilde{\bm{U}}^{\top}\nabla
f({\bm{W}})\tilde{\bm{V}}\tilde{\bm{\Sigma}}^{\frac{2(L-1-l)}{L}}\tilde{\bm{V}}$
$\displaystyle=-L\tilde{\bm{U}}\left[L^{-1}\sum_{l=0}^{L-1}\tilde{\bm{\Sigma}}^{\frac{2l}{L}}(\tilde{\bm{U}}^{\top}\nabla
f({\bm{W}})\tilde{\bm{V}})\tilde{\bm{\Sigma}}^{\frac{2(L-1-l)}{L}}\right]\tilde{\bm{V}}.$

Note that $\tilde{\bm{\Sigma}}$ is diagonal, so

$\displaystyle\tilde{\bm{\Sigma}}^{\frac{2l}{L}}(\tilde{\bm{U}}^{\top}\nabla
f({\bm{W}})\tilde{\bm{V}})\tilde{\bm{\Sigma}}^{\frac{2(L-1-l)}{L}}=(\tilde{\bm{U}}^{\top}\nabla
f({\bm{W}})\tilde{\bm{V}})\circ{\bm{H}}^{(l)},$

where
${\bm{H}}^{(l)}_{i,j}=\sigma_{i}^{\frac{2l}{L}}\sigma_{j}^{\frac{2(L-1-l)}{L}}$.
Therefore

$\displaystyle
L^{-1}\sum_{l=0}^{L-1}\tilde{\bm{\Sigma}}^{\frac{2l}{L}}(\tilde{\bm{U}}^{\top}\nabla
f({\bm{W}})\tilde{\bm{V}})\tilde{\bm{\Sigma}}^{\frac{2(L-1-l)}{L}}$
$\displaystyle=L^{-1}\sum_{l=0}^{L-1}(\tilde{\bm{U}}^{\top}\nabla
f({\bm{W}})\tilde{\bm{V}})\circ{\bm{H}}^{(l)}$
$\displaystyle=(\tilde{\bm{U}}^{\top}\nabla
f({\bm{W}})\tilde{\bm{V}})\circ{\bm{K}}^{(L)},$

where ${\bm{K}}^{(L)}=L^{-1}\sum_{l=0}^{L-1}{\bm{H}}^{(l)}$. Hence

$\displaystyle\frac{\textup{{d}}{\bm{W}}}{\textup{{d}}t}$
$\displaystyle=-L\tilde{\bm{U}}\left[(\tilde{\bm{U}}^{\top}\nabla
f({\bm{W}})\tilde{\bm{V}})\circ{\bm{K}}^{(L)}\right]\tilde{\bm{V}}.$

The entries of ${\bm{K}}^{(L)}$ can be directly calculated by

$K^{(L)}_{i,j}=L^{-1}\sum_{l=0}^{L-1}\sigma_{i}^{\frac{2l}{L}}\sigma_{j}^{\frac{2(L-1-l)}{L}}=\begin{cases}\sigma_{i}^{2-2/L},&i=j,\\\
\frac{\sigma_{i}^{2}-\sigma_{j}^{2}}{L\sigma_{i}^{2/L}-L\sigma_{j}^{2/L}},&i\neq
j.\end{cases}$

∎

######  Corollary D.2.

As $L\to\infty$, ${\bm{K}}^{(L)}$ converges to ${\bm{K}}^{*}$, where
$K^{*}_{i,i}=\sigma_{i}^{2}$,
$K^{*}_{i,j}=\frac{\sigma_{i}^{2}-\sigma_{j}^{2}}{\ln\sigma_{i}^{2}-\ln\sigma_{j}^{2}}$
for $i\neq j$.

##### Experiment details

We follow the general setting in Section C.1. The ground truth ${\bm{W}}^{*}$
is different but is generated in the same manner and has the same shape of
$20\times 20$ and $p=0.3$ is used for observation generation. We directly
apply (19), in which we compute $\tilde{\bm{V}}$ and $\tilde{\bm{U}}$ through
SVD, to simulate the trajectory together with a constant learning rate of
$\frac{10^{-3}}{L}$ for depth $L$. ${\bm{W}}(0)$ is sampled from
$10^{-3}\times\mathcal{N}(0,{\bm{I}}_{d})$.

##  Appendix E Proofs for Dynamical System

In this section, we prove Theorem 5.3 in Section 5.1. In Section E.1, we show
how to reduce Theorem 5.3 to the case where ${\bm{J}}({\bm{0}})$ is exactly a
diagonal matrix, then we prove this diagonal case in Section E.2. Finally, in
Section E.3, we discuss how to extend it to the case where
${\bm{J}}({\bm{0}})$ is non-diagonalizable.

###  E.1 Reduction to the Diagonal Case

######  Theorem E.1.

If ${\bm{J}}({\bm{0}})=\mathrm{diag}(\tilde{\mu}_{1},\dots,\tilde{\mu}_{d})$
is diagonal, then the statement in Theorem 5.3 holds.

###### Proof for Theorem 5.3.

We show how to prove Theorem 5.3 based on Theorem E.1. Let
$\frac{\textup{{d}}{\bm{\theta}}}{\textup{{d}}t}={\bm{g}}({\bm{\theta}})$ be
the dynamical system in Theorem 5.3. Let
${\bm{J}}({\bm{0}})=\tilde{{\bm{V}}}\tilde{{\bm{D}}}\tilde{{\bm{V}}}^{-1}$ be
the eigendecomposition, where $\tilde{{\bm{V}}}$ is an invertible matrix and
$\tilde{{\bm{D}}}=\mathrm{diag}(\tilde{\mu}_{1},\dots,\tilde{\mu}_{d})$. Now
we define the following new dynamics by changing the basis:

$\hat{{\bm{\theta}}}(t)=\tilde{{\bm{V}}}^{-1}{\bm{\theta}}(t).$

Then
$\frac{\textup{{d}}\hat{{\bm{\theta}}}(t)}{\textup{{d}}t}=\hat{{\bm{g}}}(\hat{{\bm{\theta}}})$
for
$\hat{{\bm{g}}}(\hat{{\bm{\theta}}}):=\tilde{{\bm{V}}}^{-1}{\bm{g}}(\tilde{{\bm{V}}}\hat{{\bm{\theta}}})$,
and the associated Jacobian matrix is
$\hat{{\bm{J}}}(\hat{{\bm{\theta}}}):=\tilde{{\bm{V}}}^{-1}{\bm{J}}(\tilde{\bm{V}}\hat{{\bm{\theta}}})\tilde{{\bm{V}}}$,
and thus
$\hat{{\bm{J}}}({\bm{0}})=\mathrm{diag}(\tilde{\mu}_{1},\dots,\tilde{\mu}_{d})$.

Now we apply Theorem E.1 to $\hat{{\bm{\theta}}}(t)$. Then
$\hat{{\bm{z}}}_{\alpha}(t):=\tilde{\bm{V}}^{-1}{\bm{z}}_{\alpha}(t)$
converges to the limit $\hat{{\bm{z}}}(t):=\lim\limits_{\alpha\to
0}\hat{{\bm{z}}}_{\alpha}(t)$. This shows that the limit
${\bm{z}}(t)=\tilde{{\bm{V}}}\hat{{\bm{z}}}(t)$ exists in Theorem 5.3. We can
also verify that ${\bm{z}}(t)$ is a solution of (6).

Given ${\bm{\delta}}_{\alpha}$ converging to ${\bm{0}}$ with positive
alignment with $\tilde{{\bm{u}}}_{1}$ as $\alpha\to 0$, we can define
$\hat{{\bm{\delta}}}_{\alpha}:=\tilde{\bm{V}}^{-1}{\bm{\delta}}_{\alpha}$,
then $\hat{{\bm{\delta}}}_{\alpha}$ converges to ${\bm{0}}$ with positive
alignment with ${\bm{e}}_{1}$, where ${\bm{e}}_{1}$ is the first vector in the
standard basis and is also the top eigenvector of $\hat{{\bm{J}}}({\bm{0}})$.
Therefore, for every $t\in(-\infty,+\infty)$, there is a constant $C>0$ such
that

$\left\|\tilde{\bm{V}}^{-1}\phi\left({\bm{\delta}}_{\alpha},t+\frac{1}{\tilde{\mu}_{1}}\log\frac{1}{\left<{\bm{\delta}}_{\alpha},\tilde{{\bm{u}}}_{1}\right>}\right)-\hat{{\bm{z}}}(t)\right\|_{2}\leq
C\cdot\|\hat{{\bm{\delta}}}_{\alpha}\|_{2}^{\frac{\tilde{\gamma}}{\tilde{\mu}_{1}+\tilde{\gamma}}}$
(20)

for every sufficiently small $\alpha$. As $\tilde{\bm{V}}$ are invertible,
this directly implies (7). ∎

###  E.2 Proof for the Diagonal Case

Now we only need to prove Theorem E.1. Let ${\bm{e}}_{1},\dots,{\bm{e}}_{d}$
be the standard basis. Then
$\tilde{{\bm{u}}}_{1}=\tilde{{\bm{v}}}_{1}={\bm{e}}_{1}$ in this diagonal
case. We only use ${\bm{e}}_{1}$ to stand for $\tilde{{\bm{u}}}_{1}$ and
$\tilde{{\bm{v}}}_{1}$ in the rest of our analysis.

Let $R>0$. Since ${\bm{g}}({\bm{\theta}})$ is $\mathcal{C}^{2}$-smooth, there
exists $\beta>0$ such that

$\|{\bm{J}}({\bm{\theta}})-{\bm{J}}({\bm{\theta}}+{\bm{h}})\|_{2}\leq\beta\|{\bm{h}}\|_{2}$
(21)

for all $\|{\bm{\theta}}\|_{2},\|{\bm{\theta}}+{\bm{h}}\|_{2}\leq R$. Then the
following can be proved by integration:

$\displaystyle{\bm{g}}({\bm{\theta}}+{\bm{h}})-{\bm{g}}({\bm{\theta}})$
$\displaystyle=\left(\int_{0}^{1}{\bm{J}}({\bm{\theta}}+\xi{\bm{h}})\textup{{d}}\xi\right){\bm{h}},$
(22)
$\displaystyle\|{\bm{g}}({\bm{\theta}}+{\bm{h}})-{\bm{g}}({\bm{\theta}})-{\bm{J}}({\bm{\theta}}){\bm{h}}\|_{2}$
$\displaystyle\leq\beta\|{\bm{h}}\|_{2}^{2}.$ (23)

By (23), we also have

$\|{\bm{g}}({\bm{\theta}})-{\bm{J}}({\bm{0}}){\bm{\theta}}\|_{2}=\|{\bm{g}}({\bm{\theta}})-{\bm{g}}({\bm{0}})-{\bm{J}}({\bm{0}}){\bm{\theta}}\|_{2}\leq\beta\|{\bm{\theta}}\|_{2}^{2}.$
(24)

Let $\kappa:=\beta/\tilde{\mu}_{1}$. We assume WLOG that $R\leq 1/\kappa$. Let
$F(x)=\log x-\log(1+\kappa x)$. It is easy to see that
$F^{\prime}(x)=\frac{1}{x+\kappa x^{2}}$ and $F(x)$ is an increasing function
with range $(-\infty,\log(1/\kappa))$. We use $F^{-1}(y)$ to denote the
inverse function of $F(x)$. Define
$T_{\alpha}(r):=\frac{1}{\tilde{\mu}_{1}}\left(F(r)-F(\alpha)\right)=\frac{1}{\tilde{\mu}_{1}}\left(\log\frac{r}{\alpha}-\log\frac{1+\kappa
r}{1+\kappa\alpha}\right)$.

Our proof only relies on the following properties of ${\bm{J}}({\bm{0}})$
(besides that $\tilde{\mu}_{1},{\bm{e}}_{1}$ are the top eigenvalue and
eigenvector of ${\bm{J}}({\bm{0}})$):

######  Lemma E.2.

For
${\bm{J}}({\bm{0}}):=\mathrm{diag}(\tilde{\mu}_{1},\dots,\tilde{\mu}_{d})$, we
have

  1. 1.

For any ${\bm{h}}\in\mathbb{R}^{d}$,
${\bm{h}}^{\top}{\bm{J}}({\bm{0}}){\bm{h}}\leq\tilde{\mu}_{1}\|{\bm{h}}\|_{2}^{2}$;

  2. 2.

For any $t\geq 0$,
$\left\|e^{t{\bm{J}}({\bm{0}})}-e^{\tilde{\mu}_{1}t}{\bm{e}}_{1}{\bm{e}}_{1}^{\top}\right\|_{2}=e^{\tilde{\mu}_{2}t}$.

###### Proof.

For Item 1,
${\bm{h}}^{\top}{\bm{J}}({\bm{0}}){\bm{h}}=\sum_{i=1}^{d}\tilde{\mu}_{i}h_{i}^{2}\leq\tilde{\mu}_{1}\|{\bm{h}}\|_{2}^{2}$.
For Item 2,
$\left\|e^{t{\bm{J}}({\bm{0}})}-e^{\tilde{\mu}_{1}t}{\bm{e}}_{1}{\bm{e}}_{1}^{\top}\right\|_{2}=\left\|\mathrm{diag}(0,e^{\tilde{\mu}_{2}t},\dots,e^{\tilde{\mu}_{d}t})\right\|_{2}=e^{\tilde{\mu}_{2}t}$.
∎

######  Lemma E.3.

For ${\bm{\theta}}(t)=\phi({\bm{\theta}}_{0},t)$ with
$\|{\bm{\theta}}_{0}\|_{2}\leq\alpha$ and $t\leq T_{\alpha}(r)$,

$\|{\bm{\theta}}(t)\|_{2}\leq\frac{1+\kappa r}{1+\kappa\alpha}\alpha\cdot
e^{\tilde{\mu}_{1}t}\leq r.$

###### Proof.

By (24) and Lemma E.2, we have

$\frac{1}{2}\frac{\textup{{d}}\|{\bm{\theta}}(t)\|_{2}^{2}}{\textup{{d}}t}=\left<{\bm{\theta}}(t),{\bm{g}}({\bm{\theta}}(t))\right>\leq\left<{\bm{\theta}}(t),{\bm{J}}({\bm{0}}){\bm{\theta}}(t)\right>+\beta\|{\bm{\theta}}(t)\|_{2}^{3}\leq\tilde{\mu}_{1}\|{\bm{\theta}}(t)\|_{2}^{2}+\beta\|{\bm{\theta}}(t)\|_{2}^{3}.$

This implies
$\frac{\textup{{d}}\|{\bm{\theta}}(t)\|_{2}}{\textup{{d}}t}\leq\tilde{\mu}_{1}(\|{\bm{\theta}}(t)\|_{2}+\kappa\|{\bm{\theta}}(t)\|_{2}^{2})$.
Since $F^{\prime}(x)=\frac{1}{x+\kappa x^{2}}$, we further have

$\frac{\textup{{d}}}{\textup{{d}}t}F(\|{\bm{\theta}}(t)\|_{2})\leq\tilde{\mu}_{1}.$

So $F(\|{\bm{\theta}}(t)\|_{2})\leq F(\alpha)+\tilde{\mu}_{1}t$. By definition
of $T_{\alpha}(r)$, we then know that $\|{\bm{\theta}}(t)\|_{2}\leq r$ for all
$t\leq T_{\alpha}(r)$. So

$\log\|{\bm{\theta}}(t)\|_{2}\leq F(\|{\bm{\theta}}(t)\|_{2})+\log(1+\kappa
r)\leq F(\alpha)+\tilde{\mu}_{1}t+\log(1+\kappa r).$

Expending $F(\alpha)$ proves the lemma. ∎

######  Lemma E.4.

For ${\bm{\theta}}(t)=\phi({\bm{\theta}}_{0},t)$ with
$\|{\bm{\theta}}_{0}\|_{2}\leq\alpha$ and $t\leq T_{\alpha}(r)$, we have

${\bm{\theta}}(t)=e^{t{\bm{J}}({\bm{0}})}{\bm{\theta}}_{0}+O(r^{2}).$

###### Proof.

Let $\hat{{\bm{\theta}}}(t)=e^{t{\bm{J}}({\bm{0}})}{\bm{\theta}}_{0}$. Then we
have

$\displaystyle\frac{1}{2}\frac{\textup{{d}}}{\textup{{d}}t}\|{\bm{\theta}}(t)-\hat{{\bm{\theta}}}(t)\|_{2}^{2}$
$\displaystyle\leq\left<{\bm{g}}({\bm{\theta}}(t))-{\bm{J}}({\bm{0}})\hat{{\bm{\theta}}}(t),{\bm{\theta}}(t)-\hat{{\bm{\theta}}}(t)\right>$
$\displaystyle=\left<{\bm{g}}({\bm{\theta}}(t))-{\bm{J}}({\bm{0}}){\bm{\theta}}(t),{\bm{\theta}}(t)-\hat{{\bm{\theta}}}(t)\right>+({\bm{\theta}}(t)-\hat{{\bm{\theta}}}(t))^{\top}{\bm{J}}({\bm{0}})({\bm{\theta}}(t)-\hat{{\bm{\theta}}}(t))$
$\displaystyle\leq\|{\bm{g}}({\bm{\theta}}(t))-{\bm{J}}({\bm{0}}){\bm{\theta}}(t)\|_{2}\cdot\|{\bm{\theta}}(t)-\hat{{\bm{\theta}}}(t)\|_{2}+\tilde{\mu}_{1}\|{\bm{\theta}}(t)-\hat{{\bm{\theta}}}(t)\|_{2}^{2},$

where the last inequality is due to Lemma E.2. By (24) and Lemma E.3, we have

$\|{\bm{g}}({\bm{\theta}}(t))-{\bm{J}}({\bm{0}}){\bm{\theta}}(t)\|_{2}\leq\beta\|{\bm{\theta}}(t)\|_{2}^{2}\leq\beta\left(\frac{1+\kappa
r}{1+\kappa\alpha}\alpha\right)^{2}\cdot e^{2\tilde{\mu}_{1}t}.$

So we have
$\frac{\textup{{d}}}{\textup{{d}}t}\|{\bm{\theta}}(t)-\hat{{\bm{\theta}}}(t)\|_{2}\leq\beta\left(\frac{1+\kappa
r}{1+\kappa\alpha}\alpha\right)^{2}\cdot
e^{2\tilde{\mu}_{1}t}+\tilde{\mu}_{1}\|{\bm{\theta}}(t)-\hat{{\bm{\theta}}}(t)\|_{2}$.
By Grönwall’s inequality,

$\|{\bm{\theta}}(t)-\hat{{\bm{\theta}}}(t)\|_{2}\leq\int_{0}^{t}\beta\left(\frac{1+\kappa
r}{1+\kappa\alpha}\alpha\right)^{2}\cdot
e^{2\tilde{\mu}_{1}\tau}e^{\tilde{\mu}_{1}(t-\tau)}d\tau.$

Evaluating the integral gives

$\|{\bm{\theta}}(t)-\hat{{\bm{\theta}}}(t)\|_{2}\leq\beta\left(\frac{1+\kappa
r}{1+\kappa\alpha}\alpha\right)^{2}e^{\tilde{\mu}_{1}t}\cdot\frac{e^{\tilde{\mu}_{1}t}-1}{\tilde{\mu}_{1}}\leq\kappa\left(\frac{1+\kappa
r}{1+\kappa\alpha}\alpha\cdot e^{\tilde{\mu}_{1}t}\right)^{2}\leq\kappa
r^{2},$

which proves the lemma. ∎

######  Lemma E.5.

Let
${\bm{\theta}}(t)=\phi({\bm{\theta}}_{0},t),\hat{{\bm{\theta}}}(t)=\phi(\hat{{\bm{\theta}}}_{0},t)$.
If
$\max\\{\|{\bm{\theta}}_{0}\|_{2},\|\hat{{\bm{\theta}}}_{0}\|_{2}\\}\leq\alpha$,
then for $t\leq T_{\alpha}(r)$,

$\|{\bm{\theta}}(t)-\hat{{\bm{\theta}}}(t)\|_{2}\leq
e^{\tilde{\mu}_{1}t+\kappa
r}\|{\bm{\theta}}_{0}-\hat{{\bm{\theta}}}_{0}\|_{2}.$

###### Proof.

For $t\leq T_{\alpha}(r)$, by (22),

$\displaystyle\frac{1}{2}\frac{\textup{{d}}}{\textup{{d}}t}\|{\bm{\theta}}(t)-\hat{{\bm{\theta}}}(t)\|_{2}^{2}$
$\displaystyle=\left<{\bm{g}}({\bm{\theta}}(t))-{\bm{g}}(\hat{{\bm{\theta}}}(t)),{\bm{\theta}}(t)-\hat{{\bm{\theta}}}(t)\right>$
$\displaystyle=({\bm{\theta}}(t)-\hat{{\bm{\theta}}}(t))^{\top}\left(\int_{0}^{1}{\bm{J}}({\bm{\theta}}_{\xi}(t))\textup{{d}}\xi\right)({\bm{\theta}}(t)-\hat{{\bm{\theta}}}(t)),$

where
${\bm{\theta}}_{\xi}(t):=\xi{\bm{\theta}}(t)+(1-\xi)\hat{{\bm{\theta}}}(t)$.
By Lemma E.3,
$\max\\{\|{\bm{\theta}}(t)\|_{2},\|\hat{{\bm{\theta}}}(t)\|_{2}\\}\leq\frac{1+\kappa
r}{1+\kappa\alpha}\alpha\cdot e^{\tilde{\mu}_{1}t}$ for all $t\leq
T_{\alpha}(r)$. So $\|{\bm{\theta}}_{\xi}(t)\|_{2}\leq\frac{1+\kappa
r}{1+\kappa\alpha}\alpha\cdot e^{\tilde{\mu}_{1}t}$. Combining these with (21)
and Lemma E.2, we have

${\bm{h}}^{\top}{\bm{J}}({\bm{\theta}}_{\xi}(t)){\bm{h}}={\bm{h}}^{\top}{\bm{J}}({\bm{0}}){\bm{h}}+{\bm{h}}^{\top}({\bm{J}}({\bm{\theta}}_{\xi}(t))-{\bm{J}}({\bm{0}})){\bm{h}}\leq\left(\tilde{\mu}_{1}+\beta\cdot\frac{1+\kappa
r}{1+\kappa\alpha}\alpha\cdot
e^{\tilde{\mu}_{1}t}\right)\|{\bm{h}}\|_{2}^{2},$

for all ${\bm{h}}\in\mathbb{R}^{d}$. Thus,
$\frac{\textup{{d}}}{\textup{{d}}t}\|{\bm{\theta}}(t)-\hat{{\bm{\theta}}}(t)\|_{2}\leq\left(\tilde{\mu}_{1}+\beta\cdot\frac{1+\kappa
r}{1+\kappa\alpha}\alpha\cdot
e^{\tilde{\mu}_{1}t}\right)\|{\bm{\theta}}(t)-\hat{{\bm{\theta}}}(t)\|_{2}$.
This implies

$\displaystyle\log\frac{\|{\bm{\theta}}(t)-\hat{{\bm{\theta}}}(t)\|_{2}}{\|{\bm{\theta}}(0)-\hat{{\bm{\theta}}}(0)\|_{2}}$
$\displaystyle\leq\int_{0}^{t}\left(\tilde{\mu}_{1}+\beta\cdot\frac{1+\kappa
r}{1+\kappa\alpha}\alpha\cdot e^{\tilde{\mu}_{1}\tau}\right)d\tau$
$\displaystyle\leq\tilde{\mu}_{1}t+\kappa\cdot\frac{1+\kappa
r}{1+\kappa\alpha}\alpha e^{\tilde{\mu}_{1}t}$
$\displaystyle\leq\tilde{\mu}_{1}t+\kappa r.$

Therefore, $\|{\bm{\theta}}(t)-\hat{{\bm{\theta}}}(t)\|_{2}\leq
e^{\tilde{\mu}_{1}t+\kappa r}\|{\bm{\theta}}(0)-\hat{{\bm{\theta}}}(0)\|_{2}$.
∎

######  Lemma E.6.

For every $t\in(-\infty,+\infty)$, ${\bm{z}}(t)$ exists and
${\bm{z}}_{\alpha}(t)$ converges to ${\bm{z}}(t)$ in the following rate:

$\left\|{\bm{z}}_{\alpha}(t)-{\bm{z}}(t)\right\|_{2}=O(\alpha),$

where $O$ hides constants depending on $g({\bm{\theta}})$ and $t$.

###### Proof.

We prove the lemma in the cases of $t\in(-\infty,F(R)/\tilde{\mu}_{1}]$ and
$t>F(R)/\tilde{\mu}_{1}$ respectively.

##### Case 1.

Fix $t\in(-\infty,F(R)/\tilde{\mu}_{1}]$. Let $\tilde{\alpha}$ be the unique
number such that $\frac{\tilde{\alpha}}{1+\kappa\tilde{\alpha}}=\alpha$ (i.e.,
$F(\tilde{\alpha})=\log\alpha$). Let $\alpha^{\prime}$ be an arbitrary number
less than $\alpha$. Let
$t_{0}:=\frac{1}{\tilde{\mu}_{1}}\log\frac{\alpha}{\alpha^{\prime}}$. Then
$t_{0}=\frac{1}{\tilde{\mu}_{1}}\left(F(\tilde{\alpha})-\log\alpha^{\prime}\right)\leq
T_{\alpha^{\prime}}(\tilde{\alpha})$. By Lemma E.4, we have

$\left\|\phi\left(\alpha^{\prime}{\bm{e}}_{1},t_{0}\right)-\alpha{\bm{e}}_{1}\right\|_{2}=\left\|\phi\left(\alpha^{\prime}{\bm{e}}_{1},t_{0}\right)-e^{t_{0}{\bm{J}}({\bm{0}})}\alpha^{\prime}{\bm{e}}_{1}\right\|_{2}=O(\tilde{\alpha}^{2}).$

Let $r:=F^{-1}(\tilde{\mu}_{1}t)\leq R$. Then
$t+\frac{1}{\tilde{\mu}_{1}}\log\frac{1}{\alpha}=T_{\tilde{\alpha}}(r)$ if
$\tilde{\alpha}<r$.

By Lemma E.3,
$\left\|\phi\left(\alpha^{\prime}{\bm{e}}_{1},t_{0}\right)\right\|_{2}\leq\tilde{\alpha}$.
Also,
$\|\alpha{\bm{e}}_{1}\|_{2}=\frac{\tilde{\alpha}}{1+\kappa\tilde{\alpha}}\leq\tilde{\alpha}$.
By Lemma E.5,

$\displaystyle\left\|{\bm{z}}_{\alpha}(t)-{\bm{z}}_{\alpha^{\prime}}(t)\right\|_{2}$
$\displaystyle=\left\|\phi\left(\alpha^{\prime}{\bm{e}}_{1},t+\frac{1}{\tilde{\mu}_{1}}\log\frac{1}{\alpha^{\prime}}\right)-\phi\left(\alpha{\bm{e}}_{1},t+\frac{1}{\tilde{\mu}_{1}}\log\frac{1}{\alpha}\right)\right\|_{2}$
$\displaystyle=\left\|\phi\left(\phi(\alpha^{\prime}{\bm{e}}_{1},t_{0}),t+\frac{1}{\tilde{\mu}_{1}}\log\frac{1}{\alpha}\right)-\phi\left(\alpha{\bm{e}}_{1},t+\frac{1}{\tilde{\mu}_{1}}\log\frac{1}{\alpha}\right)\right\|_{2}$
$\displaystyle\leq O(\tilde{\alpha}^{2}\cdot
e^{\tilde{\mu}_{1}(t+\frac{1}{\tilde{\mu}_{1}}\log\frac{1}{\alpha})+\kappa
r})$ $\displaystyle\leq O\left(\frac{\tilde{\alpha}^{2}}{\alpha}\right).$

For $\alpha$ small enough, we have $\tilde{\alpha}=O(\alpha)$, so for any
$\alpha^{\prime}\in(0,\alpha)$,

$\left\|{\bm{z}}_{\alpha}(t)-{\bm{z}}_{\alpha^{\prime}}(t)\right\|_{2}=O(\alpha).$

This implies that $\\{{\bm{z}}_{\alpha}(t)\\}$ satisfies Cauchy’s criterion
for every $t$, and thus the limit ${\bm{z}}(t)$ exists for $t\leq
F(R)/\tilde{\mu}_{1}$. The convergence rate can be deduced by taking limits
for $\alpha^{\prime}\to 0$ on both sides.

##### Case 2.

For $t=F(R)/\tilde{\mu}_{1}+\tau$ with $\tau>0$, $\phi({\bm{\theta}},\tau)$ is
locally Lipschitz with respect to ${\bm{\theta}}$. So

$\displaystyle\left\|{\bm{z}}_{\alpha}(t)-{\bm{z}}_{\alpha^{\prime}}(t)\right\|_{2}$
$\displaystyle=\left\|\phi({\bm{z}}_{\alpha}(F(R)/\tilde{\mu}_{1}),\tau)-\phi({\bm{z}}_{\alpha^{\prime}}(F(R)/\tilde{\mu}_{1}),\tau)\right\|_{2}$
$\displaystyle=O(\left\|{\bm{z}}_{\alpha}(F(R)/\tilde{\mu}_{1})-{\bm{z}}_{\alpha^{\prime}}(F(R)/\tilde{\mu}_{1})\right\|_{2})$
$\displaystyle=O(\alpha),$

which proves the lemma for $t>F(R)/\tilde{\mu}_{1}$. ∎

###### Proof for Theorem E.1.

The existence of ${\bm{z}}(t):=\lim_{\alpha\to
0}{\bm{z}}_{\alpha}(t)=\lim_{\alpha\to
0}\phi\left(\alpha{\bm{e}}_{1},t+\frac{1}{\tilde{\mu}_{1}}\log\frac{1}{\alpha}\right)$
has already been proved in Lemma E.6, where we show
$\left\|{\bm{z}}_{\alpha}(t)-{\bm{z}}(t)\right\|_{2}=O(\alpha)$.

By the continuity of $\phi(\,\cdot\,,t)$ for every $t\in\mathbb{R}$, we have

${\bm{z}}(t)=\lim_{\alpha\to
0}\phi\left(\alpha\tilde{{\bm{v}}}_{1},t+\frac{1}{\tilde{\mu}_{1}}\log\frac{1}{\alpha}\right)=\phi\left(\lim_{\alpha\to
0}\phi\left(\alpha\tilde{{\bm{v}}}_{1},\frac{1}{\tilde{\mu}_{1}}\log\frac{1}{\alpha}\right),t\right)=\phi\left({\bm{z}}(0),t\right).$

Now it is only left to prove (7). WLOG we can assume that
$\|{\bm{\delta}}_{\alpha}\|_{2}$ is decreasing and
$\frac{\alpha}{2}\leq\|{\bm{\delta}}_{\alpha}\|_{2}\leq\alpha$ (otherwise we
can do reparameterization). Then our goal becomes proving

$\left\|{\bm{\theta}}_{\alpha}(t)-{\bm{z}}(t)\right\|_{2}=O\left(\alpha^{\frac{\tilde{\gamma}}{\tilde{\mu}_{1}+\tilde{\gamma}}}\right).$
(25)

where
${\bm{\theta}}_{\alpha}(t):=\phi\left({\bm{\delta}}_{\alpha},t+\frac{1}{\tilde{\mu}_{1}}\log\frac{1}{\left<{\bm{\delta}}_{\alpha},{\bm{e}}_{1}\right>}\right)$.
We prove (25) in the cases of $t\in(-\infty,F(R)/\tilde{\mu}_{1}]$ and
$t>F(R)/\tilde{\mu}_{1}$ respectively.

##### Case 1.

Fix $t\in(-\infty,(F(R)+\log q)/\tilde{\mu}_{1}]$. Let
$\tilde{\alpha}_{1}=\alpha^{\frac{\tilde{\gamma}}{\tilde{\mu}_{1}+\tilde{\gamma}}}$.
Let
$\alpha_{1}:=e^{F(\tilde{\alpha}_{1})}=\frac{\tilde{\alpha}_{1}}{1+\kappa\tilde{\alpha}_{1}}$.
Let $t_{0}:=\frac{1}{\tilde{\mu}_{1}}(F(\tilde{\alpha}_{1})-\log\alpha)\leq
T_{\|{\bm{\delta}}_{\alpha}\|_{2}}(\tilde{\alpha}_{1})$. At time $t_{0}$, by
Lemma E.2 we have

$\left\|e^{t_{0}{\bm{J}}({\bm{0}})}-e^{\tilde{\mu}_{1}t_{0}}{\bm{e}}_{1}{\bm{e}}_{1}^{\top}\right\|_{2}=e^{\tilde{\mu}_{2}t_{0}}=e^{\frac{\tilde{\mu}_{2}}{\tilde{\mu}_{1}}(F(\tilde{\alpha}_{1})-\log\alpha)}=\left(\frac{\alpha_{1}}{\alpha}\right)^{\frac{\tilde{\mu}_{2}}{\tilde{\mu}_{1}}}.$
(26)

Let
$q_{\alpha}:=\left<\frac{{\bm{\delta}}_{\alpha}}{\alpha},{\bm{e}}_{1}\right>$.
By Definition 5.2, there exists $q>0$ such that $q_{\alpha}\geq q$ for all
sufficiently small $\alpha$. Then we have

$\displaystyle\left\|\phi\left({\bm{\delta}}_{\alpha},t_{0}\right)-\alpha_{1}q_{\alpha}{\bm{e}}_{1}\right\|_{2}$
$\displaystyle=\left\|\phi\left({\bm{\delta}}_{\alpha},t_{0}\right)-e^{t_{0}{\bm{J}}({\bm{0}})}{\bm{\delta}}_{\alpha}\right\|_{2}+\left\|\left(e^{t_{0}{\bm{J}}({\bm{0}})}-e^{\tilde{\mu}_{1}t_{0}}{\bm{e}}_{1}{\bm{e}}_{1}^{\top}\right){\bm{\delta}}_{\alpha}\right\|_{2}$
$\displaystyle=O(\tilde{\alpha}_{1}^{2})+\left(\frac{\alpha_{1}}{\alpha}\right)^{\frac{\tilde{\mu}_{2}}{\tilde{\mu}_{1}}}\|{\bm{\delta}}_{\alpha}\|_{2}$
$\displaystyle=O(\tilde{\alpha}_{1}^{2}+\alpha_{1}^{\tilde{\mu}_{2}/\tilde{\mu}_{1}}\alpha^{1-\tilde{\mu}_{2}/\tilde{\mu}_{1}})$
$\displaystyle=O(\alpha_{1}^{2}).$

Let $r:=F^{-1}(\tilde{\mu}_{1}t+\log\frac{1}{q_{\alpha}})\leq R$. Then
$t+\frac{1}{\tilde{\mu}_{1}}\log\frac{1}{\alpha_{1}q_{\alpha}}=T_{\tilde{\alpha}}(r)$
if $\tilde{\alpha}<r$. By Lemma E.3,
$\left\|\phi\left({\bm{\delta}}_{\alpha},t_{0}\right)\right\|_{2}\leq\tilde{\alpha}_{1}$.
Also,
$\|\alpha_{1}q_{\alpha}{\bm{e}}_{1}\|_{2}\leq\alpha_{1}=\frac{\tilde{\alpha}_{1}}{1+\kappa\tilde{\alpha}_{1}}\leq\tilde{\alpha}_{1}$.
By Lemma E.5,

$\displaystyle\left\|{\bm{\theta}}_{\alpha}(t)-{\bm{z}}_{\alpha_{1}}(t)\right\|_{2}$
$\displaystyle\leq\left\|\phi\left(\phi\left({\bm{\delta}}_{\alpha},t_{0}\right),t+\frac{1}{\tilde{\mu}_{1}}\log\frac{1}{\alpha_{1}q_{\alpha}}\right)-\phi\left(\alpha_{1}q_{\alpha}{\bm{e}}_{1},t+\frac{1}{\tilde{\mu}_{1}}\log\frac{1}{\alpha_{1}q_{\alpha}}\right)\right\|_{2}$
$\displaystyle=O\left(\alpha_{1}^{2}\cdot
e^{\tilde{\mu}_{1}\left(t+\frac{1}{\tilde{\mu}_{1}}\log\frac{1}{\alpha_{1}q_{\alpha}}\right)+\kappa
r}\right)$ $\displaystyle=O(\alpha_{1}).$

Combining this with the convergence rate for ${\bm{z}}_{\alpha_{1}}(t)$, we
have

$\left\|{\bm{\theta}}_{\alpha}(t)-{\bm{z}}(t)\right\|_{2}\leq\left\|{\bm{\theta}}_{\alpha}(t)-{\bm{z}}_{\alpha_{1}}(t)\right\|_{2}+\left\|{\bm{z}}_{\alpha_{1}}(t)-{\bm{z}}(t)\right\|_{2}=O(\alpha_{1}).$

##### Case 2.

For $t=(F(R)+\log q)/\tilde{\mu}_{1}+\tau$ with $\tau>0$,
$\phi({\bm{\theta}},\tau)$ is locally Lipschitz with respect to
${\bm{\theta}}$. So

$\displaystyle\left\|{\bm{\theta}}_{\alpha}(t)-{\bm{z}}(t)\right\|_{2}$
$\displaystyle=\left\|\phi({\bm{\theta}}_{\alpha}((F(R)+\log
q)/\tilde{\mu}_{1}),\tau)-\phi({\bm{z}}((F(R)+\log
q)/\tilde{\mu}_{1}),\tau)\right\|_{2}$
$\displaystyle=O(\left\|{\bm{\theta}}_{\alpha}((F(R)+\log
q)/\tilde{\mu}_{1})-{\bm{z}}((F(R)+\log q)/\tilde{\mu}_{1})\right\|_{2})$
$\displaystyle=O(\alpha_{1}),$

which proves (25) for $t>(F(R)+\log q)/\tilde{\mu}_{1}$. ∎

###  E.3 Extension to Non-Diagonalizable Case

The proof in Section E.2 can be generalized to the case where
${\bm{J}}({\bm{0}})$. Now we state the theorem formally and sketch the proof
idea. We use the notations
${\bm{g}}({\bm{\theta}}),\phi({\bm{\theta}}_{0},t),{\bm{J}}({\bm{\theta}})$ as
in Section 5.1, but we do not assume that ${\bm{J}}({\bm{0}})$ is
diagonalizable. Instead, we use
$\tilde{\mu}_{1},\tilde{\mu}_{2},\dots,\tilde{\mu}_{d}\in\mathbb{C}$ to denote
the eigenvalues of ${\bm{J}}({\bm{0}})$, repeated according to algebraic
multiplicity. We sort the eigenvalues in the descending order of the real part
of each eigenvalue, i.e.,
$\Re(\tilde{\mu}_{1})\geq\Re(\tilde{\mu}_{2})\geq\cdots\geq\Re(\tilde{\mu}_{d})$,
where $\Re(z)$ stands for the real part of a complex number $z\in\mathbb{C}$.
We call the eigenvalue with the largest real part the top eigenvalue.

######  Theorem E.7.

Assume that ${\bm{\theta}}={\bm{0}}$ is a critical point and the following
regularity conditions hold:

  1. 1.

${\bm{g}}({\bm{\theta}})$ is $\mathcal{C}^{2}$-smooth;

  2. 2.

$\phi({\bm{\theta}}_{0},t)$ exists for all ${\bm{\theta}}_{0}$ and $t$;

  3. 3.

The top eigenvalue of ${\bm{J}}({\bm{0}})$ is unique and is a positive real
number, i.e.,

$\tilde{\mu}_{1}>\max\\{\Re(\tilde{\mu}_{2}),0\\}.$

Let $\tilde{{\bm{v}}}_{1},\tilde{{\bm{u}}}_{1}$ be the left and right
eigenvectors associated with $\tilde{\mu}_{1}$, satisfying
$\tilde{{\bm{u}}}_{1}^{\top}\tilde{{\bm{v}}}_{1}=1$. Let
${\bm{z}}_{\alpha}(t):=\phi(\alpha\tilde{{\bm{v}}}_{1},t+\frac{1}{\tilde{\mu}_{1}}\log\frac{1}{\alpha})$
for every $\alpha>0$, then $\forall t\in\mathbb{R}$,
${\bm{z}}(t):=\lim\limits_{\alpha\to 0}{\bm{z}}_{\alpha}(t)$ exists and
${\bm{z}}(t)=\phi({\bm{z}}(0),t)$. If ${\bm{\delta}}_{\alpha}$ converges to
${\bm{0}}$ with positive alignment with $\tilde{{\bm{u}}}_{1}$ as $\alpha\to
0$, then for any $t\in\mathbb{R}$ and for any $\epsilon>0$, there is a
constant $C>0$ such that for every sufficiently small $\alpha$,

$\left\|\phi\left({\bm{\delta}}_{\alpha},t+\tfrac{1}{\tilde{\mu}_{1}}\log\tfrac{1}{\left<{\bm{\delta}}_{\alpha},\tilde{{\bm{u}}}_{1}\right>}\right)-{\bm{z}}(t)\right\|_{2}\leq
C\cdot\|{\bm{\delta}}_{\alpha}\|_{2}^{\frac{\tilde{\gamma}}{\tilde{\mu}_{1}+\tilde{\gamma}}-\epsilon},$
(27)

where $\tilde{\gamma}:=\tilde{\mu}_{1}-\Re(\tilde{\mu}_{2})$ is the eigenvalue
gap.

###### Proof Sketch.

Define the following two types of matrices. For $r\geq
1,a,\delta\in\mathbb{R}$, we define

${\bm{J}}^{(r)}_{a,\delta}:=\begin{bmatrix}a&\delta&&&&\\\ &a&\delta&&&\\\
&&a&\delta&&\\\ &&&\ddots&\ddots&\\\ &&&&a&\delta\\\
&&&&&a\end{bmatrix}\in\mathbb{R}^{r\times r}.$

For $r\geq 1,a,b,\delta\in\mathbb{R}$, we define

${\bm{J}}^{(r)}_{a,b,\delta}:=\begin{bmatrix}{\bm{C}}&\delta{\bm{I}}&&&&\\\
&{\bm{C}}&\delta{\bm{I}}&&&\\\ &&{\bm{C}}&\delta{\bm{I}}&&\\\
&&&\ddots&\ddots&\\\ &&&&{\bm{C}}&\delta{\bm{I}}\\\
&&&&&{\bm{C}}\end{bmatrix}\in\mathbb{R}^{2r\times 2r},$

where ${\bm{C}}=\left[\begin{smallmatrix}a&-b\\\
b&a\end{smallmatrix}\right]\in\mathbb{R}^{2\times 2}$.

By linear algebra, the real matrix ${\bm{J}}({\bm{0}})$ can be written in the
real Jordan normal form, i.e.,
${\bm{J}}({\bm{0}})=\tilde{\bm{V}}\mathrm{diag}({\bm{J}}_{[1]},\dots,{\bm{J}}_{[m]})\tilde{\bm{V}}^{-1}$,
where $\tilde{\bm{V}}\in\mathbb{R}^{d\times d}$ is an invertible matrix, and
each ${\bm{J}}_{[j]}$ is a real Jordan block. Recall that there are two types
of real Jordan blocks, ${\bm{J}}^{(r)}_{a,1}$ or ${\bm{J}}^{(r)}_{a,b,1}$. The
former one is associated with a real eigenvalue $a$, and the latter one is
associated with a pair of complex eigenvalues $a\pm bi$. The sum of sizes of
all Jordan blocks corresponding to a real eigenvalue $a$ is its algebraic
multiplicity. The sum of sizes of all Jordan blocks corresponding to a pair of
complex eigenvalues $a\pm bi$ is two times the algebraic multiplicity of
$a+bi$ or $a-bi$ (note that $a\pm bi$ have the same multiplicity).

It is easy to see that
${\bm{J}}^{(r)}_{a,\delta}={\bm{D}}{\bm{J}}^{(r)}_{a,1}{\bm{D}}^{-1}$ for
${\bm{D}}=\mathrm{diag}(\delta^{r},\delta^{r-1},\dots,\delta)\in\mathbb{R}^{r\times
r}$ and
${\bm{J}}^{(r)}_{a,b,\delta}={\bm{D}}{\bm{J}}^{(r)}_{a,b,1}{\bm{D}}^{-1}$ for
${\bm{D}}=\mathrm{diag}(\delta^{r},\delta^{r},\delta^{r-1},\delta^{r-1},\dots,\delta,\delta)\in\mathbb{R}^{2r\times
2r}$. This means for every $\delta>0$ there exists $\tilde{\bm{V}}_{\delta}$
such that
${\bm{J}}({\bm{0}})=\tilde{\bm{V}}_{\delta}{\bm{J}}_{\delta}\tilde{\bm{V}}_{\delta}^{-1}$,
where
${\bm{J}}_{\delta}:=\mathrm{diag}({\bm{J}}_{\delta[1]},\dots,{\bm{J}}_{\delta[m]})$,
${\bm{J}}_{\delta[j]}:={\bm{J}}^{(r)}_{a,\delta}$ if
${\bm{J}}_{[j]}:={\bm{J}}^{(r)}_{a,1}$, or
${\bm{J}}_{\delta[j]}:={\bm{J}}^{(r)}_{a,b,\delta}$ if
${\bm{J}}_{[j]}:={\bm{J}}^{(r)}_{a,b,1}$. Since the top eigenvalue of
${\bm{J}}({\bm{0}})$ is positive and unique, $\tilde{\mu}_{1}$ corresponds to
only one block $[\tilde{\mu}_{1}]\in\mathbb{R}^{1\times 1}$. WLOG we let
${\bm{J}}_{1}=[\tilde{\mu}_{1}]$, and thus
${\bm{J}}_{\delta[1]}=[\tilde{\mu}_{1}]$.

We only need to select a parameter $\delta>0$ and prove the theorem in the
case of ${\bm{J}}({\bm{0}})={\bm{J}}_{\delta}$ since we can change the basis
in a similar way as we have done in Section E.1. By scrutinizing the proof for
Theorem E.1, we can find that we only need to reprove Lemma E.2. However,
Lemma E.2 may not be correct since ${\bm{J}}({\bm{0}})$ is not diagonal
anymore. Instead, we prove the following:

  1. 1.

If $\delta\in(0,\tilde{\gamma})$, then
${\bm{h}}^{\top}{\bm{J}}_{\delta}{\bm{h}}\leq\tilde{\mu}_{1}\|{\bm{h}}\|_{2}^{2}$
for all ${\bm{h}}\in\mathbb{R}^{d}$;

  2. 2.

For any $\tilde{\mu}^{\prime}_{2}\in(\Re(\tilde{\mu}_{2}),\tilde{\mu}_{1})$,
if $\delta\in(0,\tilde{\mu}^{\prime}_{2}-\Re(\tilde{\mu}_{2}))$, then
$\left\|e^{t{\bm{J}}_{\delta}}-e^{\tilde{\mu}_{1}t}{\bm{e}}_{1}{\bm{e}}_{1}^{\top}\right\|_{2}\leq
e^{\tilde{\mu}^{\prime}_{2}t}$ for all $t\geq 0$.

##### Proof for Item 1.

Let ${\mathcal{K}}$ be the set of pairs $(k_{1},k_{2})$ such that $k_{1}\neq
k_{2}$ and the entry of ${\bm{J}}_{\delta}$ at the $k_{1}$-th row and the
$k_{2}$-th column is non-zero. Then we have

$\displaystyle{\bm{h}}^{\top}{\bm{J}}_{\delta}{\bm{h}}={\bm{h}}^{\top}\frac{{\bm{J}}_{\delta}+{\bm{J}}_{\delta}^{\top}}{2}{\bm{h}}$
$\displaystyle=\sum_{k=1}^{d}\Re(\tilde{\mu}_{k})h_{k}^{2}+\sum_{(k_{1},k_{2})\in{\mathcal{K}}}h_{k_{1}}h_{k_{2}}\delta$
$\displaystyle\leq\sum_{k=1}^{d}\Re(\tilde{\mu}_{k})h_{k}^{2}+\sum_{(k_{1},k_{2})\in{\mathcal{K}}}\frac{h^{2}_{k_{1}}+h^{2}_{k_{2}}}{2}\delta.$

Note that $\Re(\tilde{\mu}_{k})\leq\Re(\tilde{\mu}_{2})$ for $k\geq 2$. Also
note that there is no pair in ${\mathcal{K}}$ has $k_{1}=1$ or $k_{2}=1$, and
for every $k\geq 2$ there are at most two pairs in ${\mathcal{K}}$ has
$k_{1}=k$ or $k_{2}=k$. Combining all these together gives

${\bm{h}}^{\top}{\bm{J}}_{\delta}{\bm{h}}\leq\tilde{\mu}_{1}h_{1}^{2}+(\Re(\tilde{\mu}_{2})+\delta)\sum_{k=2}^{d}h_{k}^{2}\leq\tilde{\mu}_{1}\|{\bm{h}}\|_{2}^{2},$

which proves Item 1.

##### Proof for Item 2.

Since ${\bm{J}}_{\delta}$ is block diagonal, we only need to prove that
$\|e^{t{\bm{J}}_{\delta[j]}}\|_{2}\leq e^{\tilde{\mu}^{\prime}_{2}t}$ for
every $j\geq 2$. If
${\bm{J}}_{\delta[j]}={\bm{J}}^{(r)}_{a,\delta}=a{\bm{I}}+\delta{\bm{N}}$,
where ${\bm{N}}$ is the nilpotent matrix, then

$e^{t{\bm{J}}_{\delta[j]}}=e^{at{\bm{I}}+\delta
t{\bm{N}}}=e^{at{\bm{I}}}e^{\delta t{\bm{N}}}=e^{at}e^{\delta t{\bm{N}}},$

where the second equality uses the fact that ${\bm{I}}$ and ${\bm{N}}$ are
commutable. So we have

$\|e^{t{\bm{J}}_{\delta[j]}}\|_{2}\leq e^{at}\|e^{\delta
t{\bm{N}}}\|_{2}=e^{at}e^{\delta t\|{\bm{N}}\|_{2}}\leq e^{(a+\delta)t}.$

If
${\bm{J}}_{\delta[j]}={\bm{J}}^{(r)}_{a,\delta}={\bm{D}}+\delta{\bm{N}}^{2}$,
where ${\bm{D}}=\mathrm{diag}({\bm{C}},{\bm{C}},\dots,{\bm{C}})$ and
${\bm{N}}$ is the nilpotent matrix, then

$e^{t{\bm{J}}_{\delta[j]}}=e^{t{\bm{D}}+\delta
t{\bm{N}}^{2}}=e^{t{\bm{D}}}e^{\delta t{\bm{N}}^{2}},$

where the second equality uses the fact that ${\bm{D}}$ and ${\bm{N}}^{2}$ are
commutable. Note that
$e^{t{\bm{C}}}=e^{at}\left[\begin{smallmatrix}\cos(bt)&-\sin(bt)\\\
\sin(bt)&\cos(bt)\end{smallmatrix}\right]$, which implies
$\|e^{t{\bm{D}}}\|_{2}=\|e^{t{\bm{C}}}\|_{2}=e^{at}$. So we have

$\|e^{t{\bm{J}}_{\delta[j]}}\|_{2}\leq\|e^{t{\bm{D}}}\|_{2}\cdot\|e^{\delta
t{\bm{N}}^{2}}\|_{2}=e^{at}e^{\delta t\|{\bm{N}}^{2}\|_{2}}\leq
e^{(a+\delta)t}.$

Since $\delta\in(0,\tilde{\mu}^{\prime}_{2}-\Re(\tilde{\mu}_{2}))$, we know
that $a+\delta<\tilde{\mu}^{\prime}_{2}$, which completes the proof.

##### Proof for a fixed $\delta$.

Since Item 1 continues to hold for $\delta\in(0,\tilde{\gamma})$, Lemmas E.3,
E.4, E.5 and E.6 also hold. This proves that ${\bm{z}}(t)$ exists and
satisfies (6).

It remains to prove (27) for any $\epsilon>0$. Let
$\tilde{\gamma}^{\prime}\in(0,\tilde{\gamma})$ be a number such that
$\frac{\tilde{\gamma}^{\prime}}{\tilde{\mu}_{1}+\tilde{\gamma}^{\prime}}\geq\frac{\tilde{\gamma}}{\tilde{\mu}_{1}+\tilde{\gamma}}-\epsilon$.
Fix
$\tilde{\mu}^{\prime}_{2}=\tilde{\mu}_{1}-\tilde{\gamma}^{\prime},\delta=\tilde{\mu}^{\prime}_{2}-\Re(\tilde{\mu}_{2})$.
By Item 2, we have
$\left\|e^{t{\bm{J}}_{\delta}}-e^{\tilde{\mu}_{1}t}{\bm{e}}_{1}{\bm{e}}_{1}^{\top}\right\|_{2}\leq
e^{\tilde{\mu}^{\prime}_{2}t}$ for all $t\geq 0$. By scrutinizing the proof
for Theorem E.1, we can find that the only place we use Item 2 in Lemma E.2 is
in (26). For proving (27), we can repeat the proof while replacing all the
occurrences of $\tilde{\mu}_{2}$ by $\tilde{\mu}^{\prime}_{2}$. Then we know
that for every $t\in\mathbb{R}$, there is a constant $C>0$ such that

$\left\|\tilde{\bm{V}}_{\delta}^{-1}\phi\left({\bm{\delta}}_{\alpha},t+\tfrac{1}{\tilde{\mu}_{1}}\log\tfrac{1}{\left<{\bm{\delta}}_{\alpha},\tilde{{\bm{u}}}_{1}\right>}\right)-\tilde{\bm{V}}_{\delta}^{-1}{\bm{z}}(t)\right\|_{2}\leq
C\cdot\|\tilde{\bm{V}}_{\delta}^{-1}{\bm{\delta}}_{\alpha}\|_{2}^{\frac{\tilde{\gamma}^{\prime}}{\tilde{\mu}_{1}+\tilde{\gamma}^{\prime}}},$

for every sufficiently small $\alpha$. By definition of
$\tilde{\gamma}^{\prime}$,
$\frac{\tilde{\gamma}^{\prime}}{\tilde{\mu}_{1}+\tilde{\gamma}^{\prime}}\geq\frac{\tilde{\gamma}}{\tilde{\mu}_{1}+\tilde{\gamma}}-\epsilon$.
Since ${\bm{\delta}}_{\alpha}\to{\bm{0}}$ as $\alpha\to 0$, we have
$\|\tilde{\bm{V}}_{\delta}^{-1}{\bm{\delta}}_{\alpha}\|_{2}<1$ for
sufficiently small $\alpha$. Then we have

$\displaystyle\left\|\phi\left({\bm{\delta}}_{\alpha},t+\tfrac{1}{\tilde{\mu}_{1}}\log\tfrac{1}{\left<{\bm{\delta}}_{\alpha},\tilde{{\bm{u}}}_{1}\right>}\right)-{\bm{z}}(t)\right\|_{2}$
$\displaystyle\leq\|\tilde{\bm{V}}_{\delta}\|_{2}\cdot\left\|\tilde{\bm{V}}_{\delta}^{-1}\phi\left({\bm{\delta}}_{\alpha},t+\tfrac{1}{\tilde{\mu}_{1}}\log\tfrac{1}{\left<{\bm{\delta}}_{\alpha},\tilde{{\bm{u}}}_{1}\right>}\right)-\tilde{\bm{V}}_{\delta}^{-1}{\bm{z}}(t)\right\|_{2}$
$\displaystyle\leq\|\tilde{\bm{V}}_{\delta}\|_{2}\cdot
C\cdot\|\tilde{\bm{V}}_{\delta}^{-1}{\bm{\delta}}_{\alpha}\|_{2}^{\frac{\tilde{\gamma}^{\prime}}{\tilde{\mu}_{1}+\tilde{\gamma}^{\prime}}}$
$\displaystyle\leq
C\cdot\|\tilde{\bm{V}}_{\delta}\|_{2}\cdot\|\tilde{\bm{V}}_{\delta}^{-1}\|_{2}^{\frac{\tilde{\gamma}^{\prime}}{\tilde{\mu}_{1}+\tilde{\gamma}^{\prime}}}\cdot\|{\bm{\delta}}_{\alpha}\|_{2}^{\frac{\tilde{\gamma}}{\tilde{\mu}_{1}+\tilde{\gamma}}-\epsilon}.$

Absorbing
$\|\tilde{\bm{V}}_{\delta}\|_{2}\cdot\|\tilde{\bm{V}}_{\delta}^{-1}\|_{2}^{\frac{\tilde{\gamma}^{\prime}}{\tilde{\mu}_{1}+\tilde{\gamma}^{\prime}}}$
into $C$ proves (27). ∎

##  Appendix F Eigenvalues of Jacobians and Hessians

In this section we analyze the eigenvalues of the Jacobian
${\bm{J}}({\bm{W}})$ at critical points of (2).

For notation simplicity, we write
$\mathrm{sz}({\bm{A}}):={\bm{A}}+{\bm{A}}^{\top}$ to denote the symmetric
matrix produced by adding up ${\bm{A}}$ and its transpose, and write
$\mathrm{ac}\\{{\bm{A}},{\bm{B}}\\}={\bm{A}}{\bm{B}}+{\bm{B}}{\bm{A}}$ to
denote the anticommutator of two matrices ${\bm{A}},{\bm{B}}$. Then
${\bm{g}}({\bm{W}})$ can be written as
${\bm{g}}({\bm{W}}):=-\mathrm{ac}\\{\nabla f({\bm{W}}),{\bm{W}}\\}$.

Let ${\bm{U}}_{0}\in\mathbb{R}^{d\times r}$ be a stationary point of the
function $\mathcal{L}:\mathbb{R}^{d\times
r}\to\mathbb{R},\mathcal{L}({\bm{U}})=\frac{1}{2}f({\bm{U}}{\bm{U}}^{\top})$,
i.e., $\nabla\mathcal{L}({\bm{U}}_{0})=\nabla
f({\bm{U}}_{0}{\bm{U}}_{0}^{\top}){\bm{U}}_{0}={\bm{0}}$. By Lemma A.1, this
implies

$\nabla f({\bm{W}}_{0}){\bm{W}}_{0}={\bm{0}}$ (28)

for ${\bm{W}}_{0}:={\bm{U}}_{0}{\bm{U}}_{0}^{\top}$, and thus ${\bm{W}}_{0}$
is a critical point of (2).

For a real-valued or vector-valued function $F({\bm{\theta}})$, we use
$DF({\bm{\theta}})[{\bm{\delta}}],D^{2}F({\bm{\theta}})[{\bm{\delta}}_{1},{\bm{\delta}}_{2}]$
to denote the first- and second-order directional derivatives of
$F(\,\cdot\,)$ at ${\bm{\theta}}$.

Let $\mathcal{X}$ be a linear space, which can be $\mathbb{R}^{d\times d}$ or
$\mathbb{R}^{d\times r}$. For a function $F:\mathcal{X}\to\mathcal{X}$, we use
$DF({\bm{\theta}})$ to denote the directional derivative of $F$ at
${\bm{\theta}}$, represented by the linear operator

$DF({\bm{\theta}})[{\bm{\Delta}}]:\mathcal{X}\to\mathcal{X},{\bm{\Delta}}\mapsto
DF({\bm{\theta}})[{\bm{\Delta}}]=\lim_{t\to
0}\frac{F({\bm{\theta}}+t{\bm{\Delta}})-F({\bm{\theta}})}{t}.$

We also write
$DF({\bm{\theta}})[{\bm{\Delta}}_{1},{\bm{\Delta}}_{2}]:=\left<DF({\bm{\theta}})[{\bm{\Delta}}_{1}],{\bm{\Delta}}_{2}\right>$.

For a function $F:\mathcal{X}\to\mathbb{R}$, we use
$D^{2}F({\bm{\theta}})=D(\nabla F({\bm{\theta}}))$ to denote the second
directional derivative of $F$ at ${\bm{\theta}}$, i.e.,
$D^{2}F({\bm{\theta}})[{\bm{\Delta}}]=D(\nabla
F({\bm{\theta}}))[{\bm{\Delta}}],D^{2}F({\bm{\theta}})[{\bm{\Delta}}_{1},{\bm{\Delta}}_{2}]=D(\nabla
F({\bm{\theta}}))[{\bm{\Delta}}_{1},{\bm{\Delta}}_{2}]$.

Define ${\bm{J}}({\bm{W}}):=D{\bm{g}}({\bm{W}})$. By simple calculus, we can
compute the formula for ${\bm{J}}({\bm{W}}_{0})$:

$\displaystyle{\bm{J}}({\bm{W}}_{0})[{\bm{\Delta}}]$
$\displaystyle=-\mathrm{ac}\\{\nabla
f({\bm{W}}_{0}),{\bm{\Delta}}\\}-\mathrm{ac}\\{D^{2}f({\bm{W}}_{0})[{\bm{\Delta}}],{\bm{W}}_{0}\\},$
$\displaystyle{\bm{J}}({\bm{W}}_{0})[{\bm{\Delta}}_{1},{\bm{\Delta}}_{2}]$
$\displaystyle=-\left<\nabla
f({\bm{W}}_{0}),\mathrm{sz}({\bm{\Delta}}_{1}{\bm{\Delta}}_{2}^{\top})\right>-D^{2}f({\bm{W}}_{0})[{\bm{\Delta}}_{1},\mathrm{sz}({\bm{W}}_{0}{\bm{\Delta}}_{2}^{\top})],$

where ${\bm{\Delta}},{\bm{\Delta}}_{1},{\bm{\Delta}}_{2}\in\mathbb{R}^{d\times
d}$.

We can also compute the formula for $D^{2}\mathcal{L}({\bm{U}}_{0})$:

$\displaystyle D^{2}\mathcal{L}({\bm{U}}_{0})[{\bm{\Delta}}]$
$\displaystyle=\nabla
f({\bm{W}}_{0}){\bm{\Delta}}+D^{2}f({\bm{W}}_{0})[\mathrm{sz}({\bm{\Delta}}{\bm{U}}_{0}^{\top})]{\bm{U}}_{0},$
$\displaystyle
D^{2}\mathcal{L}({\bm{U}}_{0})[{\bm{\Delta}}_{1},{\bm{\Delta}}_{2}]$
$\displaystyle=\frac{1}{2}\left(\left<\nabla
f({\bm{W}}_{0}),\mathrm{sz}({\bm{\Delta}}_{1}{\bm{\Delta}}_{2}^{\top})\right>+D^{2}f({\bm{W}}_{0})[\mathrm{sz}({\bm{\Delta}}_{1}{\bm{U}}_{0}^{\top}),\mathrm{sz}({\bm{\Delta}}_{2}{\bm{U}}_{0}^{\top})]\right),$

where ${\bm{\Delta}},{\bm{\Delta}}_{1},{\bm{\Delta}}_{2}\in\mathbb{R}^{d\times
r}$.

###  F.1 Eigenvalues at the Origin

The eigenvalues of ${\bm{J}}({\bm{0}})$ is given in Lemma 5.4. Now we provide
the proof.

###### Proof for Lemma 5.4.

For ${\bm{W}}_{0}={\bm{0}}$, we have

$\displaystyle{\bm{J}}({\bm{0}})[{\bm{\Delta}}]$ $\displaystyle=-\nabla
f({\bm{0}}){\bm{\Delta}}-{\bm{\Delta}}\nabla f({\bm{0}})$
$\displaystyle{\bm{J}}({\bm{0}})[{\bm{\Delta}}_{1},{\bm{\Delta}}_{2}]$
$\displaystyle=-\left<\nabla
f({\bm{0}}),\mathrm{sz}({\bm{\Delta}}_{1}{\bm{\Delta}}_{2}^{\top})\right>$

It is easy to see from the second equation that ${\bm{J}}({\bm{0}})$ is
symmetric.

Let $-\nabla
f({\bm{0}})=\sum_{i=1}^{d}\mu_{i}{\bm{u}}_{1[i]}{\bm{u}}_{1[i]}^{\top}$ be the
eigendecomposition of the symmetric matrix $-\nabla f({\bm{0}})$. Then we have

$\displaystyle{\bm{J}}({\bm{0}})[{\bm{\Delta}}]$
$\displaystyle=\sum_{i=1}^{d}\mu_{i}\left({\bm{u}}_{1[i]}{\bm{u}}_{1[i]}^{\top}{\bm{\Delta}}+{\bm{\Delta}}{\bm{u}}_{1[i]}{\bm{u}}_{1[i]}^{\top}\right)$
$\displaystyle=\sum_{i=1}^{d}\sum_{j=1}^{d}\mu_{i}\left({\bm{u}}_{1[i]}{\bm{u}}_{1[i]}^{\top}{\bm{\Delta}}{\bm{u}}_{1[j]}{\bm{u}}_{1[j]}^{\top}+{\bm{u}}_{1[j]}{\bm{u}}_{1[j]}^{\top}{\bm{\Delta}}{\bm{u}}_{1[i]}{\bm{u}}_{1[i]}^{\top}\right)$
$\displaystyle=\sum_{i=1}^{d}\sum_{j=1}^{d}(\mu_{i}+\mu_{j}){\bm{u}}_{1[i]}{\bm{u}}_{1[i]}^{\top}{\bm{\Delta}}{\bm{u}}_{1[j]}{\bm{u}}_{1[j]}^{\top}$
$\displaystyle=\sum_{i=1}^{d}\sum_{j=1}^{d}(\mu_{i}+\mu_{j})\left<{\bm{\Delta}},{\bm{u}}_{1[i]}{\bm{u}}_{1[j]}^{\top}\right>{\bm{u}}_{1[i]}{\bm{u}}_{1[j]}^{\top},$

which proves (8).

For
${\bm{\Delta}}={\bm{u}}_{1[i]}{\bm{u}}_{1[j]}^{\top}+{\bm{u}}_{1[j]}{\bm{u}}_{1[i]}^{\top}$,
we have