Efficient Frequent Directions Algorithm for Sparse Matrices

Introduction

It is very common to represent data in the form of a matrix. For example, in text analysis under the bag-of-words model, a large corpus of documents can be represented as a matrix whose rows refer to the documents and columns correspond to words. A non-zero in the matrix corresponds to a word appearing in the a document. Similarly, in recommendation systems , preferences of users are represented as a matrix with rows corresponding to users and columns corresponding to items. Non-zero entires correspond to user ratings or actions.

A large set of data analytic tasks rely on obtaining a low-rank approximation of the data matrix. These include clustering, dimension reduction, principal component analysis (PCA), signal denoising, etc. Such approximations can be computed using the Singular Value Decompositions (SVD). For an $n\times d$ matrix $A$ ($d\leq n$) computing the SVD requires $O(nd^{2})$ time and $O(nd)$ space in memory on a single machine. In many scenarios, however, data matrices are extremely large and computing their SVD exactly is infeasible. Efficient approximate solutions exist for distributed setting or when data access otherwise is limited. In the row streaming model, the matrix rows are presented to the algorithm one by one in an arbitrary order. The algorithm is tasked with processing the stream in one pass while being severely restricted in its memory footprint. At the end of the stream, the algorithm must provide a sketch matrix $B$ which is a good approximation of $A$ even though it is significantly more compact. This is called matrix sketching.

Matrix sketching methods are designed to be parallelizable, space and time efficient, and easily updatable. Computing the sketch on each machine and then combining the sketches together should be as good as sketching the combined data from all the different machines. The streaming model is especially attractive since a sketch can be obtained and maintained as the data is being collected. Therefore, eliminating the need for data storage altogether.

Often matrices, as above, are sparse; most of their entries are zero. The work of argues that typical term-document matrices are sparse; documents contain no more than $5\%$ of all words. On wikipedia, most words appear on only a small constant number of pages. Similarly, in recommendation systems, in average a user rates or interacts with a small fraction of the available items: less than $6\%$ in some user-movies recommendation tasks and much fewer in physical purchases or online advertising. As such, most of these datasets are stored as sparse matrices.

There exist several techniques for producing low rank approximations of sparse matrices whose running time is $O(\operatorname{nnz}(A)\operatorname{poly}(k,1/\varepsilon))$ for some error parameter $\varepsilon\in(0,1)$. Here $\operatorname{nnz}(A)$ denotes the number of non-zeros in the matrix $A$. Examples include the power method , random projection techniques , projection-hashing , and instances of column selection techniques .

The notation $\pi_{B}(A)$ denotes the projection of the rows of $A$ on the span of the rows of $B$. In other words, $\pi_{B}(A)=AB^{\dagger}B$ where $(\cdot)^{\dagger}$ indicates taking the Moore-Penrose psuedoinverse. Alternatively, setting $[U,\Sigma,V]=\operatorname{SVD}(B)$, we have $\pi_{B}(A)=AVV^{T}$. We also denote $\pi_{B}^{k}(A)=AV_{k}V_{k}^{T}$, the right projection of $A$ on the top $k$ right singular vectors of $B$.

Matrix Sketching Prior Art

This section reviews only matrix sketching techniques that run in input sparsity time and whose output sketch is independent of the number of rows in the matrix. We categorize all known results into three main approaches (1) column/row subset selection (2) random projection based techniques and (3) iterative sketching techniques.

These techniques, which are also studied under the Column Subset Selection Problem (CSSP) in literature , form the sketch $B$ by selecting a subset of “important” columns of the input matrix $A$. They maintain the sparsity of $A$ and make the sketch $B$ to be more interpretable. These methods are not typically streaming, nor running in input sparsity time. The only method of this group which achieves both is by Drineas et al. that uses reservoir sampling to become streaming. They select $O(k/\varepsilon^{2})$ columns proportional to their squared norm and achieve the Frobenius norm error bound $\|A-\pi_{B_{k}}(A)\|_{F}^{2}\leq\|A-A_{k}\|_{F}^{2}+\varepsilon\|A\|_{F}^{2}$ with time complexity of $O((k^{2}/\varepsilon^{4})(d+k/\varepsilon^{2})+\operatorname{nnz}(A))$. In addition, they show that the spectral norm error bound $\|A-\pi_{B_{k}}(A)\|_{2}^{2}\leq\|A-A_{k}\|_{2}^{2}+\varepsilon\|A\|_{F}^{2}$ holds if one selects $O(1/\varepsilon^{2})$ columns. Rudelson et al. improved the latter error bound to $\|A-\pi_{B_{k}}(A)\|_{2}^{2}\leq\|A-A_{k}\|_{2}^{2}+\varepsilon\|A\|_{2}^{2}$ by selecting $O(r/\varepsilon^{4}\log{(r/\varepsilon^{4})})$ columns, where $r=\|A\|_{F}^{2}/\|A\|_{2}^{2}$ is the numeric rank of $A$. Note that in the result by , one would need $O(r^{2}/\varepsilon^{2})$ columns to obtain the same bound.

