A Gradient Descent Algorithm on the Grassman Manifold for Matrix Completion

Introduction

In this paper we consider the problem of reconstructing an $m\times n$ low rank matrix $M$ from a small set of observed entries. This problem is of considerable practical interest and has many applications. One example is collaborative filtering, where users submit rankings for small subsets of, say, movies, and the goal is to infer the preference of unrated movies for a recommendation system . It is believed that the movie-rating matrix is approximately low-rank, since only a few factors contribute to a user’s preferences. Other examples of matrix completion include the problem of inferring 3-dimensional structure from motion and triangulation from incomplete data of distances between wireless sensors, also known as the sensor localization problem , .

On the theoretical side, most recent work focuses on algorithms for exactly recovering the unknown low-rank matrix. They provide an upper bound on the number of observed entries that guarantee successful recovery with high probability. The main assumptions of this exact matrix completion problem are that the matrix $M$ to be recovered has rank $r\ll m,n$ and that the observed entries are known exactly. The problem is equivalent to finding a minimum rank matrix matching the observed entries . This problem is NP-hard. Adapting techniques from compressed sensing, Candès and Recht introduced a convex relaxation of this problem . They introduced the concept of incoherence and proved that for a matrix $M$ of rank $r$ which has the incoherence property, solving the convex relaxation correctly recovers the unknown matrix, with high probability, if the number of observed entries $|E|$ satisfies, $|E|\geq C(\alpha)rn^{1.2}\log n$.

Recently improved the bound to $|E|\geq C(\alpha)rn\max\{\log n,r\}$ for matrices with bounded condition number and provided an efficient algorithm called OptSpace , based on spectral methods followed by local manifold optimization. For a bounded rank $r$, it is order optimal in the sense that an $m\times n$ rank-$r$ matrix $M$ has $r(m+n-r)$ degrees of freedom and without the same number of observations it is impossible to fix them. The extra $\log n$ factor is due to a coupon-collector effect: it is necessary that $E$ contains at least one entry per row and one per column, which happens only for $|E|\geq Cn\log n$ . Candès and Tao proved a similar bound $|E|\geq C(\alpha)nr(\log n)^{6}$ with a stronger assumption on the original matrix $M$, known as the strong incoherence condition but without any assumption on the condition number of $M$. For any value of $r$, it is only suboptimal by a poly-logarithmic factor In , which appeared while we were preparing this manuscript, improved guarantees were proved for the convex relaxation algorithm. Namely, assuming only the incoherence property on the original matrix $M$, the convex relaxation correctly recovers the matrix $M$ if $|E|\geq C(\alpha)nr(\log n)^{2}$..

When the pattern of observed entries is non-random, it is interesting to ask if the solution of the rank $r$ matrix completion problem is unique . Building on the ideas from rigidity theory, Singer and Cucuringu introduce a randomized algorithm to determine whether it is possible to uniquely complete a partially observed matrix to a matrix of specific rank $r$. Furthermore, by applying their algorithm to random patterns of observed entries, one can get a lower bound on the minimum number of observed entries necessary to correctly recover the matrix $M$.

While most theoretical work focuses on proving bounds for the exact matrix completion problem, a more interesting and practical problem is when the matrix M is only approximately low rank or when the observation is corrupted by noise. The main focus of this approximate matrix completion problem is to design an algorithm to find an $m\times n$ low-rank matrix $\widehat{M}$ that best approximates the original matrix $M$ and provide a bound on the root mean squared error (RMSE) given by, ${\rm RMSE}=\frac{1}{\sqrt{mn}}||M-\widehat{M}||_{F}$. Candès and Plan introduced a generalization of the convex relaxation from to the approximate case, and provided a bound on the RMSE . More recently, a bound on the RMSE achieved by the OptSpace algorithm with noisy observations was obtained in . This bound is order optimal in a number of situations and improves over the analogous result in .

On the practical side, directly solving the convex relaxation introduced in requires solving a Semidefinite Program (SDP), the complexity of which grows proportional to $n^{3}$. In the last year, many authors have proposed efficient algorithms for solving the low-rank matrix completion problem. These include Singular Value Thresholding (SVT) , Accelerated Proximal Gradient (APG) algorithm , Fixed Point Continuation with Approximate SVD (FPCA) , Atomic Decomposition for Minimum Rank Approximation (ADMiRA) , Soft-Impute , Subspace Evolution and Transfer (SET) , Singular Value Projection (SVP) and OptSpace .

The problems each of these algorithms are trying to solve are described in Section 2. SVT is an iterative algorithm for solving the convex relaxation of the exact matrix completion problem, which minimizes the nuclear norm (a sum of the singular values) under the constraints of matching the observed entries as in (9). APG, FPCA and Soft-Impute are efficient algorithms for solving the convex relaxation of the approximate matrix completion problem, which is a nuclear norm regularized least squares problem as in (11). ADMiRA is an extension of Compressive Sampling Matching Pursuit (CoSaMP) , which is an iterative method for solving a least squares problem with bounded rank $r$ as described in (12). SVP is another approach to solve (12), which is a generalization of Iterative Hard Thresholding (IHT) .

