Fast Exact Matrix Completion with Finite Samples

Introduction

where $P_{\Omega}\left(\bm{A}\right)$ is defined as:

LRMC is by now a well studied problem with applications in several machine learning tasks such as collaborative filtering [BK07], link analysis [GL11], distance embedding [CR09] etc. Motivated by widespread applications, several practical algorithms have been proposed to solve the problem (heuristically) [RR13, HCD12].

On the theoretical front, the non-convex rank constraint implies NP-hardness in general [HMRW14]. However, under certain (by now) standard assumptions, a few algorithms have been shown to solve the problem efficiently. These approaches can be categorized into the following two broad groups:

a) The first approach relaxes the rank constraint in (1) to a trace norm constraint (sum of singular values of $\bm{X}$) and then solves the resulting convex optimization problem [CR09]. [CT09, Rec09] showed that this approach has a near optimal sample complexity (i.e. the number of observed entries of $\mathbf{M}$) of $\left\lvert\Omega\right\rvert=O\left({rn\log^{2}n}\right)$, where we abbreviate $n=n_{1}+n_{2}$. However, current iterative algorithms used to solve the trace-norm constrained optimization problem require $O\left({n^{2}}\right)$ memory and $O\left({n^{3}}\right)$ time per iteration, which is prohibitive for large-scale applications.

b) The second approach is based on an empirically popular iterative technique called Alternating Minimization (AltMin) that factorizes $\bm{X}=\bm{U}{\bm{V}}^{\top}$ where $\bm{U},\bm{V}$ have $r$ columns, and the algorithm alternately optimizes over $\bm{U}$ and $\bm{V}$ holding the other fixed. Recently, [Kes12, JNS13, Har14, HW14] showed convergence of variants of this algorithm. The best known sample complexity results for AltMin are the incomparable bounds $|\Omega|=O\left({r\kappa^{8}n{\log\frac{n}{\epsilon}}}\right)$ and $|\Omega|=O\left({\;\textrm{poly}\left(r\right)\left(\log\kappa\right){n\log\frac{n}{\epsilon}}}\right)$ due to [Kes12] and [HW14] respectively. Here, $\kappa=\sigma_{1}(\mathbf{M})/\sigma_{r}(\mathbf{M})$ is the condition number of $\mathbf{M}$ and $\epsilon$ is the desired accuracy. The computational cost of these methods is $O\left({|\Omega|r+nr^{3}}\right)$ per iteration, making these methods very fast as long as the condition number $\kappa$ is not too large.

Of the above two approaches AltMin is known to be the most practical and runs in near linear time. However, its sample complexity as well as computational complexity depend on the condition number of $\mathbf{M}$ which can be arbitrarily large. Moreover, for “exact” recovery of $\mathbf{M}$, i.e., with error $\epsilon=0$, the method requires infinitely many samples (or rather to observe the entire matrix). The dependence of sample complexity on the desired accuracy arises due to the use of independent samples in each iteration, which in turn is necessitated by the fact that using the same samples in each iteration leads to complex dependencies among iterates which are hard to analyze. Nevertheless, practitioners have been using AltMin with same samples in each iteration successfully in a wide range of applications.

Our results: In this paper, we address this issue by proposing a new algorithm called Stagewise-SVP (St-SVP) and showing that it solves the matrix completion problem exactly with a sample complexity $\left\lvert\Omega\right\rvert=O\left({nr^{5}\log^{3}n}\right)$, which is independent of both the condition number, and desired accuracy and time complexity per iteration $O\left({\left\lvert\Omega\right\rvert r^{2}}\right)$, which is near linear in $n$.

The basic block of our algorithm is a simple projected gradient descent step, first proposed by [JMD10] in the context of this problem. More precisely, given the $t^{\textrm{th}}$ iterate $\bm{X}_{t}$, [JMD10] proposed the following update rule, which they call singular value projection (SVP).

where $P_{r}$ is the projection onto the set of rank-$r$ matrices and can be efficiently computed using singular value decomposition (SVD). Note that the SVP step is just a projected gradient descent step where the projection is onto the (non-convex) set of low rank matrices. [JMD10] showed that despite involving projections onto a non-convex set, SVP solves the related problem of low-rank matrix sensing, where instead of observing elements of the unknown matrix, we observe dense linear measurements of this matrix. However, their result does not extend to the matrix completion problem and the correctness of SVP for matrix completion was left as an open question.

Our preliminary result resolves this question by showing the correctness of SVP for the matrix completion problem, albeit with a sample complexity that depends on the condition number and desired accuracy. We then develop a stage-wise variant of this algorithm, where in the $k^{\textrm{th}}$ stage, we try to recover $P_{k}\left(\mathbf{M}\right)$, there by getting rid of the dependence on the condition number. Finally, in each stage, we use independent samples for $\log n$ iterations, but use same samples for the remaining iterations, there by eliminating the dependence of sample complexity on $\epsilon$.

Paper Organization: We first present the problem setup, our main result and an overview of our techniques in the next section. We then present a “warm-up” result for the basic SVP method in Section 3. We then present our main algorithm (St-SVP) and its analysis in Section 4. We conclude the discussion in Section 5. The proofs of all the technical lemmas will follow thereafter in the appendix.

Our Results and Techniques

In this section, we will first describe the problem set up and then present our results as well as the main techniques we use.

Let $\mathbf{M}$ be an $n_{1}\times n_{2}$ matrix of rank-$r$. Let $\Omega\subseteq[n_{1}]\times[n_{2}]$ be a subset of the indices. Recall that $P_{\Omega}\left(\mathbf{M}\right)$ (as defined in (2)) is the projection of $\mathbf{M}$ on to the indices in $\Omega$. Given $\Omega,\,P_{\Omega}\left(\mathbf{M}\right)$ and $r$, the goal is to recover $\mathbf{M}$. The problem is in general ill posed, so we make the following standard assumptions on $\mathbf{M}$ and $\Omega$ [CR09].

