Global Convergence of Gradient Descent for Asymmetric Low-Rank Matrix Factorization
Year: 2022

Abstract:
We study the asymmetric low-rank factorization problem: min U∈Rm×d,V∈Rn×d 1 2 ‖UV> −Σ‖F where Σ is a given matrix of sizem×n and rank d. This is a canonical problem that admits two difficulties in optimization: 1) non-convexity and 2) non-smoothness (due to unbalancedness of U and V). This is also a prototype for more complex problems such as asymmetric matrix sensing and matrix completion. Despite being non-convex and non-smooth, it has been observed empirically that the randomly initialized gradient descent algorithm can solve this problem in polynomial time. Existing theories to explain this phenomenon all require artificial modifications of the algorithm, such as adding noise in each iteration and adding a balancing regularizer to balance the U and V. This paper presents the first proof that shows randomly initialized gradient descent converges to a global minimum of the asymmetric low-rank factorization problem with a polynomial rate. For the proof, we develop 1) a new symmetrization technique to capture the magnitudes of the symmetry and asymmetry, and 2) a quantitative perturbation analysis to approximate matrix derivatives. We believe both are useful for other related non-convex problems.

text: We study the asymmetric low-rank factorization problem:
min U∈Rm×d,V∈Rn×d
1 2 ‖UV> −Σ‖2F
where Σ is a given matrix of sizem×n and rank d. This is a canonical problem that admits two difficulties in optimization: 1) non-convexity and 2) non-smoothness (due to unbalancedness of U and V). This is also a prototype for more complex problems such as asymmetric matrix sensing and matrix completion. Despite being non-convex and non-smooth, it has been observed empirically that the randomly initialized gradient descent algorithm can solve this problem in polynomial time. Existing theories to explain this phenomenon all require artificial modifications of the algorithm, such as adding noise in each iteration and adding a balancing regularizer to balance the U and V. This paper presents the first proof that shows randomly initialized gradient descent converges to a global minimum of the asymmetric low-rank factorization problem with a polynomial rate. For the proof, we develop 1) a new symmetrization technique to capture the magnitudes of the symmetry and asymmetry, and 2) a quantitative perturbation analysis to approximate matrix derivatives. We believe both are useful for other related non-convex problems.

