Matrix Completion With Noise

I Introduction

Imagine now that we only observe a few entries of a data matrix. Then is it possible to accurately—or even exactly—guess the entries that we have not seen? For example, suppose we observe a few movie ratings from a large data matrix in which rows are users and columns are movies (we can only observe a few ratings because each user is typically rating a few movies as opposed to the tens of thousands of movies which are available). Can we predict the rating a user would hypothetically assign to a movie he/she has not seen? In general, everybody would agree that recovering a data matrix from a subset of its entries is impossible. However, if the unknown matrix is known to have low rank or approximately low rank, then accurate and even exact recovery is possible by nuclear norm minimization . This revelation, which to some extent is inspired by the great body of work in compressed sensing, is the subject of this paper.

From now on, we will refer to the problem of inferring the many missing entries as the matrix completion problem. By extension, inferring a matrix from just a few linear functionals will be called the the low-rank matrix recovery problem. Now just as sparse signal recovery is arguably of paramount importance these days, we do believe that matrix completion and, in general, low-rank matrix recovery is just as important, and will become increasingly studied in years to come. For now, we give a few examples of applications in which these problems do come up.

Collaborative filtering. In a few words, collaborative filtering is the task of making automatic predictions about the interests of a user by collecting taste information from many users . Perhaps the most well-known implementation of collaborating filtering is the Netflix recommendation system alluded to earlier, which seeks to make rating predictions about unseen movies. This is a matrix completion problem in which the unknown full matrix has approximately low rank because only a few factors typically contribute to an individual’s tastes or preferences. In the new economy, companies are interested predicting musical preferences (Apple Inc.), literary preferences (Amazon, Barnes and Noble) and many other such things.

System identification. In control, one would like to fit a discrete-time linear time-invariant state-space model

Global positioning. Finding the global positioning of points in Euclidean space from a local or partial set of pairwise distances is a problem in geometry that emerges naturally in sensor networks . For example, because of power constraints, sensors may only be able to construct reliable distance estimates from their immediate neighbors. From these estimates, we can form a partially observed distance matrix, and the problem is to infer all the pairwise distances from just a few observed ones so that locations of the sensors can be reliably estimated. This reduces to a matrix completion problem where the unknown matrix is of rank two if the sensors are located in the plane, and three if they are located are in space.

Remote sensing. The MUSIC algorithm is frequently used to determine the direction of arrival of incident signals in a coherent radio-frequency environment. In a typical application, incoming signals are being recorded at various sensor locations, and this algorithm operates by extracting the directions of wave arrivals from the covariance matrix obtained by computing the correlations of the signals received at all sensor pairs. In remote sensing applications, one may not be able to estimate or transmit all correlations because of power constraints . In this case, we would like to infer a full covariance matrix from just a few observed partial correlations. This is a matrix completion problem in which the unknown signal covariance matrix has low rank since it is equal to the number of incident waves, which is usually much smaller than the number of sensors.

There are of course many other examples including the structure-from-motion problem in computer vision, multi-class learning in data analysis , and so on.

This paper investigates whether or not one can recover low rank matrices from fewer entries, and if so, how and how well. In Section II, we will study the noiseless problem in which the observed entries are precisely those of the unknown matrix. Section III examines the more common situation in which the few available entries are corrupted with noise. We complement our study with a few numerical experiments demonstrating the empirical performance of our methods in Section IV and conclude with a short discussion (Section V).

We use the usual asymptotic notation, for instance writing $O(M)$ to denote a quantity bounded in magnitude by $CM$ for some absolute constant $C>0$.

II Exact Matrix Completion

Thus, the question is whether it is possible to recover our matrix only from the information $\mathcal{P}_{\Omega}({M})$. We will assume that the entries are selected at random without replacement as to avoid trivial situations in which a row or a column is unsampled, since matrix completion is clearly impossible in such cases. (If we have no data about a specific user, how can we guess his/her preferences? If we have no distance estimates about a specific sensor, how can we guess its distances to all the sensors?)

