Matrix completion

A number of assumptions on the sampling of the observed entries and the number of sampled entries are frequently made to simplify the analysis and to ensure the problem is not underdetermined.

To make the analysis tractable, it is often assumed that the set ${\displaystyle E}$ of observed entries and fixed cardinality is sampled uniformly at random from the collection of all subsets of entries of cardinality ${\displaystyle |E|}$. To further simplify the analysis, it is instead assumed that ${\displaystyle E}$ is constructed by Bernoulli sampling, i.e. that each entry is observed with probability ${\displaystyle p}$. If ${\displaystyle p}$ is set to ${\displaystyle {\frac {N}{mn}}}$ where ${\displaystyle N}$ is the desired expected cardinality of ${\displaystyle E}$, and ${\displaystyle m,\;n}$ are the dimensions of the matrix (let ${\displaystyle m<n}$ without loss of generality), ${\displaystyle |E|}$ is within ${\displaystyle O(n\log n)}$ of ${\displaystyle N}$ with high probability, thus Bernoulli sampling is a good approximation for uniform sampling.[1] Another simplification is to assume that entries are sampled independently and with replacement.[2]

Suppose the ${\displaystyle m}$ by ${\displaystyle n}$ matrix ${\displaystyle M}$ (with ${\displaystyle m<n}$) we are trying to recover has rank ${\displaystyle r}$. There is an information theoretic lower bound on how many entries must be observed before ${\displaystyle M}$ can be uniquely reconstructed. The set of ${\displaystyle m}$ by ${\displaystyle n}$ matrices with rank less than or equal to ${\displaystyle r}$ is an algebraic variety in ${\displaystyle {\mathbb {C} }^{m\times n}}$with dimension ${\displaystyle (n+m)r-r^{2}}$. Using this result, one can show that at least ${\displaystyle 4nr-4r^{2}}$ entries must be observed for matrix completion in ${\displaystyle {\mathbb {C} }^{n\times n}}$ to have a unique solution when ${\displaystyle r\leq n/2}$ .[3]

Secondly, there must be at least one observed entry per row and column of ${\displaystyle M}$. The Singular Value Decomposition of ${\displaystyle M}$ is given by ${\displaystyle U\Sigma V^{\dagger }}$. If row ${\displaystyle i}$ is unobserved, it is easy to see the ${\displaystyle i^{\text{th}}}$ right singular vector of ${\displaystyle M}$, ${\displaystyle v_{i}}$, can be changed to some arbitrary value and still yield a matrix matching ${\displaystyle M}$ over the set of observed entries. Similarly, if column ${\displaystyle j}$ is unobserved, the ${\displaystyle j^{\text{th}}}$ left singular vector of ${\displaystyle M}$, ${\displaystyle u_{i}}$ can be arbitrary. If we assume Bernoulli sampling of the set of observed entries, the Coupon collector effect implies that entries on the order of ${\displaystyle O(n\log n)}$ must be observed to ensure that there is an observation from each row and column with high probability.[4]

Combining the necessary conditions and assuming that ${\displaystyle r\ll m,n}$ (a valid assumption for many practical applications), the lower bound on the number of observed entries required to prevent the problem of matrix completion from being underdetermined is on the order of ${\displaystyle nr\log n}$.

The concept of incoherence arose in compressed sensing. It is introduced in the context of matrix completion to ensure the singular vectors of ${\displaystyle M}$ are not too "sparse" in the sense that all coordinates of each singular vector are of comparable magnitude instead of just a few coordinates having significantly larger magnitudes.[5][6] The standard basis vectors are then undesirable as singular vectors, and the vector ${\displaystyle {\frac {1}{\sqrt {n}}}{\begin{bmatrix}1\\1\\\vdots \\1\end{bmatrix}}}$ in ${\displaystyle \mathbb {R} ^{n}}$ is desirable. As an example of what could go wrong if the singular vectors are sufficiently "sparse", consider the ${\displaystyle m}$ by ${\displaystyle n}$ matrix ${\displaystyle {\begin{bmatrix}1&0&\cdots &0\\\vdots &&\vdots \\0&0&0&0\end{bmatrix}}}$ with singular value decomposition ${\displaystyle I_{m}{\begin{bmatrix}1&0&\cdots &0\\\vdots &&\vdots \\0&0&0&0\end{bmatrix}}I_{n}}$. Almost all the entries of ${\displaystyle M}$ must be sampled before it can be reconstructed.