$\Omega$ is generated by sampling each element of $[n_{1}]\times[n_{2}]$ independently with probability $p$.

The incoherence assumption ensures that the mass of the matrix is well spread out and a small fraction of uniformly random observations give enough information about the matrix. Both of the above assumptions are standard and are used by most of the existing results, for instance [CR09, CT09, KMO10, Rec09, Kes12]. A few exceptions include the works of [MJD09, CBSW14, BJ14].

2 Main Result

The following theorem is the main result of this paper.

Suppose $\mathbf{M}$ and $\Omega$ satisfy Assumptions 1 and 2 respectively. Also, let

where $\alpha>1$, $n:=n_{1}+n_{2}$ and $C>0$ is a global constant. Then, the output $\mathbf{\widehat{M}}$ of Algorithm 2 satisfies: $\left\lVert\mathbf{\widehat{M}}-\mathbf{M}\right\rVert_{F}\leq\epsilon,$ with probability greater than $1-n^{-10-\log\alpha}$. Moreover, the run time of Algorithm 2 is $O\left({\left\lvert\Omega\right\rvert r^{2}\log(1/\epsilon)}\right)$.

Algorithm 2 is based on the projected gradient descent update (3) and proceeds in $r$ stages where in the $k$-th stage, projections are performed onto the set of rank-$k$ matrices. See Section 4 for a detailed description and the underlying intuition behing our algorithm.

Table 1 compares our result to that for nuclear norm minimization, which is the only other polynomial time method with finite sample complexity guarantees (i.e. no dependence on the desired accuracy $\epsilon$). Note that St-SVP runs in time near linear in the ambient dimension of the matrix ($n$), where as nuclear norm minimization runs in time cubic in the ambient dimension. However, the sample complexity of St-SVP is suboptimal in its dependence on the incoherence parameter $\mu$ and rank $r$. We believe closing this gap between the sample complexity of St-SVP and that of nuclear norm minimization should be possible and leave it for future work.

3 Overview of Techniques

In this section, we briefly present the key ideas and lemmas we use to prove Theorem 1. Our proof revolves around analyzing the basic SVP step (3): $\bm{X}_{t+1}=P_{k}\left(\bm{X}_{t}+\frac{1}{p}P_{\Omega}\left(\mathbf{M}-\bm{X}_{t}\right)\right)=P_{k}(\mathbf{M}+\widehat{\mathbf{H}})$ where $p$ is the sampling probability, $\widehat{\mathbf{H}}:=\bm{X}_{t}-\mathbf{M}-\frac{1}{p}P_{\Omega}\left(\bm{X}_{t}-\mathbf{M}\right)=\mathbf{E}-\frac{1}{p}P_{\Omega}(\mathbf{E})$ and $\mathbf{E}:=\bm{X}_{t}-\mathbf{M}$ is the error matrix. Hence, $\bm{X}_{t+1}$ is given by a rank-$k$ projection of $\mathbf{M}+\widehat{\mathbf{H}}$, which is a perturbation of the desired matrix $\mathbf{M}$.

In order to carry out this argument, we write the singular vectors of $\mathbf{M}+\widehat{\mathbf{H}}$ as solutions to eigenvector equations and then use these to write $\bm{X}_{t+1}$ explicitly via Taylor series expansion. We use this technique to prove the following more general lemma.

with probability greater than $1-n^{-10-\log\alpha}$.

Proceeding in stages: If we applied Lemma 1 with $k=r$, we would require $\left\lvert\beta\right\rvert$ to be much smaller than $\sigma_{r}$. Now, $\beta$ can be thought of as $\beta\approx\sqrt{\frac{n}{p}}\left\lVert\mathbf{E}\right\rVert_{\infty}$. If we start with $\bm{X}_{0}=0$, we have $\mathbf{E}=-\mathbf{M}$, and so $\left\lVert\mathbf{E}\right\rVert_{\infty}=\|\mathbf{M}\|_{\infty}\leq\frac{\sigma_{1}\mu^{2}r}{n}$. To make $\beta\leq\sigma_{r}$, we would need the sampling probability $p$ to be quadratic in the condition number $\kappa=\sigma_{1}/\sigma_{r}$ . In order to overcome this issue, we perform SVP in $r$ stages with the $k^{th}$ stage performing projections on to the set of rank-$k$ matrices while maintaining the invariant that at the end of $(k-1)^{\textrm{th}}$ stage, $\left\lVert\mathbf{E}\right\rVert_{\infty}=O(\sigma_{k}/n)$. This lets us choose a $p$ independent of $\kappa$ while still ensuring $\beta\approx\sqrt{\frac{n}{p}}\left\lVert\mathbf{E}\right\rVert_{\infty}\leq\sigma_{k}$. Lemma 1 tells us that at the end of the $k^{\textrm{th}}$ stage, the error $\left\lVert\mathbf{E}\right\rVert_{\infty}$ is $O\left({\frac{\sigma_{k+1}}{n}}\right)$, there by establishing the invariant for the $(k+1)^{\textrm{th}}$ stage.

Using same samples: In order to reduce the error from $O\left({\frac{\sigma_{k}}{n}}\right)$ to $O\left({\frac{\sigma_{k+1}}{n}}\right)$, the $k^{\textrm{th}}$ stage would require $O\left({\log\frac{\sigma_{k}}{\sigma_{k+1}}}\right)$ iterations. Since Lemma 1 requires the elements of $\mathbf{H}$ to be independent, in order to apply it, we need to use fresh samples in each iteration. This means that the sample complexity increases with $\frac{\sigma_{k}}{\sigma_{k+1}}$, or the desired accuracy $\epsilon$ if $\epsilon<\sigma_{k+1}$. This problem is faced by all the existing analysis for iterative algorithms for matrix completion [Kes12, JNS13, Har14, HW14]. We tackle this issue by observing that when $\mathbf{M}$ is ill conditioned and $\left\lVert\mathbf{E}\right\rVert_{F}$ is very small, we can show a decay in $\left\lVert\mathbf{E}\right\rVert_{F}$ using the same samples for SVP iterations:

Let $\mathbf{M}$ and $\Omega$ be as in Theorem 1 with $\mathbf{M}$ being a symmetric matrix. Further, let $\mathbf{M}$ be ill conditioned in the sense that $\left\lVert\mathbf{M}-P_{k}(\mathbf{M})\right\rVert_{F}<\frac{\sigma_{k}}{n^{3}}$, where $\sigma_{1}\geq\cdots\geq\sigma_{r}$ are the singular values of $\mathbf{M}$. Then, the following holds for all rank-$k$ $\bm{X}$ s.t. $\left\lVert\bm{X}-P_{k}(\mathbf{M})\right\rVert_{F}<\frac{\sigma_{k}}{n^{3}}$ (w.p. $\geq 1-n^{-10-\alpha}$):

The following lemma plays a crucial role in proving Lemma 2. It is a natural extension of the Davis-Kahan theorem for singular vector subspace perturbation.

Suppose $\bm{A}$ is a matrix such that $\sigma_{k+1}(\bm{A})\leq\frac{1}{4}\sigma_{k}(\bm{A})$. Then, for any matrix $\mathbf{E}$ such that $\left\lVert\mathbf{E}\right\rVert_{F}<\frac{1}{4}\sigma_{k}(\bm{A})$, we have:

In contrast to the Davis-Kahan theorem, which establishes a bound on the perturbation of the space of singular vectors, Lemma 3 establishes a bound on the perturbation of the best rank-$k$ approximation of a matrix $\bm{A}$ with good eigen gap, under small perturbations. This is a very natural quantity while considering perturbations of low rank approximations, and we believe it may find applications in other scenarios as well. A final remark regarding Lemma 3: we suspect it might be possible to tighten the right hand side of the result to $c\min\left(\sqrt{k}\left\lVert\mathbf{E}\right\rVert_{2},\left\lVert\mathbf{E}\right\rVert_{F}\right)$, but have not been able to prove it.

Singular Value Projection

As is clear from the pseudocode in Algorithm 1, SVP is a simple projected gradient descent method for solving the matrix completion problem. Note that Algorithm 1 first splits the set $\Omega$ into $T$ random subsets and updates iterate $\bm{X}_{t}$ using $\Omega_{t}$. This step is critical for analysis as it ensures that $\Omega_{t}$ is independent of $\bm{X}_{t-1}$, allowing for the use of standard tail bounds. The following theorem is our main result for Algorithm 1:

Suppose $\mathbf{M}$ and $\Omega$ satisfy Assumptions 1 and 2 respectively with

where $n=n_{1}+n_{2},\alpha>1,\kappa=\left(\frac{\sigma_{1}}{\sigma_{r}}\right)$ with $\sigma_{1}\geq\cdots\geq\sigma_{r}$ denoting the singular values of $\mathbf{M}$, $T=\log\frac{100\mu^{2}r\left\lVert\mathbf{M}\right\rVert_{2}}{\epsilon}$ and $C>0$ is a large enough global constant. Then, the output of Algorithm 1 satisfies (w.p. $\geq 1-T{\min(n_{1},n_{2})^{-10-\log\alpha}}$): $\left\lVert\bm{X}_{T}-\mathbf{M}\right\rVert_{F}\leq\epsilon$

Stagewise-SVP

Theorem 2 is suboptimal in its sample complexity dependence on the rank, condition number and desired accuracy. In this section, we will fix two of these issues – the dependence on condition number and desired accuracy – by designing a stagewise version of Algorithm 1 and proving Theorem 1.

Our algorithm, St-SVP (pseudocode presented in Algorithm 2) runs in $r$ stages, where in the $k^{\textrm{th}}$ stage, the projection is onto the set of rank-$k$ matrices. In each stage, the goal is to obtain an approximation of $\mathbf{M}$ up to an error of $\sigma_{k+1}$. In order to do this, we use the basic SVP updates, but in a very specific way, so as to avoid the dependence on condition number and desired accuracy.

(Step II) Determine if $\sigma_{k+1}>\frac{\sigma_{k}}{n^{3}}$: Note that we can determine this, by using the $(k+1)^{\textrm{th}}$ singular value of the matrix obtained after the gradient step, i.e., $\sigma_{k+1}(\bm{X}_{k,\log n}-\frac{1}{p}P_{\Omega_{k,\log n}}(\bm{X}_{k,\log n}-\mathbf{M}))$. If true, the error $\left\lVert\bm{X}_{k,\log n}-\mathbf{M}\right\rVert_{\infty}=O\left({\frac{\sigma_{k+1}}{n}}\right)$, and so the algorithm proceeds to the $(k+1)^{\textrm{th}}$ stage.

(Step III) If not (i.e., $\sigma_{k+1}\leq\frac{\sigma_{k}}{n^{3}}$), apply SVP update for $T=\log\frac{1}{\epsilon}$ iterations with same samples: If $\sigma_{k+1}\leq\frac{\sigma_{k}}{n^{3}}$, we can use Lemma 2 to conclude that after $\log\frac{1}{\epsilon}$ iterations, the Frobenius norm of error is $\left\lVert\bm{X}_{k,\log n+T}-\mathbf{M}\right\rVert_{F}=O\left({{n}\sigma_{k+1}+\epsilon}\right)$.

We will now present a proof of Theorem 1.

Just as in Theorem 2, it suffices to prove the result for when $\mathbf{M}$ is symmetric.For every stage, we will establish the following invariant:

We will use induction. (4) clearly holds for the base case $k=1$. Now, suppose (4) holds for the $k^{\textrm{th}}$ stage, we will prove that it holds for the $(k+1)^{\textrm{th}}$ stage. The analysis follows the four step outline in the previous section:

Step I: Here, we will show that for every iteration $t$, we have:

(5) holds for $t=0$ by our induction hypothesis (4) for the $k$-th stage. Supposing it true for iteration $t$, we will show it for iteration $t+1$. The $(t+1)^{\textrm{th}}$ iterate is given by:

This proves (5). Hence, after $\log n$ steps, we have:

Step II: Let $\bm{G}:=\bm{X}_{k,\log n}-\frac{1}{p}P_{\Omega_{k,\log n}}\left(\bm{X}_{k,\log n}-\mathbf{M}\right)=\mathbf{M}+\beta\mathbf{H}$ be the gradient update with notation as above. A standard perturbation argument (Lemmas 7 and 8) tells us that:

So if $\sigma_{k+1}(\bm{G})>\frac{\sigma_{k}(\bm{G})}{n^{3}}$, then we have $\sigma_{k+1}>\frac{9\sigma_{k}}{10n^{3}}$. Since we move on to the next stage with $\bm{X}_{k+1,0}=\bm{X}_{k,\log n}$, (7) tells us that:

showing the invariant for the $(k+1)^{\textrm{th}}$ stage.

Step III: On the other hand, if $\sigma_{k+1}(\bm{G})\leq\frac{\sigma_{k}(\bm{G})}{n^{3}}$, then Lemmas 8 and 8 tell us that $\sigma_{k+1}\leq\frac{11\sigma_{k}}{10n^{3}}$. So, using Lemma 2 with $T=\log\frac{1}{\epsilon}$ iterations, we obtain:

If $\epsilon>\frac{2}{{p}}\left\lVert\mathbf{M}-P_{k}\left(\mathbf{M}\right)\right\rVert_{F}$, then we have:

On the other hand, if $\epsilon\leq\frac{2}{{p}}\left\lVert\mathbf{M}-P_{k}\left(\mathbf{M}\right)\right\rVert_{F}$, then we have:

Step IV: Using (9) and “fresh samples” analysis as in Step I (in particular (5)), we have:

which establishes the invariant for the $(k+1)^{\textrm{th}}$ stage.

Combining the invariant (4) with the exit condition as in Step III, we have: $\|\widehat{\mathbf{M}}-\mathbf{M}\|_{F}\leq\epsilon$ where $\widehat{\mathbf{M}}$ is the output of the algorithm. As there are $r$ stages, and in each stage, we need $2\log n$ sets of samples of size $O(pn^{2})$. Hence, the total samplexity is $|\Omega|=O\left({\alpha\mu^{4}r^{5}n\log^{3}n}\right)$. Similarly, total computation complexity is $O\left({\alpha\mu^{4}r^{7}n\log^{3}n\log(\|M\|_{F}/\epsilon)}\right)$.∎

Discussion and Conclusions

In this paper, we proposed a fast projected gradient descent based algorithm for solving the matrix completion problem. The algorithm runs in time $O\left({nr^{7}\log^{3}n\log 1/\epsilon}\right)$, with a sample complexity of $O\left({nr^{5}\log^{3}n}\right)$. To the best our knowledge, this is the first near linear time algorithm for exact matrix completion with sample complexity independent of $\epsilon$ and condition number of $\mathbf{M}$.

Design an efficient algorithm with information-theoretic optimal sample complexity $\left\lvert\Omega\right\rvert=O\left({nr\log n}\right)$ is still open; our result is suboptimal by a factor of ${r^{4}\log^{2}n}$ and nuclear norm approach is suboptimal by a factor of $\log n$. Another interesting direction in this area is to design optimal algorithms that can handle sampling distributions that are widely observed in practice, such as the power law distribution[MJD09].

References

Appendix A Preliminaries and Notations for Proofs

The following lemma shows that wlog we can assume $\mathbf{M}$ to be a symmetric matrix. A similar result is given in Section D of [Har14].

Define the following symmetric matrix from $\mathbf{M}$ using a dilation technique:

Note that the rank of $\widetilde{\mathbf{M}}$ is $2\cdot r$ and the incoherence of $\widetilde{\mathbf{M}}$ is bounded by $(n_{1}+n_{2})/n_{2}\mu$ (assume $n_{1}\leq n_{2}$). Note that if $n_{2}>n_{1}$, then we can split the columns of $\mathbf{M}$ in blocks of size $n_{1}$ and apply the argument separately to each block.

Now, we can split $\Omega$ to generate samples from $\mathbf{M}$ and $\mathbf{M}^{T}$, and then augment redundant samples from the part above to obtain $\widetilde{\Omega}=[n]\times[n]$.

Moreover, if we run the SVP update (3) with input $\widetilde{\mathbf{M}}$, $\widetilde{\bm{X}}$ and $\widetilde{\Omega}$, an easy calculation shows that the iterates satisfy:

where $\bm{X}_{+}$ is the output of (3) with input $\mathbf{M}$, $\bm{X}$, and $\Omega$. That is, a convergence result for $\widetilde{\bm{X}}_{+}$ would imply a convergence result for $\bm{X}_{+}$ as well. ∎

Appendix B Proof of Lemma 1

$\mathbf{H}$ is a symmetric matrix with each of its elements drawn independently, satisfying the following moment conditions:

for $i,j\in[n]$ and $2\leq k\leq 2\log n$.

That is, we wish to understand $\|\bm{X_{+}}-\mathbf{M}\|_{\infty}$ under perturbation $\mathbf{H}$. To this end, we first present a few lemmas that analyze how $\mathbf{H}$ is obtained in the context of our St-SVP algorithm and also bounds certain key quantities related to $\mathbf{H}$. We then present a few technical lemmas that are helpful for our proof of Lemma 1. The detailed proof of the lemma is given in Section B.3. See Section B.4 for proofs of the technical lemmas.