These techniques operate data-obliviously and maintain a $r\times d$ matrix $B=SA$ using a $r\times n$ random matrix $S$ which has the Johnson-Lindenstrauss Transform (JLT) property . Random projection methods work in the streaming model, are computationally efficient, and sufficiently accurate in practice . The state-of-the-art method of this approach is by Clarkson and Woodruff which was later improved slightly in . It uses a hashing matrix $S$ with only one non-zero entry in each column. Constructing this sketch takes only $O(\operatorname{nnz}(A)+n\cdot\operatorname{poly}(k/\varepsilon)+\operatorname{poly}(dk/\varepsilon))$ time, and guarantees that for any unit vector $x$ that $(1-\varepsilon)\|Ax\|\leq\|Bx\|\leq(1+\varepsilon)\|Ax\|.$ For these sparsity-efficient sketches using $r=O(d^{2}/\varepsilon^{2})$ also guarantees that $\|A-\pi_{B}(A)\|_{F}\leq(1+\varepsilon)\|A-A_{k}\|_{F}$.

1 Main Results

The expected running time of the algorithm is

Preliminaries

In this section we review some important properties about FrequentDirections and SimultaneousIteration which will be necessary for understanding and proving bounds on SparseFrequentDirections.

The FrequentDirections algorithm was introduced by Liberty and received an improved analysis by Ghashami et al.. The algorithm operates by collecting several rows of the input matrix and letting the sketch grow. Once the sketch doubles in size, a lossy DenseShrink operation reduces its size by a half. This process repeats throughout the stream. The running time of FrequentDirections and its error analysis are strongly coupled with the properties of the SVD used to perform the DenseShrink step.

An analysis of slightly generalized the one in . Let $B$ be the sketch resulting in applying FrequentDirections with a potentially different shrink operation to $A$. Then, the FrequentDirections asymptotic guarantees hold as long as the shrink operation exhibits three properties, for any positive $\Delta$ and a constant $\alpha\in(0,1)$.

For completeness, the exact guarantee is stated in Lemma 3.1.

where $\pi^{k}_{B}(\cdot)$ represents the projection operator onto $B_{k}$, the top $k$ singular vectors of $B$.

2 Simultaneous Iteration

Efficiently computing the singular vectors of matrices is one of the most well studies problems in scientific computing. Recent results give very strong approximation guarantees for block power method techniques . Several variants of this algorithm were studied under different names in the literature e.g. Simultaneous Iteration, Subspace Iteration, or Orthogonal Iteration . In this paper, we refer to this group of algorithms collectively as SimultaneousIteration. A generic version of SimultaneousIteration for rectangular matrices is described in Algorithm 1.

While this algorithm was already analyzed by , the proofs of manage to prove stable results that hold for any matrix independent of spectral gap issues. Unfortunately, an in depth discussion of these algorithms and their proof techniques is beyond the scope of this paper.

For the proof of correctness of SparseFrequentDirections, the main lemma proven by suffices. SimultaneousIteration (Algorithm 1) guarantees the three following error bounds with high probability:

Frobenius norm error bound: $\|A-ZZ^{T}A\|_{F}\leq(1+\varepsilon)\|A-A_{k}\|_{F}$

Spectral norm error bound: $\|A-ZZ^{T}A\|_{2}\leq(1+\varepsilon)\|A-A_{k}\|_{2}$

Per vector error bound: $|u_{i}^{T}AA^{T}u_{i}-z_{i}^{T}AA^{T}z_{i}|\leq\varepsilon\sigma_{k+1}^{2}$

for all $i$. Here $u_{i}$ denotes the $i$th left singular vector of $A$, and $\sigma_{k+1}$ is the ($k+1$)th singular value of $A$, and $z_{i}$ is the $i$th column of the matrix $Z$ returned by SimultaneousIteration.

