Tighter Low-rank Approximation via Sampling the Leveraged Element

Introduction

Finding a low-rank approximation to a matrix is fundamental to a wide array of machine learning techniques. The large sizes of modern data matrices has driven much recent work into efficient (typically randomized) methods to find low-rank approximations that do not exactly minimize the residual, but run much faster / parallel, with fewer passes over the data. Existing approaches involve either intelligent sampling of a few rows / columns of the matrix, projections onto lower-dimensional spaces, or sampling of entries followed by a top-$r$ SVD of the resulting matrix (with unsampled entries set to 0).

Both the sampling and the subsequent alternating minimization are naturally fast, parallelizable, and able to utilize sparsity in the input matrix. Existing literature has either focused on running in input sparsity time but approximation in (the weaker) Frobenius norm, or running in $O(n^{2})$ time with approximation in spectral norm. Our method provides the best of both worlds: it runs in input sparsity time, with just two passes over the data matrix, and yields an approximation in spectral norm. It does however have a dependence on the ratio of the first to the $r^{th}$ singular value of the matrix.

Our alternative approach also yields new methods for two related problems: directly finding the low-rank approximation of the product of two given matrices, and distributed PCA.

Our contributions are thus three new methods in this space:

Low-rank approximation of a general matrix: Our first (and main) contribution is a new method (LELA, Algorithm 1) for low-rank approximation of any given matrix: first draw a random subset of entries in a specific biased way, and then execute a weighted alternating minimization algorithm that minimizes the error on these samples over a factored form of the intended low-rank matrix. The sampling is done with only two passes over the matrix (each in input sparsity time), and both the sampling and the alternating minimization steps are highly parallelizable and compactly stored/manipulated.

For a matrix $M$, let $M_{r}$ be the best rank-$r$ approximation (i.e. the matrix corresponding to top $r$ components of SVD). Our algorithm finds a rank-$r$ matrix $\widehat{M}_{r}$ in time $O(nnz(M)+\frac{n\kappa^{2}r^{5}}{\epsilon^{2}})$, while providing approximation in spectral norm: $\|M-\widehat{M}_{r}\|\leq\|M-M_{r}\|+\epsilon\|M-M_{r}\|_{F}$, where $\kappa=\sigma_{1}(M)/\sigma_{r}(M)$ is the condition number of $M_{r}$. Existing methods either can run in input sparsity time, but provide approximations in (the weaker) Frobenius norm (i.e. with $||\cdot||$ replaced by $||\cdot||_{F}$ in the above expression), or run in $O(n^{2})$ time to approximate in spectral norm, but even then with leading constants larger than 1. Our method however does have a dependence on $\kappa$, which these do not. See Table 1 for a detailed comparison to existing results for low-rank approximation.

Direct approximation of a matrix product: We provide a new method to directly find a low-rank approximation to the product of two matrices, without having to first compute the product itself. To do so, we first choose a small set of entries (in a biased random way) of the product that we will compute, and then again run weighted alternating minimization on these samples. The choice of the biased random distribution is now different from above, as the sampling step does not have access to the product matrix. However, again both the sampling and alternating minimization are highly parallelizable.

Distributed PCA: Motivated by applications with really large matrices, recent work has looked at low-rank approximation in a distributed setting where there are $s$ servers – each have small set of rows of the matrix – each of which can communicate with a central processor charged with coordinating the algorithm. In this model, one is interested in find good approximations while minimizing both computations and the communication burden on the center.

Given a set $\Omega\subseteq[n]\times[d]$, $P_{\Omega}(M)$ is given by: $P_{\Omega}(M)(i,j)=M_{ij}$ if $(i,j)\in\Omega$ and otherwise. $R_{\Omega}(M)=w.*P_{\Omega}(M)$ denotes the Hadamard product of $w$ and $P_{\Omega}(M)$. That is, $R_{\Omega}(M)(i,j)=w_{ij}M_{ij}$ if $(i,j)\in\Omega$ and otherwise. Similarly let $R_{\Omega}^{1/2}(M)(i,j)=\sqrt{w_{ij}}M_{ij}$ if $(i,j)\in\Omega$ and otherwise.

Related results

Low rank approximation: Now we will briefly review some of the existing work on algorithms for low rank approximation. introduced the problem of computing low rank approximation of a matrix $M$ with few passes over $M$. They presented an algorithm that samples few rows and columns and does SVD to compute low rank approximation, and gave additive error guarantees. have extended these results. considered a different approach based on entrywise sampling and quantization for low rank approximation and has given additive error bounds.

have given low rank approximation algorithms with relative error guarantees in Frobenius norm. have provided guarantees on error in spectral norm which are later improved in . The main techniques of these algorithms is to use a random Gaussian or Hadamard transform matrix for projecting the matrix onto a low dimensional subspace and compute the rank-$r$ subspace. have given an algorithm based on random Hadamard transform that computes rank-$r$ approximation in time $O(\frac{n^{2}r}{\epsilon^{2}})$ and gives spectral norm bound of $\|M-\widehat{M}_{r}\|\leq c\|M-M_{r}\|+\epsilon\|M-M_{r}\|_{F}$.

One drawback of Hadamard transform is that it cannot take advantage of sparsity of the input matrix. Recently gave an algorithm using sparse subspace embedding that runs in input sparsity time with relative Frobenius norm error guarantees.

We presented some results in this area as a comparison with our results in table 1. This is a heavily subsampled set of existing results on low rank approximations. There is a lot of interesting work on very related problems of computing column/row based(CUR) decompositions, matrix sketching, low rank approximation with streaming data. Look at for more detailed discussion and comparison.

Matrix sparsification: In the matrix sparsification problem, the goal is to create a sparse sketch of a given matrix by sampling and reweighing the entries of the matrix. Various techniques for sampling have been proposed and analyzed which guarantee $\epsilon$ approximation error in Frobenius norm with $O(\frac{n}{\epsilon^{2}}\log n)$ samples . As we will see in the next section, the first step of algorithm 1 involves sampling according to a very specific distribution (similar to matrix sparsification), which has been designed for guaranteeing good error bounds for computing low rank approximation. For a comparison of various sampling distributions for the problem of low rank matrix recovery see .