2 Contributions and outline

The main contribution of this paper is to develop and implement efficient procedures, based on the OptSpace algorithm introduced in , for solving the exact and approximate matrix completion problems and add novel modifications, namely Rank Estimation and Incremental OptSpace, that allow for a broader application and a better performance.

The algorithm described in requires a knowledge of the rank of the original matrix. In the following, we introduce a procedure, Rank Estimation, that is guaranteed to correctly estimate the rank of the original matrix from a partially revealed matrix under some conditions, which in turn allows us to use the algorithm in a broader set of applications.

Next, we introduce Incremental OptSpace, a novel modification to OptSpace. We show that, empirically, Incremental OptSpace has substantially better performance than OptSpace when the underlying matrix is ill-conditioned.

Further, we carry out an extensive empirical comparison of various reconstruction algorithms. This is particularly important because performance guarantees are only “up to constants” and therefore they have limited use in comparing different algorithms. Finally, we apply our algorithm to real-world data and demonstrate that it is readily applicable to the real data.

The organization of the paper is as follows. In Section 2 we describe the low rank matrix completion problem and convex relaxations to the basic NP-hard approach, mostly to set our notation for later use. Section 3 introduces an efficient implementation of the OptSpace algorithm with novel modifications. In Section 4 we discuss the results of numerical simulations with respect to speed and accuracy and compare the performance of OptSpace with that of the other algorithms.

The model definition

In the case of exact matrix completion, we assume that the matrix $M$ has exact low rank $r\ll\min\{m,n\}$ and that the observed entries are known exactly. More precisely, we assume that there exist matrices $U$ of dimensions $m\times r$, V of dimensions $n\times r$, and a diagonal matrix $\Sigma$ of dimensions $r\times r$, such that

Notice that for a given matrix $M$, the factors $(U,V,\Sigma)$ are not unique.

Out of the $m\times n$ entries of $M$, a subset $E\subseteq[m]\times[n]$ is observed. Let $M^{E}$ be the $m\times n$ observed matrix that stores all the observed values, such that

Our goal is to find a low rank estimation $\widehat{M}(M^{E},E)$ of the original matrix $M$ from the observed matrix $M^{E}$ and the set of observed indices $E$.

If the number of observed entries $|E|$ is large enough, there is a unique rank $r$ matrix which matches the observed entries. In this case, solving the following optimization problem will recover the original matrix correctly.

where $X\in{\mathds{R}}^{m\times n}$ is the variable matrix, ${\rm rank}\,(X)$ is the rank of matrix $X$, and ${\cal P}_{E}(\cdot)$ is the projector operator defined as

The solution of this problem is the lowest rank matrix that matches the observed entries. Notice that this is optimal in the sense that if this problem does not recover the correct matrix $M$ then there exists at least one other rank-$r$ matrix that matches all the observations, and no other algorithm can do better. However, this optimization problem is NP-hard and all known algorithms require doubly exponential time in $n$ . This is especially inadequate since we are interested in cases where the dimension of the matrix $M$ is large ( eg. such as $5\cdot 10^{5}\times 2\cdot 10^{4}$ for ).

where $||X||_{*}$ denotes the nuclear norm of $X$, i.e the sum of its singular values.

In the case of approximate matrix completion problem, where the observations are contaminated by noise or the original matrix to be reconstructed is only approximately low rank, the constraint ${\cal P}_{E}(X)={\cal P}_{E}(M)$ must be relaxed. This results in either the problem

In , problem (5) is recast into the rank-$r$ matrix approximation problem of

In the following, we present an efficient algorithm, namely OptSpace, to solve the low-rank matrix completion problem which is closely related to (12), and numerically compare its performance with those of the competing algorithms in the case of exact as well as approximate matrix completion problems.

Algorithm

Algorithm 1 describes the overview of OptSpace . Each step is explained in detail in the following sections. The basic idea is to minimize the cost function $F:{\mathds{R}}^{m\times r}\times{\mathds{R}}^{n\times r}\to{\mathds{R}}$, defined as

Minimizing $F(X,Y)$ is an a priori difficult task, since $F$ is a non-convex function. The basic idea is that the singular value decomposition (SVD) of $M^{E}$ provides an excellent initial guess, and that the minimum can be found with high probability by standard gradient descent after this initialization. Two caveats must be added to this description: $(1)$ In general the matrix $M^{E}$ must be ‘trimmed’ to eliminate over-represented rows and columns; $(2)$ We need to estimate the target rank $r$.

