# Connectivity Shapes Implicit Regularization in Matrix Factorization Models
for Matrix Completion

Zhiwei Bai1, Jiajie Zhao1, Yaoyu Zhang1,2
1 School of Mathematical Sciences, Institute of Natural Sciences, MOE-LSC,
Shanghai Jiao Tong University, Shanghai 200240, P.R. China.
2 Shanghai Center for Brain Science and Brain-Inspired Technology, Shanghai
200240, P.R. China
{bai299, zjj0216<EMAIL_ADDRESS>Corresponding author:
<EMAIL_ADDRESS>
###### Abstract

Matrix factorization models have been extensively studied as a valuable test-
bed for understanding the implicit biases of overparameterized models.
Although both low nuclear norm and low rank regularization have been studied
for these models, a unified understanding of when, how, and why they achieve
different implicit regularization effects remains elusive. In this work, we
systematically investigate the implicit regularization of matrix factorization
for solving matrix completion problems. We empirically discover that the
connectivity of observed data plays a crucial role in the implicit bias, with
a transition from low nuclear norm to low rank as data shifts from
disconnected to connected with increased observations. We identify a hierarchy
of intrinsic invariant manifolds in the loss landscape that guide the training
trajectory to evolve from low-rank to higher-rank solutions. Based on this
finding, we theoretically characterize the training trajectory as following
the hierarchical invariant manifold traversal process, generalizing the
characterization of Li et al. (2020) to include the disconnected case.
Furthermore, we establish conditions that guarantee minimum nuclear norm,
closely aligning with our experimental findings, and we provide a dynamics
characterization condition for ensuring minimum rank. Our work reveals the
intricate interplay between data connectivity, training dynamics, and implicit
regularization in matrix factorization models.

##  1 Introduction

Overparameterized models have the capacity to easily fit data with random
labels (Zhang et al., 2017, 2021). However, in real-world applications, models
with more parameters than training samples still generalize well. This has led
researchers to hypothesize that overparameterized models undergo implicit
regularization, favoring certain functions as outputs. Overparameterized
matrix factorization models, $\bm{f}_{\bm{\theta}}=\bm{A}\bm{B}$ with
$\bm{\theta}=(\bm{A},\bm{B}),\bm{A},\bm{B}\in\mathbb{R}^{d\times d}$, have
served as a simplified test-bed for studying this implicit regularization. In
the context of matrix completion problems like the Netflix challenge, these
models aim to find a low-rank completion of a partially observed matrix
$\bm{M}\in\mathbb{R}^{d\times d}$. Prior works have offered seemingly
conflicting perspectives on the implicit regularization at play, with some
claiming it promotes low nuclear norm (Gunasekar et al., 2017) and others
arguing for low rank (Arora et al., 2019; Li et al., 2020; Razin and Cohen,
2020). However, a unified understanding of when, how, and why they achieve
different implicit regularization effects remains elusive.

Unlike previous works that focus on either low rank or low nuclear norm
regularization,, we systematically investigate the training dynamics and
implicit regularization of matrix factorization for matrix completion. Through
extensive experiments, we found that a certain connectivity property of
observed data (see Def. 2 for the definition of connectivity) plays a key role
in the implicit regularization effects. As shown in Fig. 1, we sample
observations randomly from a ground truth matrix
$\bm{M}^{*}\in\mathbb{R}^{d\times d}$ with $\textnormal{rank}(\bm{M}^{*})<d$
and train models
$\bm{f}_{\bm{\theta}}=\bm{A}\bm{B},\bm{A},\bm{B}\in\mathbb{R}^{d\times d}$
from small random initialization without any rank constraints. For each
observation set, we calculate the solutions with the minimum nuclear norm and
minimum rank, which serve as the ground truth benchmarks. These are then
compared with the completion matrix obtained by the model. From Fig. 1, we
observe that:

(i) Low rank bias in connected case: When the observed entries are connected,
the model consistently learns the lowest-rank solution.

(ii) Low nuclear norm bias in certain disconnected case: When the observed
entries are disconnected, the model generally does not find the minimum
nuclear norm or lowest-rank solution. However, in the special case where each
connected component is a complete bipartite subgraph, the model consistently
finds the minimum nuclear norm solution.

To understand how data connectivity modulates the implicit bias, we analyze
the loss landscape and optimization dynamics. We find a hierarchy of intrinsic
invariant manifolds $\bm{\Omega}_{k}$ of different ranks in the loss
landscape. These manifolds constrain the optimization trajectory, causing the
model to learn by incrementally ascending through higher ranks. In the
disconnected case, additional sub-$\bm{\Omega}_{k}$ invariant manifolds emerge
within the $\bm{\Omega}_{k}$ invariant manifold, preventing the model from
reaching the global lowest-rank solution. However, we prove that the minimum
nuclear norm solution is guaranteed in the disconnected with complete
bipartite subgraph case.

The contributions of our work are summarized as follows:

(i) We systematically investigate the influence of data connectivity on the
implicit regularization. Our empirical findings indicate that the connectivity
of observed data plays a key role in the implicit bias, leading to a
transition from favoring solutions with a low nuclear norm to those with a low
rank as the data becomes more connected with an increase in observations
(refer to Sec. 4).

(ii) We characterize the training dynamics of matrix factorization
theoretically, showing that the optimization trajectory follows a Hierarchical
Invariant Manifold Traversal (HIMT) process. This generalizes the
characterization of Li et al. (2020), whose proposed Greedy Low-Rank
Learning(GLRL) algorithm equivalence only corresponding to the connected case
(refer to Sec. 5 and Sec. 6.1).

(iii) Regarding the minimum nuclear norm regularization, we establish
conditions that provide guarantees closely aligned with our empirical
findings, which complement the results of Gunasekar et al. (2017). For the
minimum rank regularization, we present a dynamic characterization condition
that assures the attainment of the minimum rank solution (refer to Sec. 6.2).

Figure 1: The connectivity of observed data affects the implicit
regularization. The ground truth matrix $\bm{M}^{*}\in\mathbb{R}^{4\times 4}$
has rank ranging from 1 to 3. The sample size $n$ covers settings where $n$ is
equal to, smaller than, and larger than the $2rd-r^{2}$ threshold required for
exact reconstruction. Darker scatter points indicate a greater number of
samples, while lighter points indicate fewer samples. The positions of
observed entries are randomly chosen, and the experiment is repeated 10 times
for each sample size. (Please refer to Appendix B for additional experiments
and detailed methodology.)

##  2 Related works

##### Norm minimization and rank minimization.

Extensive research has been conducted on the implicit regularization of matrix
factorization models, focusing on norm minimization and rank minimization. For
norm minimization, Gunasekar et al. (2017) proved that gradient flow with
infinitesimal initialization converges to the minimum nuclear norm solution in
the special case of commutative observations. Ji and Telgarsky (2019);
Gunasekar et al. (2018) studied norm minimization regularization in deep
linear networks. For rank minimization, numerous works have shown that matrix
factorization models favor low-rank solutions. Arora et al. (2019); Gidel et
al. (2019); Gissin et al. (2019); Razin and Cohen (2020); Jiang et al. (2023);
Belabbas (2020) investigated how infinitesimal initialization of gradient flow
encourages low rank in specific settings. Li et al. (2020) showed that under
certain assumptions, matrix factorization dynamics are equivalent to a greedy
low-rank learning heuristic. Li et al. (2018); Stöger and Soltanolkotabi
(2021); Jin et al. (2023) established low-rank recovery guarantees for matrix
sensing problems under the Restricted Isometry Property (RIP) condition. Zhang
et al. (2022, 2023) studied a broader class of model rank minimization for
nonlinear models, of which the matrix factorization model is a special case.

##### Nonlinear dynamics.

The initialization scale can significantly influence the implicit
regularization of neural networks. Large initialization typically leads to
linear dynamics (Jacot et al., 2018) and poor generalization (Chizat et al.,
2019), while small initialization induces nonlinear dynamics (Luo et al.,
2021). In this work, we focus on the case of infinitesimal initialization,
which corresponds to highly nonlinear dynamics. An important characteristic of
nonlinear neural network dynamics is the phenomenon of condensation (Luo et
al., 2021; Zhou et al., 2022), where the network’s effective complexity is
small. The low-rank $\bm{\Omega}_{k}$ invariant manifolds we propose are
essentially a manifestation of condensation. Zhang et al. (2021, 2022); Bai et
al. (2022); Fukumizu et al. (2019); Simsek et al. (2021) established the
embedding principle of the loss landscape of neural networks and empirically
demonstrated that the training process traverses critical points embedded from
smaller subnetworks. Jacot et al. (2021) conjectured a saddle to saddle
dynamics for deep linear networks, which is conceptually analogous to the
dynamics characterization in this work.

##  3 Preliminaries

##### Matrix completion problem.

This study focuses on the matrix completion problem, which involves estimating
missing entries within a partially observed matrix. Given an incomplete matrix
$\bm{M}\in\mathbb{R}^{d\times d}$, the goal is to predict the entirety of
$\bm{M}$ based on its observed elements. The set of observed entries is
represented as $S=\\{(i_{k},j_{k}),\bm{M}_{i_{k},j_{k}}\\}_{k=1}^{n}$, where
$(i_{k},j_{k})$ indicates the row and column indices, and
$\bm{M}_{i_{k},j_{k}}$ is the corresponding value assumed non-zero in the
matrix. The set of observed indices is defined as
$S_{\bm{x}}=\\{(i_{k},j_{k})\\}_{k=1}^{n}$. Entries that are not observed,
denoted by $\star$, are considered missing or unknown. The positions of
observed elements in the matrix $\bm{M}$ are defined by a binary observation
matrix $\bm{P}$, where $\bm{P}_{ij}=1$ indicates that $\bm{M}_{ij}$ is
observed, and $\bm{P}_{ij}=0$ indicates that $\bm{M}_{ij}$ is unobserved.

##### Matrix factorization model.

Matrix factorization is a prevalent approach for addressing the matrix
completion problem. It reconstructs the matrix $\bm{W}\in\mathbb{R}^{d\times
d}$ through the product $\bm{W}=\bm{A}\bm{B}$, where
$\bm{A}\in\mathbb{R}^{d\times r}$ and $\bm{B}\in\mathbb{R}^{r\times d}$. This
work studies the overparameterized scenario with $r=d$, aiming to understand
the implicit regularization effect in the absence of explicit rank
restrictions, paralleling prior research (Gunasekar et al., 2017; Arora et
al., 2019; Li et al., 2020; Jin et al., 2023). In this work, we focus on the
asymmetric factorization, which can be represented as a parametric model:

$\bm{f}_{\bm{\theta}}=\bm{A}\bm{B},\quad\bm{A},\bm{B}\in\mathbb{R}^{d\times
d}.$ (1)

The matrix factorization model parameters are denoted by
$\bm{\theta}=(\bm{A},\bm{B})$, identified with its vectorized form
$\mathrm{vec}(\bm{\theta})\in\mathbb{R}^{2d^{2}}$. The augmented matrix is
$\bm{W}_{\mathrm{aug}}^{\top}=\left[\begin{matrix}\bm{A}^{\top}&\bm{B}\end{matrix}\right]^{\top}\in\mathbb{R}^{d\times
2d}$, and $\textnormal{row}(\bm{A})$ and $\textnormal{col}(\bm{B})$ denote the
row and column spaces of $\bm{A}$ and $\bm{B}$, respectively.

##### Loss function.

The learning process for the parameters $\bm{\theta}=(\bm{A},\bm{B})$ involves
minimizing a loss function that measures the difference between observed and
estimated entries. In this work, we focus on the mean squared error, and the
empirical risk is thus formulated as

$R_{S}(\bm{\theta})=\frac{1}{n}\|(\bm{A}\bm{B}-\bm{M})_{S_{\bm{x}}}\|_{F}^{2}:=\dfrac{1}{n}\sum\limits_{k=1}^{n}(\bm{a}_{i_{k}}\cdot\bm{b}_{\cdot,j_{k}}-\bm{M}_{i_{k},j_{k}})^{2},$
(2)

where $\bm{a}_{i}$ and $\bm{b}_{\cdot,j}$ represent the $i$-th row and $j$-th
column of matrix $\bm{A}$ and $\bm{B}$, respectively. The residual matrix
$\delta\bm{M}=(\bm{A}\bm{B}-\bm{M})_{S_{\bm{x}}}$ has elements
$\delta\bm{M}_{ij}=(\bm{A}\bm{B})_{ij}-\bm{M}_{ij}$ for $(i,j)\in S_{\bm{x}}$
and $\delta\bm{M}_{ij}=0$ for $(i,j)\notin S_{\bm{x}}$. The training dynamics
follow the gradient flow of $R_{S}(\bm{\theta})$:

$\dfrac{\mathrm{d}\bm{\theta}}{\mathrm{d}t}=-\nabla_{\bm{\theta}}R_{S}(\bm{\theta}),\quad\bm{\theta}(0)=\bm{\theta}_{0}.$
(3)

##  4 Connectivity affects implicit regularization

In this section, we define connectivity and present experimental results on
implicit regularization for connected and disconnected observational data.

######  Definition 1 (Associated Observation Graph).

Given a incomplete matrix $\bm{M}$ to be completed and its observation matrix
$\bm{P}$, the associated observation graph $G_{\bm{M}}$ is the bipartite graph
with adjacency matrix $\begin{bmatrix}\bm{0}&\bm{P}^{\top}\\\
\bm{P}&\bm{0}\end{bmatrix}$, with isolated vertices removed.

######  Definition 2 (Connectivity).

Given a incomplete matrix $\bm{M}$ to be completed, it is considered connected
if its associated observation graph $G_{\bm{M}}$ is connected; otherwise, it
is disconnected. The connected components of $\bm{M}$ are defined as the
connected components of $G_{\bm{M}}$.

The connectivity of the graph, as defined above, reflects the connectivity of
the observed data. Appendix A Sec. A.2 provides a detailed discussion on the
equivalent definition of connectivity.

In the case of disconnectivity, there is a special case where each connected
component has full observations, characterized by disconnectivity with
complete bipartite components.

######  Definition 3 (Disconnectivity with Complete Bipartite Components).

A incomplete matrix $\bm{M}$ is considered disconnected with complete
bipartite components if its associated observation graph $G_{\bm{M}}$ is
disconnected and each connected component forms a complete bipartite subgraph.

We present examples to demonstrate how connectivity influences the
characteristics of the learned solutions. Consider three matrices to be
completed, each obtained by adding one more observation to the previous
matrix: $\bm{M}_{1}$ (disconnected), $\bm{M}_{2}$ (disconnected with complete
bipartite components), and $\bm{M}_{3}$ (connected). Fig. A1 of Appendix B
illustrates the associated graphs $G_{\bm{M}}$.

$\bm{M}_{1}=\begin{bmatrix}1&2&\star\\\ 3&\star&\star\\\
\star&\star&5\end{bmatrix},\bm{M}_{2}=\begin{bmatrix}1&2&\star\\\ 3&4&\star\\\
\star&\star&5\end{bmatrix},\bm{M}_{3}=\begin{bmatrix}1&2&\star\\\ 3&4&\star\\\
6&\star&5\end{bmatrix}.$ (4)

Figs. 2(a-b) compare the learned matrices with the ground truth (GT) solutions
having the smallest nuclear norm and rank. For disconnected $\bm{M}_{1}$ (blue
bars), the learned solution achieves neither the smallest nuclear norm nor
rank. For disconnected $\bm{M}_{2}$ with complete bipartite components (green
bars), the learned matrix has the smallest nuclear norm but not rank. For
connected $\bm{M}_{3}$ (red bars), the lowest rank-2 solution is not unique;
the model identifies a particular lowest rank-2 solution, but it does not
correspond to the one with the minimum nuclear norm.

(a) Nuclear norm (b) Rank (c) $3\times 3$ matrix (d) $10\times 10$ matrix

Figure 2: (a) Nuclear norms of the learned solutions for $\bm{M}_{1}$,
$\bm{M}_{2}$, and $\bm{M}_{3}$. Dashed lines represent theoretically computed
smallest nuclear norms. (b) Singular values of the learned matrices for
$\bm{M}_{1}$, $\bm{M}_{2}$, and $\bm{M}_{3}$, with the solid line partitioning
significantly nonzero singular values (experimental rank) and the text (GT)
showing the ground truth minimum rank. Mean and standard deviation are
recorded over 100 repetitions. (c) All equivalent sampling patterns of the
$3\times 3$ matrix completion problem (see Appendix B for details). Cyan stars
marked the case learning the lowest-rank solution. (d) Reconstruction error of
the solutions for a $10\times 10$ matrix reconstruction problem with
$\bm{M}^{*}$ randomly sampled at rank $r=1$ and sample size set to the minimum
reconstruction setting $n=2rd-r^{2}$.

To thoroughly study all possible cases, we examine all sampling patterns of
the $3\times 3$ matrix completion. Fig. 2(c) shows that the model consistently
learns the lowest-rank solution for connected sampling patterns but fails to
do so for disconnected patterns. Fig. 2(d) further verifies the impact of
connectivity on low-rank matrix recovery by comparing the reconstruction error
for 100 randomly sampled rank-1 matrices using two connected sampling patterns
(red and blue dots) and one disconnected sampling pattern (green dots). The
model consistently achieves small reconstruction errors under connected
sampling patterns, while the error is significantly larger for the
disconnected pattern.

These empirical results demonstrate an implicit preference for low rank
induced by connectivity and a preference for low nuclear norm in a particular
kind of disconnection. In the following section, we will investigate the
training dynamics under both connected and disconnected scenarios.

##  5 Training dynamics in connected and disconnected cases

###  5.1 Connected case

This section empirically demonstrates the detailed dynamics of connected
observed data. Fig. 3(a) shows the connected target matrix $\bm{M}$ with a
single unknown element denoted by $\star$. The rank of $\bm{M}$ is at least
three and equals three if and only if $\star=1.2$.

(a) Matrix completion (b) Initialization scale (c) Training loss (d) Gradient
norm

(e) Singular values of $\bm{W}$

(f) Singular values of $\bm{A}$

(g) Singular values of $\bm{B}$

(h) Singular values of $\bm{W}_{\text{aug}}$

Figure 3: (a) The matrix $\bm{M}$ to be completed, with the $\star$ position
unknown. (b) The four singular values of the learned solution at different
initialization scale (Gaussian distribution, mean 0, variance from $10^{0}$ to
$10^{-16}$). (c) The training loss curves corresponding to different sampling
locations. In each sampling, we randomly cover an element and observe the
other 15 elements. (d) Evolution of the $l_{2}$-norm of the gradients
throughout the training process. The cyan crosses represent the difference
between the matrix corresponding to the saddle point and the optimal
approximation at each rank. (e-h) Evolution of singular values for matrices
$\bm{W},\bm{A},\bm{B}$, and $\bm{W}_{\mathrm{aug}}$ during training.

##### Learning lowest-rank solution.

We initialize $\bm{A}$ and $\bm{B}$ with different scales and record the
singular values of the learned matrix. As depicted in Fig. 3(b), when starting
with larger initialization, the learned solutions are almost always rank-4.
Conversely, as the initialization scale decreases, the first three singular
values of the learned solution are consistently maintained in magnitude, but
the fourth singular value keeps decreasing, resulting in the model learning
the lowest rank-3 solution.

##### Traversing progressive optima at each rank.

For a small initialization, the loss curves exhibit a steady, stepwise decline
(Fig. 3(c)). The flat periods correspond to small gradient norms, indicating
potential saddle points (Fig. 3(d)). We compare the matrices learned at these
saddle points with the optimal approximation of each rank and plot their
difference in Fig. 3(d), which is very small. These findings suggest that the
model starts near $\bm{0}$ (rank-0) and progressively finds optimal
approximations within rank-1, rank-2, and higher-rank manifolds until reaching
a global minimum.

##### Alignment of the row space of $\bm{A}$ and the column space of $\bm{B}$.

Starting with small initialization, we track the rank (number of significantly
non-zero singular values) of $\bm{W}=\bm{A}\bm{B}$, $\bm{A}$, $\bm{B}$, and
the augmented matrix $\bm{W}_{\mathrm{aug}}$ during the training process. We
observe that the rank gradually increases, with singular values growing
rapidly one after another (Fig. 3(e-h)). Throughout the entire process, we
consistently find that
$\textnormal{rank}(\bm{A})=\textnormal{rank}(\bm{B}^{\top})=\textnormal{rank}(\bm{W}_{\mathrm{aug}})$,
which implies that the row space of $\bm{A}$ and the column space of $\bm{B}$
remain aligned at all times. This alignment corresponds to a special structure
that we refer to as the “Hierarchical Intrinsic Invariant Manifold” in Sec.
6.1, which plays a crucial role in the overall dynamics of the system.

The dynamics of increasing ranks step by step aligns with the description of
Greedy Low Rank Learning (GLRL) (Li et al., 2020). However, we will show next
that when the observed data are disconnected, the learning process is not
equivalent to GLRL.

###  5.2 Disconnected case

In this section, we present a typical experiment in the disconnected
situation. As depicted in Fig. 4(a), the target matrix $\bm{M}$ contains four
unknown elements denoted by $\star$ and is disconnected. The rank of $\bm{M}$
is at least one, and there are infinitely many rank-1 solutions.

(a) Matrix completion (b) Singular values of $\bm{A}$

(c) Singular values of $\bm{B}$

(d) Singular values of $\bm{W}_{\mathrm{aug}}$

(e) Training loss (f) Output at $\star$ (g) Saddle point

(h) GLRL solution

Figure 4:  (a) The matrix to be completed, with unknown entries marked by
$\star$. (b-d) Evolution of singular values for $\bm{A}$, $\bm{B}$, and
$\bm{W}_{\mathrm{aug}}$ during training. (e) Training loss for different
disconnected sampling patterns, each covering 4 elements and observing the
remaining 5. (f) Learned values at symmetric positions $(0,1)$ and $(1,0)$
under varying initialization scales (zero mean, varying variance). Each point
represents one of ten random experiments per variance; labels show
initialization variance. Other symmetric positions exhibit similar behavior.
(g) Learned output at the saddle point corresponding to the red dot in (e).
(h) Final learned solution of the GLRL algorithm (Li et al., 2020).

##### Alignment of the row space of $\bm{A}$ and the column space of $\bm{B}$.

As shown in Fig. 4(b-e), the learning process in the disconnected case is
similar to the previous experiment: the model naturally evolves from low-rank
to high-rank, with each step increasing a singular value and satisfying
$\textnormal{rank}(\bm{A})=\textnormal{rank}(\bm{B}^{\top})=\textnormal{rank}(\bm{W}_{\mathrm{aug}})$.
Fig. 4(f) illustrates that as the initialization scale decreases, the model
tends to learn symmetric solutions. However, unlike the connected case, the
output does not approach a particular solution as the initialization
decreases. For this specific disconnected $\bm{M}$, we will show that every
symmetric solution learned is a minimal nuclear norm solution(see Sec. 6.2
Thm. 4). For fewer observations, the experimental phenomena are similar (see
Appendix B Fig. B5).

##### Lowest-rank solution is not learned.

Despite the adaptive learning behavior, the final learned solution has rank 2,
as evidenced by the two significantly non-zero singular values in Fig. 4(b-d).
Examining the dynamics (3), we find that they decouple into two independent
systems: one for the 1st and 3rd rows of $\bm{A}$ and columns of $\bm{B}$, and
another for the 2nd row of $\bm{A}$ and column of $\bm{B}$. Fig. 4(g) shows
that the model first learns the surrounding elements ${1,3,3,9}$ (rank-1
saddle point), then learns the middle element ${5}$ in the next stage. The
decoupling of dynamics is equivalent to the definition of disconnection (see
Appendix A Prop. A.4 for proof). The loss curves in Fig. 4(e) consistently
indicate that in the disconnected case, the model learns a sub-optimal
solution in the rank-1 manifold, resulting in a rank-2 solution finally.

##### Not equivalent to GLRL in disconnected case.

We compare the GLRL algorithm (Li et al., 2020) with the matrix factorization
model for solving the same matrix completion problem (Fig. 4). Li et al.
(2020) claim that the matrix factorization dynamics is mathematically
equivalent to the GLRL algorithm under reasonable assumptions. While GLRL
learns the same rank-1 saddle point shown in Fig. 4(g) in the first stage, it
then fills unobserved elements with 0, resulting in a unique rank-2 solution
(Fig. 4(h)). In contrast, the matrix factorization model learns symmetric
solutions with some degree of freedom depending on the random seed (Fig.
4(f)). The key difference is that the first critical point (Fig. 4(g)) reached
by the trajectory is a sub-optimal and not a second-order stationary point of
the rank-1 manifold as assumed by Li et al. (2020). Therefore, the equivalence
assumption between GLRL and matrix factorization does not hold in the
disconnected case.

##  6 Theoretical analysis of training dynamics and implicit regularization

###  6.1 Characterization of training dynamics

Matrix factorization models exhibit a distinctive adaptive learning behavior,
progressively evolving from low rank to high rank. Understanding this
phenomenon is rooted in grasping the global dynamics of matrix factorization
models, where the role of intrinsic invariant manifolds becomes critical.

######  Proposition 1 (Hierarchical Intrinsic Invariant Manifold (HIIM)).

(see Appendix A Prop. A.1 for Proof) Let $\bm{f}_{\bm{\theta}}=\bm{A}\bm{B}$
be a matrix factorization model and
$\\{\bm{\alpha}_{1},\cdots,\bm{\alpha}_{k}\\}$ be $k$ linearly independent
vectors. Define the manifold $\bm{\Omega}_{k}$ as
$\bm{\Omega}_{k}:=\bm{\Omega}_{k}(\bm{\alpha}_{1},\cdots,\bm{\alpha}_{k})=\\{\bm{\theta}=(\bm{A},\bm{B})\mid\textnormal{row}(\bm{A})=\textnormal{col}(\bm{B})=\textnormal{span}\\{\bm{\alpha}_{1},\cdots,\bm{\alpha}_{k}\\}\\}$.
The manifold $\bm{\Omega}_{k}$ possesses the following properties:

(i) Invariance under Gradient Flow: Given data $S$ and the gradient flow
dynamics $\dot{\bm{\theta}}=-\nabla R_{S}(\bm{\theta})$, if the initial point
$\bm{\theta}_{0}\in\bm{\Omega}_{k}$, then $\bm{\theta}(t)\in\bm{\Omega}_{k}$
for all $t\geq 0$.

(ii) Intrinsic Property: $\bm{\Omega}_{k}$ is a data-independent invariant
manifold, meaning that for any data $S$, $\bm{\Omega}_{k}$ remains invariant
under the gradient flow dynamics.

(iii) Hierarchical Structure: The manifolds $\bm{\Omega}_{k}$ form a
hierarchy:
$\bm{\Omega}_{0}\subsetneqq\bm{\Omega}_{1}\subsetneqq\cdots\subsetneqq\bm{\Omega}_{k-1}\subsetneqq\bm{\Omega}_{k}$.

Figs. 3(f-h) and Figs. 4(b-d) show that the training process with small
initialization consistently satisfies
$\textnormal{rank}(\bm{A})=\textnormal{rank}(\bm{B}^{\top})=\textnormal{rank}(\bm{W}_{\mathrm{aug}})$,
aligning with the $\bm{\Omega}_{k}$ invariant manifold. Since a non-zero
initialization in practice, the training trajectory is close to the
$\bm{\Omega}_{k}$ invariant manifold, approaches a critical point, and
transitions to the next level invariant manifold without getting trapped.

In both connected and disconnected scenarios, we observe a step-by-step
hierarchical $\bm{\Omega}_{k}$ invariant manifold traversal. In the connected
case, at each level model reaches an optimal solution. However, in the
disconnected case, we can prove that each connected component induces a
sub-$\bm{\Omega}_{k}$ invariant manifold, leading to the experimentally
observed sub-optimal solution (see Fig. 4).

######  Proposition 2 (Intrinsic Sub-$\bm{\Omega}_{k}$ Invariant Manifold).

(see Appendix A Prop. A.2 for Proof) Let $\bm{f}_{\bm{\theta}}=\bm{A}\bm{B}$
be a matrix factorization model, $\bm{M}$ be an incomplete matrix and
$\bm{\Omega}_{k}$ be an invariant manifold defined in Prop. 1. If $\bm{M}$ is
disconnected with $m$ connected components, then there exist $m$
sub-$\bm{\Omega}_{k}$ manifolds $\bm{\omega}_{k}$ such that
$\bm{\omega}_{k}\subsetneqq\bm{\Omega}_{k}$, each possessing the following
properties:

(i) Invariance under Gradient Flow: Given data $S$ and the gradient flow
dynamics $\dot{\bm{\theta}}=-\nabla R_{S}(\bm{\theta})$, if the initial point
$\bm{\theta}_{0}\in\bm{\omega}_{k}$, then $\bm{\theta}(t)\in\bm{\omega}_{k}$
for all $t\geq 0$.

(ii) Intrinsic Property: $\bm{\omega}_{k}$ is a data-value-independent
invariant manifold, meaning that for a fixed sampling pattern in $\bm{M}$ and
any observed values $S$, $\bm{\omega}_{k}$ remains invariant under the
gradient flow.

(iii) Strict Subset Relation: The output set
$\\{\bm{f}_{\bm{\theta}}\mid\bm{\theta}\in\bm{\omega}_{k}\\}$ is a proper
subset of $\\{\bm{f}_{\bm{\theta}}\mid\bm{\theta}\in\bm{\Omega}_{k}\\}$,
namely,
$\\{\bm{f}_{\bm{\theta}}\mid\bm{\theta}\in\bm{\omega}_{k}\\}\subsetneqq\\{\bm{f}_{\bm{\theta}}\mid\bm{\theta}\in\bm{\Omega}_{k}\\}$.

(a) Illustration of training trajectories in disconnected case (b) Alignment
of $\textnormal{row}(\bm{A})$ and $\textnormal{col}(\bm{B})$

Figure 5: (a) Illustrated trajectories for the experiment in Fig. 4. The blue
line represents the trajectory converging to the lowest-rank solution, and the
red line represents the actual trajectory experienced by the model. (b) The
parameter trajectory escaping from a second-order stationary point to reach
the next critical point for the experiment in Fig. 3. The 8 scatter points
represent the 4 row vectors of matrix $\bm{A}$ and the 4 column vectors of
matrix $\bm{B}$. For ease of visualization, we randomly project them onto two
dimensions and plot them in polar coordinates.

Fig. 5(a) illustrates the trajectory of the experiment in Fig. 4. In the
disconnected case, sub-$\bm{\Omega}_{k}$ invariant manifolds exist and attract
the dynamics, leading the model to learn sub-optimal solutions on the entire
$\bm{\Omega}_{k}$ invariant manifold. In fact, we can prove that these sub-
optimal solutions are necessarily strict saddle points. This loss landscape
result extends Theorem 5.10 from Li et al. (2020), which established the
findings for the specific case of symmetric matrix factorization models (see
Appendix A Sec. A.3 for a detailed discussion).

######  Theorem 1 (Loss Landscape).

(see Appendix A Thm. A.3 for Proof) Given any data $S$, the critical points of
$R_{S}(\bm{\theta})$ are either strict saddle points or global minima.

Gradient descent easily escapes saddle points (Lee et al., 2016, 2019). Fig.
5(b) shows that when the model escapes a saddle point, the parameters
initially appear chaotic but align in one direction after some time,
consistent with the “condensation” phenomenon in neural networks (Luo et al.,
2021; Zhou et al., 2022). For matrix factorization models, by meticulously
analyzing the Hessian matrix structure (see Appendix A.5), we find that this
alignment corresponds to an $\bm{\Omega}_{1}$ invariant manifold, resulting in
a rank increase of one at a time. Under reasonable assumptions, we prove that
the training trajectory follows the $\bm{\Omega}_{k}$ invariant manifold step
by step.

######  Assumption 1 (Top Eigenvalue).

Let $\delta\bm{M}=(\bm{A}_{c}\bm{B}_{c}-\bm{M})_{S_{\bm{x}}}$ be the residual
matrix at the critical point $\bm{\theta}_{c}=(\bm{A}_{c},\bm{B}_{c})$. Assume
that the top eigenvalue of the matrix $\begin{bmatrix}\bm{0}&-\delta\bm{M}\\\
-\delta\bm{M}^{\top}&\bm{0}\end{bmatrix}$ is unique.

######  Assumption 2 (Second-order Stationary Point).

Let $\bm{\Omega}$ be an $\bm{\Omega}_{k}$ invariant manifold or
sub-$\bm{\Omega}_{k}$ invariant manifold defined in Prop. 1 or 2. Assume
$\bm{\theta}_{c}$ is a second-order stationary point within $\bm{\Omega}$,
i.e., $\nabla R_{S}(\bm{\theta}_{c})=0$ and
$\bm{\theta}^{\top}\nabla^{2}R_{S}(\bm{\theta}_{c})\bm{\theta}\geq 0$ for all
$\bm{\theta}\in\bm{\Omega}$.

######  Theorem 2 (Transition to the Next Rank-level Invariant Manifold).

(see Appendix A Thm. A.4 for proof) Consider the dynamics
$\dot{\bm{\theta}}=-\nabla R_{S}(\bm{\theta})$. Let
$\varphi(\bm{\theta}_{0},t)$ denote the value of $\bm{\theta}(t)$ when
$\bm{\theta}(0)=\bm{\theta}_{0}$. Let $\bm{\Omega}$ be an $\bm{\Omega}_{k}$ or
sub-$\bm{\Omega}_{k}$ invariant manifold. Let $\bm{\theta}_{c}\in\bm{\Omega}$
be a critical point satisfying Assump. 1 and 2. Then, for randomly selected
$\bm{\theta}_{0}$, with probability 1 with respect to $\bm{\theta}_{0}$, the
limit

$\tilde{\varphi}(\bm{\theta}_{c},t):=\lim_{\alpha\to
0}\varphi\left(\bm{\theta}_{c}+\alpha\bm{\theta}_{0},t+\frac{1}{\lambda_{1}}\log\frac{1}{\alpha}\right)$
(5)

exists and falls into an invariant manifold $\bm{\Omega}_{k+1}$. Here
$\lambda_{1}$ is the top eigenvalue of $-\nabla^{2}R_{S}(\bm{\theta}_{c})$.

######  Remark.

Assump. 1 ensures that upon departing from a critical point $\bm{\theta}_{c}$,
the trajectory is constrained to escape along a single dominant
eigendirection. Assump. 2 ensures that this escape route aligns with an
$\bm{\Omega}_{1}$ invariant manifold, leading the trajectory into the
$\bm{\Omega}_{k+1}$ invariant manifold. Our Assump. 2 is more general than the
conditions posited by Li et al. (2020), as it remains valid across both
connected and disconnected configurations, and is consistent with our
empirical findings (Figs. 3 and 4). Furthermore, our framework diverges from
Li et al. (2020) by emphasizing the role of hierarchical intrinsic invariant
manifolds in directing the training trajectories, a principle potentially
extendable to a wider spectrum of models.

Thm. 2 provides a characterization of the escape trajectory. It shows that as
the point approaches a second-order stationary point
$\bm{\theta}_{c}\in\bm{\Omega}_{k}$, the trajectory generically converges to a
well-defined limit within $\bm{\Omega}_{k+1}$. Since the origin $\bm{0}$ is
always a second-order stationary point of $\bm{\Omega}_{0}$, the theorem
implies that the trajectory escaping from a small initialization will be close
to $\bm{\Omega}_{1}$. This iterative process gives rise to the phenomenon of
Hierarchical Invariant Manifold Traversal (HIMT), which involves a sequential
progression through these $\bm{\Omega}_{k}$ manifolds.

###  6.2 Implicit regularization analysis

Rank minimization is a challenging non-convex optimization problem. Li et al.
(2018); Jin et al. (2023) proved that the Restricted Isometry Property (RIP)
condition ensures a minimal rank solution. However, the RIP condition is often
too stringent for practical matrix completion. For instance, the matrix
$\bm{M}_{3}$ in Eq. (4) does not satisfy the RIP criteria, yet the model still
finds the minimum rank solution. Our empirical findings (Figs. 1, 2, 3)
suggest that a more lenient condition, specifically the connectivity of the
observed data, frequently leads to convergence towards the minimal rank
solution. Proving this result directly, however, would necessitate a
comprehensive examination of the convergence characteristics within each
$\bm{\Omega}_{k}$ invariant manifold, which is an endeavor we leave for future
work. Despite this, our insights into the system’s dynamics, i.e.,
hierarchical invariant manifold traversal, allow us to assert that if a
trajectory successfully navigates through the optimal on each rank-level
invariant manifold $\bm{\Omega}_{k}$, a solution of minimal rank can be
achieved naturally.

######  Theorem 3 (Minimum Rank).

(see Appendix A Thm. 36 for proof) Consider the dynamics
$\dot{\bm{\theta}}=-\nabla R_{S}(\bm{\theta})$, where
$\bm{\theta}(t)=(\bm{A}(t),\bm{B}(t))$, and denote
$\bm{W}_{t}=\bm{A}(t)\bm{B}(t)$. Assume $\bm{W}_{t}$ achieves an optimal
within each invariant manifold $\bm{\Omega}_{k}$. For a full rank
initialization $\bm{W}_{0}$, if the limit
$\widehat{\bm{W}}=\lim_{\alpha\rightarrow 0}\bm{W}_{\infty}(\alpha\bm{W}_{0})$
exists and is a global optimum with $\widehat{\bm{W}}_{ij}=\bm{M}_{ij}$ for
all $(i,j)\in S_{\bm{x}}$, then

$\widehat{\bm{W}}\in\operatorname{argmin}_{\bm{W}}\operatorname{rank}(\bm{W})\quad\text{s.t.}\quad\bm{W}_{ij}=\bm{M}_{ij},\forall(i,j)\in
S_{\bm{x}}.$ (6)

For a disconnected matrix $\bm{M}$, our theoretical results (Prop. 2) and
experiments (Fig. 4) confirm the existence of sub-$\bm{\Omega}_{k}$ invariant
manifolds. These manifolds attract the training trajectory, leading to sub-
optimal solutions and preventing convergence to the lowest-rank solution.

However, in a specific disconnected case, such as disconnection with complete
bipartite components, as illustrated in Figs. 1 and 2, the minimum nuclear
norm may still serve as a characterization. Gunasekar et al. (2017) proved a
special case: if the observations are commutative, then the symmetric model
will learn the minimum nuclear norm solution. Intriguingly, for the example
$\bm{M}_{2}$ in Eq. (4), even though the observations are not commutative, the
model still learns a minimum nuclear norm solution. In fact, we can prove the
following result, which aligns well with practical experiments.

######  Theorem 4 (Minimum Nuclear Norm Guarantee).

(see Appendix A Thm. 38 for proof) Consider the dynamics
$\dot{\bm{\theta}}=-\nabla R_{S}(\bm{\theta})$, where
$\bm{\theta}(t)=(\bm{A}(t),\bm{B}(t))$, and let
$\bm{W}_{t}=\bm{A}(t)\bm{B}(t)$. If the observation graph associated with the
incomplete matrix $\bm{M}$ is disconnected with complete bipartite components,
and if for a full rank initialization $\bm{W}_{0}$, the limit
$\widehat{\bm{W}}=\lim_{\alpha\rightarrow 0}\bm{W}_{\infty}(\alpha\bm{W}_{0})$
exists and is a global optimum with $\widehat{\bm{W}}_{ij}=\bm{M}_{ij}$ for
all $(i,j)\in S_{\bm{x}}$, then

$\widehat{\bm{W}}\in\operatorname{argmin}_{\bm{W}}\|\bm{W}\|_{*}\quad\text{s.t.}\quad\bm{W}_{ij}=\bm{M}_{ij},\forall(i,j)\in
S_{\bm{x}}.$ (7)

##  7 Conclusion and future work

This study presents a comprehensive experimental and theoretical investigation
of matrix factorization models. The primary objective was to develop a
cohesive framework for understanding the conditions, mechanisms, and reasons
behind the diverse implicit regularization effects exhibited by matrix
factorization models. A key finding of this research is the pivotal role of
the connectivity of observed data in shaping the implicit regularization
behavior. To elucidate this phenomenon, we identified the significance of
hierarchical invariant manifold traversal within the training dynamics.

Our experiments (Figs. 1, 2, 3) provide strong evidence that connected
observed data leads to minimum-rank solutions, as the model learns the optimal
of the $\bm{\Omega}_{k}$ invariant manifold. However, further investigation is
needed to uncover the underlying mechanisms by which connectivity facilitates
optimal attainment across different $\bm{\Omega}_{k}$ invariant manifolds.
Additionally, the trade-offs between initialization scale and training
efficiency warrant further research, as certain cases may require extremely
small initialization, potentially impacting training speed (see Appendix B
Sec. B.4).

Generalizing the insights gained from matrix factorization models to other
architectures is also an important avenue for future work. Our preliminary
experiments indicate that the learning phenomenon from low rank to high rank
persists in deep multi-layer matrix factorization and the query-key
factorization model in Transformer attention mechanisms (see Appendix B Figs.
B7, B9). These findings suggest that the hierarchical invariant manifold
traversal process uncovered in our study may have broader implications and
merit further exploration.

## Acknowledgments and Disclosure of Funding

This work is sponsored by the National Key R&D Program of China Grant No.
2022YFA1008200, the National Natural Science Foundation of China Grant No.
12101402, the Lingang Laboratory Grant No. LG-QS-202202-08, Shanghai Municipal
of Science and Technology Major Project No. 2021SHZDZX0102.

## References

  * Zhang et al. (2017) C. Zhang, S. Bengio, M. Hardt, B. Recht, O. Vinyals,  Understanding deep learning requires rethinking generalization. iclr 2017,  arXiv preprint arXiv:1611.03530 (2017).
  * Zhang et al. (2021) C. Zhang, S. Bengio, M. Hardt, B. Recht, O. Vinyals,  Understanding deep learning (still) requires rethinking generalization,  Communications of the ACM 64 (2021) 107–115.
  * Jiang et al. (2019) Y. Jiang, B. Neyshabur, H. Mobahi, D. Krishnan, S. Bengio,  Fantastic generalization measures and where to find them,  arXiv preprint arXiv:1912.02178 (2019).
  * Zhang et al. (2023) Y. Zhang, Z. Zhang, L. Zhang, Z. Bai, T. Luo, Z.-Q. J. Xu,  Optimistic estimate uncovers the potential of nonlinear models,  arXiv preprint arXiv:2307.08921 (2023).
  * Zhang et al. (2022) Y. Zhang, Z. Zhang, L. Zhang, Z. Bai, T. Luo, Z.-Q. J. Xu,  Linear stability hypothesis and rank stratification for nonlinear models,  arXiv preprint arXiv:2211.11623 (2022).
  * Gunasekar et al. (2017) S. Gunasekar, B. E. Woodworth, S. Bhojanapalli, B. Neyshabur, N. Srebro,  Implicit regularization in matrix factorization,  Advances in neural information processing systems 30 (2017).
  * Arora et al. (2019) S. Arora, N. Cohen, W. Hu, Y. Luo,  Implicit regularization in deep matrix factorization,  Advances in Neural Information Processing Systems 32 (2019).
  * Razin and Cohen (2020) N. Razin, N. Cohen,  Implicit regularization in deep learning may not be explainable by norms,  Advances in neural information processing systems 33 (2020) 21174–21187.
  * Ji and Telgarsky (2019) Z. Ji, M. Telgarsky,  Gradient descent aligns the layers of deep linear networks,  in: 7th International Conference on Learning Representations, ICLR 2019, 2019.
  * Gunasekar et al. (2018) S. Gunasekar, J. D. Lee, D. Soudry, N. Srebro,  Implicit bias of gradient descent on linear convolutional networks,  Advances in neural information processing systems 31 (2018).
  * Gidel et al. (2019) G. Gidel, F. Bach, S. Lacoste-Julien,  Implicit regularization of discrete gradient dynamics in linear neural networks,  Advances in Neural Information Processing Systems 32 (2019).
  * Gissin et al. (2019) D. Gissin, S. Shalev-Shwartz, A. Daniely,  The implicit bias of depth: How incremental learning drives generalization,  in: International Conference on Learning Representations, 2019.
  * Belabbas (2020) M. A. Belabbas,  On implicit regularization: Morse functions and applications to matrix factorization,  arXiv preprint arXiv:2001.04264 (2020).
  * Zhou et al. (2022) H. Zhou, Z. Qixuan, T. Luo, Y. Zhang, Z.-Q. Xu,  Towards understanding the condensation of neural networks at initial training,  Advances in Neural Information Processing Systems 35 (2022) 2184–2196.
  * Zhang et al. (2021) Y. Zhang, Z. Zhang, T. Luo, Z. J. Xu,  Embedding principle of loss landscape of deep neural networks,  Advances in Neural Information Processing Systems 34 (2021) 14848–14859.
  * Zhang et al. (2022) Y. Zhang, Y. Li, Z. Zhang, T. Luo, Z. J. Xu,  Embedding principle: A hierarchical structure of loss landscape of deep neural networks,  Journal of Machine Learning 1 (2022) 60–113.
  * Bai et al. (2022) Z. Bai, T. Luo, Z.-Q. J. Xu, Y. Zhang,  Embedding principle in depth for the loss landscape analysis of deep neural networks,  arXiv preprint arXiv:2205.13283 (2022).
  * Fukumizu et al. (2019) K. Fukumizu, S. Yamaguchi, Y.-i. Mototake, M. Tanaka,  Semi-flat minima and saddle points by embedding neural networks to overparameterization,  Advances in Neural Information Processing Systems 32 (2019) 13868–13876.
  * Simsek et al. (2021) B. Simsek, F. Ged, A. Jacot, F. Spadaro, C. Hongler, W. Gerstner, J. Brea,  Geometry of the loss landscape in overparameterized neural networks: Symmetries and invariances,  in: Proceedings of the 38th International Conference on Machine Learning, PMLR, 2021, pp. 9722–9732.
  * Li et al. (2020) Z. Li, Y. Luo, K. Lyu,  Towards resolving the implicit bias of gradient descent for matrix factorization: Greedy low-rank learning,  in: International Conference on Learning Representations, 2020.
  * Jin et al. (2023) J. Jin, Z. Li, K. Lyu, S. S. Du, J. D. Lee,  Understanding incremental learning of gradient descent: A fine-grained analysis of matrix sensing,  in: International Conference on Machine Learning, PMLR, 2023, pp. 15200–15238.
  * Li et al. (2018) Y. Li, T. Ma, H. Zhang,  Algorithmic regularization in over-parameterized matrix sensing and neural networks with quadratic activations,  in: Conference On Learning Theory, PMLR, 2018, pp. 2–47.
  * Luo et al. (2021) T. Luo, Z.-Q. J. Xu, Z. Ma, Y. Zhang,  Phase diagram for two-layer relu neural networks at infinite-width limit,  Journal of Machine Learning Research 22 (2021) 1–47.
  * Chizat et al. (2019) L. Chizat, E. Oyallon, F. Bach,  On lazy training in differentiable programming,  Advances in neural information processing systems 32 (2019).
  * Lee et al. (2019) J. D. Lee, I. Panageas, G. Piliouras, M. Simchowitz, M. I. Jordan, B. Recht,  First-order methods almost always avoid strict saddle points,  Mathematical programming 176 (2019) 311–337.
  * Lee et al. (2016) J. D. Lee, M. Simchowitz, M. I. Jordan, B. Recht,  Gradient descent only converges to minimizers,  in: Conference on learning theory, PMLR, 2016, pp. 1246–1257.
  * Stöger and Soltanolkotabi (2021) D. Stöger, M. Soltanolkotabi,  Small random initialization is akin to spectral learning: Optimization and generalization guarantees for overparameterized low-rank matrix reconstruction,  Advances in Neural Information Processing Systems 34 (2021) 23831–23843.
  * Du and Lee (2018) S. Du, J. Lee,  On the power of over-parametrization in neural networks with quadratic activation,  in: International conference on machine learning, PMLR, 2018, pp. 1329–1338.
  * Jacot et al. (2021) A. Jacot, F. Ged, B. Şimşek, C. Hongler, F. Gabriel,  Saddle-to-saddle dynamics in deep linear networks: Small initialization training, symmetry, and sparsity,  arXiv preprint arXiv:2106.15933 (2021).
  * Jiang et al. (2023) L. Jiang, Y. Chen, L. Ding,  Algorithmic regularization in model-free overparametrized asymmetric matrix factorization,  SIAM Journal on Mathematics of Data Science 5 (2023) 723–744.
  * Jacot et al. (2018) A. Jacot, F. Gabriel, C. Hongler,  Neural tangent kernel: Convergence and generalization in neural networks,  Advances in neural information processing systems 31 (2018).

##  Appendix A Proofs of Theoretical Results

In this section, we give all proofs for our theoretical results mentioned in
the main text.

###  A.1 Hierarchical Intrinsic Invariant Manifold and Sub Invariant Manifold

######  Proposition A.1 (Hierarchical Intrinsic Invariant Manifold (HIIM)).

Let $\bm{f}_{\bm{\theta}}=\bm{A}\bm{B}$ be a matrix factorization model and
$\\{\bm{\alpha}_{1},\cdots,\bm{\alpha}_{k}\\}$ be $k$ linearly independent
vectors. Define the manifold $\bm{\Omega}_{k}$ as
$\bm{\Omega}_{k}:=\bm{\Omega}_{k}(\bm{\alpha}_{1},\cdots,\bm{\alpha}_{k})=\\{\bm{\theta}=(\bm{A},\bm{B})\mid\textnormal{row}(\bm{A})=\textnormal{col}(\bm{B})=\textnormal{span}\\{\bm{\alpha}_{1},\cdots,\bm{\alpha}_{k}\\}\\}$.
The manifold $\bm{\Omega}_{k}$ possesses the following properties:

(i) Invariance under Gradient Flow: Given data $S$ and the gradient flow
dynamics $\dot{\bm{\theta}}=-\nabla R_{S}(\bm{\theta})$, if the initial point
$\bm{\theta}_{0}\in\bm{\Omega}_{k}$, then $\bm{\theta}(t)\in\bm{\Omega}_{k}$
for all $t\geq 0$.

(ii) Intrinsic Property: $\bm{\Omega}_{k}$ is a data-independent invariant
manifold, meaning that for any data $S$, $\bm{\Omega}_{k}$ remains invariant
under the gradient flow dynamics.

(iii) Hierarchical Structure: The manifolds $\bm{\Omega}_{k}$ form a
hierarchy:
$\bm{\Omega}_{0}\subsetneqq\bm{\Omega}_{1}\subsetneqq\cdots\subsetneqq\bm{\Omega}_{k-1}\subsetneqq\bm{\Omega}_{k}$.

###### Proof.

(i) Invariance under Gradient Flow.

By definition,
$\bm{\Omega}_{k}:=\bm{\Omega}_{k}(\bm{\alpha}_{1},\cdots,\bm{\alpha}_{k})=\\{\bm{\theta}=(\bm{A},\bm{B})\mid\textnormal{row}(\bm{A})=\textnormal{col}(\bm{B})=\textnormal{span}\\{\bm{\alpha}_{1},\cdots,\bm{\alpha}_{k}\\}\\}$.
Consider the gradient flow dynamics in (8):

$\left\\{\begin{aligned} \dot{\bm{a}}_{i}&=-\dfrac{2}{n}\sum_{j\in
I_{i}}(\bm{a}_{i}\cdot\bm{b}_{\cdot,j}-\bm{M}_{ij})\bm{b}_{\cdot,j}^{\top},\\\
\dot{\bm{b}}_{\cdot,j}&=-\dfrac{2}{n}\sum_{i\in
I_{j}}(\bm{a}_{i}\cdot\bm{b}_{\cdot,j}-\bm{M}_{ij})\bm{a}_{i}^{\top},\end{aligned}\right.$
(8)

where $I_{i}=\\{j|\exists i:(i,j)\in S_{\bm{x}}\\}$, $I_{j}=\\{i|\exists
j:(i,j)\in S_{\bm{x}}\\}$, $\bm{a}_{i}$ and $\bm{b}_{\cdot,j}$ represent the
$i$-th row and $j$-th column of $\bm{A}$ and $\bm{B}$, respectively.

For any $(i,j)\in S_{\bm{x}}$, the evolution of $\bm{a}_{i}$ is coupled with
$\bm{b}_{\cdot,j}$ for $j\in I_{i}$. The condition
$\textnormal{row}(\bm{A})=\textnormal{col}(\bm{B})=\textnormal{span}\\{\bm{\alpha}_{1},\cdots,\bm{\alpha}_{k}\\}$
ensures the existence of $k$ linearly independent vectors
$\bm{\alpha}_{1},\cdots,\bm{\alpha}_{k}\in\mathbb{R}^{d}$ such that
$\bm{a}_{i},\bm{b}_{\cdot,j}\in\mathrm{span}\\{\bm{\alpha}_{1},\cdots,\bm{\alpha}_{k}\\}$
for all $1\leq i,j\leq d$.

Consequently, if $\bm{a}_{i}$ and $\bm{b}_{\cdot,j}$ are initially in
$\mathrm{span}\\{\bm{\alpha}_{1},\bm{\alpha}_{2},\cdots,\bm{\alpha}_{k}\\}$,
they will continue to evolve within this subspace under the gradient flow
dynamics. Additionally, for $(i,j)\notin S_{\bm{x}}$, the gradients for the
corresponding $\bm{a}_{i}$ and $\bm{b}_{\cdot,j}$ will be zero, provided their
initial values are zero, maintaining this state throughout the evolution.

(ii) Intrinsic Property.

As demonstrated in part (i), $\bm{\Omega}_{k}$ is invariant under gradient
flow dynamics for any dataset $S$, confirming its status as a data-independent
invariant manifold.

(iii) Hierarchical Structure.

Th invariant manifold $\bm{\Omega}_{k}$ encompasses matrices of rank up to
$k$, including those of lower ranks. Consequently, the manifolds exhibit the
following hierarchical nesting:

$\bm{\Omega}_{0}\subsetneqq\bm{\Omega}_{1}\subsetneqq\cdots\subsetneqq\bm{\Omega}_{k-1}\subsetneqq\bm{\Omega}_{k}.$

∎

######  Proposition A.2 (Intrinsic Sub-$\bm{\Omega}_{k}$ Invariant Manifold).

Let $\bm{f}_{\bm{\theta}}=\bm{A}\bm{B}$ be a matrix factorization model,
$\bm{M}$ be an incomplete matrix and $\bm{\Omega}_{k}$ be an invariant
manifold defined in Prop. 1. If $\bm{M}$ is disconnected with $m$ connected
components, then there exist $m$ sub-$\bm{\Omega}_{k}$ manifolds
$\bm{\omega}_{k}$ such that $\bm{\omega}_{k}\subsetneqq\bm{\Omega}_{k}$, each
possessing the following properties:

(i) Invariance under Gradient Flow: Given data $S$ and the gradient flow
dynamics $\dot{\bm{\theta}}=-\nabla R_{S}(\bm{\theta})$, if the initial point
$\bm{\theta}_{0}\in\bm{\omega}_{k}$, then $\bm{\theta}(t)\in\bm{\omega}_{k}$
for all $t\geq 0$.

(ii) Intrinsic Property: $\bm{\omega}_{k}$ is a data-value-independent
invariant manifold, meaning that for a fixed sampling pattern in $\bm{M}$ and
any observed values $S$, $\bm{\omega}_{k}$ remains invariant under the
gradient flow.

(iii) Strict Subset Relation: The output set
$\\{\bm{f}_{\bm{\theta}}\mid\bm{\theta}\in\bm{\omega}_{k}\\}$ is a proper
subset of $\\{\bm{f}_{\bm{\theta}}\mid\bm{\theta}\in\bm{\Omega}_{k}\\}$,
namely,
$\\{\bm{f}_{\bm{\theta}}\mid\bm{\theta}\in\bm{\omega}_{k}\\}\subsetneqq\\{\bm{f}_{\bm{\theta}}\mid\bm{\theta}\in\bm{\Omega}_{k}\\}$.

###### Proof.

Existence.

Let us consider an incomplete matrix $\bm{M}$ whose associated observational
graph is divided into $m$ connected components, denoted by
$L_{1},L_{2},\ldots,L_{m}$. For each component $L_{p}$, we define
$S_{\bm{x}}^{L_{p}}$ as the subset of observed indices within $L_{p}$, where
$1\leq p\leq m$ and $S_{\bm{x}}$ is the set of all observed indices.

For each $L_{p}$, we can identify row indices $R_{p}$ and column indices
$C_{p}$ corresponding to the observed entries in $L_{p}$ as follows:

$R_{p}=\\{i|\exists j:(i,j)\in S_{\bm{x}}^{L_{p}}\\},\quad C_{p}=\\{j|\exists
i:(i,j)\in S_{\bm{x}}^{L_{p}}\\}.$

Here, $R_{p}$ includes the row indices and $C_{p}$ includes the column indices
of the entries observed in $L_{p}$.

Define $\bm{A}^{L_{p}}$ and $\bm{B}^{L_{p}}$ as the submatrices of $\bm{A}$
and $\bm{B}$ corresponding to $R_{p}$ and $C_{p}$, respectively, and let
$\bm{A}^{L_{p}}_{r}$ and $\bm{B}^{L_{p}}_{r}$ be the remaining rows not in
$R_{p}$ and $C_{p}$.

Let
$\bm{\Omega}_{k}:=\bm{\Omega}_{k}(\bm{\alpha}_{1},\cdots,\bm{\alpha}_{k})=\\{\bm{\theta}=(\bm{A},\bm{B})\mid\textnormal{row}(\bm{A})=\textnormal{col}(\bm{B})=\textnormal{span}\\{\bm{\alpha}_{1},\cdots,\bm{\alpha}_{k}\\}\\}$
be the given $\bm{\Omega}_{k}$ invariant manifold.

The sub-$\bm{\Omega}_{k}$ invariant manifold associated with the connected
component $L_{p}$ can be defined as

$\bm{\omega}_{k}^{L_{p}}:=\\{(\bm{\theta}=(\bm{A},\bm{B}))\mid\textnormal{row}(\bm{A}^{L_{p}})=\textnormal{col}((\bm{B}^{L_{p}}))=\textnormal{span}\\{\bm{\alpha}_{1},\cdots,\bm{\alpha}_{k}\\},\bm{A}^{L_{p}}_{r}=\bm{B}^{L_{p}}_{r}=\bm{0}\\}.$
(9)

It is easy to check $\bm{\omega}_{k}^{L_{p}}$ is a proper subset of
$\bm{\Omega}_{k}$.

(i) Invariance under Gradient Flow.

The condition
$\textnormal{row}(\bm{A}^{L_{p}})=\textnormal{col}((\bm{B}^{L_{p}}))=\textnormal{span}\\{\bm{\alpha}_{1},\cdots,\bm{\alpha}_{k}\\}$
along with $\bm{A}^{L_{p}}_{r}=\bm{B}^{L_{p}}_{r}=\bm{0}$ guarantees that
$\bm{a}_{i},\bm{b}_{\cdot,j}\in\mathrm{span}\\{\bm{\alpha}_{1},\ldots,\bm{\alpha}_{k}\\}$
for all $(i,j)\in S^{L_{p}}_{\bm{x}}$, and
$\bm{a}_{i},\bm{b}_{\cdot,j}=\bm{0}$ for all $(i,j)\notin S^{L_{p}}_{\bm{x}}$.

In other words, the sub-$\bm{\Omega}_{k}$ invariant manifold
$\bm{\omega}^{L_{p}}_{k}$ is the set of all pairs $(\bm{A},\bm{B})$ where, for
each observed position $(i,j)$ in the connected component $L_{p}$, the vectors
$\bm{a}_{i}$ and $\bm{b}_{\cdot,j}$ lie within the span of
$\\{\bm{\alpha}_{1},\cdots,\bm{\alpha}_{k}\\}$, and for any position not in
$S^{L_{p}}_{\bm{x}}$, the vectors are zero.

Considering the dynamics expressed in equation (8), it is evident that the
evolution of $\bm{a}_{i}$ is influenced by $\bm{b}_{\cdot,j}$ for $(i,j)\in
S^{L_{p}}_{\bm{x}}$. Hence, if $\bm{a}_{i}$ and $\bm{b}_{\cdot,j}$ are
initially in the span of
$\\{\bm{\alpha}_{1},\bm{\alpha}_{2},\cdots,\bm{\alpha}_{k}\\}$, they will
continue to evolve within this span under the gradient flow dynamics.
Moreover, for positions $(i,j)\notin S^{L_{p}}_{\bm{x}}$, we consider the
following scenarios:

  * •

For $(i,j)\notin S_{\bm{x}}$, since the matrix entry $\bm{M}_{ij}$ does not
contribute to the loss $R_{S}(\bm{\theta})$, the gradients for corresponding
$\bm{a}_{i}$ and $\bm{b}_{\cdot,j}$ will perpetually be zero. Thus, if their
initial values are zero, they will remain zero throughout the evolution.

  * •

For $(i,j)\in S_{\bm{x}}$ but not in $S_{\bm{x}}^{L_{p}}$, the dynamics
corresponding to different connected components are decoupled. Therefore, if
the initial values for $\bm{a}_{i}$ and $\bm{b}_{\cdot,j}$ are zero, they will
stay zero during the evolution.

(ii) Intrinsic Property.

As established in (i), the manifold $\bm{\omega}_{k}^{L_{p}}$ is invariant
under gradient flow for any data $S$ with a fixed sampling pattern, qualifying
it as a data-value-independent invariant manifold.

(iii) Strict Subset Relation.

The output set $\\{\bm{f}_{\bm{\theta}}\mid\bm{\theta}\in\bm{\Omega}_{k}\\}$
encompasses all matrices of rank $k$, whereas
$\\{\bm{f}_{\bm{\theta}}\mid\bm{\theta}\in\bm{\omega}_{k}^{L_{p}}\\}$ is
limited to rank-$k$ matrices with specific row and column indices confined to
$R_{p}$ and $C_{p}$. Consequently,
$\\{\bm{f}_{\bm{\theta}}\mid\bm{\theta}\in\bm{\omega}_{k}\\}$ forms a strict
subset of $\\{\bm{f}_{\bm{\theta}}\mid\bm{\theta}\in\bm{\Omega}_{k}\\}$, as
stated by
$\\{\bm{f}_{\bm{\theta}}\mid\bm{\theta}\in\bm{\omega}_{k}\\}\subsetneqq\\{\bm{f}_{\bm{\theta}}\mid\bm{\theta}\in\bm{\Omega}_{k}\\}$.
∎

###  A.2 Connectivity

######  Definition A.1 (Associated Observation Graph).

Given a incomplete matrix $\bm{M}$ to be completed and its observation matrix
$\bm{P}$, the associated observation graph $G_{\bm{M}}$ is the bipartite graph
with adjacency matrix $\begin{bmatrix}\bm{0}&\bm{P}^{\top}\\\
\bm{P}&\bm{0}\end{bmatrix}$, with isolated vertices removed.

######  Definition A.2 (Connectivity).

A matrix $\bm{M}$ to be completed is considered connected if its associated
observation graph $G_{\bm{M}}$ is connected, otherwise, we call it
disconnected. The connected components of $\bm{M}$ are defined as the
connected components of this graph.

######  Definition A.3 (Disconnectivity with Complete Bipartite Components).

A matrix $\bm{M}$ to be completed is considered disconnected with complete
bipartite components if its associated observation graph $G_{\bm{M}}$ is
disconnected and each connected component forms a complete bipartite subgraph.

######  Remark.

In the bipartite graph representation of the observed data, isolated vertices
correspond to entire rows or columns of the matrix $\bm{M}$ that are not
observed. These rows or columns do not contribute to the loss calculation and
have no influence on the dynamics of the matrix factorization under
infinitesimal initialization. Consequently, when analyzing the connectivity of
the observed data and its impact on the learning dynamics, these isolated
vertices can be safely disregarded.

######  Remark.

The disconnectivity of the bipartite graph representing the observed data is
equivalent to the reducibility of the adjacency matrix
$\begin{bmatrix}\bm{0}&\bm{P}^{\top}\\\ \bm{P}&\bm{0}\end{bmatrix}$, where
$\bm{P}$ is the binary observation matrix indicating the positions of the
observed entries in $\bm{M}$.

In the context of matrix completion problems, such as the Netflix problem,
connectivity has a practical interpretation. Connected components in the
bipartite graph indicate groups of users and movies that are linked by the
users’ viewing history. Users within the same connected component are related
through the movies they have watched in common. Due to this practical
significance, we prefer to use the term “connectivity” instead of
“reducibility” when discussing the structure of the observed data in matrix
completion problems.

######  Definition A.4 (Connectivity of Observed Data).

Given a matrix $\bm{M}$ to be completed, an undirected simple graph $G$ can be
induced from it: the nodes of the graph are the observed elements in the
matrix, and two nodes are adjacent if and only if they are in the same row or
column of the matrix $\bm{M}$. A matrix $\bm{M}$ to be completed is considered
connected if its induced graph $G$ is connected, otherwise, we call it
disconnected.

######  Lemma A.1.

For any simple graph $G$, if we remove all isolated vertices from $G$ to
obtain a new graph $G^{\prime}$, then $G^{\prime}$ is connected if and only if
the line graph of $G^{\prime}$, denoted as $L(G^{\prime})$, is connected.

###### Proof.

$\implies$ Assume $G^{\prime}$ is connected. Consider any two nodes in
$L(G^{\prime})$, which correspond to two edges in $G^{\prime}$, say $e_{1}$
and $e_{2}$. Since $G^{\prime}$ is connected, there exists a path connecting
the endpoints of $e_{1}$ and $e_{2}$. This path corresponds to a sequence of
edges in $G^{\prime}$, which in turn corresponds to a path connecting the
nodes representing $e_{1}$ and $e_{2}$ in $L(G^{\prime})$. Therefore,
$L(G^{\prime})$ is connected.

$\impliedby$ Conversely, assume $L(G^{\prime})$ is connected. Consider any two
vertices $v_{1}$ and $v_{2}$ in $G^{\prime}$. Since $G^{\prime}$ has no
isolated vertices, each of $v_{1}$ and $v_{2}$ is incident to at least one
edge. Let these edges be $e_{1}$ and $e_{2}$, respectively. Since
$L(G^{\prime})$ is connected, there exists a path connecting the nodes
representing $e_{1}$ and $e_{2}$ in $L(G^{\prime})$. This path corresponds to
a sequence of edges in $G^{\prime}$, which in turn corresponds to a path
connecting $v_{1}$ and $v_{2}$ in $G^{\prime}$. Therefore, $G^{\prime}$ is
connected.

In conclusion, we have proven that for any simple graph $G$, if we remove all
isolated vertices from $G$ to obtain a new graph $G^{\prime}$, then
$G^{\prime}$ is connected if and only if the line graph of $G^{\prime}$,
denoted as $L(G^{\prime})$, is connected. ∎

######  Proposition A.3.

Given a incomplete matrix $\bm{M}$, the connectivity of $\bm{M}$ defined in
Def. A.2 and Def. A.4 is equivalent.

###### Proof.

By definition, each edge of a bipartite graph corresponds to an observed data
item, and two edges in a bipartite graph are adjacent if and only if the two
corresponding observed data items are in the same row or column. Therefore,
the connectivity of the observed data is equivalent to the connectivity of the
edges of the bipartite graph, which is, in turn, equivalent to the
connectivity of the line graph of the bipartite graph.

According to Lem. A.1, for any graph $G$, if we remove all isolated vertices
from $G$ to obtain a new graph $G^{\prime}$, then $G^{\prime}$ is connected if
and only if the line graph of $G^{\prime}$, denoted as $L(G^{\prime})$, is
connected.

In the context of the bipartite graph representation of the observed data,
removing isolated vertices corresponds to removing rows and columns that
contain no observed entries. Thus, the connectivity of the bipartite graph
after removing isolated vertices is equivalent to the connectivity of the
observed data as defined in Def. A.4.

Consequently, the connectivity of the observed data as defined in Def. A.2
(based on the line graph of the bipartite graph) is equivalent to the
connectivity of the observed data as defined in Def. A.4 (based on the
connectivity of observed data). ∎

######  Definition A.5 (Decoupling of Dynamics).

Given an incomplete matrix $\bm{M}$, consider the gradient flow dynamics of
matrix factorization models, $\forall 1\leq i,j\leq d$,

$\left\\{\begin{aligned} \dot{\bm{a}}_{i}&=-\dfrac{2}{n}\sum_{j\in
I_{i}}(\bm{a}_{i}\cdot\bm{b}_{\cdot,j}-\bm{M}_{ij})\bm{b}_{\cdot,j}^{\top},\\\
\dot{\bm{b}}_{\cdot,j}&=-\dfrac{2}{n}\sum_{i\in
I_{j}}(\bm{a}_{i}\cdot\bm{b}_{\cdot,j}-\bm{M}_{ij})\bm{a}_{i}^{\top}.\end{aligned}\right.$
(10)

The dynamics are said to be decoupled if there exist disjoint subsets of
indices $R_{1},R_{2},\dots,R_{k}\subseteq\\{1,2,\dots,d\\}$ for the rows of
$\bm{A}$ and $C_{1},C_{2},\dots,C_{k}\subseteq\\{1,2,\dots,d\\}$ for the
columns of $\bm{B}$, such that for each $l\in\\{1,2,\dots,k\\}$, the dynamics
of $\\{\bm{a}_{i}:i\in R_{l}\\}$ and $\\{\bm{b}_{\cdot,j}:j\in C_{l}\\}$ form
an independent system of equations. In other words, the dynamics can be
divided into $k(k>1)$ independent subsystems, each involving a subset of rows
of $\bm{A}$ and a subset of columns of $\bm{B}$. If such a division is not
possible, the dynamics are said to be coupled.

######  Proposition A.4.

Given an incomplete matrix $\bm{M}$, if it is disconnected as defined by Def.
A.2, then the dynamics are decoupled as defined by Def. A.5; if it is
connected as defined by Def. A.2, then the dynamics are coupled as defined by
Def. A.5.

###### Proof.

Consider a matrix $\bm{M}$ to be completed, with its associated observation
graph comprising $m$ connected components, denoted as
$L_{1},L_{2},\cdots,L_{m}$. Let $S_{\bm{x}}^{L_{p}}\subseteq S_{\bm{x}}$
represent the subset of observed indices corresponding to the connected
component $L_{p}$, where $1\leq p\leq m$ and $S_{\bm{x}}$ denotes the complete
set of observed indices. If the incomplete matrix $\bm{M}$ is disconnected,
then for each connected component $L_{p}$, the subset $S_{\bm{x}}^{L_{p}}$ can
be partitioned into two subsets $R_{p}$ and $C_{p}$, $1\leq p\leq m$, such
that

$R_{p}=\\{i|\exists j:(i,j)\in S_{\bm{x}}^{L_{p}}\\},\quad C_{p}=\\{j|\exists
i:(i,j)\in S_{\bm{x}}^{L_{p}}\\}.$ (11)

In other words, $R_{p}$ contains the row indices and $C_{p}$ contains the
column indices of the observed entries in the connected component $L_{p}$.

It can be easily verified that the dynamics are decoupled in this case, as the
subsets $\\{R_{p},C_{p}\\}_{p=1}^{m}$ satisfy the conditions in Def. A.5. Each
connected component $L_{p}$ corresponds to an independent subsystem involving
the rows of $\bm{A}$ indexed by $R_{p}$ and the columns of $\bm{B}$ indexed by
$C_{p}$.

If $\bm{M}$ is connected, then its associated observation graph consists of a
single connected component, and the entire dynamics are coupled. ∎

##### Examples of connectivity and disconnectivity.

Consider three matrices to be completed, each obtained by adding one more
observation to the previous matrix: $\bm{M}_{1}$ (disconnected), $\bm{M}_{2}$
(disconnected with complete bipartite components), and $\bm{M}_{3}$
(connected).

$\bm{M}_{1}=\begin{bmatrix}1&2&\star\\\ 3&\star&\star\\\
\star&\star&5\end{bmatrix},\bm{M}_{2}=\begin{bmatrix}1&2&\star\\\ 3&4&\star\\\
\star&\star&5\end{bmatrix},\bm{M}_{3}=\begin{bmatrix}1&2&\star\\\ 3&4&\star\\\
6&\star&5\end{bmatrix}.$ (12)

Fig. A1 illustrates the associated graphs $G_{\bm{M}}$.

(a) $\bm{M}_{1}$ (b) $\bm{M}_{2}$ (c) $\bm{M}_{3}$

Figure A1:  The associated observation graphs $G_{\bm{M}}$ of the incomplete
matrices $\bm{M}_{1}$, $\bm{M}_{2}$, and $\bm{M}_{3}$ in Eq. 4. $\bm{M}_{1}$
is disconnected, with its associated observation graph consisting of two
connected components. $\bm{M}_{2}$ is also disconnected, but each connected
component of its associated observation graph forms a complete bipartite
subgraph. In contrast, $\bm{M}_{3}$ is connected, and its associated
observation graph consists of a single connected component.

###  A.3 Loss Landscape

In this paper, we focus on the problem of asymmetric matrix factorization.
Previous literature (Gunasekar et al., 2017; Li et al., 2018, 2020; Jin et
al., 2023) has predominantly concentrated on symmetric matrix factorization
problems. Although asymmetric matrix factorization models can be transformed
into symmetric cases, studying symmetric matrix factorization does not
necessarily cover all aspects of the asymmetric scenarios.

Generally, an asymmetric matrix factorization model $\bm{W}=\bm{A}\bm{B}$ can
be transformed into a symmetric situation by setting

$\bm{U}=\begin{bmatrix}\bm{A}\\\
\bm{B}^{\top}\end{bmatrix}\in\mathbb{R}^{2d\times d}.$

We then consider the model $\bm{W}^{\prime}=\bm{U}\bm{U}^{\top}$, which
corresponds to the following matrix completion problem:

$\begin{bmatrix}\bm{A}\bm{A}^{\top}&\bm{A}\bm{B}\\\
\bm{B}^{\top}\bm{A}^{\top}&\bm{B}^{\top}\bm{B}\end{bmatrix}.$

We define the loss as

$\mathcal{L}^{\prime}\left(\begin{bmatrix}\bm{A}&\bm{B}\\\
\bm{C}&\bm{D}\end{bmatrix}\right)=\frac{1}{2}\mathcal{L}(\bm{B})+\frac{1}{2}\mathcal{L}(\bm{C}^{\top}).$

Li et al. (2020) established the following results:

######  Theorem A.1 (Theorem 5.10 in Li et al. (2020)).

Let $f:\mathbb{R}^{d\times d}\rightarrow\mathbb{R}$ be a convex
$\mathcal{C}^{2}$-smooth function. (1). All stationary points of
$\mathcal{L}:\mathbb{R}^{d\times
d}\rightarrow\mathbb{R},\mathcal{L}(\bm{U})=\frac{1}{2}f\left(\bm{U}\bm{U}^{\top}\right)$
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.

The proof of this theorem relies heavily on Theorem A.2 of Du and Lee (2018),
which requires the parameter matrix $\bm{U}\in\mathbb{R}^{d\times k}$ to
satisfy the condition that $k\geq d$. In the case of symmetric matrix
factorization, where $\bm{U}\in\mathbb{R}^{d\times d}$, this condition is
naturally met. However, for asymmetric matrix factorization, where
$\bm{U}=\begin{bmatrix}\bm{A}\\\
\bm{B}^{\top}\end{bmatrix}\in\mathbb{R}^{2d\times d}$, this condition is not
satisfied, and thus the proof of Theorem A.1 is only applicable to the
symmetric case of matrix factorization.

######  Theorem A.2 (Theorem 3.1 in Du and Lee (2018))).

Let $f:\mathbb{R}^{d\times d}\rightarrow\mathbb{R}$ be a $\mathcal{C}^{2}$
convex function. Then $\mathcal{L}:\mathbb{R}^{d\times
k}\rightarrow\mathbb{R},\mathcal{L}(\bm{U})=f\left(\bm{U}\bm{U}^{\top}\right),k\geq
d$ satisfies that (1). Every local minimizer of $\mathcal{L}$ is also a global
minimizer; (2). All saddles are strict. Here saddles denote those stationary
points whose hessian are not positive semi-definite (thus including local
maximizers).

Below we give a direct proof of the loss landscape of an asymmetric matrix
factorization model.

######  Theorem A.3 (Loss Landscape).

For any data $S$, the critical points of $R_{S}(\bm{\theta})$ are either
strict saddle points or global minima.

###### Proof.

We start by recalling the definition of the loss function:

$R_{S}(\bm{\theta})=\mathcal{L}(\bm{A},\bm{B})=\dfrac{1}{2}\|\bm{A}\bm{B}-\bm{M}\|_{S_{\bm{x}}}^{2}=\dfrac{1}{2}\sum\limits_{(i,j)\in
S_{\bm{x}}}((\bm{A}\bm{B})_{ij}-\bm{M}_{ij})^{2}.$

Let $\bm{\theta}=(\bm{A},\bm{B})$ denote a critical point. We define a new
matrix, $\delta\bm{M}$, as the difference between the product of $\bm{A}$ and
$\bm{B}$ and the matrix $\bm{M}$, with this difference being computed only
over the indices in the set $S_{\bm{x}}$. More formally,
$\delta\bm{M}=(\bm{A}\bm{B}-\bm{M})_{S_{\bm{x}}}$, where the elements of
$\delta\bm{M}$ are given by:

  * •

For $(i,j)\in S_{\bm{x}}$, we have
$\delta\bm{M}_{ij}=(\bm{A}\bm{B})_{ij}-\bm{M}_{ij}$.

  * •

For $(i,j)\notin S_{\bm{x}}$, we have $\delta\bm{M}_{ij}=0$.

This definition of $\delta\bm{M}$ ensures that we only consider the
differences in the entries that belong to the set $S_{\bm{x}}$, while all
other entries are set to zero.

Consider the function:

$\displaystyle\mathcal{L}(\bm{A}+\bm{\bm{\varepsilon}},\bm{B}+\bm{\eta})$
$\displaystyle=\dfrac{1}{2}\|\delta\bm{M}+\bm{\varepsilon}\bm{B}+\bm{A}\bm{\eta}+\bm{\varepsilon}\bm{\eta}\|_{S_{\bm{x}}}^{2}$
(13)
$\displaystyle=\frac{1}{2}\|\delta\bm{M}\|_{S_{\bm{x}}}^{2}+\langle\delta\bm{M},\bm{\varepsilon}\bm{B}+\bm{A}\bm{\eta}\rangle_{S_{\bm{x}}}+\frac{1}{2}\|\bm{\varepsilon}\bm{B}+\bm{A}\bm{\eta}\|_{S_{\bm{x}}}^{2}$
$\displaystyle\quad+\langle\delta\bm{M},\bm{\varepsilon}\bm{\eta}\rangle_{S_{\bm{x}}}+o(\|\bm{\varepsilon}\|^{2},\|\bm{\eta}\|^{2}),$

where the inner product of two matrices $\bm{A},\bm{B}$ is defined as
$\left<\bm{A},\bm{B}\right>:=\mathrm{Tr}(\bm{A}\bm{B}^{\top})$.

At the critical point, the first order term
$\langle\delta\bm{M},\bm{\varepsilon}\bm{B}+\bm{A}\bm{\eta}\rangle_{S_{\bm{x}}}$
equals 0. The Hessian operator, representing the second order term, is given
by:

$h_{\bm{A},\bm{B}}(\bm{\varepsilon},\bm{\eta})=\frac{1}{2}\|\bm{\varepsilon}\bm{B}+\bm{A}\bm{\eta}\|_{S_{\bm{x}}}^{2}+\langle\delta\bm{M},\bm{\varepsilon}\bm{\eta}\rangle_{S_{\bm{x}}}.$

Our goal is to demonstrate that if $\delta\bm{M}\neq\bm{0}$, there always
exists $\bm{\varepsilon},\bm{\eta}$ such that
$h_{\bm{A},\bm{B}}(\bm{\varepsilon},\bm{\eta})<0$. To this end, we consider
the ranks of matrices $\bm{A}$ and $\bm{B}$ in two cases:

(i) $\textnormal{rank}(\bm{A})<d$ or $\textnormal{rank}(\bm{B})<d$:

Without loss of generality, we assume
$\delta\bm{M}_{ij}:=\delta\bm{M}_{ij}\neq 0$ for some $(i,j)\in S_{\bm{x}}$,
and $\textnormal{rank}\bm{A}<d$. Under these conditions, there exists a non-
zero vector $\bm{v}$ such that $\bm{A}\bm{v}=0$.

We set $\bm{\eta}^{*}_{,j}=\bm{v}$ and $\bm{\eta}^{*}_{,s}=0$ for $s\neq j$,
where $\bm{\eta}^{*}_{,j}$ denotes the $j$-th column of the matrix
$\bm{\eta}$. Let $\bm{\varepsilon}_{i}^{*}=\bm{w}^{\top}\in\mathbb{R}^{d}$ and
$\bm{\varepsilon}_{s}=0$ for $s\neq i$, where $\bm{\varepsilon}^{*}_{i}$
denotes the $i$-th row of the matrix $\bm{\varepsilon}$.

We then have:

$\displaystyle h_{\bm{A},\bm{B}}(\bm{\varepsilon}^{*},\bm{\eta}^{*})$
$\displaystyle=\frac{1}{2}\|\bm{\varepsilon}^{*}\bm{B}\|_{S_{\bm{x}}}^{2}+\delta\bm{M}_{ij}\bm{w}^{\top}\bm{v}$
$\displaystyle\leq\frac{1}{2}\|\bm{w}^{\top}\bm{B}\|_{S_{\bm{x}}}^{2}+\delta\bm{M}_{ij}\bm{w}^{\top}\bm{v}$
$\displaystyle=\frac{1}{2}\bm{w}^{\top}\bm{B}\bm{B}^{\top}\bm{w}+\delta\bm{M}_{ij}\bm{w}^{\top}\bm{v}.$

We define
$g(\bm{w},\bm{v})=\frac{1}{2}\bm{w}^{\top}\bm{B}\bm{B}^{\top}\bm{w}+\delta\bm{M}_{ij}\bm{w}^{\top}\bm{v}$
and consider:

$\displaystyle g(-\alpha\delta\bm{M}_{ij}\bm{v},\bm{v})$
$\displaystyle=\frac{1}{2}\alpha^{2}\delta\bm{M}_{ij}^{2}\bm{v}^{\top}\bm{B}\bm{B}^{\top}\bm{v}-\alpha\delta\bm{M}_{ij}^{2}\bm{v}^{\top}\bm{v}$
$\displaystyle=\frac{1}{2}\alpha^{2}\delta\bm{M}_{ij}^{2}(\bm{v}^{\top}\bm{B}\bm{B}^{\top}\bm{v}-2\bm{v}^{\top}\bm{v}).$

For $0<\alpha<\frac{2}{\lambda_{\bm{B},\max}}$, where $\lambda_{\bm{B},\max}$
represents the top eigenvalue of $\bm{B}\bm{B}^{\top}$, we find
$g(-\alpha\delta\bm{M}_{ij}\bm{v},\bm{v})<0$.

Therefore, when $\bm{w}=-\alpha\delta\bm{M}_{ij}\bm{v}$, we obtain
$h_{\bm{A},\bm{B}}(\bm{\varepsilon}^{*},\bm{\eta}^{*})<0$. This immediately
implies the critical point $\bm{\theta}=(\bm{A},\bm{B})$ is a strict saddle
point.

(ii) $\textnormal{rank}(\bm{A})=\textnormal{rank}(\bm{B})=d$:

Let $\bm{\varepsilon}=\alpha\delta\bm{M}\bm{B}^{-1}$ and $\bm{\eta}=0$. In
this scenario, the first order term
$\langle\delta\bm{M},\bm{\varepsilon}\bm{B}+\bm{A}\bm{\eta}\rangle_{S_{\bm{x}}}$
in Eq. (13) simplifies to $\alpha\|\delta\bm{M}\|_{S_{\bm{x}}}^{2}$. At a
critical point, this quantity equals zero, it implies that $\delta\bm{M}=0$.
This in turn implies that the critical point $\bm{\theta}=(\bm{A},\bm{B})$ is
a global minimum.

This concludes the proof, establishing that the critical points of
$R_{S}(\bm{\theta})$ are either strict saddle points or global minima. ∎

###  A.4 Escaping from Top Eigendirection

In this section, we focus on the dynamics of escaping from a critical point.
According to Prop. A.3, the loss landscape consists solely of strict saddle
points and a global minimum. Consequently, gradient-based methods can readily
escape from a critical point that is not a global minimum.

In the following, we will demonstrate that the escaping dynamics near a
critical point can be approximated by a linearized version of these dynamics.
For this, consider the following:

$\dot{\bm{\theta}}=-\nabla R_{S}(\bm{\theta}).$ (14)

Assume $\bm{\theta}_{c}$ is a saddle point for which $\nabla
R_{S}(\bm{\theta}_{c})=0$. We can apply a first-order Taylor expansion to the
right-hand side of Eq. (14), yielding:

$-\nabla R_{S}(\bm{\theta})=-\nabla
R_{S}(\bm{\theta}_{c})-\nabla^{2}R_{S}(\bm{\theta}_{c})(\bm{\theta}-\bm{\theta}_{c})+\mathcal{O}(\|\bm{\theta}-\bm{\theta}_{c}\|^{2}),$
(15)

where $\nabla^{2}R_{S}(\bm{\theta}_{c})$ represents the Hessian matrix. Given
that $\nabla R_{S}(\bm{\theta}_{c})=0$, the gradient flow dynamics around
$\bm{\theta}_{c}$ can be approximated as:

$\dot{\bm{\theta}}=\bm{H}(\bm{\theta}-\bm{\theta}_{c}),$ (16)

where $\bm{H}:=-\nabla^{2}R_{S}(\bm{\theta}_{c})$. Eq. (16) is a classic
linear ordinary differential equation, with the solution:

$\bm{\theta}(t)=\mathrm{e}^{t\bm{H}}(\bm{\theta}_{0}-\bm{\theta}_{c})+\bm{\theta}_{c}.$
(17)

The dynamics near a critical point can be approximated by a linearized
version. Hence, in the vicinity of a critical point, we can analyze the
linearized dynamics to understand the escape mechanism. In the following, we
will show that during this escape process, the dynamics follow a pattern
referred to as dominant eigenvalue dynamics.

Eq. (17) elucidates that the dynamics near the critical point
$\bm{\theta}_{c}$ are predominantly dictated by the properties of $\bm{H}$, a
real symmetric matrix in $\mathbb{R}^{2d^{2}\times 2d^{2}}$. Its
eigendecomposition is given by:

$\bm{H}:=-\nabla^{2}R_{S}(\bm{\theta}_{c})=\bm{Q}\bm{\Lambda}\bm{Q}^{\top},$
(18)

where $\bm{\Lambda}$ is a diagonal matrix and $\bm{Q}$ is an orthogonal
matrix. Let $\lambda_{1}>\lambda_{2}>\cdots>\lambda_{s}\in\mathbb{R}$ denote
the eigenvalues of $\bm{H}$, and let
$\bm{q}_{i1},\bm{q}_{i2},\cdots,\bm{q}_{il_{i}}$ represent the eigenvectors
corresponding to $\lambda_{i}$.

Given that $\lambda_{1}>\lambda_{2}$, the ratio
$e^{\lambda_{1}t}/e^{\lambda_{i}t}$ for $i>1$ grows exponentially fast.
Consequently, near $\bm{\theta}_{c}$, the evolution of the system is primarily
driven by the eigenvectors $\bm{q}_{11},\bm{q}_{12},\cdots,\bm{q}_{1l_{1}}$
associated with the largest eigenvalue $\lambda_{1}$.

This following proposition formalizes the intuitive idea that in the vicinity
of a saddle point, the dynamics primarily follow the direction associated with
the largest eigenvalue. This leading eigendirection becomes increasingly
dominant as time evolves, allowing for an escape from the saddle point and
facilitating a specific structured transition.

Let’s consider $\bm{\theta}_{c}$ to be a saddle point. Consider:

$\dot{\bm{\theta}}=-\nabla^{2}R_{S}(\bm{\theta}_{c})(\bm{\theta}-\bm{\theta}_{c}),\quad\bm{\theta}(0)=\bm{\theta}_{0}.$
(19)

Here, $\bm{\theta}_{0}$ and $\bm{\theta}_{c}$ are in close enough proximity
for the linearized dynamics to be valid over a sufficiently long period. We
can then establish the following proposition:

######  Proposition A.5 (Escape from Saddle Points Following a Dominant
Eigenvalue Dynamics).

Consider the linearized dynamics given by
$\dot{\bm{\theta}}=-\nabla^{2}R_{S}(\bm{\theta}_{c})(\bm{\theta}-\bm{\theta}_{c})$,
we denote $\bm{H}:=-\nabla^{2}R_{S}(\bm{\theta}_{c})$. Let
$\lambda_{1}\in\mathbb{R}$ be the largest eigenvalue of $\bm{H}$, with
corresponding eigenvectors $\bm{q}_{11},\bm{q}_{12},\cdots,\bm{q}_{1l_{1}}$.
Denote
$c_{j}=\langle\bm{\theta}_{0}-\bm{\theta}_{c},\bm{q}_{1j}\rangle,\forall 1\leq
j\leq l_{1}$. Assume there exists $j$ such that $c_{j}\neq 0$, then, given any
$\varepsilon>0$, there exists a $t_{0}>0$ such that for all $t\geq t_{0}$, the
following holds:

$\left\|\dfrac{\bm{\theta}(t)-\bm{\theta}_{c}}{\mathrm{e}^{\lambda_{1}t}\sum\limits_{j=1}^{l_{1}}c_{j}\bm{q}_{1j}}-\bm{1}\right\|<\varepsilon.$
(20)

This means that, as time progresses, the direction of the parameter evolution
increasingly aligns with the dominant eigenvectors of $\bm{H}$.

This proposition is a consequence of the fact that the solution of the
differential equation is given by $\bm{\theta}(t)=e^{\bm{H}t}\bm{\theta}(0)$,
and as $t$ tends to infinity, the term corresponding to the dominant
eigenvalue in the matrix exponential $e^{\bm{H}t}$ becomes dominant.
Therefore, near the saddle point, we have
$\bm{\theta}(t)\approx\mathrm{e}^{\lambda_{1}t}\sum\limits_{j=1}^{l_{1}}c_{j}\bm{q}_{1j}+\bm{\theta}_{c}$.

###### Proof.

The solution of the ordinary differential equation (19) is
$\bm{\theta}(t)=\mathrm{e}^{t\bm{H}}(\bm{\theta}_{0}-\bm{\theta}_{c})+\bm{\theta}_{c}$.
Here, $\bm{H}=-\nabla^{2}R_{S}(\bm{\theta}_{c})$ is a real symmetric matrix,
which can be diagonalized. Let $\lambda_{1}>\lambda_{2}>\cdots>\lambda_{s}$ be
the eigenvalues of $\bm{H}$, and let
$\bm{q}_{i1},\bm{q}_{i2},\cdots,\bm{q}_{il_{i}}$ be the eigenvectors
corresponding to $\lambda_{i}$. We can then express $\bm{\theta}(t)$ as:

$\bm{\theta}(t)=\sum\limits_{i=1}^{s}\sum\limits_{j=1}^{l_{i}}\mathrm{e}^{\lambda_{i}t}\langle\bm{\theta}_{0}-\bm{\theta}_{c},\bm{q}_{ij}\rangle\bm{q}_{ij}+\bm{\theta}_{c}.$
(21)

Next, we analyze the norm of the relative difference between $\bm{\theta}(t)$
and a term dominating its growth:

$\displaystyle\left\|\dfrac{\bm{\theta}(t)-\bm{\theta}_{c}}{\sum\limits_{j=1}^{l_{1}}\mathrm{e}^{\lambda_{1}t}\langle\bm{\theta}_{0}-\bm{\theta}_{c},\bm{q}_{1j}\rangle\bm{q}_{1j}}-\bm{1}\right\|=\left\|\dfrac{\sum\limits_{i=2}^{s}\sum\limits_{j=1}^{l_{i}}\mathrm{e}^{\lambda_{i}t}\langle\bm{\theta}_{0}-\bm{\theta}_{c},\bm{q}_{ij}\rangle\bm{q}_{ij}}{\sum\limits_{j=1}^{l_{1}}\mathrm{e}^{\lambda_{1}t}\langle\bm{\theta}_{0}-\bm{\theta}_{c},\bm{q}_{1j}\rangle\bm{q}_{1j}}\right\|$
(22)
$\displaystyle\leq\sum\limits_{i=2}^{s}\sum\limits_{j=1}^{l_{i}}\mathrm{e}^{-(\lambda_{1}-\lambda_{i})t}\left\|\dfrac{\langle\bm{\theta}_{0}-\bm{\theta}_{c},\bm{q}_{ij}\rangle\bm{q}_{ij}}{\sum\limits_{j=1}^{l_{1}}\langle\bm{\theta}_{0}-\bm{\theta}_{c},\bm{q}_{1j}\rangle\bm{q}_{1j}}\right\|$
$\displaystyle\leq\mathrm{e}^{-(\lambda_{1}-\lambda_{2})t}\sum\limits_{i=2}^{s}\sum\limits_{j=1}^{l_{i}}\left\|\dfrac{\langle\bm{\theta}_{0}-\bm{\theta}_{c},\bm{q}_{ij}\rangle\bm{q}_{ij}}{\sum\limits_{j=1}^{l_{1}}\langle\bm{\theta}_{0}-\bm{\theta}_{c},\bm{q}_{1j}\rangle\bm{q}_{1j}}\right\|.$

We define
$C=\sum\limits_{i=2}^{s}\sum\limits_{j=1}^{l_{i}}\left\|\dfrac{\langle\bm{\theta}_{0}-\bm{\theta}_{c},\bm{q}_{ij}\rangle\bm{q}_{ij}}{\sum\limits_{j=1}^{l_{1}}\langle\bm{\theta}_{0}-\bm{\theta}_{c},\bm{q}_{1j}\rangle\bm{q}_{1j}}\right\|$.
By choosing
$t_{0}=\dfrac{\log\frac{C}{\varepsilon}}{\lambda_{1}-\lambda_{2}}$, we ensure
that for all $t>t_{0}$, the following condition is met:

$\left\|\dfrac{\bm{\theta}(t)-\bm{\theta}_{c}}{\sum\limits_{j=1}^{l_{1}}\mathrm{e}^{\lambda_{1}t}\langle\bm{\theta}_{0}-\bm{\theta}_{c},\bm{q}_{1j}\rangle\bm{q}_{1j}}-\bm{1}\right\|<\varepsilon.$
(23)

∎

Prop. A.5 describes that under the linearized dynamics, the parameters will
escape from the saddle point along a specific direction. However, when
considering the original nonlinear dynamics $\dot{\bm{\theta}}(t)=-\nabla
R_{S}(\bm{\theta})$, we encounter a trade-off: we should choose $t_{0}$
sufficiently large so that the trajectory can align well with the dominant
eigendirection while escaping the saddle point, but if $t_{0}$ is too large,
the linearization approximation will fail as $\bm{\theta}(t_{0})$ moves away
from $\bm{\theta}_{c}$. Li et al. (2020) (Theorem 5.3) proved a general
dynamical result through careful analysis and error control: assuming the
eigenvector corresponding to the maximum eigenvalue is unique and the
initialization is sufficiently close to the saddle point, there always exists
a suitable $t_{0}$ such that the linear dynamics can align with the dominant
eigendirection before the linearization breaks down. We can generalize this
result to the case where the eigenvector corresponding to the largest
eigenvalue is not unique:

######  Proposition A.6.

Consider the dynamics given by $\dot{\bm{\theta}}(t)=-\nabla
R_{S}(\bm{\theta})$, we use $\varphi(\bm{\theta}_{0},t)$ to denote the value
of $\bm{\theta}(t)$ in the case of $\bm{\theta}(0)=\bm{\theta}_{0}$. At a
critical point $\bm{\theta}_{c}$, we denote the negative Hessian as
$\bm{H}:=-\nabla^{2}R_{S}(\bm{\theta}_{c})$. Let $\lambda_{1}\in\mathbb{R}$ be
the largest eigenvalue of $\bm{H}$, with corresponding eigenvectors
$\bm{q}_{11},\bm{q}_{12},\cdots,\bm{q}_{1l_{1}}$. Denote
$c_{j}=\langle\bm{\theta}_{0}-\bm{\theta}_{c},\bm{q}_{1j}\rangle,\forall 1\leq
j\leq l_{1}$, and $\bm{v}_{1}=\sum_{j=1}^{l_{1}}c_{j}\bm{q}_{1j}$. Assume
there exists $j$ such that $c_{j}\neq 0$. Let
$\bm{z}_{\alpha}(t):=\varphi\left(\bm{\theta}_{c}+\alpha\bm{v}_{1},t+\frac{1}{\lambda_{1}}\log\frac{1}{\alpha}\right)$
for every $\alpha>0$, then $\bm{z}(t):=\lim_{\alpha\rightarrow
0}\bm{z}_{\alpha}(t)$ exists and is also a solution of the given dynamics,
i.e., $\bm{z}(t)=\varphi(\bm{z}(0),t)$. Furthermore, $\forall t\in\mathbb{R}$,
there exists a constant $C>0$ such that

$\left\|\varphi\left(\bm{\theta}_{c}+\alpha\bm{\theta}_{0},t+\frac{1}{\lambda_{1}}\log\frac{1}{\alpha}\right)-\bm{z}(t)\right\|_{2}\leq
C\alpha^{\frac{\lambda_{1}-\lambda_{2}}{2\lambda_{1}-\lambda_{2}}}$

for every sufficiently small $\alpha$, where $\lambda_{1}-\lambda_{2}>0$ is
the eigenvalue gap.

###### Proof.

By Theorem 5.3 in Section 5.1 of Li et al. (2020), we know that if the
eigenspace corresponding to $\lambda_{1}$ is one-dimensional, i.e., $l_{1}=1$,
then the escaping direction will be the top eigenvector direction, and the
convergence rate is
$O(\alpha^{\frac{\lambda_{1}-\lambda_{2}}{2\lambda_{1}-\lambda_{2}}})$.
Therefore, Proposition A.6 holds in this case.

Now, if the eigenspace corresponding to $\lambda_{1}$ is not one-dimensional,
we denote
$c_{j}=\langle\bm{\theta}_{0}-\bm{\theta}_{c},\bm{q}_{1j}\rangle,\forall 1\leq
j\leq l_{1}$, and $\bm{v}_{1}=\sum_{j=1}^{l_{1}}c_{j}\bm{q}_{1j}$ will be the
escaping direction. Following the same technique as in Li et al. (2020), we
can easily verify that the convergence rate remains
$O(\alpha^{\frac{\lambda_{1}-\lambda_{2}}{2\lambda_{1}-\lambda_{2}}})$.
Therefore, Proposition A.6 holds in this case as well.

∎

###  A.5 Eigenvalues and Eigenvectors of Hessian

Suppose the dominant directions fulfill specific conditions, such as any
combination
$c_{1}\bm{q}_{11}+c_{2}\bm{q}_{12}+\cdots+c_{l_{1}}\bm{q}_{1l_{1}}$, leading
to rank 1 model parameters $(\bm{A},\bm{B})$. In such scenarios, we may
observe a phenomenon where the rank of the matrix increases incrementally.

Firstly, we analyze the eigenvector structure of the Hessian matrix at the
critical point $\bm{\theta}_{c}=(\bm{A}_{c},\bm{B}_{c})$to understand why the
parameter will enter the rank-1 invariant manifold.

##### Computation of the Hessian Matrix at a Critical Point.

To compute the Hessian matrix, we first consider the gradient:

$\displaystyle
R_{S}(\bm{\theta})=\mathbb{E}_{S}\ell\left(f(\bm{x},\bm{\theta}),f^{*}(\bm{x})\right),$
$\displaystyle\nabla_{\bm{\theta}}R_{S}(\bm{\theta})=\mathbb{E}_{S}\nabla\ell\left(\bm{f}(\bm{x},\bm{\theta}),\bm{f}^{*}(\bm{x})\right)^{\top}\nabla_{\bm{\theta}}\bm{f}_{\bm{\theta}}(\bm{x}),$
$\displaystyle=\sum_{i=1}^{d^{2}}\mathbb{E}_{S}\partial_{i}\ell\left(\bm{f}_{\bm{\theta}},\bm{f}^{*}\right)\nabla_{\bm{\theta}}\left(\bm{f}_{\bm{\theta}}\right)_{i},$
$\displaystyle=\sum_{i=1}^{d^{2}}\mathbb{E}_{S}(\bm{f}_{\bm{\theta}}-\bm{f}^{*})_{i}\nabla_{\bm{\theta}}\left(\bm{f}_{\bm{\theta}}\right)_{i},$

where $\partial_{i}\ell\left(\bm{f}_{\bm{\theta}},\bm{f}^{*}\right)$ is the
$i$-th element of
$\nabla\ell\left(\bm{f}(\bm{x},\bm{\theta}),\bm{f}^{*}(\bm{x})\right)$, and
$\left(\bm{f}_{\bm{\theta}}\right)_{i}$ is the $i$-th element of the
vectorization of $\bm{f}_{\bm{\theta}}$.

For the Hessian matrix $\bm{H}_{S}(\bm{\theta})$, we have

$\displaystyle\bm{H}(\bm{\theta})$
$\displaystyle:=\nabla_{\bm{\theta}}\nabla_{\bm{\theta}}R_{S}(\bm{\theta})=\sum_{i=1}^{d^{2}}\mathbb{E}_{S}\nabla_{\bm{\theta}}\left(\partial_{i}\ell\left(\bm{f}_{\bm{\theta}},\bm{f}^{*}\right)\right)\nabla_{\bm{\theta}}\left(\bm{f}_{\bm{\theta}}\right)_{i}+\sum_{i=1}^{d^{2}}\mathbb{E}_{S}\partial_{i}\ell\left(\bm{f}_{\bm{\theta}},\bm{f}^{*}\right)\nabla_{\bm{\theta}}\nabla_{\bm{\theta}}\left(\left(\bm{f}_{\bm{\theta}}\right)_{i}\right)$
$\displaystyle=\sum_{i,j=1}^{d^{2}}\mathbb{E}_{S}\partial_{ij}\ell\left(\bm{f}_{\bm{\theta}},\bm{f}^{*}\right)\nabla_{\bm{\theta}}\left(\bm{f}_{\bm{\theta}}\right)_{i}\left(\nabla_{\bm{\theta}}\left(\bm{f}_{\bm{\theta}}\right)_{j}\right)^{\top}+\sum_{i=1}^{d^{2}}\mathbb{E}_{S}\partial_{i}\ell\left(\bm{f}_{\bm{\theta}},\bm{f}^{*}\right)\nabla_{\bm{\theta}}\nabla_{\bm{\theta}}\left(\left(\bm{f}_{\bm{\theta}}\right)_{i}\right),$
$\displaystyle=\sum_{i,j=1}^{d^{2}}\nabla_{\bm{\theta}}\left(\bm{f}_{\bm{\theta}}\right)_{i}\left(\nabla_{\bm{\theta}}\left(\bm{f}_{\bm{\theta}}\right)_{j}\right)^{\top}+\sum_{i=1}^{d^{2}}\mathbb{E}_{S}(\bm{f_{\theta}}-\bm{f^{*}})_{i}\nabla_{\bm{\theta}}\nabla_{\bm{\theta}}\left(\left(\bm{f}_{\bm{\theta}}\right)_{i}\right),$

where $\partial_{ij}\ell\left(\bm{f}_{\bm{\theta}},\bm{f}^{*}\right)$ is the
$(i,j)$-th element of
$\nabla\nabla\ell\left(\bm{f}(\bm{x},\bm{\theta}),\bm{f}^{*}(\bm{x})\right)$.

We define matrices $\bm{H}^{(1)}(\bm{\theta})$ and $\bm{H}^{(2)}(\bm{\theta})$
as follows:

$\displaystyle\bm{H}^{(1)}(\bm{\theta}):=\sum_{i,j=1}^{d^{2}}\nabla_{\bm{\theta}}\left(\bm{f}_{\bm{\theta}}\right)_{i}\left(\nabla_{\bm{\theta}}\left(\bm{f}_{\bm{\theta}}\right)_{j}\right)^{\top},$
$\displaystyle\bm{H}^{(2)}(\bm{\theta}):=\sum_{i=1}^{d^{2}}\mathbb{E}_{S}(\bm{f_{\theta}}-\bm{f^{*}})_{i}\nabla_{\bm{\theta}}\nabla_{\bm{\theta}}\left(\left(\bm{f}_{\bm{\theta}}\right)_{i}\right),$

We further denote that
$\bm{H}(\bm{\theta}):=\bm{H}^{(1)}(\bm{\theta})+\bm{H}^{(2)}(\bm{\theta})$.

For matrix factorization model, the eigenvectors of $\bm{H}^{(2)}$ has a
special structure, as characterized by Lem. A.2.

######  Lemma A.2 (Data-Independent Interleaved Structure of Eigenvectors of
$\bm{H}^{(2)}$).

Let $\bm{\theta}_{c}=(\bm{A}_{c},\bm{B}_{c})$ be any critical point of the
matrix factorization model. If $\lambda$ is an eigenvalue of
$\bm{H}^{(2)}(\bm{\theta}_{c})\in\mathbb{R}^{2d^{2}\times 2d^{2}}$, then there
exist at least $d$ eigenvectors associated with $\lambda$. These $d$
eigenvectors take the form
$\bm{v}\otimes\bm{e}_{1},\bm{v}\otimes\bm{e}_{2},\cdots,\bm{v}\otimes\bm{e}_{d}\in\mathbb{R}^{2d^{2}}$,
where $\bm{v}\in\mathbb{R}^{2d}$ is a vector to be determined and $\bm{e}_{i}$
is the unit vector representing the $i$-th column of the identity matrix
$\bm{I}_{d}\in\mathbb{R}^{d\times d}$.

###### Proof.

Let’s denote the residual matrix at the critical point as
$\delta\bm{M}=(\bm{A}_{c}\bm{B}_{c}-\bm{M})_{S_{\bm{x}}}$, where
$(\bm{A}_{c},\bm{B}_{c})$ is a critical point. For the vectorized parameter
$\bm{\theta}_{c}$, by direct calculation the matrix
$\bm{H}^{(2)}:=-\nabla^{2}R_{S}(\bm{\theta}_{c})$ can be formulated as a block
matrix, with the diagonal blocks being $0$. The specific format is as follows:

$\bm{H}^{(2)}=\left[\begin{matrix}\bm{0}&-\bm{\delta\bm{M}}\otimes\bm{I_{d}}\\\
-\bm{\delta\bm{M}}^{\top}\otimes\bm{I_{d}}&\bm{0}\end{matrix}\right].$ (24)

Next, we compute the eigenvectors of $\bm{H}^{(2)}$. Let $\lambda$ be an
eigenvalue of $\bm{H}^{(2)}$. We need to verify that
$\bm{v}\otimes\bm{e}_{1},\bm{v}\otimes\bm{e}_{2},\cdots,\bm{v}\otimes\bm{e}_{d}\in\mathbb{R}^{2d^{2}}$
are the eigenvectors of $\bm{H}^{(2)}$ corresponding to $\lambda$, for a
particular $\bm{v}\in\mathbb{R}^{2d}$ yet to be determined. That is, we need
to ensure that for all $1\leq i\leq d$, the equation
$(\bm{H}^{(2)}-\lambda\bm{I}_{2d^{2}})(\bm{v}\otimes\bm{e}_{i})=\bm{0}$ has a
non-zero solution for $\bm{v}$. Notice that

$\displaystyle(\bm{H}^{(2)}-\lambda\bm{I}_{2d^{2}})(\bm{v}\otimes\bm{e}_{i})$
$\displaystyle=\left[\begin{matrix}-\lambda\bm{I}_{d}\otimes\bm{I}_{d}&-\bm{\delta\bm{M}}\otimes\bm{I_{d}}\\\
-\bm{\delta\bm{M}}^{\top}\otimes\bm{I_{d}}&-\lambda\bm{I}_{d}\otimes\bm{I}_{d}\end{matrix}\right](\bm{v}\otimes\bm{e}_{i})$
(25)
$\displaystyle=\left(\left[\begin{matrix}-\lambda\bm{I}_{d}&-\bm{\delta\bm{M}}\\\
-\bm{\delta\bm{M}}^{\top}&-\lambda\bm{I}_{d}\end{matrix}\right]\otimes\bm{I}_{d}\right)(\bm{v}\otimes\bm{e}_{i})$
$\displaystyle=\left(\left[\begin{matrix}-\lambda\bm{I}_{d}&-\bm{\delta\bm{M}}\\\
-\bm{\delta\bm{M}}^{\top}&-\lambda\bm{I}_{d}\end{matrix}\right]\bm{v}\right)\otimes\bm{e}_{i}.$

Since $\lambda$ is an eigenvalue of $\bm{H}^{(2)}$, the determinant of the
matrix $\bm{H}^{(2)}-\lambda\bm{I}_{2d^{2}}$ equals zero. Hence

$\textnormal{det}\left(\left[\begin{matrix}-\lambda\bm{I}_{d}&-\bm{\delta\bm{M}}\\\
-\bm{\delta\bm{M}}^{\top}&-\lambda\bm{I}_{d}\end{matrix}\right]\right)^{2d}=\textnormal{det}\left(\left[\begin{matrix}-\lambda\bm{I}_{d}&-\bm{\delta\bm{M}}\\\
-\bm{\delta\bm{M}}^{\top}&-\lambda\bm{I}_{d}\end{matrix}\right]\otimes\bm{I}_{d}\right)=0.$
(26)

Consequently, from Eq. (25), we conclude that there always exists a non-zero
vector $\bm{v}\in\mathbb{R}^{2d}$ such that
$(\bm{H}^{(2)}-\lambda\bm{I}_{2d^{2}})(\bm{v}\otimes\bm{e}_{i})=0$. Since
$\bm{v}\neq\bm{0}$, it is evident that
$\bm{v}\otimes\bm{e}_{1},\bm{v}\otimes\bm{e}_{2},\cdots,\bm{v}\otimes\bm{e}_{d}\in\mathbb{R}^{2d^{2}}$
are linearly independent, and thus they represent $d$ eigenvectors
corresponding to $\lambda$. ∎

######  Proposition A.7 (Eigenvectors Structure of $\bm{H}$ at the Origin).

Consider the dynamics given by Eq. (19), where we denote
$\bm{H}:=-\nabla^{2}R_{S}(\bm{0})$. If $\lambda$ is an eigenvalue of
$\bm{H}\in\mathbb{R}^{2d^{2}\times 2d^{2}}$, then there exist at least $d$
eigenvectors associated with $\lambda$ in $\bm{H}$. These $d$ eigenvectors
take the form
$\bm{v}\otimes\bm{e}_{1},\bm{v}\otimes\bm{e}_{2},\cdots,\bm{v}\otimes\bm{e}_{d}\in\mathbb{R}^{2d^{2}}$,
where $\bm{v}\in\mathbb{R}^{2d}$ is a vector to be determined and $\bm{e}_{i}$
is the unit vector representing the $i$-th column of the identity matrix
$\bm{I}_{d}\in\mathbb{R}^{d\times d}$.

###### Proof.

In the matrix factorization model, at the origin the gradient
$\nabla_{\bm{\theta}}\left(\bm{f}_{\bm{\theta}}\right)=\bm{0}$ and thus
$\bm{H}^{(1)}(\bm{0})=0$ and the Hessian matrix reduces to
$\bm{H}^{(2)}(\bm{0})$, making
$\bm{H}=-\nabla^{2}R_{S}(\bm{0})=-\bm{H}^{(2)}(\bm{0})$.

Let’s denote the residual matrix at the origin as
$\delta\bm{M}=(\bm{A}_{c}\bm{B}_{c}-\bm{M})_{S_{\bm{x}}}$, where at the origin
$(\bm{A}_{c},\bm{B}_{c})=(\bm{0},\bm{0})$. For the vectorized parameter
$\bm{\theta}_{c}$, the matrix $\bm{H}:=-\nabla^{2}R_{S}(\bm{0})$ can be
formulated as a block matrix, with the diagonal blocks being $0$. The specific
format is as follows:

$\bm{H}=\bm{H}^{(2)}=\left[\begin{matrix}\bm{0}&-\bm{\delta\bm{M}}\otimes\bm{I_{d}}\\\
-\bm{\delta\bm{M}}^{\top}\otimes\bm{I_{d}}&\bm{0}\end{matrix}\right].$ (27)

By Lem. A.2, the proof is completed. ∎

######  Lemma A.3 (Eigenvectors Structure of $\bm{H}$ at Second-order
Stationary Point).

