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

Introduction

This paper studies the asymmetric low-rank matrix factorization problem:

where $\eta>0$ is the learning rate and $\mathbf{U}_{0},\mathbf{V}_{0}$ are randomly initialized according to some distribution. Empirically, gradient descent with a constant learning rate can efficiently solve this problem (see, e.g., Figure 1 in Du et al. (2018)). Somehow surprisingly, there is no global convergence proof of this generic algorithm, let alone convergence rate analysis. The main difficulties are 1) the problem is non-convex and 2) this problem is not smooth with respect to $(\mathbf{U},\mathbf{V})$ because the magnitudes of them can be highly unbalanced.

To motivate the study of gradient descent for this optimization problem, we note that this is a prototypical optimization problem that illustrates the gap between practice and theory. In particular, the prediction function $\mathbf{U}\mathbf{V}^{\top}$ is homogeneous: if we multiply a factor by a scalar $c$ and divide another factor by $c$, the prediction function remains the same. This homogeneity also exists in deep learning models. Therefore, progress made in understand (1) can further help us gain understanding on other non-convex problems, such as asymmetric matrix sensing, asymmetric matrix completion, and deep learning optimization. We refer readers to Du et al. (2018) for more discussions.

For Problem (1), Du et al. (2018) showed gradient flow (gradient descent with the step size $\eta\rightarrow 0$),

However, to prove a polynomial convergence rate, the approach that solely relies on the geometry will fail because there exists a counter example (Du et al., 2017). Furthermore, for gradient descent with $\eta>0$, the key invariance no longer holds. While the invariance can still hold approximately in some way, characterizing the approximation error is highly non-trivial, and this is one of our key technical contributions.

Du et al. (2018) also studied gradient descent with decreasing step sizes $\eta_{t}=O\left(t^{-1/2}\right)$, and obtained an “approximate global optimality result": if the magnitude of the initialization is $O\left(\delta\right)$, then gradient descent converges to a $\delta$-optimal solution, i.e., this result does not establish that gradient descent converges to a global minimum. And again, there was no convergence rate. Furthermore, their result crucially relies on $\eta_{t}$ is of order $O\left(t^{-1/2}\right)$ to ensure the second order term does not diverge and thus does not apply to gradient descent with a constant learning rate.

Some previous works, e.g., Ge et al. (2015); Jin et al. (2017), modified the gradient descent algorithm to the perturbed gradient descent algorithm by adding an isotropic noise at each iteration, which can help escape strict saddle points and bypass the exponential lower bound in Du et al. (2017). To deal with the non-smooth problem, they also added a balancing regularization term (Park et al., 2017; Tu et al., 2016; Ge et al., 2017a; Li et al., 2019b), $\frac{1}{8}\|\mathbf{U}^{\top}\mathbf{U}-\mathbf{V}^{\top}\mathbf{V}\|_{F}^{2}$ to the objective function to ensure balancedness between $\mathbf{U}$ and $\mathbf{V}$ throughout the optimization process. With these two modifications, one can prove a polynomial convergence rate. However, experiments suggest that the isotropic noise and the balancing regularizer may be proof artifacts, because vanilla gradient descent applies to the original objective function (1) without any regularizer finds a global minimum efficiently. From a practical point of view, one does not want to add noise or additional regularization because it may require more hyper-parameter tuning.

The only global quantitative analysis for randomly initialized gradient is by Du et al. (2018) who proved the global convergence rate for the case where $\mathbf{\Sigma}$ has rank $1$, and $\mathbf{U}$ and $\mathbf{V}$ are two vectors. In this case, one can reduce the problem to the dynamics of $4$ variables, which can be easily analyzed. Unfortunately, it is very difficult to generalize their analysis to the general rank setting.

In this paper, we develop new techniques to overcome the technical difficulties and obtain the first polynomial convergence of randomly initialized gradient descent for solving the asymmetric low-rank matrix factorization problem. Most importantly, our analysis is completely different from existing ones: we give a thorough characterization of the entire trajectory of gradient descent.

Before presenting our main results, we emphasize that the goal of this paper is not to provide new provably efficient algorithms to solve Problem (1), but to provide a rigorous analysis of an intriguing and practically relevant phenomenon on gradient descent. This is of the same flavor as the recent breakthrough on understanding Burer-Moneiro method for solving semidefinite programs (Cifuentes and Moitra, 2019).

Here, $\sigma_{1}$ and $\sigma_{d}$ are the largest and the smallest singular values of $\mathbf{\Sigma}$, respectively. Notably, in sharp contrast to the result in Du et al. (2018), which requires the initialization depends on $\delta$, our initialization does not depend on the target accuracy. To our knowledge, this is the first global convergence result for gradient descent in solving Problem (1). Furthermore, we give a polynomial rate. The first term in $T_{\text{total}}(\delta,\eta)$ represents a warm-up phase and the second term represents the local linear convergence phase, which will be clear in the analysis sections. On the other hand, while we believe $T_{\text{total}}(\delta,\eta)$ is nearly tight, our requirement for $\eta$ is loose. An interesting future direction is further relax this requirement.