We say that a row is over-represented if its degree is more than $2|E|/m$ (twice the average degree), where degree of a row is defined as the number of observed entries in that row. Analogously, a column is over-represented if its degree is more than $2|E|/n$. Trimming is a procedure that takes $M^{E}$ and $E$ as input and outputs $\widetilde{M}^{E}$ by setting to $0$ all of the entries in over-represented rows and columns. Let $d_{l}(i)$ and $d_{r}(j)$ be the degree of $i^{th}$ row and $j^{th}$ column of $M$ respectively. Then the trimmed matrix $\widetilde{M}^{E}$ is defined as

The trimming step is essential when $|E|=\Theta(n)$, in which case there exists over-represented columns and rows of degrees $\Theta(\log n/\log\log n)$, corresponding to singular values of the order $\Theta(\sqrt{\log n/\log\log n})$. As $n$ grows large, while these spurious singular values dominate the principal components in step 3 of the Algorithm 1, the corresponding singular vectors are highly concentrated on the over-represented rows and columns (respectively for left and right singular vectors) and do not provide any useful information about the unobserved entries of $M$.

2 Estimating the rank

Define $\epsilon\equiv|E|/\sqrt{mn}$. In the case of a square matrix $M$, $\epsilon$ corresponds to the average degree per row or per column. Throughout this paper, the parameter $\epsilon$ will be frequently used as the model parameter indicating how difficult the problem instance is.

By singular value decomposition of the trimmed matrix, let

where $x_{i}$ and $y_{i}$ are the left and right singular vectors corresponding to $i$th singular value $\sigma_{i}$. Then, the following cost function is defined in terms of the singular values.

Based on the above definition, Rank Estimation consists of two steps:

The idea behind this algorithm is that, if enough entries of $M$ are revealed then there is a clear separation between the first $r$ singular values, which reveal the structure of the matrix $M$ to be reconstructed, and the spurious ones . As described in the following proposition, we can show that this simple procedure is guaranteed to reconstruct the correct rank $r$, with high probability, for $|E|$ large enough. For the proof of this proposition, we refer to Appendix A.

Assume $M$ to be a rank $r$ $m\times n$ matrix with bounded condition number $\kappa$. Then there exists a constant $C(\kappa)$ such that, if $\epsilon>C(\kappa)r$, then Rank Estimation correctly estimates the rank $r$, with high probability.

3 Rank-ρ\rho projection

Rank-$\rho$ projection consists of performing a sparse SVD on $\widetilde{M}^{E}$ and rescaling the singular values and singular vectors appropriately. From the Rank Estimation step we have the SVD of $\widetilde{M}^{E}$ in Eq. (18), namely $\widetilde{M}^{E}=\sum_{i=1}^{\min(m,n)}\sigma_{i}x_{i}y_{i}^{T}$. Define the projection :

for normalized orthogonal matrices $X_{0}\in{\mathds{R}}^{m\times\rho}$ and $Y_{0}\in{\mathds{R}}^{n\times\rho}$, and a $\rho\times\rho$ diagonal matrix $S_{0}$, defined in terms of the singular values and singular vectors in Eq. (18) as $X_{0}=\sqrt{m}[x_{1},\ldots,x_{\rho}]$, $Y_{0}=\sqrt{n}[y_{1},\ldots,y_{\rho}]$, and $S_{0}=({1}/{\epsilon})\text{diag}(\sigma_{1},\ldots,\sigma_{\rho})$. Notice that we do not need to compute the scaled singular values $S_{0}$, since we only require $X_{0}$ and $Y_{0}$ for the following local optimization step. There are a number of low complexity algorithms available for forming a sparse SVD, as well as a number of open source implementations of these algorithms.

4 Gradient descent on the Grassman manifold

The Manifold Optimization step involves gradient descent with variables $X\in{\mathds{R}}^{m\times r}$ and $Y\in{\mathds{R}}^{n\times r}$ using the cost function $F(X,Y)$ defined below. In this section, we use $r$ and ${\hat{r}}$ interchangeably to denote the estimated rank of matrix $M$.

where $f:{\mathds{R}}\times{\mathds{R}}\to{\mathds{R}}$ is an element-wise cost function. Note that compared to Eq. (14), we have additional term in Eq. (22), which is a regularization term with a regularization coefficient $\lambda\in$.

The above general formulation allows for a freedom in choosing a suitable cost function $f$ for different applications. However, a common example of the cost function $f(x,y)=\frac{1}{2}(x-y)^{2}$ works very well in practice as well as in proving performance bounds . Hence, throughout this paper, we use the squared difference as the cost function, resulting in

where the projector operator ${\cal P}_{E}$ for a given $E$ is defined in Eq. (2), and $E^{\perp}$ is the complementary set of $E$.

For the results in this paper, we choose $\lambda=0$ but we observe that using a positive $\lambda$ helps when the matrix entries are corrupted by noise. For $\lambda=0$, the gradient of $F(X,Y)$ can be written explicitly as

