Matrix Completion from a Few Entries

Introduction

Imagine that each of $m$ customers watches and rates a subset of the $n$ movies available through a movie rental service. This yields a dataset of customer-movie pairs $(i,j)\in E\subseteq[m]\times[n]$ and, for each such pair, a rating $M_{ij}\in{\mathds{R}}$. The objective of collaborative filtering is to predict the rating for the missing pairs in such a way as to provide targeted suggestions.Indeed, in 2006, Netflix made public such a dataset with $m\approx 5\cdot 10^{5}$, $n\approx 2\cdot 10^{4}$ and $|E|\approx 10^{8}$ and challenged the research community to predict the missing ratings with root mean square error below $0.8563$ [Net]. The general question we address here is: Under which conditions do the known ratings provide sufficient information to infer the unknown ones? Can this inference problem be solved efficiently? The second question is particularly important in view of the massive size of actual data sets.

A simple mathematical model for such data assumes that the (unknown) matrix of ratings has rank $r\ll m,n$. More precisely, we denote by $M$ the matrix whose entry $(i,j)\in[m]\times[n]$ corresponds to the rating user $i$ would assign to movie $j$. We assume that there exist matrices $U$, of dimensions $m\times r$, and $V$, of dimensions $n\times r$, and a diagonal matrix $\Sigma$, of dimensions $r\times r$ such that

For justification of these assumptions and background on the use of low rank matrices in information retrieval, we refer to [BDJ99]. Since we are interested in very large data sets, we shall focus on the limit $m,n\to\infty$ with $m/n=\alpha$ bounded away from and $\infty$.

We further assume that the factors $U$, $V$ are unstructured. This notion is formalized by the incoherence condition introduced by Candés and Recht [CR08], and defined in Section 2. In particular the incoherence condition is satisfied with high probability if $M=U\Sigma V^{T}$ with $U$ and $V$ uniformly random matrices with $U^{T}U=m{\mathbf{1}}$ and $V^{T}V=n{\mathbf{1}}$. Alternatively, incoherence holds if the entries of $U$ and $V$ are i.i.d. bounded random variables.

Out of the $m\times n$ entries of $M$, a subset $E\subseteq[m]\times[n]$ (the user/movie pairs for which a rating is available) is revealed. We let $M^{E}$ be the $m\times n$ matrix that contains the revealed entries of $M$, and is filled with ’s in the other positions

The set $E$ will be uniformly random given its size $|E|$.

2 Algorithm

A naive algorithm consists of the following projection operation.

Projection. Compute the singular value decomposition (SVD) of $M^{E}$ (with $\sigma_{1}\geq\sigma_{2}\geq\cdots\geq 0$)

and return the matrix ${\sf T}_{r}(M^{E})=(mn/|E|)\sum_{i=1}^{r}\sigma_{i}x_{i}y_{i}^{T}$ obtained by setting to all but the $r$ largest singular values. Notice that, apart from the rescaling factor $(mn/|E|)$, ${\sf T}_{r}(M^{E})$ is the orthogonal projection of $M^{E}$ onto the set of rank-$r$ matrices. The rescaling factor compensates the smaller average size of the entries of $M^{E}$ with respect to $M$.

It turns out that, if $|E|=\Theta(n)$, this algorithm performs very poorly. The reason is that the matrix $M^{E}$ contains columns and rows with $\Theta(\log n/\log\log n)$ non-zero (revealed) entries. The largest singular values of $M^{E}$ are of order $\Theta(\sqrt{\log n/\log\log n})$. The corresponding singular vectors are highly concentrated on high-weight column or row indices (respectively, for left and right singular vectors). Such singular vectors are an artifact of the high-weight columns/rows and do not provide useful information about the hidden entries of $M$. This motivates the definition of the following operation (hereafter the degree of a column or of a row is the number of its revealed entries).

Trimming. Set to zero all columns in $M^{E}$ with degree larger that $2|E|/n$. Set to all rows with degree larger than $2|E|/m$.

Figure 1 shows the singular value distributions of $M^{E}$ and $\widetilde{M}^{E}$ for a random rank-$3$ matrix $M$. The surprise is that trimming (which amounts to ‘throwing out information’) makes the underlying rank-$3$ structure much more apparent. This effect becomes even more important when the number of revealed entries per row/column follows a heavy tail distribution, as for real data.

In terms of the above routines, our algorithm has the following structure.

The last step of the above algorithm allows to reduce (or eliminate) small discrepancies between ${\sf T}_{r}(\widetilde{M}^{E})$ and $M$, and is described below.

Cleaning. Various implementations are possible, but we found the following one particularly appealing. Given $X\in{\mathds{R}}^{m\times r}$, $Y\in{\mathds{R}}^{n\times r}$ with $X^{T}X=m{\mathbf{1}}$ and $Y^{T}Y=n{\mathbf{1}}$, we define

The cleaning step consists in writing ${\sf T}_{r}(\widetilde{M}^{E})=X_{0}S_{0}Y_{0}^{T}$ and minimizing $F(X,Y)$ locally with initial condition $X=X_{0}$, $Y=Y_{0}$.

Notice that $F(X,Y)$ is easy to evaluate since it is defined by minimizing the quadratic function $S\mapsto{\cal F}(X,Y,S)$ over the low-dimensional matrix $S$. Further it depends on $X$ and $Y$ only through their column spaces. In geometric terms, $F$ is a function defined over the cartesian product of two Grassmann manifolds (we refer to Section 6 for background and references). Optimization over Grassmann manifolds is a well understood topic [EAS99] and efficient algorithms (in particular Newton and conjugate gradient) can be applied. To be definite, we assume that gradient descent with line search is used to minimize $F(X,Y)$.