In this paper, we show that SparseFrequentDirections can replace the computation of an exact SVD by using the results of with $\varepsilon$ being a constant. This alteration does give up the optimal asymptotic accuracy (matching that of FrequentDirections).

Sparse Frequent Directions

Our main result is stated in the next theorem. It follows from combining the proofs contained in the subsections below.

The VerifySpectral algorithm returns True if $\|C\|_{2}\leq 1$. If $\|C\|_{2}\geq 2$ it returns False with probability at least $1-\delta_{i}$.

If $\|C\|\leq 1$ than $\|C^{c\cdot\log(d/\delta_{i})}x\|\leq\|C\|^{c\cdot\log(d/\delta_{i})}\|x\|\leq 1$. If $\|C\|\geq 2$, consider execution $i$ of the method. Let $v_{1}$ denote the top singular vector of $C$. Then $\|C^{c\cdot\log(d/\delta_{i})}x\|\geq|\langle v_{1},x\rangle|2^{c\cdot\log(d/\delta_{i})}\geq 1$, for some constant $c$ as long as $|\langle v_{1},x\rangle|=\Omega(\operatorname{poly}(\delta_{i}/d))$. Let $\Phi(t^{\prime})$ denote the density function of the random variable $t^{\prime}=\langle v_{1},x\rangle$. Then $\Pr[|\langle v_{1},x\rangle|\leq t]=\int_{-t}^{t}\Phi(t^{\prime})dt^{\prime}\leq 2t\Phi(0)=O(t\sqrt{d})$. Setting the failure probability to be at most $\delta_{i}$, we conclude that $|\langle v_{1},x\rangle|=\Omega(\delta_{i}/\sqrt{d})$ with probability at least $1-\delta_{i}$. ∎

Therefore, VerifySpectral fails with probability at most $\delta_{i}$ during execution $i$. If any of VerifySpectral runs fail, BoostedSparseShrink and hence SparseFrequentDirections potentially fail. Taking the union bound over all invocations of VerifySpectral we obtain that SparseFrequentDirections fails with probability at most $\sum\delta_{i}\leq\sum_{i=1}^{\infty}\delta/2i^{2}\leq\delta$, hence it succeeds with probability at least $1-\delta$.

2 Space Usage and Runtime Analysis

Throughout this manuscript we assume the constant-word-size model. Integers and floating point numbers are represented by a constant number of bits. Random access into memory is assumed to require $O(1)$ time. In this model, multiplying a sparse matrix $A^{\prime}$ by a dense vector requires $O(\operatorname{nnz}(A^{\prime}))$ operations and storing $A^{\prime}$ requires $O(\operatorname{nnz}(A^{\prime}))$ bits of memory.

The running time of BoostedSparseShrink is dominated by those of SparseShrink and VerifySpectral, and its expected number of iterations. Note that, in expectation, they are each executed on any buffer matrix $A^{\prime}_{i}$ a small constant number of times because VerifySpectral succeeds with probability (much) greater than $1/2$. For asymptotic analysis it is identical to assuming they are each executed once.

Combining the above contributions to the total running time of the algorithm we obtain Fact 4.2.

Algorithm SparseFrequentDirections runs in expected time of

3 Error Analysis

therefore $\|A^{\prime}x\|^{2}-\|B^{\prime}x\|^{2}=(\|A^{\prime}x\|^{2}-\|Px\|^{2})+(\|Px\|^{2}-\|B^{\prime}x\|^{2})\geq 0$. ∎

We bound each term individually. The first term is bounded as

3.2 Error Analysis: BoostedSparseShrink and SparseFrequentDirections

We now consider the BoostedSparseShrink algorithm, and the looser version of Property 2 (the original version) as

Experiments

In this section we empirically validate that SparseFrequentDirections matches (and often improves upon) the accuracy of FrequentDirections, while running significantly faster on sparse real and synthetic datasets.

We do not implement SparseFrequentDirections exactly as described above. Instead we directly call SparseShrink in Algorithm 2 in place of BoostedSparseShrink. The randomized error analysis of SimultaneousIteration indicates that we may occasionally miss a subspace within a call of SimultaneousIteration and hence SparseShrink; but in practice this is not a catastrophic event, and as we will observe, does not prevent SparseFrequentDirections from obtaining small empirical error.

We ran all the algorithms under a common implementation framework to test their relative performance as accurately as possible. We ran the experiments on an Intel(R) Core(TM) 2.60 GHz CPU with 64GB of RAM running Ubuntu 14.04.3. All algorithms were coded in C, and compiled using gcc 4.8.4. All linear algebra operation on dense matrices (such as SVD) invoked those implemented in LAPACK.