where $S$ is the $r\times r$ matrix that achieves the minimum in the definition of $F(X,Y)$, Eq. (21).

One important feature of OptSpace is that $F(X,Y)$ is regarded as a function of the $r$-dimensional subspaces of ${\mathds{R}}^{m}$ and ${\mathds{R}}^{n}$ generated (respectively) by the columns of $X$ and $Y$. This interpretation is justified by the fact that $F(X,Y)=F(XA,YB)$ for any two orthogonal matrices $A$, $B\in{\mathds{R}}^{r\times r}$. The set of $r$ dimensional subspaces of ${\mathds{R}}^{m}$ is a differentiable Riemannian manifold ${\sf G}(m,r)$ (the Grassman manifold) . The gradient descent algorithm is applied to the function $F:{\sf G}(m,r)\times{\sf G}(n,r)\to{\mathds{R}}\;$ For further details on optimization by gradient descent on matrix manifolds we refer to .

In the following, we use a compact representation ${\bf x}$ for a pair $(X,Y)$, with $X$ an $n\times r$ matrix and $Y$ an $m\times r$ matrix. Similarly, the gradient is represented by : ${\rm grad}\,F({\bf x}_{k})=({\rm grad}\,F({\bf x}_{k})_{X},{\rm grad}\,F({\bf x}_{k})_{Y})$. Let ${\bf x}_{0}=(X_{0},Y_{0})$, where $X_{0}$ and $Y_{0}$ are the normalized left and right singular matrices from rank-$r$ projection. The Manifold Optimization algorithm starting at ${\bf x}_{0}$ is described below. We refer to for justifications and performance bounds of the algorithm.

For any scalar $\tau$, it is shown in that this algorithm converges to the local minimum. However, numerical experiments suggest $\tau=10^{-3}$ is a good choice. The algorithm stops when the fit error $||{\cal P}_{E}(M-\widehat{M})||_{F}/||{\cal P}_{E}(M)||_{F}$ goes below some threshold $\delta_{\rm tol}$, e.g. $10^{-6}$. The basic idea is that this is a good indicator of the relative error on the whole set, $||M-\widehat{M}||_{F}/||M||_{F}$. This stopping criterion is also used in other algorithms such as SVT in where the authors provide a convincing argument for its use.

5 A novel modification to OptSpace for ill-conditioned matrices

In this section, we describe a novel modification to the OptSpace algorithm, which has substantially better performance in the case when the matrix $M$ to be reconstructed is ill-conditioned. When the condition number $\kappa(M)$ is high, the initial guess in step 3 of OptSpace for ($u_{r},v_{r}$), the singular vectors which correspond to the smallest singular value, are often far from the correct ones. To compensate for this discrepancy, we start by first finding ($u_{1},v_{1}$), the singular vectors corresponding to the first singular value, and incrementally search for the next ones.

In the following numerical simulations, we demonstrate that Incremental OptSpace brings significant performance gains when applied to ill-conditioned matrices, cf. Section 4.

Numerical results with randomly generated matrices

The OptSpace algorithm described above was implemented in C The code is available at http://www.stanford.edu/$\sim$raghuram/optspace.html and tested on a 3.4 GHz Desktop computer with 4 GB RAM. For efficient singular value decomposition of sparse matrices, we used (a modification of) SVDLIBC Available at http://tedlab.mit.edu/$\sim$dr/svdlibc/ which is based on SVDPACKC. In this section, we compare the performance of OptSpace with other algorithms by numerical simulations. In Section 4.1, the algorithms are tested on randomly generated matrices with noiseless observations, and in Section 4.2 we compare the algorithms when we have noisy observations under different scenarios.

For exact matrix completion experiments, we use $n\times n$ test matrices $M$ of rank $r$ generated as $M=UV^{T}$. Here, $U$ and $V$ are $n\times r$ matrices with each entry being sampled independently from a standard Gaussian distribution ${\cal N}(0,1)$, unless specified otherwise. Then, each entry is revealed independently with probability $\epsilon/n$, so that on an average $n\epsilon$ entries are revealed. Numerical results show that there is no notable difference if we choose the revealed set of entries $E$ uniformly at random over all the subsets of the same size $|E|=n\epsilon$. We use $\delta_{\rm tol}=10^{-5}$ and $k_{max}=1000$ as the stopping criteria.

For approximate matrix completion, the matrices are generated as above and corrupted by additive noise $Z_{ij}$. First, in the standard scenario, $Z_{ij}$’s are independently and identically distributed according to a Gaussian distribution. For comparison, we also present numerical simulation results with different types of noise in the following subsections. Again, each entry is revealed independently with a probability $\epsilon/n$. We use $||{\cal P}_{E}(\widehat{M}-(M+Z))||_{F}^{2}\leq(1+\epsilon)|E|\sigma_{n}^{2}$ (where $\sigma_{n}^{2}$ is the noise variance per entry) as the stopping criterion.