Now by taking $\eta\rightarrow 0$, we have the following corollary for gradient flow.

Given $\delta>0$, there exists $T=O\left(\frac{1}{\sigma_{d}}\ln\frac{d\sigma_{d}}{\varepsilon}+\frac{1}{\sigma_{d}}\ln\frac{\sigma_{d}}{\delta}\right)$, such that with high probability over the initialization, for all $t\geq T$, we have $f\left(\mathbf{U}_{t},\mathbf{V}_{t}\right)\leq\delta.$ In gradient flow, $t$ is a continuous time index.

This is also the first convergence rate result of randomly initialized gradient flow for asymmetric matrix factorization. We note that our analysis on gradient flow is nearly tight. To see this, consider the ordinary differential equation $\dot{a}_{t}=(\sigma_{d}-a_{t}^{2})a_{T}$ with initial point $a_{0}>0$, then $s=a^{2}$ has analytical solution $s_{t}=\frac{\sigma_{d}e^{2\sigma_{d}t}}{e^{2\sigma_{d}t}+\frac{\sigma_{d}}{a_{0}^{2}}-1}$. Hence, to achieve a $\delta$ optimal solution, i.e. $|\sigma_{d}-a^{2}|\leq\delta$, we need $\sigma_{d}\left(\frac{\sigma_{d}}{a_{0}^{2}}-1\right)\frac{1}{\delta}\leq e^{2\sigma_{d}t}+\frac{\sigma_{d}}{a_{0}^{2}}-1$. Hence $T=\Theta\left(\frac{1}{\sigma_{d}}\ln\frac{\sigma_{d}}{\delta}\right)$ is necessary.

2 Additional Related Work

Here we discuss additional related work. First, in the symmetric setting, e.g., $\min_{\mathbf{U}}\|\mathbf{U}\mathbf{U}^{\top}-\mathbf{\Sigma}\|_{F}^{2}$, global convergence of randomly initialized gradient has been established in various settings (Jain et al., 2017; Li et al., 2018; Chen et al., 2019).In Appendix B, we show the dynamics of gradient flow actually admits a closed form, and thus can be easily analyzed. However, as has been highlighted in Li et al. (2019a, b); Park et al. (2017); Tu et al. (2016), generalization to the asymmetric case is highly non-trivial. The major technical difficulty is to deal with the unbalancedness between $\mathbf{U}$ and $\mathbf{V}$. To prevent this, additional balancing regularization is often added (Li et al., 2019b; Park et al., 2017; Tu et al., 2016; Sun and Luo, 2016), though empirically this has been shown to be unnecessary.

Another line of work showed one can first uses spectral initialization to find a near-optimal solution, then starting from there, gradient descent converges to an optimum with a linear rate (Tu et al., 2016; Zheng and Lafferty, 2016; Zhao et al., 2015; Bhojanapalli et al., 2016), though in practice random initialization often suffices. Recently, Ma et al. (2021) proved that if 1) the initialization is close to a global minimum and 2) $\mathbf{U}$ and $\mathbf{V}$ are balanced, then without adding additional balancing regularizer, gradient descent converges to a global minimum. Our stage two’s analysis is similar to theirs. However, their result cannot be directly applied to our analysis because they require a more stringent initialization than our stage two’s initial point.

Main Difficulties and Technique Overview

The starting point is the Polyak-Łojasiewicz condition: if we can establish that $\max\left\{\sigma_{d}(\mathbf{U}_{t}),\sigma_{d}(\mathbf{V}_{t})\right\}$ is lower bounded by a considerable constant $c_{\max}$, then we have $\left\|\nabla f(\mathbf{U},\mathbf{V})\right\|\geq c_{\max}\sqrt{2f(\mathbf{U},\mathbf{V})}$, which implies a linear convergence. However, the $d^{\text{th}}$ singular values of $\mathbf{U}$ and $\mathbf{V}$ are not monotonic with $t$, and they can even decrease to an extremely small value.

Hence, without loss of generality, we can assume $\mathbf{\Sigma}$ is a diagonal matrix with $\mathbf{\Sigma}_{i,i}=\sigma_{i}$, $\forall i\in[d]$, and $\mathbf{\Sigma}_{i,j}=0$ otherwise.