Finally, the implementation proposed here implicitly assumes that the rank $r$ is known. In practice this is a non-issue. Since $r\ll n$, a loop over the value of $r$ can be added at little extra cost. For instance, in collaborative filtering applications, $r$ ranges between $10$ and $30$.

3 Main results

Notice that computing ${\sf T}_{r}(\widetilde{M}^{E})$ only requires to find the first $r$ singular vectors of a sparse matrix. Our main result establishes that this simple procedure achieves arbitrarily small relative root mean square error from $O(nr)$ revealed entries. We define the relative root mean square error as

where we denote by $||A||_{F}$ the Frobenius norm of matrix $A$. Notice that the factor $(1/mn)$ corresponds to the usual normalization by the number of entries and the factor $(1/M_{\rm max}^{2})$ corresponds to the maximum size of the matrix entries where $M$ satisfies $|M_{i,j}|\leq M_{\rm max}$ for all $i$ and $j$.

Assume $M$ to be a rank $r$ matrix of dimension $n\alpha\times n$ that satisfies $|M_{i,j}|\leq M_{\rm max}$ for all $i,j$. Then with probability larger than $1-1/n^{3}$

Notice that the top $r$ singular values and singular vectors of the sparse matrix $\widetilde{M}^{E}$ can be computed efficiently by subspace iteration [Ber92]. Each iteration requires $O(|E|r)$ operations. As proved in Section 3, the $(r+1)$-th singular value is smaller than one half of the $r$-th one. As a consequence, subspace iteration converges exponentially. A simple calculation shows that $O(\log n)$ iterations are sufficient to ensure the error bound mentioned.

The ‘cleaning’ step in the above pseudocode improves systematically over ${\sf T}_{r}(\widetilde{M}^{E})$ and, for large enough $|E|$, reconstructs $M$ exactly.

Assume $M$ to be a rank $r$ matrix that satisfies the incoherence conditions A1 and A2 with $(\mu_{0},\mu_{1})$. Let $\mu=\max\{\mu_{0},\mu_{1}\}$. Further, assume $\Sigma_{\rm min}\leq\Sigma_{1},\dots,\Sigma_{r}\leq\Sigma_{\rm max}$ with $\Sigma_{\rm min},\Sigma_{\rm max}$ bounded away from and $\infty$. Then there exists a numerical constant $C^{\prime}$ such that, if

then the cleaning procedure in Spectral Matrix Completion converges, with high probability, to the matrix $M$.

This theorem is proved in Section 6. The basic intuition is that, for $|E|\geq C^{\prime}(\alpha)nr\,\max\{\log n,r\}$, $T_{r}(\widetilde{M}^{E})$ is so close to $M$ that the cost function is well approximated by a quadratic function.

Theorem 1.1 is optimal: the number of degrees of freedom in $M$ is of order $nr$, without the same number of observations is impossible to fix them. The extra $\log n$ factor in Theorem 1.2 is due to a coupon-collector effect [CR08, KMO08, KOM09]: it is necessary that $E$ contains at least one entry per row and one per column and this happens only for $|E|\geq Cn\log n$. As a consequence, for rank $r$ bounded, Theorem 1.2 is optimal. It is suboptimal by a polylogarithmic factor for $r=O(\log n)$.

4 Related work

Beyond collaborative filtering, low rank models are used for clustering, information retrieval, machine learning, and image processing. In [Faz02], the NP-hard problem of finding a matrix of minimum rank satisfying a set of affine constraints was addresses through convex relaxation. This problem is analogous to the problem of finding the sparsest vector satisfying a set of affine constraints, which is at the heart of compressed sensing [Don06, CRT06]. The connection with compressed sensing was emphasized in [RFP07], that provided performance guarantees under appropriate conditions on the constraints.

In the case of collaborative filtering, we are interested in finding a matrix $M$ of minimum rank that matches the known entries $\{M_{ij}:\,(i,j)\in E\}$. Each known entry thus provides an affine constraint. Candès and Recht [CR08] introduced the incoherent model for $M$. Within this model, they proved that, if $E$ is random, the convex relaxation correctly reconstructs $M$ as long as $|E|\geq C\,r\,n^{6/5}\log n$. On the other hand, from a purely information theoretic point of view (i.e. disregarding algorithmic considerations), it is clear that $|E|=O(n\,r)$ observations should allow to reconstruct $M$ with arbitrary precision. Indeed this point was raised in [CR08] and proved in [KMO08], through a counting argument.

The present paper describes an efficient algorithm that reconstructs a rank-$r$ matrix from $O(n\,r)$ random observations. The most complex component of our algorithm is the SVD in step $2$. We were able to treat realistic data sets with $n\approx 10^{5}$. This must be compared with the $O(n^{4})$ complexity of semidefinite programming [CR08].

Cai, Candès and Shen [CCS08] recently proposed a low-complexity procedure to solve the convex program posed in [CR08]. Our spectral method is akin to a single step of this procedure, with the important novelty of the trimming step that improves significantly its performances. Our analysis techniques might provide a new tool for characterizing the convex relaxation as well.

Theorem 1.1 can also be compared with a copious line of work in the theoretical computer science literature [FKV04, AFK+01, AM07]. An important motivation in this context is the development of fast algorithms for low-rank approximation. In particular, Achlioptas and McSherry [AM07] prove a theorem analogous to 1.1, but holding only for $|E|\geq(8\log n)^{4}n$ (in the case of square matrices).

