Guaranteed Matrix Completion via Non-convex Factorization

Introduction

There are two popular approaches to impose the low-rank structure: the nuclear norm based approach and the matrix factorization (MF) based approach. In the first approach, the whole matrix is the optimization variable and the nuclear norm (denoted as $\|\cdot\|_{*}$) of this matrix variable, which can be viewed as a convex approximation of its rank, serves as the objective function or a regularization term. For the matrix completion problem, the nuclear norm based formulation becomes either a linearly constrained minimization problem

a quadratically constrained minimization problem

On the theoretical side, it has been shown that given a rank-$r$ matrix $M$ satisfying an incoherence condition, solving (1) will exactly reconstruct $M$ with high probability provided that $O(r(m+n)\log^{2}(m+n))$ entries are uniformly randomly revealed . This result was later generalized to noisy matrix completion, whereby the optimization formulation (2) is adopted . Using a different proof framework, reference provided theoretical guarantee for a variant of the formulation (3). On the computational side, problems (1) and (2) can be reformulated as a semidefinite program (SDP) and solved to global optima by standard SDP solvers when the matrix dimension is smaller than 500. To solve problems with larger size, researchers have developed first order algorithms, including the SVT (singular value thresholding) algorithm for the formulation (1) , and several variants of the proximal gradient method for the formulation (3) . Although linear convergence of the proximal gradient method has been established for the formulation (3) under certain conditions , the per-iteration cost of computing SVD (Singular Value Decomposition) may increase rapidly as the dimension of the problem increases, making these algorithms rather slow or even useless for problems of huge size. The other major drawback is the memory requirement of storing a large $m$ by $n$ matrix.

A popular factorization based formulation for matrix completion takes the form of an unconstrained regularized square-loss minimization problem :

Despite the great empirical success, the theoretical understanding of the algorithms for the factorization based formulation is fairly limited. More specifically, the fundamental question of whether these algorithms (including many recently proposed ones) can recover the true low-rank matrix remains largely open. In this paper, we partially answer this question by showing that under similar conditions to those used in previous works, many standard optimization algorithms for a factorization based formulation (see (18)) indeed converge to the true low-rank matrix (see Theorem 3.1). Our result applies to a large class of algorithms including gradient descent, SGD and many block coordinate descent type methods such as two-block alternating minimization and block coordinate gradient descent. We also show the linear convergence of some of these algorithms (see Theorem 3.2 and Corollary 3.2).

To the best of our knowledge, our result is the first one that analyzes the geometry of matrix factorization in Euclidean space for matrix completion. In addition, our result also provides the first recovery guarantee for alternating minimization without resampling (i.e. without using independent samples in different iterations). Below we elaborate these two contributions in light of the existing works.

1) We analyze the local geometry of the matrix factorization formation (in Euclidean space). We argue that the success of many algorithms attributes mostly (or at least partially) to the geometry of the problem, rather than the specific algorithms being used. The geometrical property we establish is that the local gradient direction $-\nabla f(x)$ is aligned with the global descent direction $x^{*}-x$. For the classical matrix factorization formulation $\|M-XY^{T}\|_{F}^{2}$, we develop a novel perturbation analysis to deal with the ambiguity of the factorization. For the sampling loss $\|\mathcal{P}_{\Omega}(M-XY^{T})\|_{F}^{2}$, an incoherence regularizer (or constraint) is needed, which causes an extra difficulty of analyzing nonconvex constrained optimization. Unfortunately, projection to the constraint (or the gradient of the regularizer) is not aligned with the global direction, and we add one more regularizer to “correct” the local descent direction. A high-level lesson is that regularization may change the geometry of the problem.

2) Our result applies to the standard forms of the algorithms (though our optimization formulation is a bit different), which do not require the additional resampling scheme used in other works . We obtain a sample complexity bound that is independent of the recovery error $\epsilon$, while all previous sample complexity bounds for the matrix factorization based formulation (in Euclidean space) depend on $\epsilon$. There is a subtle theoretical issue for the resampling scheme; see more discussions in Section 1.2 and [30, Sec. 1.5.3].

2 Related works

Factorization models. The first recovery guarantee for the factorization based matrix completion is provided in , where Keshavan, Montanari and Oh considered a factorization model in Grassmannian manifold and showed that the matrix can be recovered by a proper initialization and a gradient descent method on Grassmannian manifold. Besides being quite complicated, this model is not as flexible as the factorization model in Euclidean space, and it is not easy to solve by many advanced large-scale optimization algorithms. Moreover, most algorithms in Grassmann manifold require line search, and little is known about the convergence rate.

The factorization model in Euclidean space was first analyzed in an unpublished work of Keshavan Reference is a PhD thesis that discusses various algorithms including the algorithm proposed in and alternating minimization. In this paper when we refer to , we are only referring to [17, Ch. 5] which presents resampling-based alternating minimization and the corresponding result., as well as a later work of Jain et al. . Both works considered alternating minimization with resampling scheme, a special variant of the original alternating minimization. The sample complexity bounds were later improved by Hardt and Hardt and Wooters , where in the latter work, notably, the authors devised an algorithm with a corresponding sample complexity bound independent of the condition number. However, these improvements are obtained for more sophisticated versions of resampling-based alternating minimization, not the typical alternating minimization algorithm.

Resampling. The issues of resampling have been discussed in a recent work on phase retrieval by Candès et al.. We will point out a subtle theoretical issue not mentioned in , as well as some other practical issues.

The resampling scheme (a.k.a. golfing scheme ) can be used at almost no cost for the nuclear norm approach , but for the alternating minimization it causes many issues. At first, it may seem that for both approaches resampling is a cheap way to get around a common difficulty: the dependency of the iterates on the sample set. However, there is a crucial difference: for the nuclear norm approach, resampling is just a proof technique used in a “conceptual” algorithm for constructing the dual certificate, while for the alternating minimization, resampling is used in the actual algorithm. This difference causes some issues of resampling-based alternating minimization at conceptual, practical and theoretical levels.

1) Gap between theory and algorithm. Algorithmically, an easy resampling scheme is to randomly partition the given set $\Omega$ into non-overlapping subsets $\Omega_{k},k=1,\dots,L$, as proposed in The description in has some ambiguity and it might refer to the scheme of sampling $\Omega_{k}$’s with replacement; anyhow, under this model $\Omega_{k}$’s are still dependent. See [30, Sec. 1.5.3] for more discussions. . However, the results in actually require a generative model of independent $\Omega_{k}$’s, instead of sampling $\Omega_{k}$’s based on a given $\Omega$. Therefore, the results in do not directly apply to the partition based resampling scheme that is easy to use. See [30, Sec. 1.5.3] for more discussions on this subtle issue.

This issue has been discussed by Hardt and Wooters in [20, Appendix D], and they proposed a new resampling scheme [20, Algorithm 6] to which the results in can apply, provided that the generative model of $\Omega$ is exactly known. In practice, the underlying generative model of $\Omega$ is usually unknown, in which case the scheme [20, Algorithm 6] does not work. In contrast, the classical results in and our result herein are robust to the generative model of $\Omega$: these results actually state that for an overwhelming portion of $\Omega$ with a given size, one can recover $M$ through a certain algorithm, thus for many reasonable probability distributions of $\Omega$ a high probability result holds.

2) Impracticality. As argued previously, assuming a generative model of $\Omega_{k}$’s is not practical since $\Omega$ is usually given. For given $\Omega$, the only known validated resampling scheme [20, Algorithm 6], besides not being robust to the underlying generative model of $\Omega$, might be a bit complicated to use in practice. Even the simple resampling scheme of partitioning $\Omega$ (which has not been validated yet) is rather unrealistic since each sample is used only once during the algorithm.

3) Inexact recovery. A theoretical consequence of the resampling scheme is that the required sample complexity $|\Omega|$ becomes dependent on the desired accuracy $\epsilon$, and goes to infinity as $\epsilon$ goes to zero. This is different from the classical results (and ours) where exact reconstruction only requires finite samples. While it is common to see the dependency of time complexity on the accuracy $\epsilon$, it is relatively uncommon to see the dependency of sample complexity on $\epsilon$.