To proceed, we will analyse the principle space and the complement space separately. We denote the upper $d\times d$ matrix of $\mathbf{U}$ as $U$ and denote the lower $(m-d)\times d$ matrix of $\mathbf{U}$ as $J$. Similarly, we define the upper $d\times d$ matrix of $\mathbf{V}$ as $V$ and the lower $(n-d)\times d$ matrix as $K$. Define $\Sigma:=\text{diag}(\sigma_{1},\cdots,\sigma_{d})$. We can write out the dynamics of these matrices:

Besides $A,B,J$ and $K$, there are some other special capital letters used to represent specific matrices throughout this paper. Here is a list.

We define such $S$ is because in symmetric case ($B\equiv 0$), although it is hard to find analytical solution for $A$ in continuous time case, we do find analytical form for $S$, which contains all information about the singular values of $A$.

$P$ and $Q$ are just the symmetric and skew-symmetric part of matrix $\Sigma-UV^{\top}$. Hence the linear convergence of gradient descent is equivalent the linearly diminishing of $P$ and $Q$ by Pythagorean theorem. We will mention their definitions every time we use them.

2 Symmetrization

Our key observation is that although the singular values of $U$ and $V$ may not have monotonic property, the symmetrized matrix has this property. Formally, we define

Here, $A$ represents the magnitude in the principle space and $B$ represents the magnitude of asymmetry. Empirically, we can observe that by choosing a sufficiently small learning rate $\eta$, we have two desired properties:

The smallest singular value of $A$ is almost monotonically increasing;

The norms of $B,J,K$ are almost monotonically decreasing.

The first property ensures we are learning the “signal", $\Sigma$, and the second property ensures the “noise" is disappearing. Therefore, if we can establish these two properties, we can prove the global convergence.

3 Two Stage Analysis

The analysis for asymmetric low rank case is divided into two stages. In the first stage we mainly focus on the increasing rate of $\sigma_{d}(A)$. We will prove that in gradient descent method $\sigma_{d}(A_{d})$ increases exponentially fast to $\sqrt{\frac{\sigma}{2}}$ and then $\|P\|_{op}$ drops exponentially fast to $\frac{\sigma_{d}}{4}$, while preserving $\|B\|_{F},\|J\|_{op}$ and $\|K\|_{op}$ small. In the second stage, we will use the large $\sigma_{d}(A)$ to lower bound the convergence speed of $\|\mathbf{\Sigma}-\mathbf{U}\mathbf{V}^{\top}\|_{F}^{2}$. We will prove that, once gradient descent starts at a point with small $\|P\|_{op}$, $\|B\|_{F}$, $\|J\|_{op}$ and $\|K\|_{op}$, it will converge to global optimal point exponentially fast.

Proof Sketch of Theorem 1.1

We first use a Gaussian distribution to generate matrices $U,V,J,K$ element-wisely and independentlyStrictly speaking, we cannot make any assumption on $U,V,J,K$ since they need information of singular value decomposition of $\mathbf{\Sigma}$. However, a random generation of $\mathbf{U}$ and $\mathbf{V}$ implies a random generation of $U,V,J,K$ because we use unitary transformations.. By standard random matrix theory (Corollary 2.3.5 and Theorem 2.7.5 of Tao (2012)), we know that $\exists c>0$, such that with high probability, the smallest singular value of $\frac{{U}+{V}}{2}$ is larger than $\frac{1}{c\sqrt{d}}$, the largest singular value of $\frac{{U}+{V}}{2}$ is smaller than $c\sqrt{d}$, the Frobenius norm of $B$ is less than $cd$ and the operator norms of $J$ and $K$ are less than $c\sqrt{\max\{m^{\prime},d\}}$ and $c\sqrt{\max\{n^{\prime},d\}}$, respectively, where $m^{\prime}=m-d,n^{\prime}=n-d$.

The initializations $U_{0},V_{0},J_{0},K_{0}$ are then scaled by $\varepsilon$ where $\varepsilon$ specified in Theorem 1.1.

2 Stage One: Warm-Up Phase

In this stage, we would like to prove the following theorem.

$\frac{\varepsilon^{2}}{c^{2}d}I\preceq A_{t}A_{t}^{\top}\preceq 2\Sigma$;

$\sigma_{d}(A_{T_{0}})\geq\sqrt{\frac{\sigma_{d}}{2}}$;

$\sigma_{1}(P_{T_{0}})\leq\frac{\sigma_{d}}{4}$;

$\|J_{t}\|_{op}\leq c\varepsilon\sqrt{\max\{m^{\prime},d\}}$, $\|K_{t}\|_{op}\leq c\varepsilon\sqrt{\max\{n^{\prime},d\}}$.