A short account of our results was submitted to the 2009 International Symposium on Information Theory [KOM09]. While the present paper was under completion, Cándes and Tao posted online a preprint proving a theorem analogous to 1.2 [CT09]. Once more, their approach is substantially different from ours.

5 Open problems and future directions

It is worth pointing out some limitations of our results, and interesting research directions:

1. Optimal RMSE with $O(n)$ entries. Numerical simulations with the Spectral Matrix Completion algorithm suggest that the RMSE decays much faster with the number of observations per degree of freedom $(|E|/nr)$, than indicated by Eq. (9). This improved behavior is a consequence of the cleaning step in the algorithm. It would be important to characterize the decay of RMSE with $(|E|/nr)$.

2. Threshold for exact completion. As pointed out, Theorem 1.2 is order optimal for $r$ bounded. It would nevertheless be useful to derive quantitatively sharp estimates in this regime. A systematic numerical study was initiated in [KMO08]. It appears that available theoretical estimates (including the recent ones in [CT09]) are for larger values of the rank, we expect that our arguments can be strenghtened to prove exact reconstruction for $|E|\geq C^{\prime}(\alpha)nr\log n$ for all values of $r$.

3. More general models. The model studied here and introduced in [CR08] presents obvious limitations. In applications to collaborative filtering, the subset of observed entries $E$ is far from uniformly random. A recent paper [SC09] investigates the uniqueness of the solution of the matrix completion problem for general sets $E$. In applications to fast low-rank approximation, it would be desirable to consider non-incoherent matrices as well (as in [AM07]).

Incoherence property and some notations

In order to formalize the notion of incoherence, we write $U=[u_{1},u_{2},\dots,u_{r}]$ and $V=[v_{1},v_{2},\dots,v_{r}]$ for the columns of the two factors, with $||u_{i}||=\sqrt{m}$, $||v_{i}||=\sqrt{n}$ and $u_{i}^{T}u_{j}=0$, $v_{i}^{T}v_{j}=0$ for $i\neq j$ (there is no loss of generality in this, since normalizations can be adsorbed by redefining $\Sigma$). We shall further write $\Sigma={\rm diag}(\Sigma_{1},\dots,\Sigma_{r})$ with $\Sigma_{1}\geq\Sigma_{2}\geq\cdots\geq\Sigma_{r}>0$.

The matrices $U$, $V$ and $\Sigma$ will be said to be $(\mu_{0},\mu_{1})$-incoherent if they satisfy the following properties:

For all $i\in[m]$, $j\in[n]$, we have $\sum_{k=1}^{r}{U_{i,k}^{2}}\leq\mu_{0}r$, $\sum_{k=1}^{r}{V_{i,k}^{2}}\leq\mu_{0}r$.

For all $i\in[m]$, $j\in[n]$, we have $|\sum_{k=1}^{r}{U_{i,k}(\Sigma_{k}/\Sigma_{1})V_{j,k}}|\leq\mu_{1}r^{1/2}$.

Apart from difference in normalization, these assumptions coincide with the ones in [CR08].

Probability is taken with respect to the uniformly random subset $E\subseteq[m]\times[n]$. Define $\epsilon\equiv|E|/\sqrt{mn}$. In the case when $m=n$, $\epsilon$ corresponds to the average number of revealed entries per row or column. Then, it is convenient to work with a model in which each entry is revealed independently with probability $\epsilon/\sqrt{mn}$. Since, with high probability $|E|\in[\epsilon\sqrt{\alpha}\,n-A\sqrt{n\log n},\epsilon\sqrt{\alpha}\,n+A\sqrt{n\log n}]$, any guarantee on the algorithm performances that holds within one model, holds within the other model as well if we allow for a vanishing shift in $\epsilon$.

Notice that we can assume $m\geq n$, since we can always apply our theorem to the transpose of the matrix $M$. Throughout this paper, therefore, we will assume $\alpha\geq 1$. Finally, we will use $C$, $C^{\prime}$ etc. to denote numerical constants.

Given a vector $x\in{\mathds{R}}^{n}$, $||x||$ will denote its Euclidean norm. For a matrix $X\in{\mathds{R}}^{n\times n^{\prime}}$, $||X||_{F}$ is its Frobenius norm, and $||X||_{2}$ its operator norm (i.e. $||X||_{2}=\sup_{u\neq 0}||Xu||/||u||$). The standard scalar product between vectors or matrices will sometimes be indicated by $\langle x,y\rangle$ or $\langle X,Y\rangle$, respectively. Finally, we use the standard combinatorics notation $[N]=\{1,2,\dots,N\}$ to denote the set of first $N$ integers.

Proof of Theorem 1.1 and technical results

As explained in the previous section, the crucial idea is to consider the singular value decomposition of the trimmed matrix $\widetilde{M}^{E}$ instead of the original matrix $M^{E}$, as in Eq. (5). We shall then redefine $\{\sigma_{i}\}$, $\{x_{i}\}$, $\{y_{i}\}$, by letting

Here $||x_{i}||=||y_{i}||=1$, $x_{i}^{T}x_{j}=y_{i}^{T}y_{j}=0$ for $i\neq j$ and $\sigma_{1}\geq\sigma_{2}\geq\dots\geq 0$. Our key technical result is that, apart from a trivial rescaling, these singular values are close to the ones of the full matrix $M$.

There exists a numerical constant $C>0$ such that, with probability larger than $1-1/n^{3}$

