Low-Rank Matrix and Tensor Completion via Adaptive Sampling

Introduction

Recently, the machine learning and signal processing communities have focused considerable attention toward understanding the benefits of adaptive sensing. This theme is particularly relevant to modern data analysis, where adaptive sensing has emerged as an efficient alternative to obtaining and processing the large data sets associated with scientific investigation. These empirical observations have lead to a number of theoretical studies characterizing the performance gains offered by adaptive sensing over conventional, passive approaches. In this work, we continue in that direction and study the role of adaptive data acquisition in low rank matrix and tensor completion problems.

Our study is motivated not only by prior theoretical results in favor of adaptive sensing but also by several applications where adaptive sensing is feasible. In recommender systems, obtaining a measurement amounts to asking a user about an item, an interaction that has been deployed in production systems. Another application pertains to network tomography, where a network operator is interested in inferring latencies between hosts in a communication network while injecting few packets into the network. The operator, being in control of the network, can adaptively sample the matrix of pair-wise latencies, potentially reducing the total number of measurements. In particular, the operator can obtain full columns of the matrix by measuring from one host to all others, a sampling strategy we will exploit in this paper.

Yet another example centers around gene expression analysis, where the object of interest is a matrix of expression levels for various genes across a number of conditions. There are typically two types of measurements: low-throughput assays provide highly reliable measurements of single entries in this matrix while high-throughput microarrays provide expression levels of all genes of interest across operating conditions, thus revealing entire columns. The completion problem can be seen as a strategy for learning the expression matrix from both low- and high-throughput data while minimizing the total measurement cost.

We develop algorithms with theoretical guarantees for three low-rank completion problems. The algorithms find a small subset of columns of the matrix (tensor) that can be used to reconstruct or approximate the matrix (tensor). We exploit adaptivity to focus on highly informative columns, and this enables us to do away with the usual incoherence assumptions on the row-space while achieving competitive (or in some cases better) sample complexity bounds. Specifically our results are:

In the absence of noise, we develop a streaming algorithm that enjoys both low sample requirements and computational overhead. In the matrix case, we show that $\Omega(nr^{3/2}\log r)$ adaptively chosen samples are sufficient for exact recovery, improving on the best known bound of $\Omega(nr^{2}\log^{2}n)$ in the passive setting . This also gives the first guarantee for matrix completion with coherent row space.

In the tensor case, we establish that $\Omega(nr^{T-1/2}T^{2}\log r)$ adaptively chosen samples are sufficient for recovering a $n\times\ldots\times n$ order $T$ tensor of rank $r$. We complement this with a necessary condition for tensor completion under random sampling, showing that our adaptive strategy is competitive with any passive algorithm. These are the first sample complexity upper and lower bounds for exact tensor completion.

In the noisy matrix completion setting, we modify the adaptive column subset selection algorithm of Deshpande et al. to give an algorithm that finds a rank-$r$ approximation to a matrix using $\Omega(nr^{3/2}\textrm{polylog}(n))$ samples. As before, the algorithm does not require an incoherent row space but we are no longer able to process the matrix sequentially.

Along the way, we improve on existing results for subspace detection from missing data, the problem of testing if a partially observed vector lies in a known subspace.

Related Work

The matrix completion problem has received considerable attention in recent years. A series of papers , culminating in Recht’s elegent analysis of the nuclear norm minimization program, address the exact matrix completion problem through the framework of convex optimization, establishing that $\Omega((n_{1}+n_{2})r\max\{\mu_{0},\mu_{1}^{2}\}\log^{2}(n_{2}))$ randomly drawn samples are sufficient to exactly identify an $n_{1}\times n_{2}$ matrix with rank $r$. Here $\mu_{0}$ and $\mu_{1}$ are parameters characterizing the incoherence of the row and column spaces of the matrix, which we will define shortly. Candes and Tao proved that under random sampling $\Omega(n_{1}r\mu_{0}\log(n_{2}))$ samples are necessary, showing that nuclear norm minimization is near-optimal.

The noisy matrix completion problem has also received considerable attention . The majority of these results also involve some parameter that quantifies how much information a single observation reveals, in the same vein as incoherence.

Tensor completion, a natural generalization of matrix completion, is less studied. One challenge stems from the NP-hardness of computing most tensor decompositions, pushing researchers to study alternative structure-inducing norms in lieu of the nuclear norm . Both papers derive algorithms for tensor completion, but neither provide sample complexity bounds for the noiseless case.

Our approach involves adaptive data acquisition, and consequently our work is closely related to a number of papers focusing on using adaptive measurements to estimate a sparse vector . In these problems, specifically, problems where the sparsity basis is known a priori, we have a reasonable understanding of how adaptive sampling can lead to performance improvements. As a low rank matrix is sparse in its unknown eigenbasis, the completion problem is coupled with learning this basis, which poses a new challenge for adaptive sampling procedures.

Another relevant line of work stems from the matrix approximations literature. Broadly speaking, this research is concerned with efficiently computing a structured matrix, i.e. sparse or low rank, that serves as a good approximation to a fully observed input matrix. Two methods that apply to the missing data setting are the Nystrom method and entrywise subsampling . While the sample complexity bounds match ours, the analysis for the Nystrom method has focused on positive-semidefinite kernel matrices and requires incoherence of both the row and column spaces. On the other hand, entrywise subsampling is applicable, but the guarantees are weaker than ours.

