Connectivity Shapes Implicit Regularization in Matrix Factorization Models for Matrix Completion

Authors: zhiwei  bai
jiajie  zhao
yaoyu  zhang
Venue: NeurIPS_cc_2024_Conference
Published: 2024-05-13 16:03:02

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 key 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.

Main paper: 
 introduction: 
 Overparameterized models have the capacity to easily fit data with random labels (Zhang et al., 2017 (Zhang et al., , 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, f θ = AB with θ = (A, B), A, B ∈ R d×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 M ∈ R d×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 plays a key role in the implicit regularization effects. Data connectivity, in the context of this paper, refers to the way observed entries in the matrix are linked through shared rows Figure 1 : The connectivity of observed data affects the implicit regularization. The ground truth matrix M * ∈ R 4×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 2rdr 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.) or columns. A set of observations is considered connected if there's a path between any two observed entries via other observed entries in the same rows or columns. This concept plays a crucial role in determining the behavior of matrix factorization models, as we will demonstrate throughout this paper. As shown in Fig. 1 , we sample observations randomly from a ground truth matrix M * ∈ R d×d with rank(M * ) < d and train models f θ = AB, A, B ∈ R d×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 Ω 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-Ω k invariant manifolds emerge within the Ω 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). 
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 Zhang et al. ( , 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 Ω k invariant manifolds we propose are essentially a manifestation of condensation. Zhang et al. (2021 Zhang et al. ( , 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. 
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 M ∈ R d×d , the goal is to predict the entirety of M based on its observed elements. The set of observed entries is represented as S = {(i k , j k ), M i k ,j k } n k=1 , where (i k , j k ) indicates the row and column indices, and M i k ,j k is the corresponding value assumed non-zero in the matrix. The set of observed indices is defined as S x = {(i k , j k )} n k=1 . Entries that are not observed, denoted by ⋆, are considered missing or unknown. The positions of observed elements in the matrix M are defined by a binary observation matrix P , where P ij = 1 indicates that M ij is observed, and P ij = 0 indicates that M ij is unobserved.Matrix factorization model. Matrix factorization is a prevalent approach for addressing the matrix completion problem. It reconstructs the matrix W ∈ R d×d through the product W = AB, where A ∈ R d×r and B ∈ R r×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:f θ = AB, A, B ∈ R d×d .(1)The matrix factorization model parameters are denoted by θ = (A, B), identified with its vectorized form vec(θ) ∈ R 2d 2 . The augmented matrix is W ⊤ aug = A ⊤ B ⊤ ∈ R d×2d , and row(A) and col(B) denote the row and column spaces of A and B, respectively. The augmented matrix W aug plays a crucial role in our subsequent analysis, particularly in characterizing the intrinsic invariant manifolds Ω k of the optimization process. Specifically, it allows us to establish the relationship rank(A) = rank(B ⊤ ) = rank(W aug ), which is important to understanding the invariance property under gradient flow.Loss function. The learning process for the parameters θ = (A, 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 asEQUATIONwhere a i and b •,j represent the i-th row and j-th column of matrix A and B, respectively. The residual matrix δM = (AB -M ) Sx has elements δM ij = (AB) ij -M ij for (i, j) ∈ S x and δM ij = 0 for (i, j) / ∈ S x . The training dynamics follow the gradient flow of R S (θ):EQUATIONIn all experiments, θ 0 ∼ N (0, σ 2 ) is initialized from a Gaussian distribution with mean 0 and small variance σ 2 . We use gradient descent with a small learning rate to approximate the gradient flow dynamics (Please refer to Appendix B.1 for the detailed experiment setup). 
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 M to be completed and its observation matrix P , the associated observation graph G M is the bipartite graph with adjacency matrix 0 P ⊤ P 0 , with isolated vertices removed.Definition 2 (Connectivity). Given a incomplete matrix M to be completed, it is considered connected if its associated observation graph G M is connected; otherwise, it is disconnected. The connected components of M are defined as the connected components of G 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 M is considered disconnected with complete bipartite components if its associated observation graph G 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: M 1 (disconnected), M 2 (disconnected with complete bipartite components), and M 3 (connected). Fig. A1 of Appendix B illustrates the associated graphs G M .EQUATION)Figs. 2(a-b) compare the learned matrices with the ground truth (GT) solutions having the smallest nuclear norm and rank. For disconnected M 1 (blue bars), the learned solution achieves neither the smallest nuclear norm nor rank. For disconnected M 2 with complete bipartite components (green bars), the learned matrix has the smallest nuclear norm but not rank. For connected 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.To thoroughly study all possible cases, we examine all sampling patterns of the 3×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 casesOur experimental findings indicate that when the observational data is connected, matrix factorization models often learn the lowest-rank solution starting from a small initialization. However, the required scale of initialization is not constant across different instances. We empirically observed that if the magnitude of the numerical values in the matrix to be completed varies significantly, an extremely small initialization is necessary, which, in some cases, can exceed machine precision.Consider the following two simple 2 × 2 matrix completion problems, with the only difference being that the number 3 in the first row is replaced by 20. When training begins from a small initialization, for M 4 , the fourth element only needs to learn the value 6 to be a rank-1 solution. However, for M 5 , the fourth element needs to learn the value 40 to achieve rank-1.M 4 = 1 2 3 × , M 5 = 1 2 20 × .Fig. B6 illustrates the difficulty in learning these two examples. For M 4 , an initialization variance of approximately 10 -7 is sufficient to learn the lowest-rank solution. In contrast, for M 5 , an extremely small initialization variance is required, making it challenging to learn a rank-1 solution. Yet, with an exceedingly small initialization variance of 10 -83 , we can still observe the second singular value plummeting to zero. The origin is a saddle point. The smaller the initialization, the longer it will stay at the origin. If initialization continues to decrease, training will stagnate. Therefore, if the magnitude difference is even greater, such as replacing 20 with 100 in M 5 , then with the small initialization allowed by machine precision, it is nearly impossible to learn the value 200 completely.For the matrix factorization model f θ = AB, the Hessian matrix at 0 has strictly negative eigenvalues, making the origin a strict saddle point. Under small random initialization, gradient descent escapes this saddle at an exponential rate. Our Theorem 1 ensures that only strict saddle points and global minima exist as critical points in the loss landscape. Subsequent saddle points on invariant manifolds are also strict, with exponential escape speeds, maintaining acceptable optimization process.Reaching the lowest rank solution requires parameters to escape the saddle point along the unique top eigen-direction. Fig. B6 illustrates the relationship between required initialization scale and observation magnitude differences. For lowest possible rank with large numerical magnitude differences in observations, extremely small initialization is necessary. However, for approximately low rank solutions, a relatively small initialization suffices.  
connected case: 
 This section empirically demonstrates the detailed dynamics of connected observed data. Fig. 3 (a)shows the connected target matrix M with a single unknown element denoted by ⋆. The rank of M is at least three and equals three if and only if ⋆ = 1.2.Learning lowest-rank solution. We initialize A and 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 random initialization (Gaussian distribution with mean 0 and variance 10 -16 ), 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 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 A and the column space of B. Starting with small initialization, we track the rank (number of significantly non-zero singular values) of W = AB, A, B, and the augmented matrix W 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 rank(A) = rank(B ⊤ ) = rank(W aug ), which implies that the row space of A and the column space of 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.In this section, we present a typical experiment in the disconnected situation. As depicted in Fig. 4 (a), the target matrix M contains four unknown elements denoted by ⋆ and is disconnected. The rank of M is at least one, and there are infinitely many rank-1 solutions.Alignment of the row space of A and the column space of 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 rank(A) = rank(B ⊤ ) = rank(W 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 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 A and columns of B, and another for the 2nd row of A and column of 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). In Fig. 4 (e), we fixed a rank-1 matrix and explored all nine disconnected sampling patterns with 5 observations. For each pattern, we conducted experiments with small initializations. The loss curves consistently indicate that in disconnected cases, the model learns a sub-optimal solution in the rank-1 manifold, ultimately resulting in a rank-2 solution. This demonstrates that regardless of the specific disconnected sampling pattern, the model fails to achieve the optimal low-rank solution. (Li et al., 2020) .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 := Ω k (α 1 , • • • , α k ) = {θ = (A, B) | row(A) = col(B) = span{α 1 , • • • , α k }}.The manifold Ω k possesses the following properties:(i) Invariance under Gradient Flow: Given data S and the gradient flow dynamics θ = -∇R S (θ), if the initial point θ 0 ∈ Ω k , then θ(t) ∈ Ω k for all t ≥ 0.(ii) Intrinsic Property: Ω k is a data-independent invariant manifold, meaning that for any data S, Ω k remains invariant under the gradient flow dynamics.(iii) Hierarchical Structure: The manifolds Ω k form a hierarchy: show that the training process with small initialization consistently satisfies rank(A) = rank(B ⊤ ) = rank(W aug ), aligning with the Ω k invariant manifold. Since a non-zero initialization in practice, the training trajectory is close to the Ω k invariant manifold, approaches a critical point, and transitions to the next level invariant manifold without getting trapped.Ω 0 ⫋ Ω 1 ⫋ • • • ⫋ Ω k-1 ⫋ Ω k .In both connected and disconnected scenarios, we observe a step-by-step hierarchical Ω k invariant manifold traversal. In the connected case, at each level we observe that the model reaches an optimal solution (Fig. 3 ). However, in the disconnected case, we can prove that each connected component induces a sub-Ω k invariant manifold, leading to the experimentally observed sub-optimal solution (see Fig. 4 ).Proposition 2 (Intrinsic Sub-Ω k Invariant Manifold). (see Appendix A Prop. A.2 for Proof) Let f θ = AB be a matrix factorization model, M be an incomplete matrix and Ω k be an invariant manifold defined in Prop. 1. If M is disconnected with m connected components, then there exist m sub-Ω k manifolds ω k such that ω k ⫋ Ω k , each possessing the following properties:(i) Invariance under Gradient Flow: Given data S and the gradient flow dynamics θ = -∇R S (θ), if the initial point θ 0 ∈ ω k , then θ(t) ∈ ω k for all t ≥ 0.(ii) Intrinsic Property: ω k is a data-value-independent invariant manifold, meaning that for a fixed sampling pattern in M and any observed values S, ω k remains invariant under the gradient flow.(iii) Strict Subset Relation: 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 A and the 4 column vectors of matrix B. For ease of visualization, we randomly project them onto two dimensions and plot them in polar coordinates.The output set {f θ | θ ∈ ω k } is a proper subset of {f θ | θ ∈ Ω k }, namely, {f θ | θ ∈ ω k } ⫋ {f θ | θ ∈ Ω k }.Fig. 5 (a) illustrates the trajectory of the experiment in Fig. 4 . In the disconnected case, sub-Ω k invariant manifolds exist and attract the dynamics, leading the model to learn sub-optimal solutions on the entire Ω 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 (θ) are either strict saddle points or global minima.Gradient descent easily escapes saddle points (Lee et al., 2016 (Lee et al., , 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 Ω 1 invariant manifold, resulting in a rank increase of one at a time. Under reasonable assumptions, we prove that the training trajectory follows the Ω k invariant manifold step by step.Assumption 1 (Unique Top Singular Value). Let δM = (A c B c -M ) Sx be the residual matrix at the critical point θ c = (A c , B c ). Assume that the largest singular value of δM is unique.Assumption 2 (Second-order Stationary Point). Let Ω be an Ω k invariant manifold or sub-Ω k invariant manifold defined in Prop. 1 or 2. Assume θ c is a second-order stationary point within Ω, i.e., ∇R S (θ c ) = 0 and θ ⊤ ∇ 2 R S (θ c )θ ≥ 0 for all θ ∈ Ω.Theorem 2 (Transition to the Next Rank-level Invariant Manifold). (see Appendix A Thm. A.4 for proof) Consider the dynamics θ = -∇R S (θ). Let φ(θ 0 , t) denote the value of θ(t) when θ(0) = θ 0 .Let Ω be an Ω k or sub-Ω k invariant manifold. Let θ c ∈ Ω be a critical point satisfying Assump. 1 and 2. Then, for randomly selected θ 0 , with probability 1 with respect to θ 0 , the limitφ(θ c , t) := lim α→0 φ θ c + αθ 0 , t + 1 λ 1 log 1 α (5)exists and falls into an invariant manifold Ω k+1 . Here λ 1 is the top eigenvalue of -∇ 2 R S (θ c ).Proof sketch. The main idea is to analyze the local dynamics near the critical point θ c . The nonlinear dynamics can be approximated linearly in the vicinity of θ c : dθ dt ≈ H(θ 0θ c ), where H = -∇ 2 R S (θ c ) is the negative Hessian matrix. For exact linear approximation, the solution is:θ(t) = e tH (θ 0 -θ c )+θ c . Let λ 1 > λ 2 > ... > λ s be the eigenvalues of H, with corresponding eigenvectors q ij . We can express θ(t) as: θ(t) = s i=1 li j=1 e λit ⟨θ 0θ c , q ij ⟩q ij + θ c . For sufficiently large t 0 , the dynamics follows a dominant eigenvalue dynamics: θ(t 0 ) = l1 j=1 e λ1t0 ⟨θ 0θ c , q 1j ⟩q 1j + O(e λ2t0 ). Through detailed analysis of the eigenvalues and eigenvectors of the Hessian matrix (please refer to Lems A.2-A.4 of Appendix A), we show that if the largest singular value of residual matrix δM at θ c is unique and θ c is a second-order stationary point within Ω, the first principal component l1 j=1 e λ1t0 ⟨θ 0θ c , q 1j ⟩q 1j will happen to be an Ω 1 invariant manifold. Consequently, escaping θ c increases the rank by 1, entering Ω k+1 .Remark. Assump. 1 ensures that upon departing from a critical point θ c , the trajectory is constrained to escape along a single dominant eigendirection corresponding to the largest singular value. This assumption holds for randomly generated matrix with probability 1, making it a reasonable condition in most practical scenarios. In Sec A.7 of Appendix A, we provide an special example to illustrate the situation where Assump. 1 does not hold. Remark. To ensure the escape direction falls within the Ω k+1 invariant manifold, the Hessian's top eigenvectors must satisfy rank(A) = rank(B ⊤ ) = rank(W aug ). The condition that θ c is a secondorder stationary point within Ω in Assump. 2 guarantees this Hessian structure. Our Assump. 2 is more general than conditions proposed by Li et al. (2020) , as it remains valid across both connected and disconnected configurations. Empirical findings (Figs. 3 and 4 ) indicate that this assumption consistently holds in practical scenarios.Thm. 2 provides a characterization of the escape trajectory. It shows that as the point approaches a second-order stationary point θ c ∈ Ω k , the trajectory generically converges to a well-defined limit within Ω k+1 . Since the origin 0 is always a second-order stationary point of Ω 0 , the theorem implies that the trajectory escaping from a small initialization will be close to Ω 1 . This iterative process gives rise to the phenomenon of Hierarchical Invariant Manifold Traversal (HIMT), which involves a sequential progression through these Ω k manifolds. 
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 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 Ω 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 Ω k , a solution of minimal rank can be achieved naturally.Theorem 3 (Minimum Rank). (see Appendix A Thm. A.5 for proof) Consider the dynamics θ = -∇R S (θ), where θ(t) = (A(t), B(t)), and denote W t = A(t)B(t). Assume W t achieves an optimal within each invariant manifold Ω k . For a full rank initialization W 0 , if the limit W = lim α→0 W ∞ (αW 0 ) exists and is a global optimum withEQUATIONFor a disconnected matrix M , our theoretical results (Prop. 2) and experiments (Fig. 4 ) confirm the existence of sub-Ω 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 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. A.6 for proof) Consider the dynamics θ = -∇R S (θ), where θ(t) = (A(t), B(t)), and let W t = A(t)B(t). If the observation graph associated with the incomplete matrix M is disconnected with complete bipartite components, and if for a full rank initialization W 0 , the limit W = lim α→0 W ∞ (αW 0 ) exists and is a global optimum withEQUATION 
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 Ω k invariant manifold. However, further investigation is needed to uncover the underlying mechanisms by which connectivity facilitates optimal attainment across different Ω 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. B9, B11). These findings suggest that the hierarchical invariant manifold traversal process uncovered in our study may have broader implications and merit further exploration. 
safeguards: 
 Question: Does the paper describe safeguards that have been put in place for responsible release of data or models that have a high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)? Answer: [NA] Justification: The paper poses no such risks. Guidelines:• The answer NA means that the paper poses no such risks.• Released models that have a high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere to usage guidelines or restrictions to access the model or implementing safety filters. • Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images. • We recognize that providing effective safeguards is challenging, and many papers do not require this, but we encourage authors to take this into account and make a best faith effort. 12. Licenses for existing assets Question: Are the creators or original owners of assets (e.g., code, data, models), used in the paper, properly credited and are the license and terms of use explicitly mentioned and properly respected? Answer: [NA] Justification: The paper does not use existing assets. Guidelines:• The answer NA means that the paper does not use existing assets.• The authors should cite the original paper that produced the code package or dataset.• The authors should state which version of the asset is used and, if possible, include a URL. • The name of the license (e.g., CC-BY 4.0) should be included for each asset.• For scraped data from a particular source (e.g., website), the copyright and terms of service of that source should be provided. • If assets are released, the license, copyright information, and terms of use in the package should be provided. For popular datasets, paperswithcode.com/datasets has curated licenses for some datasets. Their licensing guide can help determine the license of a dataset. 
 
 Appendix: In this section, we give all proofs for our theoretical results mentioned in the main text.Proposition A.1 (Hierarchical Intrinsic Invariant Manifold (HIIM)). Let f θ = AB be a matrix factorization model and{α 1 , • • • , α k } be k linearly independent vectors. Define the manifold Ω k as Ω k := Ω k (α 1 , • • • , α k ) = {θ = (A, B) | row(A) = col(B) = span{α 1 , • • • , α k }}.The manifold Ω k possesses the following properties:(i) Invariance under Gradient Flow: Given data S and the gradient flow dynamics θ = -∇R S (θ), if the initial point θ 0 ∈ Ω k , then θ(t) ∈ Ω k for all t ≥ 0.(ii) Intrinsic Property: Ω k is a data-independent invariant manifold, meaning that for any data S, Ω k remains invariant under the gradient flow dynamics.(iii) Hierarchical Structure: The manifolds Ω k form a hierarchy:Ω 0 ⫋ Ω 1 ⫋ • • • ⫋ Ω k-1 ⫋ Ω k .Proof. (i) Invariance under Gradient Flow.By definition,EQUATIONwhereI i = {j|∃i : (i, j) ∈ S x }, I j = {i|∃j : (i, j) ∈ Sx }, a i and b •,j represent the i-th row and j-th column of A and B, respectively.For any (i, j) ∈ S x , the evolution of a i is coupled with b •,j for j∈ I i . The condition row(A) = col(B) = span{α 1 , • • • , α k } ensures the existence of k linearly independent vectors α 1 , • • • , α k ∈ R d such that a i , b •,j ∈ span{α 1 , • • • , α k } for all 1 ≤ i, j ≤ d.Consequently, if a i and b •,j are initially in span{α 1 , α 2 , • • • , α k }, they will continue to evolve within this subspace under the gradient flow dynamics. Additionally, for (i, j) / ∈ S x , the gradients for the corresponding a i and b •,j will be zero, provided their initial values are zero, maintaining this state throughout the evolution.(ii) Intrinsic Property.As demonstrated in part (i), Ω 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 Ω k encompasses matrices of rank up to k, including those of lower ranks. Consequently, the manifolds exhibit the following hierarchical nesting:Ω 0 ⫋ Ω 1 ⫋ • • • ⫋ Ω k-1 ⫋ Ω k . Proposition A.2 (Intrinsic Sub-Ω k Invariant Manifold).Let f θ = AB be a matrix factorization model, M be an incomplete matrix and Ω k be an invariant manifold defined in Prop. 1. If M is disconnected with m connected components, then there exist m sub-Ω k manifolds ω k such that ω k ⫋ Ω k , each possessing the following properties:(i) Invariance under Gradient Flow: Given data S and the gradient flow dynamics θ = -∇R S (θ), if the initial point θ 0 ∈ ω k , then θ(t) ∈ ω k for all t ≥ 0.(ii) Intrinsic Property: ω k is a data-value-independent invariant manifold, meaning that for a fixed sampling pattern in M and any observed values S, ω k remains invariant under the gradient flow.(iii) Strict Subset Relation: The output set {f θ | θ ∈ ω k } is a proper subset of {f θ | θ ∈ Ω k }, namely, {f θ | θ ∈ ω k } ⫋ {f θ | θ ∈ Ω k }. Proof. Existence.Let us consider an incomplete matrix M whose associated observational graph is divided into m connected components, denoted by L 1 , L 2 , . . . , L m . For each component L p , we define S Lp x as the subset of observed indices within L p , where 1 ≤ p ≤ m and S 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: LetR p = {i|∃j : (i, j) ∈ S Lp x }, C p = {j|∃i : (i, j) ∈ S Lp x }.Ω k := Ω k (α 1 , • • • , α k ) = {θ = (A, B) | row(A) = col(B) = span{α 1 , • • • , α k }} be the given Ω k invariant manifold.The sub-Ω k invariant manifold associated with the connected component L p can be defined asω Lp k := {(θ = (A, B)) | row(A Lp ) = col((B Lp )) = span{α 1 , • • • , α k }, A Lp r = B Lp r = 0}. (9) It is easy to check ω Lp k is a proper subset of Ω k . (i) Invariance under Gradient Flow. The condition row(A Lp ) = col((B Lp )) = span{α 1 , • • • , α k } along with A Lp r = B Lp r = 0 guarantees that a i , b •,j ∈ span{α 1 , . . . , α k } for all (i, j) ∈ S Lp x , and a i , b •,j = 0 for all (i, j) / ∈ S Lp x .In other words, the sub-Ω k invariant manifold ω Lp k is the set of all pairs (A, B) where, for each observed position (i, j) in the connected component L p , the vectors a i and b •,j lie within the span of {α 1 , • • • , α k }, and for any position not in S Lp x , the vectors are zero. Considering the dynamics expressed in equation ( 8), it is evident that the evolution of a i is influenced by b •,j for (i, j) ∈ S Lp x . Hence, if a i and b •,j are initially in the span of {α 1 , α 2 , • • • , α k }, they will continue to evolve within this span under the gradient flow dynamics. Moreover, for positions (i, j) / ∈ S Lp x , we consider the following scenarios:• For (i, j) / ∈ S x , since the matrix entry M ij does not contribute to the loss R S (θ), the gradients for corresponding a i and b •,j will perpetually be zero. Thus, if their initial values are zero, they will remain zero throughout the evolution.• For (i, j) ∈ S x but not in S Lp x , the dynamics corresponding to different connected components are decoupled. Therefore, if the initial values for a i and b •,j are zero, they will stay zero during the evolution.(ii) Intrinsic Property.As established in (i), the manifold ω Lp k is invariant under gradient flow for any data S with a fixed sampling pattern, qualifying it as a data-value-independent invariant manifold. {f θ | θ ∈ ω k } forms a strict subset of {f θ | θ ∈ Ω k }, as stated by {f θ | θ ∈ ω k } ⫋ {f θ | θ ∈ Ω k }.Definition A.1 (Associated Observation Graph). Given a incomplete matrix M to be completed and its observation matrix P , the associated observation graph G M is the bipartite graph with adjacency matrix 0 P ⊤ P 0 , with isolated vertices removed.Definition A.2 (Connectivity). A matrix M to be completed is considered connected if its associated observation graph G M is connected, otherwise, we call it disconnected. The connected components of M are defined as the connected components of this graph. Definition A.3 (Disconnectivity with Complete Bipartite Components). A matrix M to be completed is considered disconnected with complete bipartite components if its associated observation graph G 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 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 0 P ⊤ P 0 , where P is the binary observation matrix indicating the positions of the observed entries in 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 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 M .A matrix 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 ′ , then G ′ is connected if and only if the line graph of G ′ , denoted as L(G ′ ), is connected.Proof. =⇒ Assume G ′ is connected. Consider any two nodes in L(G ′ ), which correspond to two edges in G ′ , say e 1 and e 2 . Since G ′ 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 ′ , which in turn corresponds to a path connecting the nodes representing e 1 and e 2 in L(G ′ ). Therefore, L(G ′ ) is connected.⇐= Conversely, assume L(G ′ ) is connected. Consider any two vertices v 1 and v 2 in G ′ . Since G ′ 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 ′ ) is connected, there exists a path connecting the nodes representing e 1 and e 2 in L(G ′ ). This path corresponds to a sequence of edges in G ′ , which in turn corresponds to a path connecting v 1 and v 2 in G ′ . Therefore, G ′ 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 ′ , then G ′ is connected if and only if the line graph of G ′ , denoted as L(G ′ ), is connected.Proposition A.3. Given a incomplete matrix M , the connectivity of 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 ′ , then G ′ is connected if and only if the line graph of G ′ , denoted as L(G ′ ), 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 M , consider the gradient flow dynamics of matrix factorization models,EQUATION)The dynamics are said to be decoupled if there exist disjoint subsets of indices R 1 , R 2 , . . . , R k ⊆ {1, 2, . . . , d} for the rows of A and C 1 , C 2 , . . . , C k ⊆ {1, 2, . . . , d} for the columns of B, such that for each l ∈ {1, 2, . . . , k}, the dynamics of {a i : i ∈ R l } and {b •,j : j ∈ 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 A and a subset of columns of B. If such a division is not possible, the dynamics are said to be coupled. If 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: M 1 (disconnected), M 2 (disconnected with complete bipartite components), and M 3 (connected).EQUATIONThe observation matrix P is:P 1 = 1 1 0 1 0 0 0 0 1 , P 2 = 1 1 0 1 1 0 0 0 1 , P 3 = 1 1 0 1 1 0 1 0 1 (13)And the adjacency matrix is:A 1 = 0 P ⊤ 1 P 1 0 , A 2 = 0 P ⊤ 2 P 2 0 , A 3 = 0 P ⊤ 3 P 3 0 (14)Given the adjacency matrix A, we can obtain a graph G M . Fig. A1 illustrates the associated graphs G M , from which we can see that M 1 is disconnected, with its associated observation graph consisting of two connected components. M 2 is also disconnected, but each connected component of its associated observation graph forms a complete bipartite subgraph. In contrast, M 3 is connected, and its associated observation graph consists of a single connected component. M 1 is disconnected, with its associated observation graph consisting of two connected components. M 2 is also disconnected, but each connected component of its associated observation graph forms a complete bipartite subgraph. In contrast, M 3 is connected, and its associated observation graph consists of a single connected component.In this paper, we focus on the problem of asymmetric matrix factorization. Previous literature (Gunasekar et al., 2017; Li et al., 2018 Li et al., , 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 W = AB can be transformed into a symmetric situation by settingU = A B ⊤ ∈ R 2d×d .We then consider the model W ′ = U U ⊤ , which corresponds to the following matrix completion problem:AA ⊤ AB B ⊤ A ⊤ B ⊤ B .We define the loss as (1). All stationary points of L : R d×d → R, L(U ) = 1 2 f U U ⊤ are either strict saddles or global minimizers; (2). For any random initialization, GF (1) converges to strict saddles of L(U ) with probability 0.L ′ A B C D = 1 2 L(B) + 1 2 L(C ⊤ ).The proof of this theorem relies heavily on Theorem A.2 of Du and Lee (2018) , which requires the parameter matrix U ∈ R d×k to satisfy the condition that k ≥ d. In the case of symmetric matrix factorization, where U ∈ R d×d , this condition is naturally met. However, for asymmetric matrix factorization, where U = A B ⊤ ∈ R 2d×d , this condition is not satisfied, and thus the proof of Theorem A.1 is only applicable to the symmetric case of matrix factorization. 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 (θ) are either strict saddle points or global minima.Proof. We start by recalling the definition of the loss function:R S (θ) = L(A, B) = 1 2 ∥AB -M ∥ 2 Sx = 1 2 (i,j)∈Sx ((AB) ij -M ij ) 2 .Let θ = (A, B) denote a critical point. We define a new matrix, δM , as the difference between the product of A and B and the matrix M , with this difference being computed only over the indices in the set S x . More formally, δM = (AB -M ) Sx , where the elements of δM are given by:• For (i, j) ∈ S x , we haveδM ij = (AB) ij -M ij .• For (i, j) / ∈ S x , we have δM ij = 0.This definition of δM ensures that we only consider the differences in the entries that belong to the set S x , while all other entries are set to zero.Consider the function:EQUATIONwhere the inner product of two matrices A, B is defined as ⟨A, B⟩ := Tr(AB ⊤ ).At the critical point, the first order term ⟨δM , εB + Aη⟩ Sx equals 0. The Hessian operator, representing the second order term, is given by:h A,B (ε, η) = 1 2 ∥εB + Aη∥ 2 Sx + ⟨δM , εη⟩ Sx .Our goal is to demonstrate that if δM ̸ = 0, there always exists ε, η such that h A,B (ε, η) < 0. To this end, we consider the ranks of matrices A and B in two cases:(i) rank(A) < d or rank(B) < d:Without loss of generality, we assume δM ij := δM ij ̸ = 0 for some (i, j) ∈ S x , and rankA < d. Under these conditions, there exists a non-zero vector v such that Av = 0.We set η * ,j = v and η * ,s = 0 for s ̸ = j, where η * ,j denotes the j-th column of the matrix η. Let ε * i = w ⊤ ∈ R d and ε s = 0 for s ̸ = i, where ε * i denotes the i-th row of the matrix ε. We then have:h A,B (ε * , η * ) = 1 2 ∥ε * B∥ 2 Sx + δM ij w ⊤ v ≤ 1 2 ∥w ⊤ B∥ 2 F + δM ij w ⊤ v = 1 2 w ⊤ BB ⊤ w + δM ij w ⊤ v.We define g(w, v) = 1 2 w ⊤ BB ⊤ w + δM ij w ⊤ v and consider:g(-αδM ij v, v) = 1 2 α 2 δM 2 ij v ⊤ BB ⊤ v -αδM 2 ij v ⊤ v = 1 2 α 2 δM 2 ij (v ⊤ BB ⊤ v -2v ⊤ v).For 0 < α < 2 λ B,max , where λ B,max represents the top eigenvalue of BB ⊤ , we find g(-αδM ij v, v) < 0.Therefore, when w = -αδM ij v, we obtain h A,B (ε * , η * ) < 0. This immediately implies the critical point θ = (A, B) is a strict saddle point.(ii) rank(A) = rank(B) = d:Let ε = αδM B -1 and η = 0. In this scenario, the first order term ⟨δM , εB + Aη⟩ Sx in Eq. ( 15) simplifies to α∥δM ∥ 2 Sx . At a critical point, this quantity equals zero, it implies that δM = 0. This in turn implies that the critical point θ = (A, B) is a global minimum. This concludes the proof, establishing that the critical points of R S (θ) are either strict saddle points or global minima.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:EQUATIONAssume θ c is a saddle point for which ∇R S (θ c ) = 0. We can apply a first-order Taylor expansion to the right-hand side of Eq. ( 16), yielding:EQUATIONwhere ∇ 2 R S (θ c ) represents the Hessian matrix. Given that ∇R S (θ c ) = 0, the gradient flow dynamics around θ c can be approximated as:EQUATIONwhere H := -∇ 2 R S (θ c ). Eq. ( 18) is a classic linear ordinary differential equation, with the solution:EQUATIONThe 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. ( 19) elucidates that the dynamics near the critical point θ c are predominantly dictated by the properties of H, a real symmetric matrix in R 2d 2 ×2d 2 . Its eigendecomposition is given by:EQUATION)where Λ is a diagonal matrix and Q is an orthogonal matrix. Let λ 1 > λ 2 > • • • > λ s ∈ R denote the eigenvalues of H, and let q i1 , q i2 , • • • , q ili represent the eigenvectors corresponding to λ i .Given that λ 1 > λ 2 , the ratio e λ1t /e λit for i > 1 grows exponentially fast. Consequently, near θ c , the evolution of the system is primarily driven by the eigenvectors q 11 , q 12 , • • • , q 1l1 associated with the largest eigenvalue λ 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 θ c to be a saddle point. Consider:EQUATIONHere, θ 0 and θ 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 θ= -∇ 2 R S (θ c )(θ -θ c ), we denote H := -∇ 2 R S (θ c ).Let λ 1 ∈ R be the largest eigenvalue of H, with corresponding eigenvectors q 11 , q 12 , • • • , q 1l1 . Denote c j = ⟨θ 0θ c , q 1j ⟩, ∀1 ≤ j ≤ l 1 . Assume there exists j such that c j ̸ = 0, then, given any ε > 0, there exists a t 0 > 0 such that for all t ≥ t 0 , the following holds:EQUATION)This means that, as time progresses, the direction of the parameter evolution increasingly aligns with the dominant eigenvectors of H.This proposition is a consequence of the fact that the solution of the differential equation is given by θ(t) = e Ht θ(0), and as t tends to infinity, the term corresponding to the dominant eigenvalue in the matrix exponential e Ht becomes dominant. Therefore, near the saddle point, we have θ(t) ≈ e λ1t l1 j=1 c j q 1j + θ c .Proof. The solution of the ordinary differential equation ( 21) is θ(t) = e tH (θ 0θ c ) + θ c . Here, H = -∇ 2 R S (θ c ) is a real symmetric matrix, which can be diagonalized. Let λ 1 > λ 2 > • • • > λ s be the eigenvalues of H, and let q i1 , q i2 , • • • , q ili be the eigenvectors corresponding to λ i . We can then express θ(t) as:EQUATIONNext, we analyze the norm of the relative difference between θ(t) and a term dominating its growth:EQUATION)We define C = s i=2 li j=1 ⟨θ 0θ c , q ij ⟩q ij l1 j=1 ⟨θ 0θ c , q 1j ⟩q 1j. By choosing t 0 = log C ε λ 1λ 2 , we ensure that for all t > t 0 , the following condition is met:EQUATIONProp. 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 θ(t) = -∇R S (θ), 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 θ(t 0 ) moves away from θ 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 θ(t) = -∇R S (θ), we use φ(θ 0 , t) to denote the value of θ(t) in the case of θ(0) = θ 0 . At a critical point θ c , we denote the negative Hessian asH := -∇ 2 R S (θ c ).Let λ 1 ∈ R be the largest eigenvalue of H, with corresponding eigenvectors q 11 , q 12 , • • • , q 1l1 . Denote c j = ⟨θ 0θ c , q 1j ⟩, ∀1 ≤ j ≤ l 1 , and v 1 = l1 j=1 c j q 1j . Assume there exists j such that c j ̸ = 0. Let z α (t) := φ θ c + αv 1 , t + 1 λ1 log 1 α for every α > 0, then z(t) := lim α→0 z α (t) exists and is also a solution of the given dynamics, i.e., z(t) = φ(z(0), t). Furthermore, ∀t ∈ R, there exists a constant C > 0 such thatφ θ c + αθ 0 , t + 1 λ 1 log 1 α -z(t) 2 ≤ Cα λ 1 -λ 2 2λ 1 -λ 2for every sufficiently small α, where λ 1λ 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 λ 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(α λ 1 -λ 2 2λ 1 -λ 2 ). Therefore, Proposition A.6 holds in this case. Now, if the eigenspace corresponding to λ 1 is not one-dimensional, we denote c j = ⟨θ 0θ c , q 1j ⟩, ∀1 ≤ j ≤ l 1 , and v 1 = l1 j=1 c j 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(αλ 1 -λ 2 2λ 1 -λ 2 ).Therefore, Proposition A.6 holds in this case as well.Suppose the dominant directions fulfill specific conditions, such as any combination c 1 q 11 + c 2 q 12 + • • • + c l1 q 1l1 , leading to rank 1 model parameters (A, 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 θ c = (A c , 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:R S (θ) = E S ℓ (f (x, θ), f * (x)) , ∇ θ R S (θ) = E S ∇ℓ (f (x, θ), f * (x)) ⊤ ∇ θ f θ (x), = d 2 i=1 E S ∂ i ℓ (f θ , f * ) ∇ θ (f θ ) i , = d 2 i=1 E S (f θ -f * ) i ∇ θ (f θ ) i ,where∂ i ℓ (f θ , f * ) is the i-th element of ∇ℓ (f (x, θ), f * (x)), and (f θ ) i is the i-th element of the vectorization of f θ .For the Hessian matrix H S (θ), we haveH(θ) := ∇ θ ∇ θ R S (θ) = d 2 i=1 E S ∇ θ (∂ i ℓ (f θ , f * )) ∇ θ (f θ ) i + d 2 i=1 E S ∂ i ℓ (f θ , f * ) ∇ θ ∇ θ ((f θ ) i ) = d 2 i,j=1 E S ∂ ij ℓ (f θ , f * ) ∇ θ (f θ ) i ∇ θ (f θ ) j ⊤ + d 2 i=1 E S ∂ i ℓ (f θ , f * ) ∇ θ ∇ θ ((f θ ) i ) , = d 2 i,j=1 ∇ θ (f θ ) i ∇ θ (f θ ) j ⊤ + d 2 i=1 E S (f θ -f * ) i ∇ θ ∇ θ ((f θ ) i ) ,where∂ ij ℓ (f θ , f * ) is the (i, j)-th element of ∇∇ℓ (f (x, θ), f * (x)).We define matrices H (1) (θ) and H (2) (θ) as follows:H (1) (θ) := d 2 i,j=1 ∇ θ (f θ ) i ∇ θ (f θ ) j ⊤ , H (2) (θ) := d 2 i=1 E S (f θ -f * ) i ∇ θ ∇ θ ((f θ ) i ) ,We further denote that H(θ) := H (1) (θ) + H (2) (θ).For matrix factorization model, the eigenvectors of H (2) has a special structure, as characterized by Lem. A.2.Lemma A.2 (Data-Independent Interleaved Structure of Eigenvectors of H (2) ). Let θ c = (A c , B c ) be any critical point of the matrix factorization model. If λ is an eigenvalue of H (2) (θ c ) ∈R 2d 2 ×2d 2 , then there exist at least d eigenvectors associated with λ. These d eigenvectors take theform v ⊗ e 1 , v ⊗ e 2 , • • • , v ⊗ e d ∈ R 2d 2, where v ∈ R 2d is a vector to be determined and e i is the unit vector representing the i-th column of the identity matrixI d ∈ R d×d .Proof. Let's denote the residual matrix at the critical point as δM = (A c B c -M ) Sx , where (A c , B c ) is a critical point. For the vectorized parameter θ c , by direct calculation the matrix H (2) := -∇ 2 R S (θ c ) can be formulated as a block matrix, with the diagonal blocks being 0. The specific format is as follows:EQUATION)Next, we compute the eigenvectors of H (2) . Let λ be an eigenvalue of H (2) . We need to verify thatv ⊗ e 1 , v ⊗ e 2 , • • • , v ⊗ e d ∈ R 2d 2are the eigenvectors of H (2) corresponding to λ, for a particular v ∈ R 2d yet to be determined. That is, we need to ensure that for all 1 ≤ i ≤ d, the equation (H (2) -λI 2d 2 )(v ⊗ e i ) = 0 has a non-zero solution for v. Notice thatEQUATION)Since λ is an eigenvalue of H (2) , the determinant of the matrix H (2) -λI 2d 2 equals zero. HenceEQUATION)Consequently, from Eq. ( 27), we conclude that there always exists a non-zero vector v ∈ R 2d such that (H (2) -λI 2d 2 )(v ⊗ e i ) = 0. Since v ̸ = 0, it is evident that v ⊗ e 1 , v ⊗ e 2 , • • • , v ⊗ e d ∈ Rv ⊗ e 1 , v ⊗ e 2 , • • • , v ⊗ e d ∈ R 2d 2, where v ∈ R 2d is a vector to be determined and e i is the unit vector representing the i-th column of the identity matrixI d ∈ R d×d .Proof. In the matrix factorization model, at the origin the gradient ∇ θ (f θ ) = 0 and thus H (1) (0) = 0 and the Hessian matrix reduces toH (2) (0), making H = -∇ 2 R S (0) = -H (2) (0).Let's denote the residual matrix at the origin as δM = (A c B c -M ) Sx , where at the origin (A c , B c ) = (0, 0). For the vectorized parameter θ c , the matrix H := -∇ 2 R S (0) can be formulated as a block matrix, with the diagonal blocks being 0. The specific format is as follows:EQUATION)By Lem. A.2, the proof is completed.Lemma A.3 (Eigenvectors Structure of H at Second-order Stationary Point). Let Ω denote an Ω k invariant manifold or sub-Ω k invariant manifold defined in Prop. A.1 and A.2, and consider a second-order stationary point θ c within Ω, i.e., ∇R S (θ c ) = 0 and θ ⊤ ∇ 2 R S (θ c )θ ≥ 0 for all θ ∈ Ω. Then, the eigenvectors corresponding to the negative eigenvalues of the Hessian matrix H(θ c ) are contained within the span of the eigenvectors corresponding to the negative eigenvalues of H (2) (θ c ).Proof. Recall the definitions of H (1) (θ) and H (2) (θ) given by:H (1) (θ) := d 2 i,j=1 ∇ θ (f θ ) i ∇ θ (f θ ) j ⊤ , H (2) (θ) := d 2 i=1 E S (f θ -f * ) i ∇ 2 θ (f θ ) i .The Hessian matrixH(θ) at θ is H(θ) := H (1) (θ) + H (2) (θ).The manifold Ω is an affine subspace with orthogonal complement denoted by Ω ⊥ . Let H have the following block representation in the bases of Ω and Ω ⊥ :EQUATION)Since Ω is an invariant subspace under the gradient flow, we have Hθ ∈ Ω for all θ ∈ Ω, which implies that H 12 = 0. Since H is symmetry, we have H 21 = 0.Let λ < 0 be a negative eigenvalue of H := H(θ c ) with v as the corresponding eigenvector. Since H 11 is positive semi-definite, v must lie in Ω ⊥ .At a critical point θ c = (A c , B c ), by direct calculation, the gradient ∇ θ f θc can be structured as:EQUATION)where ⊗ denotes the Kronecker product.Note that ∇ θ (f θ * ) j is the j-th column of matrix ∇ θ f θc and it falls precisely within the defined Ω k invariant manifold or sub-Ω k invariant manifold Ω. Therefore, we have:EQUATION)which implies that v is orthogonal to the image of ∇ θ (f θc ), placing it in the null space of H (1) (θ c ).As a result, we have:H(θ c )v = H (1) (θ c ) + H (2) (θ c ) v = H (2) (θ c )v = λv.Thus, the eigenvector v of the Hessian H(θ c ), corresponding to the negative eigenvalue λ, is also an eigenvector of H (2) (θ c ), confirming that v is within the span of the eigenvectors of H (2) (θ c ).Proposition A.8. The linear combination of the eigenvectors ofEQUATIONBy Definition A.1, the data-independent invariant manifold generated by c isEQUATION) for some β = [β 1 , β 2 , • • • , β 2d ] ⊤ ∈ R 2d, and the vectorized parameter θ = vec((A, B)) ∈ R 2d takes the form β ⊗ c. Let β = v, and the proof is complete.Lemma A.4. Suppose α 1 , α 2 , • • • , α k+1 ∈ R d are linearly independent, the data- independent invariant manifold exhibits the property Ω k (α 1 , α 2 , • • • , α k ) + Ω 1 (α k+1 ) = Ω k+1 (α 1 , α 2 , • • • , α k+1 ). Proof. Assume that θ = (A, B) ∈ Ω k+1 (α 1 , α 2 , • • • , α k+1 ).Then A, B should adopt the form:EQUATION) Denote c i = β i γ i ∈ R 2d, then the vectorized parameters θ := vec(θ) can be expressed as:EQUATION)Since we know that Let φ(θ 0 , t) denote the value of θ(t) when θ(0) = θ 0 . Let Ω be a Ω k invariant manifold or sub-Ω k invariant manifold. Let θ c ∈ Ω be a critical point satisfying Assump. A.1 and A.2. Then, for randomly selected θ 0 , with probability 1 with respect to θ 0 , the limitk i=1 c i ⊗ α i ∈ Ω k (α 1 , α 2 , • • • , α k ) and c k+1 ⊗ α k+1 ∈ Ω 1 (α), it is straight- forward to validate that Ω k+1 (α 1 , α 2 , • • • , α k+1 ) = Ω k (α 1 , α 2 , • • • , α k ) + Ω 1 (α k+1 ).φ(θ c , t) := lim α→0 φ θ c + αθ 0 , t + 1 λ 1 log 1 α (37)exists and falls into an invariant manifold Ω k+1 . Here λ 1 is the top eigenvalue of negative Hessian -∇ 2 R S (θ c ).Proof. At the critical point θ c , we denote the negative Hessian as H := -∇ 2 R S (θ c ). Let λ 1 ∈ R be the largest eigenvalue of H, with corresponding eigenvectors q 11 , q 12 , • • • , q 1l1 . Denote c j = ⟨θ 0θ c , q 1j ⟩, ∀1 ≤ j ≤ l 1 , and v 1 = l1 j=1 c j q 1j . For a randomly selected θ 0 , with probability 1, there exists at least one j such that c j ̸ = 0.Consider the path z α (t) := φ θ c + αv 1 , t + 1 λ1 log 1 α for every α > 0. By Prop. A.6, the limit z(t) := lim α→0 z α (t) exists and satisfies the dynamics z(t) = φ(z(0), t).Furthermore, ∀t ∈ R, there exists a constant C > 0 such thatφ θ c + αθ 0 , t + 1 λ 1 log 1 α -z(t) 2 ≤ Cα λ 1 -λ 2 2λ 1 -λ 2for every sufficiently small α, where λ 1λ 2 > 0 is the eigenvalue gap.This implies that the limit lim α→0 φ θ c + αθ 0 , t + 1 λ1 log 1 α exists andφ(θ c , t) := lim α→0 φ θ c + αθ 0 , t + 1 λ 1 log 1 α = lim α→0 φ θ c + αv 1 , t + 1 λ 1 log 1 α .Assuming θ c ∈ Ω k and satisfies Assumps. A.1 and A.2, we aim to show the existence of a rank-(k+1) invariant manifold Ω k+1 containing θ c + αv 1 .Define the following matrices:H (1) (θ c ) := d 2 i,j=1 ∇ θc (f θc ) i ∇ θc (f θc ) j ⊤ , H (2) (θ c ) := d 2 i=1 E S (f θc -f * ) i ∇ 2 θc (f θc ) i .The Hessian matrix H(θ c ) at θ c can be expressed as H(θ c ) := H (1) (θ c ) + H (2) (θ c ), and we have H = -H(θ c ).Assump. A.1 and Lem. A.2 imply that there exist exactly d eigenvectors associated with the top eigenvalue λ 1 of -H (2) (θ c ). These eigenvectors are of the form v ⊗ e i for i = 1, . . . , d, where v ∈ R 2d is a vector to be determined and e i is the i-th standard basis vector in R d . By Assump. A.2 and Lem. A.3, the eigenvectors corresponding to λ 1 of H are contained within the span of the eigenvectors associated with the negative eigenvalues of H (2) (θ c ).Prop. A.8 ensures that the escaping direction v 1 lies within a rank-1 invariant manifold Ω 1 . Lem. A.4 then guarantees the existence of an invariant manifold Ω k+1 that includes θ c + αv 1 . Since Ω k+1 is invariant under the gradient flow, the trajectory φ θ c + αv 1 , t + 1 λ1 log 1 α remains within Ω k+1 .Finally, since Ω k+1 is a closed subspace, the limit φ(θ c , t) lies in Ωk+1 , concluding the proof.Consider the 2 × 2 matrix completion problem: M = 2 ⋆ ⋆ 2 . In this case, the two numbers on the diagonal are identical, which causes the maximum singular value of the residual matrix at the origin to be non-unique, violating Assump. 1. Consequently, the training process will jump directly from the rank 0 to the rank 2 invariant manifold, thereby missing the lowest rank solution of rank 1. This behavior is demonstrated in Fig. A2 , which shows experimental results for this scenario. Theorem A.5 (Minimum Rank). Let Ω denote an invariant as defined previously. Assume W t achieves a global minimum within each invariant manifold Ω k . If the limit W = lim α→0 W ∞ (αI) exists and is a global optimum with W (i, j) = M (i, j) for all (i, j) ∈ S x , thenEQUATIONProof. Consider the invariant manifold Ω k , which is defined as follows:Ω k := Ω k (α 1 , • • • , α k ) = {(A, B)|a i , b •,j ∈ span{α 1 , • • • , α k }, ∀1 ≤ i, j ≤ d},where a i denotes the i-th row of A, b •,j denotes the j-th column of B, and α 1 , . . . , α k are independent vectors that span the invariant subspace associated with Ω k .According to Thm. A.4, the training trajectory adheres to a hierarchical traversal across invariant manifolds. For any matrix C with rank(C) ≤ k, we will show that there always exists θ = (A, B) ∈ Ω k such that AB = C. Therefore, Ω k contains all matrices of rank k.In fact, Since rank(C) ≤ k, we can express C as a sum of k rank-one matrices:C = k i=1 u i v ⊤ iwhere u i and v i are column vectors. By the definition of Ω k , for any θ = (A, B) ∈ Ω k , each row of A and each column of B can be expressed as a linear combination of {α 1 , • • • , α k }:a i = k j=1 c ij α j b •,j = k i=1 d ij α iwhere c ij and d ij are scalars. We can write:A =     k j=1 c 1j α j . . . k j=1 c dj α j     , B = k i=1 d i1 α i • • • k i=1 d id α i .Now, we can express the product AB as:AB = k i=1 k j=1 d l=1 c li d jl α i α ⊤ jBy choosing appropriate values for c ij and d ij , we can make AB = C. This is possible because the outer products α i α ⊤ j span the same subspace as the rank-one matrices u i v ⊤ i in the expression of C. Therefore, for any matrix C with rank(C) ≤ k, there always exists θ = (A, B) ∈ Ω k such that AB = C.If the output matrix W t attains optimums within each Ω 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.9 Minimum Nuclear Norm Guarantee Lemma A.5 (Minimal Nuclear Norm Computation). Given a matrix M to be completed with observed diagonal entries, i.e., diag(M ) = v, the minimal nuclear norm solution among all possible completions is ∥v∥ 1 .Proof. The nuclear norm of a matrix is the dual of the spectral norm ∥ • ∥ 2 , defined as:∥A∥ * = max ∥X∥2≤1 ⟨A, X⟩.Given that ∥ diag(sign(v))∥ 2 ≤ 1, for any matrix A with diag(A) = v, it follows that:∥A∥ * ≥ ⟨A, diag(sign(v))⟩ = ⟨v, sign(v)⟩ = ∥v∥ 1 .Specifically, the nuclear norm of the diagonal matrix with v on its diagonal is ∥ diag(v)∥ * = ∥v∥ 1 , which establishes that the diagonal matrix with v is indeed a minimizer for the nuclear norm.Proposition A.9 (Minimum Nuclear Norm Guarantee in Diagonal Case). Consider the dynamics θ = -∇R S (θ), where θ(t) = (A(t), B(t)) and denote W t = A(t)B(t). If the observation data is diagonal, and if for a full rank initialization W 0 , the limit W = lim α→0 W ∞ (αW 0 ) exists and is a global optimum withEQUATIONProof. Without loss of generality, assume that M is a diagonal matrix given by:M =     µ 1 µ 2 . . . µ d     .By Lem. A.5, the minimal nuclear norm among all possible completions is |µ1 | + |µ 2 | + • • • + |µ d |.When the matrix to be completed is diagonal, the evolution of the i-th row of A is influenced only by the i-th column of B. Hence, the dynamics decouple into d independent parts, each equivalent to learning a scalar µ i . The learning process thus unfolds in d stages, with each stage passing through a critical point to learn a respective µ i .By Lem. A.2, the second term of the Hessian matrix can be expressed as:H (2) = 0 -δM ⊗ I d -δM ⊤ ⊗ I d 0 .According to Lem. A.2, the d eigenvectors of H (2) take the form v ⊗ e 1 , v ⊗ e 2 , . . . , v ⊗ e d ∈ R 2d 2 , where v ∈ R 2d is a vector to be determined and e i is the unit vector corresponding to the i-th column of the identity matrixI d ∈ R d×d . (i) Suppose µ 1 > µ 2 > • • • > µ 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 δM at any critical point remains a diagonal matrix. While starting to learn µ i from a critical point, direct calculation confirms that vector v = [e i , e i ] ⊤ ∈ R 2d . Escaping from each saddle point θ c , the trajectory θ(t)θ c approximates d i=1 c i (v ⊗ e i ), which satisfies B ,i = (A ⊤ ) ,i . Thus, learning a diagonal matrix M using the asymmetric model AB is equivalent to using a symmetric model AA ⊤ . The final outcome ensures that diag(AB) = diag(AA ⊤ ) = diag(M ).The nuclear norm of AA ⊤ equals the sum of its eigenvalues, which is precisely the trace of the matrix, and tr(AA⊤ ) = tr(M ) = µ 1 + µ 2 + • • • + µ d . Therefore, the nuclear norm of the learned matrix W = AB = AA ⊤ remains |µ 1 | + |µ 2 | + • • • + |µ d |. (ii) If some µ i < 0, assume without loss of generality that |µ 1 | > |µ 2 | > • • • > |µ n | > 0.While starting to learn µ i from a critical point, direct calculation confirms that v = [e i , sign(µ i )e i ] ⊤ ∈ R 2d . Escaping from each saddle point θ c , the trajectory θ(t)θ c approximates d i=1 c i (v ⊗ e i ), satisfying B ,i = sign(µ i )(A ⊤ ) ,i . Hence, AB = AA ⊤ Q, where Q is an orthogonal matrix given by:Q =     sign(µ 1 ) sign(µ 2 ) . . . sign(µ d )     .The final result ensures that diag(AB) = diag(AA ⊤ Q) = diag(M ), meaning diag(AA ⊤ ) = diag(QM ). The nuclear norm of AA ⊤ equals the sum of its eigenvalues, which is the trace of the matrix, and tr(AA⊤ ) = tr(QM ) = |µ 1 | + |µ 2 | + • • • + |µ d |.Since an orthogonal transformation does not change the nuclear norm of a matrix, the nuclear norm of the final learned matrixW = AB = AA ⊤ Q is still |µ 1 | + |µ 2 | + • • • + |µ d |.Theorem A.6 (Minimum Nuclear Norm Guarantee). Consider the dynamics θ = -∇R S (θ), where θ(t) = (A(t), B(t)), and let W t = A(t)B(t). If the observation graph associated with the matrix M to be completed is disconnected with complete bipartite components, and if for a full rank initialization W 0 , the limit W = lim α→0 W ∞ (αW 0 ) exists and is a global optimum withEQUATIONProof. Consider a matrix M ∈ R d×d composed of m connected components, with each component forming a complete bipartite subgraph. Since M is disconnected, it can be represented in a block diagonal form without loss of generality:M =     M 1 M 2 . . . M m     ,where each block M i ∈ R di×d ′ i , andm i=1 d i = d, m i=1 d ′ i = d, representing the sum of the dimensions of the blocks.Each block 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 M forms a complete bipartite subgraph, the block M i is fully observed. We can do singular value decomposition (SVD) on each sub-block M i asM i = U i Σ i V ⊤i , where U i and V i are orthogonal matrices, and Σ i is a diagonal matrix with the singular values of M i .Construct block diagonal matrices U and V as follows:U =     U 1 U 2 . . . U m     , V =     V 1 V 2 . . . V m     .This leads to the diagonal matrix:U M V ⊤ =     Σ 1 Σ 2 . . . Σ m     =     µ 1 µ 2 . . . µ d     .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.,|µ 1 | + |µ 2 | + • • • + |µ d |.Consider an incomplete matrix M whose associated observational graph is divided into m connected components, denoted by L 1 , L 2 , . . . , L m . For each component L p , we define Sx as the subset of observed indices within L p , where 1 ≤ p ≤ m and S 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|∃j : (i, j) ∈ S Lp x }, C p = {j|∃i : (i, j) ∈ S Lp x }.Here, R p includes the row indices and C p includes the column indices of the entries observed in L p .Define A Lp and B Lp as the submatrices of A and B corresponding to R p and C p , respectively. The evolution of A Lp is influenced only by B Lp . Thus, the dynamics decouple into m independent parts, each equivalent to learning a fully observed matrix M i . Accordingly, we partition A and B into m blocks:A =    A 1 . . . A m    , B = [B 1 • • • B m ] ,whereA i = A Li and B j = B Lj . Denote L i = ∥A i B i -M i ∥ 2 2 . The overall loss function L = m i=1 L i can be decomposed into m independent parts. Performing orthogonal transformations Ãi = U ⊤ i A i and Bi = B i V i , we obtain a diagonal loss Li = ∥ Ãi Bi -Σ i ∥ 2 2 for 1 ≤ i ≤ C.Since gradient descent is the steepest descent in the l 2 norm and orthogonal transformations preserve this norm, the dynamics of optimizing Li are equivalent to those of optimizing L i .Without loss of generality, assumeµ 1 > µ 2 > • • • > µ d > 0.Otherwise, as with Prop. A.9, a sign orthogonal transformation Q can be applied without changing the nuclear norm.By Prop. A.9, the learning result for a diagonal matrix implies Bi = Ã⊤ i , which meansB i V i = U ⊤ i A ⊤ i . The final learning result Ã =    Ã1 . . . Ãm    , B = B1 • • • Bm , satisfies B = Ã⊤ .The result ensures that diag( Ã B) = diag( Ã Ã⊤ ) = diag(Σ). The nuclear norm of Ã Ã⊤ equals the sum of its eigenvalues, which is the trace of the matrix, and Tr(Ã Ã⊤ ) = Tr(Σ) = |µ 1 | + |µ 2 | + • • • + |µ d |.Since orthogonal transformations do not alter the nuclear norm of a matrix, the nuclear norm ofW = AB is also |µ 1 | + |µ 2 | + • • • + |µ d |, concluding the proof.In this section, we present the supplementary experiments mentioned in the main text and the details of experiments.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.In the connectivity experiments corresponding to Fig. 1 , we explore the behavior of randomly generated 4 × 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: 2rdr 2 , which meets the threshold for exact reconstruction, 2rdr 2 -1, which is just below the threshold, and 2rdr 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. In the connectivity experiments depicted in Figs. 2(c-d ), we examine the behavior of randomly generated matrices of size 4 × 4 and 10 × 10 with a rank of 1. The matrices are sampled at a size of 2rdr 2 , which corresponds to the threshold for exact reconstruction. We evaluate two connected and one disconnected sampling patterns. A similar approach is taken for the 10 × 10 matrices. x i , and the standard deviation is calculated by:EQUATIONEQUATIONσ = 1 n -1 n i=1 (x i -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 × 4 matrix with rank 2, we selected a sample count less than the threshold of 2rdr 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. × × × × × × × × × × × × × × ×                                For a given sample size, there are different sampling models corresponding to connected or disconnected sampling. As shown in Fig. B1 , 7 observations are sampled, but different sampling positions affect connectivity or disconnection. To thoroughly study all possible cases, we examine all sampling cases of a 3 × 3 matrix completion, as illustrated in Figure 2(c) .For a 3 × 3 matrix, the sample size varies from 1 to 9. When using the matrix decomposition model f θ = AB for matrix completion, the dynamics obtained by exchanging rows or columns or transposing the matrix to be completed are equivalent. These three operations allow us to divide all sampling patterns equally.For sample size = 1, there is only one sampling pattern in the equivalent sense, and the observation matrix P is:P 1 = 1 0 0 0 0 0 0 0 0 .where 1 indicates that the position is observed and non-zero, and 0 means that the position is not observed or is 0. For sample size = 2, there are only 2 sampling patterns in the equivalent sense:P 1 = 1 1 0 0 0 0 0 0 0 , P 2 = 1 0 0 0 1 0 0 0 0For sample size = 3, there are only 4 sampling patterns in the equivalent sense:P 1 = 1 1 1 0 0 0 0 0 0 , P 2 = 1 1 0 1 0 0 0 0 0 , P 3 = 1 1 0 0 0 1 0 0 0 , P 4 = 1 0 0 0 1 0 0 0 1For sample size = 4, there are only 5 sampling patterns in the equivalent sense:P 1 = 1 1 1 1 0 0 0 0 0 , P 2 = 1 1 0 1 1 0 0 0 0 , P 3 = 1 1 0 0 1 1 0 0 0 , P 4 = 1 1 0 1 0 0 0 0 1 , P 5 = 1 1 0 0 0 1 0 0 1For sample size = 5, there are only 3 sampling patterns in the equivalent sense:P 1 = 1 1 1 1 1 0 0 0 0 , P 2 = 1 1 1 1 0 0 1 0 0 , P 3 = 1 1 0 1 1 0 0 0 1For sample size = 6, there are only 4 sampling patterns in the equivalent sense:P 1 = 1 1 1 1 1 1 0 0 0 , P 2 = 1 1 1 1 1 0 0 0 1 , P 3 = 1 1 1 1 1 0 1 0 0 , P 4 = 1 1 0 0 1 1 1 0 1For sample size = 7, there are only 2 sampling patterns in the equivalent sense:P 1 = 1 1 1 1 1 1 1 0 0 , P 2 = 1 1 1 1 1 0 0 1 1For sample size = 8, there is only 1 sampling pattern in the equivalent sense:P 1 = 1 1 1 1 1 1 1 1 0For sample size = 9, there is only 1 sampling pattern in the equivalent sense:P 1 = 1 1 1 1 1 1 1 1 1 B.To validate the scalability of our findings, we extended our experiments to higher-dimensional matrices. We conducted tests on 20 × 20 matrices, employing both connected (Fig. B7 ) and disconnected (Fig. B8 ) sampling patterns, while monitoring rank evolution during training. Our results consistently corroborated the main findings:(i) Connected observations converged to optimal low-rank solutions.(ii) Disconnected observations yielded higher-rank solutions.(iii) The Hierarchical Invariant Manifold Traversal (HIMT) process was observed in both connected and disconnected scenarios. In the context of depth-3 matrix factorization models, we consider the functional form:f θ = ABC, where A, B, C ∈ R d×d .Figs. B9 and B10 suggest that even for a depth-3 model, the learning process exhibits a progression from low rank to high rank structures.Within the Transformer architecture, the matrix factorization component retains its significance. The attention mechanism is formalized as follows:f θ (X) = h i=1 softmax row XW Qi W ⊤ Ki X ⊤ √ d k XW Vi W Oi ,where the row-wise softmax operation is applied to the attention scores, and the sum is over the h different attention heads, with W Qi , W Ki , W Vi , W Oi representing the learnable weight matrices The attention module's ability to capture low-rank representations is reflected in the depth-2 matrix factorization model. As illustrated in Fig. B11 , the attention models consistently learns representations that evolve from lower to higher ranks. • The factors of variability that the error bars are capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, or overall run with given experimental conditions). • The method for calculating the error bars should be explained (closed form formula, call to a library function, bootstrap, etc.) • The assumptions made should be given (e.g., Normally distributed errors).• It should be clear whether the error bar is the standard deviation or the standard error of the mean. • It is OK to report 1-sigma error bars, but one should state it. The authors should preferably report a 2-sigma error bar than state that they have a 96% CI, if the hypothesis of Normality of errors is not verified. • For asymmetric distributions, the authors should be careful not to show in tables or figures symmetric error bars that would yield results that are out of range (e.g. negative error rates). Guidelines:• The answer NA means that there is no societal impact of the work performed.• If the authors answer NA or No, they should explain why their work has no societal impact or why the paper does not address societal impact.• Examples of negative societal impacts include potential malicious or unintended uses (e.g., disinformation, generating fake profiles, surveillance), fairness considerations (e.g., deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations. • The conference expects that many papers will be foundational research and not tied to particular applications, let alone deployments. However, if there is a direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that an improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that a generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster. • The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly, harms that could arise when the technology is being used as intended but gives incorrect results, and harms following from (intentional or unintentional) misuse of the technology. • If there are negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition to attacks, mechanisms for monitoring misuse, mechanisms to monitor how a system learns from feedback over time, improving the efficiency and accessibility of ML). 


Decision: The paper investigates the implicit regularization of matrix factorization for solving matrix completion problems, revealing that the connectivity of observed data has a strong influence on the impact of implicit regularization. It received a positive response from the reviewers, who unanimously found it novel and original, with sound numerical and theoretical analysis and a high-quality presentation.

Reviews: review: paper_summary: paper_summary: This paper systematically investigates the implicit regularization of matrix factorization for solving matrix completion problems. The authors find by empirical results 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 becomes more connected with increased observations. They 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. The optimization trajectory follows a Hierarchical Invariant Manifold Traversal (HIMT) process, generalizing the characterization of Li et al. (2020), whose proposed Greedy Low-Rank Learning (GLRL) algorithm equivalence only corresponding to the connected case. Regarding the minimum nuclear norm regularization, the authors establish conditions that provide guarantees closely aligned with the empirical findings, and they present a dynamic characterization condition that assures the attainment of the minimum rank solution.
strength_weakness: strength_weakness: - In some cases, an extremely small initialization is required, which may potentially impact the training speed. 
- The proofs of the theoretical results are quite complex and may be difficult to follow for readers who are not familiar with the mathematical background. A more intuitive explanation or proof sketch could help readers better understand the key ideas.
- It is suggested to provide justifications to demonstrate the reasonableness of Assumption 1.
questions: questions: No
limitations: limitations: N/A
paper_summary: paper_summary: This paper attempts to present a unified understanding of when and how matrix factorization models have different implicit regularization effects. Their key finding is that connectivity of the observed data plays an important role: (i) in the connected cases, the model learns the lowest-ranked solution and (ii) in the disconnected case, it seems to sort of depend, but generally does not find the minimum nuclear norm solution. They identify invariant manifolds in the loss landscape that guides the training trajectory from these lower ranked to higher ranked solutions. They support a few of their findings with theoretical results.
strength_weakness: strength_weakness: - I am unsure of how restrictive the assumptions used for the theoretical results are which are used to prove Theorem 2. The authors do make a remark about the assumptions, but it is still not immediately clear to me that these are fair assumptions to make.
- As the authors mention, it would be interesting (and important) future work to show how the connectivity of the observed data holds attainment across the invariant manifolds.
questions: questions: - I am trying to understand Figure 4(e): do the different dashed lines correspond to different variances of the initialization? In that case, I am assuming that the figure is trying to show that despite the initialization scales, the model always learns a sub-optimal solution in this disconnected case.
limitations: limitations: I have addressed a few limitations in the weaknesses section.
paper_summary: paper_summary: This paper studies the training dynamics of matrix factorisation for matrix completion, optimising the vanilla mean-squared-error loss function via gradient descent with small initialisation. The authors characterise the observation pattern of the underlying matrix via the connectivity of its associated bipartite graph, and show both empirically and theoretically that this plays an important role in the implicit regularisation of the learned solution.

The paper provides empirical evidence that in both the connected and disconnected cases, the row/column spaces of the factor matrices remain aligned at all times, and that the loss curves exhibit steady step-wise decline as the solutions traverse solutions of increasing rank. In the connected case, it is observed that the learned solutions of each rank are near optimal, however in the disconnected case, suboptimal solutions are learned.

These empirical observations and explained by theoretical results which examine the training dynamics of the continuous dynamical system following the gradient flow of the loss function with infinitesimal initialisation. The authors begin by defining a Hierarchical Inrinsic Invariant Manifold and Proposition 1 describes how the gradient flow behaves on these manifolds and that they form a hierarchy. Proposition 2 shows that in the special case that when the observation graph is disconnected, the gradient flow dynamics are restricted to sub-manifolds. Theorem 1 (which extends an analgous result for symmetric matrix factorization models in Li et al. (2020)) shows that the critical points of the gradient flow are either strict saddle points or global minima and Theorem 2 shows that the gradient flow traverses the hierarchy of Intrinsic Invariant Manifolds (or sub-manifolds in the disconnected case), a phenomenon coined "Hierarchical Invariant Manifold Traversal". The flagship results of the paper are Theorems 3 and 4 which show that in the connected case, the gradient flow converges to the minimum rank solution, and in the a special disconnected case (disconnected with complete bipartite components), the gradient flow converges to the minimum nuclear norm solution.
strength_weakness: strength_weakness: I find this work hard to fault. While the theoretical results could be criticised for being fairly complex, having spend some time understanding them, the payoff makes it absolutely worthwhile. I don't see that they could be made simpler.
questions: questions: - What the three bars for each matrix in Fig 2b? What is the thick vertical line?
- In Fig 3d, is there reason to expect the norms of the gradients to be on the same scale as the difference between the matrix at the saddle points and the optimal approximation for the corresponding rank? They seem to be close to the gradient norms at the saddle points. Is this a coincidence or is there a reason to expect this?
- In Fig 3c and 4e, loss trajectories for different sampling patterns are plotted. If I understand correctly, in Fig 3a for example, each line corresponds to a different entry of the matrix being unobserved. In this case is $\star = 1.2$, so the underlying rank is three, or does it take a different value? Is it the case that the matrix has rank three if and only if the unobserved entry takes a specific value?
- In Fig 4f, are you counting position indexes from zero? In Section 3 (lines 97 and 99), you count position indexes from one. 
- The red dot in Fig 4e would perhaps be best a little larger and plotted on top of the other lines. At present it is quite hard to see.
- Assumption 1 is equivalent to the top singular value of $-\delta \mathbf{M}$​ being unique.  This seems like a more intuitive way to write this assumption. Is there a reason you chose to write the assumption in terms of the eigenvalues of the symmetric dilation matrix?
- In 231, it is claimed that "in the connected case, at each level [the] model reaches an optimal solution". Is this implied by one of the theorems, by previous work, or is this an empirical observation?
limitations: limitations: The authors are upfront about the limitations of their work in the conclusions, which present exciting directions for future research.
paper_summary: paper_summary: This study empirical reveals that the connectivity of observed data significantly influences the implicit bias and identifies a hierarchy of intrinsic invariant manifolds in the loss landscape, providing a preliminary framework for understanding the mechanisms behind implicit regularisation in matrix factorisation models.
strength_weakness: strength_weakness: Please refer to the "Questions".
questions: questions: 1, Have you verified your findings on high dimensional matrix with experiments?\
2, What's the intuition of the augmented matrix? See $\textit{Line 110}$?\
3, What exactly "small initialisation" you defined in the experiments? See $\textit{Line 164}$\
Minor suggestion(s):\
$\textbf{A.}$ $Eq.3$ after $\textit{Line 117}$: You have mentioned $\theta(0) = \theta_0$, which indicates the initialisation (if I understand correctly). However, the initialisation method was mentioned far from here at $\textit{Line 36-37}$. It could be easier to follow if you can optimise this.\
$\textbf{B.}$ I expect a detailed example of connectivity and disconnectivity in the appendix (see $\textit{Line 561 - 564 and Fig. A1}$), which will make it easier for a wider audience to understand.
limitations: limitations: Yes.
paper_summary: paper_summary: The paper studies implicit regularization in matrix completion, a problem that estimates missing entries from a partially observed matrix. In particular, in this paper, matrix completion is formulated as an optimization problem of the form: $\min_{A, B} f(A,B) = \|M - AB\|_S^2$ where $S$ is the subset of observed entries and the function $\|\cdot\|_S$ is the Frobenius norm of a matrix after masking entries outside $S$ to zero. Using this formulation, the authors study the implicit regularization of the gradient flow associated with $f$ whose initial condition is close to zero. Their main results characterize the influence of connectivity - a property defined based on $S$ - on the dynamics and solutions obtained by the gradient flow. The results are supported both theoretically and numerically.
strength_weakness: strength_weakness: - The main concern is the interest of the problem: why should we care about data connectivity (in the sense of this paper)? I think the motivation of this study needs to be further justified.
questions: questions: 1. Major questions/suggestions:
- I wonder under which regime of $n$ (the size of the matrix) and $m$ (the number of observed entries) do we have more than/exactly one connected component? I think this question is important to justify the motivation of this work.
2. Minor questions/suggestions:
- The experiments in Section 4 need to be provided with more details: What is the algorithm used? If it is gradient descent, what is the learning rate? Or the solution is obtained by solving an ODE.
- Throughout the introduction, the term "connectivity" is repeatedly used but its description or definition is not provided (it is delayed until Section 3). The authors might consider describing/defining this term earlier to improve the readability.
- Line 612: The second inequality should be $\|w^T B\|\_F^2$ instead of $\|w^T B\|\_{S_x}^2$.
limitations: limitations: There are no particular limitations or potential negative impacts of this paper.

Comments: title: Thank you for your detailed response
Response to rebuttals
comment: Thank you for clearing up my confusion on the figure and the detailed response to the assumptions. I believe that my current assessment of the paper is correct and will keep my score, which recommends acceptance.
Thank you for the rebuttals. The authors' responses to my main concerns are generally satisfactory to me. Alongside the reviews of other reviewers, I am inclined to maintain my positive assessment of this submission.
I appreciate the authors' thorough response. After careful consideration, I have decided to raise my score.
Thanks the authors for their responses. I am quite satisfied with their answers, except for the second point: the threshold $(n - 1)^2 + 1$ seems too large for a graph to be connected. The classical random graph model yields a much better result, see [this link](https://mathoverflow.net/questions/60075/connectivity-of-the-erd%C5%91s-r%C3%A9nyi-random-graph). This is, however, just for information. I keep my score as it is and recommend acceptance.
I thank the authors for their clarifications and I continue to strongly support this paper.