where it is understood that $\Sigma_{q}=0$ for $q>r$.

This result generalizes a celebrated bound on the second eigenvalue of random graphs [FKS89, FO05] and is illustrated in Fig. 1: the spectrum of $\widetilde{M}^{E}$ clearly reveals the rank-$3$ structure of $M$.

As shown in Section 5, Lemma 3.1 is a direct consequence of the following estimate.

There exists a numerical constant $C>0$ such that, with probability larger than $1-1/n^{3}$

The proof of this lemma is given in Section 4.

where we used Lemma 3.2 for the second inequality and Lemma 3.1 for the last inequality. Now, for any matrix $A$ of rank at most $2r$, $||A||_{F}\leq\sqrt{2r}||A||_{2}$, whence

The result follows by using $|E|=\epsilon\sqrt{mn}$.

Proof of Lemma 3.2

We want to show that $|x^{T}(\widetilde{M}^{E}-\frac{\epsilon}{\sqrt{mn}}M)y|\leq CM_{\rm max}\sqrt{\alpha\epsilon}$ for each $x\in{\mathds{R}}^{m}$, $y\in{\mathds{R}}^{n}$ such that $||x||=||y||=1$. Our basic strategy (inspired by [FKS89]) will be the following: $(1)$ Reduce to $x$, $y$ belonging to discrete sets $T_{m}$, $T_{n}$; $(2)$ Bound the contribution of light couples by applying union bound to these discretized sets, with a large deviation estimate on the random variable $Z$, defined as $Z\equiv\sum_{L}{x_{i}\widetilde{M}^{E}_{i,j}y_{j}}-\frac{\epsilon}{\sqrt{mn}}x^{T}My$; $(3)$ Bound the contribution of heavy couples using bound on the discrepancy of corresponding graph.

The technical challenge is that a worst-case bound on the tail probability of $Z$ is not good enough, and we must keep track of its dependence on $x$ and $y$. The definition of light and heavy couples is provided in the following section.

Notice that $T_{n}\subseteq S_{n}\equiv\{x\in{\mathds{R}}^{n}:\,||x||\leq 1\}$. Next remark is proved in [FKS89, FO05], and relates the original problem to the discretized one.

Let $R\in{\mathds{R}}^{m\times n}$ be a matrix. If $|x^{T}Ry|\leq B$ for all $x\in T_{m}$ and $y\in T_{n}$, then $|x^{\prime T}Ry^{\prime}|\leq(1-\Delta)^{-2}B$ for all $x^{\prime}\in S_{m}$ and $y^{\prime}\in S_{n}$.

Hence it is enough to show that, with high probability, $|x^{T}(\widetilde{M}^{E}-\frac{\epsilon}{\sqrt{mn}}M)y|\leq CM_{\rm max}\sqrt{\alpha\epsilon}$ for all $x\in T_{m}$ and $y\in T_{n}$.

A naive approach would be to apply concentration inequalities directly to the random variable $x^{T}(\widetilde{M}^{E}-\frac{\epsilon}{\sqrt{mn}}M)y$. This fails because the vectors $x$, $y$ can contain entries that are much larger than the typical size $O(n^{-1/2})$. We thus separate two contributions. The first contribution is due to light couples $L\subseteq[m]\times[n]$, defined as

The second contribution is due to its complement $\overline{L}$, which we call heavy couples. We have

In the next two subsections, we will prove that both contributions are upper bounded by $CM_{\rm max}\sqrt{\alpha\epsilon}$ for all $x\in T_{m}$, $y\in T_{n}$. Applying Remark 4.1 to $|x^{T}(\widetilde{M}^{E}-\frac{\epsilon}{\sqrt{mn}}M)y|$, this proves the thesis.

2 Bounding the contribution of light couples

Let us define the subset of row and column indices which have not been trimmed as ${\cal A}_{l}$ and ${\cal A}_{r}$:

where $\deg(\cdot)$ denotes the degree (number of revealed entries) of a row or a column. Notice that ${\cal A}=({\cal A}_{l},{\cal A}_{r})$ is a function of the random set $E$. It is easy to get a rough estimate of the sizes of ${\cal A}_{l}$, ${\cal A}_{r}$.

There exists $C_{1}$ and $C_{2}$ depending only on $\alpha$ such that, with probability larger than $1-1/n^{4}$, $|{\cal A}_{l}|\geq m-\max\{e^{-C_{1}\epsilon}m,C_{2}\alpha\}$, and $|{\cal A}_{r}|\geq n-\max\{e^{-C_{1}\epsilon}n,C_{2}\}$.

For the proof of this claim, we refer to Appendix A. For any $E\subseteq[m]\times[n]$ and $A=(A_{l},A_{r})$ with $A_{l}\subseteq[m]$, $A_{r}\subseteq[n]$, we define $M^{E,A}$ by setting to zero the entries of $M$ that are not in $E$, those whose row index is not in $A_{l}$, and those whose column index not in $A_{r}$. Consider the event

with $\delta\equiv\max\{e^{-C_{1}\epsilon},C_{2}\alpha\}$ and $H(x)$ the binary entropy function.

Let $x\in S_{m}$, $y\in S_{n}$, $Z=\sum_{(i,j)\in L}x_{i}M^{E,A}_{ij}y_{j}-\frac{\epsilon}{\sqrt{mn}}x^{T}My$, and assume $|A_{l}|\geq m(1-\delta)$, $|A_{r}|\geq n(1-\delta)$ with $\delta$ small enough. Then

We begin by bounding the mean of $Z$ as follows (for the proof of this statement we refer to Appendix B).