This example shows that one cannot hope to complete the matrix if some of the singular vectors of the matrix are extremely sparse—above, one cannot recover $M$ without sampling all the entries in the first row, see for other related pathological examples. More generally, if a row (or column) has no relationship to the other rows (or columns) in the sense that it is approximately orthogonal, then one would basically need to see all the entries in that row to recover the matrix $M$. Such informal considerations led the authors of to introduce a geometric incoherence assumption, but for the moment, we will discuss an even simpler notion which forces the singular vectors of $M$ to be spread across all coordinates. To express this condition, recall the singular value decomposition (SVD) of a matrix of rank $r$,

If the singular vectors of $M$ are sufficiently spread, the hope is that there is a unique low-rank matrix which is consistent with the observed entries. If this is the case, one could, in principle, recover the unknown matrix by solving

A popular alternative is the convex relaxation

where $I_{n}$ is the $n\times n$ identity matrix. Hence, (II.4) is an SDP, which one can express by writing $\|X\|_{*}$ as the optimal value of the SDP dual to (II.5).

In , it is proven that nuclear-norm minimization succeeds nearly as soon as recovery is possible by any method whatsoever.

then $M$ is the unique solution to (II.4) with probability at least $1-n^{-3}$. In other words: with high probability, nuclear-norm minimization recovers all the entries of ${M}$ with no error.

As a side remark, one can obtain a probability of success at least $1-n^{-\beta}$ for $\beta$ by taking $C$ in (II.6) of the form $C^{\prime}\beta$ for some universal constant $C^{\prime}$.

An $n_{1}\times n_{2}$ matrix of rank $r$ depends upon $r(n_{1}+n_{2}-r)$ degrees of freedom This can be seen by counting the degrees of freedom in the singular value decomposition.. When $r$ is small, the number of degrees of freedom is much less than $n_{1}n_{2}$ and this is the reason why subsampling is possible. (In compressed sensing, the number of degrees of freedom corresponds to the sparsity of the signal; i.e. the number of nonzero entries.) What is remarkable here, is that exact recovery by nuclear norm minimization occurs as soon as the sample size exceeds the number of degrees of freedom by a couple of logarithmic factors. Further, observe that if $\Omega$ completely misses one of the rows (e.g. one has no rating about one user) or one of the columns (e.g. one has no rating about one movie), then one cannot hope to recover even a matrix of rank $1$ of the form $M=xy^{*}$. Thus one needs to sample every row (and also every column) of the matrix. When $\Omega$ is sampled at random, it is well established that one needs at least on the order $O(n\log n)$ for this to happen as this is the famous coupon collector’s problem. Hence, (II.6) misses the information theoretic limit by at most a logarithmic factor.

To obtain similar results for all values of the rank, introduces the strong incoherence property with parameter $\mu$ stated below.

Let $P_{U}$ (resp. $P_{V}$) be the orthogonal projection onto the singular vectors $u_{1},\ldots,u_{r}$ (resp. $v_{1},\ldots$, $v_{r}$). For all pairs $(a,a^{\prime})\in[n_{1}]\times[n_{1}]$ and $(b,b^{\prime})\in[n_{2}]\times[n_{2}]$,

These conditions do not assume anything about the singular values. As we will see, incoherent matrices with a small value of the strong incoherence parameter $\mu$ can be recovered from a minimal set of entries. Before we state this result, it is important to note that many model matrices obey the strong incoherence property with a small value of $\mu$.

Suppose the singular vectors obey (II.2) with $\mu_{B}=O(1)$ (which informally says that the singular vectors are not spiky), then with the exception of a very few peculiar matrices, $M$ obeys the strong incoherence property with $\mu=O(\sqrt{\log n})$.

Assume that the column matrices $[u_{1},\ldots,u_{r}]$ and $[v_{1},\ldots,v_{r}]$ are independent random orthogonal matrices, then with high probability, $M$ obeys the strong incoherence property with $\mu=O(\sqrt{\log n})$, at least when $r\geq\log n$ as to avoid small samples effects.