We first give some intuitions about the five conditions in Theorem 3.1. The first condition represents the “signal" is properly bounded from below and above throughout stage one. The second condition shows the magnitude of asymmetry is small throughout stage one. We note that it is crucial to study the Frobenius norm of $B$ instead of operator norm, because Frobenius norm admits a nice expansion for analysis. The third condition is an important one, which guarantees after $T_{0}$ iterations, we have enough “signal" strength in the principal space. The fourth condition is a technical one, which represents the symmetric error is small after $T_{0}$ iterations. The fifth condition represents the magnitude of the complement space remains small.

The proof of Theorem 3.1 is quite challenging and require new technical ideas and careful calculations, which we explain below.

If we only consider a differential equation $\dot{S}=(\Sigma-S)S+S(\Sigma-S)$, a well-known theorem (Theorem 12 in Lax (2007)) shows that if the singular values of $S$ are different from each other, and $\xi$ is the singular vector that $S\xi=\sigma_{d}(S)\xi$, then the derivative of $\sigma_{d}(S)$ is exactly $\xi^{\top}\dot{S}\xi$, which is lower bounded by $2(\sigma_{d}-\sigma_{d}(S))\sigma_{d}(S)$. To adapt it to discrete case, we prove the following lemma.

This lemma shows if we ignore perturbations from $B,J$ and $K$, then for small $\eta$ (when $\eta^{2}$ is of smaller order than the first term), the least eigenvalue of $S$ increases at a geometric rate.

However, there are also some small perturbation terms about $B,J$ and $K$ while doing analysis. $\|J\|_{op}$ and $\|K\|_{op}$ are easy to give an upper bound, since by (8) and (9), we know that by choosing small enough $\eta$, they are monotonically decreasing. However, the dynamic of $B$ is highly non-trivial. After some careful calculations (cf. (3.2.5)), we find that the increasing rate of $\|B\|_{F}^{2}$ is related to the smallest eigenvalue of $P:=\Sigma-AA^{\top}+BB^{\top}$: if $\max\{0,-\lambda_{d}(P)\}$ is small, then $\|B\|_{F}^{2}$ increases slowly.

Now we would like to give a lower bound on $\lambda_{d}(P)$. Inspired by gradient flow case, $P$ and $S$ are almost complementary of each other, and their dynamic behaves similarly. Hence we have $\dot{P}\approx-(\Sigma-P)P-P(\Sigma-P)$ with some small perturbation terms about $B,J$ and $K$. Hence we can use lemma 3.3 to give a lower bound in discrete case.

Notice that we use $B$ while analyzing $P$ and use $P$ while analyzing $B$. Hence, during the whole process, we need to bound both of them inductively.

Finally, once $\sigma_{d}(A)$ increases to a relatively large amount, we can use it to prove that $\|P\|_{op}$ will decrease exponentially fast to $\frac{\sigma_{d}}{4}$. One cannot simply prove that $P$ converges to zero in this stage, since the perturbation term $B$ will never converge to zero.

2.1 Assumptions

We make some assumptions on $A$ and $B$ in iterations $t\leq T_{0}$, where $T_{0}$ will be defined at the end of subsubsection 3.2.4, and we will verify the assumptions in the end.

$\frac{\varepsilon^{2}}{c^{2}d}I\preceq AA^{\top}\preceq 2\Sigma$.

The Frobenius norm of $B$ is bounded by $e_{b}d\varepsilon$ for some $e_{b}\geq c$, where $e_{b}$ will be determined laterWe will show later that it is appropriate to choose $e_{b}=2c$.. Hence its operator norm is also bounded by $e_{b}d\varepsilon$.

2.2 Dynamics on A, B and P

The dynamics on $J$ and $K$ is trivial, since by equations (8) and (9), i.e.

we know that if we choose $\eta\leq\frac{1}{3\sigma_{1}}$, one can inductively proved that $0\preceq V_{t}^{\top}V_{t}+K_{t}^{\top}K_{t}\preceq 3\sigma_{1}I$ and $0\preceq U_{t}^{\top}U_{t}+J_{t}^{\top}J_{t}\preceq 3\sigma_{1}I$ by using the first two assumptions in subsection 3.2.1. And then it follows that the operator norms of $J$ and $K$ are monotonically decreasing in this stage.

However, it is non-trivial to prove that $\|B\|_{op}$ keeps small. We will analyze the dynamics of $A,B$ and $P:=\Sigma-AA^{\top}+BB^{\top}$ together inductively.

First of all, from equations (6) and (7), we can write down the dynamics of $A:=\frac{U+V}{2}$ and $B:=\frac{U-V}{2}$ as following.

2.3 Dynamics on A

Given (13), we can give a lower bound for the minimal singular value of $A_{t+1}$.