For $A=(A_{l},A_{r})$, let $M^{A}$ be the matrix obtained from $M$ by setting to zero those entries whose row index is not in $A_{l}$, and those whose column index not in $A_{r}$. Define the potential contribution of the light couples $a_{ij}$ and independent random variables $Z_{ij}$ as

Let $Z_{1}=\sum_{i,j}{Z_{ij}}$ so that $Z=Z_{1}-\frac{\epsilon}{\sqrt{mn}}x^{T}My$. Note that $\sum_{i,j}{a_{ij}^{2}}\leq\sum_{i,j}{\left(x_{i}M^{A}_{ij}y_{j}\right)^{2}}\leq M_{\rm max}^{2}$. Fix $\lambda=\sqrt{mn}/2M_{\rm max}\sqrt{\epsilon}$ so that $|\lambda a_{i,j}|\leq 1/2$, whence $e^{\lambda a_{ij}}-1\leq\lambda a_{ij}+2(\lambda a_{ij})^{2}$. It then follows that

Hence, assuming $\alpha\geq 1$, there exists a numerical constant $C^{\prime}$ such that, for $C>C^{\prime}\sqrt{\alpha}$, the first term is of order $e^{-\Theta(n)}$, and this finishes the proof.

3 Bounding the contribution of heavy couples

Let $Q$ be an $m\times n$ matrix with $Q_{ij}=1$ if $(i,j)\in E$ and $i\not\in{\cal A}_{r}$, $j\not\in{\cal A}_{l}$ (i.e. entry $(i,j)$ is not trimmed by our algorithm), and $Q_{ij}=0$ otherwise. Since $|M_{ij}|\leq M_{\rm max}$, the heavy couples satisfy $|x_{i}y_{j}|\geq\sqrt{\epsilon/mn}$. We then have

Notice that $Q$ is the adjacency matrix of a random bipartite graph with vertex sets $[m]$ and $[n]$ and maximum degree bounded by $2\epsilon\max(\alpha^{1/2},\alpha^{-1/2})$. The following remark strengthens a result of [FO05].

Given vectors $x$, $y$, let $\overline{L}^{\prime}=\{(i,j):|x_{i}y_{j}|\geq C\sqrt{\epsilon/mn}\}$. Then there exist a constant $C^{\prime}$ such that, $\sum_{(i,j)\in\overline{L}^{\prime}}Q_{ij}|x_{i}y_{j}|\leq C^{\prime}(\sqrt{\alpha}+\frac{1}{\sqrt{\alpha}})\sqrt{\epsilon}$, for all $x\in T_{m}$, $y\in T_{n}$ with probability larger than $1-1/2n^{3}$.

For the reader’s convenience, a proof of this fact is proposed in Appendix C. The analogous result in [FO05] (for the adjacency matrix of a non-bipartite graph) is proved to hold only with probability larger than $1-e^{-C\epsilon}$. The stronger statement quoted here can be proved using concentration of measure inequalities. The last remark implies that for all $x\in T_{m}$, $y\in T_{n}$, and $\alpha\geq 1$, the contribution of heavy couples is bounded by $CM_{\rm max}\sqrt{\alpha\epsilon}$ for some numerical constant $C$ with probability larger than $1-1/2n^{3}$.

Proof of Lemma 3.1

Recall the variational principle for the singular values.

Here $H$ is understood to be a linear subspace of ${\mathds{R}}^{n}$.

Using Eq. (19) with $H$ the orthogonal complement of ${\rm span}(v_{1},\ldots,v_{q-1})$, we have, by Lemma 3.2,

The lower bound is proved analogously, by using Eq. (20) with $H={\rm span}(v_{1},\ldots,v_{q})$.

Minimization on Grassmann manifolds and proof of Theorem 1.2

The function $F(X,Y)$ defined in Eq. (6) and to be minimized in the last part of the algorithm can naturally be viewed as defined on Grassmann manifolds. Here we recall from [EAS99] a few important facts on the geometry of Grassmann manifold and related optimization algorithms. We then prove Theorem 1.2. Technical calculations are deferred to Sections 7, 8, and to the appendices.

We recall that, for the proof of Theorem 1.2, it is assumed that $\Sigma_{\rm min}$, $\Sigma_{\rm max}$ are bounded away from and $\infty$. Numerical constants are denoted by $C,C^{\prime}$ etc. Finally, throughout this section, we use the notation $X^{(i)}\in{\mathds{R}}^{r}$ to refer to the $i$-th row of the matrix $X\in{\mathds{R}}^{m\times r}$ or $X\in{\mathds{R}}^{n\times r}$.

Denote by ${\sf O}(d)$ the orthogonal group of $d\times d$ matrices. The Grassmann manifold is defined as the quotient ${\sf G}(n,r)\simeq{\sf O}(n)/{\sf O}(r)\times{\sf O}(n-r)$. In other words, a point in the manifold is the equivalence class of an $n\times r$ orthogonal matrix $A$

For consistency with the rest of the paper, we will assume the normalization $A^{T}A=n\,{\mathbf{1}}$. To represent a point in ${\sf G}(n,r)$, we will use an explicit representative of this form. More abstractly, ${\sf G}(n,r)$ is the manifold of $r$-dimensional subspaces of ${\mathds{R}}^{n}$.

It is easy to see that $F(X,Y)$ depends on the matrices $X$, $Y$ only through their equivalence classes $[X]$, $[Y]$. We will therefore interpret it as a function defined on the manifold ${\sf M}(m,n)\equiv{\sf G}(m,r)\times{\sf G}(n,r)$:

In the following, a point in this manifold will be represented as a pair ${\bf x}=(X,Y)$, with $X$ an $n\times r$ orthogonal matrix and $Y$ an $m\times r$ orthogonal matrix. Boldface symbols will be reserved for elements of ${\sf M}(m,n)$ or of its tangent space, and we shall use ${\bf u}=(U,V)$ for the point corresponding to the matrix $M=U\Sigma V^{T}$ to be reconstructed.

Given ${\bf x}=(X,Y)\in{\sf M}(m,n)$, the tangent space at ${\bf x}$ is denoted by ${\sf T}_{{\bf x}}$ and can be identified with the vector space of matrix pairs ${\bf w}=(W,Z)$, $W\in{\mathds{R}}^{m\times r}$, $Z\in{\mathds{R}}^{n\times r}$ such that $W^{T}X=Z^{T}Y=0$. The ‘canonical’ Riemann metric on the Grassmann manifold corresponds to the usual scalar product $\langle W,W^{\prime}\rangle\equiv{\rm Tr}(W^{T}W^{\prime})$. The induced scalar product on ${\sf T}_{{\bf x}}$ between ${\bf w}=(W,Z)$ and ${\bf w}^{\prime}=(W^{\prime},Z^{\prime})$ is $\langle{\bf w},{\bf w}^{\prime}\rangle=\langle W,W^{\prime}\rangle+\langle Z,Z^{\prime}\rangle$.

This metric induces a canonical notion of distance on ${\sf M}(m,n)$ which we denote by $d({\bf x}_{1},{\bf x}_{2})$ (geodesic or arc-length distance). If ${\bf x}_{1}=(X_{1},Y_{1})$ and ${\bf x}_{2}=(X_{2},Y_{2})$ then

where the arc-length distances $d(X_{1},X_{2})$, $d(Y_{1},Y_{2})$ on the Grassmann manifold can be defined explicitly as follows. Let $\cos\theta=(\cos\theta_{1},\dots,\cos\theta_{r})$, $\theta_{i}\in[-\pi/2,\pi/2]$ be the singular values of $X_{1}^{T}X_{2}/m$. Then

The $\theta_{i}$’s are called the ‘principal angles’ between the subspaces spanned by the columns of $X_{1}$ and $X_{2}$. It is useful to introduce two equivalent notions of distance:

Notice that $d_{\rm c}$ and $d_{\rm p}$ do not depend on the specific representatives $X_{1}$, $X_{2}$, but only on the equivalence classes $[X_{1}]$ and $[X_{2}]$. Distances on ${\sf M}(m,n)$ are defined through Pythagorean theorem, e.g. $d_{\rm c}({\bf x}_{1},{\bf x}_{2})=\sqrt{d_{\rm c}(X_{1},X_{2})^{2}+d_{\rm c}(Y_{1},Y_{2})^{2}}$.

The geodesic, chordal and projection distance are equivalent, namely

For the reader’s convenience, a proof of this fact is proposed in Appendix D.

An important remark is that geodesics with respect to the canonical Riemann metric admit an explicit and efficiently computable form. Given ${\bf u}\in{\sf M}(m,n)$, ${\bf w}\in{\sf T}_{{\bf u}}$ the corresponding geodesic is a curve $t\mapsto{\bf x}(t)$, with ${\bf x}(t)={\bf u}+{\bf w}t+O(t^{2})$ which minimizes arc-length. If ${\bf u}=(U,V)$ and ${\bf w}=(W,Z)$ then ${\bf x}(t)=(X(t),Y(t))$ where $X(t)$ can be expressed in terms of the singular value decomposition $W=L\Theta R^{T}$ [EAS99]:

which can be evaluated in time of order $O(nr)$. An analogous expression holds for $Y(t)$.

2 Gradient and incoherence

The gradient of $F$ at ${\bf x}$ is the vector ${\rm grad}\,F({\bf x})\in{\sf T}_{{\bf x}}$ such that, for any smooth curve $t\mapsto{\bf x}(t)\in{\sf M}(m,n)$ with ${\bf x}(t)={\bf x}+{\bf w}\,t+O(t^{2})$, one has

In order to write an explicit representation of the gradient of our cost function $F$, it is convenient to introduce the projector operator

The two components of the gradient are then

where $Q_{X},Q_{Y}\in{\mathds{R}}^{r\times r}$ are determined by the condition ${\rm grad}\,F({\bf x})\in{\sf T}_{{\bf x}}$. This yields

3 Algorithm

At this point the gradient descent algorithm is fully specified. It takes as input the factors of ${\sf T}_{r}(\widetilde{M}^{E})$, to be denoted as ${\bf x}_{0}=(X_{0},Y_{0})$, and minimizes a regularized cost function

where $X^{(i)}$ denotes the $i$-th row of $X$, and $Y^{(j)}$ the $j$-th row of $Y$. The role of the regularization is to force ${\bf x}$ to remain incoherent during the execution of the algorithm.

We will take $\rho=n\epsilon$. Notice that $G(X,Y)$ is again naturally defined on the Grassmann manifold, i.e. $G(X,Y)=G(XQ,YQ^{\prime})$ for any $Q,Q^{\prime}\in{\sf O}(r)$.

We have $G(X,Y)=0$ on ${\cal K}(3\mu_{0})$. Notice that ${\bf u}\in{\cal K}(\mu_{0})$ by the incoherence property. Also, by the following remark proved in Appendix D, we can assume that ${\bf x}_{0}\in{\cal K}(3\mu_{0})$.