Recall that the SVP update (3) is given by: $\bm{X_{+}}=P_{k}(\bm{X}-\frac{1}{p}P_{\Omega}(\bm{X}-\mathbf{M}))=P_{k}(\mathbf{M}+\mathbf{H})$ where $\mathbf{H}=\bm{E}-\frac{1}{p}P_{\Omega}(\bm{E})$ and $\bm{E}=\bm{X}-\mathbf{M}$. Our first lemma shows that matrices of the form $\bm{E}-\frac{1}{p}P_{\Omega}(\bm{E})$, scaled appropriately, satisfy Definition 1, i.e., satisfies the assumption of Lemma 1.

Let $\bm{A}$ be a symmetric $n\times n$ matrix. Suppose $\Omega\subseteq[n]\times[n]$ is obtained by sampling each element with probability $p\in\left[\frac{1}{4n},0.5\right]$. Then the matrix

We now present a critical lemma for our proof which bounds $\|H^{a}u\|_{\infty}$ for $2\leq a\leq\log n$. Note that the entries of $H^{a}$ can be dependent on each other, hence we cannot directly apply standard tail bounds. Our proof follows along very similar lines to Lemma 6.5 of [EKYY13]; see Appendix D for a detailed proof.

Suppose $\widehat{\mathbf{H}}$ satisfies Definition 1. Fix $1\leq a\leq\log n$. Let $\bm{e}_{r}$ denote the $r^{\textrm{th}}$ standard basis vector. Then, for any fixed vector $\bm{u}$, we have:

with probability greater than $1-n^{1-2\log\frac{c}{4}}$.

Next, we bound $\|H\|_{2}$ using matrix Bernstein inequality by [Tro12]; see Appendix B.4 for a proof.

Suppose $\mathbf{H}$ satisfies Definition 1. Then, w.p. $\geq 1-1/n^{10+\log\alpha}$, we have: $\left\lVert\mathbf{H}\right\rVert_{2}\leq 3\sqrt{\alpha}.$

B.2 Technical Lemmas useful for Proof of Lemma 1

In this section, we present the technical lemmas used by our proof of Lemma 1.

First, we present the well known Weyl’s perturbation inequality [Bha97]:

Suppose $\bm{B}=\bm{A}+\bm{N}$. Let $\lambda_{1},\cdots,\lambda_{n}$ and $\sigma_{1},\cdots,\sigma_{n}$ be the eigenvalues of $\bm{B}$ and $\bm{A}$ respectively. Then we have:

Next, we present a natural perturbation lemma that bounds the spectral norm distance of $A$ to $AB^{-1}A$ where $B=P_{k}(A+E)$ and $E$ is a perturbation to $A$.

B.3 Detailed Proof of Lemma 1

We are now ready to present a proof of Lemma 1. Recall that $\bm{X_{+}}=P_{k}(\mathbf{M}+\beta\mathbf{H})$, hence,

where $(\bm{u}_{i},\lambda_{i})$ is the $i^{\textrm{th}}\,(i\leq k)$ top eigenvector-eigenvalue pair (in terms of magnitude).

Now, as $\mathbf{H}$ satisfies conditions of Definition 1, we can apply Lemma 7 to obtain:

Using (10), we have: $\left(\mathbf{I}-\frac{\beta}{\lambda_{i}}\mathbf{H}\right)\bm{u}_{i}=\frac{1}{\lambda_{i}}\mathbf{M}\bm{u}_{i}$. Moreover, using (12), $\mathbf{I}-\frac{\beta}{\lambda_{i}}\mathbf{H}$ is invertible. Hence, using Taylor series expansion, we have:

Letting $\bm{U}\bm{\Lambda}{\bm{U}}^{\top}$ denote the eigenvalue decomposition (EVD) of $\bm{X_{+}}$, we obtain:

Using Lemma 9, we have the following bound for the first term above:

Let $\mathbf{M}=\mathbf{U}^{*}\mathbf{\Sigma}\left(\mathbf{U}^{*}\right)^{\top}$ denote the EVD of $\mathbf{M}$. We now bound the terms in the summation in (13) for $1\leq a+b<\log n$.

where $(\zeta_{1})$ follows from Lemma 6 and $(\zeta_{2})$ follows from (16).

where we used Lemma 10 to bound $\left\lVert\mathbf{M}\bm{U}{\bm{\Lambda}^{-(a+b+1)}}{\bm{U}}^{\top}{\mathbf{M}}\right\rVert_{2}$ and Lemma 7 to bound $\left\lVert\mathbf{H}\right\rVert_{2}$. The last inequality follows from using $(1/2)^{a+b}\leq 1/n\leq\frac{\mu^{2}r^{2}}{n}$ as $a+b>\log n$.

Plugging (17), (18) and (19) in (13) gives us:

B.4 Proofs of Technical Lemmas from Section B.1, Section B.2

The lemma now follows using matrix Bernstein inequality (Lemma 16). ∎

Let $\mathbf{M}=\bm{U}^{*}\mathbf{\Sigma}{\bm{U}^{*}}^{\top}$ be the eigenvalue decomposition $\mathbf{M}$. We have:

where $(\zeta_{1})$ follows from the incoherence of $\mathbf{M}$. ∎

Let $\bm{W}=\bm{U}\bm{\Lambda}{\bm{U}}^{\top}+\widetilde{\bm{U}}\widetilde{\bm{\Lambda}}{\widetilde{\bm{U}}}^{\top}$ be the eigenvalue decomposition of $\bm{W}$. Since $P_{k}(\bm{W})=\bm{U}\bm{\Lambda}{\bm{U}}^{\top}$, we see that $\left\lvert\lambda_{k}\right\rvert\geq\left\lvert\widetilde{\lambda}_{i}\right\rvert$.