It is also worth briefly mentioning the vast body of literature on column subset selection, the task of approximating a matrix by projecting it onto a few of its columns. While the best algorithms, namely volume sampling and sampling according to statistical leverages , do not seem to be readily applicable to the missing data setting, some algorithms are. Indeed our procedure for noisy matrix completion is an adaptation of an existing column subset selection procedure .

Our techniques are also closely related to ideas employed for subspace detection – testing whether a vector lies in a known subspace – and subspace tracking – learning a time-evolving low-dimensional subspace from vectors lying close to that subspace. Balzano et al. prove guarantees for subspace detection with known subspace and a partially observed vector, and we will improve on their result en route to establishing our guarantees. Subspace tracking from partial information has also been studied , but little is known theoretically about this problem.

Definitions and Preliminaries

where $e_{j}$ denotes the $j$th standard basis element.

In previous analyses of matrix completion, the incoherence assumption is that both the row and column spaces of the matrix have coherences upper bounded by $\mu_{0}$. When both spaces are incoherent, each entry of the matrix reveals roughly the same amount of information, so there is little to be gained from adaptive sampling, which typically involves looking for highly informative measurements. Thus the power of adaptivity for these problems should center around relaxing the incoherence assumption, which is the direction we take in this paper. Unfortunately, even under adaptive sampling, it is impossible to identify a rank one matrix that is zero in all but one entry without observing the entire matrix, implying that we cannot completely eliminate the assumption. Instead, we will retain incoherence on the column space, but remove the restrictions on the row space.

Exact Completion Problems

The pseudocode of the algorithm is given in Algorithm 1. Our first main result characterizes the performance of the tensor completion algorithm. We defer the proof to the appendix.

In the special case of a $n\times\ldots\times n$ tensor of order $T$, the algorithm succeeds with high probability using $\Omega(nr^{T-1/2}\mu_{0}^{T-1}T^{2}\log(Tr/\delta))$ samples, exhibiting a linear dependence on the tensor dimensions. In comparison, the only guarantee we are aware of shows that $\Omega\left(\left(\prod_{t=2}^{T_{1}}n_{t}\right)r\right)$ samples are sufficient for consistent estimation of a noisy tensor, exhibiting a much worse dependence on tensor dimension . In the noiseless scenario, one can unfold the tensor into a $n_{1}\times\prod_{t=2}^{T}n_{t}$ matrix and apply any matrix completion algorithm. Unfortunately, without exploiting the additional tensor structure, this approach will scale with $\prod_{t=2}^{T}n_{t}$, which is similarly much worse than our guarantee. Note that the naïve procedure that does not perform the recursive step has sample complexity scaling with the product of the dimensions and is therefore much worse than the our algorithm.

The most obvious specialization of Theorem 2 is to the matrix completion problem:

observations. The algorithm runs in $O(n_{1}n_{2}r+r^{3}m)$ time.

A few comments are in order. Recht guaranteed exact recovery for the nuclear norm minimization procedure as long as the number of observations exceeds $32(n_{1}+n_{2})r\max\{\mu_{0},\mu_{1}^{2}\}\beta\log^{2}(2n_{2})$ where $\beta$ controls the probability of failure and $||UV^{T}||_{\infty}\leq\mu_{1}\sqrt{r/(n_{1}n_{2})}$ with $\mu_{1}$ as another coherence parameter. Without additional assumptions, $\mu_{1}$ can be as large as $\mu_{0}\sqrt{r}$. In this case, our bound improves on his in its the dependence on $r,\mu_{0}$ and logarithmic terms.

The Nystrom method can also be applied to the matrix completion problem, albeit under non-uniform sampling. Given a PSD matrix, one uses a randomly sampled set of columns and the corresponding rows to approximate the remaining entries. Gittens showed that if one samples $O(r\log r)$ columns, then one can exactly reconstruct a rank $r$ matrix . This result requires incoherence of both row and column spaces, so it is more restrictive than ours. Almost all previous results for exact matrix completion require incoherence of both row and column spaces.

The one exception is a recent paper by Chen et al. that we became aware of while preparing the final version of this work . They show that sampling the matrix according to statistical leverages of the rows and columns can eliminate the need for incoherence assumptions. Specifically, when the matrix has incoherent column space, they show that by first estimating the leverages of the columns, sampling the matrix according to this distribution, and then solving the nuclear norm minimization program, one can recover the matrix with $\Omega(nr\mu_{0}\log^{2}n)$ samples. Our result improves on theirs when $r$ is small compared to $n$, specifically when $\sqrt{r}\log r\leq\log^{2}n$, which is common.

Our algorithm is also very computationally efficient. Existing algorithms involve successive singular value decompositions ($O(n_{1}n_{2}r)$ per iteration), resulting in much worse running times.