Let $\bm{\Omega}$ denote an $\bm{\Omega}_{k}$ invariant manifold or
sub-$\bm{\Omega}_{k}$ invariant manifold defined in Prop. A.1 and A.2, and
consider a second-order stationary point $\bm{\theta}_{c}$ within
$\bm{\Omega}$, i.e., $\nabla R_{S}(\bm{\theta}_{c})=0$ and
$\bm{\theta}^{\top}\nabla^{2}R_{S}(\bm{\theta}_{c})\bm{\theta}\geq 0$ for all
$\bm{\theta}\in\bm{\Omega}$. Then, the eigenvectors corresponding to the
negative eigenvalues of the Hessian matrix $\bm{H}(\bm{\theta}_{c})$ are
contained within the span of the eigenvectors corresponding to the negative
eigenvalues of $\bm{H}^{(2)}(\bm{\theta}_{c})$.

###### Proof.

Recall the definitions of $\bm{H}^{(1)}(\bm{\theta})$ and
$\bm{H}^{(2)}(\bm{\theta})$ given by:

$\displaystyle\bm{H}^{(1)}(\bm{\theta})$
$\displaystyle:=\sum_{i,j=1}^{d^{2}}\nabla_{\bm{\theta}}\left(\bm{f}_{\bm{\theta}}\right)_{i}\left(\nabla_{\bm{\theta}}\left(\bm{f}_{\bm{\theta}}\right)_{j}\right)^{\top},$
$\displaystyle\bm{H}^{(2)}(\bm{\theta})$
$\displaystyle:=\sum_{i=1}^{d^{2}}\mathbb{E}_{S}\left[\left(\bm{f}_{\bm{\theta}}-\bm{f}^{*}\right)_{i}\nabla_{\bm{\theta}}^{2}\left(\bm{f}_{\bm{\theta}}\right)_{i}\right].$