Since $\left\lVert\bm{E}\right\rVert_{2}\leq\frac{\beta_{k}}{2}$, we see that

Using the eigenvalue decomposition of $\bm{W}$, we have the following expansion:

Applying triangle inequality and using $\left\lVert\bm{B}\bm{C}\right\rVert_{2}\leq\left\lVert\bm{B}\right\rVert_{2}\left\lVert\bm{C}\right\rVert_{2}$, we get:

Using the above inequality with (21), we obtain:

This proves the second claim of the lemma.

The last claim of the lemma follows by using triangle inequality and (21) in the above equation. ∎

Appendix C Proof of Lemma 2

We now present a proof of Lemma 2 that show decrease in the Frobenius norm of the error matrix, despite using same samples in each iteration. In order to state our proof, we will first introduce certain notations and provide a few perturbation results that might be of independent interest. Then, in next subsection, we will present a detailed proof of Lemma 2. Finally, in Section C.3, we present proofs of the technical lemmas given below.

In order to state our first supporting lemma, we will introduce the concept of tangent spaces of matrices [Bha97].

Let $\bm{A}$ be a matrix with EVD (eigenvalue decomposition) $\bm{U}^{*}\mathbf{\Sigma}{\bm{U}^{*}}^{\top}$. The following space of matrices is called the tangent space of $\bm{A}$:

That is, if $\bm{A}=\bm{U}^{*}\mathbf{\Sigma}{\bm{U}^{*}}^{\top}$ is the EVD of $\bm{A}$, then any matrix $\bm{B}$ can be decomposed into four mutually orthogonal terms as

where $\bm{U}^{*}_{\perp}$ is a basis of the orthogonal space of $\bm{U}^{*}$. The first three terms above are in ${\mathcal{T}}(\bm{A})$ and the last term is in ${{\mathcal{T}}(\bm{A})}^{\perp}$. We let ${\mathcal{P}}_{{\mathcal{T}}(\bm{A})}$ and ${\mathcal{P}}_{{\mathcal{T}}(\bm{A})^{\perp}}$ denote the projection operators onto ${\mathcal{T}}(\bm{A})$ and ${\mathcal{T}}(\bm{A})^{\perp}$ respectively.

Let $\bm{A}$ and $\bm{B}$ be two symmetric matrices. Suppose further that $\bm{B}$ is rank-$k$. Then, we have:

Next, we present a few technical lemmas related to norm of $M-P_{\Omega}(M)$:

Let $M$, $\Omega$ be as given in Lemma 2 and let $p=|\Omega|/n^{2}$ be the sampling probability. Then, For every $r\times r$ matrix $\mathbf{\widehat{\Sigma}}$, we have $(w.p.\geq 1-n^{-10-\alpha})$:

Let $M$, $\Omega$, $p$ be as given in Lemma 2. Then, for every $i,j\in[r]$, we have $(w.p.\geq 1-n^{-10-\alpha})$:

Let $M$, $\Omega$, $p$ be as given in Lemma 2. Then, for every $i,j\in[r]$ and $s\in[n]$, we have $(w.p.\geq 1-n^{-10-\alpha})$:

C.2 Detailed Proof of Lemma 2

Let $\mathbf{E}:={\bm{X}}-{P_{k}\left(\mathbf{M}\right)}$, $\mathbf{H}:=\mathbf{E}-\frac{1}{p}P_{\Omega}\left(\mathbf{E}\right)$ and $\bm{G}:=\bm{X}-\frac{1}{p}P_{\Omega}\left(\bm{X}-\mathbf{M}\right)=P_{k}\left(\mathbf{M}\right)+\mathbf{H}-\frac{1}{p}P_{\Omega}\left(\mathbf{M}-P_{k}\left(\mathbf{M}\right)\right)$. That is, $\bm{X_{+}}=P_{k}\left(\bm{G}\right)$.

For simplicity, in this section, we let $\mathbf{M}=\bm{U}^{*}\mathbf{\Sigma}{\bm{U}^{*}}^{\top}+\bm{U}^{*}_{\perp}\underline{\mathbf{\Sigma}}{\bm{U}^{*}_{\perp}}^{\top}$ denote the eigenvalue decomposition (EVD) of $\mathbf{M}$ with $P_{k}\left(\mathbf{M}\right)=\bm{U}^{*}\mathbf{\Sigma}{\bm{U}^{*}}^{\top}$, and also let $\underline{\mathbf{M}}=\bm{U}^{*}_{\perp}\underline{\mathbf{\Sigma}}{\bm{U}^{*}_{\perp}}^{\top}$. We also use the shorthand notation ${\mathcal{T}}:={\mathcal{T}}\left(P_{k}\left(\mathbf{M}\right)\right)$.

Representing $\bm{X}$ in terms of its projection onto ${\mathcal{T}}$ and its complement, we have:

where the last conclusion follows from Lemma 11 and the hypothesis that $\left\lVert\bm{X}-P_{k}\left(\mathbf{M}\right)\right\rVert_{F}<\frac{\left\lvert\sigma_{k}\right\rvert}{n^{2}}$. Using $\|\mathbf{E}\|_{F}\leq\sigma_{k}/n^{2}$, we have:

where we used the hypothesis that $\left\lVert\mathbf{M}-P_{k}\left(\mathbf{M}\right)\right\rVert_{F}<\frac{\sigma_{k}}{n^{2}}$ in the second inequality.

Since $\bm{X_{+}}=P_{k}\left(P_{k}\left(\mathbf{M}\right)+\mathbf{H}-\frac{1}{p}P_{\Omega}\left(\underline{\mathbf{M}}\right)\right)$, using Lemma 3 with (25), (26), we have:

Now, using Claim 1, we have $\left\lVert{\mathcal{P}}_{{\mathcal{T}}}\left({\mathbf{H}-\frac{1}{p}P_{\Omega}\left(\underline{\mathbf{M}}\right)}\right)\right\rVert_{F}<\frac{1}{10}\left\lVert P_{k}\left(\mathbf{M}\right)-\bm{X}\right\rVert_{F}+\frac{2}{\sqrt{p}}\left\lVert\underline{\mathbf{M}}\right\rVert_{F}$, which along with the above equation establishes the lemma. We now state and prove the claim bounding $\left\lVert{\mathcal{P}}_{{\mathcal{T}}}\left({\mathbf{H}-\frac{1}{p}P_{\Omega}\left(\underline{\mathbf{M}}\right)}\right)\right\rVert_{F}$ that we used above to finish the proof.

Assume notation defined in the section above. Then, we have:

We first bound $\left\lVert{\mathcal{P}}_{{\mathcal{T}}}\left({\mathbf{H}}\right)\right\rVert_{F}$. Recalling that $P_{k}\left(\mathbf{M}\right)=\bm{U}^{*}\mathbf{\Sigma}{\bm{U}^{*}}^{\top}$ is the EVD of $P_{k}\left(\mathbf{M}\right)$, we have:

Step I: To bound the first term in (27), we use Lemma 12 to obtain:

Step II: To bound the second term, we let $\bm{U}:=\bm{U}^{*}_{\perp}{\bm{\Lambda}_{1}}^{\top}$, and proceed as follows:

where $(\zeta_{1})$ follows from Lemma 13. This means that we can bound the second term as:

Step III: We now let $\bm{U}:=\bm{U}^{*}_{\perp}\bm{\Lambda}_{2}$ and turn to bound the third term in (27). We have:

Step IV: To bound the last term in (27), we use Lemma 11 to conclude

Combining (28), (29), (30) and (31), we have:

Claim now follows by combining (32) and (33). ∎

C.3 Proofs of Technical Lemmas from Section C.1

Let $B=\bm{U}\bm{\Lambda}{\bm{U}}^{\top}$ be EVD of $B$. Then, we have:

We will now prove Lemma 3, which is a natural extension of the Davis-Kahan theorem. In order to do so, we will first recall the Davis-Kahan theorem:

Let $\bm{A}=\bm{U}^{*}\mathbf{\Sigma}{\bm{U}^{*}}^{\top}+\bm{U}^{*}_{\perp}\mathbf{\widehat{\Sigma}}{\bm{U}^{*}_{\perp}}^{\top}$ be the EVD of $\bm{A}$ with $P_{k}\left(\bm{A}\right)=\bm{U}^{*}\mathbf{\Sigma}{\bm{U}^{*}}^{\top}$. Similarly, let $\bm{A}+\mathbf{E}=\bm{U}\bm{\Lambda}{\bm{U}}^{\top}+{\bm{U}}_{\perp}\bm{\widehat{\Lambda}}{{\bm{U}}_{\perp}}^{\top}$ denote the EVD of $\bm{A}+\mathbf{E}$ with $P_{k}\left(\bm{A}+\mathbf{E}\right)=\bm{U}\bm{\Lambda}{\bm{U}}^{\top}$. Expanding $P_{k}\left(\bm{A}+\mathbf{E}\right)$ into components along $\bm{U}^{*}$ and orthogonal to it, we have:

Before going on to bound the terms in (34), let us make some observations. We first use Lemma 8 to conclude that

Applying Theorem 3 with $S_{1}=\left[\frac{-\left\lvert\sigma_{k}\right\rvert}{2},\frac{\left\lvert\sigma_{k}\right\rvert}{2}\right]$ and $S_{2}=\left(-\infty,\frac{-3\left\lvert\sigma_{i}\right\rvert}{4}\right]\cup\left[\frac{3\left\lvert\sigma_{i}\right\rvert}{4},\infty\right)$, with separation parameter $\nu=\frac{\left\lvert\sigma_{i}\right\rvert}{4}$, we see that

We are now ready to bound the last two terms in the right hand side of (34). Firstly, we have:

where the last step follows from (36) and the assumption on $\left\lVert\mathbf{E}\right\rVert_{F}$. For the other term, we have:

where we used (35). Combining the above two inequalities with (34) proves the lemma. ∎

Finally, we present proofs for Lemma 12, Lemma 13, Lemma 14.

Using Theorem 1 by [BJ14], the followings $\forall\widehat{\Sigma}$ (w.p. $\geq 1-n^{-10-\alpha}$):

Lemma now follows by using the assumed value of $p$ in the above bound along with the fact that $\left(\bm{U}^{*}\mathbf{\widehat{\Sigma}}{\bm{U}^{*}}^{\top}-\frac{1}{p}P_{\Omega}\left(\bm{U}^{*}\mathbf{\widehat{\Sigma}}{\bm{U}^{*}}^{\top}\right)\right)\bm{U}^{*}$ is a rank-$r$ matrix. ∎

Let $H=\frac{1}{\beta}\left(\bm{u}^{*}_{j}{\bm{u}^{*}_{i}}^{\top}-\frac{1}{p}P_{\Omega}\left(\bm{u}^{*}_{j}{\bm{u}^{*}_{i}}^{\top}\right)\right)$, where $\beta=\frac{2\mu^{2}r}{\sqrt{n\cdot p}}$. Then, using Lemma 5, $H$ satisfies the conditions of Definition 1. Lemma now follows by using Lemma 7 and using $p$ as given in the lemma statement. ∎

Let $b_{i}$ be a set of independent bounded random variables, then the following holds $\forall\ t>0$:

Appendix D Proof of Lemma 6

We will prove the statement for $r=1$. The lemma can be proved by taking a union bound over all $r$. In order to prove the lemma, we will calculate a high order moment of the random variable

and then use Markov inequality. We use the following notation which is mostly consistent with Lemma $6.5$ of [EKYY13]. We abbreviate $(i,j)$ as $\alpha$ and denote $\widehat{h}_{ij}$ by $\widehat{h}_{\alpha}$. We further let

We now split the matrix $\widehat{\mathbf{H}}$ into two parts $\mathbf{H}$ and $\mathbf{H}^{\prime}$ which correspond to the upper triangular and lower triangular parts of $\widehat{\mathbf{H}}$. This means

The above summation has $2^{a}$ terms, of which we consider only

The resulting factor of $2^{a}$ does not change the result.

Abbreviating $\boldsymbol{\alpha}:=(\alpha_{1},\cdots,\alpha_{a})$, and

where the summation runs only over those $\boldsymbol{\alpha}$ such that $\alpha_{1}(1)=1$.

Calculating the $k^{\textrm{th}}$ moment expansion of $X_{a}$ for some even number $k$, we obtain:

For each valid $\boldsymbol{\alpha}=(\boldsymbol{\alpha}^{s})=(\alpha_{l}^{s})$, we define the partition $\Gamma(\boldsymbol{\alpha})$ of the index set $\left\{(s,l):s\in[k];l\in[a]\right\}$, where $(s,l)$ and $(s^{\prime},l^{\prime})$ are in the same equivalence class if $\alpha_{l}^{s}=\alpha_{l^{\prime}}^{s^{\prime}}$. We first bound the contribution of all $\boldsymbol{\alpha}$ corresponding to a partition $\Gamma$ in the summation (39) and then bound the total number of partitions $\Gamma$ possible. Since each $h_{\alpha}$ is centered, we can conclude that any partition $\Gamma$ that has a non-zero contribution to the summation in (39) satisfies:

each equivalence class of $\Gamma$ contains at least two elements.

We further bound the summation in (39) by taking absolute values of the summands

where the summation runs over $(\boldsymbol{\alpha}^{1},\cdots,\boldsymbol{\alpha}^{k})$ that correspond to valid partitions $\Gamma$. Fixing one such partition $\Gamma$, we bound the contribution to (40) of all the terms $\boldsymbol{\alpha}$ such that $\Gamma(\boldsymbol{\alpha})=\Gamma$.

We denote $G\equiv G(\Gamma)$ to be the graph constructed from $\Gamma$ as follows. The vertex set $V(G)$ is given by the equivalence classes of $\Gamma$. For every $(s,l)$, we have an edge between the equivalence class of $(s,l)$ and the equivalence class of $(s,l+1)$.

Each term in (40) can be bounded as follows:

where the last step follows from property $(*)$ above and Definition 1.

Using the above, we can bound (40) as follows:

where $v:=\left\lvert V(G)\right\rvert$ denotes the number of vertices in $G$.

Factorizing the above summation over different components of $G$, we obtain

where $l$ denotes the number of connected components of $G$, $G_{j}$ denotes the $j^{\textrm{th}}$ component of $G$, and $v_{j}$ denotes the number of vertices in $G_{j}$. We will now bound terms corresponding to one connected component at a time. Pick a connected component $G_{j}$. Since $\alpha_{1}^{s}(1)=1$ for every $s\in[a]$, we know that there exists a vertex $\alpha_{\gamma}\in G_{j}$ such that $\alpha_{\gamma}(1)=1$. Pick one such vertex as a root vertex and create a spanning tree $T_{j}$ of $G_{j}$. We use the bound $B_{\alpha_{\gamma}\alpha_{\gamma^{\prime}}}\leq 1$ for every $\left\{\gamma,\gamma^{\prime}\right\}\in E_{j}\setminus T_{j}$. The remaining summation $\sum_{\alpha_{1},\cdots,\alpha_{v_{j}}}\left(\prod_{\left\{\gamma,\gamma^{\prime}\right\}\in T_{j}}B_{\alpha_{\gamma}\alpha_{\gamma^{\prime}}}\right)$ can be calculated bottom up from leaves to the root. Since

Noting that the number of partitions $\Gamma$ is at most $(ka)^{ka}$, we obtain the bound

Choosing $k=2\lceil{\frac{\log n}{a}}\rceil$ and applying $k^{\textrm{th}}$ moment Markov inequality, we obtain

Applying a union bound now gives us the result.

Appendix E Empirical Results

In this section, we compare the performance of St-SVP with SVP on synthetic examples. We do not however include comparison to other matrix completion methods like nuclear norm minimization or alternating minimization; see [JMD10] for a comparison of SVP with those methods.

We implemented both the methods in Matlab and all the results are averaged over 5 random trials. In each trial we generate a random low rank matrix and observe $|\Omega|=5(n_{1}+n_{2})r\log(n_{1}+n_{2})$ entries from it uniformly at random.

In the first experiment, we fix the matrix size ($n_{1}=n_{2}=5000$) and generate random matrices with varying rank $r$. We choose the first singular value to be $1$ and the remaining ones to be $1/r$, giving us a condition number of $\kappa=r$. Figure 1 (a) & (b) show the error in recovery and the run time of the two methods, where we define the recovery error as $\left\lVert\widehat{\mathbf{M}}-\mathbf{M}\right\rVert_{2}/\left\lVert\mathbf{M}\right\rVert_{2}$. We see that St-SVP recovers the underlying matrix much more accurately as compared to SVP. Moreover, St-SVP is an order of magnitude faster than SVP.

In the next experiment, we vary the condition number of the generated matrices. Interestingly, for small $\kappa$, both SVP and St-SVP recover the underlying matrix in similar time. However, for larger $\kappa$, the running time of SVP increases significantly and is almost two orders of magnitude larger than that of St-SVP. Finally, we study the two methods with varying matrix sizes while keeping all the other parameters fixed ($r=10$, $\kappa=1/r$). Here again, St-SVP is much faster than SVP.