In a recent work the authors have managed to remove the dependency of the required sample size on $\epsilon$ by using a singular value projection algorithm. However, considers a matrix variable of the same size as the original matrix, which requires significantly more memory than the matrix factorization approach considered in this paper. Moreover, it requires resampling at a number of iterations (though not all), which may suffer from the same issues we mentioned earlier. The resampling is also required in the recent work of ; see [30, Sec. 1.5.3] for more discussions.

Other works on non-convex formulations. Non-convex formulation has also been studied for the phase retrieval problem in some recent works . These works provide theoretical guarantee for some algorithms specially tailored to certain non-convex formulations and with specific initializations. The major difference between and is that the former requires independent samples in each iteration, while the latter uses the same samples throughout in the proposed algorithm. As mentioned earlier, such a difference also exists between all previous works on alternating minimization for matrix completion and our work.

Finally, we note that there is a growing list of works on the theoretical guarantee of non-convex formulations for various problems, such as sparse regression (e.g. ), sparse PCA , robust PCA and EM (Expected-Maximization) algorithm . We emphasize several aspects that distinguish our paper from other recent works on non-convex optimization. First, our paper is one of the first to analyze the (local) geometry of the problem. Second, we deal with non-symmetric matrix factorization which has a more bizarre geometry than symmetric matrix factorization and some other models. Third, one difficulty of our problem essentially lies in nonconvex constrained optimization (though we consider the closely related regularized form).

3 Proof Overview and Techniques

Basic idea: local geometry. The very first question is what kind of property can ensure global convergence for non-convex optimization. We will establish a local geometrical property of a regularized objective such that any stationary point in a local region is globally optimal. This is achieved in three steps: (i) study the local geometry of the fully observed objective $\|M-XY^{T}\|_{F}^{2}$; (ii) study the local geometry of the matrix completion objective $\|\mathcal{P}_{\Omega}(M-XY^{T})\|_{F}^{2}$; (iii) study the local geometry of a regularized objective. Next, we will discuss the difficulties involved in each step and describe how we address these difficulties.

Local geometry of $\|M-XY^{T}\|_{F}^{2}$. We start by considering a simple case that $M$ is fully observed and the objective function is $f(X,Y)=\|M-XY^{T}\|_{F}^{2}$. What is the geometrical landscape of this function? In the simplest case $m=n=r=1$ and $f(x,y)=(xy-1)^{2}$, the set of stationary points is $\{(x,y)\mid xy=1\}\cup\{(0,0)\}$, in which $(0,0)$ is a saddle point and the curve $xy=1$ consists of global optima. We plot the function around the curve $xy=1$ in the positive orthant in Figure 1.

An interpretation is that the negative gradient direction $-\nabla f$ should be aligned with the global direction $(u,v)-(x,y)$; a convex function has a similar property, but the difference is that here the global direction is adjusted according to the position of $(x,y)$.

For general $m,n,r$, the geometrical landscape is probably much more complicated than the scalar case. Nevertheless, we can still prove that the convexity of $\|M-Z\|^{2}$ is partially preserved when reparameterizing $Z$ as $Z=XY$. The exact expression is a variant of (6) which we will discuss in more detail later. Technically, we need to connect the Euclidean space and the quotient manifold via “coupled perturbation analysis”: given $X,Y$ such that $\|XY^{T}-M\|_{F}$ is small, find decomposition $M=UV^{T}$ such that $U,V$ are close to $X$ and $Y$ respectively (a simpler version of Proposition 4.1). The difference from traditional perturbation analysis of Wedin (i.e. if two matrices are close then their row/column spaces are close) is that in the row/column spaces are fixed while in our problem $U,V$ are up to our choice.

Local geometry of $\|\mathcal{P}_{\Omega}(M-XY^{T})\|_{F}^{2}$. Let us come back to the original matrix completion problem, in which an additional sampling operator $\mathcal{P}_{\Omega}$ is introduced. Similarly, we hope that $f_{\Omega}(Z)=\frac{1}{2}\|\mathcal{P}_{\Omega}(M-Z)\|^{2}$ is strongly convex and this strong convexity can be partially preserved after reparametrization $Z=XY^{T}$. However, one issue is that the function $f_{\Omega}(Z)$ is possibly non-strongly-convex (though still convex). In fact, if $f_{\Omega}$ is locally strongly convex around $M$, then we should have

Assuming $Z$ is rank-$r$, this inequality can be rewritten as

where $K(\delta)$ is a neighborhood of $M$ defined as $\{(X,Y)\mid\|XY^{T}-M\|_{F}\leq\delta\}$ and $C$ is a numerical constant. We wish (5) to hold with high probability (w.h.p.) for random $\Omega$ in which each position in $M$ is chosen with probability $p$. This inequality is closely related to matrix RIP (restricted isometry property) in (see equation (III.4) therein). If $X,Y$ are independent of $\Omega$, then (5) follows easily from the concentration inequalities. Unfortunately, if $X,Y$ are chosen arbitrarily instead of independently from $\Omega$, the bound (5) may fail to hold.

A solution, as employed in , is to utilize a random graph lemma in which provides a bound on $\|\mathcal{P}_{\Omega}(A)\|_{F}$ for any rank-$1$ matrix $A$ (possibly dependent on $\Omega$). This lemma, combined with another probability result in , implies a bound on $\|\mathcal{P}_{\Omega}(M-XY^{T})\|_{F}$. However, this bound is not good enough since it only leads to (5) when $\delta=O(1/n)$. The underlying reason is that the bound given by the random graph lemma is actually quite loose if $X$ or $Y$ have unbalanced rows, i.e. certain row has large norm. One solution is to force the iterates to have bounded row norms (a.k.a. incoherent), by adding a constraint or regularizer. With the incoherence requirement on $X,Y$, now (5) can be shown to be hold for $\delta=O(1)$, or more precisely, $\delta=O(\Sigma_{\min})$, where $\Sigma_{\min}$ is the minimum eigenvalue of $M$. With such a $\delta$, it is possible to find an initial point in the region $K(\delta)$.

In summary, although $f_{\Omega}(Z)=\frac{1}{2}\|\mathcal{P}_{\Omega}(Z-M)\|_{F}^{2}$ is possibly non-strongly-convex, by restricting to an incoherent neighborhood of $M$ it is “relative” strongly convex (called “relative” since we fix $M$ in (5)). More specifically, we have that w.h.p.

where $K_{1}$ denotes the set of $(X,Y)$ with bounded row norms. Note that this inequality also implies that global optimally in $\mathcal{B}$ leads to exact recovery; or equivalently, zero training error leads to zero generalization error.

Having established the geometry of $f_{\Omega}(Z)$, we can use the same technique for the fully observed case to show the local geometry For illustration purpose, we present a two-step approach: first establish a geometrical property of $f_{\Omega}(Z)$, then extend the property to $f_{\Omega}(XY^{T})$. However, our current proof does not follow the two-step approach but directly establish the property of $f_{\Omega}(XY^{T})$. In fact, although we establish the property of $f_{\Omega}(Z)$ in Claim 3.1, the proof of this claim is very similar to the proof of (7). of

Denoting $\bm{x}=(X,Y),\bm{x}^{*}=(U,V)$ and utilizing $\nabla F(\bm{x}^{*})=0$, (7) becomes

It links the local optimality measure $\|\nabla F(\bm{x})\|$ with the global optimality measure $\text{dist}(\bm{x},\mathcal{X}^{*})=\min_{\bm{x}^{*}\in\mathcal{X}^{*}}\|\bm{x}-\bm{x}^{*}\|$, and implies that any stationary point of $F$ in $\mathcal{B}$ is a global minimum.

If (8) holds for arbitrary $\bm{x},\bm{x}^{*}$ then $F$ would be strongly convex in $\bm{x}$. Let us emphasize again two differences of (8) with local strong convexity: i) since $\bm{x}^{*}$ is not arbitrary but has to be one global minimum, (8) indicates local “relative convexity” of $F$; ii) due to the ambiguity of factorization, $\bm{x}^{*}$ should be chosen according to $\bm{x}$, thus (8) indicates local relative convexity up to a group transformation (it might be conceptually helpful to view it as a property in the quotient manifold, but we do not explicitly exploit its structure).