Matrix completion: Matrix completion problem is to recover a $n\times n$ rank-$r$ matrix from observing small number of ($O(nr\log(n))$) random entries. Nuclear norm minimization is shown to recover the matrix from uniform random samples if the matrix is incoherentA $n\times d$ matrix $A$ of rank-$r$ with SVD $U^{*}\Sigma^{*}(V^{*})^{T}$ is incoherent if $\|(U^{*})^{i}\|^{2}\leq\frac{\mu_{0}r}{n},\forall i$ and $\|(V^{*})^{j}\|^{2}\leq\frac{\mu_{0}r}{d},\forall j$ for some constant $\mu_{0}$. . Similar results are shown for other algorithms like OptSpace and alternating minimization . Recently has given guarantees for recovery of any matrix under leverage score sampling from $O(nr\log^{2}(n))$ entries.

Distributed PCA: In distributed PCA, one wants to compute rank-$r$ approximation of a $n\times d$ matrix that is stored across $s$ servers with small communication between servers. One popular model is row partition model where subset of rows are stored at each server. Algorithms in achieve $O(\frac{dsr}{\epsilon})$ communication complexity with relative error guarantees in Frobenius norm, under this model. Recently have considered the scenario of arbitrary splitting of a $n\times d$ matrix and given an algorithm that has $O(\frac{dsr}{\epsilon})$ communication complexity with relative error guarantees in Frobenius norm.

Low-rank Approximation of Matrices

Sampling: A crucial ingredient of our approach is using the correct sampling distribution. Recent results in matrix completion indicate that a small number ($O(nr\log^{2}(n))$) of samples drawn in a way biased by leverage scoresIf SVD of $M_{r}=U^{*}\Sigma^{*}(V^{*})^{T}$ then leverage scores of $M_{r}$ are $||(U^{*})^{i}||^{2}$ and $||(V^{*})^{j}||^{2}$ for all $i,j$. can capture all the information in any exactly low-rank matrix. While this is indicative, here we have neither access to the leverage scores, nor is our matrix exactly low-rank. We approximate the leverage scores with the row and column norms ($||M^{i}||^{2}$ and $||M_{j}||^{2}$), and account for the arbitrary high-rank nature of input by including an $L_{1}$ term in the sampling; the distribution is given in eq. (2). Computationally, our sampling procedure can be done in two passes and $O(nnz(M)+m\log n)$ time.

Weighted alternating minimization: In our second step, we directly optimize over the factored form of the intended low-rank matrix, by minimizing a weighted squared error over the sampled elements from stage 1. That is, we first express the low-rank approximation $\widehat{M}_{r}$ as $UV^{T}$ and then iterate over $U$ and $V$ alternatingly to minimize the weighted $L_{2}$ error over the sampled entries (see Sub-routine 2). Note that this is different from taking principal components of a 0-filled version of the sampled matrix. The weights give higher emphasis to elements with smaller sampling probabilities. In particular, the goal is to minimize the following objective function:

That is, if $T=\log(\frac{\|M\|}{\epsilon\|M-M_{r}\|_{F}})$, we have:

Note that our time and sample complexity depends quadratically on $\kappa$. Recent results in the matrix completion literature shows that such a dependence can be improved to $\log(\kappa)$ by using a slightly more involved analysis . We believe a similar analysis can be combined with our techniques to obtain tighter bounds; we leave a similar tighter analysis for future research as such a proof would be significantly more tedious and would take away from the key message of this paper.

In the first step we take 2 passes over the matrix to compute the sampling distribution (2) and sampling the entries based on this distribution. It is easy to show that this step would require $O(nnz(M)+m\log(n))$ time. Next, the initialization step of WAltMin procedure requires computing rank-$r$ SVD of $R_{\Omega_{0}}(M)$ which has at most $m$ non-zero entries. Hence, the procedure can be completed in $O(mr)$ time using standard techniques like power method. Note that by Lemma 3.2 we need top-$r$ singular vectors of $R_{\Omega_{0}}(M)$ only upto constant approximation. Further $t$-th iteration of alternating minimization takes $O(2|\Omega_{2t+1}|r^{2})$ time. So, the total time complexity of our method is $O(nnz(M)+mr^{2})$. As shown in Theorem 3.1, our method requires $m=O(\frac{nr^{3}}{\epsilon^{2}}\kappa^{2}\log(n)\log^{2}(\frac{\|M\|}{\epsilon\|M-M_{r}\|_{F}}))$ samples. Hence, the total run-time of our algorithm is: $O(nnz(M)+\frac{nr^{5}}{\epsilon^{2}}\kappa^{2}\log(n)\log^{2}(\frac{\|M\|}{\epsilon\|M-M_{r}\|_{F}}))$.

Remarks: Now we will discuss how to sample entries of $M$ using sampling method (2) in $O(nnz(M)+m\log(n))$ time. Consider the following multinomial based sampling model: sample the number of elements per row (say $m_{i}$) by doing $m$ draws using a multinomial distribution over the rows, given by $\{0.5(\frac{d\|M^{i}\|^{2}}{(n+d)\|M\|_{F}^{2}}+\frac{1}{n+d})+0.5\frac{\|M^{i}\|_{1}}{\|M\|_{1,1}}\}$. Then, sample $m_{i}$ elements of the row-$i$, using $\{0.5\frac{\|M_{j}\|^{2}}{\|M\|_{F}^{2}}+0.5\frac{|M_{ij}|}{\|M\|_{1,1}}\}$ over $j\in[d]$, with replacement.

The failure probability in this model is bounded by 2 times the failure probability if the elements are sampled according to (2) . Hence, we can instead use the above mentioned multinomial model for sampling. Moreover, $\|M^{i}\|$, $\|M^{i}\|_{1}$ and $\|M_{j}\|$ can be computed in time $O(nnz(M)+n)$, so $m_{i}$’s can be sampled efficiently. Moreover, the multinomial distribution for all the rows can be computed in time $O(d+nnz(M))$, $O(d)$ work for setting up the first $\|M_{j}\|$ term and $nnz(M)$ term for changing the base distribution wherever $M_{ij}$ is non-zero. Hence, the total time complexity is $O(nnz(M)+m\log n)$.