The key ingredient in our proofs is a result pertaining to subspace detection, the task of testing if a subsampled vector lies in a subspace. This result, which improves over the results of Balzano et al. , is crucial in obtaining our sample complexity bounds, and may be of independent interest.

Where $\alpha=\sqrt{2\frac{\mu(v)}{m}\log(1/\delta)}+2\frac{\mu(v)}{3m}\log(1/\delta)$, $\beta=6\log(d/\delta)+\frac{4}{3}\frac{d\mu(v)}{m}\log^{2}(d/\delta)$, $\gamma=\sqrt{\frac{8d\mu(U)}{3m}\log(2d/\delta)}$ and $\mu(v)=n||v||_{\infty}^{2}/||v||_{2}^{2}$.

This theorem shows that if $m=\Omega(\max\{\mu(v),d\mu(U),d\sqrt{\mu(U)\mu(v)}\}\log d)$ then the orthogonal projection from missing data is within a constant factor of the fully observed one. In contrast, Balzano et al. give a similar result that requires $m=\Omega(\max\{\mu(v)^{2},d\mu(U),d\mu(U)\mu(v)\}\log d)$ to get a constant factor approximation. In the matrix case, this improved dependence on incoherence parameters brings our sample complexity down from $nr^{2}\mu_{0}^{2}\log r$ to $nr^{3/2}\mu_{0}\log r$. We conjecture that this theorem can be further improved to eliminate another $\sqrt{r}$ factor from our final bound.

We adapt the proof strategy of Candes and Tao to the tensor completion problem and establish the following lower bound for uniform sampling:

Fix $1\leq m,r\leq\min_{t}n_{t}$ and $\mu_{0}>1$. Fix $0<\delta<1/2$ and suppose that we do not have the condition:

Theorem 5 implies that as long as the right hand side of Equation 6 is at most $\epsilon<1$, and:

then with probability at least $\delta$ there are infinitely many matrices that agree on the observed entries. This gives a necessary condition on the number of samples required for tensor completion. Note that when $T=2$ we recover the known lower bound for matrix completion.

Theorem 5 gives a necessary condition under uniform sampling. Comparing with Theorem 2 shows that our procedure outperforms any passive procedure in its dependence on the tensor dimensions. However, our guarantee is suboptimal in its dependence on $r$. The extra factor of $\sqrt{r}$ would be eliminated by a further improvement to Theorem 5, which we conjecture is indeed possible.

For adaptive sampling, one can obtain a lower bound via a parameter counting argument. Observing the $(i_{1},\ldots,i_{T})$th entry leads to a polynomial equation of the form $\sum_{k}\prod_{t}a_{k}^{(t)}(i_{t})=M_{i_{1},\ldots,i_{T}}$. If $m<r(\sum_{t}n_{t})$, this system is underdetermined showing that $\Omega((\sum_{t}n_{t})r)$ observations are necessary for exact recovery, even under adaptive sampling. Thus, our algorithm enjoys sample complexity with optimal dependence on matrix dimensions.

Noisy Matrix Completion

Our algorithm for noisy matrix completion is an adaptation of the column subset selection (CSS) algorithm analyzed by Deshpande et al. . The algorithm builds a candidate column space in rounds; at each round it samples additional columns with probability proportional to their projection on the orthogonal complement of the candidate column space.

To concretely describe the algorithm, suppose that at the beginning of the $l$th round we have a candidate subspace $U_{l}$. Then in the $l$th round, we draw $s$ additional columns according to the distribution where the probability of drawing the $i$th column is proportional to $||\mathcal{P}_{U_{l}^{\perp}}c_{i}||_{2}^{2}$. Observing these $s$ columns in full and then adding them to the subspace $U_{l}$ gives the candidate subspace $U_{l+1}$ for the next round. We initialize the algorithm with $U_{1}=\emptyset$. After $L$ rounds, we approximate each column $c$ with $\hat{c}=U_{L}(U_{L\Omega}^{T}U_{L\Omega})^{-1}U^{T}_{L\Omega}c_{\Omega}$ and concatenate these estimates to form $\hat{M}$.

The challenge is that the algorithm cannot compute the sampling probabilities without observing entries of the matrix. However, our results show that with reliable estimates, which can be computed from few observations, the algorithm still performs well.

Let $\Omega$ be the set of all observations over the course of the algorithm, let $U_{L}$ be the subspace obtained after $L=log(n_{1}n_{2})$ rounds and $\hat{M}$ be the matrix whose columns ${\hat{c}_{i}=U_{L}(U_{L\Omega}^{T}U_{L\Omega})^{-1}U_{L\Omega}^{T}c_{\Omega i}}$. Then there are constants $c_{1},c_{2}$ such that:

$\hat{M}$ can be computed from $\Omega((n_{1}+n_{2})r^{3/2}\mu(U)\textrm{polylog}(n_{1}n_{2}))$ observations. In particular, if $||A||_{F}^{2}=1$ and $R_{ij}\sim\mathcal{N}(0,\sigma^{2}/(n_{1}n_{2}))$, then there is a constant $c_{\star}$ for which:

The main improvement in the result is in relaxing the assumptions on the underlying matrix $A$. Existing results for noisy matrix completion require that the energy of the matrix is well spread out across both the rows and the columns (i.e. incoherence), and the sample complexity guarantees deteriorate significantly without such an assumption . As a concrete example, Negahban and Wainwright use a notion of spikiness, measured as $\sqrt{n_{1}n_{2}}\frac{||A||_{\infty}}{||A||_{F}}$ which can be as large as $\sqrt{n_{2}}$ in our setup, e.g. when the matrix is zero except for on one column and constant across that column.

The choices of $||A||_{F}^{2}=1$ and noise variance rescaled by $\frac{1}{n_{1}n_{2}}$ enable us to compare our results with related work . Thinking of $n_{1}=n_{2}=n$ and the incoherence parameter as a constant, our results imply consistent estimation as long as $\sigma^{2}=\omega\left(\frac{n}{r^{2}\textrm{polylog}(n)}\right)$. On the other hand, thinking of the spikiness parameter as a constant, show that the error is bounded by $\frac{\sigma^{2}nr\log n}{m}$ where $m$ is the total number of observations. Using the same number of samples as our procedure, their results implies consistency as long as $\sigma^{2}=\omega(r\textrm{polylog}(n))$. For small $r$ (i.e. $r=O(1)$), our noise tolerance is much better, but their results apply even with fewer observations, while ours do not.

Simulations

We verify Corollary 3’s linear dependence on $n$ in Figure 1, where we empirically compute the success probability of the algorithm for varying values of $n$ and $p=m/n$, the fraction of entries observed per column. Here we study square matrices of fixed rank $r=5$ with $\mu(U)=1$. Figure 1 shows that our algorithm can succeed with sampling a smaller and smaller fraction of entries as $n$ increases, as we expect from Corollary 3. In Figure 1, we instead plot success probability against total number of observations per column. The fact that the curves coincide suggests that the samples per column, $m$, is constant with respect to $n$, which is precisely what Corollary 3 implies. Finally, in Figure 1, we rescale instead by $n/\log^{2}n$, which corresponds to the passive sample complexity bound . Empirically, the fact that these curves do not line up demonstrates that our algorithm requires fewer than $\log^{2}n$ samples per column, outperforming the passive bound.

The second row of Figure 1 plots the same probability of success curves for the Singular Value Thresholding (SVT) algorithm . As is apparent from the plots, SVT does not enjoy a linear dependence on $n$; indeed Figure 1 confirms the logarithmic dependency that we expect for passive matrix completion, and establishes that our algorithm has empirically better performance.

In the third row, we study the algorithm’s dependence on $r$ on $500\times 500$ square matrices. In Figure 1 we plot the probability of success of the algorithm as a function of the sampling probability $p$ for matrices of various rank, and observe that the sample complexity increases with $r$. In Figure 1 we rescale the $x$-axis by $r^{-3/2}$ so that if our theorem is tight, the curves should coincide. In Figure 1 we instead rescale the $x$-axis by $r^{-1}$ corresponding to our conjecture about the performance of the algorithm. Indeed, the curves line up in Figure 1, demonstrating that empirically, the number of samples needed per column is linear in $r$ rather than the $r^{3/2}$ dependence in our theorem.

To confirm the computational improvement over existing methods, we ran our matrix completion algorithm on large-scale matrices, recording the running time and error in Table 3. To contrast with SVT, we refer the reader to Table 5.1 in . As an example, recovering a $10000\times 10000$ matrix of rank $100$ takes close to 2 hours with the SVT, while it takes less than 5 minutes with our algorithm.

For the noisy algorithm, we study the dependence on row-space incoherence. In Figure 3, we plot the reconstruction error as a function of the row space coherence for our procedure and the semidefinite program of Negahban and Wainwright , where we ensure that both algorithms use the same number of observations. It’s readily apparent that the SDP decays in performance as the row space becomes more coherent while the performance of our procedure is unaffected.

Conclusions and Open Problems

In this work, we demonstrate how sequential active algorithms can offer significant improvements in time, and measurement overhead over passive algorithms for matrix and tensor completion. We hope our work motivates further study of sequential active algorithms for machine learning.

Several interesting theoretical questions arise from our work:

Can we tighten the dependence on rank for these problems? In particular, can we bring the dependence on $r$ down from $r^{3/2}$ to linear? Simulations suggest this is possible.

Can one generalize the nuclear norm minimization program for matrix completion to the tensor completion setting while providing theoretical guarantees on sample complexity?

We hope to pursue these directions in future work.

References

Appendix A Proof of Corollary 3

Corollary 3 is considerably simpler to prove than Theorem 2, so we prove the former in its entirety before proceeding to the latter. To simplify the presentation, a number of technical lemmas regarding incoherence and concentration of measure are deferred to sections E and F, respectively.

When $m\geq 32r\mu_{0}\log(1/\delta)$. Here we used that $\mu(v)\leq r\mu(U)$ since $v\in\textrm{span}(U)$. For $\gamma$:

Which certainly holds when $m\geq 36r^{3/2}\mu_{0}\log(r/\delta)$, concluding the proof. ∎

So these columns are all recovered exactly. This step only adds a factor of $\delta$ to the failure probability, leading to the final term in the failure probability of the theorem.