Local geometry with regularizers/constraints. The property (8) is still not desirable. The original purpose of studying geometry is to show there is no spurious “1st order local-min” (point that satisfies 1st order optimality conditions). To establish the geometrical property with sampling, we restrict to an incoherent set $K_{1}$, but this restriction changes the meaning of the 1st order local-min. In fact, to ensure the iterates stay in the incoherent region $K_{1}$, we need to solve a constrained optimization problem $\min_{\bm{x}\in K_{1}}F(\bm{x})$ or a regularized problem $\min_{\bm{x}}F(\bm{x})+G_{1}(\bm{x})$ where $G_{1}$ is a regularizer forcing $\bm{x}$ to be in $K_{1}$. Standard optimization algorithms converge to the KKT points of $\min_{\bm{x}\in K_{1}}F(\bm{x})$ or the stationary points of $F+G_{1}$, which may not be the stationary points of $F$. The property (8) only implies any stationary point of $F$ in $\mathcal{B}$ is globally optimal.

We shall focus on the regularized problem $\min F+G_{1}$; the constrained problem $\min_{\bm{x}\in K_{1}}F$ is similar. Because of the extra regularizer, the property (8) is not enough. We need to prove a result similar to (8), but with $\nabla F$ replaced by $\nabla F+\nabla G_{1}$:

then combining with the existing result (8) we are done; unfortunately, we do not know how to prove (10). Intuitively, (10) means that $-\nabla G_{1}(\bm{x})$, which is almost the same direction as the projection to the incoherent region $K_{1}$, is positively correlated with the global direction $\bm{x}^{*}-\bm{x}$. At first sight, this seems trivially true because for any point $\bar{\bm{x}}\in K_{1}$ we have $\langle\nabla G_{1}(\bm{x}),\bm{x}-\bar{\bm{x}}\rangle\geq 0$ (as illutrated in Fig. 2). However, a rather strange issue is that $\bm{x}^{*}$ is chosen to be a point in $\{(U,V)\mid UV^{T}=M\}$ that is close to $\bm{x}$, thus there is no guarantee that $\bm{x}^{*}$ lies in $K_{1}$. An underlying reason is that the global optimum set $\{(U,V)\mid UV^{T}=M\}$ is unbounded and thus not a subset of $K_{1}$. If we enforce $(U,V)$ to be in $K_{1}$, we may not be able to find $(U,V)$ that is close enough to $(X,Y)$.

Technically, the issue is that $(U,V)$ chosen in Proposition 4.1 have row-norms bounded above by quantities proportional to the norms of $X,Y$, and can be higher than the row-norms of $X,Y$ (threshold of $K_{1}$). To resolve this issue, we add an extra regularizer $G_{2}(X,Y)$ to force $(X,Y)$ to lie in $K_{2}$, a set of matrix pairs with bounded norms. This extra bound makes $\langle\nabla G_{1}(\bm{x}),\bm{x}-\bm{x}^{*}\rangle\geq 0$ straightforward to prove, but a similar issue arises: now we need to prove (8) for $F+G_{1}+G_{2}$ instead of $F$. Again, it suffices to prove that for any $\bm{x}\in K(\delta)\cap K_{1}\cap K_{2}$ there exists $\bm{x}^{*}$ such that

Constrained perturbation analysis. The desired inequality (11) is implied by the following condition on $U,V$: $\|U\|_{F}\leq\|X\|_{F},\|V\|_{F}\leq\|Y\|_{F}$ when $\|X\|_{F},\|Y\|_{F}$ are large. Recall that previously we try to find $U,V$ that are close to $X,Y$; see Proposition 4.1. Now we need to impose extra constraints on $U,V$, giving rise to Proposition 4.2. The extra constraints make the perturbation analysis significantly more involved; in fact, we apply a sophisticated iterative procedure to construct the factorization $M=UV^{T}$. The main steps of the proof are briefly given in Appendix C.2.

One crucial component of our proof can be viewed as the perturbation analysis for “preconditioning”. Roughly speaking, the basic problem is: given an $r\times r$ matrix $\hat{X}$ with a large condition number, find another matrix $\hat{U}$ with the same Frobenius norm as $\hat{X}$ but smaller inverse Frobenious norm (i.e. $\|\hat{U}^{-1}\|_{F}\leq\frac{1}{1-\delta}\|\hat{X}^{-1}\|_{F}$). In other words, we want to reduce $\sum_{i=1}^{r}\frac{1}{\sigma_{i}^{2}}$ with $\sum_{i=1}^{r}\sigma_{i}^{2}$ fixed, where $\sigma_{i}$’s are all singular values. Intuitively, by reducing $\sum_{i=1}^{r}\frac{1}{\sigma_{i}^{2}}$ we reduce the discrepancy of singular values. This process is somewhat similar to preconditioning in numerical algebra that reduces the gap between the largest and smallest eigenvalue. The precise statement of the basic problem and its relation with the key technical result Proposition 4.2 are provided in Appendix C.2.1.

Algorithm requirements. We provide three conditions and show that if an algorithm satisfies either of them, then with specific initialization the iterates will stay in the desired basin (see Proposition 5.1). A special case of the third condition has been used in for Grassmann manifold optimization. Together, these three conditions cover a wide spectrum of algorithms including GD, SGD and block coordinate descent type methods.

Proof outline. The overall proof can be divided into two parts: the geometrical property (Lemma 3.1) and the algorithm property (Lemma 3.2). For the geometrical property, Lemma 3.1 states that the regularized objective function $F+G_{1}+G_{2}$ enjoys some nice geometrical property in a certain local region around the global optima, thus there is no other stationary point in this region. For the algorithm property, Lemma 3.2 states that starting from an easily computable initial point, many standard algorithms generate a sequence that are inside the desired region and these algorithms also converge to stationary points. Since these stationary points must be global optima by Lemma 3.1, we obtain that these algorithms converge to the global optima.

4 Other Remarks

Difference with previous works. As discussed earlier, one major challenge is to bound $\mathcal{P}_{\Omega}(A)$ when $A$ may be dependent on $\Omega$. One simple strategy as adopted in is to use a resampling scheme to decouple $A$ and the observation set. This strategy artificially avoids this difficulty, and causes a few issues discussed earlier in Section 1.2. Another strategy, as employed in , is to use a random graph lemma in .

We apply the random graph lemma of when extending the local geometry of $\|M-XY^{T}\|_{F}^{2}$ to $\|\mathcal{P}_{\Omega}(M-XY^{T})\|_{F}^{2}$. The difference of our work with is that we study the local geometry in Euclidean space (and, indirectly, the geometry of the quotient manifold), which is quite different from the local geometry in Grassmann manifold studied in . Technically, the complications of the proof in are mostly due to heavy computation of various quantities in Grassmann manifold; in addition, much effort is spent in estimating the terms related to the extra factor $S$ which enables the decoupling of $X$ and $Y$ ( actually uses a three-factor decomposition $XSY^{T}$). For our problem, one difficulty is to “pull back” the distance in the quotient manifold to the Euclidean space, by the coupled perturbation analysis. Another difficulty is to align the gradient of the regularizer with the global direction (this is not an issue for Grassman manifold), which requires a more sophisticated perturbation analysis. The difficulties have been discussed in detail in Section 1.3.

Symmetric PSD or rank-1 case. The symmetric PSD (positive semi-definite) case or the rank-1 case are easier to deal with, because in the 3-step study of the local geometry the third step is not necessary. When $M$ is rank-1 (possibly non-symmetric), the regularizer $G_{2}(\cdot)$ may still be needed, but Proposition 4.2 is trivial since its assumptions cannot hold for $r=1$. When $M$ is symmetric PSD, a popular approach is to use a symmetric factorization $M=XX^{T}$ instead of the non-symmetric factorization, and the loss function becomes $\|\mathcal{P}_{\Omega}(M-XX^{T})\|_{F}^{2}$. The same proof in our paper can be translated to this symmetric PSD case, except that the third step is not necessary. In fact, it is possible to show that (10) holds without any additional requirement on $\bm{x}$. As a result, the regularizer $G_{2}$ and a major technical result Proposition 4.2 are not needed. In both the symmetric PSD and rank-1 case, we only need to establish the intermediate result (7) and the proof can be greatly simplified. Stronger sample complexity and time complexity bounds may be established in these two cases.

Simulation Results The regularizers are introduced due to theoretical purposes; interestingly, they turn out to be helpful in the numerical experiments (the comments below are extracted from the thesis [30, Chapter 2]).