The Hessian matrix $\bm{H}(\bm{\theta})$ at $\bm{\theta}$ is
$\bm{H}(\bm{\theta}):=\bm{H}^{(1)}(\bm{\theta})+\bm{H}^{(2)}(\bm{\theta})$.

The manifold $\bm{\Omega}$ is an affine subspace with orthogonal complement
denoted by $\bm{\Omega}^{\perp}$. Let $\bm{H}$ have the following block
representation in the bases of $\bm{\Omega}$ and $\bm{\Omega}^{\perp}$:

$\bm{H}=\begin{bmatrix}\bm{H}_{11}&\bm{H}_{12}\\\
\bm{H}_{21}&\bm{H}_{22}\end{bmatrix}.$ (28)

Since $\bm{\Omega}$ is an invariant subspace under the gradient flow, we have
$\bm{H}\bm{\theta}\in\bm{\Omega}$ for all $\bm{\theta}\in\bm{\Omega}$, which
implies that $\bm{H}_{12}=\bm{0}$. Since $\bm{H}$ is symmetry, we have
$\bm{H}_{21}=\bm{0}$.

Let $\lambda<0$ be a negative eigenvalue of $\bm{H}:=\bm{H}(\bm{\theta}_{c})$
with $\bm{v}$ as the corresponding eigenvector. Since $\bm{H}_{11}$ is
positive semi-definite, $\bm{v}$ must lie in $\bm{\Omega}^{\perp}$.