Appendix B Proof of Theorem 2

We first focus on the recovery of the tensor in total, expressing this in terms of failure probabilities in the recursion. Then we inductively bound the failure probability of the entire algorithm. Finally, we compute the total number of observations. For now, define $\tau_{T}$ to be the failure probability of recovering a $T$-order tensor.

By Lemma 13, the subspace spanned by the mode-$T$ tensors has incoherence at most $r^{T-2}\mu_{0}^{T-1}$ and rank at most $r$ and each slice has incoherence at most $r^{T-1}\mu_{0}^{T-1}$. By the same argument as Lemma 7, we see that with $m\geq 36r^{T-1/2}\mu_{0}^{T-1}\log(2r/\delta)$ the projection test succeeds in identifying informative subtensors (those not in our current basis) with probability $\geq 1-4\delta$. With a union bound over these $r$ subtensors, the failure probability becomes $\leq 4r\delta+\delta$, not counting the probability that we fail in recovering these subtensors, which is $r\tau_{T-1}$.

For each order $T-1$ tensor that we have to recover, the subspace of interest has incoherence at most $r^{T-3}\mu^{T-2}$ and with probability $\geq 1-4r\delta$ we correctly identify each informative subtensor as long as $m\geq 36r^{T-3/2}\mu^{T-2}\log(2r/\delta)$. Again the failure probability is $\leq 4r\delta+\delta+r\tau_{T-2}$.

To compute the total failure probability we proceed inductively. $\tau_{1}=0$ since we completely observe any one-mode tensor (vector). The recurrence relation is:

We also compute the sample complexity inductively. Let $m_{T}$ denote the number of samples needed to complete a $T$-order tensor. Then $m_{1}=n_{1}$ and:

The running time is computed in a similar way to the matrix case. Assume that the running time to complete an order $t$ tensor is:

Note that this is exactly the running time of our Algorithm in the matrix case.

Per order $T-1$ subtensor, the projection and reconstructions take $O(r\prod_{t=1}^{T-1}n_{t})$, which in total contributes a factor of $O(r\prod_{t=1}^{T}n_{t})$. At most $r$ times, we must complete an order $T-1$ subtensor, and invert the matrix $U_{\Omega}^{T}U_{\Omega}$. These two together take in total:

Finally the cost of the Gram-schmidt process is $r^{2}\prod_{t=1}^{T-1}n_{t}$ which is dominated by the other costs. In total the running time is:

Appendix C Proof of Theorem 6

We will first prove a more general result and obtain Theorem 6 as a simple consequence.

Let $M=A+R$ where $A=U\Sigma V^{T}$ and $R_{ij}\sim\mathcal{N}(0,\sigma^{2})$. Let $M_{r}$ denote the best rank $r$ approximation to $M$. Assume that $A$ is rank $r$ and $\mu(U)\leq\mu_{0}$. For every $\delta,\epsilon\in(0,1)$ sample a set of size $s=\frac{5Lr}{2\delta\epsilon}$ at each of the $L$ rounds of the algorithm and compute $\hat{M}$ as prescribed. Then with probability $\geq 1-9\delta$:

and the algorithm has expected sample complexity:

The proof of this result involves some modifications to the analysis in . We will follow their proof, allowing for some error in the sampling probabilities, and arrive at a recovery guarantee. Then we will show how these sampling probabilities can be well-approximated from limited observations.

The first Lemma analyzes a single round of the algorithm, while the second gives an induction argument to chain the first across all of the rounds. These are extensions of Theorems 2.1 and Theorems 1.2, respectively, from .

Then with probability $\geq 1-\delta$ we have:

Where $\mathcal{P}_{H,r}$ denotes a projection on to the best $r$-dimensional subspace of $H$ and $M_{r}$ is the best rank $r$ approximation to $M$.

The proof closely mirrors that of Theorem 2.1 in . The main difference is that we are using an estimate of the correct distribution, and this will result in some additional error.

For each $i=1,\ldots,n_{2}$ and for each $l=1,\ldots s$ define:

That is the $i$th column of the residual $E$, scaled by the $i$th entry of the $j$th right singular vector, and the sampling probability. Defining $X^{(j)}=\frac{1}{s}\sum_{l=1}^{s}X_{l}^{(j)}$, we see that:

We will now proceed to bound the second central moment of $w^{(j)}$.

Now we use the probabilities $\hat{p}_{i}$ to evaluate each term in the summation:

This gives us an upper bound on the second central moment:

To complete the proof, let $y^{(j)}=w^{(j)}/\sigma_{j}$ and define the matrix $F=(\sum_{j=1}^{k}y^{(j)}u^{(j)T})M$. Since $y^{(j)}\in W$, the column space of $F$ is contained in $W$ so $||M-\mathcal{P}_{W}(M)||_{F}^{2}\leq||M-F||_{F}^{2}$.

We now use Markov’s inequality on the second term. Specifically, with probability $\geq 1-\delta$ we have:

Suppose that $(1+\alpha_{2})/(1-\alpha_{1})\leq c$ for some constant $c$ and for each of $L$ rounds of sampling. Let $S_{1},\ldots,S_{L}$ denote the sets of columns selected at each round and set $s=\frac{Lcr}{\delta\epsilon}$. Then with probability $\geq 1-\delta$ we have:

The proof is by induction on the number of rounds $L$. We will have each round of the algorithm fail with probability $\delta/L$ so that the total failure probability will be at most $\delta$. The base case follows from Lemma 9. At the $l$th round, the same lemma tells us:

Plugging in our choice of $s$ and the definition of $E$:

and applying the induction hypothesis we have:

To complete the proof, we just need to compute how many observations are necessary to ensure that $(1+\alpha_{2})/(1-\alpha_{1})\leq c$. We can do this by manipulating Theorem 4 and upper bounding the incoherences of the subspaces throughout the execution of the algorithm.

with probability $\geq 1-6\delta$ as long as the expected number of samples observed per column $m$ satisfies:

To establish the result, we will use the concentration results from Section F and the incoherence results form Section E. The goal will be to apply Theorem 4 with a union bound across all rounds and all columns, but we first need to bound the incoherences.

With a union bound, Lemma 14 shows that each column (once projected onto the orthogonal complement of one of the subspaces) has incoherence $O(r\mu(U)\log(n_{1}n_{2}L/\delta))$ with probability $\geq 1-\delta$. At the same time, Lemma 15 reveals that with probability $\geq 1-\delta$ all of the subspaces in the algorithm have incoherence at most $O(\mu(U)\log(n_{1}L/\delta))$.

By boosting the size of $m$ by a constant, we can make $\alpha\leq 1/4$. For $\gamma$ we have:

if we choose the constants correctly. Finally we have:

again using our definition of $m$. In particular, if we make this bound $\leq 1/4$ we then have that:

We are essentially done proving the theorem. The total number of samples used is:

We also completely observe $\Omega(L^{2}r/\delta\epsilon)$ columns. In total this gives us the sample complexity bound in Theorem 8. The failure probability is $\leq 7\delta$ ($6\delta$ from Lemma 11 and $\delta$ from Lemma 10).

So far we have recovered a subspace that can be used to approximate $M$. Unfortunately, we cannot actually compute $\mathcal{P}_{U_{L}}M$ given limited samples. Instead, for each column $c$, we compute $\hat{c}=U_{L}(U_{L\Omega}^{T}U_{L\Omega})^{-1}U_{\Omega L}c_{\Omega}$ and use $\hat{c}$ as our estimate of the column. This is similar to another projection operation, and the error will only be a constant factor worse than before.

Let $c_{i}$ denote a column of the matrix $M$ and let $\hat{U}$ denote the subspace at the end of the adaptive algorithm. Write $\hat{c}=\hat{U}(\hat{U}_{\Omega}^{T}\hat{U}_{\Omega})^{-1}\hat{U}_{\Omega}c$ Then with probability $\geq 1-2\delta$:

With $\beta$ and $\gamma$ defined as in Theorem 4.

Decompose $c=x+y$ where $x\in\hat{U}$ and $y\in\hat{U}^{\perp}$. It’s easy to see that $x=\hat{U}(\hat{U}_{\Omega}^{T}\hat{U}_{\Omega})^{-1}\hat{U}_{\Omega}x_{\Omega}$ so we are left with:

Because $y\in U^{\perp}$ so the cross term is zero. The second term here is equivalant to:

By Lemma 3 in the first term is upper bounded by $\frac{n_{1}^{2}}{(1-\gamma)^{2}m^{2}}$ while Lemma 17 reveals that the second term is upper bounded by $\beta\frac{m}{n_{1}^{2}}r\mu(\hat{U})||y||^{2}$. Combining these two yields the result. ∎

We already showed that with our choice of $m$, the expression in the above Lemma is smaller than $5/4$. Moreover the probability of failure simply adds $2\delta$ to the total failure probability. Thus:

and the last expression we bounded previously.

To prove the main theorem, it is best to view $M$ as equal to $A$ on all of the unobserved entries. In other words, if $\Omega$ is the set of all observations over the course of the algorithm, the random matrix $R$ is zero on $\Omega^{C}$. Since we never observed $M$ on $\Omega^{C}$, we have no way of knowing whether $M$ was equal to $A$ on those coordinates. It is therefore fair to write $M=A+R_{\Omega}$ where $R$ is zero on $\Omega^{C}$.

We expand the norm and then apply the main theorem:

Now since $M_{r}$ is the best rank $r$ approximation to $M$ (in Frobenius norm) and since $A$ is rank $r$, we know that $||M-M_{r}||_{F}\leq||M-A||_{F}$. With this substitution and setting $\epsilon=1/2,L=\log_{2}(n_{1}n_{2})$ we will arrive at the result (below constants are denoted by $c$ and they change from line to line):

which holds as long as $n_{1}n_{2}$ is sufficiently large.

Appendix D Proof of Theorem 5

We start by giving a proof in the matrix case, which is a slight variation of the proof by Candes and Tao . Then we turn to the tensor case, where only small adjustments are needed to establish the result. We work in the Bernoulli model, noting that Candes’ and Tao’s arguments demonstrate how to adapt these results to the uniform-at-random sampling model.