First, the simulation suggests that the imbalance of the rows of $X$ or $Y$ is an important issue for matrix completion in practice, a phenomenon not reported before to our knowledge. The table in Figure 2.10 of shows that when $|\Omega|$ is small, in all successful instances the iterates are balanced, while in all failed instances the iterates are unbalanced. This contrast occurs for many standard algorithms such as AltMin,GD and SGD.

Second, adding only the regularizer $G_{1}$ helps, but not too much. Adding an extra regularizer $G_{2}$ can push the sample complexity to be very close to the fundamental limit, at least for the synthetic Gaussian data. These experiments seem to indicate that the new regularizers do change the geometry of the problem.

Necessity of incoherence? While our regularizers are helpful when $|\Omega|$ is small, an open question is whether the row-norm requirement is needed for the local geometry when $|\Omega|$ is large. We observe that the row-norms can be automatically controlled by standard algorithms for the synthetic Gaussian data when there are, say, $5rn$ samples for $n\times n$ matrices. There are two possible explanations (assuming a large $|\Omega|$): (i) the local geometrical property (7) holds without the incoherence requirement; (ii) (7) still requires incoherence, but there is an unknown mechanism for many algorithms to control the row-norms.

To exclude the first possibility, we need to find $(X,Y)\in K(\delta)$ such that $\nabla F(X,Y)=0$ but $XY^{T}\neq M$; since (7) holds, such $(X,Y)$ must have unbalanced row-norms. Such an example would validate the necessity of the incoherence restriction for the local geometry. Note that the necessity of incoherence for the local geometry is different from the necessity of an incoherence regularizer/constraint for a specific algorithm. Even if the local geometry requires incoherence, it remains an interesting question why many algorithms can automatically control row-norms when $|\Omega|$ is large.

5 Notations and organization

Organization. The rest of the paper is organized as follows. In Section 2 we introduce the problem formulation and four typical algorithms. In Section 3, we present the main results and the main lemmas used in the proofs of these results. The proof of the two lemmas used in proving Theorem 3.1 are given in Section 4 and Section 5 respectively. The proof of the first lemma depends on two “coupled perturbation analysis” results Proposition 4.1 and Proposition 4.2, the proofs of which are given in Appendix B and Appendix C respectively. The proof of a lemma used in proving Theorem 3.2 is given in Appendix E.

Problem Formulation and Algorithms

Incoherence condition. The incoherence condition for the matrix completion problem is first introduced by Candès and Recht in and has become a standard assumption for low-rank matrix recovery problems (except a few recent works such as ). We will define an incoherence condition for an $m\times n$ matrix $M$ which is the same as that in .

We say a matrix $M=\hat{U}\Sigma\hat{V}^{T}$ (compact SVD of $M$) is $\mu$-incoherent if:

It can be shown that $\mu\in[1,\frac{\max\{m,n\}}{r}]$. For some popular random models for generating $M$, the incoherence condition holds with a parameter scaling as $\sqrt{r\log n}$ (see ). In this paper, we just assume that $M$ is $\mu$-incoherent. Note that the incoherence condition implies that $\hat{U},\hat{V}$ have bounded row norm. Throughout the paper, we also use the terminology “incoherent” to (imprecisely) describe $m\times r$ or $n\times r$ matrices that have bounded row norm (see the definition of set $K_{1}$ in (30)).

Random sampling model. In the statement of the results in this paper, the probability is taken with respect to the uniform random model of $\Omega\subseteq[m]\times[n]$ with fixed size $|\Omega|=S$ (i.e. $\Omega$ is generated uniformly at random from set $\{\Omega^{\prime}\subseteq[m]\times[n]:\text{ the size of }\Omega^{\prime}\text{ is }S\}$ ). We remark that this model is “equivalent to” a Bernolli model that each entry of $M$ is included into $\Omega$ independently with probability $p=\frac{S}{mn}$ in the sense that if the success of an algorithm holds for the Bernolli model with a certain $p$ with high probability, then the success also holds for the uniform random model with $|\Omega|=pmn$ with high probability (see or [31, Sec. 1D] for more details). Thus in the proofs we will instead use the Bernolli model.

2 Problem formulation

We consider a variant of (P0) with incoherence-control regularizers. In particular, we introduce two types of regularization terms besides the square loss function: the first type is designed to force the iterates $X_{k},Y_{k}$ to be incoherent (i.e. with bounded row norm), and the second type is designed to upper bound the norm of $X_{k}$ and $Y_{k}$. Note that (P0) is related to the Lagrangian method, while our regularizer is based on the penalty function method for constrained optimization problems. We can also view the regularizer $\lambda(\|X\|_{F}^{2}+\|Y\|_{F}^{2})$ as a “soft regularizer”, and our new regularizer as a “hard regularizer”. The advantage of the hard regularizer is that it does not distort the optimal solution.

Our regularizers are smooth functions with simple gradients, thus the algorithms for our formulation have similar per-iteration computation cost as the algorithms for the formulation without regularizers. In the numerical experiments, we find that when $|\Omega|$ is large, the iterates are always incoherent and bounded, and our algorithms are the same as the traditional algorithms for the unregularized formulation; when $|\Omega|$ is relatively small, the traditional algorithms may produce high error, and our regularizer becomes active and significantly reduce the error. In some sense, our algorithms for the new formulation are “better” versions of the traditional algorithms, and our theoretical results can be viewed as a validation of the traditional algorithms in the “large-$|\Omega|$ regime” and a validation of the modified algorithm in the “small-$\Omega$” regime. Preliminary simulation results show that many algorithms for the proposed formulation can recover the matrix when $|\Omega|$ is very close to the fundamental limit, significantly improving upon the traditional algorithms; see [30, Chapter 3].

The regularization function $G$ is defined as follows:

where $A^{(i)}$ denotes the $i$th row of a matrix A,

Here, $I_{\mathcal{C}}$ is the indicator function of a set $\mathcal{C}$, i.e. $I_{\mathcal{C}}(z)$ equals $1$ when $z\in\mathcal{C}$ and otherwise. $\rho$ is a constant specified shortly. Throughout the paper, $\delta$ and $\delta_{0}$ are defined as

where $C_{d}$ is some numerical constant. The coefficient $\rho$ is defined as (a larger $\rho$ also works)

The numerical constant $C_{T}>5$ will be specified in the proof of our main result. The parameter $\beta_{T}$ is chosen to be of the same order as $\|\hat{U}\Sigma^{1/2}\|_{F}$ and $\|\hat{V}\Sigma^{1/2}\|_{F}$, and $\beta_{1},\beta_{2}$ are chosen to be of the same order as $\sqrt{r}\|(\hat{U}\Sigma^{1/2})^{(i)}\|,\sqrt{r}\|(\hat{V}\Sigma^{1/2})^{(j)}\|$. The additional factor $\sqrt{3r}$ is due to technical consideration (to prove (256)). Our regularizer $G$ involves $\Sigma_{\max}$ and $\mu$ which depend on the unknown matrix $M$; in practice, we can estimate $\Sigma_{\max}$ by $c_{1}\sqrt{\frac{\|\mathcal{P}_{\Omega}(M)\|_{F}^{2}}{pr}}$, and estimate $\mu$ by $c_{2}\frac{\sqrt{mn}}{r\Sigma_{\max}}\max_{(i,j)\in\Omega}|M_{ij}|$ (according to (203)) where $c_{1},c_{2}$ are numerical constants to tune.

It is easy to verify that $G_{0}$ is continuously differentiable. The choice of function $G_{0}$ is not unique; in fact, we can choose any $G_{0}$ that satisfies the following requirements: a) $G_{0}$ is convex and continuously differentiable; b) $G_{0}(z)=0,z\in$. In , $G_{0}$ is chosen as $G_{0}(z)=I_{[1,\infty]}(z)(e^{(z-1)^{2}}-1)$, which also satisfies these two requirements. Choosing different $G_{0}$ does not affect the proof except the change of numerical constants (which depend on $G_{0}(3/2),G_{0}^{\prime}(3/2),G_{0}^{\prime\prime}(3/2)$). Note that the requirement of $G_{0}$ being non-decreasing and convex guarantees the convexity of $G(X,Y)$. In fact, according to the well-known result that the composition of a non-decreasing convex function and a convex function is a convex function, and notice that $\|X^{(i)}\|^{2},\|Y^{(j)}\|^{2},\|X\|_{F}^{2},\|Y\|_{F}^{2}$ are convex, we have that each component of $G$ is convex and thus $G$ is convex.