For the first part, we could define $S_{t}:=A_{t}A_{t}^{\top}$, and $\overline{S}_{t+1}:=(I+\eta(\Sigma-S_{t}))S_{t}(I+\eta(\Sigma-S_{t}))$. Then according to lemma 3.2, by choosing $\beta=\frac{1}{2}$ and $\eta\leq\frac{1}{16\sigma_{1}}$, we havePlease see (23) for the full steps for this inequality.

For simplicity, we denote $\sigma_{d}(A_{t})$ by $a_{t}$, and define $s_{t}=\sigma_{d}(S_{t})=a_{t}^{2}$.

After some routine computationsPlease see section D for details., we can prove that it takes at most $T_{1}:=O\left(\frac{1}{\eta\sigma_{d}}\ln\frac{d\sigma_{d}}{\varepsilon^{2}}\right)$ iterations to make $a_{t}$ to at least $\sqrt{\frac{\sigma_{d}}{2}}$, and additional computations show that, if $a_{t}$ is always bounded by $\sqrt{2\sigma_{1}}$, then once $a_{t}$ becomes larger than $\sqrt{\frac{\sigma_{d}}{2}}$, it is always larger than $\sqrt{\frac{\sigma_{d}}{2}}$.

2.4 Dynamics on P

To bound $P_{t}$ by equation (15), we need to first bound the norms of $C_{t}$ and $D_{t}$ by (16) and (17).

By simple triangle inequalities we havePlease see (24) for full steps for this inequality.

where the last inequality holds when choosing $\varepsilon\leq\frac{\sqrt{\sigma_{1}}e_{b}d}{c^{2}(m+n)}$. Then we can conclude that

where $E_{t}$ is a matrix with operator norm less than $O(\eta^{2}\sigma_{1}^{3}+\eta e_{b}^{2}\varepsilon^{2}(m+n)d\sigma_{1}+\eta^{2}\sigma_{1}^{2}e_{b}^{2}d^{2}\varepsilon^{2})$. By choosing $\varepsilon\leq\frac{\sqrt{\sigma_{1}}}{e_{b}d}$ and $\eta\leq\frac{c^{2}\varepsilon^{2}(m+n)d}{\sigma_{1}^{2}}$, we have $\|E_{t}\|_{op}\leq O(\eta e_{b}^{2}\varepsilon^{2}(m+n)d\sigma_{1})$. Further more, by choosing $\beta=\frac{1}{2}$ and $\eta\leq\frac{1}{16\sigma_{1}}$ in lemma 3.3, we have

Because $P_{0}$ is initially positive, we know that

This lower bound verifies the assumption that $A_{t}A_{t}^{\top}\preceq 2\Sigma$, since $A_{t}A_{t}^{\top}=\Sigma-P+B_{t}B_{t}^{\top}\preceq\Sigma+O(e_{b}^{2}\varepsilon^{2}(m+n)d\kappa)I+e_{b}^{2}\varepsilon^{2}dI\preceq 2\Sigma$ by choosing $e_{b}^{2}\varepsilon^{2}=O\left(\frac{\sigma_{d}}{(m+n)d\kappa}\right)$.

On the other hand, we can also analyze the operator norm of $P$ by using formula (19), since $\sigma_{d}(A_{t})\geq\sqrt{\frac{\sigma_{d}}{2}}$ for $t\geq T_{1}$. This implies that

Inequality (21) shows that we only need at most $T_{2}:=O\left(\frac{1}{\eta\sigma_{d}}\ln\kappa\right)$ iterations after $T_{1}$ to make $\sigma_{1}(P_{t})\leq\frac{\sigma_{d}}{4}$. Because $\varepsilon^{2}\leq\sigma_{d}$, we have the total number of iteration $T_{0}:=T_{1}+T_{2}=O\left(\frac{1}{\eta\sigma_{d}}\ln\frac{d\sigma_{d}}{\varepsilon^{2}}\right)$.

2.5 Dynamics on B

To verify the assumption about $\|B\|_{F}$ made in subsection 3.2.1, we cannot simply use the equation (14), since the error term $\|(AB^{\top}-BA^{\top})A\|_{op}$ is approximately $O(\sigma_{1}\|B\|_{op})$, which will perturb the analysis seriously. Inspired by the continuous case that $\dot{\|B\|_{F}^{2}}=2\left\langle B,\dot{B}\right\rangle=\text{Tr}(B^{\top}PB)-\frac{1}{2}\|Q\|_{F}^{2}\leq\text{Tr}(B^{\top}PB)$, where $Q=AB^{\top}-BA^{\top}$ if we assume $J=K=0$. In this inequality, we hide the term $(AB^{\top}-BA^{\top})AB^{\top}$ in $-\frac{1}{2}\|Q\|_{F}^{2}$ and wipe it completely in our analysis.

Hence, for discrete case, we have the following inequalityPlease see (25) for full steps for this inequality.,

where the last equation is because we have chosen $\eta=O\left(\frac{\sigma_{d}\varepsilon^{2}}{d\sigma_{1}^{3}}\right)$ and $e_{b}^{2}\varepsilon^{2}=O\left(\frac{\sigma_{d}}{(m+n)d\kappa}\right)$.

3 Stage Two: Local Convergence Phase

We have proved in theorem 3.1 that the gradient descent achieved a pretty good point at $T_{0}$, i.e. $\|B_{T_{0}}\|_{F}\leq 2cd\varepsilon$ and $\sigma_{1}(P_{T_{0}})\leq\frac{\sigma_{d}}{4}$. In this subsection, we will prove that start from this point, the gradient descent will converge linearly to the global optimal point. Then theorem 1.1 follows.

We will prove inductively on the following conditions:

$\|B\|_{F}=O(\frac{\sigma_{d}}{\sqrt{\sigma_{1}}})$;

$\Delta_{t}:=\|\Sigma-U_{T_{0}+t}V_{T_{0}+t}^{\top}\|_{op}\leq\left(1-\frac{\eta\sigma_{d}}{2}\right)^{t}\frac{2}{5}\sigma_{d}$;

$\sigma_{d}(U),\sigma_{d}(V)\geq\sqrt{\frac{\sigma_{d}}{2}}$.

Intuitively, the (1) guarantees the magnitude of asymmetry remains small; (2) guarantees that in the principal space, the error converges to with a geometric rate; and (3) guarantees the “signal" in the principal space remains lower bounded.

First of all, it is easy to prove linear convergence of $J$ and $K$ by using assumption (3). Now we can verify the assumptions inductively.

$(1)+(2)\Rightarrow(3)$: We can prove $\sigma_{d}(UV^{\top})=\Theta(\sigma_{d})$. Because $U-V$ is small, (3) follows by triangle inequality.

$(1)+(3)\Rightarrow(2)$: Consider continuous-time case, if we assume $J=K=0$, the time derivative of $\Sigma-UV^{\top}$ is $-(\Sigma-UV^{\top})VV^{\top}-UU^{\top}(\Sigma-UV^{\top})$. Hence the convergence rate is lower bounded by $\sigma_{d}(U)$ and $\sigma_{d}(V)$. Because the perturbation term $J$ and $K$ decreases exponentially, assumption (2) follows naturally. We transform this intuition to the discrete-time case.

$(2)+(3)\Rightarrow(1)$: Again, we use (3.2.5) to show that the increasing rate of $\|B\|_{F}^{2}$ is bounded by $\Delta$. Because $\Delta$ decreases exponentially, $\|B\|_{F}^{2}$ cannot diverge to infinity, but increase by a $\text{poly}(m,n,\kappa)$ factor. Then by taking $\varepsilon$ sufficiently small can we verify the assumption (1).

To sum up, we have $\|\mathbf{\Sigma}-\mathbf{U}_{t}\mathbf{V}_{t}^{\top}\|_{F}^{2}=\|\Sigma-U_{t}V_{t}^{\top}\|_{F}^{2}+\|U_{t}K_{t}^{\top}\|_{F}^{2}+\|J_{t}V_{t}^{\top}\|_{F}^{2}+\|J_{t}K_{t}^{\top}\|_{F}^{2}$, which can be further bounded by

for some universal constant $C$. Hence one only needs $T_{f}:=O\left(\frac{\ln\frac{\sigma_{d}}{\delta}}{\eta\sigma_{d}}\right)$ iterations after $T_{0}$ to achieve an $\delta$-optimal point.

Conclusion

This paper proved that randomly initialized gradient descent converges to a global minimum of the asymmetric low-rank matrix factorization problem with a polynomial convergence rate. This result explains the empirical phenomena observed in prior work, and confirms that gradient descent with a constant learning rate still enjoys the auto-balancing property as argued in Du et al. (2018).

We believe our requirement of the step size $\eta$ is loose and a tighter analysis may improve the running time of gradient descent. Another interesting direction is to apply our techniques to other related problems such as asymmetric matrix sensing, asymmetric matrix completion and linear neural networks.

References

Appendix A Omitted Derivations of Formulas

We have omitted a number of complicated formulas in the main text to provide clear intuition and concise proof sketch. We will list all mentioned formulas here for readers’ reference.

Appendix B Dynamics in the Symmetric and Full-Rank Case

We consider the case where $U=V=A$ and $\Sigma$ is symmetric and full-rank, and we use gradient flow. We can derive the dynamics of $S=AA^{\top}$ as $\dot{S}:=(\Sigma-S)S+S(\Sigma-S)$, which is a quadratic ordinary differential equation and it is hard to solve directly.