Candès and Recht[1] define the coherence of a matrix ${\displaystyle U}$ with column space an ${\displaystyle r-}$dimensional subspace of ${\displaystyle \mathbb {R} ^{n}}$ as ${\displaystyle \mu (U)={\frac {n}{r}}\max _{i<n}\|P_{U}e_{i}\|^{2}}$, where ${\displaystyle P_{U}}$ is the orthogonal projection onto ${\displaystyle U}$. Incoherence then asserts that given the singular value decomposition ${\displaystyle U\Sigma V^{\dagger }}$ of the ${\displaystyle m}$ by ${\displaystyle n}$ matrix ${\displaystyle M}$,

1. ${\displaystyle \mu (U),\;\mu (V)\leq \mu _{0}}$
2. The entries of ${\displaystyle \sum _{k}u_{k}v_{k}^{\dagger }}$ have magnitudes upper bounded by ${\displaystyle \mu _{1}{\sqrt {\frac {r}{mn}}}}$

for some ${\displaystyle \mu _{0},\;\mu _{1}}$.

The rank minimization problem is NP-hard. One approach, proposed by Candès and Recht, is to form a convex relaxation of the problem and minimize the nuclear norm ${\displaystyle \|M\|_{*}}$ (which gives the sum of the singular values of ${\displaystyle M}$) instead of ${\displaystyle {\text{rank}}(M)}$ (which counts the number of non zero singular values of ${\displaystyle M}$).[1] This is analogous to minimizing the L1-norm rather than the L0-norm for vectors. The convex relaxation can be solved using semidefinite programming (SDP) by noticing that the optimization problem is equivalent to

${\displaystyle {\begin{aligned}&{\underset {W_{1},W_{2}}{\text{min}}}&&{\text{trace}}(W_{1})+{\text{trace}}(W_{2})\\&{\text{subject to}}&&X_{ij}=M_{ij}\;\;\forall i,j\in E\\&&&{\begin{bmatrix}W_{1}&X\\X^{T}&W_{2}\end{bmatrix}}\succeq 0\end{aligned}}}$

The complexity of using SDP to solve the convex relaxation is ${\displaystyle O({\text{max}}(m,n)^{4})}$. State of the art solvers like SDP3 can only handle matrices of size up to 100 by 100 [11] An alternative first order method that approximately solves the convex relaxation is the Singular Value Thresholding Algorithm introduced by Cai, Candès and Shen.[11]

Candès and Recht show, using the study of random variables on Banach spaces, that if the number of observed entries is on the order of ${\displaystyle \max {\{\mu _{1}^{2},{\sqrt {\mu _{0}}}\mu _{1},\mu _{0}n^{0.25}\}}nr\log n}$ (assume without loss of generality ${\displaystyle m<n}$), the rank minimization problem has a unique solution which also happens to be the solution of its convex relaxation with probability ${\displaystyle 1-{\frac {c}{n^{3}}}}$ for some constant ${\displaystyle c}$. If the rank of ${\displaystyle M}$ is small (${\displaystyle r\leq {\frac {n^{0.2}}{\mu _{0}}}}$), the size of the set of observations reduces to the order of ${\displaystyle \mu _{0}n^{1.2}r\log n}$. These results are near optimal, since the minimum number of entries that must be observed for the matrix completion problem to not be underdetermined is on the order of ${\displaystyle nr\log n}$.

This result has been improved by Candès and Tao.[4] They achieve bounds that differ from the optimal bounds only by polylogarithmic factors by strengthening the assumptions. Instead of the incoherence property, they assume the strong incoherence property with parameter ${\displaystyle \mu _{3}}$. This property states that:

1. ${\displaystyle |\langle e_{a},P_{U}e_{a'}\rangle -{\frac {r}{m}}1_{a=a'}|\leq \mu _{3}{\frac {\sqrt {r}}{m}}}$ for ${\displaystyle a,a'\leq m}$ and ${\displaystyle |\langle e_{b},P_{U}e_{b'}\rangle -{\frac {r}{n}}1_{b=b'}|\leq \mu _{3}{\frac {\sqrt {r}}{n}}}$ for ${\displaystyle b,b'\leq n}$
2. The entries of ${\displaystyle \sum _{i}u_{i}v_{i}^{\dagger }}$ are bounded in magnitude by ${\displaystyle \mu _{3}{\sqrt {\frac {r}{mn}}}}$

Intuitively, strong incoherence of a matrix ${\displaystyle U}$ asserts that the orthogonal projections of standard basis vectors to ${\displaystyle U}$ has magnitudes that have high likelihood if the singular vectors were distributed randomly.[5]

Candès and Tao find that when ${\displaystyle r}$ is ${\displaystyle O(1)}$ and the number of observed entries is on the order of ${\displaystyle \mu _{3}^{4}n(\log n)^{2}}$, the rank minimization problem has a unique solution which also happens to be the solution of its convex relaxation with probability ${\displaystyle 1-{\frac {c}{n^{3}}}}$ for some constant ${\displaystyle c}$. For arbitrary ${\displaystyle r}$, the number of observed entries sufficient for this assertion hold true is on the order of ${\displaystyle \mu _{3}^{2}nr(\log n)^{6}}$

Keshavan, Montanari and Oh[9] consider a variant of matrix completion where the rank of the ${\displaystyle m}$ by ${\displaystyle n}$ matrix ${\displaystyle M}$, which is to be recovered, is known to be ${\displaystyle r}$. They assume Bernoulli sampling of entries, constant aspect ratio ${\displaystyle {\frac {m}{n}}}$, bounded magnitude of entries of ${\displaystyle M}$ (let the upper bound be ${\displaystyle M_{\text{max}}}$), and constant condition number ${\displaystyle {\frac {\sigma _{1}}{\sigma _{r}}}}$ (where ${\displaystyle \sigma _{1}}$ and ${\displaystyle \sigma _{r}}$ are the largest and smallest singular values of ${\displaystyle M}$ respectively). Further, they assume the two incoherence conditions are satisfied with ${\displaystyle \mu _{0}}$ and ${\displaystyle \mu _{1}{\frac {\sigma _{1}}{\sigma _{r}}}}$ where ${\displaystyle \mu _{0}}$ and ${\displaystyle \mu _{1}}$ are constants. Let ${\displaystyle M^{E}}$ be a matrix that matches ${\displaystyle M}$ on the set ${\displaystyle E}$ of observed entries and is 0 elsewhere. They then propose the following algorithm:

1. Trim ${\displaystyle M^{E}}$ by removing all observations from columns with degree larger than ${\displaystyle {\frac {2|E|}{n}}}$ by setting the entries in the columns to 0. Similarly remove all observations from rows with degree larger than ${\displaystyle {\frac {2|E|}{n}}}$.
2. Project ${\displaystyle M^{E}}$ onto its first ${\displaystyle r}$ principal components. Call the resulting matrix ${\displaystyle {\text{Tr}}(M^{E})}$.
3. Solve ${\displaystyle \min _{X,Y}\min _{S\in \mathbb {R} ^{r\times r}}{\frac {1}{2}}\sum _{i,j\in E}(M_{ij}-(XSY^{\dagger })_{ij})^{2}+\rho G(X,Y)}$ where ${\displaystyle G(X,Y)}$ is some regularization function by gradient descent with line search. Initialize ${\displaystyle X,\;Y}$ at ${\displaystyle X_{0},\;Y_{0}}$ where ${\displaystyle {\text{Tr}}(M_{E})=X_{0}S_{0}Y_{0}^{\dagger }}$. Set ${\displaystyle G(X,Y)}$ as some function forcing ${\displaystyle X,\;Y}$ to remain incoherent throughout gradient descent if ${\displaystyle X_{0}}$ and ${\displaystyle Y_{0}}$ are incoherent.
4. Return the matrix ${\displaystyle XSY^{\dagger }}$.

Steps 1 and 2 of the algorithm yield a matrix ${\displaystyle {\text{Tr}}(M^{E})}$ very close to the true matrix ${\displaystyle M}$ (as measured by the root mean square error (RMSE) with high probability. In particular, with probability ${\displaystyle 1-{\frac {1}{n^{3}}}}$, ${\displaystyle {\frac {1}{mnM_{\text{max}}^{2}}}\|M-{\text{Tr}}(M^{E})\|_{F}^{2}\leq C{\frac {r}{m|E|}}{\sqrt {\frac {m}{n}}}}$ for some constant ${\displaystyle C}$. ${\displaystyle \|\cdot \|_{F}}$ denotes the Frobenius norm. Note that the full suite of assumptions is not needed for this result to hold. The incoherence condition, for example, only comes into play in exact reconstruction. Finally, although trimming may seem counter intuitive as it involves throwing out information, it ensures projecting ${\displaystyle M^{E}}$ onto its first ${\displaystyle r}$ principal components gives more information about the underlying matrix ${\displaystyle M}$ than about the observed entries.

In Step 3, the space of candidate matrices ${\displaystyle X,\;Y}$ can be reduced by noticing that the inner minimization problem has the same solution for ${\displaystyle (X,Y)}$ as for ${\displaystyle (XQ,YR)}$ where ${\displaystyle Q}$ and ${\displaystyle R}$ are orthonormal ${\displaystyle r}$ by ${\displaystyle r}$ matrices. Then gradient descent can be performed over the cross product of two Grassman manifolds. If ${\displaystyle r\ll m,\;n}$ and the observed entry set is in the order of ${\displaystyle nr\log n}$, the matrix returned by Step 3 is exactly ${\displaystyle M}$. Then the algorithm is order optimal, since we know that for the matrix completion problem to not be underdetermined the number of entries must be in the order of ${\displaystyle nr\log n}$.

Alternating minimization represents a widely applicable and empirically successful approach for finding low-rank matrices that best fit the given data. For example, for the problem of low-rank matrix completion, this method is believed to be one of the most accurate and efficient, and formed a major component of the winning entry in the Netflix problem. In the alternating minimization approach, the low-rank target matrix is written in a bilinear form:

${\displaystyle X=UV^{T}}$;

the algorithm then alternates between finding the best ${\displaystyle U}$ and the best ${\displaystyle V}$. While the overall problem is non-convex, each sub-problem is typically convex and can be solved efficiently. Jain, Netrapalli and Sanghavi [10] have given one of the first guarantees for performance of alternating minimization for both matrix completion and matrix sensing.

The alternating minimization algorithm can be viewed as an approximate way to solve the following non-convex problem:

${\displaystyle {\begin{aligned}&{\underset {U,V\in \mathbb {R} ^{n\times k}}{\text{min}}}&\|P_{\Omega }(UV^{T})-P_{\Omega }(M)\|_{F}^{2}\\\end{aligned}}}$

The AltMinComplete Algorithm proposed by Jain, Netrapalli and Sanghavi is listed here:[10]

1. Input: observed set ${\displaystyle \Omega }$, values ${\displaystyle P_{\Omega }(M)}$
2. Partition ${\displaystyle \Omega }$ into ${\displaystyle 2T+1}$ subsets ${\displaystyle \Omega _{0},\cdots ,\Omega _{2T}}$ with each element of ${\displaystyle \Omega }$ belonging to one of the ${\displaystyle \Omega _{t}}$ with equal probability (sampling with replacement)
3. ${\displaystyle {\hat {U}}^{0}=SVD({\frac {1}{p}}P_{\Omega _{0}}(M),k)}$ i.e., top-${\displaystyle k}$ left singular vectors of ${\displaystyle {\frac {1}{p}}P_{\Omega _{0}}(M)}$
4. Clipping: Set all elements of ${\displaystyle {\hat {U}}^{0}}$ that have magnitude greater than ${\displaystyle {\frac {2\mu {\sqrt {k}}}{\sqrt {n}}}}$ to zero and orthonormalize the columns of ${\displaystyle {\hat {U}}^{0}}$
5. for ${\displaystyle t=0,\cdots ,T-1}$ do
6. ${\displaystyle \quad {\hat {V}}^{t+1}\leftarrow {\text{argmin}}_{V\in \mathbb {R} ^{n\times k}}\|P_{\Omega _{t+1}}({\hat {U}}V^{T}-M)\|_{F}^{2}}$
7. ${\displaystyle \quad {\hat {U}}^{t+1}\leftarrow {\text{argmin}}_{U\in \mathbb {R} ^{m\times k}}\|P_{\Omega _{T+t+1}}(U({\hat {V}}^{t+1})^{T}-M)\|_{F}^{2}}$
8. end for
9. Return ${\displaystyle X={\hat {U}}^{T}({\hat {V}}^{T})^{T}}$

They showed that by observing ${\displaystyle |\Omega |=O(({\frac {\sigma _{1}^{*}}{\sigma _{k}^{*}}})^{6}k^{7}\log n\log(k\|M\|_{F}/\epsilon ))}$ random entries of an incoherent matrix ${\displaystyle M}$, AltMinComplete algorithm can recover ${\displaystyle M}$ in ${\displaystyle O(\log(1/\epsilon ))}$ steps. In terms of sample complexity (${\displaystyle |\Omega |}$), theoretically, Alternating Minimization may require a bigger ${\displaystyle \Omega }$ than Convex Relaxation. However empirically it seems not the case which implies that the sample complexity bounds can be further tightened. In terms of time complexity, they showed that AltMinComplete needs time

${\displaystyle O(|\Omega |k^{2}\log(1/\epsilon ))}$.

It is worth noting that, although convex relaxation based methods have rigorous analysis, alternating minimization based algorithms are more successful in practice.

Collaborative filtering is the task of making automatic predictions about the interests of a user by collecting taste information from many users. Companies like Apple, Amazon, Barnes and Noble, and Netflix are trying to predict their user preferences from partial knowledge. In these kind of matrix completion problem, the unknown full matrix is often considered low rank because only a few factors typically contribute to an individual's tastes or preference.

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

${\displaystyle {\begin{aligned}x(t+1)&=Ax(t)+Bu(t)\\y(t)&=Cx(t)+Du(t)\end{aligned}}}$

to a sequence of inputs ${\displaystyle u(t)\in \mathbb {R} ^{m}}$ and outputs ${\displaystyle y(t)\in \mathbb {R} ^{p},t=0,\ldots ,N}$. The vector ${\displaystyle x(t)\in \mathbb {R} ^{n}}$ is the state of the system at time ${\displaystyle t}$ and ${\displaystyle n}$ is the order of the system model. From the input/output pair, one would like to recover the matrices ${\displaystyle A,B,C,D}$ and the initial state ${\displaystyle x(0)}$. This problem can also be viewed as a low-rank matrix completion problem.

The localization (or global positioning) problem emerges naturally in IoT sensor networks. The problem is to recover the sensor map in Euclidean space from a local or partial set of pairwise distances. Thus it is a matrix completion problem with rank two if the sensors are located in a 2-D plane and three if they are in a 3-D space.[12]