The sampling result below is general, nonasymptotic and optimal up to a few logarithmic factors.

With the same notations as in Theorem 1, there is a numerical constant $C$ such that if

$M$ is the unique solution to (II.4) with probability at least $1-n^{-3}$.

In other words, if a matrix is strongly incoherent and the cardinality of the sampled set is about the number of degrees of freedom times a few logarithmic factors, then nuclear-norm minimization is exact. This improves on an earlier result of Candès and Recht who proved—under slightly different assumptions—that on the order of $n^{6/5}r\log n$ samples were sufficient, at least for values of the rank obeying $r\leq n^{1/5}$.

We would like to point out a result of a broadly similar nature, but with a completely different recovery algorithm and with a somewhat different range of applicability, which was recently established by Keshavan, Oh, and Montanari . Their conditions are related to the incoherence property introduced in , and are also satisfied by a number of reasonable random matrix models. There is, however, another condition which states that the singular values of the unknown matrix cannot be too large or too small (the ratio between the top and lowest value must be bounded). This algorithm 1) trims each row and column with too few entries; i.e replaces the entries in those rows and columns by zero and 2) computes the SVD of the trimmed matrix and truncate it as to only keep the top $r$ singular values (note that the value of $r$ is needed here). The result is that under some suitable conditions discussed above, this recovers a good approximation to the matrix $M$ provided that the number of samples be on the order of $nr$. The recovery is not exact but only approximate although the authors have announced that one could add extra steps to the algorithm to provide an exact recovery if one has more samples (on the order of $nr\log n$). At the time of this writing, such results are not yet publicly available.

We cannot possibly rehash the proof of Theorem 2 from in this paper, or even explain the main technical steps, because of space limitations. We will, however, detail sufficient and almost necessary conditions for the low-rank matrix $M$ to be the unique solution to the SDP (II.4). This will be useful to establish stability results.

Recall the SVD (II.1) of $M$ and the “sign matrix” $E$ (II.7). It is is well-known that $Z\in\partial\|M\|_{*}$ if and only if $Z$ is of the form,

In English, $Z$ is a subgradient if it can be decomposed as the sign matrix plus another matrix with spectral norm bounded by one, whose column (resp. row) space is orthogonal to the span of $u_{1},\ldots,u_{r}$, (resp. of $v_{1},\ldots,v_{r}$). Another way to put this is by using notations introduced in . Let $T$ be the linear space spanned by elements of the form ${u}_{k}{x}^{*}$ and ${y}{v}_{k}^{*}$, $k\in[r]$, and let $T^{\perp}$ be the orthogonal complement to $T$. Note that $T^{\perp}$ is the set of matrices obeying $P_{U}W=0$ and $WP_{V}=0$. Then, $Z\in\partial\|M\|_{*}$ if and only if

We say that $\Lambda$ is a dual certificate if $\Lambda$ is supported on $\Omega$ ($\Lambda=\mathcal{P}_{\Omega}(\Lambda)$), $\mathcal{P}_{T}(\Lambda)=E$ and $\|\mathcal{P}_{T^{\perp}}(\Lambda)\|\leq 1$.

Hence, there is a clear analogy and one can think of $T$ defined above as playing the role of the support set in the sparse recovery problem.

With this in place, we shall make use of the following lemma from :

Suppose there exists a dual certificate $\Lambda$ and consider any $H$ obeying $\mathcal{P}_{\Omega}(H)=0$. Then

With $\Lambda=E+\mathcal{P}_{T^{\perp}}(\Lambda)$ and $Z=E+\mathcal{P}_{T^{\perp}}(Z)$, we have

since $\mathcal{P}_{\Omega}(H)=0$. Now we use the fact that the nuclear and spectral norms are dual to one another. In particular, there exists $\|Z\|\leq 1$ such that $\langle Z,\mathcal{P}_{T^{\perp}}(H)\rangle=\|\mathcal{P}_{T^{\perp}}(H)\|_{*}$ and $|\langle\Lambda,\mathcal{P}_{T^{\perp}}(H)\rangle|=|\langle\mathcal{P}_{T^{\perp}}(\Lambda),\mathcal{P}_{T^{\perp}}(H)\rangle|\leq\|\mathcal{P}_{T^{\perp}}(\Lambda)\|\|\mathcal{P}_{T^{\perp}}(H)\|_{*}$. Therefore,