We remark that (P1) can be interpreted as the penalized version of the following constrained problem (see, e.g. )

To illustrate this, note that the constraint $f_{1}(X)\triangleq\frac{3\|X\|_{F}^{2}}{2\beta_{T}^{2}}-1\leq 0$ corresponds to the penalty term $\rho G_{0}(f_{1}(X)+1)=\rho\max\{0,f_{1}(X)\}^{2}$ which appears as the third term in $G(X,Y)$, and similarly other constraints correspond to other terms in $G(X,Y)$. In other words, the regularization function $G(X,Y)$ is just a penalty function for the constraints of the problem (19). The function $\max\{0,\cdot\}^{2}$ is a popular choice for the penalty function in optimization (see, e.g. ), which motivates our choice of $G_{0}$ in (14). Our result can be extended to cover the algorithms for the constrained version (19), or a partially regularized formulation (e.g. only penalize the violation of the constraint $\|X\|_{F}^{2}\leq\frac{2}{3}\beta_{T}^{2},\|Y\|_{F}^{2}\leq\frac{2}{3}\beta_{T}^{2}$).

One commonly used assumption in the optimization literature is that the gradient of the objective function is Lipschitz continuous. For any positive number $\beta$, define a bounded set

The following result shows that this assumption (Lipschitz continuous gradients) holds for our objective function within a bounded set.

where $\|(X,Y)-(U,V)\|_{F}=\sqrt{\|X-U\|_{F}^{2}+\|Y-V\|_{F}^{2}}$.

The proof of Claim 2.1 is given in Appendix A.1.

3 Row-scaled Spectral Initialization

Our results require the initial point to be close enough to the global optima. To be more precise, we want the initial point to be in an incoherent neighborhood of the original matrix $M$ (this neighborhood will be specified later). Special initialization is also required in other works on non-convex formulations .

The initialization procedure is given in Table 1. The property of the initial point generated by this procedure will be presented in Claim 5.2.

In the numerical experiments, we find that the proposed initialization is not better than random initialization if we use the proposed formulation with the incoherence-control regularizer. In contrast, for traditional formulations (either unregularized or with a regularizer $\lambda(\|X\|_{F}^{2}+\|Y\|_{F}^{2})$) the proposed initialization does lead to better recovery performance (lower sample complexity). We also notice that the row-scaling step is crucial for this improvement since simply initializing via the spectral method does not help too much. See [30, Chapter 3] for the simulation results and discussions.

4 Algorithms

Our result applies to many standard algorithms such as gradient descent, SGD and block coordinate descent type methods (including alternating minimization, block coordinate gradient descent, block successive upper bound minimization, etc.). We will describe several typical algorithms in this subsection.

where $G_{0}^{\prime}(z)=I_{[1,\infty]}(z)2(z-1)$, and $\bar{X}^{(i)}$ (resp. $\bar{Y}^{(j)}$) denotes a matrix with the $i$-th (resp. $j$-th) row being $X^{(i)}$ (resp. $Y^{(j)}$) and the other rows being zero.

We first present a gradient descent algorithm in Table 2. There are many choices of stepsizes such as constant stepsize, exact line search, limited line search, diminishing stepsize and Armijo rule . We present three stepsize rules here: constant stepsize, restricted Armijo rule and restricted line search (the latter two are the variants of Armijo rule and exact line search). Note that the restricted line search rule is similar to that used in for the gradient descent method over Grassmannian manifolds. To simplify the notations, we denote $\bm{x}_{k}(\eta)\triangleq(X_{k}(\eta),Y_{k}(\eta))$ and $d(\bm{x}_{k}(\eta),\bm{x}_{0})\triangleq\sqrt{\|X_{k}(\eta)-X_{0}\|_{F}^{2}+\|Y_{k}(\eta)-Y_{0}\|_{F}^{2}}.$

AltMin (alternating minimization) belongs to the class of block coordinate descent (BCD) type methods. One can update the blocks in different orders (e.g. cyclic , randomized or parallel) and solve the subproblem inexactly. Commonly used inexact BCD type algorithms include BCGD (block coordinate gradient descent, which updates each variable by a single gradient step ) and BSUM (block successive upper bound minimization, which updates each variable by minimizing an upper bound of the objective function ). BCD-type methods have been widely used in engineering (e.g. ). In the context of matrix completion, Hastie et al. proposed an algorithm that could be viewed as a BSUM algorithm. Just considering different choices of the blocks will lead to different algorithms for the matrix completion problem . Our result applies to many BCD type methods, including the two-block alternating minimization, BCGD and BSUM. While it is not very interesting to list all possible algorithms to which our results are applicable, we just present two specific algorithms for illustration.

The first BCD type algorithm we present is (two-block) AltMin, which, in the context of matrix completion, usually refers to the algorithm that alternates between $X$ and $Y$ by updating one factor at a time with the other factor fixed. Although the overall objective function is non-convex, each subproblem of $X$ or $Y$ is convex and thus can be solved efficiently. The details are given in Table 3.

Theoretically speaking, AltMin for our formulation (P1) is not as efficient as the vanilla AtlMin for (P0) since an extra inner loop is needed to solve the subproblem. However, we remark that in the regimes of $|\Omega|$ that the vanilla AltMin works, the least square solution $X$ (resp. $Y$) is always bounded and incoherent (empirical observation), in which case the regularizer $G$ is inactive; therefore, the gradient updates in Table 4 do not happen. In the regimes of $|\Omega|$ that the vanilla AltMin fails, $G$ is active and the gradient updates do happen; however, instead of solving the subproblem exactly, one could perform one gradient step and the algorithm becomes the popular variant BCGD . Our main result of exact recovery still holds for BCGD (the proof for Algorithm 3 in Claim 5.3 can be applied to BCGD since BCGD is a special case of BSUM).

In the second BCD type algorithm called row BSUM, we update the rows of $X$ and $Y$ cyclically by minimizing an upper bound of the objective function; see Table 5. The extra terms $\frac{\lambda_{0}}{2}\|X^{(i)}-X_{k-1}^{(i)}\|^{2}$ or $\frac{\lambda_{0}}{2}\|Y^{(j)}-Y_{k-1}^{(j)}\|^{2}$ are added to make the subproblems strongly convex, which help prove convergence to stationary points. Such a technique has also been used in the alternating least square algorithm for tensor decomposition . Note that for the two-block BCD algorithm, convergence to stationary points can be guaranteed even when the subproblems are not strongly convex , thus in Algorithm 2 we do not add the extra terms. The benefit of cyclically updating the rows is that each subproblem can be solved efficiently using a simple binary search; see Appendix A.2 for the details. We remark again that instead of solving the subproblem exactly, one could just perform one gradient step to update each row of $X$ and $Y$ (with $\lambda=0$) and our result still holds.

The fourth algorithm we present is SGD (stochastic gradient descent) tailored for our problem (P1). In the optimization literature, this algorithm for minimizing the sum of finitely many functions is more commonly referred to as “incremental gradient method”, while SGD represents the algorithm for minimizing the expectation of a function; nevertheless, in this paper we follow the convention in the computer science literature and still call it “SGD”. In SGD, at each iteration we pick a component function and perform a gradient update. Similar to the BCD type methods where the blocks can be chosen in different orders, one can pick the component functions in a cyclic order, in an essentially cyclic order, or in a random order (either sampling with replacement or without replacement). In practice, the version of sampling without replacement converges much faster than the version of sampling with replacement (see [30, Chapter 2] for simulation results). In general, the understanding of sampling without replacement for optimization algorithms is quite limited (see, e.g., for one example of such analysis).

In this paper we only consider the cyclic order, and use a standard stepsize rule for SGD which requires the stepsizes $\{\eta_{k}\}$ to go to zero as $k\rightarrow\infty$, but neither too fast nor too slow (this choice guarantees convergence to stationary points even for nonconvex problems). One such choice of stepsizes is $\eta_{k}=O(1/k)$. We remark that our results also apply to other versions of SGD with different update orders or stepsize rules as long as they converge to stationary points.

and $\{f_{k}(X,Y)\}_{k=1}^{|\Omega|+m+n+2}$ denotes the collection of all component functions. With these definitions, the SGD algorithm is given in Table 6.