heading: 1 Introduction
text: This paper studies the asymmetric low-rank matrix factorization problem:
min U∈Rm×d,V∈Rn×d
f (U,V) := 1
2 ‖UV> −Σ‖2F . (1)
where Σ ∈ Rm×n is a given matrix of rank d. While solving this optimization problem is not hard (e.g., using power method), in this paper, we are interested in using randomly initialized gradient descent to solve this problem:
Ut+1 = Ut + η(Σ−UtV>t )Vt; (2) Vt+1 = Vt + η(Σ−UtV>t )>Ut, (3)
35th Conference on Neural Information Processing Systems (NeurIPS 2021).
where η > 0 is the learning rate and U0,V0 are randomly initialized according to some distribution. Empirically, gradient descent with a constant learning rate can efficiently solve this problem (see, e.g., Figure 1 in Du et al. [2018]). Somehow surprisingly, there is no global convergence proof of this generic algorithm, let alone convergence rate analysis. The main difficulties are 1) the problem is non-convex and 2) this problem is not smooth with respect to (U,V) because the magnitudes of them can be highly unbalanced.
To motivate the study of gradient descent for this optimization problem, we note that this is a prototypical optimization problem that illustrates the gap between practice and theory. In particular, the prediction function UV> is homogeneous: if we multiply a factor by a scalar c and divide another factor by c, the prediction function remains the same. This homogeneity also exists in deep learning models. Therefore, progress made in understand (1) can further help us gain understanding on other non-convex problems, such as asymmetric matrix sensing, asymmetric matrix completion, and deep learning optimization. We refer readers to Du et al. [2018] for more discussions.
For Problem (1), Du et al. [2018] showed gradient flow (gradient descent with the step size η → 0), U̇ = ( Σ−UV> ) V and V̇ = ( Σ−UV> )> U,
converges to the global minimum but no rate was given. Here U : R → Rm×d,U = U(t) and V : R → Rn×d,V = V(t) are an integral curve over manifold Rm×d × Rn×d, and U̇ := dUdt , V̇ := dVdt . Key in their proof is an invariance maintained by gradient flow: d dt ( U>U−V>V ) = 0. This invariance implies that if initially the difference between the magnitudes of U0 and V0 is small, then the difference remains small. This in turn guarantees the smoothness on the gradient flow trajectory. Du et al. [2018] further uses a geometric result (all saddle points in the objective function are strict and all local minima are global minima [Ge et al., 2015, 2017b, 2016, Li et al., 2019b]), and then invokes the stable manifold theorem to show the global convergence of gradient flow [Lee et al., 2016, Panageas and Piliouras, 2016].
However, to prove a polynomial convergence rate, the approach that solely relies on the geometry will fail because there exists a counter example [Du et al., 2017]. Furthermore, for gradient descent with η > 0, the key invariance no longer holds.1
Du et al. [2018] also studied gradient descent with decreasing step sizes ηt = O ( t−1/2 ) , and obtained an “approximate global optimality result": if the magnitude of the initialization is O (δ), then gradient descent converges to a δ-optimal solution, i.e., this result does not establish that gradient descent converges to a global minimum. And again, there was no convergence rate. Furthermore, their result crucially relies on ηt is of order O ( t−1/2 ) to ensure the second order term does not diverge and thus does not apply to gradient descent with a constant learning rate.
Some previous works, e.g., Ge et al. [2015], Jin et al. [2017], modified the gradient descent algorithm to the perturbed gradient descent algorithm by adding an isotropic noise at each iteration, which can help escape strict saddle points and bypass the exponential lower bound in Du et al. [2017]. To deal with the non-smooth problem, they also added a balancing regularization term [Park et al., 2017, Tu et al., 2016, Ge et al., 2017a, Li et al., 2019b], 18‖U
>U−V>V‖2F to the objective function to ensure balancedness between U and V throughout the optimization process. With these two modifications, one can prove a polynomial convergence rate. However, experiments suggest that the isotropic noise and the balancing regularizer may be proof artifacts, because vanilla gradient descent applies to the original objective function (1) without any regularizer finds a global minimum efficiently. From a practical point of view, one does not want to add noise or additional regularization because it may require more hyper-parameter tuning.
The only global quantitative analysis for randomly initialized gradient is by Du et al. [2018] who proved the global convergence rate for the case where Σ has rank 1, and U and V are two vectors. In this case, one can reduce the problem to the dynamics of 4 variables, which can be easily analyzed. Unfortunately, it is very difficult to generalize their analysis to the general rank setting. Recent work also studies this case with noise injected to the gradient update [Liu et al., 2021].
In this paper, we develop new techniques to overcome the technical difficulties and obtain the first polynomial convergence of randomly initialized gradient descent for solving the asymmetric low-rank
1While the invariance can still hold approximately in some way, characterizing the approximation error is highly non-trivial, and this is one of our key technical contributions.
matrix factorization problem. Most importantly, our analysis is completely different from existing ones: we give a thorough characterization of the entire trajectory of gradient descent.
Before presenting our main results, we emphasize that the goal of this paper is not to provide new provably efficient algorithms to solve Problem (1), but to provide a rigorous analysis of an intriguing and practically relevant phenomenon on gradient descent. This is of the same flavor as the recent breakthrough on understanding Burer-Moneiro method for solving semidefinite programs [Cifuentes and Moitra, 2019].