We first illustrate the rate of convergence of OptSpace . In Figure 1, we plot the fit error, $||{\cal P}_{E}(\widehat{M}-M)||_{F}/n$ and the prediction error $||\widehat{M}-M||_{F}/n$, with respect to the number of iterations of the Manifold Optimization step. These plots are obtained for matrices with $n=1000$ and $r=10$ and averaged over $10$ instances. The results are shown for two values of $\epsilon$: $100$ and $200$. We can see that the prediction error decays exponentially with the number of iterations in both cases. Also, the prediction error is very close to the fit error, thus lending support to the validity of the chosen stopping criterion.

We next study the reconstruction rate of the algorithm. We declare a matrix to be reconstructed if $||M-\widehat{M}||_{F}/||M||_{F}\leq 10^{-4}$. The reconstruction rate is the fraction of instances for which the matrix was reconstructed.

In Figure 2, we plot the reconstruction rate as function of $|E|/n$ for OptSpace on randomly generated rank-$4$ matrices for different matrix sizes $n$. As predicted by Theorem 1.2 of , threshold of the reconstruction rate of OptSpace is upper bounded by $|E|=Cn(\log n)^{2}$, for fixed rank $r=4$. Here, an extra factor of $\log n$ comes from the fact that if we generate random factors $U$ and $V$ from a Gaussian distribution, then the incoherence parameter $\mu_{0}$ scales like $\log n$. However, the location of the threshold is surprisingly close to the lower bound proved in which scales as $|E|=Cn\log n$. The lower bound provides a threshold below which the problem admits more than one solution. Note that the lower bound is displayed only for the case when $n=1000$.

In Figure 3, we plot the reconstruction rate for randomly generated matrices with dimensions $m=n=500$ using OptSpace . The resulting reconstruction rate is plotted for different ranks $r$ as a function of $|E|/n$. As rank increases and for fixed $n$, the reconstruction rate has a sharp threshold at $|E|=Crn\log n$. This indicates that in practice the dependence of the threshold on the rank scales like $r$ rather than $r^{2}$ as predicted by Theorem 1.2 of . Also, for all values of rank, the location of the threshold is surprisingly close to the lower bound proved in , below which the problem admits more than one solution.

In Figure 4, we plot the reconstruction rate of OptSpace as a function of $|E|/n$ for rank $10$ matrices of dimension $m=n=1000$. Also plotted are the reconstruction rates obtained for the convex relaxation approach of solved using the Singular Value Thresholding algorithm , the FPCA algorithm from and ADMiRA . We compare these with a theoretical lower bound on the reconstruction rate described in . Various algorithms exhibit threshold at different values of $|E|/n$, and the threshold depends on the problem size $n$ and the rank $r$. This figure clearly illustrates that OptSpace outperforms the other algorithms on random data, and this was consistent for various values of $n$ and $r$.

In the following Tables 1 and 2, we present numerical results obtained using these algorithms for different values of $n$ and $r$. Table 1 presents results for smaller values of $\epsilon$ and hence for hard problems, whereas Table 2 presents results for larger values of $\epsilon$ which are relatively easy problems. Note that the values of $\epsilon$ used in Table 1 all correspond to $|E|\leq 2.6d(n,r)$ where $d(n,r)=2nr-r^{2}$ is the number of degrees of freedom. We ran into Out of Memory problems for the FPCA algorithm for $n\geq 20000$ and hence we omit these problems from the table. All the results presented in tables are averaged over 5 instances. On the easy problems, all the algorithms achieved similar performances, whereas on the hard problems, OptSpace outperforms other algorithms on most of instances.

2 Approximate matrix completion

In this section we compare the performance of different algorithms for matrix completion with noisy observations. As a metric, we use the relative root mean squared error defined as

For direct comparison we start with an example taken from . In this example, $M$ is a square matrix of dimensions $n\times n$ and rank $r$ generated as $M=UV^{T}$ with fixed $n=600$. $U$ and $V$ are $n\times r$ matrices with each entry being sampled independently from a standard Gaussian distribution ${\cal N}(0,\sigma_{s}^{2}=20/\sqrt{n})$. As before, each entry is revealed independently with probability $\epsilon/n$. Each entry is corrupted by added noise matrix $Z$, so that the observation for the index $(i,j)$ is $M_{ij}+Z_{ij}$. Further, $Z$ has each entries drawn from i.i.d. standard Gaussian distribution ${\cal N}(0,1)$. In the following we refer to this noise model as the standard scenario. We also refer to for the data for the convex relaxation approach and the information theoretic lower bound.

Figure 5 compares the average root mean squared error achieved by the different algorithms for a fixed rank $r=2$ as a function of $|E|/n$. After one iteration, for most values of $\epsilon$, OptSpace has a smaller root mean square error than the convex relaxation approach and in about $10$ iterations, it becomes indistinguishable from the information theoretic lower bound. In Figure 6, we compare the average root mean squared error obtained for a fixed sample size $\epsilon=120$ as a function of the rank. Again, for most values of $r$, after one iteration OptSpace has a smaller root mean square error than the convex relaxation based algorithm.