Main Results

The main result of this paper is that Algorithms 1-4 (standard optimization algorithms) will converge to the global optima of problem (P1) given in (18) and reconstruct $M$ exactly with high probability, provided that the number of revealed entries is large enough. Similar to the results for nuclear norm minimization , the probability is taken with respect to the random choice of $\Omega$, and the result also applies to a uniform random model of $\Omega$.

then with probability at least $1-2/n^{4}$, each of Algorithms 1-4 reconstructs $M$ exactly. Here, we say an algorithm reconstructs $M$ if each limit point $(X^{*},Y^{*})$ of the sequence $\{X_{k},Y_{k}\}$ generated by this algorithm satisfies $X^{*}(Y^{*})^{T}=M$.

This result shows that although (18) is a non-convex optimization problem, many standard algorithms can converge to the global optima with certain initialization. Different from all previous works on alternating minimization for matrix completion, our result does not require the algorithm to use independent samples in different iterations. To the best of our knowledge, our result is the first one that provides theoretical guarantee for alternating minimization without resampling. In addition, this result also provides the first exact recovery guarantee for many algorithms such as gradient descent, SGD and BSUM.

As demonstrated in (and proved in [5, Theorem 1.7]), $O(nr\log n)$ entries are the minimum requirement to recover the original matrix: $O(nr)$ is the number of degrees of freedom of a rank $r$ matrix $M$, and the additional $\log n$ factor is due to the coupon collector effect . For $r=O(1)$ and $\kappa$ bounded, Theorem 3.1 is order optimal in terms of the sample complexity since only $O(n\log n)$ entries are needed to exactly recover $M$. For $r=\mathcal{O}(\log n)$, however, our result is suboptimal by a polylogarithmic factor. The initialization has contributed $r^{4}\kappa^{4}$ to the sample complexity bound, and we expect that using other initialization procedures (e.g. the one proposed in ) can reduce the exponents of $r$ and $\kappa$.

(Linear convergence) Under the same condition of Theorem 3.1, with probability at least $1-2/n^{4}$, Algorithm 1a (gradient descent with constant stepsize) converges linearly; more precisely, the sequence $\{X_{k},Y_{k}\}$ generated by Algorithm 1a satisfies

where $\xi=\frac{1}{C_{g}r^{5}\kappa^{3}}p\Sigma_{\min}$ (here $C_{g}$ is a numerical constant), $\eta_{1}$ is the stepsize and $\eta_{1}\xi<1$.

Under the same condition of Theorem 3.1, with probability at least $1-1/(2n^{4})$, we have

The proof of this claim is given in Appendix D.2. This result is a simple corollary of several intermediate bounds established in the proof of Lemma 3.1.

To prove Theorem 3.1, we only need to prove two lemmas which describe the local geometry of the regularized objective in (P1) and the properties of the algorithms respectively. Roughly speaking, the first lemma shows that any stationary point of (P1) in a certain region is globally optimal, and the second lemma shows that each of Algorithms 1-4 converges to stationary points in that region. This region can be viewed as an “incoherent neighborhood” of $M$, and can be formally defined as $K_{1}\cap K_{2}\cap K(\delta)$, where $K_{1},K_{2}$ are defined as

Note that $K_{2}=\Gamma(\beta_{T})$ by our definition of $\Gamma$ in (20). As mentioned in Section 2.1, we only need to consider a Bernolli model of $\Omega$ where each entry is included into $\Omega$ with probability $p=\frac{S}{mn}$, where $S$ satisfies (27).

The first lemma describes the local geometry and implies that any stationary point $(X,Y)$ in $K_{1}\cap K_{2}\cap K(\delta)$ satisfies $XY^{T}=M$. The main steps to derive this geometrical property is described in Section 1.3. The formal proof will be given in Section 4.

The second lemma describes the properties of the algorithms we presented. Throughout the paper, “under the same condition of Lemma 3.1” means “assume $\delta$ is defined by (16) and $\Omega$ is generated by a Bernolli model with expected cardinality $S$ satisfying (27), where $C_{0},C_{d}$ are the same numerical constants as those in Lemma 3.1”. The proof of Lemma 3.2 will be given in Section 5.

Under the same conditions of Lemma 3.1, with probability at least $1-1/n^{4}$, the sequence $(X_{k},Y_{k})$ generated by either of Algorithms 1-4 has the following properties: (a) Each limit point of $(X_{k},Y_{k})$ is a stationary point of (P1). (b) $(X_{k},Y_{k})\in K_{1}\cap K_{2}\cap K(\delta),\;\forall k\geq 0$.

Intuitively, $\|X_{k}^{(i)}\|,\|Y_{k}^{(j)}\|,\|X_{k}\|_{F},\|Y_{k}\|_{F}$ are bounded because of the regularization terms we introduced and that the objective function is decreasing, and $\|M-X_{k}Y_{k}^{T}\|_{F}$ is bounded because the objective function is decreasing (however, the intuition is not enough and the proof requires some extra effort). In Section 5 we provide some easily verifiable conditions for Property (b) to hold (see Proposition 5.1), so that Lemma 3.2 and Theorem 3.1 can be extended to other algorithms.

With these two lemmas, the proof of Theorem 3.1 is quite straightforward and presented below.

Remark: Note that $X_{*}Y_{*}^{T}=M$ does not necessarily imply the global optimality of $(X_{*},Y_{*})$ since we have not proved $G(X_{*},Y_{*})=0$. Nevertheless, the global optimality can be easily proved using a different version of Lemma 3.1 (see the discussion before Lemma 3.3); in other words, Theorem 3.1 can be slightly strengthened to “Algorithm 1-4 converge to the global optima of problem (P1)”, instead of “Algorithm 1-4 recover $M$”.

The same argument can be used to show a more general result than Theorem 3.1, as stated in the following corollary.

Under the same conditions of Theorem 3.1, any algorithm satisfying Properties (a) and (b) in Lemma 3.2 reconstructs $M$ exactly with probability at least $1-2/n^{4}$.

2 Proof of Theorem 3.2

The proof of Theorem 3.2 applies a standard framework for first order methods: the convergence rate (or iteration complexity) can be derived from the “cost-to-go estimate” and the “sufficient descent” condition. For instance, the linear convergence $f(\bm{x}_{k})-f^{*}\leq(1-c_{1}c_{2})^{k}$ is a direct corollary of the cost-to-go estimate $\|\nabla f(\bm{x}_{k})\|^{2}\geq c_{1}[f(\bm{x}_{k})-f^{*}]$ and the sufficient descent condition $f(\bm{x}_{k})-f(\bm{x}_{k+1})\geq c_{2}\|\nabla f(\bm{x}_{k})\|^{2}$, where $f^{*}$ is the minimum value of $f$, and $c_{1},c_{2}$ are certain constants. We remark that using other optimization frameworks may lead to stronger time complexity bounds; this is left as future work.

(Cost-to-go estimate) Under the same conditions of Lemma 3.1, with probability at least $1-1/n^{4}$, the following holds:

where $\xi=\frac{1}{C_{g}r^{5}\kappa^{3}}p\Sigma_{\min}$ (here $C_{g}\geq 1$ is a numerical constant).

The following claim shows that Algorithm 1a satisfies the sufficient descent condition. It is easy to prove: it is well known that for minimizing a function (possibly non-convex) with Lipschitz continuous gradient, the gradient descent method with constant step-size satisfies the sufficient decrease condition.

(Sufficient descent) For the sequence $\bm{x}_{k}=(X_{k},Y_{k})$ generated by Algorithm 1a (gradient descent with constant stepsize), we have

where $\eta_{1}$ is the stepsize bounded above by $\bar{\eta_{1}}$ defined in (238).

The linear convergence can be easily derived from Lemma 3.1 and Claim 3.2. For completeness, we present the proof below.

Proof of Theorem 3.2: According to Property (b) of Lemma 3.2, with probability at least $1-1/n^{4}$, $(X_{k},Y_{k})\in K_{1}\cap K_{2}\cap K(\delta)$ for all $k$. According to Lemma 3.3 and Claim 3.2, we have (with probability at least $1-2/n^{4}$)