We compare the performance of the two algorithms on both synthetic and real datasets. Each dataset is an $n\times d$ matrix $A$ containing $n$ datapoints in $d$ dimensions.

The real dataset is part of the $20$ Newsgroups dataset , that is a collection of approximately $20{,}000$ documents, partitioned across $20$ different newsgroups. However we use the ‘by date’ version of the data, where features (columns) are tokens and rows correspond to documents. This data matrix is a zero-one matrix with $11{,}314$ rows and $117{,}759$ columns. In our experiment, we use the transpose of the data and picked the first $d=3000$ columns, hence the subset matrix has $n=117{,}759$ rows and $d=3000$ columns; roughly $0.15\%$ of the subset matrix is non-zeros.

The synthetic data generates $n$ rows i.i.d. Each row receives exactly $z\ll d$ non-zeros (with default $z=100$ and $d=1000$), with the remaining entries as . The non-zeros are chosen as either $1$ or $-1$ at random. Each non-zero location is chosen without duplicates among the columns. The first $1.5z$ columns (e.g., 150), the “head”, have a higher probability of receiving a non-zero than the last $d-1.5z$ columns, the “tail”. The process to place a non-zero first chooses the head with probability $0.9$ or the tail with probability $0.1$. For whichever set of columns it chooses (head or tail), it places the non-zero uniformly at random among those columns.

Projection Error: proj-err $=\|A-\pi_{B_{k}}(A)\|_{F}^{2}/\|A-A_{k}\|_{F}^{2}$,

Covariance Error: cov-err $=\|A^{T}A-B^{T}B\|_{2}/\|A\|_{F}^{2}$,

1 Observations

By considering Table 2 on synthetic data and Figure 1 on the real data, we can vary and learn many aspects of the runtime and accuracy of SparseFrequentDirections and FrequentDirections.

Consider the last row of Table 2, the “Run Time” row, and the last column of Figure 1. SparseFrequentDirections is clearly faster than FrequentDirections for all datasets, except when the synthetic data becomes dense in the last column of the “Run Time” row, where $d=1000$ and $\operatorname{nnz}\text{per row}=500$ in the right-most data point. For the default values the improvement is between about a factor of $1.5$x and $2x$, but when the matrix is very sparse the improvement is $10$x or more. Very sparse synthetic examples are seen in the left data points of the last column, and in the right data points of the second column, of the “Run Time” row.

In particular, these two plots (the second and fourth columns of the “Run Time” row) really demonstrate the dependence of SparseFrequentDirections on $\operatorname{nnz}(A)$ and of FrequentDirections on $n\cdot d$. In the last column, we fix the matrix size $n$ and $d$, but increase the number of non-zeros $\operatorname{nnz}(A)$; the runtime of FrequentDirections is basically constant, while for SparseFrequentDirections it grows linearly. In the second column, we fix $n$ and $\operatorname{nnz}(A)$, but increase the number of columns $d$; the runtime of FrequentDirections grows linearly while the runtime for SparseFrequentDirections is basically constant.

These algorithms are designed for datasets with extremely large values of $n$; yet we only run on datasets with $n$ up to $60{,}000$ in Table 2, and $117{,}759$ in Figure 1. However, both FrequentDirections and SparseFrequentDirections have runtime that grows linearly with respect to the number of rows (assuming the sparsity is at an expected fixed rate per row for SparseFrequentDirections). This can also be seen empirically in the first column of the “Run Time” row where, after a small start-up cost, both FrequentDirections and SparseFrequentDirections grow linearly as a function of the number of data points $n$. Hence, it is valid to directly extrapolate these results for datasets of increased $n$.

We will next discuss the accuracy, as measured in Projection Error in the top row of Table 2 and left plot of Figure 1, and in Covariance Error in the middle row of Table 2 and middle plot of Figure 1. We observe that both FrequentDirections and SparseFrequentDirections obtain very small error (much smaller than upper bounded by the theory), as has been observed elsewhere . Moreover, the error for SparseFrequentDirections always nearly matches, or improves over FrequentDirections. We can likely attribute this improvement to being able to process more rows in each batch, and hence needing to perform the shrinking operation fewer overall times. The one small exception to SparseFrequentDirections having less Covariance Error than FrequentDirections is for extreme sparse datasets in the leftmost data points of Table 2, last column – we attribute this to some peculiar orthogonality of columns with near equal norms due to extreme sparsity.

References