heading: 1.1 Main Results
text: Our main result is below. Theorem 1.1. Suppose each entry of U0 and V0 are initialized using Gaussian distribution with mean 0 and variance ε2, where ε = Õ ( σd√
d3σ1(m+n)
) .2 Then there exists Ttotal(δ, η) =
O (
1 ησd ln dσdε + 1 ησd ln σdδ
) such that for any δ > 0 and learning rate η = O ( σdε 2
dσ31
) , we have that
with high probability over the initialization, when t > Ttotal(δ, η),f (Ut,Vt) ≤ δ.
Here, σ1 and σd are the largest and the smallest singular values of Σ, respectively. Notably, in sharp contrast to the result in Du et al. [2018], which requires the initialization depends on δ, our initialization does not depend on the target accuracy. To our knowledge, this is the first global convergence result for gradient descent in solving Problem (1). Furthermore, we give a polynomial rate. The first term in Ttotal(δ, η) represents a warm-up phase and the second term represents the local linear convergence phase, which will be clear in the analysis sections. On the other hand, while we believe Ttotal(δ, η) is nearly tight, our requirement for η is loose. An interesting future direction is further relax this requirement.
Now by taking η → 0, we have the following corollary for gradient flow. Corollary 1.2. Given δ > 0, there exists T = O (
1 σd ln dσdε + 1 σd ln σdδ
) , such that with high
probability over the initialization, for all t ≥ T , we have f (Ut,Vt) ≤ δ.3
This is also the first convergence rate result of randomly initialized gradient flow for asymmetric matrix factorization. We note that our analysis on gradient flow is nearly tight. To see this, consider the ordinary differential equation ȧt = (σd − a2t )aT with initial point a0 > 0, then s = a2 has analytical solution st = σde 2σdt
e2σdt+ σd a20 −1 . Hence, to achieve a δ optimal solution, i.e. |σd − a 2| ≤ δ, we need σd ( σd a20 − 1 ) 1 δ ≤ e 2σdt + σd a20 − 1. Hence T = Θ ( 1 σd ln σdδ ) is necessary. Remark 1.2. We note that the weights will also converge. In section 3.3, we will prove that U and V are always bounded, which implies the gradient norm is bounded in terms of the loss. Since the loss function converges linearly, we have that the gradient norm converges to 0 linearly. The convergence of U and V follows from that the trajectory of (Ut,Vt) forms a Cauchy sequence.

heading: 1.3 Additional Related Work
text: Here we discuss additional related work. First, in the symmetric setting, e.g., minU ‖UU> −Σ‖2F , global convergence of randomly initialized gradient has been established in various settings [Jain et al., 2017, Li et al., 2018, Chen et al., 2019].4 However, as has been highlighted in Li et al. [2019a,b], Park et al. [2017], Tu et al. [2016], generalization to the asymmetric case is highly non-trivial. The major technical difficulty is to deal with the unbalancedness between U and V. To prevent this, additional balancing regularization is often added [Li et al., 2019b, Park et al., 2017, Tu et al., 2016, Sun and Luo, 2016], though empirically this has been shown to be unnecessary.
Another line of work showed one can first uses spectral initialization to find a near-optimal solution, then starting from there, gradient descent converges to an optimum with a linear rate [Tu et al.,
2Õ hides logarithmic terms. 3In gradient flow, t is a continuous time index. 4In Appendix B, we show the dynamics of gradient flow actually admits a closed form, and thus can be easily
analyzed.
2016, Zheng and Lafferty, 2016, Zhao et al., 2015, Bhojanapalli et al., 2016], though in practice random initialization often suffices. Recently, Ma et al. [2021] proved that if 1) the initialization is close to a global minimum and 2) U and V are balanced, then without adding additional balancing regularizer, gradient descent converges to a global minimum. Our stage two’s analysis is similar to theirs. However, their result cannot be directly applied to our analysis because they require a more stringent initialization than our stage two’s initial point. There is also other work studying the balancing effect in terms of the learning rate [Wang et al., 2021].
Notations. Throughout the paper, bold letters, e.g., U,V,Σ, are reserved for matrices running in the algorithm, non-bold letters, e.g., U, V,Σ are for our analysis. For a matrix W with rank r, denote σi(W ) as the ith largest singular value of W , ∀i ∈ [r]. Furthermore, if W is symmetric, denote λi(W ) as the ith largest eigenvalue of W . Let Σ ∈ Rm×n be a rank-d matrix with singular value σ1 ≥ · · · ≥ σd > 0, and define its conditional as κ := σ1σd . Our goal is to factorize Σ into UV >.

heading: 2 Main Difficulties and Technique Overview