Table 4 illustrate how the performance changes with different noise power for fixed $n=1000$. We present the results of our experiments with different ranks and noise ratios defined as

Next, in the following series of examples, we illustrate how the performances change under different noise models. In the following, $M$ is a square matrix generated as $UV^{T}$ like above, but $U$ and $V$ now have each entry sampled independently from a standard Gaussian distribution ${\cal N}(0,1)$, unless specified otherwise. As before, each entry is revealed independently with probability $\epsilon/n$ and the observation is corrupted by added noise matrix $Z$. We compare the resulting RMSE of FPCA, ADMiRA and OptSpace on this randomly generated matrices with noisy observations and missing entries. Since ADMiRA requires a target rank, we use the rank estimated using Rank Estimation described in Section 3.2. For FPCA we choose $\mu=\sqrt{2np}\sigma$, where $p=|E|/n^{2}$ and $\sigma^{2}$ is the variance of each entry in $Z$. A convincing argument for this choice of $\mu$ is given in .

In the standard scenario, we typically make the following three assumptions on the noise matrix $Z$. (1) The noise $Z_{ij}$ does not depend on the value of the matrix $M_{ij}$. (2) The entries of $Z$, $\{Z_{ij}\}$, are independent. (3) The distribution of each entries of $Z$ is Gaussian. The matrix completion algorithms described in Section 2 are expected to be especially effective under this standard scenario for the following two reasons. First, the squared error objective function that the algorithms minimize is well suited for the Gaussian noise. Second, the independence of $Z_{ij}$’s ensure that the noise matrix is almost full rank and the singular values are evenly distributed. This implies that for a given noise power $||Z||_{F}$, the spectral norm $||Z||_{2}$ is much smaller than $||Z||_{F}$. In the following, we fix $m=n=500$ and $r=4$, and study how the performance changes with different noise. Each of the simulation results is averaged over 10 instances and is shown with respect to two basic parameters, the average number of revealed entries per row $|E|/n$ and the noise ratios $N$, defined as Eq. (24).

In this standard scenario, the noise $Z_{ij}$’s are distributed as i.i.d. Gaussian N(0,$\sigma^{2}$). Note that the noise ratio is equal to $N=\sigma/2$. The accuracy of the estimation is measured using RMSE. We compare the resulting RMSE of FPCA, ADMiRA and OptSpace to the RMSE of the oracle estimate, which is $\sigma\sqrt{(2nr-r^{2})/|E|}$ .

Figure 7 shows the performance for each of the algorithms with respect to $|E|/n$ under the standard scenario for fixed $N=1/2$. For most values of $|E|$, the simple rank-$r$ projection has the worst performance. However, when all the entries are revealed and the noise is i.i.d. Gaussian, the simple rank-$r$ projection coincides with the oracle bound, which in this simulation corresponds to the value $|E|/n=500$. Note that the behavior of the performance curves of FPCA, ADMiRA, and OptSpace with respect to $|E|$ is similar to the oracle bound, which is proportional to $1/\sqrt{|E|}$.

Among the three algorithms, FPCA has the largest RMSE, and OptSpace is very close to the oracle bound for all values of $|E|$. Note that when all the values are revealed, ADMiRA is an efficient way of implementing rank-$r$ projection, and the performances are expected to be similar. This is confirmed by the observation that for $|E|/n\geq 400$ the two curves are almost identical. One of the reasons why the RMSE of FPCA does not decrease with $|E|$ for large values of $|E|$ is that FPCA overestimates the rank and returns estimated matrices with rank much higher than $r$, whereas the rank estimation algorithm used for ADMiRA and OptSpace always returned the correct rank $r$ for $|E|/n\geq 80$.

2.2 Multiplicative Gaussian noise

In sensor network localization , where the entries of the matrix corresponds to the pair-wise distances between the sensors, the observation noise is oftentimes assumed to be multiplicative. In formulae, $Z_{ij}=\xi_{ij}M_{ij}$, where $\xi_{ij}$’s are distributed as i.i.d. Gaussian with zero mean. The variance of $\xi_{ij}$’s are chosen to be $1/r$ so that the resulting noise ratio is $N=1/2$. Note that in this case, $Z_{ij}$’s are mutually dependent through $M_{ij}$’s and the values of the noise also depend on the value of the matrix entry $M_{ij}$.

Figure 9 shows the RMSE with respect to $|E|/n$ under multiplicative Gaussian noise. The RMSE of the rank-$r$ projection for $|E|/n=40$ is larger than $1.5$ and is omitted in the figure. The bottommost line corresponds to the oracle performance under standard scenario, and is displayed here, and all of the following figures, to serve as a reference for comparison. The main difference with respect to Figure 7 is that most of the performance curves are larger under multiplicative noise. For the same value of the noise ratio $N$, it is more difficult to distinguish the noise from the original matrix, since the noise is now correlated with the matrix $M$.