which concludes the proof. ∎A consequence of this lemma are the sufficient conditions below.

Consider any feasible perturbation $M+H$ obeying $\mathcal{P}_{\Omega}(H)=0$. Then by assumption, Lemma 4 gives

unless $\mathcal{P}_{T^{\perp}}(H)=0$. Assume then that $\mathcal{P}_{T^{\perp}}({H})=0$; that is to say, ${H}\in T$. Then $\mathcal{P}_{\Omega}({H})=0$ implies that ${H}=0$ by the injectivity assumption. The conclusion is that $M$ is the unique minimizer since any nontrivial perturbation increases the nuclear norm. ∎

The methods for proving that matrix completion by nuclear minimization is exact, consist in constructing a dual certificate.

Under the assumptions of either Theorem 1 or Theorem 2, there exists a dual certificate obeying $\|\mathcal{P}_{T^{\perp}}(\Lambda)\|\leq 1/2$. In addition, if $p=m/(n_{1}n_{2})$ is the fraction of observed entries, the operator $\mathcal{P}_{T}\mathcal{P}_{\Omega}\mathcal{P}_{T}:T\rightarrow T$ is one-to-one and obeys

where $\mathcal{I}:T\rightarrow T$ is the identity operator.

III Stable Matrix Completion

In any real world application, one will only observe a few entries corrupted at least by a small amount of noise. In the Netflix problem, users’ ratings are uncertain. In the system identification problem, one cannot determine the locations $y(t)$ with infinite precision. In the global positioning problem, local distances are imperfect. And finally, in the remote sensing problem, the signal covariance matrix is always modeled as being corrupted by the covariance of noise signals. Hence, to be broadly applicable, we need to develop results which guarantee that reasonably accurate matrix completion is possible from noisy sampled entries. This section develops novel results showing that this is, indeed, the case.

where $\{Z_{ij}:(i,j)\in\Omega\}$ is a noise term which may be stochastic or deterministic (adversarial). Another way to express this model is as

where $Z$ is an $n\times n$ matrix with entries $Z_{ij}$ for $(i,j)\in\Omega$ (note that the values of $Z$ outside of $\Omega$ are irrelevant). All we assume is that $\|\mathcal{P}_{\Omega}(Z)\|_{F}\leq\delta$ for some $\delta>0$. For example, if $\{Z_{ij}\}$ is a white noise sequence with standard deviation $\sigma$, then $\delta^{2}\leq(m+\sqrt{8m})\sigma^{2}$ with high probability, say. To recover the unknown matrix, we propose solving the following optimization problem:

Among all matrices consistent with the data, find the one with minimum nuclear norm. This is also an SDP, and let $\hat{M}$ be the solution to this problem.

Our main result is that this reconstruction is accurate.

With the notations of Theorem 6, suppose there exists a dual certificate obeying $\|\mathcal{P}_{T^{\perp}}(\Lambda)\|\leq 1/2$ and that $\mathcal{P}_{T}\mathcal{P}_{\Omega}\mathcal{P}_{T}\succeq\frac{p}{2}\mathcal{I}$ (both these conditions are true with very large probability under the assumptions of the noiseless recovery Theorems 1 and 2). Then $\hat{M}$ obeys

For small values of $p$ (recall this is the fraction of observed entries), the error is of course at most just about $4\sqrt{\frac{2\min(n_{1},n_{2})}{p}}\,\delta$. As we will see from the proof, there is nothing special about 1/2 in the condition $\|\mathcal{P}_{T^{\perp}}(\Lambda)\|\leq 1/2$. All we need is that there is a dual certificate obeying $\|\mathcal{P}_{T^{\perp}}(\Lambda)\|\leq a$ for some $a<1$ (the value of $a$ only influences the numerical constant in (III.3)). Further, when $Z$ is random, (III.3) holds on the event $\|\mathcal{P}_{\Omega}(Z)\|_{F}\leq\delta$.