heading: 2.1 A Reduction to Principal and Complement Spaces.
text: The starting point is the Polyak-Łojasiewicz condition: if we can establish that max {σd(Ut), σd(Vt)} is lower bounded by a considerable constant cmax, then we have ‖∇f(U,V)‖ ≥ cmax √ 2f(U,V), which implies a linear convergence. However, the dth singular values of U and V are not monotonic with t, and they can even decrease to an extremely small value.
To deal with this issue, we consider the following transformation. Let the singular value decomposition of Σ is Σ ≡ ΦΣ′Ψ>, where Φ ∈ Rm×m and Ψ ∈ Rn×n are unitary matrices, and Σ′ is diagonal matrix. Define U′t := Φ −1Ut and V′t = Ψ −1Vt. Then we can rewrite equations (2) and (3) as
U′t+1 = U ′ t + η(Σ ′ −U′tV′t > )V′t; (4)
V′t+1 = V ′ t + η(Σ ′ −U′tV′t > )>U′t. (5)
Hence, without loss of generality, we can assume Σ is a diagonal matrix with Σi,i = σi, ∀i ∈ [d], and Σi,j = 0 otherwise.
To proceed, we will analyse the principal space and the complement space separately. We denote the upper d × d matrix of U as U and denote the lower (m − d) × d matrix of U as J . Similarly, we define the upper d × d matrix of V as V and the lower (n − d) × d matrix as K. Define Σ := diag(σ1, · · · , σd). We can write out the dynamics of these matrices:
Ut+1 = Ut + η(Σ− UtV >t )Vt − ηUtK>t Kt; (6) Vt+1 = Vt + η(Σ− UtV >t )>Ut − ηVtJ>t Jt; (7) Jt+1 = Jt − ηJt(V >t Vt +K>t Kt); (8) Kt+1 = Kt − ηKt(U>t Ut + J>t Jt). (9)
Additional Notations Throughout this paper, we have some notation conventions. First of all, if we omit the subscript (iteration number) of a matrix, then it represents that this matrix at any iteration t. If some matrices without subscripts appear in the same equation, it means the equation holds for arbitrary iteration t, and the subscript for each matrix should be the same. For instance, if we define A := U+V2 , it means we define At := Ut+Vt 2 , ∀t ≥ 0. Similarly, we define B := U−V 2 .
Besides A,B, J and K, there are some other special capital letters used to represent specific matrices throughout this paper. Here is a list.
S = AA>; P = Σ−AA> +BB>; Q = AB> −BA>.
We define such S is because in symmetric case (B ≡ 0), although it is hard to find analytical solution for A in continuous time case, we do find analytical form for S, which contains all information about the singular values of A.
P and Q are just the symmetric and skew-symmetric part of matrix Σ − UV >. Hence the linear convergence of gradient descent is equivalent the linearly diminishing of P and Q by Pythagorean theorem. We will mention their definitions every time we use them.

heading: 2.2 Symmetrization
text: Our key observation is that although the singular values of U and V may not have monotonic property, the symmetrized matrix has this property. Formally, we define
A := U + V
2 and B := U − V 2 .
Here,A represents the magnitude in the principal space andB represents the magnitude of asymmetry. Empirically, we can observe that by choosing a sufficiently small learning rate η, we have two desired properties:
1. The smallest singular value of A is almost monotonically increasing; 2. The norms of B, J,K are almost monotonically decreasing.
The first property ensures we are learning the “signal", Σ, and the second property ensures the “noise" is disappearing. Therefore, if we can establish these two properties, we can prove the global convergence.

heading: 2.3 Two Stage Analysis
text: The analysis for asymmetric low rank case is divided into two stages. In the first stage we mainly focus on the increasing rate of σd(A). We will prove that in gradient descent method σd(Ad) increases exponentially fast to √ σ 2 and then ‖P‖op drops exponentially fast to σd 4 , while preserving ‖B‖F , ‖J‖op and ‖K‖op small. In the second stage, we will use the large σd(A) to lower bound the convergence speed of ‖Σ−UV>‖2F . We will prove that, once gradient descent starts at a point with small ‖P‖op, ‖B‖F , ‖J‖op and ‖K‖op, it will converge to global optimal point exponentially fast.

heading: 3 Proof Sketch of Theorem 1.1

heading: 3.1 Initialization
text: We first use a Gaussian distribution to generate matrices U, V, J,K element-wisely and independently5. By standard random matrix theory (Corollary 2.3.5 and Theorem 2.7.5 of Tao [2012]), we know that ∃c > 0, such that with high probability, the smallest singular value of U+V2 is larger than
1 c √ d , the largest singular value of U+V2 is smaller than c √ d, the Frobenius norm of B is less than cd and the operator norms of J and K are less than c √ max{m′, d} and c √
max{n′, d}, respectively, where m′ = m− d, n′ = n− d. The initializations U0, V0, J0,K0 are then scaled by ε where ε specified in Theorem 1.1.