2.3 Outliers

In structure from motion , the entries of the matrix corresponds to the position of points of interest in $2$-dimensional images captured by cameras in different angles and locations. However, due to failures in the feature extraction algorithm, some of the observed positions are corrupted by large noise where as most of the observations are noise free. To account for such outliers, we use the following model.

The value of $a$ is chosen according to the target noise ratio $N=a/20$. The noise is independent of the matrix entries and $Z_{ij}$’s are mutually independent, but the distribution is now non-Gaussian.

Figure 10 shows the performance of the algorithms with respect to $|E|/n$ and the noise ratio $N$ with outliers. Comparing the first figure to Figure 7, we can see that the performance for large value of $|E|$ is less affected by outliers compared to the small values of $|E|$. The second figure clearly shows how the performance degrades for non-Gaussian noise when the number of samples is small. The algorithms minimize the squared error $||{\cal P}_{E}(X)-{\cal P}_{E}(N)||_{F}^{2}$ as in (11) and (12). For outliers, a suitable algorithm would be to minimize the $l_{1}$-norm of the errors instead of the $l_{2}$-norm . Hence, for this simulation with outliers, we can see that the performance of the rank-$r$ projection, ADMiRA and OptSpace is worse than the Gaussian noise case. However, the performance of FPCA is almost the same as in the standard scenario.

2.4 Quantization noise

One common model for noise is the quantization noise. For a regular quantization, we choose a parameter $a$ and quantize the matrix entries to the nearest value in $\{\ldots$, $-a/2$, $a/2$, $3a/2$, $5a/2$, $\ldots\}$. The parameter $a$ is chosen carefully such that the resulting noise ratio is $1/2$. The performance for this quantization is expected to be worse than the multiplicative noise case. The reason is that now the noise is deterministic and completely depends on the matrix entries $M_{ij}$, whereas in the multiplicative noise model it was random.

Figure 11 shows the performance against $|E|/n$ within quantization noise. The overall behavior of the performance curves is similar to Figure 7, but most of the curves are shifted up. Note that the bottommost line is the oracle performance in the standard scenario which is the same in all the figures. Compared to Figure 9, for the same value of $N=1/2$, quantization is much more detrimental than the multiplicative noise as expected.

2.5 Ill conditioned matrices

Numerical results with real data matrices

In this section, we consider low-rank matrix completion problems in the context of recommender systems, based on two real data sets: the Jester joke data set and the Movielens data set . The Jester joke data set contains $4.1\times 10^{6}$ ratings for 100 jokes from 73,421 users. The dataset is available at http://www.ieor.berkeley.edu/$\sim$goldberg/jester-data/ Since the number of users is large compared to the number of jokes, we randomly select $n_{u}\in\{100,1000,2000,4000\}$ users for comparison purposes. As in , we randomly choose two ratings for each user as a test set, and this test set, which we denote by $T$, is used in computing the prediction error in Normalized Mean Absolute Error (NMAE). The Mean Absolute Error (MAE) is defined as in .

where $M_{ij}$ is the original rating in the data set and $\widehat{M}_{ij}$ is the predicted rating for user $i$ and item $j$. The Normalized Mean Absolute Error (NMAE) is defined as

where $M_{\rm max}$ and $M_{\rm min}$ are upper and lower bounds for the ratings. In the case of Jester joke, all the ratings are in $$which implies that$M_{\rm max}=10$and$M_{\rm min}=-10$.

The numerical results for Jester joke data set using Incremental OptSpace, FPCA and ADMiRA are presented in the first four columns of Table 5. In the table, $rank$ indicates the rank used to estimate the matrix and $time$ is the running time of each matrix completion algorithm. To get an idea of how good the predictions are, consider the case where each missing entries is predicted with a random number drawn uniformly at random in $$and the actual rating is also a random number with same distribution. After a simple computation, we can see that the resulting NMAE of the random prediction is 0.333. As another comparison, for the same data set with$n_{u}=18000$, simple nearest neighbor algorithm and Eigentaste both yield NMAE of 0.187 . The NMAE of Incremental OptSpace is lower than these simple algorithms even for$n_{u}=100$and tends to decrease with$n_{u}$.

Numerical simulation results on the Movielens data set is also shown in the last row of the above table. The data set contains $100,000$ ratings for $1,682$ movies from $942$ users. The dataset is available at http://www.grouplens.org/node/73 We use $80,000$ randomly chosen ratings to estimate the $20,000$ ratings in the test set, which is called $u1.base$ and $u1.test$, respectively, in the movielens data set. In the last column of Table 5, we compare the resulting NMAE using Incremental OptSpace , FPCA and ADMiRA.