Roughly speaking, our theorem states the following: when perfect noiseless recovery occurs, then matrix completion is stable vis a vis perturbations. To be sure, the error is proportional to the noise level $\delta$; when the noise level is small, the error is small. Moreover, improving conditions under which noiseless recovery occurs, has automatic consequences for the more realistic recovery from noisy samples.

A significant novelty here is that there is just no equivalent of this result in the compressed sensing or statistical literature for our matrix completion problem does not obey the restricted isometry property (RIP) . For matrices, the RIP would assume that the sampling operator obeys

for all matrices $X$ with sufficiently small rank and $\delta<1$ sufficiently small . However, the RIP does not hold here. To see why, let the sampled set $\Omega$ be arbitrarily chosen and fix $(i,j)\notin\Omega$. Then the rank-1 matrix $e_{i}e_{j}^{*}$ whose $(i,j)$th entry is 1, and vanishes everywhere else, obeys $\mathcal{P}_{\Omega}(e_{i}e_{j}^{*})=0$. Clearly, this violates (III.4).

It is nevertheless instructive to compare (III.3) with the bound one would achieve if the RIP (III.4) were true. In this case, would give

for some numerical constant $C_{0}$. That is, an estimate which would be better by a factor proportional to $1/\sqrt{\min(n_{1},n_{2})}$. It would be interesting to know whether or not estimates, which are as good as what is achievable under the RIP, hold for the RIPless matrix completion problem. We will return to such comparisons later (Section III-B).

We close this section by emphasizing that our methods are also applicable to sparse signal recovery problems in which the RIP does not hold.

We use the notation of the previous section, and begin the proof by observing two elementary properties. The first is that since $M$ is feasible for (III.2), we have the cone constraint

The second is that the triangle inequality implies the tube constraint

since $M$ is feasible. We will see that under our hypotheses, (III.5) and (III.6) imply that $\hat{M}$ is close to $M$. Set $\hat{M}=M+H$ and put $H_{\Omega}:=\mathcal{P}_{\Omega}(H)$, $H_{\Omega^{c}}:=\mathcal{P}_{\Omega^{c}}(H)$ for short. We need to bound $\|H\|^{2}_{F}=\|H_{\Omega}\|_{F}^{2}+\|H_{\Omega^{c}}\|_{F}^{2}$, and since (III.6) gives $\|H_{\Omega}\|_{F}\leq 2\delta$, it suffices to bound $\|H_{\Omega^{c}}\|_{F}$. Note that by the Pythagorean identity, we have

and it is thus sufficient to bound each term in the right hand-side.

We start with the second term. Let $\Lambda$ be a dual certificate obeying $\|\mathcal{P}_{T^{\perp}}(\Lambda)\|\leq 1/2$, we have

The second inequality follows from Lemma 4. Therefore, with $\|\mathcal{P}_{T^{\perp}}(\Lambda)\|\leq 1/2$, the cone constraint gives

Since the nuclear norm dominates the Frobenius norm, $\|\mathcal{P}_{T^{\perp}}(H_{\Omega^{c}})\|_{F}\leq\|\mathcal{P}_{T^{\perp}}(H_{\Omega^{c}})\|_{*}$, we have

where the second inequality follows from the Cauchy-Schwarz inequality, and the last from (III.6).

To develop a bound on $\|\mathcal{P}_{T}(H_{\Omega})\|_{F}$, observe that the assumption $\mathcal{P}_{T}\mathcal{P}_{\Omega}\mathcal{P}_{T}\succeq\frac{p}{2}\mathcal{I}$ together with $\mathcal{P}_{T}^{2}=\mathcal{P}_{T}$, $\mathcal{P}_{\Omega}^{2}=\mathcal{P}_{\Omega}$ give

But since $\mathcal{P}_{\Omega}(H_{\Omega^{c}})=0=\mathcal{P}_{\Omega}\mathcal{P}_{T}(H_{\Omega^{c}})+\mathcal{P}_{\Omega}\mathcal{P}_{T^{\perp}}(H_{\Omega^{c}})$, we have

As a consequence of this and (III.7), we have

The theorem then follows from this inequality together with (III.8).

III-B Comparison with an oracle

We would like to return to discussing the best possible accuracy one could ever hope for. For simplicity, assume that $n_{1}=n_{2}=n$, and suppose that we have an oracle informing us about $T$. In many ways, going back to the discussion from Section II-A, this is analogous to giving away the support of the signal in compressed sensing . With this precious information, we would know that $M$ lives in a linear space of dimension $2nr-r^{2}$ and would probably solve the problem by the method of least squares:

That is, we would find the matrix in $T$, which best fits the data in a least-squares sense. Let $\mathcal{A}:T\rightarrow\Omega$ (we abuse notations and let $\Omega$ be the range of $\mathcal{P}_{\Omega}$) defined by $\mathcal{A}:=\mathcal{P}_{\Omega}\mathcal{P}_{T}$. Then assuming that the operator $\mathcal{A}^{*}\mathcal{A}=\mathcal{P}_{T}\mathcal{P}_{\Omega}\mathcal{P}_{T}$ mapping $T$ onto $T$ is invertible (which is the case under the hypotheses of Theorem 7), the least-squares solution is given by

Let $Z^{\prime}$ be the minimal (normalized) eigenvector of $\mathcal{A}^{*}\mathcal{A}$ with minimum eigenvalue $\lambda_{\text{min}}$, and set $Z=\delta\lambda_{\text{min}}^{-1/2}\mathcal{A}(Z^{\prime})$ (note that by definition $\mathcal{P}_{\Omega}(Z)=Z$ since $Z$ is in the range of $\mathcal{A}$). By construction, $\|Z\|_{F}=\delta$, and

since by assumption, all the eigenvalues of $\mathcal{A}^{*}\mathcal{A}=\mathcal{P}_{T}\mathcal{P}_{\Omega}\mathcal{P}_{T}$ lie in the interval $[p/2,3p/2]$. The matrix $Z$ defined above also maximizes $\|(\mathcal{A}^{*}\mathcal{A})^{-1}\mathcal{A}^{*}(Z)\|_{F}$ among all matrices bounded by $\delta$ and so the oracle achieves

with adversarial noise. Consequently, our analysis looses a $\sqrt{n}$ factor vis a vis an optimal bound that is achievable via the help of an oracle.

The diligent reader may argue that the least-squares solution above may not be of rank $r$ (it is at most of rank $2r$) and may thus argue that this is not the strongest possible oracle. However, as explained below, if the oracle gave $T$ and $r$, then the best fit in $T$ of rank $r$ would not do much better than (III.12). In fact, there is an elegant way to understand the significance of this oracle which we now present. Consider a stronger oracle which reveals the row space of the unknown matrix $M$ (and thus the rank of the matrix). Then we would know that the unknown matrix is of the form

Using our previous notations, the oracle gives away $T_{0}\subset T$ where $T_{0}$ is the span of elements of the form $yv_{k}^{*}$, $k\in[r]$, and is more precise. If $\mathcal{A}_{0}:T_{0}\rightarrow\Omega$ is defined by $\mathcal{A}_{0}:=\mathcal{P}_{\Omega}\mathcal{P}_{T_{0}}$, then the least-squares solution is now

Because all the eigenvalues of $\mathcal{A}_{0}^{*}\mathcal{A}_{0}$ belong to $[\lambda_{\text{min}}(\mathcal{A}^{*}\mathcal{A}),\lambda_{\text{max}}(\mathcal{A}^{*}\mathcal{A})]$, the previous analysis applies and this stronger oracle would also achieve an error of size about $p^{-1/2}\delta$. In conclusion, when all we know is $\|\mathcal{P}_{\Omega}(Z)\|_{F}\leq\delta$, one cannot hope for a root-mean squared error better than $p^{-1/2}\delta$.

since all the $2nr-r^{2}$ eigenvalues of $(\mathcal{A}^{*}\mathcal{A})^{-1}$ are just about equal to $p^{-1}$. When $nr\ll m$, this is better than (III.12).

IV Numerical Experiments

For a peek at the results, consider Table I.

In order to stably recover $M$ from a fraction of noisy entries, the following regularized nuclear norm minimization problem was solved using the FPC algorithm from ,

It is a standard duality result that (IV.1) is equivalent to (III.2), for some value of $\mu$, and thus one could use (IV.1) to solve (III.2) by searching for the value of $\mu(\delta)$ giving $\|\mathcal{P}_{\Omega}(\hat{M}-Y)\|_{F}=\delta$ (assuming $\|\mathcal{P}_{\Omega}(Y)\|_{F}>\delta$). We use (IV.1) because it works well in practice, and because the FPC algorithm solves (IV.1) nicely and accurately. We also remark that a variation on our stability proof could also give a stable error bound when using the SDP (IV.1).

In order to interpret our numerical results, they are compared to those achieved by the oracle, see Section III-B. To this end, Figure 2 plots three curves for varying values of $n,p,$ and $r$: 1) the RMS error introduced above, 2) the RMS error achievable when the oracle reveals $T$, and the problem is solved using least squares, 3) the estimated oracle root expected MS error derived in Section III-B, i.e. $\sqrt{\text{df}/[n^{2}p]}=\sqrt{\text{df}/m}$, where $\text{df}=r(2n-r)$. In our experiments, as $n$ and $m/\text{df}$ increased, with $r=2$, the RMS error of the nuclear norm problem appeared to be fit very well by $1.68\,\sqrt{\text{df}/m}$. Thus, to compare the oracle error to the actual recovered error, we plotted the oracle errors times 1.68. We also note that in our experiments, the RMS error was never greater than $2.25\,\sqrt{\text{df}/m}$.

No one can predict the weather. We conclude the numerical section with a real world example. We retrieved from the website a $366\times 1472$ matrix whose entries are daily average temperatures at 1472 different weather stations throughout the world in 2008. Checking its SVD reveals that this is an approximately low rank matrix as expected. In fact, letting $M$ be the temperature matrix, and calling $M_{2}$ the matrix created by truncating the SVD after the top two singular values gives $\|M_{2}\|_{F}/\|M\|_{F}=.9927$.

To test the performance of our matrix completion algorithm, we subsampled 30% of $M$ and then recovered an estimate, $\hat{M}$, using (IV.1). Note that this is a much different problem than those proposed earlier in this section. Here, we attempt to recover a matrix that is not exactly low rank, but only approximately. The solution gives a relative error of $\|\hat{M}-M\|_{F}/\|M\|_{F}=.166$. For comparison The number 2 is somewhat arbitrary here, although we picked it because there is a large drop-off in the size of the singular values after the second. If, for example, $M_{10}$ is the best rank-10 approximation, then $\|M_{10}-M\|_{F}/\|M\|_{F}=.081$., exact knowledge of the best rank-2 approximation achieves $\|M_{2}-M\|_{F}/\|M\|_{F}=.121$. Here $\mu$ has been selected to give a good cross-validated error and is about 535.

V Discussion

This paper reviewed and developed some new results about matrix completion. By and large, low-rank matrix recovery is a field in complete infancy abounding with interesting and open questions, and if the recent avalanche of results in compressed sensing is any indication, it is likely that this field will experience tremendous growth in the next few years.

At a computational level, one would like to have available a suite of efficient algorithms for minimizing the nuclear norm under convex constraints and, in general, for finding low-rank matrices obeying convex constraints. Algorithms with impressive performance in some situations have already been proposed but the computational challenges of solving problems with millions if not billions of unknowns obviously still require much research.

E. C. is supported by ONR grants N00014-09-1-0469 and N00014-08-1-0749 and by the Waterman Award from NSF. E. C. would like to thank Terence Tao and Stephen Becker for some very helpful discussions.

References