heading: 3.2 Stage One: Warm-Up Phase
text: In this stage, we would like to prove the following theorem. Theorem 3.1. By choosing ε = Õ (
σd√ d3σ1(m+n)
) and η = O ( σdε 2
dσ31
) , we have that there exists
T0 = O (
1 ησd ln dσdε2
) , such that ∀t ≤ T0,
• ε 2 c2dI AtA > t 2Σ;
• ‖Bt‖F ≤ 2cdε; 5Strictly speaking, we cannot make any assumption on U, V, J,K since they need information of singular value decomposition of Σ. However, a random generation of U and V implies a random generation ofU, V, J,K because we use unitary transformations.
• σd(AT0) ≥ √ σd 2 ;
• σ1(PT0) ≤ σd4 ; • ‖Jt‖op ≤ cε √ max{m′, d}, ‖Kt‖op ≤ cε √ max{n′, d}.
We first give some intuitions about the five conditions in Theorem 3.1. The first condition represents the “signal" is properly bounded from below and above throughout stage one. The second condition shows the magnitude of asymmetry is small throughout stage one. We note that it is crucial to study the Frobenius norm of B instead of operator norm, because Frobenius norm admits a nice expansion for analysis. The third condition is an important one, which guarantees after T0 iterations, we have enough “signal" strength in the principal space. The fourth condition is a technical one, which represents the symmetric error is small after T0 iterations. The fifth condition represents the magnitude of the complement space remains small.
Proof Sketch. The proof of Theorem 3.1 is quite challenging and require new technical ideas and careful calculations, which we explain below.
To analyze the dynamic of σd(A), let us recall how it behaves in the continuous-time case. Our main idea is that, instead of analyzing A itself, we consider the symmetric matrix S := AA>. Then dS dt ≈ (Σ− S)S + S(Σ− S) plus some small perturbation terms about B, J and K.
If we only consider a differential equation Ṡ = (Σ − S)S + S(Σ − S), a well-known theorem (Theorem 12 in Lax [2007]) shows that if the singular values of S are different from each other, and ξ is the singular vector that Sξ = σd(S)ξ, then the derivative of σd(S) is exactly ξ>Ṡξ, which is lower bounded by 2(σd − σd(S))σd(S). To adapt it to discrete case, we prove the following lemma. Lemma 3.2. Suppose S,Σ ∈ Rd×d are two definite positive matrices, η > 0, and S′ = (I + η(Σ− S))S(I + η(Σ− S)). Suppose σ1(S) ≤ 2σ1, σd(Σ) ≥ σd and σ1(Σ) ≤ σ1. Define s = σd(S) and s′ = σd(S
′). Then ∀β ∈ (0, 1) and η ≤ β8σ1 ,
s′ ≥ (1 + η(σd − s))2s− 8 + 6β
1− β σ31η 2.
This lemma shows if we ignore perturbations from B, J and K, then for small η (when η2 is of smaller order than the first term), the least eigenvalue of S increases at a geometric rate.
However, there are also some small perturbation terms aboutB, J andK while doing analysis. ‖J‖op and ‖K‖op are easy to give an upper bound, since by (8) and (9), we know that by choosing small enough η, they are monotonically decreasing. However, the dynamic of B is highly non-trivial. After some careful calculations (cf. (19)), we find that the increasing rate of ‖B‖2F is related to the smallest eigenvalue of P := Σ−AA> +BB>: if max{0,−λd(P )} is small, then ‖B‖2F increases slowly. Now we would like to give a lower bound on λd(P ). Inspired by gradient flow case, P and S are almost complementary of each other, and their dynamic behaves similarly. Hence we have Ṗ ≈ −(Σ− P )P − P (Σ− P ) with some small perturbation terms about B, J and K. Hence we can use lemma 3.3 to give a lower bound in discrete case. Lemma 3.3. Suppose P,Σ ∈ Rd×d are two symmetric matrices, η > 0, and P ′ = (I − η(Σ − P ))P (I − η(Σ − P )). Suppose σ1(P ) ≤ 2σ1 and σdI Σ σ1I . Define p = λd(P ) and p′ = λd(P
′). Then ∀β ∈ (0, 1) and η ≤ β8σ1 , p′ ≥ { (1− ησd)2p− 8+6β1−β σ 3 1η
2, if p < 0; 0, if p ≥ 0.
Notice that we use B while analyzing P and use P while analyzing B. Hence, during the whole process, we need to bound both of them inductively.
Finally, once σd(A) increases to a relatively large amount, we can use it to prove that ‖P‖op will decrease exponentially fast to σd4 . One cannot simply prove that P converges to zero in this stage, since the perturbation term B will never converge to zero.
Below we give more details.

heading: 3.2.1 Assumptions
text: We make some assumptions on A and B in iterations t ≤ T0, where T0 will be defined at the end of subsubsection 3.2.4, and we will verify the assumptions in the end.
(1) ε 2 c2dI AA > 2Σ.
(2) The Frobenius norm of B is bounded by ebdε for some eb ≥ c, where eb will be determined later6. Hence its operator norm is also bounded by ebdε.

heading: 3.2.2 Dynamics on A, B and P
text: The dynamics on J and K is trivial, since by equations (8) and (9), i.e.
Jt+1 = Jt − ηJt(V >t Vt +K>t Kt); Kt+1 = Kt − ηKt(U>t Ut + J>t Jt),
we know that if we choose η ≤ 13σ1 , one can inductively proved that 0 V > t Vt +K > t Kt 3σ1I and 0 U>t Ut + J>t Jt 3σ1I by using the first two assumptions in subsection 3.2.1. And then it follows that the operator norms of J and K are monotonically decreasing in this stage.
However, it is non-trivial to prove that ‖B‖op keeps small. We will analyze the dynamics of A,B and P := Σ−AA> +BB> together inductively. First of all, from equations (6) and (7), we can write down the dynamics ofA := U+V2 andB := U−V 2 as following.
At+1 = At + η(Σ−AtA>t +BtB>t )At − η(AtB>t −BtA>t )Bt
−ηAt K>t Kt + J > t Jt
2 − ηBt K>t Kt − J>t Jt 2 ; (10)
Bt+1 = Bt − η(Σ−AtA>t +BtB>t )Bt + η(AtB>t −BtA>t )At
−ηAt K>t Kt − J>t Jt
2 − ηBt
K>t Kt + J > t Jt
2 . (11)
We can further calculate
Pt+1 = Pt − ηPt(Σ− Pt)− η(Σ− Pt)Pt + η2(PtPtPt − PtΣPt)− 2ηBtB>t Pt −2ηPtBtB>t − η(At + ηPtAt)C>t − ηCt(At + ηPtAt)> − η2CtC>t +η(Bt + ηPtBt)D > t + ηDt(Bt + ηPtBt) > + η2DtD > t (12)
where
Ct := −AtB>t Bt +BtA>t Bt −At K>t Kt + J > t Jt
2 −Bt K>t Kt − J>t Jt 2 ; (13)
Dt := +AtB > t At −BtA>t At −At K>t Kt − J>t Jt 2 −Bt K>t Kt + J > t Jt 2 , (14)
are two small perturbation terms.

heading: 3.2.3 Dynamics on A
text: Given (10), we can give a lower bound for the minimal singular value of At+1.
σd(At+1) ≥ σd(At + η(Σ−AtA>t )At) −η ∥∥∥∥BtB>t At −AtB>t Bt +BtA>t Bt −AtK>t Kt + J>t Jt2 −BtK>t Kt − J>t Jt2 ∥∥∥∥ op .
6We will show later that it is appropriate to choose eb = 2c.
For the first part, we could define St := AtA>t , and St+1 := (I + η(Σ − St))St(I + η(Σ − St)). Then according to lemma 3.2, by choosing β = 12 and η ≤ 1 16σ1 , we have7
σd(At+1) ≥ √ (1 + η(σd − σd(At)2))2σd(At)2 − 22σ31η2
−1.5 √
2σ1η(e 2 b + c 2)ε2(m+ n)d. (15)
For simplicity, we denote σd(At) by at, and define st = σd(St) = a2t . After some routine computations8, we can prove that it takes at most T1 := O (
1 ησd ln dσdε2
) iterations
to make at to at least √ σd 2 , and additional computations show that, if at is always bounded by √ 2σ1,
then once at becomes larger than √ σd 2 , it is always larger than √ σd 2 .

heading: 3.2.4 Dynamics on P
text: To bound Pt by equation (12), we need to first bound the norms of Ct and Dt by (13) and (14).
By simple triangle inequalities we have9
‖Ct‖op ≤ √ 2σ1(e 2 b + c 2)(m+ n)dε2; ‖Dt‖op ≤ 8σ1ebdε,
where the last inequality holds when choosing ε ≤ √ σ1ebd
c2(m+n) . Then we can conclude that
Pt+1 = (I − η(Σ− Pt))Pt(I − η(Σ− Pt)) + Et, (16)
where Et is a matrix with operator norm less than O(η2σ31 + ηe 2 bε 2(m+ n)dσ1 + η 2σ21e 2 bd 2ε2). By choosing ε ≤
√ σ1 ebd and η ≤ c 2ε2(m+n)d σ21 , we have ‖Et‖op ≤ O(ηe2bε2(m + n)dσ1). Further more,
by choosing β = 12 and η ≤ 1 16σ1 in lemma 3.3, we have
λd(Pt+1) ≥ max { (1− ησd)2λd(Pt)−O(ηe2bε2(m+ n)dσ1),−O(ηe2bε2(m+ n)dσ1) } .
Because P0 is initially positive, we know that
λd(Pt) ≥ −O(e2bε2(m+ n)dκ). (17)
This lower bound verifies the assumption that AtA>t 2Σ, since AtA>t = Σ − P + BtB>t Σ +O(e2bε 2(m+ n)dκ)I + e2bε 2dI 2Σ by choosing e2bε2 = O
( σd
(m+n)dκ
) .
On the other hand, we can also analyze the operator norm of P by using formula (16), since σd(At) ≥ √ σd 2 for t ≥ T1. This implies that
σ1(Pt+1) ≤ (
1− ησd 2
)2 σ1(Pt) +O(ηe 2 bε 2(m+ n)dσ1),
and it follows that σ1(Pt+T1) ≤ (
1− ησd 2
)2t σ1(PT1) +O(e 2 bε 2(m+ n)dκ). (18)
Inequality (18) shows that we only need at most T2 := O (
1 ησd
lnκ ) iterations after T1 to make
σ1(Pt) ≤ σd4 . Because ε 2 ≤ σd, we have the total number of iteration T0 := T1 + T2 = O (
1 ησd ln dσdε2
) .
7Please see (20) for the full steps for this inequality. 8Please see section D for details. 9Please see (21) for full steps for this inequality.

heading: 3.2.5 Dynamics on B
text: To verify the assumption about ‖B‖F made in subsection 3.2.1, we cannot simply use the equation (11), since the error term ‖(AB> −BA>)A‖op is approximately O(σ1‖B‖op), which will perturb the analysis seriously. Inspired by the continuous case that ˙‖B‖2F = 2 〈 B, Ḃ 〉 = Tr(B>PB) −
1 2‖Q‖ 2 F ≤ Tr(B>PB), where Q = AB> −BA> if we assume J = K = 0. In this inequality, we hide the term (AB> −BA>)AB> in − 12‖Q‖ 2 F and wipe it completely in our analysis.
Hence, for discrete case, we have the following inequality10,
‖Bt+1‖2F − ‖Bt‖2F ≤ −2ηλd(Pt)‖Bt‖2F + η‖B>t At‖F ‖K>t Kt − J>t Jt‖F +η2‖(Σ−AtA>t +BtB>t )Bt + (AtB>t −BtA>t )At
−At K>t Kt − J>t Jt
2 −Bt
K>t Kt + J > t Jt
2 ‖2F (19)
≤ O(ηe2bε2(m+ n)dκ)‖Bt‖2F +O(η √ σ1eb(m+ n)d 2ε3),
where the last equation is because we have chosen η = O ( σdε 2
dσ31
) and e2bε 2 = O ( σd (m+n)dκ ) .
By some routine calculations in section E, we have that by choosing ε = Õ (
σd√ σ1eb(m+n)
) , it is
appropriate to choose eb = 2c, so that ‖BT ‖2F ≤ e2bd2ε2, and induction holds.

heading: 3.3 Stage Two: Local Convergence Phase
text: We have proved in theorem 3.1 that the gradient descent achieved a pretty good point at T0, i.e. ‖BT0‖F ≤ 2cdε and σ1(PT0) ≤ σd4 . In this subsection, we will prove that start from this point, the gradient descent will converge linearly to the global optimal point. Then theorem 1.1 follows.
We will prove inductively on the following conditions:
(1) ‖B‖F = O( σd√σ1 );
(2) ∆t := ‖Σ− UT0+tV >T0+t‖op ≤ ( 1− ησd2 )t 2 5σd;
(3) σd(U), σd(V ) ≥ √ σd 2 .
Intuitively, the (1) guarantees the magnitude of asymmetry remains small; (2) guarantees that in the principal space, the error converges to 0 with a geometric rate; and (3) guarantees the “signal" in the principal space remains lower bounded.
Proof Sketch. First of all, it is easy to prove linear convergence of J and K by using assumption (3). Now we can verify the assumptions inductively.
(1) + (2)⇒ (3): We can prove σd(UV >) = Θ(σd). Because U −V is small, (3) follows by triangle inequality.
(1) + (3)⇒ (2): Consider continuous-time case, if we assume J = K = 0, the time derivative of Σ−UV > is −(Σ−UV >)V V >−UU>(Σ−UV >). Hence the convergence rate is lower bounded by σd(U) and σd(V ). Because the perturbation term J and K decreases exponentially, assumption (2) follows naturally. We transform this intuition to the discrete-time case.
(2) + (3) ⇒ (1): Again, we use (19) to show that the increasing rate of ‖B‖2F is bounded by ∆. Because ∆ decreases exponentially, ‖B‖2F cannot diverge to infinity, but increase by a poly(m,n, κ) factor. Then by taking ε sufficiently small can we verify the assumption (1).
The full proof is deferred to appendix.
10Please see (22) for full steps for this inequality.
Summary of Stage 2. To sum up, we have ‖Σ − UtV>t ‖2F = ‖Σ − UtV >t ‖2F + ‖UtK>t ‖2F + ‖JtV >t ‖2F + ‖JtK>t ‖2F , which can be further bounded by
‖Σ−UT0+tV>T0+t‖ 2 F ≤
( 1− ησd
2 )t 2 5 σd + 2(c 2 + c4)ε2σ1(m+ n)d ( 1− ησd 2 )2t ≤ C ( 1− ησd
2
)t σd,
for some universal constant C. Hence one only needs Tf := O ( ln σd δ
ησd
) iterations after T0 to achieve
an δ-optimal point.

heading: 4 Conclusion
text: This paper proved that randomly initialized gradient descent converges to a global minimum of the asymmetric low-rank matrix factorization problem with a polynomial convergence rate. This result explains the empirical phenomena observed in prior work, and confirms that gradient descent with a constant learning rate still enjoys the auto-balancing property as argued in Du et al. [2018].
We believe our requirement of the step size η is loose and a tighter analysis may improve the running time of gradient descent. Another interesting direction is to apply our techniques to other related problems such as asymmetric matrix sensing, asymmetric matrix completion and linear neural networks.11

heading: Acknowledgement
text: Simon S. Du gratefully acknowledges funding from NSF Award’s IIS-2110170 and DMS-2134106.

Review summaries:
This paper studies the problem of decomposing an m x n matrix of rank r into the product of an m x r matrix and an r x n matrix. While various methods exits for solving this problem, the goal of the paper is to provide theoretical guarantees for gradient descent with random initialization and fixed learning rate. The main result gives sufficient conditions, in terms of the dimensions (m,n), the rank r, and the condition number of the underlying matrix, under which the gradient descent method converges to a global optimum with a polynomial rate.   This work improves upon prior work which either: focused only the symmetric setting, does not give convergence results, or considered different initializations. The analysis and the conditions are based on two different stages. In the first stages,  the iterates converge slowly to get to a point `````near to an optimal solution. In the second stage, the iterate convergence quickly.   No numerical results are provided. 

This paper studies the performance of the gradient descent (GD) algorithm for a prototypical instance of the asymmetric matrix factorization (AMF) problem. The main result of the paper is establishing the polynomial convergence of GD for an ideal and noiseless, but unregularized instance of the asymmetric matrix factorization. To the best of my knowledge, the technical contributions of the paper are sound and correct.

This paper provides a proof of polynomial convergence rate whp of randomly initialized gradient descent on the non-convex problem of matrix factorization with the Frobenius norm. Notably, this problem was considered in previous works and only weaker results have been shown, e.g. no convergence rate, or polynomial convergence rate with the addition of noise and regularization.