Next, to get some insight on the structure of real data, we look at a complete sub matrix where all the entries are known. With Jester joke data set, we deleted all users containing missing entries, and generated a complete matrix $M$ with $14,116$ users and $100$ jokes. The distribution of the singular values of $M$ is shown in Figure 13. We must point out that this rating matrix is not low-rank or even approximately low-rank, although it is common to make such assumptions. This is one of the difficulties in dealing with real data. The other aspect is that the samples are not drawn uniformly at random as commonly assumed in .

Finally we test the incoherence assumption for the Netflix dataset in Figure 14 and Figure 15. For a $m\times n$ matrix whose singular value decomposition is given by $M=U\Sigma V^{T}$, $M$ is said to be $(\mu_{0},\mu_{1})$-incoherent if it satisfies the following properties :

There exists a constant $\mu_{0}>0$ such that for all $i\in[m]$, $j\in[n]$ we have $\sum_{k=1}^{r}{U_{i,k}^{2}}\leq\mu_{0}r/m$, $\sum_{k=1}^{r}{V_{i,k}^{2}}\leq\mu_{0}r/n$.

There exists $\mu_{1}$ such that $|\sum_{k=1}^{r}{U_{i,k}V_{j,k}}|\leq\mu_{1}\sqrt{r/mn}$.

To check if A1 holds for the Netflix movie ratings matrix, we run OptSpace on the Netflix dataset and plot cumulative sum of the sorted row norms of the left and right factors defined as follows. Let the output of OptSpace be $X\in{\mathds{R}}^{m\times r}$, $Y\in{\mathds{R}}^{n\times r}$ and $S\in{\mathds{R}}^{r\times r}$. Here $m=480,189$ is the number of users and $n=17,770$ is the number of movies. For the target rank we used $r=5$. Let $x_{i}=\frac{m}{r}||X^{(i)}||^{2}$ and $y_{i}=\frac{n}{r}||Y^{(i)}||^{2}$ where $X^{(i)}$ and $Y^{(i)}$ denote the $i$th row of the left factor $X$ and the right factor $Y$ respectively. Define a permutation $\pi_{l}:[m]\rightarrow[m]$ which sorts $x_{i}$’s in a non-decreasing order such that $x_{\pi_{l}(1)}\leq x_{\pi_{l}(2)}\leq\ldots\leq x_{\pi_{l}(m)}$. Here, we used the standard combinatorics notation $[k]=\{1,2,\ldots,k\}$ for an integer $k$. Similarly, we can define $\pi_{r}:[n]\rightarrow[n]$ for $y_{i}$’s.

In Figure 14, we plot $\sum_{i=1}^{k}x_{i}$ vs. $k$. For comparison, we also plot the corresponding results for a randomly generated matrix $X_{G}$. Generate $U\in{\mathds{R}}^{m\times r}$ by sampling its entries $U_{ij}$ independently and distributed as ${\cal N}(0,1)$ and let $X_{G}$ be the left singular vectors of $U$. Since $x_{i}$’s are scaled by $m/r$, when $k=m$ we have $\sum_{i=1}^{m}x_{i}=m$. This is also true for the random matrix $X_{G}$. Figure 15 shows the corresponding plots for $Y$. For a given matrix, if A1 holds with a small $\mu_{0}$ then the corresponding curve would be close to a straight line. The curvature in the plots is indicative of the disparity among the row weigths of the factors. We can see that a randomly generated matrix would satisfy A1 with a smaller value of $\mu_{0}$ compared to the movie ratings matrix, hence can be said to be more incoherent. The factor corresponding to movies has a larger disparity than the factor corresponding to users, and hence challenges the incoherence assumption.

Appendix A Proof of Proposition 3.1

The matrix $M$ to be reconstructed is factorized as Eq. (1), where $\Sigma={\rm diag}(\Sigma_{1},\ldots,\Sigma_{r})$ is a diagonal matrix of the singular values. We start from following key lemma.

There exists a numerical constant $C>0$ such that, with high probability

where it is understood that $\Sigma_{q}=0$ for $q>r$, and ${M_{\rm max}}=\max\{M_{ij}\}$.

The proof of this lemma is provided in . Applying this result to the cost function $R(i)=\frac{\sigma_{i+1}+\sigma_{1}\sqrt{i/\epsilon}}{\sigma_{i}}$, we get the following bounds :

Let, $\beta=\Sigma_{1}/\Sigma_{r}$ and $\xi=({M_{\rm max}}\sqrt{\alpha})/(\Sigma_{1}\sqrt{r})$. After some calculus, we establish that for

we have the desired inequality $R(r)<R(i)$ for all $i\neq r$. This proves the remark.

Acknowledgment

We thank Andrea Montanari for stimulating discussions and helpful comments on the subject of this paper. This work was partially supported by a Terman fellowship, the NSF CAREER award CCF-0743978 and the NSF grant DMS-0806211.

References