Let $U,X\in{\mathds{R}}^{n\times r}$ with $U^{T}U=X^{T}X=n{\mathbf{1}}$ and $U\in{\cal K}(\mu_{0})$ and $d(X,U)\leq\delta\leq\frac{1}{16}$. Then there exists $X^{\prime\prime}\in{\mathds{R}}^{n\times r}$ such that $X^{\prime\prime T}X^{\prime\prime}=n{\mathbf{1}}$, $X^{\prime\prime}\in{\cal K}(3\mu_{0})$ and $d(X^{\prime\prime},U)\leq 4\delta$. Further, such an $X^{\prime\prime}$ can be computed in a time of $O(nr^{2})$.

In the above, $\gamma$ must be set in such a way that $d({\bf u},{\bf x}_{0})\leq\gamma$. The next remark determines the correct scale.

Let $U,X\in{\mathds{R}}^{m\times r}$ with $U^{T}U=X^{T}X=m{\mathbf{1}}$, $V,Y\in{\mathds{R}}^{n\times r}$ with $V^{T}V=Y^{T}Y=n{\mathbf{1}}$, and $M=U\Sigma V^{T}$, $\widehat{M}=XSY^{T}$ for $\Sigma={\rm diag}(\Sigma_{1},\dots,\Sigma_{r})$ and $S\in{\mathds{R}}^{r\times r}$. If $\Sigma_{1},\dots,\Sigma_{r}\geq\Sigma_{\rm min}$, then

As a consequence of this remark and Theorem 1.1, we can assume that $d({\bf u},{\bf x}_{0})\leq C(\frac{\Sigma_{\rm max}}{\Sigma_{\rm min}})\,\frac{\mu_{1}r\sqrt{\alpha}}{\sqrt{\epsilon}}$. We shall then set $\gamma=C^{\prime}(\frac{\Sigma_{\rm max}}{\Sigma_{\rm min}})\,\frac{\mu_{1}r\sqrt{\alpha}}{\sqrt{\epsilon}}$ (the value of $C^{\prime}$ is set in the course of the proof).

Before passing to the proof of Theorem 1.2, it is worth discussing a few important points concerning the gradient descent algorithm.

The appropriate choice of $\gamma$ might seem to pose a difficulty. In reality, this parameter is introduced only to simplify the proof. We will see that the constraint $d({\bf x}_{k}(t),{\bf x}_{0})\leq\gamma$ is, with high probability, never saturated.

Indeed, the line minimization instruction 5 (which might appear complex to implement) can be replaced by a standard step selection procedure, such as the one in [Arm66].

Similarly, there is no need to know the actual value of $\mu_{0}$ in the regularization term. One can start with $\mu_{0}=1$ and then repeat the optimization doubling it at each step.

The Hessian of $F$ can be computed explicitly as well. This opens the way to quadratically convergent minimization algorithms (e.g. the Newton method).

4 Proof of Theorem 1.2

The proof of Theorem 1.2 breaks down in two lemmas. The first one implies that, in a sufficiently small neighborhood of ${\bf u}$, the function ${\bf x}\mapsto F({\bf x})$ is well approximated by a parabola.

There exists numerical constants $C_{0},C_{1},C_{2}$ such that the following happens. Assume $\epsilon\geq C_{0}\mu_{0}\sqrt{\alpha}\,r\max\{\log n;\mu_{0}r\sqrt{\alpha}(\Sigma_{\rm max}/\Sigma_{\rm min})^{4}\}$ and $\delta\leq\Sigma_{min}/C_{0}\Sigma_{\rm max}$. Then

for all ${\bf x}\in{\sf M}(m,n)\cap{\cal K}(4\mu_{0})$ such that $d({\bf x},{\bf u})\leq\delta$, with probability at least $1-1/n^{4}$. Here $S\in{\mathds{R}}^{r\times r}$ is the matrix realizing the minimum in Eq. (6).

The second Lemma implies that ${\bf x}\mapsto F({\bf x})$ does not have any other stationary point (apart from ${\bf u}$) within such a neighborhood.

There exists numerical constants $C_{0},C$ such that the following happens. Assume $\epsilon\geq C_{0}\mu_{0}r\sqrt{\alpha}(\Sigma_{\rm max}/\Sigma_{\rm min})^{2}\max\{\log n;\mu_{0}r\sqrt{\alpha}(\Sigma_{\rm max}/\Sigma_{\rm min})^{4}\}$ and $\delta\leq\Sigma_{min}/C_{0}\Sigma_{\rm max}$. Then

for all ${\bf x}\in{\sf M}(m,n)\cap{\cal K}(4\mu_{0})$ such that $d({\bf x},{\bf u})\leq\delta$, with probability at least $1-1/n^{4}$.

(Theorem 1.2) Let $\delta>0$ be such that Lemma 6.4 and Lemma 6.5 are verified, and $C_{1}$, $C_{2}$ be defined as in Lemma 6.4. We further assume $\delta\leq\sqrt{(e^{1/9}-1)/C_{2}}$. Take $\epsilon$ large enough such that, $d({\bf u},{\bf x}_{0})\leq\min(1,(C_{1}/C_{2})^{1/2}(\Sigma_{\rm min}/\Sigma_{\rm max}))\delta/10$. Further, set the algorithm parameter to $\gamma=\delta/4$.

${\bf x}_{k}\in{\cal K}(4\mu_{0})$ for all $k$.

Indeed ${\bf x}_{0}\in{\cal K}(3\mu_{0})$ whence $\widetilde{F}({\bf x}_{0})=F({\bf x}_{0})\leq C_{2}\sqrt{\alpha}n\epsilon\Sigma_{\rm max}^{2}\,\delta^{2}$. The claim follows because $\widetilde{F}({\bf x}_{k})$ is non-increasing and $\widetilde{F}({\bf x})\geq\rho\,G(X,Y)\geq n\epsilon\sqrt{\alpha}\Sigma_{\rm max}^{2}(e^{1/9}-1)$ for ${\bf x}\not\in{\cal K}(4\mu_{0})$, where we choose $\rho$ to be $n\epsilon\sqrt{\alpha}\Sigma_{\rm max}^{2}$.

$d({\bf x}_{k},{\bf u})\leq\delta/10$ for all $k$.

Since we set $\gamma=\delta/4$, by triangular inequality, we can assume to have $d({\bf x}_{k},{\bf u})\leq\delta/2$. Since $d({\bf x}_{0},{\bf u})^{2}\leq(C_{1}\Sigma_{\rm min}^{2}/C_{2}\Sigma_{\rm max}^{2})(\delta/10)^{2}$, we have $\widetilde{F}({\bf x})\geq F({\bf x})\geq F({\bf x}_{0})$ for all ${\bf x}$ such that $d({\bf x},{\bf u})\in[\delta/10,\delta]$. Since $\widetilde{F}({\bf x}_{k})$ is non-increasing and $\widetilde{F}({\bf x}_{0})=F({\bf x}_{0})$, the claim follows.

Notice that, by the last observation, the constraint $d({\bf x}_{k}(t),{\bf x}_{0})\leq\gamma$ is never saturated, and therefore our procedure is just gradient descent with exact line search. Therefore by [Arm66] this must converge to the unique stationary point of $\widetilde{F}$ in ${\cal K}(4\mu_{0})\cap\{{\bf x}:\,d({\bf x},{\bf u})\leq\delta/10\}$, which, by Lemma 6.5, is ${\bf u}$. ∎

Proof of Lemma 6.4

The following Lemma will be used several times in the following.

There exist two numerical constants $C_{1},C_{2}$ suct that the following happens. If $\epsilon\geq C_{1}\log n$ then, with probability larger than $1-1/n^{5}$,

for all $x\in{\mathds{R}}^{m}$, $y\in{\mathds{R}}^{n}$.

Write $x_{i}=x_{0}+x^{\prime}_{i}$ where $\sum_{i}x_{i}^{\prime}=0$. Then

where we recall that $\deg(j)=\{i\in[m]:\,$ such that $(i,j)\in E\}$. Further $|x_{0}|=|\sum_{i}x_{i}/m|\leq||x||_{1}/m$. The first term is upper bounded by

For $\epsilon\geq C_{1}\log n$, with probability larger than $1-1/2n^{5}$, the maximum degree is bounded by $(9/C_{1})\sqrt{\alpha}\epsilon$ which is of same order as the average degree. Therefore this term is at most $C_{2}\sqrt{\alpha}\epsilon||x||_{1}||y||_{1}/m$.

The second term is upper bounded by $C_{2}\sqrt{\alpha\epsilon}||x^{\prime}||_{2}||y||_{2}$ using Theorem 1.1 in [FO05] or, equivalently, Theorem 3.1 in the case $r=1$ and $M_{\rm max}=1$. It can be shown to hold with probability larger than $1-1/2n^{5}$ with a large enough numerical constant $C_{2}$. The thesis follows because $||x^{\prime}||_{2}\leq||x||_{2}$. ∎

2 Preliminary facts and estimates

This subsection contains some remarks that will be useful in the proof of Lemma 6.5 as well.

Let ${\bf w}=(W,Z)\in{\sf T}_{{\bf u}}$, and $t\mapsto(X(t),Y(t))$ be the geodesic such that $(X(t),Y(t))=(U,V)+(W,Z)t+O(t^{2})$. By setting $(X,Y)=(X(1),Y(1))$, we establish a one-to-one correspondence between the points ${\bf x}$ as in the statement and a neighborhood of the origin in ${\sf T}_{{\bf u}}$. If we let $W=L\Theta R^{T}$ be the singular value decomposition of $W$ (with $L^{T}L=m{\mathbf{1}}$ and $R^{T}R={\mathbf{1}}$), the explicit expression for geodesics in Eq. (29) yields

An analogous expression can obviously be written for $Y=V+\overline{Z}$. Notice that, by the equivalence between chordal and canonical distance, Remark 6.1, we have

If ${\bf u}\in{\cal K}(\mu_{0})$ and ${\bf x}\in{\cal K}(4\mu_{0})$, then $(\overline{W},\overline{Z})\in{\cal K}(10\mu_{0})$ and ${\bf w}=(W,Z)\in{\cal K}(5\pi^{2}\mu_{0}/2)$.

The first fact follows from $||\overline{W}^{(i)}||^{2}\leq 2||X^{(i)}||^{2}+2||U^{(i)}||^{2}$. In order to prove ${\bf w}\in{\cal K}(5\pi^{2}\mu_{0}/2)$, we notice that

The claim follows by showing a similar bound for $||Z^{(i)}||^{2}$. ∎

We next prove a simple a priori estimate.

There exist numerical constants $C_{1},C_{2}$ such that the following holds with probability larger than $1-1/n^{5}$. If $\epsilon\geq C_{1}\log n$, then for any $(X,Y)\in{\cal K}(\mu)$ and $S\in{\mathds{R}}^{r\times r}$,