At a critical point $\bm{\theta}_{c}=(\bm{A}_{c},\bm{B}_{c})$, by direct
calculation, the gradient $\nabla_{\bm{\theta}}\bm{f}_{\bm{\theta}_{c}}$ can
be structured as:

$\nabla_{\bm{\theta}}\bm{f}_{\bm{\theta}_{c}}=\left[\begin{array}[]{cccc}\bm{B}_{c}&&&\\\
&\bm{B}_{c}&&\\\ &&\ddots&\\\ a_{11}\bm{I}&a_{21}\bm{I}&\cdots&a_{d1}\bm{I}\\\
a_{12}\bm{I}&a_{22}\bm{I}&\cdots&a_{d2}\bm{I}\\\
\vdots&\vdots&\ddots&\vdots\\\
a_{1d}\bm{I}&a_{2d}\bm{I}&\cdots&a_{dd}\bm{I}\end{array}\right]=\left[\begin{array}[]{c}\bm{I}\otimes\bm{B}_{c}\\\
\bm{A}_{c}^{\top}\otimes\bm{I}\end{array}\right]_{2d^{2}\times d^{2}},$ (29)

where $\otimes$ denotes the Kronecker product.

Note that $\nabla_{\bm{\theta}}\left(\bm{f}_{\bm{\theta}^{*}}\right)_{j}$ is
the $j$-th column of matrix $\nabla_{\bm{\theta}}f_{\bm{\theta}_{c}}$ and it
falls precisely within the defined $\bm{\Omega}_{k}$ invariant manifold or
sub-$\bm{\Omega}_{k}$ invariant manifold $\bm{\Omega}$. Therefore, we have:

$\left(\nabla_{\bm{\theta}}\left(\bm{f}_{\bm{\theta}_{c}}\right)_{j}\right)^{\top}\bm{v}=0\quad\forall
1\leq j\leq d^{2},$ (30)

which implies that $\bm{v}$ is orthogonal to the image of
$\nabla_{\bm{\theta}}\left(\bm{f}_{\bm{\theta}_{c}}\right)$, placing it in the
null space of $\bm{H}^{(1)}(\bm{\theta}_{c})$.

As a result, we have:

$\bm{H}(\bm{\theta}_{c})\bm{v}=\left(\bm{H}^{(1)}(\bm{\theta}_{c})+\bm{H}^{(2)}(\bm{\theta}_{c})\right)\bm{v}=\bm{H}^{(2)}(\bm{\theta}_{c})\bm{v}=\lambda\bm{v}.$

Thus, the eigenvector $\bm{v}$ of the Hessian $\bm{H}(\bm{\theta}_{c})$,
corresponding to the negative eigenvalue $\lambda$, is also an eigenvector of
$\bm{H}^{(2)}(\bm{\theta}_{c})$, confirming that $\bm{v}$ is within the span
of the eigenvectors of $\bm{H}^{(2)}(\bm{\theta}_{c})$. ∎

###  A.6 Transition to the Next Rank-level Invariant Manifold

######  Proposition A.8.

The linear combination of the eigenvectors of $\bm{H}^{(2)}$:
$c_{1}(\bm{v}\otimes\bm{e}_{1})+c_{2}(\bm{v}\otimes\bm{e}_{2})+\cdots+c_{d}(\bm{v}\otimes\bm{e}_{d})$
falls within the invariant manifold $\bm{\Omega}_{1}(\bm{c})$, where
$\bm{c}=(c_{1},c_{2},\cdots,c_{d})^{\top}$.

###### Proof.

Notice that

$\displaystyle
c_{1}(\bm{v}\otimes\bm{e}_{1})+c_{2}(\bm{v}\otimes\bm{e}_{2})+\cdots+c_{d}(\bm{v}\otimes\bm{e}_{d})$
(31)
$\displaystyle=\bm{v}\otimes[c_{1}\bm{e}_{1}+c_{2}\bm{e}_{2}+\cdots+c_{d}\bm{e}_{d}]$
$\displaystyle=\bm{v}\otimes\bm{c}.$

By Definition A.1, the data-independent invariant manifold generated by
$\bm{c}$ is
$\bm{\Omega}_{1}(\bm{c})=\\{(\bm{A},\bm{B})|\bm{a}_{i},\bm{b}_{j}\in\mathrm{span}\\{\bm{c}\\},\forall
1\leq i,j\leq d\\}$. If
$\bm{\theta}=(\bm{A},\bm{B})\in\bm{\Omega}_{1}(\bm{c})$, then $\bm{A},\bm{B}$
must take the form

$\bm{A}=\left[\begin{matrix}\beta_{1}\bm{c}\\\ \beta_{2}\bm{c}\\\ \vdots\\\
\beta_{d}\bm{c}\end{matrix}\right],\quad\bm{B}=\left[\begin{matrix}\beta_{d+1}\bm{c}\\\
\beta_{d+2}\bm{c}\\\ \vdots\\\ \beta_{2d}\bm{c}\end{matrix}\right],$ (32)

for some
$\bm{\beta}=[\beta_{1},\beta_{2},\cdots,\beta_{2d}]^{\top}\in\mathbb{R}^{2d}$,
and the vectorized parameter
$\bm{\theta}=\mathrm{vec}((\bm{A},\bm{B}))\in\mathbb{R}^{2d}$ takes the form
$\bm{\beta}\otimes\bm{c}$. Let $\bm{\beta}=\bm{v}$, and the proof is complete.
∎

######  Lemma A.4.

Suppose
$\bm{\alpha}_{1},\bm{\alpha}_{2},\cdots,\bm{\alpha}_{k+1}\in\mathbb{R}^{d}$
are linearly independent, the data-independent invariant manifold exhibits the
property
$\bm{\Omega}_{k}(\bm{\alpha}_{1},\bm{\alpha}_{2},\cdots,\bm{\alpha}_{k})+\bm{\Omega}_{1}(\bm{\alpha}_{k+1})=\bm{\Omega}_{k+1}(\bm{\alpha}_{1},\bm{\alpha}_{2},\cdots,\bm{\alpha}_{k+1})$.

###### Proof.

Assume that
$\bm{\theta}=(\bm{A},\bm{B})\in\bm{\Omega}_{k+1}(\bm{\alpha}_{1},\bm{\alpha}_{2},\cdots,\bm{\alpha}_{k+1})$.
Then $\bm{A},\bm{B}$ should adopt the form:

$A=\sum\limits_{i=1}^{k}\bm{\beta}_{i}\bm{\alpha}_{i}^{\top}+\bm{\beta}_{k+1}\bm{\alpha}_{k+1}^{\top},B=\sum\limits_{i=1}^{k}\bm{\gamma}_{i}\bm{\alpha}_{i}^{\top}+\bm{\gamma}_{k+1}\bm{\alpha}_{k+1}^{\top}.$
(33)

Denote $\bm{c}_{i}=\left[\begin{matrix}\bm{\beta}_{i}\\\
\bm{\gamma}_{i}\end{matrix}\right]\in\mathbb{R}^{2d}$, then the vectorized
parameters $\bm{\theta}:=\textnormal{vec}(\bm{\theta})$ can be expressed as:

$\bm{\theta}=\sum\limits_{i=1}^{k}\bm{c}_{i}\otimes\bm{\alpha}_{i}+\bm{c}_{k+1}\otimes\bm{\alpha}_{k+1}.$
(34)

Since we know that
$\sum\limits_{i=1}^{k}\bm{c}_{i}\otimes\bm{\alpha}_{i}\in\bm{\Omega}_{k}(\bm{\alpha}_{1},\bm{\alpha}_{2},\cdots,\bm{\alpha}_{k})$
and $\bm{c}_{k+1}\otimes\bm{\alpha}_{k+1}\in\bm{\Omega}_{1}(\bm{\alpha})$, it
is straightforward to validate that
$\bm{\Omega}_{k+1}(\bm{\alpha}_{1},\bm{\alpha}_{2},\cdots,\bm{\alpha}_{k+1})=\bm{\Omega}_{k}(\bm{\alpha}_{1},\bm{\alpha}_{2},\cdots,\bm{\alpha}_{k})+\bm{\Omega}_{1}(\bm{\alpha}_{k+1})$.
∎

######  Assumption A.1 (Unique Top Eigenvalue).

Let $\delta\bm{M}=(\bm{A}_{c}\bm{B}_{c}-\bm{M})_{S_{\bm{x}}}$ be the residual
matrix at the critical point $\bm{\theta}_{c}=(\bm{A}_{c},\bm{B}_{c})$. Assume
that the top eigenvalue of the matrix $\begin{bmatrix}\bm{0}&-\delta\bm{M}\\\
-\delta\bm{M}^{\top}&\bm{0}\end{bmatrix}$ is unique.

######  Assumption A.2 (Second-order Stationary Point).

Let $\bm{\Omega}$ be an $\bm{\Omega}_{k}$ invariant manifold or
sub-$\bm{\Omega}_{k}$ invariant manifold defined in Prop. A.1 or A.2. Assume
$\bm{\theta}_{c}$ is a second-order stationary point within $\bm{\Omega}$,
i.e., $\nabla R_{S}(\bm{\theta}_{c})=0$ and
$\bm{\theta}^{\top}\nabla^{2}R_{S}(\bm{\theta}_{c})\bm{\theta}\geq 0$ for all
$\bm{\theta}\in\bm{\Omega}$.

######  Theorem A.4 (Transition to the Next Rank-level Invariant Manifold).

Consider the dynamics $\dot{\bm{\theta}}=-\nabla R_{S}(\bm{\theta})$. Let
$\varphi(\bm{\theta}_{0},t)$ denote the value of $\bm{\theta}(t)$ when
$\bm{\theta}(0)=\bm{\theta}_{0}$. Let $\bm{\Omega}$ be a $\bm{\Omega}_{k}$
invariant manifold or sub-$\bm{\Omega}_{k}$ invariant manifold. Let
$\bm{\theta}_{c}\in\bm{\Omega}$ be a critical point satisfying Assump. A.1 and
A.2. Then, for randomly selected $\bm{\theta}_{0}$, with probability 1 with
respect to $\bm{\theta}_{0}$, the limit

$\tilde{\varphi}(\bm{\theta}_{c},t):=\lim_{\alpha\to
0}\varphi\left(\bm{\theta}_{c}+\alpha\bm{\theta}_{0},t+\frac{1}{\lambda_{1}}\log\frac{1}{\alpha}\right)$
(35)

exists and falls into an invariant manifold $\bm{\Omega}_{k+1}$. Here
$\lambda_{1}$ is the top eigenvalue of negative Hessian
$-\nabla^{2}R_{S}(\bm{\theta}_{c})$.

###### Proof.

At the critical point $\bm{\theta}_{c}$, we denote the negative Hessian as
$\bm{H}:=-\nabla^{2}R_{S}(\bm{\theta}_{c})$. Let $\lambda_{1}\in\mathbb{R}$ be
the largest eigenvalue of $\bm{H}$, with corresponding eigenvectors
$\bm{q}_{11},\bm{q}_{12},\cdots,\bm{q}_{1l_{1}}$.

Denote
$c_{j}=\langle\bm{\theta}_{0}-\bm{\theta}_{c},\bm{q}_{1j}\rangle,\forall 1\leq
j\leq l_{1}$, and $\bm{v}_{1}=\sum_{j=1}^{l_{1}}c_{j}\bm{q}_{1j}$. For a
randomly selected $\bm{\theta}_{0}$, with probability 1, there exists at least
one $j$ such that $c_{j}\neq 0$.

Consider the path
$\bm{z}_{\alpha}(t):=\varphi\left(\bm{\theta}_{c}+\alpha\bm{v}_{1},t+\frac{1}{\lambda_{1}}\log\frac{1}{\alpha}\right)$
for every $\alpha>0$. By Prop. A.6, the limit $\bm{z}(t):=\lim_{\alpha\to
0}\bm{z}_{\alpha}(t)$ exists and satisfies the dynamics
$\bm{z}(t)=\varphi(\bm{z}(0),t)$.

Furthermore, $\forall t\in\mathbb{R}$, there exists a constant $C>0$ such that

$\left\|\varphi\left(\bm{\theta}_{c}+\alpha\bm{\theta}_{0},t+\frac{1}{\lambda_{1}}\log\frac{1}{\alpha}\right)-\bm{z}(t)\right\|_{2}\leq
C\alpha^{\frac{\lambda_{1}-\lambda_{2}}{2\lambda_{1}-\lambda_{2}}}$

for every sufficiently small $\alpha$, where $\lambda_{1}-\lambda_{2}>0$ is
the eigenvalue gap.

This implies that the limit $\lim_{\alpha\to
0}\varphi\left(\bm{\theta}_{c}+\alpha\bm{\theta}_{0},t+\frac{1}{\lambda_{1}}\log\frac{1}{\alpha}\right)$
exists and

$\tilde{\varphi}(\bm{\theta}_{c},t):=\lim_{\alpha\to
0}\varphi\left(\bm{\theta}_{c}+\alpha\bm{\theta}_{0},t+\frac{1}{\lambda_{1}}\log\frac{1}{\alpha}\right)=\lim_{\alpha\to
0}\varphi\left(\bm{\theta}_{c}+\alpha\bm{v}_{1},t+\frac{1}{\lambda_{1}}\log\frac{1}{\alpha}\right).$

Assuming $\bm{\theta}_{c}\in\bm{\Omega}_{k}$ and satisfies Assumps. A.1 and
A.2, we aim to show the existence of a rank-$(k+1)$ invariant manifold
$\bm{\Omega}_{k+1}$ containing $\bm{\theta}_{c}+\alpha\bm{v}_{1}$.

Define the following matrices:

$\displaystyle\bm{H}^{(1)}(\bm{\theta}_{c})$
$\displaystyle:=\sum_{i,j=1}^{d^{2}}\nabla_{\bm{\theta}_{c}}\left(\bm{f}_{\bm{\theta}_{c}}\right)_{i}\left(\nabla_{\bm{\theta}_{c}}\left(\bm{f}_{\bm{\theta}_{c}}\right)_{j}\right)^{\top},$
$\displaystyle\bm{H}^{(2)}(\bm{\theta}_{c})$
$\displaystyle:=\sum_{i=1}^{d^{2}}\mathbb{E}_{S}\left[\left(\bm{f}_{\bm{\theta}_{c}}-\bm{f}^{*}\right)_{i}\nabla_{\bm{\theta}_{c}}^{2}\left(\bm{f}_{\bm{\theta}_{c}}\right)_{i}\right].$

The Hessian matrix $\bm{H}(\bm{\theta}_{c})$ at $\bm{\theta}_{c}$ can be
expressed as
$\bm{H}(\bm{\theta}_{c}):=\bm{H}^{(1)}(\bm{\theta}_{c})+\bm{H}^{(2)}(\bm{\theta}_{c})$,
and we have $\bm{H}=-\bm{H}(\bm{\theta}_{c})$.

Assump. A.1 and Lem. A.2 imply that there exist exactly $d$ eigenvectors
associated with the top eigenvalue $\lambda_{1}$ of
$-\bm{H}^{(2)}(\bm{\theta}_{c})$. These eigenvectors are of the form
$\bm{v}\otimes\bm{e}_{i}$ for $i=1,\ldots,d$, where $\bm{v}\in\mathbb{R}^{2d}$
is a vector to be determined and $\bm{e}_{i}$ is the $i$-th standard basis
vector in $\mathbb{R}^{d}$. By Assump. A.2 and Lem. A.3, the eigenvectors
corresponding to $\lambda_{1}$ of $\bm{H}$ are contained within the span of
the eigenvectors associated with the negative eigenvalues of
$\bm{H}^{(2)}(\bm{\theta}_{c})$.

Prop. A.8 ensures that the escaping direction $\bm{v}_{1}$ lies within a
rank-1 invariant manifold $\bm{\Omega}_{1}$. Lem. A.4 then guarantees the
existence of an invariant manifold $\bm{\Omega}_{k+1}$ that includes
$\bm{\theta}_{c}+\alpha\bm{v}_{1}$. Since $\bm{\Omega}_{k+1}$ is invariant
under the gradient flow, the trajectory
$\varphi\left(\bm{\theta}_{c}+\alpha\bm{v}_{1},t+\frac{1}{\lambda_{1}}\log\frac{1}{\alpha}\right)$
remains within $\bm{\Omega}_{k+1}$.

Finally, since $\bm{\Omega}_{k+1}$ is a closed subspace, the limit
$\tilde{\varphi}(\bm{\theta}_{c},t)$ lies in $\bar{\bm{\Omega}}_{k+1}$,
concluding the proof. ∎

###  A.7 Minimum Rank

######  Theorem A.5 (Minimum Rank).

Let $\bm{\Omega}$ denote an invariant as defined previously. Assume
$\bm{W}_{t}$ achieves a global minimum within each invariant manifold
$\bm{\Omega}_{k}$. If the limit $\widehat{\bm{W}}=\lim_{\alpha\rightarrow
0}\bm{W}_{\infty}(\alpha I)$ exists and is a global optimum with
$\widehat{\bm{W}}(i,j)=\bm{M}(i,j)$ for all $(i,j)\in S_{\bm{x}}$, then

$\widehat{\bm{W}}\in\operatorname{argmin}_{\bm{W}}\operatorname{rank}(\bm{W})\quad\text{s.t.}\quad\bm{W}(i,j)=\bm{M}(i,j),\forall(i,j)\in
S_{\bm{x}}.$ (36)

###### Proof.

Consider the invariant manifold $\bm{\Omega}_{k}$, which is defined as
follows:

$\bm{\Omega}_{k}:=\bm{\Omega}_{k}(\bm{\alpha}_{1},\cdots,\bm{\alpha}_{k})=\\{(\bm{A},\bm{B})|\bm{a}_{i},\bm{b}_{\cdot,j}\in\mathrm{span}\\{\bm{\alpha}_{1},\cdots,\bm{\alpha}_{k}\\},\forall
1\leq i,j\leq d\\},$

where $\bm{a}_{i}$ denotes the $i$-th row of $\bm{A}$, $\bm{b}_{\cdot,j}$
denotes the $j$-th column of $\bm{B}$, and
$\bm{\alpha}_{1},\ldots,\bm{\alpha}_{k}$ are independent vectors that span the
invariant subspace associated with $\bm{\Omega}_{k}$.

According to Thm. A.4, the training trajectory adheres to a hierarchical
traversal across invariant manifolds. For any matrix $\bm{C}$ with
$\textnormal{rank}(\bm{C})\leq k$, we will show that there always exists
$\bm{\theta}=(\bm{A},\bm{B})\in\bm{\Omega}_{k}$ such that
$\bm{A}\bm{B}=\bm{C}$. Therefore, $\bm{\Omega}_{k}$ contains all matrices of
rank $k$.

In fact, Since $\textnormal{rank}(\bm{C})\leq k$, we can express $\bm{C}$ as a
sum of $k$ rank-one matrices:
$\bm{C}=\sum_{i=1}^{k}\bm{u}_{i}\bm{v}_{i}^{\top}$ where $\bm{u}_{i}$ and
$\bm{v}_{i}$ are column vectors. By the definition of $\bm{\Omega}_{k}$, for
any $\bm{\theta}=(\bm{A},\bm{B})\in\bm{\Omega}_{k}$, each row of $\bm{A}$ and
each column of $\bm{B}$ can be expressed as a linear combination of
$\\{\bm{\alpha}_{1},\cdots,\bm{\alpha}_{k}\\}$:
$\bm{a}_{i}=\sum_{j=1}^{k}c_{ij}\bm{\alpha}_{j}$
$\bm{b}_{\cdot,j}=\sum_{i=1}^{k}d_{ij}\bm{\alpha}_{i}$ where $c_{ij}$ and
$d_{ij}$ are scalars. We can write:

$\bm{A}=\begin{bmatrix}\sum_{j=1}^{k}c_{1j}\bm{\alpha}_{j}\\\ \vdots\\\
\sum_{j=1}^{k}c_{dj}\bm{\alpha}_{j}\end{bmatrix},\bm{B}=\begin{bmatrix}\sum_{i=1}^{k}d_{i1}\bm{\alpha}_{i}&\cdots&\sum_{i=1}^{k}d_{id}\bm{\alpha}_{i}\end{bmatrix}.$

Now, we can express the product $\bm{A}\bm{B}$ as:
$\bm{A}\bm{B}=\sum_{i=1}^{k}\sum_{j=1}^{k}\left(\sum_{l=1}^{d}c_{li}d_{jl}\right)\bm{\alpha}_{i}\bm{\alpha}_{j}^{\top}$
By choosing appropriate values for $c_{ij}$ and $d_{ij}$, we can make
$\bm{A}\bm{B}=\bm{C}$. This is possible because the outer products
$\bm{\alpha}_{i}\bm{\alpha}_{j}^{\top}$ span the same subspace as the rank-one
matrices $\bm{u}_{i}\bm{v}_{i}^{\top}$ in the expression of $\bm{C}$.

Therefore, for any matrix $\bm{C}$ with $\textnormal{rank}(\bm{C})\leq k$,
there always exists $\bm{\theta}=(\bm{A},\bm{B})\in\bm{\Omega}_{k}$ such that
$\bm{A}\bm{B}=\bm{C}$.

If the output matrix $\bm{W}_{t}$ attains optimums within each
$\bm{\Omega}_{k}$, it suggests that the optimization process is selecting the
best approximation from the set of all possible rank-$k$ matrices. Provided
that each step in the optimization is optimal, the resulting solution will
naturally be the matrix with the lowest feasible rank that satisfies the
matrix completion criteria, thereby completing the proof. ∎

###  A.8 Minimum Nuclear Norm Guarantee

######  Lemma A.5 (Minimal Nuclear Norm Computation).

Given a matrix $\bm{M}$ to be completed with observed diagonal entries, i.e.,
$\mathrm{diag}(\bm{M})=\bm{v}$, the minimal nuclear norm solution among all
possible completions is $\|\bm{v}\|_{1}$.

###### Proof.

The nuclear norm of a matrix is the dual of the spectral norm $\|\cdot\|_{2}$,
defined as:

$\|\bm{A}\|_{*}=\max_{\|\bm{X}\|_{2}\leq 1}\langle\bm{A},\bm{X}\rangle.$

Given that $\|\operatorname{diag}(\operatorname{sign}(\bm{v}))\|_{2}\leq 1$,
for any matrix $\bm{A}$ with $\operatorname{diag}(\bm{A})=\bm{v}$, it follows
that:

$\|\bm{A}\|_{*}\geq\langle\bm{A},\operatorname{diag}(\operatorname{sign}(\bm{v}))\rangle=\langle\bm{v},\operatorname{sign}(\bm{v})\rangle=\|\bm{v}\|_{1}.$

Specifically, the nuclear norm of the diagonal matrix with $\bm{v}$ on its
diagonal is $\|\operatorname{diag}(\bm{v})\|_{*}=\|\bm{v}\|_{1}$, which
establishes that the diagonal matrix with $\bm{v}$ is indeed a minimizer for
the nuclear norm. ∎

#### Diagonal Observations

######  Proposition A.9 (Minimum Nuclear Norm Guarantee in Diagonal Case).

Consider the dynamics $\dot{\bm{\theta}}=-\nabla R_{S}(\bm{\theta})$, where
$\bm{\theta}(t)=(\bm{A}(t),\bm{B}(t))$ and denote
$\bm{W}_{t}=\bm{A}(t)\bm{B}(t)$. If the observation data is diagonal, and if
for a full rank initialization $\bm{W}_{0}$, the limit
$\widehat{\bm{W}}=\lim_{\alpha\rightarrow 0}\bm{W}_{\infty}(\alpha\bm{W}_{0})$
exists and is a global optimum with $\widehat{\bm{W}}_{ij}=\bm{M}_{ij}$ for
all $(i,j)\in S_{\bm{x}}$, then

$\widehat{\bm{W}}\in\operatorname{argmin}_{\bm{W}}\|\bm{W}\|_{*}\quad\text{s.t.}\quad\bm{W}_{ij}=\bm{M}_{ij},\forall(i,j)\in
S_{\bm{x}}.$ (37)

###### Proof.

Without loss of generality, assume that $\bm{M}$ is a diagonal matrix given
by:

$\bm{M}=\begin{bmatrix}\mu_{1}&&&\\\ &\mu_{2}&&\\\ &&\ddots&\\\ &&&\mu_{d}\\\
\end{bmatrix}.$

By Lem. A.5, the minimal nuclear norm among all possible completions is
$|\mu_{1}|+|\mu_{2}|+\cdots+|\mu_{d}|$. When the matrix to be completed is
diagonal, the evolution of the $i$-th row of $\bm{A}$ is influenced only by
the $i$-th column of $\bm{B}$. Hence, the dynamics decouple into $d$
independent parts, each equivalent to learning a scalar $\mu_{i}$. The
learning process thus unfolds in $d$ stages, with each stage passing through a
critical point to learn a respective $\mu_{i}$.

By Lem. A.2, the second term of the Hessian matrix can be expressed as:

$\bm{H}^{(2)}=\left[\begin{matrix}\bm{0}&-\bm{\delta\bm{M}}\otimes\bm{I_{d}}\\\
-\bm{\delta\bm{M}}^{\top}\otimes\bm{I_{d}}&\bm{0}\end{matrix}\right].$

According to Lem. A.2, the $d$ eigenvectors of $\bm{H}^{(2)}$ take the form
$\bm{v}\otimes\bm{e}_{1},\bm{v}\otimes\bm{e}_{2},\ldots,\bm{v}\otimes\bm{e}_{d}\in\mathbb{R}^{2d^{2}}$,
where $\bm{v}\in\mathbb{R}^{2d}$ is a vector to be determined and $\bm{e}_{i}$
is the unit vector corresponding to the $i$-th column of the identity matrix
$\bm{I}_{d}\in\mathbb{R}^{d\times d}$.

(i) Suppose $\mu_{1}>\mu_{2}>\cdots>\mu_{d}>0$.

With a infinitesimal initialization, the training dynamics first focus on the
element with the largest singular value, then proceed sequentially to blocks
with smaller singular values. This pattern of learning is consistent with the
concept of “sequential learning” as reported in the literature (Gidel et al.,
2019; Gissin et al., 2019; Jiang et al., 2023).

For a diagonal observation matrix, the residual matrix $\delta\bm{M}$ at any
critical point remains a diagonal matrix. While starting to learn $\mu_{i}$
from a critical point, direct calculation confirms that vector
$\bm{v}=[\bm{e}_{i},\bm{e}_{i}]^{\top}\in\mathbb{R}^{2d}$. Escaping from each
saddle point $\bm{\theta}_{c}$, the trajectory
$\bm{\theta}(t)-\bm{\theta}_{c}$ approximates
$\sum_{i=1}^{d}c_{i}(\bm{v}\otimes\bm{e}_{i})$, which satisfies
$\bm{B}_{,i}=(\bm{A}^{\top})_{,i}$. Thus, learning a diagonal matrix $\bm{M}$
using the asymmetric model $\bm{A}\bm{B}$ is equivalent to using a symmetric
model $\bm{A}\bm{A}^{\top}$. The final outcome ensures that
$\mathrm{diag}(\bm{A}\bm{B})=\mathrm{diag}(\bm{A}\bm{A}^{\top})=\mathrm{diag}(\bm{M})$.

The nuclear norm of $\bm{A}\bm{A}^{\top}$ equals the sum of its eigenvalues,
which is precisely the trace of the matrix, and
$\operatorname{tr}(\bm{A}\bm{A}^{\top})=\operatorname{tr}(\bm{M})=\mu_{1}+\mu_{2}+\cdots+\mu_{d}$.
Therefore, the nuclear norm of the learned matrix
$\bm{W}=\bm{A}\bm{B}=\bm{A}\bm{A}^{\top}$ remains
$|\mu_{1}|+|\mu_{2}|+\cdots+|\mu_{d}|$.

(ii) If some $\mu_{i}<0$, assume without loss of generality that
$|\mu_{1}|>|\mu_{2}|>\cdots>|\mu_{n}|>0$.

While starting to learn $\mu_{i}$ from a critical point, direct calculation
confirms that
$\bm{v}=[\bm{e}_{i},\operatorname{sign}(\mu_{i})\bm{e}_{i}]^{\top}\in\mathbb{R}^{2d}$.
Escaping from each saddle point $\bm{\theta}_{c}$, the trajectory
$\bm{\theta}(t)-\bm{\theta}_{c}$ approximates
$\sum_{i=1}^{d}c_{i}(\bm{v}\otimes\bm{e}_{i})$, satisfying
$\bm{B}_{,i}=\operatorname{sign}(\mu_{i})(\bm{A}^{\top})_{,i}$. Hence,
$\bm{A}\bm{B}=\bm{A}\bm{A}^{\top}\bm{Q}$, where $\bm{Q}$ is an orthogonal
matrix given by:

$\bm{Q}=\begin{bmatrix}\operatorname{sign}(\mu_{1})&&&\\\
&\operatorname{sign}(\mu_{2})&&\\\ &&\ddots&\\\
&&&\operatorname{sign}(\mu_{d})\\\ \end{bmatrix}.$

The final result ensures that
$\mathrm{diag}(\bm{A}\bm{B})=\mathrm{diag}(\bm{A}\bm{A}^{\top}\bm{Q})=\mathrm{diag}(\bm{M})$,
meaning $\mathrm{diag}(\bm{A}\bm{A}^{\top})=\mathrm{diag}(\bm{Q}\bm{M})$. The
nuclear norm of $\bm{A}\bm{A}^{\top}$ equals the sum of its eigenvalues, which
is the trace of the matrix, and
$\operatorname{tr}(\bm{A}\bm{A}^{\top})=\operatorname{tr}(\bm{Q}\bm{M})=|\mu_{1}|+|\mu_{2}|+\cdots+|\mu_{d}|$.

Since an orthogonal transformation does not change the nuclear norm of a
matrix, the nuclear norm of the final learned matrix
$\bm{W}=\bm{A}\bm{B}=\bm{A}\bm{A}^{\top}\bm{Q}$ is still
$|\mu_{1}|+|\mu_{2}|+\cdots+|\mu_{d}|$. ∎

#### Disconnected with Complete Bipartite Components

######  Theorem A.6 (Minimum Nuclear Norm Guarantee).

Consider the dynamics $\dot{\bm{\theta}}=-\nabla R_{S}(\bm{\theta})$, where
$\bm{\theta}(t)=(\bm{A}(t),\bm{B}(t))$, and let
$\bm{W}_{t}=\bm{A}(t)\bm{B}(t)$. If the observation graph associated with the
matrix $\bm{M}$ to be completed is disconnected with complete bipartite
components, and if for a full rank initialization $\bm{W}_{0}$, the limit
$\widehat{\bm{W}}=\lim_{\alpha\rightarrow 0}\bm{W}_{\infty}(\alpha\bm{W}_{0})$
exists and is a global optimum with $\widehat{\bm{W}}_{ij}=\bm{M}_{ij}$ for
all $(i,j)\in S_{\bm{x}}$, then

$\widehat{\bm{W}}\in\operatorname{argmin}_{\bm{W}}\|\bm{W}\|_{*}\quad\text{s.t.}\quad\bm{W}_{ij}=\bm{M}_{ij},\forall(i,j)\in
S_{\bm{x}}.$ (38)

###### Proof.

Consider a matrix $\bm{M}\in\mathbb{R}^{d\times d}$ composed of $m$ connected
components, with each component forming a complete bipartite subgraph. Since
$\bm{M}$ is disconnected, it can be represented in a block diagonal form
without loss of generality:

$\bm{M}=\begin{bmatrix}\bm{M}_{1}&&&\\\ &\bm{M}_{2}&&\\\ &&\ddots&\\\
&&&\bm{M}_{m}\\\ \end{bmatrix},$

where each block $\bm{M}_{i}\in\mathbb{R}^{d_{i}\times d_{i}^{\prime}}$, and
$\sum_{i=1}^{m}d_{i}=d,\sum_{i=1}^{m}d_{i}^{\prime}=d$, representing the sum
of the dimensions of the blocks.

Each block $\bm{M}_{i}$ corresponds some singular values of the corresponding
Hessian matrix at a critical point. With a infinitesimal initialization, the
training dynamics first focus on the block with the largest singular value,
then proceed sequentially to blocks with smaller singular values. This pattern
of learning is consistent with the concept of “sequential learning” as
reported in the literature (Gidel et al., 2019; Gissin et al., 2019; Jiang et
al., 2023).

Since each connected component of $\bm{M}$ forms a complete bipartite
subgraph, the block $\bm{M}_{i}$ is fully observed. We can do singular value
decomposition (SVD) on each sub-block $\bm{M}_{i}$ as
$\bm{M}_{i}=\bm{U}_{i}\Sigma_{i}\bm{V}_{i}^{\top}$, where $\bm{U}_{i}$ and
$\bm{V}_{i}$ are orthogonal matrices, and $\Sigma_{i}$ is a diagonal matrix
with the singular values of $\bm{M}_{i}$.

Construct block diagonal matrices $\bm{U}$ and $\bm{V}$ as follows:

$\bm{U}=\begin{bmatrix}\bm{U}_{1}&&&\\\ &\bm{U}_{2}&&\\\ &&\ddots&\\\
&&&\bm{U}_{m}\\\ \end{bmatrix},\quad\bm{V}=\begin{bmatrix}\bm{V}_{1}&&&\\\
&\bm{V}_{2}&&\\\ &&\ddots&\\\ &&&\bm{V}_{m}\\\ \end{bmatrix}.$

This leads to the diagonal matrix:

$\bm{U}\bm{M}\bm{V}^{\top}=\begin{bmatrix}\bm{\Sigma}_{1}&&&\\\
&\bm{\Sigma}_{2}&&\\\ &&\ddots&\\\ &&&\bm{\Sigma}_{m}\\\
\end{bmatrix}=\begin{bmatrix}\mu_{1}&&&\\\ &\mu_{2}&&\\\ &&\ddots&\\\
&&&\mu_{d}\\\ \end{bmatrix}.$

Orthogonal transformations preserve the nuclear norm, so by Lem. A.5, the
minimal nuclear norm among all possible completions is the sum of the absolute
values of the diagonal entries, i.e., $|\mu_{1}|+|\mu_{2}|+\cdots+|\mu_{d}|$.

Consider an incomplete matrix $\bm{M}$ whose associated observational graph is
divided into $m$ connected components, denoted by $L_{1},L_{2},\ldots,L_{m}$.
For each component $L_{p}$, we define $S_{\bm{x}}^{L_{p}}$ as the subset of
observed indices within $L_{p}$, where $1\leq p\leq m$ and $S_{\bm{x}}$ is the
set of all observed indices.

For each $L_{p}$, we can identify row indices $R_{p}$ and column indices
$C_{p}$ corresponding to the observed entries in $L_{p}$ as follows:

$R_{p}=\\{i|\exists j:(i,j)\in S_{\bm{x}}^{L_{p}}\\},\quad C_{p}=\\{j|\exists
i:(i,j)\in S_{\bm{x}}^{L_{p}}\\}.$

Here, $R_{p}$ includes the row indices and $C_{p}$ includes the column indices
of the entries observed in $L_{p}$.

Define $\bm{A}^{L_{p}}$ and $\bm{B}^{L_{p}}$ as the submatrices of $\bm{A}$
and $\bm{B}$ corresponding to $R_{p}$ and $C_{p}$, respectively. The evolution
of $\bm{A}^{L_{p}}$ is influenced only by $\bm{B}^{L_{p}}$. Thus, the dynamics
decouple into $m$ independent parts, each equivalent to learning a fully
observed matrix $\bm{M}_{i}$.

Accordingly, we partition $\bm{A}$ and $\bm{B}$ into $m$ blocks:

$\bm{A}=\begin{bmatrix}\bm{A}_{1}\\\ \vdots\\\
\bm{A}_{m}\end{bmatrix},\quad\bm{B}=\begin{bmatrix}\bm{B}_{1}&\cdots&\bm{B}_{m}\end{bmatrix},$

where $\bm{A}_{i}=\bm{A}^{L_{i}}$ and $\bm{B}_{j}=\bm{B}^{L_{j}}$.

Denote $L_{i}=\|\bm{A}_{i}\bm{B}_{i}-\bm{M}_{i}\|_{2}^{2}$. The overall loss
function $L=\sum_{i=1}^{m}L_{i}$ can be decomposed into $m$ independent parts.
Performing orthogonal transformations
$\tilde{\bm{A}}_{i}=\bm{U}_{i}^{\top}\bm{A}_{i}$ and
$\tilde{\bm{B}}_{i}=\bm{B}_{i}\bm{V}_{i}$, we obtain a diagonal loss
$\tilde{L}_{i}=\|\tilde{\bm{A}}_{i}\tilde{\bm{B}}_{i}-\bm{\Sigma}_{i}\|_{2}^{2}$
for $1\leq i\leq C$.

Since gradient descent is the steepest descent in the $l_{2}$ norm and
orthogonal transformations preserve this norm, the dynamics of optimizing
$\tilde{L}_{i}$ are equivalent to those of optimizing $L_{i}$.

Without loss of generality, assume $\mu_{1}>\mu_{2}>\cdots>\mu_{d}>0$.
Otherwise, as with Prop. 37, a sign orthogonal transformation $\bm{Q}$ can be
applied without changing the nuclear norm.

By Prop. 37, the learning result for a diagonal matrix implies
$\tilde{\bm{B}}_{i}=\tilde{\bm{A}}_{i}^{\top}$, which means
$\bm{B}_{i}\bm{V}_{i}=\bm{U}_{i}^{\top}\bm{A}_{i}^{\top}$.

The final learning result

$\tilde{\bm{A}}=\begin{bmatrix}\tilde{\bm{A}}_{1}\\\ \vdots\\\
\tilde{\bm{A}}_{m}\end{bmatrix},\quad\tilde{\bm{B}}=\begin{bmatrix}\tilde{\bm{B}}_{1}&\cdots&\tilde{\bm{B}}_{m}\end{bmatrix},$

satisfies $\tilde{\bm{B}}=\tilde{\bm{A}}^{\top}$.

The result ensures that
$\mathrm{diag}(\tilde{\bm{A}}\tilde{\bm{B}})=\mathrm{diag}(\tilde{\bm{A}}\tilde{\bm{A}}^{\top})=\mathrm{diag}(\bm{\Sigma})$.
The nuclear norm of $\tilde{\bm{A}}\tilde{\bm{A}}^{\top}$ equals the sum of
its eigenvalues, which is the trace of the matrix, and
$\textnormal{Tr}(\tilde{\bm{A}}\tilde{\bm{A}}^{\top})=\textnormal{Tr}(\Sigma)=|\mu_{1}|+|\mu_{2}|+\cdots+|\mu_{d}|$.

Since orthogonal transformations do not alter the nuclear norm of a matrix,
the nuclear norm of $\bm{W}=\bm{A}\bm{B}$ is also
$|\mu_{1}|+|\mu_{2}|+\cdots+|\mu_{d}|$, concluding the proof. ∎

##  Appendix B Experimental Setup and Supplementary Experiments

In this section, we present the supplementary experiments mentioned in the
main text and the details of experiments.

###  B.1 Experimental Setup

For all our experiments, we employ gradient descent with a carefully chosen
small learning rate. A learning rate is deemed suitable when it yields a
smooth, monotonically decreasing training trajectory for the loss function,
free from any abrupt fluctuations or oscillations. We initialize all model
parameters using a Gaussian distribution with a mean of zero and a variance
that is detailed for each specific experiment. Because of the small size of
the experiment, the experiment can be completed on a single CPU.

The criterion for the sufficiency of training in all cases is a training loss
that falls below $10^{-10}$. To ascertain the rank of the matrix produced by
the learning process, we utilize a technique of extrapolation with an
infinitesimally small initialization. As depicted in Fig. 3(b), if a singular
value persistently diminishes in response to decreasing initialization
magnitudes, it is then inferred that such a singular value will not contribute
to the rank in the context of an infinitesimal initialization.

We have included code in the Supplementary Material that determines the
connectivity of a partially observed matrix and provides specific examples
illustrating the implicit regularization effects. This code can be used to
reproduce our results and explore the relationship between data connectivity
and the implicit biases of matrix factorization models in various matrix
completion scenarios.

###  B.2 Connectivity Experiments

In the connectivity experiments corresponding to Fig. 1, we explore the
behavior of randomly generated $4\times 4$ matrices with intrinsic ranks of 1,
2, and 3. To investigate the impact of sampling density on matrix
reconstruction, we sample matrices at three different levels: $2rd-r^{2}$,
which meets the threshold for exact reconstruction, $2rd-r^{2}-1$, which is
just below the threshold, and $2rd-r^{2}+1$, which exceeds the threshold.

For each sampling size, we randomly generate 10 sets of sampling positions. We
then assess the connectivity of the sampled positions and compute both the
rank and the nuclear norm of the solutions obtained through gradient descent.
As an illustration, in Fig. B1, panel (a) presents a scenario with connected
sampling positions, panel (b) shows disconnected sampling positions, and panel
(c) depicts disconnected sampling with each disconnected component forming a
complete bipartite graph.

(a) Connected (b) Disconnected (c) Disconnected with complete bipartite
components

Figure B1: Examples of connected sampling and disconnected sampling patterns
in Fig. 1.

In the connectivity experiments depicted in Figs. 2(c-d), we examine the
behavior of randomly generated matrices of size $4\times 4$ and $10\times 10$
with a rank of 1. The matrices are sampled at a size of $2rd-r^{2}$, which
corresponds to the threshold for exact reconstruction. We evaluate two
connected and one disconnected sampling patterns.

Fig. B2(a) displays the first connected sampling pattern, where all entries in
the first row and the first column are sampled. Fig. B2(b) illustrates the
second connected sampling pattern, which forms a “Z” shape across the matrix.
Fig. B2(c) shows the disconnected sampling pattern, where the samples are
split into two unconnected blocks, one in the top-left and the other in the
bottom-right of the matrix. A similar approach is taken for the $10\times 10$
matrices.

(a) Connected sampling 1 (b) Connected sampling 2 (c) Disconnected sampling

Figure B2: Examples of connected sampling and disconnected sampling patterns
in Fig. 2(c).

For Figs. 2 (a-b), we performed 100 random initializations for each
initialization scale, recorded the mean and standard deviation, and plotted
them on the figure. For a sequence $x_{1},x_{2},\cdots,x_{n}$, the mean is
calculated by $\bar{x}=\frac{1}{n}\sum_{i=1}^{n}x_{i}$, and the standard
deviation is calculated by:

$\sigma=\sqrt{\frac{1}{n-1}\sum_{i=1}^{n}(x_{i}-\bar{x})^{2}}.$

Figs. 2(c-d) demonstrate that when the target matrix has a rank of 1 and the
number of samples meets the minimum requirement for reconstruction with
connected sampling positions, the matrix factorization model is capable of
accurately reconstructing the original target matrix.

In scenarios where the target matrix has a higher rank, we have extended our
experiments accordingly. For a randomly chosen $4\times 4$ matrix with rank 2,
we selected a sample count less than the threshold of $2rd-r^{2}=12$,
specifically 10 samples, while ensuring that the sampling pattern is
connected, as shown in Fig. B3. The resulting solution from the matrix
completion has a rank of 2, which is the minimal rank that fits the sampled
data.

Figure B3: The sample pattern in Fig. B4.

Fig. B4 reveals distinct behaviors of the matrix completion depending on the
scale of initialization. With a larger initialization, the third and fourth
singular values of the completed matrix remain relatively significant,
suggesting that the model does not converge to the lowest rank solution. On
the other hand, with a smaller initialization, the third and fourth singular
values are uniformly small, indicating that the model successfully converges
to the lowest rank solution.

(a) Initialization scale $\sigma=10^{-1}$ (b) Initialization scale
$\sigma=10^{-7}$ (c) Reconstruction error

Figure B4: For a randomly selected $4\times 4$ matrix with rank 2, we chose 10
samples, fewer than the threshold $2rd-r^{2}=12$, ensuring connected sampling
positions, as shown in Fig. B3. The figures show the experimental results for
the four singular values of the matrix learned under Gaussian initialization
with mean 0 and standard deviations of $10^{-1}$ (a) and $10^{-7}$ (b),
respectively. (c) Reconstruction error of the solutions for a $4\times 4$
matrix reconstruction problem with $\bm{M}^{*}$ randomly sampled at rank $r=1$
and sample size set to the minimum reconstruction setting $n=2rd-r^{2}$. Red
and blue scatter points represent two connected sampling patterns, while green
points represent a disconnected pattern.

(a) Matrix completion (b) Singular values of $\bm{A}$

(c) Singular values of $\bm{B}$

(d) Singular values of $\bm{W}_{\text{aug}}$

(e) Training loss (f) Output at $\star$ (g) Saddle point

(h) GLRL solution

Figure B5:  (a) The matrix to be completed, with unknown entries marked by
$\star$. (b-d) Evolution of singular values for $\bm{A}$, $\bm{B}$, and
$\bm{W}_{\mathrm{aug}}$ during training. (e) Training loss for the
disconnected sampling pattern. (f) Learned values at symmetric positions
$(0,1)$ and $(1,0)$ under varying initialization scales (zero mean, varying