2 Proof Overview:

We now present the key steps in our proof of Theorem 3.1. As mentioned in the previous section, our algorithm proceeds in two steps: entry-wise sampling of the given matrix $M$ and then weighted alternating minimization (WAltMin) to obtain a low-rank approximation of $M$.

Hence, the goal is to analyze the WAltMin procedure, with input samples obtained using (2), to obtain the bounds in Theorem 3.1. Now, WAltMin is an iterative procedure solving an inherently non-convex problem, $\min_{U,V}\sum_{(i,j)\in\Omega}w_{ij}(\mathbf{e}_{i}^{T}UV^{T}\mathbf{e}_{j}-M_{ij})^{2}$. Hence, it is prone to local minimas or worse, might not even converge. However, recent results for low-rank matrix completion have shown that alternating minimization (with appropriate initialization) can indeed be analyzed to obtain exact matrix completion guarantees.

Our proof also follows along similar lines, where we show that the initialization procedure (step 4 of Sub-routine 2) provides an accurate enough estimate of $M_{r}$ and then at each step, we show a geometric decrease in distance to $M_{r}$. However, our proof differs from the previous works in two key aspects: a) existing proof techniques of alternating minimization assume that each element is sampled uniformly at random, while we can allow biased and approximate sampling, b) existing techniques crucially use the assumption that $M_{r}$ is incoherent, while our proof avoids this assumption using the weighted version of AltMin.

We now present our bounds for initialization as well as for each step of the WAltMin procedure. Theorem 3.1 follows easily from the two bounds.

Let the set of entries $\Omega$ be generated according to (2). Also, let $m\geq C\frac{n}{\delta^{2}}\log(n)$. Then, the following holds (w.p. $\geq 1-\frac{2}{n^{10}}$):

Also, if $\|M-M_{r}\|_{F}\leq\frac{1}{576\kappa r^{1.5}}\|M_{r}\|_{F}$, then the following holds (w.p. $\geq 1-\frac{2}{n^{10}}$):

where $\widehat{U}^{(0)}$ is the initial iterate obtained using Steps 4, 5 of Sub-Procedure 2. $\kappa=\sigma_{1}^{*}/\sigma_{r}^{*}$, $\sigma_{i}^{*}$ is the $i$-th singular value of $M$, $M_{r}=U^{*}\Sigma^{*}(V^{*})^{T}$.

Let $\mathcal{P}_{r}(A)$ be the best rank-$r$ approximation of $A$. Then, Lemma 3.2 and Weyl’s inequality implies that:

Now, we can have two cases: 1) $\|M-M_{r}\|_{F}\geq\frac{1}{576\kappa r^{1.5}}\|M_{r}\|_{F}$: In this case, setting $\delta=\epsilon/(\kappa r^{1.5})$ in (4) already implies the required error bounds of Theorem 3.1There is a small technicality here: alternating minimization can potentially worsen this bound. But the error after each step of alternating minimization can be effectively checked using a small cross-validation set and we can stop if the error increases.. 2) $\|M-M_{r}\|_{F}\leq\frac{1}{576\kappa r^{1.5}}\|M_{r}\|_{F}$. In this regime, we will show now that alternating minimization reduces the error from initial $\delta\|M\|_{F}$ to $\delta\|M-M_{r}\|_{F}$.

Let hypotheses of Theorem 3.1 hold. Also, let $\|M-M_{r}\|_{F}\leq\frac{1}{576\kappa r\sqrt{r}}\|M_{r}\|_{F}$. Let $\widehat{U}^{(t)}$ be the $t$-th step iterate of Sub-Procedure 2 (called from Algorithm 1), and let $\widehat{V}^{(t+1)}$ be the $(t+1)$-th iterate (for $V$). Also, let $\|(U^{(t)})^{i}\|\leq 8\sqrt{r}\kappa\sqrt{\frac{\|M_{j}\|^{2}}{\|M\|_{F}^{2}}+\frac{|M_{ij}|}{\|M\|_{F}}}$ and $dist({U}^{(t)},U^{*})\leq\frac{1}{2}$, where $U^{(t)}$ is a set of orthonormal vectors spanning $\widehat{U}^{(t)}$. Then, the following holds (w.p. $\geq 1-\gamma/T$):

and $\|(V^{(t+1)})^{j}\|\leq 8\sqrt{r}\kappa\sqrt{\frac{\|M_{j}\|^{2}}{\|M\|_{F}^{2}}+\frac{|M_{ij}|}{\|M\|_{F}}}$, where $V^{(t+1)}$ is a set of orthonormal vectors spanning $\widehat{V}^{(t+1)}$.

The above lemma shows that distance between $\widehat{V}^{(t+1)}$ and $V^{*}$ (and similarly, $\widehat{U}^{(t+1)}$ and $U^{*}$) decreases geometrically up to $\epsilon\|M-M_{r}\|_{F}/\sigma^{*}_{r}$. Hence, after $\log(\|M\|_{F}\|/\zeta)$ steps, the first error term in the bounds above vanishes and the error bound given in Theorem 3.1 is obtained.

Note that, the sampling distribution used for our result is a “hybrid” distribution combining leverage scores and the $L_{1}$-style sampling. However, if $M$ is indeed a rank-$r$ matrix, then our analysis can be extended to handle the leverage score based sampling itself ($q_{ij}=m\cdot\frac{\|M^{i}\|^{2}+\|M_{j}\|^{2}}{2n\|M\|_{F}^{2}}$). Hence our results also show that weighted alternating minimization can be used to solve the coherent-matrix completion problem introduced in .

3 Direct Low-rank Approximation of Matrix Product

In this section we present a new pass efficient algorithm for the following problem: suppose we are given two matrices, and desire a low-rank approximation of their product $AB$; in particular, we are not interested in the actual full matrix product itself (as this may be unwieldy to store and use, and thus wasteful to produce in its entirety). One example setting where this arises is when one wants to calculate the joint counts between two very large sets of entities; for example, web companies routinely come across settings where they need to understand (for example) how many users both searched for a particular query and clicked on a particular advertisement. The number of possible queries and ads is huge, and finding this co-occurrence matrix from user logs involves multiplying two matrices – query-by-user and user-by-ad respectively – each of which is itself large.

We give a method that directly produces a low-rank approximation of the final product, and involves storage and manipulation of only the efficient factored form (i.e. one tall and one fat matrix) of the final intended low-rank matrix. Note that as opposed to the previous section, the matrix does not already exist and hence we do not have access to its row and column norms; so we need a new sampling scheme (and a different proof of correctness).

Algorithm: Suppose we are given an $n_{1}\times d$ matrix $A$ and another $d\times n_{2}$ matrix $B$, and we wish to calculate a rank-$r$ approximation of the product $A\cdot B$. Our algorithm proceeds in two stages:

Choose a biased random set $\Omega\subset[n_{1}]\times[n_{2}]$ of elements as follows: choose an intended number $m$ (according to Theorem 3.4 below) of sampled elements, and then independently include each $(i,j)\in[n_{1}]\times[n_{2}]$ in $\Omega$ with probability given by $\hat{q}_{ij}=\min\{1,q_{ij}\}$ where

Then, find $P_{\Omega}(A\cdot B)$, i.e. only the elements of the product $AB$ that are in this set $\Omega$.

Run the alternating minimization procedure WAltMin$(P_{\Omega}(A\cdot B),\Omega,r,\hat{q},T)$, where $T$ is the number of iterations (again chosen according to Theorem 3.4 below). This produces the low-rank approximation in factored form.

Remarks: Note that the sampling distribution now depends only on the row norms $\|A^{i}\|^{2}$ of $A$ and the column norms $\|B_{j}\|^{2}$ of $B$; each of these can be found completely in parallel, with one pass over each row/column of the matrices $A$ / $B$. A second pass, again parallelizable, calculates the element $(A\cdot B)_{ij}$ of the product, for $(i,j)\in\Omega$. Once this is done, we are again in the setting of doing weighted alternating minimization over a small set of samples – the setting we had before, and as already mentioned this too is highly parallelizable and very fast overall. In particular, the computation complexity of the algorithm is $O(|\Omega|\cdot(d+r^{2}))=O(m(d+r^{2}))=O(\frac{nr^{3}\kappa^{2}}{\epsilon^{2}}\cdot(d+r^{2}))$ (suppressing terms dependent on norms of $A$ and $B$ ), where $n=\max\{n_{1},n_{2}\}$.

We now present our theorem on the number of samples and iterations needed to make this procedure work with at least a constant probability.

Further when $\|Y-Y_{r}\|_{F}\leq\|Y_{r}\|_{F}$ we get:

Distributed Principal Component Analysis

Modern large-scale systems have to routinely compute PCA of data matrices with millions of data points embedded in similarly large number of dimensions. Now, even storing such matrices on a single machine is not possible and hence most industrial scale systems use distributed computing environment to handle such problems. However, performance of such systems depend not only on computation and storage complexity, but also on the required amount of communication between different servers.

In contrast, a distributed setting extension of our LELA algorithm 1 has linear communication complexity $O(ds+\frac{nr^{5}}{\epsilon^{2}})$ and computes rank-$r$ approximation $\widehat{M}_{r}$, with $||M-\widehat{M}_{r}||\leq||M-M_{r}||+\epsilon||M-M_{r}||_{F}$. Now note that if $n\approx d$ and if $s$ scales with $n$ (which is a typical requirement), then our communication complexity can be significantly better than that of . Moreover, our method provides spectral norm bounds as compared to relatively weak Frobenius bounds mentioned above.

Algorithm: The distributed version of our LELA algorithm depends crucially on the following observation: given $V$, each row of $U$ can be updated independently. Hence, servers need to communicate rows of $V$ only. There also, we can use the fact that each server requires only $O(nr/s\log n)$ rows of $V$ to update their corresponding $U_{\mathbf{r_{k}}}$. $U_{\mathbf{r_{k}}}$ denote restriction of ${U}$ to rows in set $\{r_{k}\}$ and outside and similarly $\widehat{V}_{\mathbf{c_{k}}},\widehat{Y}_{\mathbf{c_{k}}}$ denote restriction of $V,\widehat{Y}$ to rows $\{c_{k}\}$.

We now describe the distributed version of each of the critical step of LELA algorithm. See Algorithm 3 for a detailed pseudo-code. For simplicity, we dropped the use of different set of samples in each iteration. Correspondingly the algorithm will modify to distributing samples $\Omega^{k}$ into $2T+1$ buckets and using one in each iteration. This simplification doesn’t change the communication complexity.

Sampling: For sampling, we first compute column norms $\|M_{j}\|,\forall 1\leq j\leq d$ and communicate to each server. This operation would require $O(ds)$ communication. Next, each server (server $k$) samples elements from its rows $\{\mathbf{r_{k}}\}$ and stores $R_{\Omega^{k}}(M_{\mathbf{r_{k}}})$ locally. Note that because of independence over rows, the servers don’t need to transmit their samples to other servers.

Initialization: In the initialization step, our algorithm computes top $r$ right singular vector of $R_{\Omega}(M)$ by iterations $\widehat{Y}^{(t+1)}=R_{\Omega}(M)^{T}R_{\Omega}(M)Y^{t}=\sum_{k}R_{\Omega^{k}}(M)^{T}R_{\Omega^{k}}(M)Y^{t}$. Now, note that computing $R_{\Omega^{k}}(M)Y^{t}$ requires server $k$ to access atmost $|\Omega^{k}|$ columns of $Y^{t}$. Hence, the total communication from the CP to all the servers in this round is $O(|\Omega|r)$. Similarly, each column of $R_{\Omega^{k}}(M)^{T}R_{\Omega^{k}}(M)Y^{t}$ is only $|\Omega^{k}|$ sparse. Hence, total communication from all the servers to CP in this round is $O(|\Omega|r)$. Now, we need constant many rounds to get a constant factor approximation to SVD of $R_{\Omega}(M)$, which is enough for good initialization in WAltMin procedure. Hence, total communication complexity of the initialization step would be $O(|\Omega|r)$.

Alternating Minimization Step: For alternating minimization, update to rows of $U$ is computed at the corresponding servers and the update to $V$ is computed at the CP. For updating $\widehat{U}^{(t+1)}_{\mathbf{r_{k}}}$ at server $k$, we use the following observation: updating $\widehat{U}^{(t+1)}_{\mathbf{r_{k}}}$ requires atmost $|\Omega^{k}|$ rows of $\widehat{V}^{(t)}$. Hence, the total communication from CP to all the servers in the $t$-th iteration is $O(|\Omega|r)$. Next, we make a critical observation that update $\widehat{V}^{(t+1)}$ can be computed by adding certain messages from each server (see Algorithm 3 for more details). Message from server $k$ to CP is of size $O(|\Omega^{k}|r^{2})$. Hence, total communication complexity in each round is $O(|\Omega|r^{2})$ and total number of rounds is $O(\log(\|M\|_{F}/\zeta))$.

We now combine the above given observations to provide error bounds and communication complexity of our distributed PCA algorithm:

Let the $n\times d$ matrix $M$ be distributed over $s$ servers according to the row-partition model. Let $m\geq\frac{C}{\gamma}\frac{nr^{3}}{\epsilon^{2}}\kappa^{2}\log(n)\log^{2}(\frac{||M_{r}||}{\zeta})$. Then, the algorithm 3 on completion will leave matrices $\widehat{U}_{\mathbf{r_{k}}}^{(t+1)}$ at server $k$ and $\widehat{V}^{(t+1)}$ at CP such that the following holds (w.p. $\geq 1-\gamma$): $||M-\widehat{U}^{(t+1)}(\widehat{V}^{(t+1)})^{T}||\leq||M-M_{r}||+\epsilon||M-M_{r}||_{F}+\zeta$, where $\widehat{U}^{(t+1)}=\sum_{k}\widehat{U}_{\mathbf{r_{k}}}^{(t+1)}$. This algorithm has a communication complexity of $O(ds+|\Omega|r^{2})=O(ds+\frac{nr^{5}\kappa^{2}}{\epsilon^{2}}\log^{2}(\frac{||M_{r}||}{\zeta}))$ real numbers.

As discussed above, each update to $\widehat{V}^{(t)}$ and $\widehat{U}^{(t)}$ are computed exactly as given in the WAltMin procedure (Sub-routine 2). Hence, error bounds for the algorithm follows directly from Theorem 3.1. Communication complexity bounds follows by observing that $|\Omega|\leq 2m$ w.h.p.

Remark: The sampling step given above suggests another simple algorithm where we can compute $P_{\Omega}(M)$ in a distributed fashion and communicate the samples to CP. All the computation is performed at CP afterwards. Hence, the total communication complexity would be $O(ds+|\Omega|)=O(ds+\frac{nr^{3}\kappa^{2}}{\epsilon^{2}}\log(\frac{||M_{r}||}{\zeta})$, which is lesser than the communication complexity of Algorithm 3. However, such an algorithm is not desirable in practice, because it is completely reliant on one single server to perform all the computation. Hence it is slower and is fault-prone. In contrast, our algorithm can be implemented in a peer-peer scenario as well and is more fault-tolerant.

Also, the communication complexity bound of Theorem 4.1 only bounds the total number of real numbers transferred. However, if each of the number requires several bits to communicate then the real communication can still be very large. Below, we bound each of the real number that we transfer, hence providing a bound on the number of bits transferred.

Bit complexity: First we will bound $w_{ij}M_{ij}$. Note that we need to bound this only for $(i,j)\in\Omega$. Now, $|w_{ij}M_{ij}|\leq||R_{\Omega}(M)||_{\infty}\leq||R_{\Omega}(M)||\leq||M||+\epsilon||M||_{F}\leq 2*ndM_{max}$, where the third inequality follows from Lemma 3.2. Hence if the entries of the matrix $M$ are being represented using $b$ bits initially then the algorithm needs to use $O(b+\log(nd))$ bits. By the same argument we get a bound of $O(b+\log(nd))$ bits for computing $||M^{i}||^{2},\forall i;||M_{j}||^{2},\forall j;||M||_{F}^{2}$ and $||M||_{1,1}$.

Further at any stage of the WAltMin iterations $||\widehat{U}^{(t)}(\widehat{V}^{(t+1)})^{T}||_{\infty}\leq||\widehat{U}^{(t)}(\widehat{V}^{(t+1)})^{T}||\leq 2||M||_{F}$. So this stage also needs $O(b+\log(n))$ bits for computation. Hence overall the bit complexity of each of the real numbers of the algorithm 3 is $O(b+\log(nd))$ , if $b$ bits are needed for representing the matrix entries. That is, overall communication complexity of the algorithm is $O((b+\log(nd))\cdot(ds+\frac{nr^{3}\kappa^{2}}{\epsilon^{2}}\log(\frac{||M_{r}||}{\zeta}))$.

Simulations

In this section we present some simulation results on synthetic data to show the error performance of the algorithm 1. First we consider the setting of finding low rank approximation of a given matrix $M$. Later we consider the setting of computing low rank approximation of $A\cdot B$, given $A$ and $B$ without computing the product.

For simulations we consider random matrices of size 1000 by 1000 and rank-5. $M_{r}$ is a rank 5 matrix with all singular values 1. We consider two cases, one in which $M_{r}$ is incoherent and other in which $M_{r}$ is coherent. Recall that a $n\times d$ rank-$r$ matrix $M_{r}$ is an incoherent matrix if $\|(U^{*})^{i}\|^{2}\leq\frac{\mu_{0}r}{n},\forall i$ and $\|(V^{*})^{j}\|^{2}\leq\frac{\mu_{0}r}{d},\forall j$, where SVD of $M_{r}$ is $U^{*}\Sigma^{*}(V^{*})^{T}$. Intuitively incoherent matrices have mass spread over almost all entries whereas coherent matrices have mass concentrated on only few entries.

To generate matrices with varying incoherence parameter $\mu_{0}$, we use the power law matrices model . $M_{r}=DUV^{T}D$, where $U$ and $V$ are random $n\times r$ orthonormal matrices and $D$ is a diagonal matrix with $D_{ii}\propto\frac{1}{i^{\alpha}}$. For $\alpha=0$ $M_{r}$ is an incoherent matrix with $\mu_{0}=O(1)$ and for $\alpha=1$ $M_{r}$ is a coherent matrix with $\mu_{0}=O(n)$.

The input to algorithms is the matrix $M=M_{r}+Z$, where $Z$ is a Gaussian noise matrix with $||Z||=0.01,0.05$ and $0.1$. Correspondingly Frobenius norm of $Z$ is $||Z||*\sqrt{1000}/2$, which is $0.16,0.79$ and $1.6$ respectively.

Each simulation is averaged over 20 different runs. We run the WAltMin step of the algorithm for 15 iterations. Note that using different set of samples in each iteration of WAltMin subroutine 2 is generally observed to be not required in practice. Hence we use the same set of samples for all iterations.

In the first plot we compare the error $||M_{r}-\widehat{M}_{r}||$ of our algorithm LELA 1 with the random projection based algorithm . We use the matrix with each entry an independent Gaussian random variable as the sketching matrix, for the random projection algorithm. Other choices are Walsh-Hadamard based transform and sparse embedding matrices .

We compare the error of both algorithms as we vary number of samples $m$ for algorithm 1, equivalently varying the dimension of random projection $l=m/n$ for the random projection algorithm. In figure 1 we plot the error $||M_{r}-\widehat{M}_{r}||$ with varying number of samples $m$ for both the algorithms. For incoherent matrices we see that LELA algorithm has almost the same error as the random projection algorithm 1(a). But for coherent matrices we notice that in figure 1(b) LELA has significantly smaller error.

Now we consider the setting of computing low rank approximation of $YY^{T}$ given $Y$ using algorithm LELA direct discussed in section 3.3 with sampling (5). In figure 2 we compare this algorithm with a stagewise algorithm, which computes low rank approximation $\widehat{Y}_{r}$ from $Y$ first and sets the rank-$r$ approximation of $YY^{T}$ as $\widehat{Y}_{r}\widehat{Y}^{T}_{r}$. As discussed in section 3.3 direct approximation of $YY^{T}$ has less error than that of computing $\widehat{Y}_{r}\widehat{Y}^{T}_{r}$. Again plot 2(a) is for incoherent matrices and plot 2(b) is for coherent matrices.

Finally in figure 3 we consider the case where $A$ and $B$ are two rank $2r$ matrices with $AB$ being a rank $r$ matrix. Here the top $r$ dimensional row space of $A$ is orthogonal to the top $r$ dimensional column space of $B$. Hence simple algorithms that compute rank $r$ approximation of $A$ and $B$ first and then multiply will have high error as compared to that of LELA direct.

References

Appendix A Concentration Inequalities

In this section we will review couple of concentration inequalities we use in the proofs.

Let $X_{1},...X_{n}$ be independent scalar random variables. Let $|X_{i}|\leq L,\forall i~{}w.p.~{}1$. Then,

Recall the Shatten-$p$ norm of a matrix $X$ is

$\sigma_{i}(X)$ is the $i$th singular value of $X$. In particular for $p=2$, Shatten-$2$ norm is the Frobenius norm of the matrix.

[Matrix Chebyshev Inequality ] Let $X$ be a random matrix. For all $t>0$,

Appendix B Proofs of section 3

In this section we will present proof for Theorem 3.1. For simplicity we will only discuss proofs for the case when matrix is square. Rectangular case is a simple extension. We will provide proofs of the supporting lemmas first.

Let $\hat{q}_{ij}=\min\{q_{ij},1\}.$ This is to make sure the probabilities are all less than 1. Recall the definition of weights $w_{ij}=1/\hat{q}_{ij}$ when $\hat{q}_{ij}>0$ and else. Note that $\sum_{ij}\hat{q}_{ij}\leq m$. Also let $m\geq\beta nr\log(n)$.

Let $\{\delta_{ij}\}$ be the indicator random variables and $\delta_{ij}=1$ with probability $\hat{q}_{ij}$. Define $\Omega$ to be the sampling operator with $\Omega_{ij}=\delta_{ij}$. Define the weighted sampling operator $R_{\Omega}$ such that, $R_{\Omega}(M)_{ij}=\delta_{ij}w_{ij}M_{ij}$.

First we will abstract out the properties of the sampling distribution (2) that we use in the rest of the proof.

For $\Omega$ generated according to (2) and under the assumptions of Lemma 3.2 the following holds, for all $(i,j)$ such that $q_{ij}\leq 1$.

The proof of the lemma B.1 is straightforward from the definition of $q_{ij}$.

Now we will provide proof of the initialization lemma 3.2. Proof of lemma 3.2:

The proof of this lemma has two parts. 1) We show that

2) We show that the trimming step of algorithm 2 gives the required row norm bounds on $\widehat{U}^{(0)}$.

Proof of the first step: We prove the proof of the first part using the matrix Bernstein inequality. Note that the $L1$ term in the sampling distribution will help in getting good bounds on absolute magnitude of random variables $X_{ij}$ in this proof.

First we will bound $\|X_{ij}\|$. When $q_{ij}\geq 1$, $\hat{q}_{ij}=1$ and $\delta_{ij}=1$, and $X_{ij}=0$ with probability 1. Hence we only need to consider cases when $\hat{q}_{ij}=q_{ij}\leq 1$. We will assume this in all the proofs without explicitly mentioning it any more.

$\zeta_{1}$ follows from $\hat{q}_{ij}\leq 1$.

Hence, $\|X_{ij}\|$ is bounded by $L=\frac{2n}{m}\|M\|_{F}$. Now we will bound the variance.

Hence we get $\|M-\mathcal{P}_{r}(R_{\Omega}(M))\|\leq\|M-R_{\Omega}(M)\|+\|R_{\Omega}(M)-\mathcal{P}_{r}(R_{\Omega}(M))\|\leq\|M-M_{r}\|+2\delta\|M\|_{F}$, which implies $\|M_{r}-\mathcal{P}_{r}(R_{\Omega}(M))\|\leq 2\|M-M_{r}\|+2\delta\|M\|_{F}$.

Let SVD of $\mathcal{P}_{r}(R_{\Omega}(M))$ be $U^{(0)}\Sigma^{(0)}(V^{(0)})^{T}$. Hence,

This implies $dist(U^{(0)},U^{*})\leq\frac{2\|M-M_{r}\|+2\delta\|M\|_{F}}{\sigma^{*}_{r}}\leq\frac{1}{144r}.$ This follows from the assumption $\|M-M_{r}\|_{F}\leq\frac{1}{576\kappa r^{1.5}}\|M_{r}\|_{F}$ and $\delta\leq\frac{1}{576\kappa r^{1.5}}$. $\kappa=\frac{\sigma^{*}_{1}}{\sigma^{*}_{r}}$ is the the condition number of $M_{r}$.

Proof of the trimming step: From previous step we know that $\|R_{\Omega}(M)-M\|\leq\delta\left\|M\right\|_{F}$ and consequently $dist(U^{(0)},U^{*})\leq\delta_{2}$. Let,

be the estimates for the left leverages scores of the matrix $M$. Since $\|M-M_{r}\|_{F}\leq\|M_{r}\|_{F}$, $l_{i}^{2}\geq\frac{\sum_{k=1}^{r}(\sigma^{*}_{k})^{2}(U^{*}_{ik})^{2}}{\sum_{k=1}^{r}(\sigma^{*}_{k})^{2}}$.

since $\left|(U^{(0)}_{j})_{i}-(\bar{u}_{j})_{i}\right|\geq 2l_{i}-l_{i}=\sqrt{\frac{\sum_{k=1}^{r}(\sigma^{*}_{k})^{2}(U^{*}_{ik})^{2}}{\sum_{k=1}^{r}(\sigma^{*}_{k})^{2}}}$.

First we will show that this trimming step will not increase the distance to $U^{*}$ by much. To bound the $dist(\widehat{U}^{(0)},U^{*})$ consider, $\|(u^{*}_{\perp})^{T}U\|$, where $u^{*}_{\perp}$ is some vector perpendicular to $U^{*}$.

for $\delta_{2}\leq\frac{1}{144r}$. Second we will bound $\|(\widehat{U}^{(0)})^{i}\|$.

Hence we finish the proof of the second part of the lemma. ∎

B.2 Weighted AltMin analysis

We first provide proof of Lemma 3.3 for rank-1 case to explain the main ideas and in the next section we will discuss rank-$r$ case. Hence $M_{1}=\sigma^{*}u^{*}(v^{*})^{T}$. Before the proof of the lemma we will prove couple of supporting lemmas.

The weighted alternating minimization updates at the $t+1$ iteration are,

Writing in terms of power method updates we get,

where $B$ and $C$ are diagonal matrices with $B_{jj}=\sum_{i}\delta_{ij}w_{ij}(u^{t}_{i})^{2}$ and $C_{jj}=\sum_{i}\delta_{ij}w_{ij}u^{t}_{i}u^{*}_{i}$ and $y$ is the vector $R_{\Omega}(M-M_{1})^{T}u^{t}$ with entries $y_{j}=\sum_{i}\delta_{ij}w_{ij}u^{t}_{i}(M-M_{1})_{ij}$.

Now we will bound the error caused by the $M-M_{1}$ component in each iteration.

For $\Omega$ generated according to (2) and under the assumptions of Lemma 3.3 the following holds:

with probability greater that $1-\frac{\gamma}{T\log(n)}$, for $m\geq\beta n\log(n)$, $\beta\geq\frac{4c_{1}^{2}T}{\gamma\delta^{2}}$. Hence,$\left\|(u^{t})^{T}R_{\Omega}(M-M_{1})\right\|\leq dist(u^{t},u^{*})\|M-M_{1}\|+\delta\left\|M-M_{1}\right\|_{F},$ for constant $\delta$.

$\zeta_{1}$ follows from the fact that $X_{ij}$ are zero mean independent random variables. $\zeta_{2}$ follows from (20). Hence applying the matrix Chebyshev inequality for $p=2$ and $t=\delta\|M-M_{1}\|_{F}$ gives the result. ∎

For $\Omega$ sampled according to (2) and under the assumptions of Lemma 3.3 the following holds:

with probability greater that $1-\frac{2}{n^{2}}$, for $m\geq\beta n\log(n)$, $\beta\geq\frac{16}{\delta_{1}^{2}}$ and $\delta_{1}\leq 3$.

$\zeta_{1}$ follows from (11). Also it is easy to check that $|X_{j}|\leq\frac{16n}{m}$. Now applying Bernstein inequality gives the result. ∎

For $\Omega$ sampled according to (2) and under the assumptions of Lemma 3.3 the following holds:

with probability greater than $1-\frac{2}{n^{2}}$, for $m\geq\beta n\log(n),\beta\geq\frac{48c_{1}^{2}}{\delta_{1}^{2}}$ and $\delta_{1}\leq 3$.

Hence the $jth$ coordinate of the error term in equation (16) is

Recall that $\alpha_{i}=u^{t}_{i}(\langle u^{*},u^{t}\rangle u^{t}_{i}-u^{*}_{i}).$ Let $X_{ij}=\delta_{ij}w_{ij}\alpha_{i}v^{*}_{j}e_{j}e_{1}^{T}$, for $i,j$ in $[1,..n]$. Note that $X_{ij}$ are independent random matrices. Then $(\langle u^{t},u^{*}\rangle B-C)v^{*}$ is the first and the only column of the matrix $\sum_{i,j=1}^{n}X_{ij}$. We will bound $\|\sum_{ij=1}^{n}X_{ij}\|$ using matrix Bernstein inequality.

because $\sum_{i}\alpha_{i}=0$. Now we will give a bound on $\|X_{ij}\|$.

The lemma follows from applying matrix Bernstein inequality. ∎

Now we will provide proof of lemma 3.3. Proof of lemma 3.3:[Rank-1 case]

Now we will prove that the distance between $u^{t}$, $u^{*}$ and $v^{t+1},v^{*}$ decreases with each iteration. Recall that from the assumptions of the lemma we have the following row norm bounds for $u^{t}$;

First we will prove that $dist(u^{t},u^{*})$ decreases in each iteration and second we will prove that $v^{t+1}$ satisfies similar bound on its row norms.

Using Lemma B.3, Lemma B.4 and equation (16) we get,

Hence by applying the noise bounds Lemma B.2 we get,

$\zeta_{1}$ follows from $\delta_{1}\leq\frac{1}{2}$. $\zeta_{2}$ follows from using $\langle u^{t},u^{*}\rangle\geq\langle u^{0},u^{*}\rangle$. $\zeta_{3}$ follows from $(\langle u^{*},u^{0}\rangle-2\delta_{1}\sqrt{1-\langle u^{*},u^{0}\rangle^{2}}\geq\frac{1}{2}$, $\delta\leq\frac{1}{20}$ and $\delta_{1}\leq\frac{1}{20}$. Hence

From Lemma B.3 and (20) we get that $\left|\sum_{i}\delta_{ij}w_{ij}(u^{t}_{i})^{2}-1\right|\leq\delta_{1}$ and $\left|\sum_{i}\delta_{ij}w_{ij}u^{*}_{i}u^{t}_{i}-\langle u^{*},u^{t}\rangle\right|\leq\delta_{1}$, when $\beta\geq\frac{16c_{1}^{2}}{\delta_{1}^{2}}.$ Hence,

We will bound using $\sum_{i}\delta_{ij}w_{ij}u^{t}_{i}M_{ij}$ using Bernstein inequality.

$\sum_{i}\operatorname{Var}(X_{i})=\sum_{i}\hat{q}_{ij}(1-\hat{q}_{ij})(w_{ij})^{2}(u^{t}_{i})^{2}M_{ij}^{2}\leq\sum_{i}w_{ij}(u^{t}_{i})^{2}M_{ij}^{2}\leq\frac{4nc_{1}^{2}}{m}\|M_{j}\|^{2}$. Finally $|X_{ij}|\leq\left|w_{ij}u^{t}_{i}M_{ij}\right|\leq\frac{4nc_{1}}{m}|M_{ij}|/\sqrt{\frac{|M_{ij}|}{\|M\|_{F}}}\leq\frac{4nc_{1}}{m}\sqrt{|M_{ij}|\|M\|_{F}}$.

Hence applying Bernstein inequality with $t=\delta\sqrt{\|M_{j}\|^{2}+|M_{ij}|\|M\|_{F}}$ gives, $\sum_{i}\delta_{ij}w_{ij}u^{t}_{i}M_{ij}\leq(1+\delta_{1})\sqrt{\|M_{j}\|^{2}+|M_{ij}|\|M\|_{F}}$ with probability greater than $1-\frac{2}{n^{3}}$ when $m\geq\frac{24c_{1}^{2}}{\delta_{1}^{2}}n\log(n)$. For $\delta_{1}\leq\frac{1}{20}$, we get, $\|\widehat{u}^{t}\|\left|\widehat{v}_{j}^{t+1}\right|\leq\frac{21}{19}\sqrt{\|M_{j}\|^{2}+|M_{ij}|\|M\|_{F}}$.

Now we will bound $\|\widehat{v}^{t+1}\|$.

$\zeta_{1}$ follows from Lemma B.4 and equations (16) and (23). $\zeta_{2}$ follows from using $\langle u^{*},\widehat{u}^{0}\rangle\leq\langle u^{*},u^{t}\rangle$ and $\delta_{1}\leq\frac{1}{20}$. $\zeta_{3}$ follows from the argument: for $\delta_{1}\leq\frac{1}{16}$, $\langle\widehat{u}^{0},u^{*}\rangle-2\delta_{1}\sqrt{1-\langle u^{*},\widehat{u}^{0}\rangle^{2}}$ is greater than $\frac{1}{2}$, if $\langle\widehat{u}^{0},u^{*}\rangle\geq\frac{3}{5}$. This holds because $dist(u^{*},\widehat{u}^{0})\leq\frac{4}{5}$ from Lemma 3.2.

Hence we have shown that $v^{t+1}$ satisfies the row norm bounds. From Lemma B.2 we have, in each iteration with probability greater than $1-\frac{\gamma}{T\log(n)}$ we have $\left\|(u^{t})^{T}R_{\Omega}(M-M_{1})\right\|\leq dist(u^{t},u^{*})\|M-M_{1}\|+\delta\left\|M-M_{1}\right\|_{F}$. Hence the probability of failure in $T$ iterations is less than $\gamma.$ Lemma now follows from assumption on $m$. ∎

Now we have all the elements needed for proof of the Theorem 3.1. Proof of Theorem 3.1:[Rank-1 case]

Lemma 3.2 has shown that $\widehat{u}^{0}$ satisfies the row norm bounds condition. From Lemma 3.3 we get $dist(v^{t+1},v^{*})\leq\frac{1}{2}dist(u^{t},u^{*})+5\delta\|M-M_{1}\|_{F}/\sigma^{*}$. Hence $dist(v^{t+1},v^{*})\leq\frac{1}{4^{t}}dist(\widehat{u}^{0},u^{*})+10\delta\|M-M_{1}\|_{F}/\sigma^{*}$. After $t=O(\log(\frac{1}{\zeta}))$ iterations we get $dist(v^{t+1},v^{*})\leq\zeta+10\delta\|M-M_{1}\|_{F}/\sigma^{*}$ and $dist(u^{t},u^{*})\leq\zeta+10\delta\|M-M_{1}\|_{F}/\sigma^{*}$.