However, if we define $\overline{X}:=S^{-1}$, we have $S\overline{X}\equiv I$. Taking the derivative implies $\dot{S}\overline{X}+S\dot{\overline{X}}=0$. Hence, $\dot{\overline{X}}=-S^{-1}\dot{S}S^{-1}$. Substitute $\dot{S}=(\Sigma-S)S+S(\Sigma-S)$ in it, we have

which is a linear ordinary differential equation.

For simplicity, define $X:=\overline{X}-\Sigma^{-1}$. Then

Similarly, because $P$’s dynamic is $\dot{P}=-(\Sigma-P)P-P(\Sigma-P)$, we have

And it is interesting to verify that $S(t)+P(t)\equiv\Sigma$ by using the following lemma.

Appendix C Proof of Lemmas

Since $\Sigma$ is invertible, we only need to verify the equation after right multiplying both side by $\Sigma^{-1}$. We have

where (31) is because $\Sigma$ commutes with $E$, (32) is because $\Sigma=S+P$ and finally (C) is because $\left(E(PS^{-1})E\right)^{-1}=E^{-1}(SP^{-1})E^{-1}$. ∎

First of all, we can expand the expression of $S^{\prime}$ and split it in the following terms.

For the first term $\beta S-2\eta S^{2}+\eta^{2}S^{3}$, its eigenvalues are $\beta s_{i}-2\eta s_{i}^{2}+\eta^{2}s_{i}^{3}$ since $S$ is commutable with itself, where $s_{i}$ is the $i^{\text{th}}$ largest singular value of $S$. By the assumptions $s_{i}\leq 2\sigma_{1}$ and $\eta\leq\frac{\beta}{8\sigma_{1}}$, we see the smallest eigenvalue of $\beta S-2\eta S^{2}+\eta^{2}S^{3}$ is exactly $\beta s-2\eta s^{2}+\eta^{2}s^{3}$.

For the second term, it can be rewritten as

Hence, the minimal singular value can be bounded by $\left(\sqrt{1-\beta}+\frac{\eta\sigma_{d}}{\sqrt{1-\beta}}\right)^{2}s$.

Finally, the last term can be lower bounded by $-\eta^{2}\sigma_{1}\left(-\frac{\beta}{1-\beta}\Sigma S\Sigma-\Sigma SS-SS\Sigma\right)\geq-\frac{8+6\beta}{1-\beta}\eta^{2}\sigma_{1}^{3}$. Summing up all three terms and we get

If $p\geq 0$, it suggests that $P$ is positive semi-definite, and $P^{\prime}$ is positive semi-definite, too. Hence $p^{\prime}\geq 0$ if $p\geq 0$.

If $p\leq 0$, we can expand the expression of $P^{\prime}$ and split it in the following terms.

For the first term $\beta P+2\eta P^{2}+\eta^{2}P^{3}$, its eigenvalues are $\beta p_{i}+2\eta p_{i}^{2}+\eta^{2}p_{i}^{3}$ since $P$ is commutable with itself, where $p_{i}$ is the $i^{\text{th}}$ largest eigenvalue of $P$. By the assumptions $|p_{i}|\leq 2\sigma_{1}$ and $\eta\leq\frac{\beta}{8\sigma_{1}}$, we see the smallest eigenvalue of $\beta P+2\eta P^{2}+\eta^{2}P^{3}$ is exactly $\beta p+2\eta p^{2}+\eta^{2}p^{3}$.

For the second term, it can be rewritten as

Hence, the minimal eigenvalue can be bounded by $\left(\sqrt{1-\beta}-\frac{\eta\sigma_{d}}{\sqrt{1-\beta}}\right)^{2}p$ if $p\leq 0$.

Finally, the last term can be lower bounded by $-\eta^{2}\sigma_{1}\left(-\frac{\beta}{1-\beta}\Sigma P\Sigma-\Sigma PP-PP\Sigma\right)\geq-\frac{8+6\beta}{1-\beta}\eta^{2}\sigma_{1}^{3}$. Summing up all three terms and we get that when $p\leq 0$,

Appendix D Solving the Iteration Formula of a

In this section we analyze the iteration formula (18).

We first consider the case when $a_{t}\leq\sqrt{\frac{\sigma_{d}}{2}}$. Notice that $a_{t}\geq\frac{\varepsilon}{c\sqrt{d}}$, we have

where we choose $\eta$ so small that $22\sigma_{1}^{3}\eta^{2}\leq\frac{\varepsilon^{2}}{c^{2}d}$.

By taking $\varepsilon=O\left(\frac{\sigma_{d}}{\sqrt{d^{3}\sigma_{1}}e_{b}^{2}(m+n)}\right)$ and $\eta=O\left(\frac{\sigma_{d}\varepsilon^{2}}{d\sigma_{1}^{3}}\right)$, we have $\frac{1}{2}\eta(\sigma_{d}-a_{t}^{2})a_{t}\geq\frac{1}{2}\eta\frac{\sigma_{d}}{2}\frac{\varepsilon}{c\sqrt{d}}\geq 22\frac{c\sqrt{d}}{\varepsilon}\sigma_{1}^{3}\eta^{2}+1.5\sqrt{2\sigma_{1}}\eta(e_{b}^{2}+c^{2})\varepsilon^{2}(m+n)d$, hence,

Subtracting $\sigma_{d}$ by (36), we have

Hence, $\frac{s_{T}}{\sigma_{d}-s_{T}}\geq(1+\eta\sigma_{d})^{T}\frac{s_{0}}{\sigma_{d}}$. So, it takes at most $T_{1}:=O\left(\frac{1}{\eta\sigma_{d}}\ln\frac{d\sigma_{d}}{\varepsilon^{2}}\right)$ iterations to bring $a_{t}$ to at least $\sqrt{\frac{\sigma_{d}}{2}}$.

Appendix E Solving the Iteration Formula on B

The iteration formula can be summarized as

where $p=O(\eta e_{b}^{2}\varepsilon^{2}(m+n)d\kappa)$ and $q=O(\eta\sqrt{\sigma_{1}}e_{b}(m+n)d^{2}\varepsilon^{3})$. Moreover, we have

Appendix F Proof of Stage Two

Here is the full version of the proof. Initially, $\|\Delta_{0}\|_{op}=\|P_{T_{0}}+Q_{T_{0}}\|_{op}$ where $Q=AB^{\top}-BA^{\top}$. Hence $\|\Delta_{0}\|_{op}\leq\sigma_{1}(P_{T_{0}})+\sigma_{1}(Q_{T_{0}})\leq\frac{\sigma_{d}}{4}+\sqrt{2\sigma_{1}}\sigma_{1}(B_{T_{0}})\leq\frac{\sigma_{d}}{3}$. Then for $U_{T_{0}}$ we have $\frac{2\sigma_{d}}{3}\leq\sigma_{d}(\Sigma)-\sigma_{1}(\Delta_{0})\leq\sigma_{d}(U_{T_{0}}V_{T_{0}}^{\top})\leq\sigma_{d}(U_{T_{0}}U_{T_{0}}^{\top})-2\sigma_{1}(U_{T_{0}}B_{T_{0}}^{\top})\leq\sigma_{d}(U_{T_{0}}U_{T_{0}}^{\top})-4\sqrt{2\sigma_{1}}O\left(\frac{\sigma_{d}}{\sqrt{\sigma_{1}}}\right)$. Hence $\sigma_{d}(U_{T_{0}})\geq\sqrt{\frac{\sigma_{d}}{2}}$. We can do the same thing on $V_{T_{0}}$.

First of all, by equations (8) and (9), we have

Expanding $\Sigma-U_{t+1+T_{0}}V_{t+1+T_{0}}^{\top}$ by brute forcePlease see (26) for the result of the expanding., we get

Thus we can now verify that $\Delta_{t}\leq\left(1-\frac{\eta\sigma_{d}}{2}\right)^{t}\frac{2}{5}\sigma_{d}$. Together with the linear convergence of $J$ and $K$, we know the gradient descent converge linearly. Notice that by using the operator norm of $\Delta_{t}$, we can easily prove that $\sigma_{d}(U)$ and $\sigma_{d}(V)$ in the next iteration is at least $\sqrt{\frac{\sigma_{d}}{2}}$ once given $\|B_{T_{0}+t}\|_{F}$ is small.

To give an upper bound on $\|B\|_{F}$, we still use equation (3.2.5).

First of all, we have $\|P\|_{F}^{2}+\|Q\|_{F}^{2}=\|\Sigma-UV^{\top}\|_{F}^{2}$, since $P+Q=\Sigma-UV^{\top}$, and $\left\langle P,Q\right\rangle=0$. Hence, $\|P_{t+T_{0}}\|_{F}\leq\sqrt{d}\Delta_{t}$ and $\|Q_{t+T_{0}}\|_{F}\leq\sqrt{d}\Delta_{t}$.

To solve this iteration formula, we first notice that the product of the main coefficient is bounded by a universal constant,

we can then write it into an iteration formula about $\frac{\|B_{t+T_{0}}\|_{F}^{2}}{\Xi_{t}}$,

By taking $\varepsilon=O\left(\frac{\sigma_{d}}{\sqrt{\sigma_{1}(m+n)d}}\right)$ and $\eta=O\left(\frac{\sigma_{d}}{d\sigma_{1}^{2}}\right)$, induction on $\|B\|_{F}$ holds.