In the matrix case, suppose that $l_{1}=\frac{n_{1}}{r}$ and $l_{2}=\frac{n_{2}}{\mu_{0}r}$ are both integers. Define the following blocks $R_{1},\ldots R_{r}\subset[n_{1}]$ and $C_{1},\ldots C_{r}\subset[n_{2}]$ as:

Now consider the $n_{1}\times n_{2}$ family of matrices defined by:

$\mathcal{M}$ is a family of block-diagonal matrices where the blocks have size $l_{1}\times l_{2}$. Each block has constant rows whose entries may take arbitrary values in $[1,\sqrt{\mu_{0}}]$. For any $M\in\mathcal{M}$, the incoherence of the column space can be computed as:

A similar calculation reveals that the row space is also incoherent with parameter $\mu_{0}$.

Unique identification of $M$ is not possible unless we observe at least one entry from each row of each diagonal block. If we did not, then we could vary that corresponding coordinate in the appropriate $u_{k}$ and find infinitely many matrices $M^{\prime}\in\mathcal{M}$ that agree with our observations, have rank and incoherence at most $r$ and $\mu_{0}$ respectively. Thus, the probability of successful recovery is no larger than the probability of observing one entry of each row of each diagonal block.

The probability that any row of any block is unsampled is $\pi_{1}=(1-p)^{l_{2}}$ and the probability that all rows are sampled is $(1-\pi_{1})^{n_{1}}$. This must upper bound the success probability $1-\delta$. Thus:

or $\pi_{1}\leq 2\delta/n_{1}$ as long as $\delta<1/2$. Substituting $\pi_{1}=(1-p)^{l_{2}}$ we obtain:

as a necessary condition for unique identification of $M$.

Exponentiating both sides, writing $p=\frac{m}{n_{1}n_{2}}$ and the fact that $1-e^{-x}>x-x^{2}/2$ gives us:

when $\mu_{0}r/n_{2}\log(\frac{n_{1}}{2\delta})\leq\epsilon<1$.

D.2 Tensor Case

Fix $T$, the order of the tensor and suppose that $l_{1}=\frac{n_{1}}{r}$ is an integer. Moreover, suppose that $l_{t}=\frac{n_{t}}{\mu_{0}r}$ is an integer for $1<t\leq T$. Define a set of blocks, one for each mode and the family

This is a family of block-diagonal tensors and just as before, straightforward calculations reveal that each subspace is incoherent with parameter $\mu_{0}$. Again, unique identification is not possible unless we observe at least one entry from each row of each diagonal block. The difference is that in the tensor case, there are $\prod_{i\neq 1}l_{i}$ entries per row of each diagonal block so the probability that any single row is unsampled is $\pi_{1}=(1-p)^{\prod_{i\neq 1}l_{i}}$. Again there are $n_{1}$ rows and any algorithm that succeeds with probability $1-\delta$ must satisfy:

Which implies $\pi_{1}\leq 2\delta/n_{1}$ (assuming $\delta<1/2$). Substituting in the definition of $\pi_{1}$ we have:

The same approximations as before yield the bound (as long as $\frac{\mu_{0}^{T-1}r^{T-1}}{\prod_{i\neq j}n_{i}}\log(\frac{n_{1}}{2\delta})\leq\epsilon<1$):

Appendix E Properties about Incoherence

A significant portion of our proofs revolve around controlling incoherences of various subspaces used throughout the execution of the algorithms The following technical lemmas will enable us to work with these quantities.

$\mu(W_{1})\leq\frac{\textrm{dim}(U_{1})}{d^{\prime}}\mu(U_{1})$.

For the first property, since $W_{1}$ is a subspace of $U_{1}$, $\mathcal{P}_{W_{1}}e_{j}=\mathcal{P}_{W_{1}}\mathcal{P}_{U_{1}}e_{j}$ so $||\mathcal{P}_{W_{1}}e_{j}||_{2}^{2}\leq||\mathcal{P}_{U_{1}}e_{j}||_{2}^{2}$. The result now follows from the definition of incoherence.

For the second property, we instead compute the incoherence of:

Let $U$ be the column space of $M$ and let $V$ be some other subspace of dimension at most $n_{1}/2-k$. Let $v_{i}=\mathcal{P}_{V^{\perp}}c_{i}$ for each column $c_{i}$. Then with probability $\geq 1-\delta$:

Decompose $v_{i}=x_{i}+r_{i}$ where $x_{i}\in U\cap V^{\perp}$ and $r_{i}\in U^{\perp}\cap V^{\perp}$. Since each column is composed of a deterministic component living in $U$ and a random component, it must be the case that $r_{i}$ is a random gaussian vector living in $U^{\perp}\cap V^{\perp}$, which is a subspace of dimension at least $n_{1}-d\geq n_{1}-\textrm{dim}(U)-\textrm{dim}(V)$. We can now proceed with the bound:

For the second line, we used that $\frac{\sum_{i}a_{i}}{\sum_{i}b_{i}}\leq\sum_{i}\frac{a_{i}}{b_{i}}$ whenever $a_{i},b_{i}\geq 0$ which is the case here. Finally we use Lemma 19 on the denominator, 20 on the numerator, and a union bound over all $n_{2}$ columns. For $(n_{1}-d)$ sufficiently large (as long as $\sqrt{(n_{1}-d)\log n_{2}/\delta}\leq(n_{1}-d)/4$) and if $d\leq n_{1}/2$ we can bound as:

Let $I_{l}=\bigcup_{i=1}^{l}S_{i}$ and let $U_{l}=\textrm{span}(\{c_{i}\}_{i\in I_{l}})$ as in the execution of the noisy algorithm. If $|I_{l}|\leq n_{1}/2$ then with probability $\geq 1-\delta$, for all $l\in[L]$, we have:

It is clear that $U_{l}\subset\textrm{span}(\{c_{i}\}_{i\in I_{l}})\bigcup\textrm{span}(\{r_{i}\}_{i\in I_{l}})$ which will make things much easier to analyze. Note that $\textrm{span}(\{c_{i}\}_{i\in I_{l}})\subset U$ the original incoherent subspace and let $R_{I_{l}}$ denote the random matrix of columns corresponding to $I_{l}$. We then have:

Now if $|I_{l}|\leq n_{1}/2$ and $\delta$ is not exponentially small, the contribution from the random matrix is:

So the total incoherence will be (note that $\textrm{dim}(U_{l})=|I_{l}|$ with probability $1$ since $|I_{l}|\leq n_{1}/2$):

The failure probability here is $\delta n_{1}L$ if there are $L$ rounds so the incoherence is:

Appendix F A Collection of Concentration Results

We enumerate several concentration of measure lemmas that we use throughout our proofs. Many of these are well known results and we provide the references to their proofs.

We improve on the result of Balzano et al. to establish Theorem 4. The proof parallels theirs but with improvements to two key Lemmas. The improvement stems from using Bernstein’s inequality in lieu of standard Chernoff bounds in the concentration arguments and carries over into our sample complexity guarantees. Here we state and prove the two lemmas and then sketch the overal proof.

With the same notations as Theorem 4, with probability $\geq 1-2\delta$.

The difference between Lemma 16 and Lemma 1 from is in the definition of $\alpha$. Here we have reduced the relationship between $\mu(v)$ and $m$ from $\mu(v)^{2}/m$ to $\mu(v)/m$. The proof is an application of Bernstein’s inequality.

Let $X_{i}=v_{\Omega(i)}^{2}$ so that $\sum_{i=1}^{m}X_{i}=||v_{\Omega}||_{2}^{2}$. We can compute the variance and bound for $X_{i}$ as:

Finally plugging in the definition of $\alpha$ from the theorem shows that the right hand side is $\leq 2\delta$. ∎

In similar spirit to Lemma 16 we can also improve Lemma 2 from using Bernstein’s inequality:

With the same notations as Theorem 4, with probability at least $1-\delta$:

Again the improvement in our Lemma is in the expression $\beta$ where we have an improved dependence between $m$ and $\mu(y)$. The proof is an application of Bernstein’s inequality. Note that:

Where we have defined $X_{ji}=\sum_{k=1}^{n}u_{jk}v_{j}\mathbf{1}_{\Omega(i)=k}$. We have:

We apply Bernstein’s inequality and take a union bound, so that with probability $\geq 1-\delta$:

Notice also that $||u_{j}||_{\infty}^{2}\leq d\mu(U)/n$. Plugging in these bounds, with probability $\geq 1-\delta$:

Where we used that $||v||_{\infty}^{2}\leq||v||_{2}^{2}\mu(v)/n$ via the definition of incoherence. ∎

It will also be essential for these projections matrices to be invertible even with missing observations, as this will allow us to reconstruct columns of the matrix.

Let $\delta>0$ and $m\geq\frac{8}{3}r\mu_{0}\log(2r/\delta)$, Then:

with probability $\geq 1-\delta$, provided that $\gamma<1$. In particular $U_{\Omega}^{T}U_{\Omega}$ is invertible.

Let $W_{\Omega}^{T}W_{\Omega}=(U_{\Omega}^{T}U_{\Omega})^{-1}$. If Lemma 18 holds, $U_{\Omega}^{T}U_{\Omega}$ is invertible so

yields the lower bound. The upper bound follows from the same decomposition and Lemma 16. ∎

F.2 Concentration for Gaussian Vectors and Matrices

Let $X\sim\chi_{d}^{2}$. Then with probability $\geq 1-2\delta$:

A gaussian random vector $r$ has incoherence that depends on $||r||^{2}_{\infty}$ so it is crucial that we can control the maximum of gaussian random variables.

Let $X_{1},\ldots,X_{n}\sim\mathcal{N}(0,\sigma^{2})$. Then with probability $\geq 1-\delta$:

Finally, we will be projecting on to perturbed subspaces so we will need to control the coherence of these subspaces. The spectrum of the perturbation, will play a role in the coherence calculations.

Let $R$ be a $n\times t$ whose entries are independent standard normal random variables. Then for every $\epsilon\geq 0$, with probability $1-2\exp\{-\epsilon^{2}/2\}$, one has: