LARGE LEARNING RATE TAMES HOMOGENEITY: CONVERGENCE AND BALANCING EFFECT
Year: 2022

Abstract:
Recent empirical advances show that training deep models with large learning rate often improves generalization performance. However, theoretical justifications on the benefits of large learning rate are highly limited, due to challenges in analysis. In this paper, we consider using Gradient Descent (GD) with a large learning rate on a homogeneous matrix factorization problem, i.e., minX,Y ‖A−XY ‖F. We prove a convergence theory for constant large learning rates well beyond 2/L, where L is the largest eigenvalue of Hessian at the initialization. Moreover, we rigorously establish an implicit bias of GD induced by such a large learning rate, termed ‘balancing’, meaning that magnitudes of X and Y at the limit of GD iterations will be close even if their initialization is significantly unbalanced. Numerical experiments are provided to support our theory.

text: Recent empirical advances show that training deep models with large learning rate often improves generalization performance. However, theoretical justifications on the benefits of large learning rate are highly limited, due to challenges in analysis. In this paper, we consider using Gradient Descent (GD) with a large learning rate on a homogeneous matrix factorization problem, i.e., minX,Y ‖A−XY >‖2F. We prove a convergence theory for constant large learning rates well beyond 2/L, where L is the largest eigenvalue of Hessian at the initialization. Moreover, we rigorously establish an implicit bias of GD induced by such a large learning rate, termed ‘balancing’, meaning that magnitudes of X and Y at the limit of GD iterations will be close even if their initialization is significantly unbalanced. Numerical experiments are provided to support our theory.

heading: 1 INTRODUCTION
text: Training machine learning models such as deep neural networks involves optimizing highly nonconvex functions. Empirical results indicate an intimate connection between training algorithms and the performance of trained models (Le et al., 2011; Bottou et al., 2018; Zhang et al., 2021; Soydaner, 2020; Zhou et al., 2020). Especially for widely used first-order training algorithms (e.g., GD and SGD), the learning rate is of essential importance and has received extensive focus from researchers (Smith, 2017; Jastrzebski et al., 2017; Smith, 2018; Gotmare et al., 2018; Liu et al., 2019; Li & Arora, 2019). A recent perspective is that large learning rates often lead to improved testing performance compared to the counterpart trained with small learning rates (Smith & Topin, 2019; Yue et al., 2020). Towards explaining the better performance, a common belief is that large learning rates encourage the algorithm to search for flat minima, which often generalize better and are more robust than sharp ones (Seong et al., 2018; Lewkowycz et al., 2020).
Despite abundant empirical observations, theoretical understandings of the benefits of large learning rate are still limited for non-convex functions, partly due to challenges in analysis. For example, the convergence (of GD or SGD) under large learning rate is not guaranteed. Even for globally smooth functions, very few general results exist if the learning rate exceeds certain threshold (Kong & Tao, 2020). Besides, popular regularity assumptions such as global smoothness for simplified analyses are often absent in homogeneous models, including commonly used ReLU neural networks.
This paper theoretically studies the benefits of large learning rate in a matrix factorization problem
min X,Y
1
2
∥∥A−XY >∥∥2 F , where A ∈ Rn×n, X, Y ∈ Rn×d. (1)
We consider Gradient Descent (GD) for solving (1): at the k-th iteration, we have
Xk+1 = Xk + h(A−XkY >k )Yk and Yk+1 = Yk + h(A> − YkX>k )Xk,
where h is the learning rate. Despite its simple formula, problem (1) serves as an important foundation of a variety of problems, including matrix sensing (Chen & Wainwright, 2015; Bhojanapalli et al., 2016; Tu et al., 2016), matrix completion (Keshavan et al., 2010; Hardt, 2014), and linear neural networks (Ji & Telgarsky, 2018; Gunasekar et al., 2018).
Problem (1) possesses several intriguing properties. Firstly, the objective function is non-convex, and critical points are either global minima or saddles (see e.g., Baldi & Hornik (1989); Li et al. (2019b); Valavi et al. (2020a); Chen et al. (2018)). Secondly, problem (1) is homogeneous in X and Y , meaning that rescaling X,Y to aX, a−1Y for any a 6= 0 will not change the objective’s value. This property is shared by commonly used ReLU neural networks. A direct consequence of homogeneity is that global minima of (1) are non-isolated and can be unbounded. The curvatures at these global minima are highly dependent on the magnitudes of X,Y . When X,Y have comparable magnitudes, the largest eigenvalue of Hessian is small, and this corresponds to a flat minimum; on the contrary, unbalanced X and Y give a sharp minimum. Last but not the least, the homogeneity impairs smoothness conditions of (1), rendering the gradient being not Lipschitz continuous unless X,Y are bounded. See a formal discussion in Section 2.
Existing approaches for solving (1) often uses explicit regularization (Ge et al., 2017; Tu et al., 2016; Cabral et al., 2013; Li et al., 2019a), or infinitesimal (or diminishing) learning rates for controlling the magnitudes of X,Y (Du et al., 2018; Ye & Du, 2021). In this paper, we go beyond the scope of aforementioned works, and analyze GD with a large learning rate for solving (1). In particular, we allow the learning rate h to be as large as approximately 4/L (see more explanation in Section 2), where L denotes the largest eigenvalue of Hessian at GD initialization. In connection to empirical observations, we provide positive answers to the following two questions:
Does GD with large learning rate converge at least for some cases of (1)? Does larger learning rate biases toward flatter minima (i.e., X,Y with comparable magnitudes)?
We theoretically show the convergence of GD with large learning rate for the two situations n = 1, d ∈ N+ or d = 1, n ∈ N+ with isotropic A. We also observe a, perhaps surprising, “balancing effect” for general matrix factorization (i.e., any d, n, and A), meaning that when h is sufficiently large, the difference between X and Y shrinks significantly at the convergence of GD compared to its initial, even if the initial point is close to an unbalanced global minimum. In fact, with a proper large learning rate h, ‖Xk − Yk‖2F may decrease by an arbitrary factor at its limit. The following is a simple example of our theory for n = 1 (i.e. scalar factorization), and more general results will be presented later with a precise bound for h depending on the initial condition and A.
Theorem 1.1 (Informal version of Thm.3.1 & 3.2). Given scalar A and initial condition X0, Y0 ∈ R1×d chosen almost everywhere, with learning rate h . 4/L, GD converges to a global minimum (X,Y ) satisfying ‖X‖2F + ‖Y ‖2F ≤ 2h . Consequently, its extent of balancing is quantified by ‖X − Y ‖2F ≤ 2h − 2A.
We remark that having a learning rate h ≈ 4/L is far beyond the commonly analyzed regime in optimization. Even for globally L-smooth objective, traditional theory requires h < 2/L for GD convergence and h = 1/L is optimal for convex functions (Boyd et al., 2004), not to mention that our problem (1) is never globally L-smooth due to homogeneity. Modified equation provides a tool for probing intermediate learning rates (see Hairer et al. (2006, Chapter 9) for a general review, and Kong & Tao (2020, Appendix A) for the specific setup of GD), but the learning rate here is too large for modified equation to work (see Appendix C). In fact, besides blowing up, GD with large learning rate may have a zoology of limiting behaviors (see e.g., Appendix B for convergence to periodic orbits under our setup, and Kong & Tao (2020) for convergence to chaotic attractors).
Our analyses (of convergence and balancing) leverage various mathematical tools, including a proper partition of state space and its dynamical transition (specifically invented for this problem), stability theory of discrete time dynamical systems (Alligood et al., 1996), and geometric measure theory (Federer, 2014).
The rest of the paper is organized as: Section 2 provides the background of studying (1) and discusses related works; Section 3 presents convergence and balancing results for scalar factorization problems; Section 4 generalizes the theory to rank-1 matrix approximation; Section 5 studies problem (1) with arbitrary A and its arbitrary-rank approximation; Section 6 summarizes the paper and discusses broadly related topics and future directions.

heading: 2 BACKGROUND AND RELATED WORK
text: Notations. ‖v‖ is `2 norm of a column or row vector v. ‖M‖F is the Frobenius norm of a matrix M .
Sharp and flat minima in (1) We discuss the curvatures at global minima of (1). To ease the presentation, consider simplified versions of (1) with either n = 1 or d = 1. In this case, X and Y become vectors and we denote them as x, y, respectively. We show the following proposition characterizing the spectrum of Hessian at a global minimum. Proposition 2.1. When n = 1, d ∈ N+ or d = 1, n ∈ N+ in (1), the largest eigenvalue of Hessian at a global minimum (x, y) is ‖x‖2 + ‖y‖2 and the smallest eigenvalue is 0.
Homogenity implies global minimizers of (1) are not isolated, which is consistent with the 0 eigenvalue. On the other hand, if the largest eigenvalue ‖x‖2 + ‖y‖2 is large (or small), then the curvature at such a global minimum is sharp (or flat), in the direction of the leading eigenvector. Meanwhile, note this sharpness/flatness is an indication of the balancedness between magnitudes of x, y at a global minimum. To see this, singular value decomposition (SVD) yields that at a global minimum, (x, y) satisfies ∥∥xy>∥∥2 F = σ2max(A). Therefore, large ‖x‖ 2
+ ‖y‖2 is obtained when |‖x‖ − ‖y‖| is large, i.e., x and y magnitudes are unbalanced, and small ‖x‖2 + ‖y‖2 is obtained when balanced.
Large learning rate We study smoothness properties of (1) and demonstrate that our learning rate is well beyond conventional optimization theory. We first define the smoothness of a function. Definition 2.2 (L-smooth). A function f ∈ C1 defined on RN is L-smooth if for all u1, u2 ∈ RN ,
‖∇f(u1)−∇f(u2)‖ ≤ L‖u1 − u2‖. (2)
If further f ∈ C2, then∇2f LI . Moreover, if we have (2) for u1, u2 ∈ X ⊆ RN , we call it locally L-smooth.
In traditional optimization (Nesterov, 2003; Polyak, 1987; Nesterov, 1983; Polyak, 1964; Beck & Teboulle, 2009), most analyzed objective functions often satisfy (i) (some relaxed form of) convexity or strong convexity, and (ii) L-smoothness. Choosing a step size h < 2/L guarantees the convergence of GD to a minimum by the existing theory (reviewed in Appendix F.3). Our choice of learning rate h ≈ 4/L (more precisely, 4/L0; see below) goes beyond the classical analyses.
Besides, in our problem (1), the regularity is very different. Even simplified versions of (1), i.e., with either n = 1 or d = 1, suffer from (i) non-convexity and (ii) unbounded eigenvalues of Hessian, i.e., no global smoothness (see Appendix F.2 for more details). As shown in Du et al. (2018); Ye & Du (2021); Ma et al. (2021), decaying or infinitesimal learning rate ensures that the GD trajectory stays in a locally smooth region. However, the gap between the magnitudes of X,Y can only be maintained in that case. We show, however, that larger learning rate can shrink this gap. More precisely, if initial condition is well balanced, there is no need to use large learning rate; otherwise, we can use learning rate as large as approximately 4/L, and within this range, larger h provides smaller gap between X and Y at the limit (i.e. infinitely many GD iterations).
Related work Matrix factorization problems in various forms have been extensively studied in the literature. The version of (1) is commonly known as the low-rank factorization, although here d can arbitrary and we consider both d ≤ rank(A) and d > rank(A) cases. Baldi & Hornik (1989); Li et al. (2019b); Valavi et al. (2020b) provide landscape analysis of (1). Ge et al. (2017); Tu et al. (2016) propose to penalize the Frobenius norm of ∥∥X>X − Y >Y ∥∥2 F
to mitigate the homogeneity and establish global convergence guarantees of GD for solving (1). Cabral et al. (2013); Li et al. (2019a) instead penalize individual Frobenius norms of ‖X‖2F + ‖Y ‖ 2 F. We remark that penalizing ‖X‖2F + ‖Y ‖ 2 F is closely related to nuclear norm regularization, since the variational formula for nuclear norm ‖Z‖2∗ = minZ=XY > ‖X‖ 2 F + ‖Y ‖ 2 F. More recently, Liu et al. (2021) consider using injected noise as regularization to GD and establish a global convergence (see also Zhou et al. (2019); Liu et al. (2022). More specifically, by perturbing the GD iterate (Xk, Yk) with Gaussian noise, GD will converge to a flat global optimum. On the other hand, Du et al. (2018); Ye & Du (2021); Ma et al. (2021) show that even without explicit regularization, when the learning rate of GD is infinitesimal,
i.e., GD approximating gradient flow, X and Y maintain the gap in their magnitudes. Such an effect is more broadly recognized as implicit bias of learning algorithms (Neyshabur et al., 2014; Gunasekar et al., 2018; Soudry et al., 2018; Li et al., 2020; 2021). Built upon this implicit bias, Du et al. (2018) further prove that GD with diminishing learning rates converges to a bounded global minimum of (1), and this conclusion is recently extended to the case of a constant small learning rate (Ye & Du, 2021). Our work goes beyond the scopes of these milestones and considers matrix factorization with much larger learning rates.
Additional results exist that demonstrate large learning rate can improve performance in various learning problems. Most of them involve non-constant learn rates. Specifically, Li et al. (2019c) consider a two-layer neural network setting, where using learning rate annealing (initially large, followed by small ones) can improve classification accuracy compared to training with small learning rates. Nakkiran (2020) shows that the observation in Li et al. (2019c) even exists in convex problems. Lewkowycz et al. (2020) study constant large learning rates, and demonstrate distinct algorithmic behaviors of large and small learning rates, as well as empirically illustrate large learning rate yields better testing performance on neural networks. Their analysis is built upon the neural tangent kernel perspective (Jacot et al., 2018), with a focus on the kernel spectrum evolution under large learning rate. Worth noting is, Kong & Tao (2020) also study constant large learning rates, and show that large learning rate provides a mechanism for GD to escape local minima, alternative to noisy escapes due to stochastic gradients.

heading: 3 OVERPARAMETERIZED SCALAR FACTORIZATION
text: In order to provide intuition before directly studying the most general problem, we begin with a simple special case, namely factorizing a scalar by two vectors. It corresponds to (1) with n = 1 and d ∈ N+, and this overparameterized problem is written as
min x,y∈R1×d
1 2 (µ− xy>)2, (3)
where µ is assumed without loss of generality to be a positive scalar. Problem (3) can be viewed as univariate regression using a linear two-layer neural network with the quadratic loss, which is studied in Lewkowycz et al. (2020) with atomic data distribution. Yet our analysis in the sequel can be used to study arbitrary univariate data distributions; see details in Section 6.
Although simplified, problem (3) is still nonconvex and exhibits the same homogeneity as (1). The convergence of its large learning rate GD optimization was previously not understood, let alone balancing. Many results that we will obtain for (3) will remain true for more general problems.
We first prove that GD converges despite of > 2/L learning rate and for almost all initial conditions: Theorem 3.1 (Convergence). Given (x>0 , y>0 ) ∈ (Rd×Rd)\B where B is some Lebesgue measure-0 set, when the learning rate h satisfies
h ≤ min {
4
‖x0‖2 + ‖y0‖2 + 4µ ,
1
3µ
} ,

heading: GD converges to a global minimum.
text: Theorem 3.1 says that choosing a constant learning rate depending on GD initialization guarantees the convergence for almost every starting point in the whole space. This result is even stronger than the already nontrivial convergence under small learning rate with high probability over random initialization (Ye & Du, 2021). Furthermore, the upper bound on h is sufficiently large: on the one hand, suppose GD initialization (x0, y0) is close to an unbalanced global minimum. By Proposition 2.1, we can check that the largest eigenvalue L(x0, y0) of Hessian∇2f(x0, y0) is approximately ‖x0‖2 + ‖y0‖2. Consequently, our upper bound of h is almost 4/L, which is beyond 2/L (see Section 2 for more details). On the other hand, we observe numerically that the 4/L upper bound is actually very close to the stability limit of GD when initialized away from the origin (see more details in Appendix E).
The convergence in Theorem 3.1 has an interesting searching-to-converging transition as depicted in Figure 1, where we observe two phases. In Phase 1, large learning rate drives GD to search for flat
regions, escaping from the attraction of sharp minima. After some iterations, the algorithm enters the vicinity of a global minimum with more balanced magnitudes in x, y. Then in Phase 2, GD converges to the found balanced global minimum. We remark that the searching-to-converging transition also appears in Lewkowycz et al. (2020). However, the algorithmic behaviors are not the same. In fact, in our searching phase (phase 1), the objective function does not exhibit the blow-up phenomenon. In our convergence phase (phase 2), the analysis relies on a detailed state space partition (see Line 192) due to nonconvex nature of (3), while the analysis in Lewkowycz et al. (2020) is akin to monotone convergence in a convex problem.
In comparison with the dynamics of small learning rate, we note that the searching phase (Phase 1) is vital to the convergence analysis. Meanwhile, the searching phase induces a balancing effect of large learning rate. The following theorem explicitly quantifies the extent of balancing. Theorem 3.2 (Balancing). Under the same initial condition and learning rate h as Theorem 3.1, GD for (3) converges to a global minimizer (x, y) satisfying
‖x‖2 + ‖y‖2 ≤ 2 h .
Consequently, its extent of balancing is quantified by
‖x− y‖2 ≤ 2 h − 2µ.
One special case of Theorem 3.2 is the following theorem, which states that no matter how close to a global minimum does GD start, if this minimum does not correspond to well-balanced norms, a large learning rate will take the iteration to a more balanced limit. We also demonstrate a sharp shrinkage in the distance between x and y. Corollary 3.3 (From ‘unbalanced’ to ‘balanced’). For any δ ∈ (0, µ), let the GD initialization satisfy
(x0, y0) ∈ { (u, v) : |uv> − µ| < δ, ‖u‖2 + ‖v‖2 > 8µ } \B,
where B is some Lebesgue measure-0 set. When the learning rate h satisfies h = 4‖x0‖2+‖y0‖2+4µ , the extent of balancing at the limiting point (x, y) of GD obeys
‖x− y‖2 < 1 2 ‖x0 − y0‖2 + 2µ.
Both Theorem 3.2 and Corollary 3.3 suggest that larger learning rate yields better balancing effect, as ‖x− y‖2 at the limit of GD may decrease a lot and is controlled by the learning rate. We remark that the balancing effect is a consequence of large learning rate, as small learning rate can only maintain the difference in magnitudes of x, y (Du et al., 2018).
In addition, the actual balancing effect can be quite strong with ‖xk − yk‖2 decreasing to be almost 0 at its limit under a proper choice of large learning rate. Figure 2 illustrates an almost perfect balancing case when h = 4‖x0‖2+‖y0‖2+4 ≈ 0.0122 is chosen as the upper bound. The difference ‖xk − yk‖ decreases from approximately 17.9986 to 0.0154 at its limit. Additional experiments with various learning rates and initializations can be found in Appendix A.
Technical Overview We sketch the main ideas behind Theorem 3.1, which lead to the balancing effect in Theorem 3.2. Full proof is deferred to Appendix G.
The convergence is proved by handling Phase 1 and 2 separately (see Fig.1). In Phase 1, we prove that ‖xk‖2 +‖yk‖2 has a decreasing trend as GD searches for flat minimum. We show that ‖xk‖2 +‖yk‖2 may not be monotone, i.e., it either decreases every iteration or decreases every other iteration.
In Phase 2, we carefully partition the state space and show GD at each partition will eventually enter a monotone convergence region. Note that the partition is based on detailed understanding of the dynamics and is highly nontrivial. Attentive readers may refer to Appendix G for more details. The combination of Phase 1 & 2 is briefly summarized as a proof flow chart in Figure 3.

heading: 4 RANK-1 APPROX. OF ISOTROPIC A (AN UNDER-PARAMETERIZED CASE)
text: Given insights from scalar factorization, we consider rank-1 factorization of an isotropic matrix A, i.e., A = µIn×n with µ > 0, d = 1, and n ∈ N+. The corresponding optimization problem is
min x,y∈Rn×1
1
2 ∥∥µIn×n − xy>∥∥2F . (4) Although similar at an uncareful glance, Problems (4) and (3) are rather different unless n = d = 1. First of all, Problem (4) is under-parameterized for n > 1, while (3) is overparameterized. More importantly, we’ll show that, when (x, y) is a global minimum of (4), x, y must be aligned, i.e., x = `y for some ` > 0. In the scalar factorization problem, however, no such alignment is required. As a result, the set of global minima of (4) is an n-dimensional submanifold embedded in a 2n-dimensional space, while in the scalar factorization problem the set of global minimum is a (2d− 1)-dimensional
submanifold — one rank deficient — in a 2d-dimensional space. We expect the convergence in (4) is more complicated than that in (3), since searching in (4) is demanding.
To prove the convergence of large learning rate GD for (4), our theory consists of two steps: (i) show the convergence of the alignment between x and y (this is new); (ii) use that to prove the convergence of the full iterates (i.e., x & y). Step (i) first: Theorem 4.1 (Alignment). Given (x0, y0) ∈ (Rn × Rn)\B, where B is some Lebesgue measure-0 set, when learning rate h ≤ min { 4
‖x0‖2+‖y0‖2+4 √ 7µ , 1 2 √ 7µ
} , the iterator (xk, yk) of GD at the
k-th iteration satisfies | cos(∠(xk, yk))| → 1 as k →∞.
Proof sketch. A sufficient condition for the convergence of | cos(∠(xk, yk))| is ‖xk‖2 ‖yk‖2 − (x>k yk)
2 → 0. To ease the presentation, let Uk = x>k yk, Vk = x>k xk, and Wk = y>k yk. By the GD update and some algebraic manipulation, we derive
Vk+1Wk+1 − U2k+1 = rk · (VkWk − U2k ),
where rk = (1 − h(Vk + Wk) + h2(VkWk − µ2))2. When k is sufficiently large, we can show a uniform upper bound on rk < 1− c for some constant c > 0. In this way, we deduce that VkWk−U2k will exponentially decay and converge to 0. More details are provided in Appendix H.
Theorem 4.1 indicates that GD iterations will converge to the neighbourhood of {(x, y) : x = `y, for some ` ∈ R\{0}}. This helps establish the global convergence as stated in Step (ii). Theorem 4.2 (Convergence). Under the same initial conditions and learning rate h as Theorem 4.1, GD for (4) converges to a global minimum.
Similar to the over-parametrized scalar case, this convergence can also be split into two phases where phase 1 motivates the balancing behavior with the decrease of ‖xk‖2 + ‖yk‖2, and phase 2 ensures the convergence. The following balancing theorem is thus obtained. Theorem 4.3 (Balancing). Under the same initial conditions and learning rate h as Theorem 4.1, GD for (4) converges to a global minimizer that obeys
‖x‖2 + ‖y‖2 ≤ 2 h ,
and its extent of balancing is quantified by
‖x− y‖2 ≤ 2 h − 2µ.
This conclusion is the same as the one in Section 3. A quantitatively similar corollary like Corollary 3.3 can also be obtained from the above theorem, namely, if x0 and y0 start from an unbalanced point near a minimum, the limit will be a more balanced one.

heading: 5 GENERAL MATRIX FACTORIZATION
text: In this section, we consider problem (1) with an arbitrary matrix A ∈ Rn×n. We replicate the problem formulation here for convenience,
min X,Y ∈Rn×d
1
2
∥∥A−XY >∥∥2 F . (5)
Note this is the most general case with n, d ∈ N+ and any square matrixA. Due to this generalization, we no longer utilize the convergence analysis and instead, establish the balancing theory via stability analysis of GD as a discrete time dynamical system.
Let µ1 ≥ µ2 ≥ · · · ≥ µn ≥ 0 be the singular values of A. Assume for technical convenience ‖A‖2F = ∑n i=1 µ 2 i being independent of d and n. We denote the singular value decomposition of A as A = UDV >, where U,D, V ∈ Rn×n, U, V are orthogonal matrices and D is diagonal. Then we establish the following balancing effect.
Theorem 5.1. Given almost all the initial conditions, for any learning rate h such that GD for (5) converges to a point (X,Y ), there exists c = c(n, d) > c0 with constant c0 > 0 independent of h, n, and d, such that (X,Y ) satisfies
c(‖X‖2F + ‖Y ‖ 2 F) <
2 h ,
and the extent of balancing is quantified by ∥∥X − (UV >)Y ∥∥2
F < 2ch −2 ∑min{d,n} i=1 µi, which means
∣∣ ‖X‖F − ‖Y ‖F ∣∣2 < 2ch − 2 √√√√min{d,n}∑
i=1
µ2i .
In particular, when d = 1, i.e., rank-1 factorization of an arbitrary A, the constant c equals 1.
We observe that the extent of balancing can be quantified under some rotation of Y . This is necessary, since for factorizing a general matrix A (which can be asymmetric), at a global minimum, X,Y may only align after a rotation (which is however fixed by A, independent of initial or final conditions). Figure 4 illustrates an example of the balancing effect under different learning rates. Evidently, larger learning rate leads to a more balanced global minimizer. Additional experiments with various dimensions, learning rates, and initializations can be found in Appendix A; a similar balancing effect is also shown for additional problems including matrix sensing and matrix completion there.
Different from previous sections, Theorem 5.1 builds on stability analysis by viewing GD as a dynamical system in discrete time. More precisely, the proof of Theorem 5.1 (see Appendix I for details) consists of two parts: (i) the establishment of an easier but equivalent problem via the rotation of X and Y , (ii) stability analysis of the equivalent problem.
For (i), by singular value decomposition (SVD), A = UDV >, where U,D, V ∈ Rn×n, U and V are orthogonal matrices, and D is a non-negative diagonal matrix. Let Xk = URk, Yk = V Sk. Then{
Xk+1 = Xk + h(A−XkY >k )Yk Yk+1 = Yk + h(A−XkY >k )>Xk
⇔ { URk+1 = URk + h(UDV > − URkS>k V >)V Sk
V Sk+1 = V Sk + h(UDV > − URkS>k V >)>URk ⇔ {
Rk+1 = Rk + h(D −RkS>k )Sk Sk+1 = Sk + h(D −RkS>k )>Rk .
Therefore, GD for problem (5) is equivalent to GD for the following problem
min R,S∈Rn×d
1
2
∥∥D −RS>∥∥2 F
and it thus suffices to work with diagonal non-negative A.
For (ii), here is a brief description of the idea of stability analysis: consider each iteration of GD as a mapping ψ from uk ∼ (Xk, YK) to uk+1 ∼ (Xk+1, Yk+1), where matrices X and Y are flattened
and concatenated into a vector so that ψ is a closed map on vector space R2dn. GD iteration is thus a discrete time dynamical system on state space R2dn given by
uk+1 = ψ(uk) = uk − h∇f(uk), where f is the objective function f(uk) = 12 ∥∥A−XkY >k ∥∥2F, and gradient returns a vector that collects all component-wise partial derivatives.
It’s easy to see that any stationary point of f , denoted by u∗, is a fixed point of ψ, i.e., u∗ = ψ(u∗). What fixed point will the iterations of ψ converge to? For this, the following notions are helpful: Proposition 5.2. Consider a fixed point u∗ of ψ. If all the eigenvalues of Jacobian matrix∇ψ(u∗) are of complex modulus less than 1, it is a stable fixed point. Proposition 5.3. Consider a fixed point u∗ of ψ. If at least one eigenvalue of Jacobian matrix ∇ψ(u∗) is of complex modulus greater than 1, it is an unstable fixed point.
Roughly put, the stable set of an unstable fixed point is of negligible size when compared to that of a stable fixed point, and thus what GD converges to is a stable fixed point for almost all initial conditions (Alligood et al., 1996). Thus, we investigate the stability of each global minimum of f (each saddle of f is an unstable fixed point of ψ and thus is irrelevant). By a detailed evaluation of ∇ψ’s eigenvalues (Appendix I), we see that a global minimum (X,Y ) of f corresponds to a stable fixed point of GD iteration if
∣∣∣1− ch(‖X‖2F + ‖Y ‖2F)∣∣∣ < 1, i.e., it is balanced as in Thm. 5.1.

heading: 6 CONCLUSION AND DISCUSSION
text: In this paper, we demonstrate an implicit regularization effect of large learning rate on the homogeneous matrix factorization problem solved by GD. More precisely, a phenomenon termed as “balancing” is theoretically illustrated, which says the difference between the two factors X and Y may decrease significantly at the limit of GD, and the extent of balancing can increase as learning rate increases. In addition, we provide theoretical analysis of the convergence of GD to the global minimum, and this is with large learning rate that can exceed the typical limit of 2/L, where L is the largest eigenvalue of Hessian at GD initialization.
For the matrix factorization problem analyzed here, large learning rate avoids bad regularities induced by the homogeneity between X and Y . We feel it is possible that such balancing behavior can also be seen in problems with similar homogeneous properties, for example, in tensor decomposition (Kolda & Bader, 2009), matrix completion (Keshavan et al., 2010; Hardt, 2014), generalized phase retrieval (Candes et al., 2015; Sun et al., 2018), and neural networks with homogeneous activation functions (e.g., ReLU). Besides the balancing effect, the convergence analysis under large learning rate may be transplanted to other non-convex problems and help discover more implicit regularization effects.
In addition, factorization problems studied here are closely related to two-layer linear neural networks. For example, one-dimensional regression via a two-layer linear neural network can be formulated as the scalar factorization problem (3): Suppose we have a collection of data points (xi, yi) ∈ R× R for i = 1, . . . , n. We aim to train a linear neural network y = (u>v)x with u, v ∈ Rd for fitting the data. We optimize u, v by minimizing the quadratic loss,
(u∗, v∗) ∈ arg min u,v
1
n n∑ i=1 ( yi − (u>v)xi )2 = arg min u,v 1 n (∑n i=1 xiyi∑n i=1 x 2 i − u>v )2 . (6)
As can be seen, taking µ = ∑n
i=1 xiyi∑n i=1 x 2 i
recovers (3). In this regard, our theory indicates that training of u, v by GD with large learning rate automatically balances u, v, and the obtained minimum is flat. This may provide some initial understanding of the improved performance brought by large learning rates in practice. Note that (6) generalizes to arbitrary data distribution of training a two-layer linear network with atomic data (i.e., x = 1 and y = 0) in Lewkowycz et al. (2020).
It is important to clarify, however, that there is a substantial gap between this demonstration and extensions to general neural networks, including deep linear and nonlinear networks. Although we suspect that large learning rate leads to similar balancing effect of weight matrices in the network, rigorous theoretical analysis is left as a future direction.

heading: ACKNOWLEDGMENTS
text: We thank anonymous reviewers and area chair for suggestions that improved the quality of this paper. The authors are grateful for partial supports from NSF DMS-1847802 (YW and MT) and ECCS-1936776 (MT).

heading: A ADDITIONAL EXPERIMENTS

heading: A.1 MORE RESULTS FOR MATRIX FACTORIZATION
text: In this section, we present more experiments with different choices of n, d, various initializations and scalings, and a broader range of learning rates. All these experiments verify our claim on the balancing effect of large learning rate, i.e., the shrinkage of the gap between the magnitudes of X and Y exists for general matrix factorization problems, i.e., for any choice of n, d ∈ N+.
Then we experiment with the general matrix factorization as a supplement of our theory in Section 5 where we only rigorously prove the balancing effect given the convergence of GD. In the following examples, we show that there indeed exist large learning rates that trigger the shrinkage of the gap between ‖Xk‖F and ‖Yk‖F and at the same time guarantee the convergence of GD. Figure 8, 9, and 10 correspond to the over-parameterized version of matrix factorization problem (5) with A ∈ R6×6 (asymmetric, generated by i.i.d. Gaussian) and X,Y ∈ R6×100. The initial conditions are randomly generated with (‖X0‖F , ‖Y0‖F) = (9, 1), (‖X0‖F , ‖Y0‖F) = (19, 1), and
(‖X0‖F , ‖Y0‖F) = (99, 1) respectively. Similarly, the learning rates are chosen from large to small as h0, 67h0, 5 7h0, 4 7h0, 3 7h0, 2 7h0 for the 1st-6th columns respectively where 6 7h0 (the 2nd column) is picked near the stability limit. We also similarly provide two examples of the under-parameterized version in Figure 11 and 12. The two examples correspond to problem (5) with A ∈ R100×100 (asymmetric, similarly generated as the previous one) and X,Y ∈ R100×3. Here we use a shifted loss which is to subtract the global minimum error.
As is shown in Figure 8, 9, 10, 11, and 12, we can observe the similar phenomenon as the scalar case, in the sense that (1) larger learning rate gives a smaller the gap between ‖Xk‖F and ‖Yk‖F in the limit, except for the overly large h that causes GD to diverge; (2) when the learning rate becomes sufficiently big, a two-phase pattern of convergence appears, which does not manifest in traditional optimization analysis and thus indicates that the learning rate is already larger than that permitted by traditional theory, and yet one still has convergence.
More precisely, the 2nd-4th columns are the large learning rate cases, where one observes a shrinkage of the unbalancedness and non-monotonicity of the loss at the beginning, while the right two columns corresponds to the small learning rates for which ‖Xk‖F and ‖Yk‖F barely change as k increases, and the decrease of loss is monotone. These are all evidence of the consistency between the most general case of matrix factorization in Section 5 and the special cases in Section 3 and 4.

heading: A.2 MATRIX SENSING AND MATRIX COMPLETION
text: As we demonstrated, the balancing effect is a nontrivial implicit bias created by large learning rate in GD, and this was rigorously established for the problem of matrix factorization. Matrix factorization already corresponds to a class of important problems as we consider arbitrary n and d; however, we feel balancing is an effect even more general, and thus we now demonstrate it empirically on two related additional problems, namely matrix sensing and matrix completion.
In Figure 13, we consider matrix sensing corresponding to the problem
min X,Y ∈Rn×d
1
2m m∑ i=1 (bi − 〈Ai, XY >〉)2.
Here Ai ∈ R100×100 are generated from element-wise i.i.d. Gaussian; bi ∈ R are generated from uniform distribution [0, 1]; m = 10; X,Y ∈ R100×6; 〈U, V 〉 = tr(V >U). The problem is solved via GD. Experiments in the figure correspond to learning rates chosen from large to small as h0, 45h0, 3 5h0, 2 5h0 where h0 (the 1st column) is picked near the stability limit. We can see from the figure that matrix sensing exhibits a similar balancing effect with matrix factorization that larger learning rate leads to more balanced norms between X and Y .
In Figure 14, we consider matrix completion corresponding to the problem
min X,Y ∈Rn×d
1
2 ∥∥PΩ(A−XY T )∥∥2F , where A ∈ R10×10 is a low-rank matrix with rank 2; X,Y ∈ R10×2; PΩ(U) = (Uij)(i,j)∈Ω + (0)(i,j)/∈Ω with sparsity−0.6 where sparsity=(number of non-zero elements)/(number of all elements). The problem is solved via GD. In the above mentioned figure, the learning rates are similarly chosen as the above matrix sensing example. Likewise, the decrease of learning rate results in the increase in the gap between the norms of X and Y .
Both the matrix sensing and matrix completion above hold the same homogeneity property between X and Y . The “balancing effect” that we proved for matrix factorization is also observed when and only when the learning rate h is large, in which case the norms of X and Y become close at the convergence of GD (and yes, GD still converges even though h is large enough such that the convergence is not monotone.
We feel the techniques invented and employed in this paper can extend to these two cases, but that is beyond the scope of this paper.

heading: B GD CONVERGING TO PERIODIC ORBIT
text: Consider the objective (1−xy)2/2. Take step size h = 1.9. Then GD can converge to periodic orbits with period 2, 3 and 4 respectively in Figure 15.

heading: C MODIFIED EQUATION FAILS FOR LARGE LEARNING RATES
text: Consider the objective function f(x, y) = (1− xy)2/2. Then the GD update is the following{ xk+1 = xk + h(1− xkyk)yk yk+1 = yk + h(1− ykxk)xk ⇒ ( xk+1 yk+1 ) = ( xk yk ) + h(−∇f(xk, yk)). (7)
By backward error analysis (Hairer et al., 2006, Chapter 9), the modified equation can better approximate GD than gradient flow and is defined as follows(
ẋ ẏ
) = −∇f(x, y)− h
2 ∇2f(x, y)∇f(x, y) +O(h2). (8)
Figure 16 shows the trajectories of the 1st order modified equation of (8) and GD (7) with initial condition x = 4, y = 10 and h = 0.026. The x-axis represents the time t and for GD, the time point for kth step is kh. We compare the absolute values of x and y of both methods due to the symmetry of the global minima xy = 1. As is shown in the figure, even if GD almost converges to the most balanced minimizer, the solutions of modified equation are still far away from each other, x ≈ 0.1 and y ≈ 9. Actually, large learning rate h fall outside the convergence domain of the modified equation which thus is not an appropriate tool for the analysis of norm balancing.

heading: D PROOF OF PROPOSITION 2.1
text: Proof of Proposition 2.1. This proposition is a direct consequence of Theorem F.2 and I.3.

heading: E OVER-PARAMETRIZED SCALAR DECOMPOSITION: STABILITY LIMIT
text: Consider the objective (1− xy)2/2. Choose the initial condition to be x0 = 20, y0 = 0.07 and use GD update. The upper bound of h in Theorem 3.1 is 4
x20+y 2 0+4µ
, where µ = 1. In Figure 17, when
h = 4 x20+y 2 0+4 (the left one), GD converges; however, when h is slightly larger than this bound, it blows up. Hence our restriction for h is very close to the stability limit.

heading: F OVERPARAMETRIZED OBJECTIVE: LARGE LEARNING RATE
text: In this section, we use column vector instead of row vector for sake of better understanding, i.e., our objective function is (µ− x>y)2/2.

heading: F.1 EIGENVALUES OF HESSIAN
text: Lemma F.1 (Matrix determinant lemma). Suppose A is an invertible n× n matrix and u, v ∈ Rn are column vectors. Then
det(A+ uv>) = (1 + v>A−1u) detA.
Theorem F.2. The eigenvalues of the Hessian of (µ − x>y)2/2 are ±(µ − x>y) repeated n − 1 times and 12 ( ‖x‖2 + ‖y‖2 ± √ (‖x‖2 + ‖y‖2)2 + 4(µ− x>y)2 − 8(µ− x>y)x>y ) . Especially, at x>y = µ, the eigenvalues are ‖x‖2 + ‖y‖2 and 0 repeated 2n− 1 times.
Proof. Consider the objective f(x, y) = (µ− x>y)2/2. Its Hessian is the following
H =
( yy> (x>y − µ)In + yx> (x>y − µ)In + xy> xx> )
=
( 0 (x>y − µ)In
(x>y − µ)In 0
) + ( y x )( y x )> .
Then to calculate the eigenvalues, we also need λI2n −H = (
λIn (µ− x>y)In (µ− x>y)In λIn ) − ( y x )( y x )> ∆ = B − ( y x )( y x )> .
By Lemma F.1, we have for invertible B
det ( B − ( y x )( y x )>) = ( 1− ( y x )> B−1 ( y x )) detB.
Since (µ− x>y)In and λIn commute, we have
detB = det((λ2 − (µ− x>y)2)In) = (λ2 − (µ− x>y)2)n.
By the formula for inversion of block matrix, we have for λ2 − (µ− x>y)2 6= 0,
B−1 =
( λ
λ2−(µ−x>y)2 In − µ−x>y λ2−(µ−x>y)2 In
− µ−x >y
λ2−(µ−x>y)2 In λ λ2−(µ−x>y)2 In
) .
Then combining all these and by the continuity of characteristic polynomial, we obtain the following expression for all λ
det(λI2n −H) = (λ2 − (µ− x>y)2)n−1(λ2 − λ(x>x+ y>y)− (µ− x>y)2 + 2(µ− x>y)x>y).
When x>y = µ, it becomes
det(λI2n −H) = λ2n−1(λ− (x>x+ y>y)).

heading: F.2 LOCAL NON-CONVEXITY NEAR THE GLOBAL MINIMUM
text: Consider the regionD = {δ > x>y−µ > 0}, a small neighbourhood of the global minima µ = x>y for some δ > 0. Consider two points [x, y], [x̄, ȳ] ∈ D with x = cx̄ and y = ȳ/c for c > 0 a constant.
Then f(x, y) + 〈∇f(x, y), [x̄− x, ȳ − y]〉 = (µ− x>y)2/2 + (µ− x>y) ( y>(x− x̄) + x>(y − ȳ) ) = (µ− x̄>ȳ)2/2 + (2− c− 1/c)(µ− x̄>ȳ)x̄>ȳ ≥ (µ− x̄>ȳ)2/2 = f(x̄, ȳ).
This contradicts the definition of convexity. Therefore the objective is not locally convex.
F.3 A REVIEW OF TRADITIONAL CONVERGENCE ANALYSIS OF GD UNDER L-SMOOTHNESS
For general function f , GD is defined as follows
xk+1 = xk − h∇f(xk). (9)
Theorem F.3. If f : RN → R is L−smooth, then with h < 2L , GD converges to a stationary point.
Proof. Let min f = f∗.
f(xk+1) ≤ f(xk) + 〈∇f(xk), xk+1 − xk〉+ L
2 ‖xk+1 − xk‖2
= f(xk)− h(1− L
2 h)‖∇f(xk)‖2.
Then N∑ k=1 ‖∇f(xk)‖2 ≤ 1 h(1− L2 h) (f(x0)− f(xN ))
≤ 1 h(1− L2 h) (f(x0)− f∗).
Therefore limk→∞ ‖∇f(xk)‖2 = 0, i.e., GD converges to a stationary point.

heading: G PROOF OF THEOREM 3.1, THEOREM 3.2, AND COROLLARY 3.3
text: In this section, we use column vectors for x and y instead of row vectors for better understanding. Then the objective function is 12 (µ− x >y)2.
The following Theorem G.1 and G.2 are the main theorems of the over-parametrized scalar case. Next, for sake of convenience, let u2k = ‖xk‖2 + ‖yk‖2. Theorem G.1 (Main 1). Let h = 4
u20+4cµ . Assume c ≥ 1 and u20 > 8µ. Then GD converges to a
point in {(x, y) : ‖x‖2 + ‖y‖2 ≤ 2h , x >y = µ} except for a Lebesgue measure-0 set of (x0, y0).
Proof. By Lemma G.7 and Lemma G.11, the theorem holds.
Theorem G.2 (Main 2). Let h = 48µ+4cµ = 1 (2+c)µ . Assume c ≥ 1 and u 2 0 ≤ 8µ. Then GD converges to a point in {(x, y) : ‖x‖2 + ‖y‖2 ≤ 2h , x >y = µ} except for a Lebesgue measure-0 set of (x0, y0).
Proof. For the measure-0 set, by Lemma G.5, the set of points converging to {‖x‖2 + ‖y‖2 > 2h} is measure-0; by the proof of Lemma G.11, the set of points converging to the origin is measure-0; also {u20 = ‖x‖ 2 + ‖y‖2 = 0} is a hyperplane and thus is measure-0. Hence the set of all the initial conditions not converging to {‖x‖2 + ‖y‖2 > 2h , x >y = µ} is measure-0.
Since u20 ≤ 8µ, if u20 > 2h , by Lemma G.7, it will decrease to u 2 k ≤ 2h for some k. Then by Lemma G.11, we have the convergence.
With the above two theorems, we can prove all the theorems and corollary in Section 3.
Proof of Theorem 3.1 and Theorem 3.2. From Theorem G.1, h = 4 u20+4cµ for c ≥ 1 which implies h ≤ 4
u20+4µ for all u20 > 8µ. Similarly, for Theorem G.2, h ≤ 13µ . When u 2 0 = 8µ, 4 u20+4µ = 4
8µ+4µ = 1 3µ ; when u 2 0 > 8µ, 4 u20+4µ < 13µ ; when u 2 0 < 8µ, 4 u20+4µ > 13µ . Also, all the limit points are in {‖x‖2 + ‖y‖2 ≤ 2h}. Then we can get Theorem 3.1 and 3.2 where the second inequality of Theorem 3.2 is because at the global minimum x>y = µ and also ‖x − y‖2 = ‖x‖2 + ‖y‖2 − 2x>y.
Proof of Corollary 3.3. For |x>0 y0 − µ| < δ, ‖x0‖2 + ‖y0‖2 > 8µ, we have
x>0 y0 < µ+ δ ⇒ ‖x0 − y0‖ 2 = ‖x0‖2 + ‖y0‖2 − 2x>0 y0 > 4 h − 4µ− 2(µ+ δ) = 4 h − 6µ− 2δ.
Also, from Theorem 3.2, ‖x− y‖2 ≤ 2h − 2µ and 2 h − 2µ < 4 h − 6µ− 2δ < 4 h − 8µ. We then obtain ‖x− y‖2 ≤ 12 ‖x0 − y0‖ 2 + 2µ.
We will divide the proof of the two main theorems into two phases: (1) when 2h < x 2 k + y 2 k < 4 h , we would like to show that GD escapes to a smaller ball x2k + y 2 k ≤ 2hexcept for a measure-0 set; (2) once GD enters x2k + y 2 k ≤ 2h , it will converge to the global minimum inside this region except for a measure-0 set.
G.1 PHASE 1: 2h < u 2 k < 4 h
We first deal with the situation where GD just converges in this region by showing that such points form a null set. This is stated in Theorem G.5.
Theorem G.3 and Corollary G.4 are the preliminary of proving Theorem G.5. Also Corollary G.4 is a direct result of Theorem G.3.
Theorem G.3. Let f : RN → RM and f ∈ C1. If the set of critical points of f is a null-set, i,e.,
L({x ∈ RN : ∇f(x) is not invertible}) = 0,
then L(f−1(B)) = 0 for any null-set B.
Proof. Let G = {|∇f | 6= 0} and G is an open set, where |∇f | denotes the determinant of ∇f . By implicit function theorem, we have G ∩ f−1(z) is a (N − 1)-submanifold in C1. Consider G∩ {f ∈ B} = ⋃ z∈B G∩ f−1(z), where L(G∩ f−1(z)) = 0 from the above discussion. Hence, if B is countable, L(G ∩ {f ∈ B}) = 0.
When B is uncountable, using co-area formula from geometric measure theory, we have for any bounded open ball Br,∫
G∩Br g|∇f |dL = ∫ R ∫ f−1(z)∩G∩Br g(x)dLN−1(x)dz = ∫ B ∫ f−1(z)∩G∩Br g(x)dLN−1(x)dz.
Let g be the indicator of ⋃ z∈B G ∩Br ∩ f−1(z). Then since B is a null-set, we have∫
B ∫ f−1(z)∩G∩Br g(x)dLN−1(x)dz = ∫ B ∫ f−1(z)∩G∩Br dLN−1(x)dz = 0
⇒ ∫ G∩Br g|∇f |dL = 0.
Then g|∇f | = 0 a.e. in G ∩Br. Also, |∇f | 6= 0 a.e.. Therefore, g = 0 a.e. in G ∩Br, i.e., L ( ⋃ z∈B G ∩Br ∩ f−1(z) ) = ∫ gdL = 0.
Since we can find a sequence of bounded open ball Br, s.t., RN = ⋃∞ r=1Br, and also L(Gc) = 0, we have
L(f−1(B)) ≤ L ( ⋃ z∈B G ∩ f−1(z) ) + L(Gc) = L ( ∞⋃ r=1 ⋃ z∈B G ∩Br ∩ f−1(z) ) + L(Gc) = 0.
Corollary G.4. Let ψ : R2d → R2d be the GD iteration map, i.e., ψ(x, y) = [x+ h(µ− x>y)y, y+ h(µ− x>y)x]>. Then if
L({det(Dψ) = 0}) = 0,
then L(ψ−1(B)) = 0 for any null-set B. Lemma G.5. Given h = 4
u20+4cµ and u2k > 2 h , GD will not converge to the points in {‖x‖ 2 +‖y‖2 > 2 h , x >y = µ}, except for a measure-0 set.
Proof. For d = 1,
det(Dψ) = 1− hx2 − hy2 − 3h2x2y2 + 4h2xyµ− h2µ2 = p1(x, y) · · · pm(x, y) = 0,
where pi(x, y) is irreducible polynomial and pi(x, y) = 0 is a co-dimensional-1 manifold. Then {(x, y) : det(Dψ) = 0} = ⋃m i=1{pi(x, y) = 0} are measure zero, i.e., L({det(Dψ) = 0}) = 0. Similarly for d > 1, we also have L({det(Dψ) = 0}) = 0.
From the GD iteration, we have
x>k+1yk+1 − µ = (x>k yk − µ) · (1− h(‖xk‖2 + ‖yk‖2)− h2x>k yk(µ− x>k yk)).
Then let B = {(x, y) : 1 − h(‖x‖2 + ‖y‖2) − h2x>y(µ − x>y) = 0}. For any (x, y) ∈ B, let [x+, y+]
> = ψ(x, y). Then x+y+ = µ. Similarly, L(B) = 0. By Corollary G.4 L(ψ−1(B)) = 0 and then L(ψ−n(B)) = L(ψ−1 ◦ · · · ◦ ψ−1(B)) = 0. Let ψ−0(B) = B and G = ⋃∞ i=0 ψ
−n(B). Hence
L(G) = L ( ∞⋃ i=0 ψ−n(B) ) ≤ ∞∑ i=0 L(ψ−n(B)) = 0.
Moreover, for any 0 < < µ, assume |x>k yk − µ| < . When u2k ≥ 2h + h(µ+ ),
|x>k+1yk+1 − µ| = |x>k yk − µ| · |1− hu2k − h2x>k yk(µ− x>k yk)| = |x>k yk − µ| · (−1 + hu2k + h2x>k yk(µ− x>k yk))
≥ |x>k yk − µ| · ( − 1 + h ( 2 h + h(µ+ ) ) − h2(µ+ ) ) > |x>k yk − µ|.
Hence, {‖x‖2 + ‖y‖2 > 2h , x >y = µ} is not the limit of this GD map ψ, except for the measure-0 set G.
Next we show that GD will be bounded inside {‖x‖2 + ‖y‖2 < 4h}.
Lemma G.6. Given h = 4 u20+4cµ , then for 0 ≤ k < min{k : u2k ≤ 2h}, we have u 2 k ≤ 4h − 3µ < 4 h for all k.
Proof. First u20 = 4 h − 4cµ ≤ 4 h − 4µ < 4 h . Then if x > 0 y0 > µ or x > 0 y0 < − hµu20 4−hu20 , from Lemma G.8, u21 < u 2 0 < 4 h . If − hµu20 4−hu20
≤ x>0 y0 < µ, from Lemma G.9 and its proof, u22 < u20 < 4h and u1 ≤ u20 + µ c ≤ u 2 0 + µ ≤ 4h − 3µ < 4 h . Therefore, iteratively we have u 2 k < 4 h .
Therefore, without loss of generality, we can just assume u2k ≤ u20 for a fixed kth iteration that we need to analyze because for every two step there exists an ith iteration such that u2i ≤ u20 and we can choose k = i to do the analysis.
Lemma G.8 and G.9 describe one-step or two-step decay of u2k, which lead to the primary result, Lemma G.7, in phase one. Moreover, the proof of Lemma G.7 contains a finer characterization of u2k making it possible to end phase one and enter phase two. Lemma G.7. Given h = 4
u20+4cµ for c ≥ 1, GD will enter {‖x‖2 + ‖y‖2 ≤ 2h} except for a
measure-0 set of initial conditions.
Proof. By Lemma G.6 and its discussion, assume without loss of generality 2h < u 2 k ≤ u20. From Lemma G.8 and G.9, the region where the decrease of u2k may be small is when xk, yk are close to x>k yk = µ or x > k yk = − hµu2k 4−hu2k .
Since
x>k+1yk+1 − µ = (x>k yk − µ)(1− hu2k − h2x>k yk(µ− x>k yk)),
consider sk = 1 − hu2k − h2x>k yk(µ − x>k yk). From the proof of Lemma G.9, we know when x>k yk = − hµu2k 4−hu2k
, sk < −1+2h2µ2 < 0. Therefore, there exists δ > 0, s.t., for all xk, yk ∈ {|x>k yk+ hµu2k
4−hu2k | < δ, 2h < u 2 k ≤ u20}, sk < 0. Then x>k+1yk+1 > µ. Hence, when −δ ≤ x>k yk + hµu2k 4−hu2k < 0, we will skip this step and consider the decrease of the next step with x>k+1yk+1 > µ. Also, there exist β1 = β1(δ) > 0, s.t., when x>k yk + hµu2k
4−hu2k < −δ,
u2k+1 − u2k = h(µ− x>k yk)((4− hu2k)x>k yk + hu2kµ) < −β1.
From the proof of Lemma G.9, when − hµu 2 k
4−hu2k ≤ x>k yk ≤ 0, for β2 = max{hµ2(4 − hu20)(1 −
2h2µ2)2, hµ2(8− h(2u20 + µ c ))} > 0,
u2k+2 − u2k < max{−h(µ− x>k yk)µ(4− hu2k+1)(1− 2h2µ2)2,−h(µ− x>k yk)2(8− h(u2k + u2k+1))} ≤ −β2.
Fix an small in 0 < < µ. When x>k yk ≥ µ + , from the proof of Lemma G.8, we have for β3 = 4hµ > 0,
u2k+1 − u2k ≤ 4hµ(µ− x>k yk) ≤ −4hµ = −β3.
When 0 < x>k yk ≤ µ − with > 0, from the proof of Lemma G.9, we have for β4 = h 2(8 − h(2u20 + µ c )) > 0,
u2k+2 − u2k ≤ −h(µ− x>k yk)2(8− h(u2k + u2k+1)) ≤ −β4.
When |x>k yk−µ| < , assume |x>k yk−µ| = k. In this case, we have sk < 0 meaning GD oscillates around x>y = µ. Also, if 0 < x>k yk < µ and u 2 k > 2 h , sk = 1 − hu 2 k − h2x>y(µ − x>y) < −1, i.e., if GD is in 0 < x>y < µ, then we have |x>k+1yk+1 − µ| > |x>k yk − µ|. Hence we only need to focus on the other side which is x>k yk > µ. From Lemma G.5, if u 2 k ≥ 2h + h(µ + ), |x>k+1yk+1−µ| > |x>k yk −µ|. Then within finite steps (note all these steps satisfy either one-step or two-step decrease of u2k; we ignore these steps only because the decrease maybe small), we will have either |x>k yk − µ| keeps increasing or u2K < 2h + Kh(µ+ K) (here we also consider x > KyK > µ). For the former one, the decrease of u2k for each step will be lower bounded away from 0; for the latter one, from the discussion above, u2K+1 − u2K ≤ −4hµ K ⇒ u2K+1 ≤ 2h .
From Lemma G.5, we know GD will not terminate in finite steps except for measure-0 set. From all the discussion above, we have that u2k decreases by a constant for either one-step or two-step except and hence GD will enter ‖x‖2 + ‖y‖2 ≤ 2h except for a measure-0 set.
Lemma G.8. Given h = 4 u20+4cµ and u2k < 4 h , when x > k yk > µ or x > k yk < − hµu2k 4−hu2k , we have u2k+1 − u2k < 0.

heading: Proof.
text: u2k+1 − u2k = h(µ− x>k yk)((4− hu2k)x>k yk + hu2kµ) = h(µ− x>k yk)(4x>k yk + hu2k(µ− x>k yk))
If x>k yk > µ, by u 2 k < 4 h ,
u2k+1 − u2k = h(µ− x>k yk)(4x>k yk + hu2k(µ− x>k yk)) ≤ 4hµ(µ− x>k yk) < 0.
If x>k yk < − hµu2k 4−hu2k ⇒ x>k yk < µ, then,
u2k+1 − u2k = h(µ− x>k yk)((4− hu2k)x>k yk + hu2kµ) < 0.
Lemma G.9. Given h = 4 u20+4cµ , c > max{ 12 , 2hµ}, c > 0, and u 2 k ≤ u20, when − hµu2k 4−hu2k
≤ x>k yk < µ, then u 2 k+2 − u2k < 0.
Proof. For every u2k, there exist a constant ck > 0, s.t. h = 4 u2k+4ckµ . Since 2h < u 2 k ≤ u20 < 4h , we have c ≤ ck < 12hµ . Since
x>k+1yk+1 − µ = (x>k yk − µ)(1− hu2k − h2x>k yk(µ− x>k yk)) = (x>k yk − µ)sk,
then sk = 1−hu2k−h2x>k yk(µ−x>k yk) is bounded by the value at x>k yk = − hµu2k 4−hu2k or x>k yk = µ, i.e.,
sk = 1− hu2k − h2x>k yk(µ− x>k yk) ≤ max { 1− hu2k + 4h3u2kµ 2
(4− hu2k)2 , 1− hu2k } Since u2k > 2 h , we have 1− hu 2 k < −1. For the other one, since u2k > 2h and ck < 1 2hµ ,
1− hu2k + 4h3u2kµ 2
(4− hu2k)2 = (1− 3c2k)u2k + 4c3kµ c2k(u 2 k + 4ckµ) < −1 + 2h2µ2,
where c > hµ+ √
2−h2µ2 2(1−h2µ2) . This is because either u 2 0 > 8µ ⇒ hµ < 13µ or h ≤ 1 3µ and then
hµ+ √
2−h2µ2 2(1−h2µ2) < 1 ≤ c. Also note here this is c not ck; this value −1 + 2h 2µ2 is achieved when
u2k = 2/h and ck = hµ+ √ 2−h2µ2 2(1−h2µ2) or 1 2hµ .
Also, when 0 ≤ x>k yk < µ,
sk = 1− hu2k − h2x>k yk(µ− x>k yk) ≤ 1− hu2k < −1.
We then prove 4− hu2k+1 > 0. When x>k yk < µ,
u2k+1 − u2k = h(µ− x>k yk)((4− hu2k)x>k yk + hu2kµ).
The maximum is achieved at x>k yk = 2−hu2k 4−hu2k µ, i.e., u2k+1 − u2k ≤ µ ck ≤ µc . Hence 4− hu 2 k+1 > 0 when c > 1/2.
Since x>k yk < µ and sk ≤ 0,
u2k+2 − u2k = u2k+2 − u2k+1 + u2k+1 − u2k = h(µ− x>k+1yk+1)((4− hu2k+1)x>k+1yk+1 + hu2k+1µ)
+ h(µ− x>k yk)((4− hu2k)x>k yk + hu2kµ) = h(µ− x>k yk)(((4− hu2k)x>k yk + hu2kµ) + sk(hu2k+1µ+ (4− hu2k+1)(sk(x>k yk − µ) + µ))) = h(µ− x>k yk) ( x>k yk(4− hu2k) + s2k(4− hu2k+1)(x>k yk − µ) + hu2kµ+ sk(hu2k+1µ+ (4− hu2k+1)µ) ) .
When 0 ≤ x>k yk < µ, s < −1, then
u2k+2 − u2k ≤ −h(µ− x>k yk)2(8− h(u2k + u2k+1)) < 0
When − hµu 2 k
4−hu2k ≤ x>k yk < 0, sk ≤ 0, then for c ≥ 2hµ,
u2k+2 − u2k < h(µ− x>k yk)µ(−4 + hu2k + 8h2µ2 − (4− hu2k+1)(1− 2h2µ2)2) < −h(µ− x>k yk)µ(4− hu2k+1)(1− 2h2µ2)2 < 0
where this bound is achieved by taking sk = −1 + 2h2µ2 and x>k yk = 0.
G.2 PHASE 2: u2k ≤ 2h
In this part, we will show the convergence of GD in Lemma G.11.
There are two convergence patterns: (i) transversal convergence, i.e., oscillating around the valley; (ii) unilateral convergence, i.e., converging from one side of the valley.
The key point of the proof of pattern (i) is to analyze the change of sk = 1−hu2k−h2x>k yk(µ−x>k yk). We know x>k+1yk+1 − µ = sk(x>k yk − µ), i.e., sk measures the change of the loss. If for some constant K > 0 we have |sk| < 1 ∀k ≥ K such that |x>k yk − µ| is guaranteed to decrease to 0, then the convergence follows. Moreover, we will need to analyze s2k and s2k+1 separately because in this case sk < 0 and x>k yk − µ changes sign at each step.
For pattern (ii), we will show that the trajectory of GD is bounded in a subset of this region |sk| < 1 that guarantees the decrease of the loss.
Before presenting our main result in phase 2, we first show a boundedness theorem when GD enters this phase.
Lemma G.10. Once GD enters {‖x‖2 + ‖y‖2 ≤ 2h}, it will stay bounded inside {‖x‖ 2 + ‖y‖2 ≤ 2 h + 2µ} and re-enter {‖x‖ 2 + ‖y‖2 ≤ 2h} within finite steps where the number of such steps do not depend on the number of iteration.
Proof. If at step K, u2K ≤ 2h , then from Lemma G.15, u 2 K+1 ≤ 2h + µ and it returns to phase 1. From the proof of Lemma G.9, u2K+2 ≤ u2K+1 + µ ≤ 2h + 2µ and we know either u 2 K+2 < u 2 K+1 or u2K+3 < u 2 K+1 and so on until it re-enters u 2 i ≤ 2h for some i ≥ K. Therefore, all the u 2 k ≤ 2h + 2µ once GD enters {‖x‖2 + ‖y‖2 ≤ 2h}.
Then by Lemma G.10, we can always pick a kth iteration such that u2k ≤ 2h . Also, when u 2 0 > 8µ, we have hµ < 13 . Together with the choice of h = 1 (2+c)µ when 0 < u 2 0 ≤ 8µ, we have hµ ≤ 13 . Then we obtain the following convergence of GD inside phase 2. Lemma G.11. Given hµ ≤ 13 , if GD enters {‖x‖ 2 + ‖y‖2 ≤ 2h}, it converge to x >y = µ inside this region except for a measure-0 set of initial conditions.
Proof. First, from Corollary G.4, we know the set of points converging to (0, 0) in finite steps is measure 0. For all the other points, we have the following discussion.
If sk = 0, then x>k+1yk+1 = µ, i.e., it converges.
If x>k yk < −µ and 1h + 2hµ 2 < u2k ≤ 2h , then by Lemma G.14,
u2k+1 − u2k = h(µ− x>k yk)((4− hu2k)x>k yk + hu2kµ) ≤ 2h(µ2 − (x>k yk)2) < 0.
Hence for δ1 > 0, when −µ− δ1 < x>k yk < −µ, we have x>k+1yk+1 > µ and we will skip this step and directly consider the (k+1)th iteration; when x>k yk ≤ −µ−δ1, u2k+1−u2k ≤ −2h((µ+δ1)2−µ2). Namely when x>k yk < −µ, either u2k+1 is some constant away from u2k or we can ignore the decrease in this step and directly look at the next one with x>k+1yk+1 > µ.
If x>k yk < −µ and u2k ≤ 1h + 2hµ 2, then
u2k+1 − u2k = h(µ− x>k yk)(4x>k yk + hu2k(µ− x>k yk)) ≤ h(µ− x>k yk)(4x>k yk + (1 + 2h2µ2)(µ− x>k yk)) ≤ −4hµ2(1− 2h2µ2) < 0.
If x>k yk > 3µ, then u 2 k+1 − u2k ≤ 2h(µ2 − x2ky2k) < −16hµ2. Also, when u2k ≤ 6µ, x>k yk ≤ 3µ.
Next, for the rest of the region, we first consider 1h + 2hµ 2 < u2k ≤ 2h (by Lemma G.10, we are always able to find u2k ≤ 2h ). In this region, by Lemma G.18, we know sk < 0, i.e., (x > k yk − µ)(x>k+1yk+1 − µ) < 0. We divide the whole region into two parts: x>k yk > µ and x>k yk < µ. Therefore all the xk+2i should be on the same side for i such that 1h + 2hµ 2 < u2k+2i ≤ 2h .
If −µ ≤ x>k yk < µ and sk ≤ −1, by Lemma G.16, we have u2k+2 − u2k ≤ −4h(1 − hµ)(µ − x>k yk)
2 < 0. Hence when GD does not converge, i.e., {(xk+2n, yk+2n)}∞n=0 does not converge , u2k will keep decreasing with u2k+2i ≤ 2h for allK+2iwith i ≥ 0 such that sk+2i ≤ −1. Moreover, there existsN > 0 s.t. sk+2N > −1 because sk = 1−hu2k−h2x>k yk(µ−x>k yk) and if xk+2i, yk+2i is the first iteration to leave this region, then it has to satisfy x>k+2iyk+2i < −µ which implies sk+2i > −1, and then it follows from Lemma G.13 and G.12 just as the following two paragraphs of discussion; if it never leaves this region and GD is not converging, i.e., |µ− x>k yk| has a lower bound, then u2k+2i will keep decreasing until sk+2N > −1 (from the proof of Lemma G.12, if µ− x>k yk is very small, then sk+1 < sk meaning both sides are blowing up and therefore |µ− x>k yk| has a lower bound if sk ≤ −1).
If −µ ≤ x>k yk < µ and sk > −1, from Lemma G.12 and Lemma G.13, similar to the previous discussion, there exists N1 > 0, s.t., for all k > N1, sk > min{−1 + c0,−1 + c1(µ− x>k yk)2} for some c0, c1 > 0, namely |µ− x>k yk| strictly decreases. Then it will either converge in 1h + 2hµ 2 < ‖x‖2 + ‖y‖2 ≤ 2h or enter ‖x‖ 2 + ‖y‖2 ≤ 1h + 2hµ 2.
If µ < x>k yk ≤ 3µ, from the definition of sk, we have sk > −1. Also by Lemma G.13, we know sm > min{−1 + c0,−1 + c1(µ− x>k yk)2} for all m > k. Hence similarly, it will either converge or enter u2k ≤ 1h + 2hµ 2.
Then we consider u2k ≤ 1h + 2hµ 2. From Lemma G.17 and Lemma G.20, we know |sk| < 1 and will be away from 1 by a constant in this area. We will show that GD stays bounded in the converging domain. If sk ≤ 0, it follows from the previous discussion and the proof of Lemma G.17 and Lemma G.20. If sk > 0, when the trajectory is in {|x>y| > µ}, from Lemma G.14, we have u2k+1 < u 2 k; when it is in {|x>y| < µ}, then |xk+1 − yk+1| < |xk − yk| from Lemma G.19. Hence the trajectory will stay inside the monotone decreasing region. Next, since |µ− x>k yk| monotonically decreases and is lower bounded by 0, we have that |µ− x>k yk| converges. If |µ− x>k yk| converges to C > 0, then min{−1 + c0,−1 + c1(µ− x>k yk)2} < sk < 1− c2 for some c0, c1, c2 > 0 meaning |µ− x>k yk| does not converges to C. Contradiction. Therefore, |µ− x>k yk| converges to 0.
Lemma G.12. When −µ ≤ x>k yk < µ and 1h + 2hµ 2 < u2k ≤ 2h + 2µ, if sk > −1, then sk+2 > min{− 34 − h2µ2 4 ,− 5h2µ2 2 ,−1 + h 2(µ− x>k yk)2} > −1, for 92h 2µ2 < 1.
Proof. From previous statement, sk < 0 when u2k > 1 h + 2hµ 2. For sk = 1− hu2k − h2x>k yk(µ− x>k yk), let = µ− x>k yk. Then x>k+1yk+1 − µ = −sk .
sk+1 = 1− hu2k+1 − h2x>k+1yk+1(µ− x>k+1yk+1) = sk − h2 µ(3 + sk) + h2 2(3 + s2k − hu2k).
Hence if is very small, we have sk+1 < sk. Also, xk+2yk+2 − µ = sk+1(x>k+1yk+1 − µ) = −sk+1sk . Moreover, u2k+1 ≤ 2h + 2µ by Lemma G.6.
sk+2 = 1− hu2k+2 − h2x>k+2yk+2(µ− x>k+2yk+2) = sk − h2 µ(3 + 4sk + sksk+1) + h2 2(3 + s2k − hu2k + s2k(3 + s2k+1 − hu2k+1))) ≥ sk − h2 µ(3 + 4sk + sksk+1) + h2 2(1 + (1 + s2k+1)s2k) sk+1=µ/(2sk )
≥ sk − h2µ2
4 − h2 (3 + 4sk) + h2 2(1 + s2k)
If 3 + 4sk ≥ 0, then
sk+2 ≥ sk − h2µ2
4 − h2 (3 + 4sk) + h2 2(1 + s2k)
=(3+4sk)µ/(2(1+s 2 k)) ≥ sk − h2µ2 4 − h 2µ2(3 + 4sk) 2
4(1 + s2k)
sk=0 or sk=−3/4 ≥ min{−3 4 − h 2µ2 4 ,−5h 2µ2 2 } > −1.
If 3 + 4sk < 0, i.e., −1 < sk < − 34 , then
sk+2 = sk − h2 µ(3 + 4sk + sksk+1) + h2 2(3 + s2k − hu2k + s2k(3 + s2k+1 − hu2k+1)) = sk − h2 µ(3 + 4sk + sk(sk − h2 µ(3 + sk) + h2 2(3 + s2 − hu2k)))
+ h2 2(3 + s2k − hu2k + s2k(3 + s2k+1 − hu2k+1)) = sk − h2µ (sk(3 + sk) + 3 + sk)− h4µ 3sk(3 + s2k − hu2k)
+ h2 2(3 + s2k − hu2k + s2k(3 + s2k+1 − hu2k+1) + h2µ2sk(3 + sk)).
When −1 < sk < 1−3h 2 µ h2 µ , sk − h 2µ (sk(3 + sk) + 3 + sk) > −1. Since ≤ 2µ, we have 1−3h2 µ h2 µ ≥ 1−6h2µ2 2h2µ2 ≥ − 3 4 for 9 2h 2µ2 < 1. Also,
− h4µ 3sk(3 + s2k − hu2k) > 3
4 h4µ 3 > 0,
h2 2(3 + s2k − hu2k + s2k(3 + s2k+1 − hu2k+1) + h2µ2sk(3 + sk)) > h2 2 > 0. Hence
sk+2 > −1− h4µ 3sk(3 + s2k − hu2k) + h2 2(3 + s2k − hu2k + s2k(3 + s2k+1 − hu2k+1) + h2µ2sk(3 + sk)) > −1 + h2 2
Lemma G.13. When x>k yk > µ and 1 h + 2hµ 2 < u2k ≤ 2h + 2µ, if sk > −1, then sk+1 > sk and sk+2 > min{−1 + h2(µ− x>k yk)2( 12 − 3 2h 2µ2),− 12 − 3 2h
4(µ− x>k yk)2µ2} > −1, for 8h2µ2 < 1. Moreover, if u2k+1 ≤ 2h , sk ≤ −1, and sk+1 > −1, then sk+2 > sk + h 2(µ− x>k yk)2(1− 8h2µ2).
Proof. From previous statement, sk < 0 when u2k > 1 h + 2hµ 2. Without loss of generality, we can just assume u2k < 3 h (otherwise u 2 k+2 < 3 h and we can use this as our starting point). For x > k yk > µ, sk = 1− hu2k − h2x>k yk(µ− x>k yk). Let x>k yk − µ = . Then
sk+1 = 1− hu2k+1 − h2x>k+1yk+1(µ− x>k+1yk+1) = sk + h2 µ(3 + sk) + h2 2(3 + s2k − hu2k) > sk,
where sk > −3 in this region. Also, x>k+2yk+2 − µ = sk+1sk .
Then for sk ≤ −1, u2k+1 ≤ 2h , and sk+1 > −1,
sk+2 = sk + h 2 µ(3 + 4sk + sksk+1) + h 2 2(3 + s2k − hu2k + s2k(3 + s2k+1 − hu2k+1)) = sk + h 2 µ(3 + 4sk+1 + s 2 k+1) + h 2 2 [ (3 + s2k − hu2k)(1− 4h2 µ− h2 µsk+1)
+ s2k(2− hu2k+1) + s2k(1 + s2k+1)− h2µ2(4 + sk+1)(3 + sk) ]
> sk + h 2 2(1− 8h2µ2)
where the inequality is achieved by sk+1 > −1, u2k ≤ 4h , u 2 k+1 ≤ 2h , 8h 2µ2 < 1, ≤ 2µ.
For sk > −1,
sk+2 = sk + h 2 µ(3 + 4sk + s 2 k) + h 2 2 [ (3 + s2k − hu2k)(1 + skh2 µ)
+ s2k(3 + s 2 k+1 − hu2k+1) + h2µ2sk(3 + sk) ] > sk + h 2 µ(3 + 4sk + s 2 k) + h 2 2 [ s2k(1 + 2skh 2µ2) + s2k + h 2µ2sk(3 + sk)
] = sk + h 2 µ(3 + 4sk + s 2 k) + h 2 2 [ 2s2k + h 2µ2sk(3 + 2s 2 k + sk)
] When sk > −1, we have sk + h2 µ(3 + 4sk + s2k) > −1 and h2 µ(3 + 4sk + s2k) > 0. Since 8h2µ2 < 1, 2s2k + h 2µ2sk(3 + 2s 2 k + sk) ≥ 12 − 3 2h 2µ2 > 0 for −1 < sk ≤ − 12 . For − 1 2 < sk ≤ 0, 2s2k + h2µ2sk(3 + 2s2k + sk) ≥ h2µ2sk(3 + 2s2k + sk) ≥ − 32h 2µ2 > − 12 . Hence sk+2 > min{−1 + h2 2( 12 − 3 2h 2µ2),− 12 − 3 2h 4 2µ2} > −1.
Lemma G.14. Given u2k ≤ 2h , when |x > k yk| > µ, u2k+1 − u2k < 0.
Proof.
u2k+1 − u2k = h(µ− x>k yk)(4x>k yk + hu2k(µ− x>k yk)) = h(µ− x>k yk)4x>k yk + h2u2k(µ− x>k yk)2
≤ h(µ− x>k yk)4x>k yk + 2h(µ− x>k yk)2 = 2h(µ2 − (x>k yk)2) < 0
Lemma G.15. Given u2k ≤ 2h , u 2 k+1 − u2k ≤ µ.

heading: Proof.
text: u2k+1 − u2k = h(µ− x>k yk)(4x>k yk + hu2k(µ− x>k yk)) ≤ 4hµ2
4− hu2k ,
where the maximum is achieved at x>k yk = 2−hu2k 4−hu2k µ. When h = 4 u20+4cµ , from the proof of Lemma G.9, uk+1 − u2k ≤ µ. When h = 1(2+c)µ and in this case u 2 k ≤ 2h = (1 + c 2µ), uk+1 − u 2 k ≤
8 7(2+c)µ < µ.
Lemma G.16. When−µ ≤ x>k yk < µ, if sk = 1−hu2k−h2x>k yk(µ−x>k yk) ≤ −1, u2k+2−u2k ≤ −h(4− 4hµ)(µ− x>k yk)2 < 0.

heading: Proof.
text: u2k+2 − u2k = h(µ− x>k yk)(((4− hu2k)x>k yk + hu2kµ) + s(hu2k+1µ+ (4− hu2k+1)(s(x>k yk − µ) + µ))) = h(µ− x>k yk) ( −(4− hu2k + s2(4− hu2k+1))(µ− x>k yk) + 4(1 + s)µ ) ≤ −h(8− h(u2k + u2k+1))(µ− x>k yk)2
≤ −h(4− 4hµ)(µ− x>k yk)2 < 0
where the first inequality is because when sk = −1 it has the maximum value; the second is from Lemma G.10.
Lemma G.17. When u2k ≤ 2h , then sk < 1 and will be away from 1 by a constant.

heading: Proof.
text: sk = 1− hu2k − h2x>k yk(µ− x>k yk) ≤ 1− hu2k + h2 u2k 2
( µ+
u2k 2 ) where the equality is achieved when x>k yk = − u2k 2 .
Since 0 < u2k ≤ 2h ,
sk ≤ 1− hu2k + h2 u2k 2
( µ+
u2k 2
) ≤ max{1, hµ} < 1.
Actually this sk will be away from 1 by a constant. This is because when u2k are close to 0, sk will be close to 1; however, u2k will be increasing because of the decrease of |µ − x>k yk| and thus sk will decrease accordingly.
Lemma G.18. When u2k > 1 h + 2hµ 2 and −µ ≤ x>k yk ≤ 3µ, (x>k+1yk+1 − µ)(x>k yk − µ) < 0, meaning it oscillates around x>y = µ.

heading: Proof.
text: 1− hu2k − h2x>k yk(µ− x>k yk) < 1− h( 1
h + 2hµ2) + 2h2µ2 = 0
Lemma G.19. When −µ ≤ x>k yk < µ and hµ ≤ 13 , we have |xk+1 − yk+1| < |xk − yk|.
Proof.
|xk+1 − yk+1| = |xk − yk| |1− h(µ− x>k yk)|.
Since hµ < 1, then
−1 < 1− 2hµ ≤ 1− h(µ− x>k yk) < 1.
Lemma G.20. When u2k ≤ 1h + 2hµ 2 and hµ ≤ 13 , sk ≥ − 1 4 .
Proof. By hµ ≤ 13 ,
1− hu2k − h2x>k yk(µ− x>k yk) ≥ −2h2µ2 − h2µ2 4 ≥ −9h 2µ2 4 ≥ −1 4 .

heading: H PROOF OF THEOREM 4.1, THEOREM 4.2, AND THEOREM 4.3
text: The proof in this section mainly follows the strategy of the scalar case. More precisely, all the expressions, e.g., u2k+1 and x > k+1yk+1 − µ, in this rank-1 approximation case contains two parts: the one that can be handled using the technique in scalar case, and the other one measuring the extent of alignment between xk and yk. Again, let u2k = ‖xk‖2 + ‖yk‖2.
Theorem H.1 (Main 1). Let h = 4 u20+4cµ
. Assume c ≥ √ 7 and u20 > 4 √ 7µ. Then GD converges to
a point in {(x, y) : ‖x‖2 + ‖y‖2 ≤ 2h , x >y = µ} except for a Lebesgue measure-0 set of (x0, y0).
Proof. The proof follows from Lemma H.3 and Lemma H.8.
Theorem H.2 (Main 2). Let h = 4 4 √ 7µ+4cµ = 1 ( √ 7+c)µ . Assume c ≥ √ 7 and u20 ≤ 4 √ 7µ. Then GD converges to a point in {(x, y) : ‖x‖2 + ‖y‖2 ≤ 2h , x >y = µ} except for a Lebesgue measure-0 set of (x0, y0).
Proof. Since c ≥ √ 7, in this case we have u20 ≤ 4 √ 7µ ≤ 2( √
7 + c)µ = 2h . Therefore, the convergence follows from Lemma H.8.
Hence we can show the proof of the theorems in Section 4.
Proof of Theorem 4.1. By Lemma H.3 and H.7, the proof follows from Lemma H.9.
Proof of Theorem 4.2 and 4.3. Similar to scalar case, h = 4 u20+4cµ
for c ≥ √
7 implies h ≤ 4 u20+4 √ 7µ
and h = 1 ( √ 7+c)µ implies h ≤ 1 2 √ 7 . Therefore, the convergence and balancing of GD follow from Theorem H.1 and H.2. The second inequality of Theorem 4.3 follows from x>y = µ at the global minimum because of diagonal and non-negative A and the best rank-1 approximation in SVD.

heading: H.1 PHASE 1
text: The main theorem in phase 1 is the following.
Lemma H.3. Given h = 4 u20+4cµ
for c ≥ √
7, GD will enter {‖x‖2 + ‖y‖2 ≤ 2h} except for a measure-0 set of initial conditions.
Proof. Since the expressions, u2k+2−u2k and u2k+2−u2k, of the rank-1 version have an additional term upper bounded by c0((x>k yk)
2 − x>k xky>k yk) for some constant c0 > 0 which is non-positive, we have u2k decreases by a constant for either one- or two-step if (x > k yk)
2−x>k xky>k yk < −c1 for some constant c1 > 0 when 2h < u 2 k < 4 h . Therefore, we can pick a K > 0 s.t. (x > k yk)
2 − x>k xky>k yk < −c1 for k ≤ K. If GD enters {‖x‖2 + ‖y‖2 ≤ 2h} within K steps, then we are done. Otherwise, we can just assume that for all the k > K, −c1 ≤ (x>k yk)2 − x>k xky>k yk ≤ 0. Then again, like the scalar case, we only need to consider when xk, yk are close to x>k yk = µ and x > k yk = − hµu2k 4−hu2k . For the region near x>k yk = − hµu2k
4−hu2k , we can always choose a c1 s.t. there exists a δ > 0, when
−δ ≤ x>k yk + hµu2k 4−hu2k < 0, we have both x>k+1yk+1 > µ and sk < 0. Then all the discussion just follows the scalar case. For the region near x>k yk = µ, i.e., |x>k yk − µ| < for some 0 < < µ, consider
x>k+1yk+1 − µ = (x>k yk − µ)sk + h2(2µ− x>k yk)((x>k yk)2 − x>k xky>k yk).
Similar to scalar decomposition, we only need to consider the case when sk ≤ −1 because if sk > −1 which can only happen when x>k yk > µ, we have u 2 k ≤ 2h + kh(µ + ) and then u 2 k+1 ≤ 2h . If sk ≥ −1 for all the following iterations, together with the bounded second term in x>k+1yk+1 − µ, i.e., −2h2µc1 ≤ h2(2µ− x>k yk)((x>k yk)2 − x>k xky>k yk) ≤ 0, we have that there exists K > 0 s.t. |x>KyK − µ| ≥ 1 for some 1 > 0 unless it leaves this region |x>k yk − µ| < . Therefore, by the same discussion in scalar case, u2k will decrease either in one step or two steps by at least a constant.
Moreover, by Lemma H.4, GD will not terminate inside { 2h < ‖x‖ 2 + ‖y‖2 < 4h} in finite step. Thus, GD is guaranteed to enter {‖x‖2 + ‖y‖2 ≤ 2h} in finite steps.
Lemma H.4. GD will not converge to any fixed points in {‖x‖2 +‖y‖2 > 2h} except for a measure-0 set.
Proof. Similar to the discussion in scalar case, the set of points converging to x>y = µ in finite step is measure-0.
Also, since µx>y = x>xy>y = µ2, the fixed points of GD map requires x>xy>y = (x>y)2. This can be seen in the following expression
x>k+1yk+1 − µ = (x>k yk − µ)(1− hu2k − h2x>k yk(µ− x>k yk)) + h2(2µ− x>k yk)((x>k yk)2 − x>k xky>k yk),
where the convergence of x>k yk − µ requires the convergence of x>k xky>k yk − (x>k yk)2.
Hence if x>k xky > k yk − (x>k yk)2 does not converge, then GD will not converge. If x>k xky>k yk − (x>k yk) 2 converges, then c ≥ x>k xky>k yk − (x>k yk)2 ≥ 0 for all the rest of the iterations with some c > 0, and we further assume |x>k yk − µ| < for any 0 < < µ. As is shown in scalar case, when u2k ≥ 2h + h(µ+ ), we have sk = 1− hu 2 k − h2x>k yk(µ− x>k yk) ≥ −1. Then
x>k+1yk+1 − µ = (x>k yk − µ)sk + h2(2µ− x>k yk)((x>k yk)2 − x>k xky>k yk)
will never converge (because sk ≤ −1 and the second term is bounded).
Lemma H.5. When x>k yk > µ or x>k yk < − hµu2k
4−hu2k , we have u2k+1 − u2k < 0.

heading: Proof.
text: u2k+1 − u2k = h(µ− x>k yk)((4− hu2k)x>k yk + hu2kµ) + h(4− hu2k)((x>k yk)2 − x>k xky>k yk) ≤ h(µ− x>k yk)((4− hu2k)x>k yk + hu2kµ).
Hence the proof follows 1D case.
Lemma H.6. When − hµu 2 k
4−hu2k ≤ x>k yk < µ, we have u2k+2 − u2k < 0.
Proof. Let sk = 1− hu2k − h2x>k yk(µ− x>k yk), the same as 1D.
x>k+1yk+1 − µ = (x>k yk − µ)(1− hu2k − h2x>k yk(µ− x>k yk)) + h2(2µ− x>k yk)((x>k yk)2 − x>k xky>k yk) ∆ = (x>k yk − µ)sk +M.
u2k+2 − u2k = u2k+1 − u2k + u2k+2 − u2k+1 = h(µ− x>k yk)((4− hu2k)x>k yk + hu2kµ) + h(4− hu2k)((x>k yk)2 − x>k xky>k yk)
+ h(µ− x>k+1yk+1)((4− hu2k+1)x>k+1yk+1 + hu2k+1µ) + h(4− hu2k+1)((x>k+1yk+1)2 − x>k+1xk+1y>k+1yk+1)
= [ h(µ− x>k yk)(((4− hu2k)x>k yk + hu2kµ) + sk(hu2k+1µ+ (4− hu2k+1)(sk(x>k yk − µ) + µ))) ] + [ h(4− hu2k)((x>k yk)2 − x>k xky>k yk) + hM2(−4 + hu2k+1)
+ hM(−hu2k+1µ+ (−4 + hu2k+1)µ+ 2(4− hu2k+1)(µ− x>k yk)sk) ]
∆ = I + II.
I is the same as 1D u2k+2 − u2k and hence (1) < 0. For II,
II = hM2(−4 + hu2k+1) + h((x>k yk)2 − x>k xky>k yk) × (4− hu2 + h2(2µ− x>k yk)(−hu2k+1µ+ (−4 + hu2k+1)µ+ 2(4− hu2k+1)(µ− x>k yk)sk)︸ ︷︷ ︸
III
).
From 1D results, sk < 0. Then, III achieves its lower bound at x>k yk = µ or x > k yk = − hµu2k 4−hu2k . For x>k yk = µ,
III = 4− hu2k − 4h2µ2 ≥ 4hµ(c− hµ) > 0,
where the inequality is from u2k ≤ u20 = 4h − 4cµ and c > hµ.
For x>k yk = − hµu2k
4−hu2k ,
III = 4− hu2k − 4h2(8− hu2k)(4− hu2k − 2sk(4− hu2k+1))µ2
(4− hu2k)2
≥ 4− hu2k − 4h2(8− hu2k)(4− hu2k − 2sk(4− hu2k))µ2
(4− hu2k)2
= 4− hu2k − 4h2(8− hu2k)(1− 2sk)µ2 4− hu2k u2k≤u 2 0
≥ 4− hu2k − 4h2(4 + 4chµ)(1− 2sk)µ2
4chµ
u2k≤u 2 0
≥ 4hµ c (−1 + c2 + 2sk + ch(−1 + 2sk)µ)
∗ ≥ c2(1 + 8h2µ2) + c(hµ− h 3µ3 2 )− 7− h 2µ2
2 ∗∗ ≥ 0,
where ∗ follows from sk = 1− hu2k − h2x>k yk(µ− x>k yk)
x>k yk=µ/2
≥ 1− hu2k − h2µ2
4 u2k≤u 2 0
≥ −3 + 4chµ− h 2µ2
4 ,
and ∗∗ follows from c ≥ √ 7.
From 1D discussion, −4 + hu2k+1 < 0. Also, (x>k yk)2 − x>k xky>k yk < 0. Hence we have (3) ≥ 0.
Overall, (2) ≤ 0 and therefore the whole discussion of (1) is the same as 1D.

heading: H.2 PHASE 2
text: The proof in phase 2 is partly different from the scalar case due to the alignment. Before presenting the main convergence lemma, we first state the boundedness of GD in the following lemma. Lemma H.7. Once GD enters {‖x‖2 + ‖y‖2 ≤ 2h}, it will stay bounded inside {‖x‖
2 + ‖y‖2 ≤ 2 h + 2hµ 2(1 + 11−h2µ2 )} and re-enter {‖x‖ 2 + ‖y‖2 ≤ 2h} within finite steps where the number of such steps do not depend on the number of iteration.
Proof. Assume u2k ≤ 2h . Then u2k+1 − u2k = h(µ− x>k yk)((4− hu2k)x>k yk + hu2kµ) + h(4− hu2k)((x>k yk)2 − x>k xky>k yk) ≤ 2h(µ2 − x>k xky>k yk) ≤ 2hµ2, where the first inequality follows from u2k ≤ 2h . If further u 2 k+1 ≤ 2h + 2hµ
2, then by x>k+1yk+1 ≥ − √ x>k+1xk+1y > k+1yk+1 and −1 + h2µ2 < 0, we have
u2k+2 − u2k+1 = h(µ− x>k+1yk+1)((4− hu2k+1)x>k+1yk+1 + hu2k+1µ) + h(4− hu2k+1)((x>k+1yk+1)2 − x>k+1xk+1y>k+1yk+1)
≤ 2h(x>k+1xk+1y>k+1yk+1(−1 + h2µ2) + µ2(1− 2x>k+1yk+1h2µ+ h2µ2)) ≤ 2h(x>k+1xk+1y>k+1yk+1(−1 + h2µ2) + µ2(1 + 2 √ x>k+1xk+1y > k+1yk+1h 2µ+ h2µ2))
≤ 2hµ 2
1− h2µ2 .
The first inequality is from u2k+1 ≤ 2h + 2hµ 2; the second inequality follows from the two conditions mentioned above; the third inequality follows from choosing √ x>k+1xk+1y > k+1yk+1 = h2µ3 1−h2µ2 .
From previous proof, we know when u2k > 2 h , by Lemma H.3, we have GD will enter {‖x‖ 2 + ‖y‖2 ≤ 2h} in finite steps with the decrease of u 2 k in either one step or two steps. Then, u 2 i ≤ 2 h + 2hµ 2(1 + 11−h2µ2 ) for all i ≥ k.
Let Uk = x>k yk, Vk = x > k xk, Wk = y > k yk. Then we can rewrite the iteration for the three variables
Uk+1 = h(1− hVk)µWk + hµVk(1− hWk) + Uk(1 + h2VkWk − h(Vk +Wk) + h2µ2) Vk+1 = 2hµUk + Vk − 2hVkWk + h2(µ2Wk − 2µUkWk + VkW 2k ) Wk+1 = 2hµUk +Wk − 2hVkWk + h2(µ2Vk − 2µUkVk + V 2kWk) .
(10)
When h = 4 u20+4cµ or h = 1 2 √ 7µ , we have hµ ≤ 1 2 √ 7 . The following is the main lemma in phase 2.
Lemma H.8. Given hµ ≤ 1 2 √ 7 , if GD enters {‖x‖2 +‖y‖2 ≤ 2h}, it converges to x >y = ‖x‖ ‖y‖ = µ, i.e., the global minimum, inside this region except for a measure-0 set of initial conditions.
Proof. This proof follows the same method as the scalar decomposition. We need to prove that all the conclusions in scalar case also holds in this rank one decomposition.
Consider u2k+1 − u2k when x>k yk < −µ or x>k yk > µ,
u2k+1 − u2k = h(µ− x>k yk)((4− hu2k)x>k yk + hu2kµ) + h(4− hu2k)((x>k yk)2 − x>k xky>k yk) = h(µ− x>k yk)((4− hu2k)x>k yk + hu2kµ) + h(4− hu2k)(U2k − VkWk) ≤ h(µ− x>k yk)((4− hu2k)x>k yk + hu2kµ),
where the last expression is the same as scalar case and therefore under the conditions such that u2k+1 − u2k < 0 in scalar case, we also have u2k+1 − u2k < 0 in rank-one approximation.
As in scalar case, we also denote sk = 1− hu2k − h2x>k yk(µ− x>k yk). First consider 1h + 2hµ 2 < u2k ≤ 2h . If −µ ≤ x > k yk < µ, u 2 k+1 ≥ u2k, and sk ≤ −1, consider u2k+2 − u2k. From Lemma H.6, the extra terms of u2k+2 − u2k compared to scalar case is
hM2(−4 + hu2k+1) + h((x>k yk)2 − x>k xky>k yk) × (4− hu2 + h2(2µ− x>k yk)(−hu2k+1µ+ (−4 + hu2k+1)µ+ 2(4− hu2k+1)(µ− x>k yk)sk)︸ ︷︷ ︸
III
)
where M = h2(2µ − x>k yk)((x>k yk)2 − x>k xky>k yk). Then we will prove that III is positive and then all the extra terms will be less than 0. By rewriting III in terms of x>k yk, we get
III = 2(x>k yk) 2h2sk(4− hu2k+1) + x>k yk(h3u2k+1µ− 6h2sk(4− hu2k+1)µ− h2(−4 + hu2k+1)µ)
− 2h3u2k+1µ2 + 4h2sk(4− hu2k+1)µ2 + 2h2(−4 + hu2k+1)µ2 + 4− hu2k ≥ 4− hu2k + 12h2(−1 + 4sk)µ2 − 12h3sku2k+1µ2
≥ 4− hu2k + 12h2(−1 + 4sk)µ2 − 12h3sku2kµ2
= 4− hu2k + 12h2(−1 + sk(4− hu2k))µ2 > 0,
where the first inequality is because it obtains the minimum at x>k yk = −µ; the second inequality follows from u2k+1 ≥ u2k; the last inequality is from the range of u2k and sk. Therefore the extra terms of u2k+2 − u2k is always negative and then we have the same conclusion as scalar case, i.e., u2k+2 − u2k < −C(µ− x>k yk)2 for some constant C > 0.
Next, consider sk+1 when −µ ≤ x>k yk ≤ 3µ and 1h + 2hµ 2 < u2k ≤ 2h which implies −1−
h2µ2
4 ≤ sk < 0. The extra terms of sk+1 compared to scalar case is
h2(U2k − VkWk) ( − (4− hu2k) + (−x>k yk + 2µ)(2h2sk(x>k yk − µ) + h2µ) + h4x(−x>k yk + 2µ)2 ) ︸ ︷︷ ︸
IV
,
where U2k − VkWk ≤ 0. Then consider the upper bound of IV
IV = −4− 2(x>k yk)2h2sk + x>k yk(−h2µ+ 6h2sµ) + hu2k + 2h2µ2 − 4h2skµ2
≤ −4 + hu2k + (−h2 − 4h2sk)µ2 < 0,
where the first inequality follows from x>k yk ≤ 3µ. Therefore the extra terms is non-negative, namely all the lower bound for sk+1 in the scalar case also holds for rank-one approximation.
Also we need to consider sk+2 when sk > −1, −µ ≤ x>k yk ≤ 3µ, and 1h + 2hµ 2 < u2k ≤ 2h . For scalar case, from Lemma G.13 with = x>k+1yk+1 − µ, we have
sk+2 = sk+1 + (h 2(3 + s2k+1)− h3u2k) 2 + h2 (−3µ− sk+1µ)
= 3h2 2 + h2s2k+1 2 − h3u2k 2 − 3h2 µ+ sk+1(1− h2 µ).
Since the minimum point w.r.t. sk+1 is −(1−h2 µ) 2h2 2 < −1− h2µ2
4 , we have if sk+1 is larger, then sk+2 is also larger. Given the same value for sk for both scalar and rank-one cases, sk+1 is larger from the above discussion and thus sk+2 is also larger, namely, the lower bound for sk+2 also holds in rank-one cases.
Next consider u2k ≤ 1h + 2hµ 2 where |sk| < 1 and we need the iteration to have certain restriction such that it will be bounded in a region. If sk < 0, then it follows the same from the above discussion. If sk > 0 and x>k yk > µ or x > k yk < −µ, then by u2k+1 − u2k < 0, it will stay in {x>k yk > µ, u2k ≤ 1h + 2hµ 2} until convergence. If sk > 0 and −µ ≤ x>k yk < µ, similar to scalar case, we consider |x>k xk − y>k yk| = |Vk −Wk| and
|Vk+1 −Wk+1| = |Vk −Wk| · |1− h2VkWk + 2h2Ukµ− h2µ2|,
where Uk = x>k yk and U 2 k ≤ VkWk. Then
1− h2VkWk + 2h2Ukµ− h2µ2 ≤ 1− h2U2k + 2h2Ukµ− h2µ2 < 1,
1− h2VkWk + 2h2Ukµ− h2µ2 ≥ 1− h2 u4k 4 − 2h2µ− h2µ2 ≥ 1− 1 4 (1 + 2h2µ2)2 − 3h2µ2 > 0.
Therefore if −µ ≤ x>k yk < µ and u2k ≤ 1h + 2hµ 2, we have |Vk+1 − Wk+1| < |Vk − Wk|, meaning together with 0 ≤ µ− x>k+1yk+1 < µ− x>k yk, xk+1 and yk+1 is bounded in the monotone convergence region.
So far we retrieve all the conditions for convergence. Therefore, from the proof of the convergence in scalar case, we have that there exists a constant K > 0, s.t., |sk| < 1 and 0 < 2µ− x>k yk < 2µ for all k > K. Then let Sk = h2(2µ− x>k yk) and we have
x>k+1yk+1 − µ = (x>k yk − µ)(1− hu2k − h2x>k yk(µ− x>k yk)) + h2(2µ− x>k yk)((x>k yk)2 − x>k xky>k yk)
= (x>k yk − µ)sk + Sk(U2k − VkWk)
In fact x>k yk − µ → 0 as k → 0. First, since Sk(U2k − VkWk) → 0 as k → 0 by Lemma H.9 and |sk| < 1, we have x>k yk − µ is bounded. Also, from previous discussion, |sk| ≤ max{C2, 1 − C3(x > k yk − µ)2} for some constant C2, C3 > 0. For any 0 > 0, there exists K1 > 0 and 1 0 such that |Sk(U2k − VkWk)| < 1 for all k > K1 and C2 ≤ 1 − C3 21. Therefore x>k yk − µ will
decreases to O( 1). This is because
|x>k+1yk+1 − µ| ≤ (1− C3 21)|x>k yk − µ|+ 2µ|U2K1 − VK1WK1 |C 2(k−K1) 0
≤ |(x>K1yK1 − µ)|(1− C3 2 1) k−K1
+ 2h2µ(VK1WK1 − U2K1) k∑
i=K1
C 2(i−K1) 0 (1− C3 21)k−i
≤ |(x>K1yK1 − µ)|(1− C3 2 1) k−K1 + 2h2µ(VK1WK1 − U2K1)(k −K1)C k−K1 4
where C4 = max{C0, 1 − C3 21} < 1 and we have kCk4 → 0 as k → ∞. Therefore there exists K2 ≥ K1 s.t. x>K2yK2 − µ = O( 1). Then we will show for all k > K2, |x > k yk − µ| < 0. As is shown before there exists K3 > K2 and 2 1 s.t. |Sk(U2k − VkWk)| < 2. Since |sk| < 1, the increase of |x>k yk − µ| is also O( 1) and thus |x>k yk − µ| = O( 1) < 0 for k > K2. Then we have the convergence of |x>k yk − µ| → 0.
In general, we have x>k yk → µ and | cos(∠(xk, yk))| → 0, which implies GD converges to the global minimum x>y = ‖x‖ ‖y‖ = µ.
For the proof of alignment, we only need to consider when ‖xk‖2 +‖yk‖2 ≤ 2h + 2hµ 2(1 + 11−h2µ2 ). Also, we have VkWk − U2k ≥ 0 and
Vk+1Wk+1 − U2k+1 = (VkWk − U2k )(1− h(Vk +Wk) + h2(VkWk − µ2))2.
Thus we denote lk = 1 − h(Vk + Wk) + h2(VkWk − µ2) which characterizes the change of VkWk − U2k .
Lemma H.9. If ‖xk‖2 + ‖yk‖2 ≤ 2h + 2hµ 2(1 + 11−h2µ2 ), then VkWk − U 2 k converges to 0 and also | cos(∠(xk, yk))| converges to 1.
Proof. By 0 ≤ VkWk ≤ (Vk+Wk2 ) 2 and the boundedness theorem, we have
Vk +Wk ≤ 2
h + 2hµ2(1 +
1
1− h2µ2 ),
and
lk = 1− h(Vk +Wk) + h2(VkWk − µ2) ≥ −1− h2µ2(3 + 2
1− h2µ2 ), (11)
lk = 1− h(Vk +Wk) + h2(VkWk − µ2) ≤ h2
4 (Vk +Wk)
2 − h(Vk +Wk) + 1− h2µ2
≤ 1− h2µ2 < 1.
We will show that there exists K > 0 and a constant C > 0 only depending on hµ, s.t. when k > K, |lk| < 1− C.
If Vk +Wk < 2h − hµ 2, then
lk = 1− h(Vk +Wk) + h2(VkWk − µ2) > −1 + h2VkWk > −1.
(If VkWk = 0, then VkWk − U2k = 0 meaning it just converge.)
Therefore when Vk + Wk < 2h − hµ 2, we have |lk| < 1. We only need to consider lk+1 when 2 h − hµ 2 ≤ Vk +Wk ≤ 2h + 2hµ 2(1 + 11−h2µ2 ). Next we replace lk, Vk,Wk, Uk by l, V,W,U for
simplicity.
lk+1 = 1− h(Vk+1 +Wk+1) + h2(Vk+1Wk+1 − µ2) = 1− h(Vk+1 +Wk+1) + h2(Vk+1Wk+1 − U2k+1) + h2(U2k+1 − µ2) = 1− h(Vk+1 +Wk+1) + h2l2(VW − U2) + h2(U2k+1 − µ2) (10) = l + h2(3 + l2 − h(1 + 4h2µ2)(V +W ))(VW − U2)
+ h2U2(3 + (l + 2h2µ2)2 − h(1 + 4h2µ2)(V +W )) + h2µ2(h2(W + V )2 − h(W + V ) + 4h4V 2W 2) + 2h2Uµ((1− (l + 2h2µ2))h(W + V )− 2 + 2(l + 2h2µ2)(1− l − h2µ2)) (12)
If l ≥ −4h2µ2, we have lk+1 = l + h
2(3 + l2 − h(1 + 4h2µ2)(V +W ))(VW − U2) + h2U2(3 + (l + 2h2µ2)2 − h(1 + 4h2µ2)(V +W )) + h2µ2(h2(W + V )2 − h(W + V ) + 4h4V 2W 2) + 2h2Uµ((1− (l + 2h2µ2))h(W + V )− 2 + 2(l + 2h2µ2)(1− l − h2µ2))
≥ l + 2h2Uµ((1− (l + 2h2µ2))h(W + V )− 2 + 2(l + 2h2µ2)(1− l − h2µ2)) ≥ l − h2µ(V +W )(5h2µ2 + (9h4µ4)/8)
≥ −4h2µ2 − hµ(2 + 2h2µ2(1 + 1 1− h2µ2 ))(5h2µ2 + (9h4µ4)/8)
> −1,
where the first inequality is because all the removed terms are positive when 2h − hµ 2 ≤ V +W ≤ 2 h + 2hµ
2(1 + 11−h2µ2 ); the second inequality follows from recollecting l and take l from quadratic formula with the h(V +W ) taken at the bound, and U > −V+W2 ; the third inequality follows from the bounds of h(V +W ) and l ; the last inequality is from hµ ≤ 1
2 √ 7 .
If l < −4h2µ2 ⇒ l + h2µ2 < 0. By rewriting lk+1 in the following way lk+1 = l + h
2(4h4V 2W 2 − h(V +W ) + h2(V +W )2)µ2 + 4h4U2µ2(l + h2µ2) + 2h2Uµ(−2− h(V +W )(−1 + l + 2h2µ2)− 2(−1 + l + h2µ2)(l + 2h2µ2)) + h2VW (3 + l2 − h(V +W )(1 + 4h2µ2)),
we have that the minimum of lk+1 is achieved when U = ± √ VW . Thus we will consider two cases: U ≥ 0 and U < 0. Note we always have U2 < VW (otherwise VW − U2 = 0 just converges) and hence lk+1 is larger than the corresponding value at U = ± √ VW . Let q = h|U |, v = hµ, q = cv for c ≥ 0. Then we have
q2 = c2v2 ≤ h2VW = h(V +W )− 1 + l − v2 < 1, i.e., c < 1 v . (13)
Next, let
l′k+1 = lk+1−h2(3 + l2−h(1 + 4h2µ2)(V +W ))(VW −U2) < lk+1−h2( 1
2 + l2)(VW −U2).
First let U ≥ 0 and q = hU < h √ VW . Choose q = hU = h √ VW . Then by
h(V +W ) = 1− l + h2(VW − µ2) = 1− l + c2v2 − v2, we have the minimum value of lk+1 when U ≥ 0
lk+1 ≥ l′k+1 = l + (−1 + c)(−((−1 + l)l) + c(2 + l + l2))v2 (14) − (−1 + c)2(1 + c2 − 2l + 2cl)v4 + (−1 + c)4v6.
1) When −1 < l < −4h2µ2 and c = 1, we have lk+1 = l′k+1 = l = lk.
2) When −1 < l < −4h2µ2 and c > 1, the derivative of (14) with respect to c is d
dc l′k+1 = (1− l)lv2 + (−2− l − l2)v2 + 2v4 − 6lv4 − 4v6 + 3c2(2v4 − 2lv4 − 4v6)
+ 4c3(−v4 + v6) + 2c((2 + l + l2)v2 − 2v4 + 6lv4 + 6v6), which is a third order polynominal with negative coefficient of the highest degree term and has a root in [0, 1]. Also, when c = 1, ddc l ′ k+1 = 2v 2(1 + l) > 0. Therefore, when 1 < c < 1v , l ′ k+1 either increases or first increases then decreases, namely, the minimum of l′k+1 is achieved at c = 1 or c = 1v . If c = 1, by the previous discussion, l ′ k+1 = l. If c = 1 v ,
l′k+1 = 1 + l 2(−1 + v)2 − v2 − 2v3 + 5v4 − 4v5 + v6 + l(2− 2v + 5v2 − 6v3 + 2v4)
≥ v 2(−24 + 44v − 25v2 + 4v3)
4(−1 + v)2 > −1,
where the first inequality follows from the property of quadratic function by taking l = − (2−2v+5v
2−6v3+2v4) 2(−1+v)2 . Therefore when −1 < l < −4h 2µ2 and c > 1, lk+1 > −1.
3) When −1 < l < −4h2µ2 and c < 1, rewrite l′k+1 with respect to l and we have
l′k+1 = 2(−1 + c)cv2 − (−1 + c)2v4 − (−1 + c)2c2v4 + (−1 + c)4v6 + l2((1− c)v2
+ (−1 + c)cv2) + l(1 + (−1 + c)v2 + (−1 + c)cv2 + 2(−1 + c)2v4 − 2(−1 + c)2cv4), where the minimum point is
l = − 1 2(−1 + c)2v2 − 1 2 − (1− c)v2 < − 1 2v2 − 1 2 < −1.
Therefore the minimum of l′k+1 is obtained at l = −1, i.e.,
l′k+1 ≥ −1 + (−1 + c)2v2(2− (3− 2c+ c2)v2 + (−1 + c)2v4) ≥ −1 + (−1 + c)2v2(2− 3v2 + (−1 + c)2v4) > −1,
where the inequality is from 0 ≤ c < 1. Therefore when −1 < l < −4h2µ2 and c < 1, lk+1 > −1.
4) When l ≤ −1 and c ≥ 54 , consider l ′ k+1 − l (15) and take the derivative ddc (l ′ k+1 − l) = ddc l ′ k+1. From previous discussion, consider ddc l ′ k+1 at c = 5 4 , i.e.,
d
dc l′k+1 = 1/16v 2(48 + 8l2 − 23v2 + v4 + l(40− 6v2)) > v2 − (161v4)/16 + (325v6)/16 > 0,
where the first inequality follows from l > −1− 6v2 by (11). Therefore the minimum of l′k+1 − l is at c = 54 or c = 1 v . If c = 5 4 ,
l′k+1 − l = v2 256 (160 + 16l2 − 41v2 + v4 − 8l(−18 + v2)) > v 2 256 (32− 705v2 + 625v4) > 0
where the first inequality is from l > −1− 6v2 by (11). If c = 1v ,
l′k+1 − l = l2(−1 + v)2 + l(−1 + v)(−1 + v − 4v2 + 2v3) + (−1 + v)(−1− v + 2v3 − 3v4 + v5) ≥ (1− v)3(1 + 3v + v2 − v3) > 0
where the first inequality is from l ≤ −1. Therefore when l ≤ −1 and c ≥ 54 , lk+1 − lk > min{ v 2
256 (32− 705v 2 + 625v4), (1− v)3(1 + 3v + v2 − v3)} > 0.
In the next two cases c ≤ 1 and 1 < c < 54 , we no longer assume q = h √ VW but want the same
lower bound of l′k+1. Since q = hU < h √ VW , we have h(V +W ) > 1− l + c2v2 − v2 > 0. By rewriting l′k+1 to be
l′k+1 =l − 4cv2 + 3c2v2 + 4c4v6 + 4cv2(1− l − v2)(l + 2v2) + v2h2(V +W )2
+ c2v2(l + 2v2)2 + (−v2 + 2cv2(1− l − 2v2)− c2v2(1 + 4v2))h(V +W ),
it can be seen from the minimum point w.r.t. h(V +W ) that when h(V +W ) > 0 and 0 ≤ c < 54 , l′k+1 decreases as h(V +W ) decreases. Therefore by h(V +W ) > 1− l + c2v2 − v2, we have the same expression for the lower bound of l′k+1 but with U 2 < VW , i.e.,
lk+1 − l′k+1 = h2(3 + l2 − h(1 + 4h2µ2)(V +W ))(VW − U2) > h2( 1
2 + l2)(VW − U2),
l′k+1 > l + (−1 + c)(−((−1 + l)l) + c(2 + l + l2))v2
− (−1 + c)2(1 + c2 − 2l + 2cl)v4 + (−1 + c)4v6.
5) When l ≤ −1 and c ≤ 1, consider l′k+1 − l
l′k+1 − l > (−1 + c)(−((−1 + l)l) + c(2 + l + l2))v2 (15) − (−1 + c)2(1 + c2 − 2l + 2cl)v4 + (−1 + c)4v6
= (−1 + c)2l2v2 + (−1 + c)lv2(1 + c− 2(1− c)v2 + 2(1− c)cv2) + (−1 + c)v2(2c+ (1− c)v2 + (1− c)c2v2 + (−1 + c)3v4)
where the minimum point is l = − 12 + 1−(1−c)v2+(1−c)cv2 1−c > −1. By l ≤ −1 we have
l′k+1 − l > v2(2(−1 + c)2 − (−1 + c)2(3− 2c+ c2)v2 + (−1 + c)4v4) = v2((−1 + c)2(2− (3− 2c+ c2)v2) + (−1 + c)4v4) ≥ v2((−1 + c)2(2− 3v2) + (−1 + c)4v4) ≥ 0.
Therefore when l ≤ −1 and c ≤ 1, lk+1 − lk > h2( 12 + l 2)(VW − U2) ≥ 32h 2(VW − U2).
6) When l ≤ −1 and 1 < c < 54 , consider Uk+1. By h(V +W )− h 2VW = 1− l − v2, we have
Uk+1 = h(1− hV )µW + hµV (1− hW ) + U(1 + h2VW − h(V +W ) + h2µ2) = µ((c− 1)(l + v2) + 1 + cv2 − h2VW ).
Thus denote hUk+1 = ck+1hµ = ck+1v and we have
ck+1 = (c− 1)(l + v2) + 1 + cv2 − h2VW < (c− 1)(l + v2) + 1 + cv2 − c2v2 < 1.
Then if lk+1 − l ≥ 0, then lk+2 − l = lk+2 − lk+1 + lk+1 − l > lk+2 − lk+1 > 32h 2(VW − U2). Otherwise, we have lk+1 < l and consider lk+2. Also, the distance between 1 and ck+1 is larger than that of 1 and c, i.e.,
1− ck+1 − (c− 1) > (1− (c− 1)(l + v2) + 1 + cv2 − c2v2)− (c− 1) (16) = 1 + l + v2 + c2v2 + c(−1− l − 2v2) > 0.
From previous discussion, l′k+1 − l increases as l decreases. Also ddc l ′ k+1 < 0 when 0 ≤ c < 1, so l′k+1 − l increases as c decreases. By (16), we have ck+1 < 2− c. Then
l′k+2 − l = l′k+2 − lk+1 + lk+1 − l > l′k+2 − lk+1 + l′k+1 − l > (−1 + c)(−((−1 + l)l) + c(2 + l + l2))v2
− (−1 + c)2(1 + c2 − 2l + 2cl)v4 + (−1 + c)4v6
+ (1− c)v2(−((−1 + l)l)− (−2 + c)(2 + l + l2) + (−1 + c)(5 + c2 + 2l − 2c(2 + l))v2 − (−1 + c)3v4)
= 2v2(−1 + c)2((2 + l + l2) + (−2− (−1 + c)2)v2 + (−1 + c)2v4) > 2v2(−1 + c)2((2− 3v2) + (−1 + c)2v4) > 0,
where the second inequality follows from ck+1 < 2 − c and lk+1 < l; the last inequality is by l ≤ −1 and 1 < c < 54 . Therefore when l ≤ −1 and 1 < c < 5 4 , lk+2 − lk > 3 2h
2(VW − U2) or 2v2(−1 + c)2((2− 3v2) + (−1 + c)2v4).
Second, let U < 0 and q = −hU < h √ VW . Choose q = −hU = h √ VW . Then similarly we have the minimum value of lk+1 when U < 0
l′k+1 = l + (1 + c)((−1 + l)l + c(2 + l + l2))v2 − (1 + c)2(1 + c2 (17) − 2l − 2cl)v4 + (1 + c)4v6
= (1 + c)2l2v2 + l(1 + (−1 + c2)v2 + 2(1 + c)3v4) + (1 + c)v2(v2(−1 + v2) + c3v2(−1 + v2) + c2v2(−1 + 3v2) + c(2− v2 + 3v4))
≥ −1 + (−4 + 4c+ 3c2)v2 + (15 + 16c− 5c2 − 8c3 − 2c4)v4 (18) + (−5− 4c+ 2c2 + c3)2v6,
where the inequality is because when l > −1− 6v2 + c2v2 (from (11)), l′k+1 increases as l increases.
1) When c ≥ 1, since (−5− 4c+ 2c2 + c3)2v6 ≥ 0, we consider (−4 + 4c+ 3c2)v2 + (15 + 16c− 5c2 − 8c3 − 2c4)v4 denoted to be lpartk+1
d
dc lpartk+1 = (4 + 6c)v 2 + (16− 10c− 24c2 − 8c3)v4.
It has a root between -2 and 0. Also when c = 1, ddc l part k+1 > 0; when c = 1 v , d dc l part k+1 < 0. Thus, the minimum of lpartk+1 is achieved at c = 1 or c = 1 v , i.e., l part k+1 ≥ v2(3 + 16v2). Then
l′k+1 ≥ −1 + v2(3 + 16v2) + (−5− 4c+ 2c2 + c3)2v6 > 0.
Therefore when c ≥ 1, lk+1 > −1.
2) When 0 ≤ c < 1 and l > −1, from previous discussion of l′k+1 we have
l′k+1 ≥ −1 + 2(1 + c)2v2 − (1 + c)2(3 + 2c+ c2)v4 + (1 + c)4v6
≥ −1 + (1 + c)2v2(2− 6v2) + (1 + c)4v6
≥ −1 + v2(2− 6v2) + v6 > −1.
Therefore when 0 ≤ c < 1 and l > −1, lk+1 > −1.
3) When 0 ≤ c < 1 and l ≤ −1, consider l′k+1 − l. By (17),
l′k+1 − l = (1 + c)2l2v2 + l((−1 + c2)v2 + 2(1 + c)3v4) + (1 + c)v2(v2(−1 + v2) + c3v2(−1 + v2) + c2v2(−1 + 3v2) + c(2− v2 + 3v4))
≥ 2(1 + c)2v2 − (1 + c)2(3 + 2c+ c2)v4 + (1 + c)4v6
≥ v2(2− 6v2) + v6,
where the first inequality is because the minimum point of l is greater than -1. Therefore when 0 ≤ c < 1 and l ≤ −1, lk+1 − lk > v2(2− 6v2) + v6.
Actually all the proofs of lk > −1 ⇒ lk+1 > −1 + C for some C > 0 are valid for all 0 < Vk +Wk ≤ 2h + 2hµ 2(1 + 11−h2µ2 ).
In general, if lk > −1, then there exists a constant C > 0 only depending on hµ s.t. li > −1 +C for all i > k. If lk ≤ −1, then either (1) lk+1 increases at least by a fixed value only depending on hµ, or (2) lk+2 or lk+1 increases by 32h
2(VkWk − U2k ). If it keeps staying in case (2), then since lk ≤ −1, VkWk − U2k will be larger and larger until li > −1 for some i. Therefore, there exists a K > 0, s.t. when k > K, lk > −1 +C. Then because lk ≤ 1−h2µ2, let 0 < C0 = max{1−C, 1−h2µ2} < 1 and we have |lk| ≤ C0 < 1 for all k > K. Further,
VkWk − U2k ≤ C 2(k−K) 0 (VKWK − U2K)→ 0 as k →∞,
i.e., 1− U 2 k
VkWk → 0⇒ | cos(∠(xk, yk))| → 1.
Otherwise, there exists a K > 0, s.t., VKWK − U2K = 0, then VkWk − U2k = 0 for all k ≥ K. There are two cases. (i) xK and yK are already aligned. (ii) one of VK and WK is 0, WOLG assume
VK = 0. Then from the GD update,{ xk+1 = xk + h(µI − xky>k )yk yk+1 = yk + h(µI − ykx>k )xk ,
xK+1 and yK+1 will be aligned for both cases. Then for all k ≥ K + 1, xk and yk is aligned, i.e., | cos(∠(xk, yk))| = 1.

heading: I PROOF OF THEOREM 5.1

heading: I.1 GENERAL RANK 1 APPROXIMATION FOR NON-NEGATIVE DIAGONAL MATRIX
text: We first discuss an easy case where A is a non-negative diagonal matrix (later on we will see this is actually the canonical case) and consider
min x,y∈Rn
‖A− xy>‖2F /2. (19)
Assume A = µ1In1 . . . µmInm , where Ini is an identity matrix of dimension ni × ni, ∑m i=1 ni = n, µi ≥ 0, µi 6= µj for i 6= j. Let Si = {1, 2, · · · , ni} + ∑i−1 j=1 nj , i.e., an index set describing the positional indices of Ini . Let µmax = max{µ1, · · · , µm} and Smax be the corresponding index set. Let [x∞, y∞] be a fixed point of the GD map. Let xs,∞ and ys,∞ be the sth element of x∞ and y∞.
The following theorem characterizes the fixed points of the GD map. Theorem I.1. The fixed points of GD map satisfy ∑ s∈Si x 2 s,∞· ∑ s∈Si y 2 s,∞ = µ 2 i , xs,∞ = ys,∞ = 0 for s /∈ Si, ∀i = 1, · · · , n. Thus ‖x∞‖‖y∞‖ = µi, ∀i = 1, · · · , n.
Proof. To solve the fixed points, we have ∀i = 1, · · · , n and x ∈ Si{ (A− x∞y>∞)y∞ = Ay∞ − (y>∞y∞)x∞ = 0 (A− x∞y>∞)>x∞ = A>x∞ − (x>∞x∞)y∞ = 0 i.e., { µiys,∞ = (y > ∞y∞)xs,∞
µixs,∞ = (x > ∞x∞)ys,∞
⇒ µ2ixs,∞ys,∞ = (x>∞x∞)(y>∞y∞)xs,∞ys,∞
Obviously, x∞ = y∞ = 0 is a fixed point. If xs,∞ 6= 0, then ys,∞ 6= 0 and we have µ2i = (x>∞x∞)(y > ∞y∞). Since µi 6= µj for i 6= j, there exists at most one i, s.t. µ2i = (x>∞x∞)(y>∞y∞) and all the xs,∞, ys,∞ s.t. s ∈ Sj , j 6= i is zero.
This theorem identifies that each of these fixed points is concentrated in one of the positions of the eigenvalue blocks. Note the fixed points contain both stable and unstable ones, i.e., all the global minima and saddles are of the above form. Moreover, the global minima are the points obeying ∑ s∈Smax x 2 s,∞ · ∑ s∈Smax y 2 s,∞ = µ 2 max, xs,∞ = ys,∞ = 0 for s /∈ Smax and also ‖x∞‖‖y∞‖ = µmax. The saddles are the ones with µi 6= µmax, ∀i = 1, · · · , n. Remark I.2 (The alignment of x∞ and y∞ at the global minimum). At the global minimum of the objective (19) with diagonal and non-negative A defined above, x∞y>∞ is the rank-1 approximation of A with the largest eigenvalue, i.e., we have µmax = tr(x∞y>∞) = x > ∞y∞. Also, from the above theorem, ‖x∞‖‖y∞‖ = µmax. Therefore for each global minimum, there exists a scalar l > 0, s.t., x∞ = ly∞, i.e., x∞ and y∞ are aligned at the global minimum.
Based on the analytical form of eigenvalues of the Jacobian of GD map (19) in Theorem I.1, we can establish the following relation between x∞ and y∞ at the global minimum via stability analysis.
Theorem I.3. For almost all initial conditions, if GD converges to a minimum, then x>∞y∞ = µmax and
‖x∞‖2 + ‖y∞‖2 < 2
h ; (20)
the extent of balancing is quantified by
‖x∞ − y∞‖2 < 2
h − 2µmax. (21)
Proof. Since each unstable fixed point has its stable set with negligible size when compared to the stable set of stable fixed point, then for almost all initial conditions, if GD converges to a global minimum, this minimum is a stable fixed point of the GD map.
From the above theorem, we know the fixed points are µ2i = (x > ∞x∞)(y > ∞y∞). Assume µ 2 = (x>∞x∞)(y > ∞y∞).
For the Jacobian J = I2n + ( −h(y>∞y∞)In hA− 2hx∞y>∞ hA> − 2hy∞x>∞ −h(x>∞x∞)In ) = I2n + Ã, since the lower two blocks commute, we have
det(λI − Ã) = det(λ2In + λh(x>∞x∞ + y>∞y∞)In + h2µ2In − h2A2
+ h2(2Ay∞x > ∞ + 2x∞y > ∞A− 4x∞y>∞y∞x>∞))
= det(λ2In + λh(x > ∞x∞ + y > ∞y∞)In + h 2µ2In − h2A2
+ h2(2(A− x∞y>∞)y∞x>∞ + 2x∞y>∞(A− y∞x>∞))) = det(λ2In + λh(x > ∞x∞ + y > ∞y∞)In + h 2µ2In − h2A2)
= n∏ i=1 (λ2 + h(x>∞x∞ + y > ∞y∞)λ+ h 2(µ2 − µ2i )).
Hence the two values, 1 and 1 − h(x>∞x∞ + y>∞y∞), are eigenvalues of J at each nonzero fixed point.
If µ = max{µi}, then |1− h(x>∞x∞ + y>∞y∞)| < 1⇒all the eigenvalues are bounded by 1. Hence x>∞x∞ + y > ∞y∞ < 2 h .
If µ 6= max{µi}, then 1+ 12 (−h(x > ∞x∞+y > ∞y∞)+ √ h2(x>∞x∞ + y > ∞y∞)
2 + 4h2((µ2i − µ2))) > 1 for µi ≥ µ. It has at least one unstable direction. It’s stable set is measure 0.
The second inequality of the theorem follows from x>∞y∞ = µ at the global minimum because of diagonal and non-negative A and SVD.

heading: I.2 GENERAL MATRIX FACTORIZATION FOR NON-NEGATIVE DIAGONAL MATRIX
text: In this section, we consider problem (5) via GD iteration but with A ∈ Rn×n to be an arbitrary non-negative diagonal matrix.
Let µ1 ≥ µ2 ≥ · · · ≥ µn ≥ 0 be the diagonal elements of A. Let r = rank(A) which of course satifies r ≤ n. Also assume that ‖A‖F does not scale with n or d (without loss of generality we can assume ‖A‖F = O(1)). Then we will have our balancing result in the following. Theorem I.4. For almost all initial conditions, if GD for (5) converges to a global minimum, then there exists c = c(n, d) > c0 with constant c0 > 0 independent of n, d, and learning rate h, such that, (X,Y ) satisfies
c(‖X∞‖2F + ‖Y∞‖ 2 F) <
2 h ,
and the extent of balancing is quantified by
‖X∞ − Y∞‖2F < 2
ch − 2 min{d,r}∑ i=1 µi.
Proof. Let r =rank(A) ≤ n. Since each unstable fixed point has its stable set with size negligible when compared to the stable set of stable fixed point, then for almost all initial conditions, if GD converges to a global minimum, this minimum is a stable fixed point of the GD map.
For matrix X = (x1 x2 · · · xd) ∈ Rn×d, where xi ∈ Rn is the ith column vector, we vectorize the
variable and get vec(X) =  x1 x2 ... xd . Then we can rewrite the GD iteration, {
Xk+1 = Xk + h(A−XkY >k )Yk Yk+1 = Yk + h(A > − YkX>k )Xk
⇒ {
vec(Xk+1) = vec(Xk) + h(I ⊗A− T (vec(Yk)vec(Xk)>))vec(Yk) vec(Yk+1) = vec(Yk) + h(I ⊗A− T (vec(Xk)vec(Yk)>))vec(Xk) ,
where ⊗ is the Kronecker product, and T : Rnd×nd → Rnd×nd is a linear operator, s.t.
T  M11 · · · M1d... ... Md1 · · · Mdd   = M > 11 · · · M>1d ... ... M>d1 · · · M>dd  ,
then T (vec(Y )vec(X)>) =  x1y > 1 x2y > 1 · · · xdy>1 x1y > 2 x2y > 2 · · · xdy>2 ... ...
... x1y > d x2y > d · · · xdy>d
 .
The Jacobian of the vectorized GD map is (replace Xk and Yk by X and Y ) I2nd − hM , where
M =
( Y >Y ⊗ In T (vec(Y )vec(X)>) + I ⊗ (XY > −A)
T (vec(X)vec(Y )>) + I ⊗ (Y X> −A) X>X ⊗ In
) .
We would like to obtain an estimation of the eigenvalues of I −hM for stability analysis. Since M is the Hessian of the objective at the global minimum defined, all the eigenvalues are non-negative. Also, since ‖A‖F is independent of h, n and d, we can just assume without loss of generality ‖A‖F = O(1) and then take ‖X‖F and ‖Y ‖F to be O(1) at each global minimum. Then
tr(M) = n(‖X‖2F + ‖Y ‖ 2 F) = O(n).
Next consider tr(M2). We only need the diagonal blocks of M2. The two diagonal blocks are the following
M1 = (Y >Y ⊗ In)2
+ [T (vec(Y )vec(X)>) + I ⊗ (XY > −A)][T (vec(X)vec(Y )>) + I ⊗ (Y X> −A)], M2 = (X >X ⊗ In)2
+ [T (vec(X)vec(Y )>) + I ⊗ (Y X> −A)][T (vec(Y )vec(X)>) + I ⊗ (XY > −A)].
Since we evaluate this matrix at the minimum, we have the gradient equals 0, i.e.,
(A−XY >)Y = 0, (A− Y X>)X = 0.
By the mixed-product property of Kronecker product, we have
M1 = (Y >Y )2 ⊗ In + T (vec(Y )vec(X)>)T (vec(X)vec(Y )>)
+ Id ⊗ [(A−XY >)(A− Y X>)], M2 = (X
>X)2 ⊗ In + T (vec(X)vec(Y )>)T (vec(Y )vec(X)>) + Id ⊗ [(A− Y X>)(A−XY >)].
Then
tr(M2) = 2 ‖X‖2F ‖Y ‖ 2 F + n( ∥∥X>X∥∥2 F + ∥∥Y >Y ∥∥2 F ) + d ∥∥A−XY >∥∥2 F .
Also, this is the global minimum, meaning∥∥A−XY >∥∥2 F = tr((A−XY >)(A− Y X>)) = {∑r i=d+1 µ 2 i , d < r
0, d ≥ r
Therefore, when d ≥ r, we have
tr(M2) = 2 ‖X‖2F ‖Y ‖ 2 F + n( ∥∥X>X∥∥2 F + ∥∥Y >Y ∥∥2 F ) = O(n).
When d < r, we have
tr(M2) = 2 ‖X‖2F ‖Y ‖ 2 F + n( ∥∥X>X∥∥2 F + ∥∥Y >Y ∥∥2 F ) + d ∥∥A−XY >∥∥2 F
< 2 ‖X‖2F ‖Y ‖ 2 F + n( ∥∥X>X∥∥2 F + ∥∥Y >Y ∥∥2 F ) + n ∥∥A−XY >∥∥2 F = O(n).
Therefore, the number of non-zero eigenvalues is N = O(n). Let λmax be the largest eigenvalue of M . Then there exist a constant c0 > 0
λmax ≥ tr(M)/N ≥ c0(‖X‖2F + ‖Y ‖ 2 F).
From stability analysis, we need the absolute values of the eigenvalues of I − hM to be bounded by 1, i.e.,
1− hc0(‖X‖2F + ‖Y ‖ 2 F) ≥ 1− hλmax ≥ −1⇒ c0(‖X‖ 2 F + ‖Y ‖ 2 F) ≤
2 h .
Obviously, if we pick a constant c = c(n, d) for each n and d instead of c0, the above inequality also holds with c(n, d) > c0. For the derivation of the second inequality, since at the global minimum tr(X∞Y >∞) =∑min{d,r} i=1 µi = ∑min{d,n} i=1 µi because for non-negative diagonal A, its SVD of the non-zero part satisfies Unon−zero = Vnon−zero , we have
‖X∞ − Y∞‖2F = ‖X∞‖ 2 F + ‖Y∞‖ 2 F − 2 tr(X∞Y > ∞) ≤
2
ch − 2 min{d,n}∑ i=1 µi.

heading: I.3 FROM DIAGONAL MATRICES TO GENERAL MATRICES
text: Proof of Theorem 5.1. First all the minima of this problem are global minima and all the saddles are unstable fixed points. Since each unstable fixed point has its stable set with size negligible when compared to the stable set of stable fixed point, then for almost all initial conditions, if GD converges to a point, for almost all situations, this point is a stable fixed point of the GD map.
By singular value decomposition (SVD), A = UDV >, where U,D, V ∈ Rn×n, U and V are orthogonal matrices, and D is a non-negative diagonal matrix. Let Xk = URk, Yk = V Sk. Then{
Xk+1 = Xk + h(A−XkY >k )Yk Yk+1 = Yk + h(A−XkY >k )>Xk
⇔ { URk+1 = URk + h(UDV > − URkS>k V >)V Sk
V Sk+1 = V Sk + h(UDV > − URkS>k V >)>URk ⇔ {
Rk+1 = Rk + h(D −RkS>k )Sk Sk+1 = Sk + h(D −RkS>k )>Rk .
Therefore, problem (5) is equivalent to the following problem solved by GD
min R,S∈Rn×d
∥∥D −RS>∥∥2 F /2. (22)
From Theorem I.3 and I.4, we obtain: if GD for (22) converges to a global minimum (R,S), then there exists a constant c > 0 independent of n, r, d, and h, s.t., the limiting minimum satisfies
c(‖R‖2F + ‖S‖ 2 F) ≤
2 h ,
and the extent of balancing is quantified by
‖R− S‖2F ≤ 2
ch − 2 min{d,r}∑ i=1 µi.
Especially, when d = 1, n ∈ N+, c = 1.
Since X = UR, Y = V S ⇒ R = U>X, S = V >Y , we have ‖R‖F = ‖X‖F, ‖S‖F = ‖Y ‖F, and ‖R− S‖F = ∥∥U>X − V >Y ∥∥ F = ∥∥X − (UV >)Y ∥∥ F . The we obtain the first and second
inequalities in Theorem 5.1. Also, by SVD, ∥∥XY >∥∥2 F = ∑min{d,n} i=1 µ 2 i ≤ ‖X‖ 2 F ‖Y ‖ 2 F, we obtain the third inequality.

Review summaries:
This paper studies the properties of gradient descent with large learning rate in the matrix factorization problem. The goal is to understand when gradient descent converges to a global minimum where the two factors are roughly balanced in norm, which is a shallow minimum that may be hoped to generalize better. For small learning rates (i.e. where classical optimization results guarantee convergence to a critical point), this does not happen. Hence, this paper considers learning rates beyond this regime.  For two special cases (when the matrices are 1xd or dx1), it is shown that under a large learning rate, gradient descent still converges to a global minimum, and this minimum is roughly "balanced". For the general case, it is shown that for almost any initial point, if the learning rate is such that gradient descent converges to a global minimum, then this minimum is roughly balanced.

The authors study the problem of minimizing matrix factorization squared error with gradient descent. They show convergence with larger learning rates than previously considered. They further show that gradient descent has an implicit bias towards finding matrix factors which are magnitude-balanced in a certain sense.  

This paper is concerned with the non-convex matrix factorization problem. This paper considers the square loss, which is optimized via gradient descent. Due to symmetries of the problem, there are infinitely many different global minima. This paper aims to understand to which minima gradient descent with large learning rate (=step size) converges to.   The main claim of this paper is that gradient descent with large learning rate has an implicit regularization effect, meaning that it converges to factors of the underlying matrix, which are implicitly balanced in norm. To substantiate the claim, the paper provides several theoretical results as well as some numerical experiments.  

This paper studies gradient descent optimization of matrix factorization (various forms, including the most generic case) under large learning rate. The authors prove convergence results that apply when the learning rate is larger than 2/L and up to ~4/L, and also show that optimization with large learning rates leads to dynamic balancing between the two matrix factors in the factorization, an effect that does not arise when using small learning rates. In proving some of the theorems, the authors observe the optimization can have two phases: one in which gradient descent is driven to search for flat regions, and a second phase near the global minimum where convergence occurs towards a balanced solution for the factorization problem.