The stepsize $\eta_{1}$ can be bounded as $0<\eta_{1}\leq\bar{\eta}_{1}\overset{\eqref{eta 1 first bound}}{\leq}\frac{1}{4\beta_{T}^{2}}=\frac{1}{4C_{T}r\Sigma_{\max}}\leq\frac{1}{\Sigma_{\max}}$. Since $0<\xi=\frac{1}{C_{g}r^{5}\kappa^{3}}p\Sigma_{\min}\leq\Sigma_{\min}$, we have $0<\eta_{1}\xi\leq\frac{\Sigma_{\min}}{\Sigma_{\max}}\leq 1$, which implies $0<1-\frac{1}{2}\eta_{1}\xi<1$. Then the relation (35) leads to

The same argument can be used to show a more general result than Theorem 3.2, as stated in the following corollary.

Under the same conditions of Theorem 3.1, any algorithm satisfying Properties (a) and (b) in Lemma 3.2 and the sufficient decrease condition (34) has the linear convergence property, i.e. generates a sequence $(X_{k},Y_{k})$ that satisfies (28).

Proof of Lemma 3.1

In Section 4.1, we will show that to prove Lemma 3.1, we only need to construct $U,V$ to satisfy three inequalities that $\|\mathcal{P}_{\Omega}((U-X)(V-Y)^{T})\|_{F}$ and $\|((U-X)(V-Y)^{T})\|_{F}$ are bounded above and $\langle\nabla_{X}G,X-U\rangle+\langle\nabla_{Y}G,Y-V\rangle$ is bounded below. In Section 4.2 we describe two propositions that specify the choice of $U,V$, and then we show that such $U,V$ satisfy the three desired inequalities in Section 4.2 and subsequent subsections.

To ensure (32) holds, we only need to ensure that the following two inequalities hold:

Using the expressions of $\nabla_{X}F,\nabla_{Y}F$ in (24), we bound $\phi_{F}$ as follows:

The reason to decompose $M-XY^{T}$ as $a+b$ is the following. In order to bound $\|\mathcal{P}_{\Omega}(M-XY^{T})\|_{F}$, we notice $E(\mathcal{P}_{\Omega}(M-XY^{T}))=p(M-XY^{T})$ and wish to prove $\|\mathcal{P}_{\Omega}(M-XY^{T})\|_{F}^{2}\approx O(pd^{2}).$ However, $\|\mathcal{P}_{\Omega}(A)\|_{F}$ could be as large as $\|A\|_{F}$ if the matrix $A$ is not independent of the random subset $\Omega$ (e.g. choose $A$ s.t. $A=\mathcal{P}_{\Omega}(A)$). This issue can be resolved by decomposing $XY^{T}-M$ as $a+b$ and bounding $\|\mathcal{P}_{\Omega}(a)\|_{F}$ and $\|\mathcal{P}_{\Omega}(b)\|_{F}$ separately. In fact, $\|\mathcal{P}_{\Omega}(a)\|_{F}$ can be bounded because $a$ lies in a space spanned by the matrices with the same row space or column space as $M$, which is independent of $\Omega$ (Theorem 4.1 in ). $\|\mathcal{P}_{\Omega}(b)\|_{F}$ can be bounded according to a random graph lemma of , which requires $U,V,X,Y$ to be incoherent (i.e. have bounded row norm).

We claim that (37a) is implied by the following two inequalities:

In fact, assume (40a) and (40b) are true, we prove $\phi_{F}\geq pd^{2}/4$ as follows. By $XY^{T}-M=a+b$ we have

By Theorem 4.1 in , for $|\Omega|$ satisfying (27) with large enough $C_{0}$, we have that with probability at least $1-1/(2n^{4})$, $\|\mathcal{P}_{\mathcal{T}}\mathcal{P}_{\Omega}\mathcal{P}_{\mathcal{T}}(a)-p\mathcal{P}_{\mathcal{T}}(a)\|_{F}\leq\frac{1}{6}p\|a\|_{F}$ (note that this bound holds uniformly for all $a\in\mathcal{T}$, thus also holds when $a$ is dependent on $\Omega$). Since $a\in\mathcal{T}$, this inequality can be simplified to

Following the analysis of [4, Corollary 4.3], we have

The absolute value of the second term can be bounded as

which implies $-\frac{1}{6}p\|a\|_{F}^{2}\leq\langle a,\mathcal{P}_{\mathcal{T}}\mathcal{P}_{\Omega}(a)-pa\rangle\leq\frac{1}{6}p\|a\|_{F}^{2}$. Substituting into (44), we obtain that with probability at least $1-1/(2n^{4})$,

The first inequality of the above relation implies

According to (39) and the bounds (46) and (40a), we have $\phi_{F}/(pd^{2})\geq\frac{27}{40}+2(\frac{1}{5})^{2}-\frac{3}{5}\sqrt{\frac{27}{40}}\geq\frac{1}{4}$, which proves (37a).

In summary, to find a factorization $M=UV^{T}$ such that (32) holds, we only need to ensure that the factorization satisfies (40b), (40a) and (37b). In the following three subsections, we will show that such a factorization $M=UV^{T}$ exists. Specifically, $U,V$ will be defined in Table 7 and the three desired inequalities will be proved in Corollary 4.2, Proposition 4.3 and Claim 4.1 respectively.

2 Definitions of U,V𝑈𝑉U,V and key technical results

We construct $U,V$ according to two propositions, which will be stated in this subsection and proved in the appendix. The first proposition states that if $XY^{T}$ is close to $M$, then there exists a factorization $M=UV^{T}$ such that $U$ (resp. $V$) is close to $X$ (resp. $Y$), and $U,V$ are incoherent. Roughly speaking, this proposition shows the continuity of the factorization map $Z=XY^{T}\mapsto(X,Y)$ near a low-rank matrix $M$. The condition $X,Y\in K_{1}\cap K_{2}\cap K(\delta)$ and (16) implies that $d\triangleq\|M-XY^{T}\|_{F}\leq\delta=\frac{\Sigma_{\min}}{C_{d}r^{1.5}\kappa}$ and $\|X\|_{F}\leq\beta_{T},\|Y\|_{F}\leq\beta_{T}$, thus for large enough $C_{d}$, the assumptions of Proposition 4.1 hold. Similarly, the assumptions of the other results in this subsection also hold.

The proof of Proposition 4.1 is given in Appendix B.

Remark 1: A symmetric result that switches $X,U$ and $Y,V$ in the above proposition holds: under the conditions of Proposition (4.1), there exist $U,V$ satisfying (48) with $U,V$ reversed, i.e. $UV^{T}=M$, $\|V\|_{F}(1-\frac{d}{\Sigma_{\min}})\leq\|Y\|_{F}$, $\|U-X\|_{F}\leq\frac{3\beta_{T}}{\Sigma_{\rm min}}d,\|V-Y\|_{F}\leq\frac{6\beta_{T}}{5\Sigma_{\rm min}}d$, and $\|U^{(i)}\|^{2}\leq\frac{3r\mu}{2m}\beta_{T}^{2},\|V^{(j)}\|^{2}\leq\frac{r\mu}{n}\beta_{T}^{2}$.

Remark 2: To prove Theorem 3.1 (convergence), we only need $\|U\|_{F}\leq\|X\|_{F}$; here the slightly stronger requirement $\|U\|_{F}\leq(1-\frac{d}{\Sigma_{\min}})\|X\|_{F}$ is for the purpose of proving Theorem 3.2 (linear convergence).

Remark 3: Without the incoherence assumption on $M$, by the same proof we can show that there still exist $U,V$ satisfying (48a) and (48c), i.e. $M=UV^{T}$ and $U,V$ are close to $X,Y$ respectively. Such a result bears some similarity with the classical perturbation theory for singular value decomposition . In particular, proved that for two low-rank matricesThe result in also covered the case of two approximately low-rank matrices, but we only consider the case of exact low-rank matrices here. that are close, the spaces spanned by the left (resp. right) singular vectors of the two matrices are also close. Note that the singular vectors themselves may be very sensitive to perturbations and no such perturbation bounds can be established (see [63, Sec. 6]). The difference of our work with the classical perturbation theory is that we do not consider SVD of two matrices; instead, we allow one matrix to have an arbitrary factorization, and the factorization of the other matrix can be chosen accordingly. Since we do not have any restriction on the factorization $XY^{T}$ (except the dimensions) and the norms of $X$ and $Y$ can be arbitrarily large, the distance between two corresponding factors has to be proportional to the norm of one single factor, which explains the coefficient $\beta_{T}$ in (48c).

Unfortunately, Proposition 4.1 is not strong enough to prove $\phi_{G}\geq 0$ when both $\|X\|_{F}$ and $\|Y\|_{F}$ are large (see an analysis in Section 4.4). To resolve this issue, we need to prove the second proposition in which there is an additional assumption that both $\|X\|_{F}$ and $\|Y\|_{F}$ are large, and an additional requirement that both $\|U\|_{F}$ and $\|V\|_{F}$ are bounded (by the norms of original factors $\|X\|_{F}$ and $\|Y\|_{F}$ respectively). More specifically, the proposition states that if $M$ is close to $XY^{T}$, and both $\|X\|_{F}$ and $\|Y\|_{F}$ are large, then there is a factorization $M=UV^{T}$ such that $U$ (resp.$\ V$) is close to $X$ (resp.$\ Y$), and $\|U\|_{F}\leq\|X\|_{F},\|V\|_{F}\leq\|Y\|_{F}$. For the purpose of proving linear convergence, we prove a slightly stronger result that $\|V\|_{F}\leq(1-d/\Sigma_{\min})\|Y\|_{F}$ . The previous result Proposition 4.1 can be viewed as a perturbation analysis for an arbitrary factorization, while Proposition 4.2 can be viewed as an enhanced perturbation analysis for a constrained factorization. Although Proposition 4.2 is just a simple variant of Proposition 4.1, it seems to require a much more involved proof than Proposition 4.1. See the formal proof of Proposition 4.2 in Appendix C.

Remark: A symmetric result that switches $X,U$ and $Y,V$ in the above proposition still holds; the only change is that (50b) will become $\|U\|_{F}\leq(1-\frac{d}{\Sigma_{\min}})\|X\|_{F},\;\|V\|_{F}\leq\|Y\|_{F}$. It is easy to prove a variant of the above proposition in which (50b) is changed to $\|U\|_{F}\leq(1-\frac{d}{2\Sigma_{\min}})\|X\|_{F},\|V\|_{F}\leq(1-\frac{d}{2\Sigma_{\min}})\|Y\|_{F}$; in other words, the asymmetry of $X,U$ and $Y,V$ in (50b) is artificial. Nevertheless, Proposition 4.2 is enough for our purpose.

Throughout the proof of Lemma 3.1, $U,V$ are defined in Table 4.2.

According to Proposition 4.1 and Proposition 4.2 (and their symmetric results), the properties of $U,V$ defined in Tabel 7 are summarized in the following corollary. For simplicity, we only present the case that $\|X\|_{F}\leq\|Y\|_{F}$; in the other case that $\|X\|_{F}>\|Y\|_{F}$, a symmetric result of Corollary 4.1 holds.

Suppose $d\triangleq\|XY^{T}-M\|_{F}\leq\frac{\Sigma_{\min}}{C_{d}r}$ and $\|X\|_{F}\leq\|Y\|_{F}$, then $U,V$ defined in Table 7 satisfy:

In (51b), we bound $\|U-X\|_{F}\|V-Y\|_{F}$ by $O(d^{2})$ with a rather complicated coefficient, but to prove (40b) we need a bound $O(d)$ with a coefficient $1/10$. Under a slightly stronger condition on $d$ than that of Corollary 4.1, which still holds for $(X,Y)\in K(\delta)$ with $\delta$ defined in (16), we can prove the bound (40b) by (51b).

There exists a numerical constant $C_{d}$ such that if

then $U,V$ defined in Table 7 satisfy (40b).

Proof of Corollary 4.2: According to (51b) , we have

where the last inequliaty follows from (52) with $C_{d}\geq 650C_{T}.$ $\Box$

In the next two subsections, we will use the properties in Corollary 4.1 to prove (40a) and (37b).

The following result states that for $U,V$ defined in Table 7, (40a) holds.

Under the same conditions as Lemma 3.1, with probability at least $1-1/(2n^{4})$, the following is true. For any $(X,Y)\in K_{1}\cap K_{2}\cap K(\delta)$ and $U,V$ defined in Table 7, we have

Proof of Proposition 4.3: We need the following random graph lemma [31, Lemma 7.1].

Let $Z=U-X,W=V-Y$ and $z_{i}=\|Z^{(i)}\|^{2}$, $w_{j}=\|W^{(j)}\|^{2}$. We have

Analogous to the proof of (40b) in Corollary 4.2, we can prove that $\|U-X\|_{F}\|V-Y\|_{F}\leq d/(10\sqrt{C_{1}})$ for large enough $C_{d}$ (in fact, $C_{d}\geq 650C_{T}\sqrt{C_{1}}$ suffices). Therefore, we have

We still need to bound $\|z\|_{2}$ and $\|w\|_{2}.$ We have

Here, the third inequliaty follows from the property (51c) in Corollary 4.1 and the condition $(X,Y)\in K_{1}$ (which implies $\|X^{(i)}\|\leq\beta_{1}$), and the fourth inequliaty follows from the definition of $\beta_{1}$ in (15). Similarly,

Thus the second term in (56) can be bounded as

where the last inequality is equivalent to $520^{2}C_{1}^{2}C_{T}^{4}\alpha^{\frac{3}{2}}\mu^{2}r^{7}\kappa^{4}\leq\frac{9}{100^{2}}|\Omega|/n$, which holds due to (27) with large enough numerical constant $C_{0}$. Plugging (57) and (60) into (56), we get $\|\mathcal{P}_{\Omega}((U-X)(V-Y)^{T})\|_{F}^{2}\leq\frac{p}{25}d^{2}=\frac{p}{25}\|M-XY^{T}\|_{F}^{2}$. $\Box$

In this subsection, we prove the following claim.

$U,V$ defined in Table 7 satisfy (37b), i.e. $\phi_{G}=\langle\nabla_{X}G,X-U\rangle+\langle\nabla_{Y}G,Y-V\rangle\geq 0$.

By the expressions of $\nabla_{X}G,\nabla_{Y}G$ in (24), we have

where $G_{0}^{\prime}(z)=I_{[1,\infty]}(z)2(z-1)$.

We only need to prove (62a); the proof of (62b) is similar. We consider two cases.

Case 1: $\|X^{(i)}\|^{2}\leq\frac{2\beta_{1}^{2}}{3}.$ Note that $\frac{3\|X^{(i)}\|^{2}}{2\beta_{1}^{2}}\leq 1$ implies $G_{0}^{\prime}(\frac{3\|X^{(i)}\|^{2}}{2\beta_{1}^{2}})=0$, thus $h_{1i}=0$.

Case 2: $\|X^{(i)}\|^{2}>\frac{2\beta_{1}^{2}}{3}.$ By Corollary 4.1 and the fact that $\beta_{1}^{2}=\beta_{T}^{2}\frac{3\mu r}{m}$, we have

As a result, $\langle X^{(i)},X^{(i)}\rangle=\|X^{(i)}\|\|X^{(i)}\|>\|X^{(i)}\|\|U^{(i)}\|\geq\langle X^{(i)},U^{(i)}\rangle,$ which implies $\langle X^{(i)},X^{(i)}-U^{(i)}\rangle\geq 0$. Combining this inequality with the fact that $G_{0}^{\prime}(\frac{3\|X^{(i)}\|^{2}}{2\beta_{1}^{2}})\geq 0,$ we get $h_{1i}\geq 0.$

Without loss of generality, we can assume $\|X\|_{F}\leq\|Y\|_{F},$ and we will apply Corollary 4.1 to prove (64). If $\|Y\|_{F}<\|X\|_{F}$, we can apply a symmetric result of Corollary 4.1 to prove (64). We further consider three cases.

Case 1: $\|X\|_{F}\leq\|Y\|_{F}\leq\sqrt{\frac{2}{3}}\beta_{T}.$ In this case $G_{0}^{\prime}(\frac{3\|X\|_{F}^{2}}{2\beta_{T}^{2}})=G_{0}^{\prime}(\frac{3\|Y\|_{F}^{2}}{2\beta_{T}^{2}})=0$, which implies $h_{2}=h_{4}=0$, thus $\eqref{h_2, h_4 >=0}$ holds.