Frequent Directions : Simple and Deterministic Matrix Sketching

Introduction

The data streaming paradigm considers computation on a large data set $A$ where data items arrive in arbitrary order, are processed, and then never seen again. It also enforces that only a small amount of memory is available at any given time. This small space constraint is critical when the full data set cannot fit in memory or disk. Typically, the amount of space required is traded off with the accuracy of the computation on $A$. Usually the computation results in some summary $S(A)$ of $A$, and this trade-off determines how accurate one can be with the available space resources.

Modern large data sets are often viewed as large matrices. For example, textual data in the bag-of-words model is represented by a matrix whose rows correspond to documents. In large scale image analysis, each row in the matrix corresponds to one image and contains either pixel values or other derived feature values. Other large scale machine learning systems generate such matrices by converting each example into a list of numeric features. Low rank approximations for such matrices are used in common data mining tasks such as Principal Component Analysis (PCA), Latent Semantic Indexing (LSI), and k-means clustering. Regardless of the data source, the optimal low rank approximation for any matrix is obtained by its truncated Singular Value Decompositions (SVD).

In large matrices as above, one processor (and memory) is often incapable of handling all of the dataset $A$ in a feasible amount of time. Even reading a terabyte of data on a single processor can take many hours. Thus this computation is often spread among some set of machines. This renders standard SVD algorithms infeasible. Given a very large matrix $A$, a common approach is to compute in the streaming paradigm a sketch matrix $B$ that is significantly smaller than the original. A good sketch matrix $B$ is such that $A\approx B$ or $A^{T}A\approx B^{T}B$ and so computations can be performed on $B$ rather than on $A$ without much loss in precision.

In this paper, we propose a fourth approach: frequent directions. It is deterministic and draws on the similarity between the matrix sketching problem and the item frequency estimation problem. We provide additive and relative error bounds for it; we show how to merge summaries computed in parallel on disjoint subsets of data; we show it achieves the optimal tradeoff between space and accuracy, up to constant factors, for any row-update based summary; and we empirically demonstrate that it outperforms exemplars from all of the above described approaches.

2 Item Frequency Approximation

Our algorithm for sketching matrices is an extension of a well known algorithm for approximating item frequencies in streams. The following section shortly overviews the frequency approximation problem. Here, a stream $A=\{a_{1},\cdots,a_{n}\}$ has $n$ elements where each $a_{i}\in[d]$. Let $f_{j}=|\{a_{i}\in A\mid a_{i}=j\}|$ be the frequency of item $j$ and stands for number of times item $j$ appears in the stream. It is trivial to produce all item frequencies using $O(d)$ space simply by keeping a counter for each item. Although this method computes exact frequencies, it uses space linear to the size of domain which might be huge. Therefore, we are interested in using less space and producing approximate frequencies $\hat{f}_{j}$.

3 Connection to Matrix Sketching

4 Main results

Assuming a constant word size, any randomized matrix approximation streaming algorithm in the row-wise-updates model, which guarantees $\|A-\pi_{B}^{k}(A)\|^{2}_{F}\leq(1+\varepsilon)\|A-A_{k}\|^{2}_{F}$ and that succeeds with probability at least $2/3$, must use $\Omega(kd/\varepsilon)$ space.

Theorem 1.4 claims that FrequentDirections is optimal with respect to the tradeoff between sketch size and resulting accuracy. On the other hand, in terms of running time, FrequentDirections is not known to be optimal.

5 Practical Implications

Frequent Directions

This section proves our main results for Algorithm 1 which is our simplest and most space efficient algorithm. The reader will notice that we occasionally use inequalities instead of equalities at different parts of the proof to obtain three Properties. This is not unintentional. The reason is that we want the same exact proofs to hold also for Algorithm 1 which we describe in Section 3. Algorithm 1 is conceptually identical to Algorithm 1, it requires twice as much space but is far more efficient. Moreover, any algorithm which produces an approximate matrix $B$ which satisfies the following facts (for any choice of $\Delta$) will achieve the error bounds stated in Lemma 1.1 and Lemma 1.2.

In what follows, we denote by $\delta_{i}$, $B_{[i]}$, $C_{[i]}$ the values of $\delta$, $B$ and $C$ respectively after the $i$th row of $A$ was processed. Let $\Delta=\sum_{i=1}^{n}\delta_{i}$, be the total mass we subtract from the stream during the algorithm. To prove our result we first prove three auxiliary properties.

For any vector $x$ we have $\|Ax\|^{2}-\|Bx\|^{2}\geq 0$.

Use the observations that $\langle a_{i},x\rangle^{2}+\|B_{[i-1]}x\|^{2}=\|C_{[i]}x\|^{2}$.

To see this, first note that $\|C_{[i]}x\|^{2}-\|B_{[i]}x\|^{2}\leq\|C_{[i]}^{T}C_{[i]}-B_{[i]}^{T}B_{[i]}\|\leq\delta_{i}$. Now, consider the fact that $\|C_{[i]}x\|^{2}=\|B_{[i-1]}x\|^{2}+\|a_{i}x\|^{2}$. Substituting for $\|C_{[i]}x\|^{2}$ above and taking the sum yields

Combining this with $\sum_{i}\|C_{[i]}x\|^{2}-\|B_{[i]}x\|^{2}\leq\sum_{i}\delta_{i}=\Delta$ yields that $\|Ax\|^{2}-\|Bx\|^{2}\leq\Delta$. ∎

Now we can show that projecting $A$ onto $B_{k}$ provides a relative error approximation. Here, $y_{i}$ correspond to the singular vectors of $A$ as above and $v_{i}$ to the singular vectors of $B$ in a similar fashion.

Running Time Analysis

In extremely large datasets, the processing is often distributed among several machines. Each machine receives a disjoint input of raw data and is tasked with creating a small space summary. Then to get a global summary of the entire data, these summaries need to be combined. The core problem is illustrated in the case of just two machines, each process a data set $A_{1}$ and $A_{2}$, where $A=[A_{1};A_{2}]$, and create two summaries $B_{1}$ and $B_{2}$, respectively. Then the goal is to create a single summary $B$ which approximates $A$ using only $B_{1}$ and $B_{2}$. If $B$ can achieve the same formal space/error tradeoff as each $B_{1}$ to $A_{1}$ in a streaming algorithm, then the summary is called a mergeable summary .

First note that $B^{\prime}$ satisfies all facts with $\Delta^{\prime}=\Delta_{1}+\Delta_{2}$, by additivity of squared spectral norm along any direction $x$ (e.g. $\|B_{1}x\|^{2}+\|B_{2}x\|^{2}=\|B^{\prime}x\|^{2}$) and squared Frobenious norms (e.g. $\|B_{1}\|_{F}^{2}+\|B_{2}\|_{F}^{2}=\|B^{\prime}\|_{F}^{2}$), but has space twice as large as desired. Property 1 is straight forward for $B$ since $B$ only shrinks in all directions in relation to $B^{\prime}$. For Property 2 follows by considering any unit vector $x$ and expanding

This property trivially generalizes to any number of partitions of $A$. It is especially useful when the matrix (or data) is distributed across many machines. In this setting, each machine can independently compute a local sketch. These sketches can then be combined in an arbitrary order using FrequentDirections.

2 Worst case update time

Space Lower Bounds

In this section we show that FrequentDirections is space optimal with respect to the guaranteed accuracy. We present nearly-matching lower bounds for each case. We show the number of bits needed is equivalent to $d$ times the number of rows FrequentDirections requires. We first prove Theorem 1.3 showing the covariance error bound in Theorem 1.1 is nearly tight, regardless of streaming issues.

We start with a simple intuitive lemma showing an $\Omega(kd)$ lower bound, which we will refer to. We then prove our main $\Omega(kd/\epsilon)$ lower bound.

Any streaming algorithm which, for every input $A$, with constant probability (over its internal randomness) succeeds in outputting a matrix $R$ for which $\|A-AR^{\dagger}R\|_{F}\leq(1+\varepsilon)\|A-A_{k}\|_{F}$ must use $\Omega(kd)$ bits of space.

Let $\mathcal{S}$ be the set of $k$-dimensional subspaces over the vector space $GF(2)^{d}$, where $GF(2)$ denotes the finite field of $2$ elements with the usual modulo $2$ arithmetic. The cardinality of $\mathcal{S}$ is known to be

where the inequalities assume that $k\leq d/3$.

where $\alpha,\beta_{1},\ldots,\beta_{k}$ are integers. We can assume that the greatest common divisor (gcd) of $\alpha,\beta_{1},\ldots,\beta_{k}$ is $1$, otherwise the same conclusion holds after we divide $\alpha,\beta_{1},\ldots,\beta_{k}$ by the gcd. Note that (1) implies that $\alpha v=\sum_{i}\beta_{i}A^{2}_{i}\bmod 2$, i.e., when we take each of the coordinates modulo $2$. Since the $\beta_{i}$ cannot all be divisible by $2$ (since $\alpha$ would then be odd and so by the gcd condition the left hand side would contain a vector with at least one odd coordinate, contradicting that the right hand side is a vector with even coordinates), and the rows of $A^{2}$ form a basis over $GF(2^{d})$, the right hand side must be non-zero, which implies that $\alpha=1\bmod 2$. This implies that $v$ is in the span of the rows of $A^{2}$ over $GF(2^{d})$, a contradiction.

where $|R|$ denotes the expected length of the encoding of $R$, $H$ is the entropy function, and $I$ is the mutual information. For background on information theory, see . This completes the proof. ∎

2 Intuition for Main Lower Bound

The only other lower bounds for streaming algorithms for low rank approximation that we know of are due to Clarkson and Woodruff . As in their work, we use the Index problem in communication complexity to establish our bounds, which is a communication game between two players Alice and Bob, holding a string $x\in\{0,1\}^{r}$ and an index $i\in[r]\mathrel{\mathop{=}}:\{1,2,\ldots,r\}$, respectively. In this game Alice sends a single message to Bob who should output $x_{i}$ with constant probability. It is known (see, e.g., ) that this problem requires Alice’s message to be $\Omega(r)$ bits long. If Alg is a streaming algorithm for low rank approximation, and Alice can create a matrix $A_{x}$ while Bob can create a matrix $B_{i}$ (depending on their respective inputs $x$ and $i$), then if from the output of Alg on the concatenated matrix $[A_{x};B_{i}]$ Bob can output $x_{i}$ with constant probability, then the memory required of Alg is $\Omega(r)$ bits, since Alice’s message is the state of Alg after running it on $A_{x}$.

The main technical challenges are thus in showing how to choose $A_{x}$ and $B_{i}$, as well as showing how the output of Alg on $[A_{x};B_{i}]$ can be used to solve Index. This is where our work departs significantly from that of Clarkson and Woodruff . Indeed, a major challenge is that in Theorem 1.4, we only require the output to be the matrix $R$, whereas in Clarkson and Woodruff’s work from the output one can reconstruct $AR^{\dagger}R$. This causes technical complications, since there is much less information in the output of the algorithm to use to solve the communication game.

The intuition behind the proof of Theorem 1.4 is that given a $2\times d$ matrix $A=[1,x;1,0^{d}]$, where $x$ is a random unit vector, then if $P=R^{\dagger}R$ is a sufficiently good projection matrix for the low rank approximation problem on $A$, then the second row of $AP$ actually reveals a lot of information about $x$. This may be counterintuitive at first, since one may think that $[1,0^{d};1,0^{d}]$ is a perfectly good low rank approximation. However, it turns out that $[1,x/2;1,x/2]$ is a much better low rank approximation in Frobenius norm, and even this is not optimal. Therefore, Bob, who has $[1,0^{d}]$ together with the output $P$, can compute the second row of $AP$, which necessarily reveals a lot of information about $x$ (e.g., if $AP\approx[1,x/2;1,x/2]$, its second row would reveal a lot of information about $x$), and therefore one could hope to embed an instance of the Index problem into $x$. Most of the technical work is about reducing the general problem to this $2\times d$ primitive problem.

3 Proof of Main Lower Bound for Preserving Subspaces

Now let $c>0$ be a small constant to be determined. We consider the following two player problem between Alice and Bob: Alice has a $ck/\varepsilon\times d$ matrix $A$ which can be written as a block matrix $[I,R]$, where $I$ is the $ck/\varepsilon\times ck/\varepsilon$ identity matrix, and $R$ is a $ck/\varepsilon\times(d-ck/\varepsilon)$ matrix in which the entries are in $\{-1/(d-ck/\varepsilon)^{1/2},+1/(d-ck/\varepsilon)^{1/2}\}$. Here $[I,R]$ means we append the columns of $I$ to the left of the columns of $R$; see Figure 1. Bob is given a set of $k$ standard unit vectors $e_{i_{1}},\ldots,e_{i_{k}}$, for distinct $i_{1},\ldots,i_{k}\in[ck/\varepsilon]=\{1,2,\ldots,ck/\varepsilon\}$. Here we need $c/\varepsilon>1$, but we can assume $\varepsilon$ is less than a sufficiently small constant, as otherwise we would just need to prove an $\Omega(kd)$ lower bound, which is established by Lemma 4.1.

Denote this problem by $f$. We will show the randomized $1$-way communication complexity of this problem $R^{1-way}_{1/4}(f)$, in which Alice sends a single message to Bob and Bob fails with probability at most $1/4$, is $\Omega(kd/\varepsilon)$ bits. More precisely, let $\mu$ be the following product distribution on Alice and Bob’s inputs: the entries of $R$ are chosen independently and uniformly at random in $\{-1/(d-ck/\varepsilon)^{1/2},+1/(d-ck/\varepsilon)^{1/2}\}$, while $\{i_{1},\ldots,i_{k}\}$ is a uniformly random set among all sets of $k$ distinct indices in $[ck/\varepsilon]$. We will show that $D^{1-way}_{\mu,1/4}(f)=\Omega(kd/\varepsilon)$, where $D^{1-way}_{\mu,1/4}(f)$ denotes the minimum communication cost over all deterministic $1$-way (from Alice to Bob) protocols which fail with probability at most $1/4$ when the inputs are distributed according to $\mu$. By Yao’s minimax principle (see, e.g., ), $R^{1-way}_{1/4}(f)\geq D^{1-way}_{\mu,1/4}(f)$.

We use the following two-player problem Index in order to lower bound $D^{1-way}_{\mu,1/4}(f)$. In this problem Alice is given a string $x\in\{0,1\}^{r}$, while Bob is given an index $i\in[r]$. Alice sends a single message to Bob, who needs to output $x_{i}$ with probability at least $2/3$. Again by Yao’s minimax principle, we have that $R^{1-way}_{1/3}({\sf Index})\geq D^{1-way}_{\nu,1/3}({\sf Index})$, where $\nu$ is the distribution for which $x$ and $i$ are chosen independently and uniformly at random from their respective domains. The following is well-known.

$D^{1-way}_{\nu,1/3}({\sf Index})=\Omega(r).$

For $c$ a small enough positive constant, and $d\geq k/\varepsilon$, we have $D^{1-way}_{\mu,1/4}(f)=\Omega(dk/\varepsilon)$.

We proceed to lower bound $\|B-BP\|_{F}^{2}$ in a certain way, which will allow our reduction to Index to be carried out. We need the following fact:

((2.4) of ) Let $A$ be an $m\times n$ matrix with i.i.d. entries which are each $+1/\sqrt{n}$ with probability $1/2$ and $-1/\sqrt{n}$ with probability $1/2$, and suppose $m/n<1$. Then for all $t>0$,

where $\alpha,\alpha^{\prime}>0$ are absolute constants. Here $\|A\|_{2}$ is the operator norm $\sup_{x}\frac{\|Ax\|}{\|x\|}$ of $A$.

We apply Fact 4.2 to the matrix $R$, which implies,

and using that $d\geq k/\varepsilon$ and $c>0$ is a sufficiently small constant, this implies

where $\beta>0$ is an absolute constant (depending on $c$). Note that for $c>0$ sufficiently small, $(1+3\sqrt{c})^{2}\leq 1+7\sqrt{c}$. Let $\mathcal{E}$ be the event that $\|R\|_{2}^{2}\leq 1+7\sqrt{c}$, which we condition on.

We partition the rows of $B$ into $B_{1}$ and $B_{2}$, where $B_{1}$ contains those rows whose projection onto the first $ck/\varepsilon$ coordinates equals $e_{i}$ for some $i\notin\{i_{1},\ldots,i_{k}\}$. Note that $B_{1}$ is $(ck/\varepsilon-k)\times d$ and $B_{2}$ is $2k\times d$. Here, $B_{2}$ is $2k\times d$ since it includes the rows in $A$ indexed by $i_{1},\ldots,i_{k}$, together with the rows $e_{i_{1}},\ldots,e_{i_{k}}$. Let us also partition the rows of $R$ into $R_{T}$ and $R_{S}$, so that the union of the rows in $R_{T}$ and in $R_{S}$ is equal to $R$, where the rows of $R_{T}$ are the rows of $R$ in $B_{1}$, and the rows of $R_{S}$ are the non-zero rows of $R$ in $B_{2}$ (note that $k$ of the rows are non-zero and $k$ are zero in $B_{2}$ restricted to the columns in $R$).

For any unit vector $u$, write $u=u_{R}+u_{S}+u_{T}$, where $S=\{i_{1},\ldots,i_{k}\},T=[ck/\varepsilon]\setminus S$, and $R=[d]\setminus[ck/\varepsilon]$, and where $u_{A}$ for a set $A$ is on indices $j\notin A$. Then, conditioned on $\mathcal{E}$ occurring, $\|Bu\|^{2}\leq(1+7\sqrt{c})(2-\|u_{T}\|^{2}-\|u_{R}\|^{2}+2\|u_{S}+u_{T}\|\|u_{R}\|).$

Let $C$ be the matrix consisting of the top $ck/\varepsilon$ rows of $B$, so that $C$ has the form $[I,R]$, where $I$ is a $ck/\varepsilon\times ck/\varepsilon$ identity matrix. By construction of $B$, $\|Bu\|^{2}=\|u_{S}\|^{2}+\|Cu\|^{2}$. Now, $Cu=u_{S}+u_{T}+Ru_{R}$, and so

We will also make use of the following simple but tedious fact.

For $x\in$, the function $f(x)=2x\sqrt{1-x^{2}}-x^{2}$ is maximized when $x=\sqrt{1/2-\sqrt{5}/10}$. We define $\zeta$ to be the value of $f(x)$ at its maximum, where $\zeta=2/\sqrt{5}+\sqrt{5}/10-1/2\approx.618$.

Setting $y^{2}=1-x^{2}$, we can equivalently maximize $f(y)=-1+2y\sqrt{1-y^{2}}+y^{2}$, or equivalently $g(y)=2y\sqrt{1-y^{2}}+y^{2}$. Differentiating this expression and equating to , we have

Multiplying both sides by $\sqrt{1-y^{2}}$ one obtains the equation $4y^{2}-2=2y\sqrt{1-y^{2}}$, and squaring both sides, after some algebra one obtains $5y^{4}-5y^{2}+1=0$. Using the quadratic formula, we get that the maximizer satisfies $y^{2}=1/2+\sqrt{5}/10$, or $x^{2}=1/2-\sqrt{5}/10$. ∎

Conditioned on $\mathcal{E}$ occurring, $\|B\|_{2}^{2}\leq(1+7\sqrt{c})(2+\zeta).$

Suppose we replace the vector $u_{S}+u_{T}$ with an arbitrary vector supported on coordinates in $S$ with the same norm as $u_{S}+u_{T}$. Then the right hand side of this expression cannot increase, which means it is maximized when $\|u_{T}\|=0$, for which it equals $(1+7\sqrt{c})(2-\|u_{R}\|^{2}+2\sqrt{1-\|u_{R}\|^{2}}\|u_{R}\|),$ and setting $\|u_{R}\|$ to equal the $x$ in Fact 4.3, we see that this expression is at most $(1+7\sqrt{c})(2+\zeta)$. ∎

Write the projection matrix $P$ output by the streaming algorithm as $UU^{T}$, where $U$ is $d\times k$ with orthonormal columns $u^{i}$ (note that $R^{\dagger}R=P$). Applying Lemma 4.2 and Fact 4.3 to each of the columns $u^{i}$, we have:

Using the matrix Pythagorean theorem, we thus have,

Therefore, if Bob succeeds in solving $f$ on input $B$, then,

Comparing (4) and (5), we arrive at, conditioned on $\mathcal{E}$:

where $c_{1}>0$ is a constant that can be made arbitrarily small by making $c>0$ arbitrarily small.

Since $P$ is a projector, $\|BP\|_{F}=\|BU\|_{F}$. Write $U=\hat{U}+\bar{U}$, where the vectors in $\hat{U}$ are supported on $T$, and the vectors in $\bar{U}$ are supported on $[d]\setminus T$. We have,

where the first inequality uses $\|B\hat{U}\|_{F}\leq\|B\|_{2}\|\hat{U}\|_{F}$ and (6), the second inequality uses that event $\mathcal{E}$ occurs, and the third inequality holds for a constant $c_{2}>0$ that can be made arbitrarily small by making the constant $c>0$ arbitrarily small.

Combining with (5) and using the triangle inequality,

where $c_{3}>0$ is a constant that can be made arbitrarily small for $c>0$ an arbitrarily small constant (note that $c_{2}>0$ also becomes arbitrarily small as $c>0$ becomes arbitrarily small). Hence, $\|B\bar{U}\|_{F}^{2}\geq(2+\zeta)k-c_{3}k$, and together with Corollary 4.1, that implies $\|\bar{U}\|_{F}^{2}\geq k-c_{4}k$ for a constant $c_{4}$ that can be made arbitrarily small by making $c>0$ arbitrarily small.

Our next goal is to show that $\|B_{2}\bar{U}\|_{F}^{2}$ is almost as large as $\|B\bar{U}\|_{F}^{2}$. Consider any column $\bar{u}$ of $\bar{U}$, and write it as $\bar{u}_{S}+\bar{u}_{R}$. Hence,

Suppose $\|R_{S}\bar{u}_{R}\|=\tau\|\bar{u}_{R}\|$ for a value $0\leq\tau\leq 1+7\sqrt{c}$. Then

and hence, letting $\tau_{1},\ldots,\tau_{k}$ denote the corresponding values of $\tau$ for the $k$ columns of $\bar{U}$, we have

Comparing the square of (7) with (9), we have

where $c_{5}>0$ is a constant that can be made arbitrarily small by making $c>0$ an arbitrarily small constant. Now, $\|\bar{U}\|_{F}^{2}\geq k-c_{4}k$ as shown above, while since $\|R_{s}\bar{u}_{R}\|=\tau_{i}\|\bar{u}_{R}\|$ if $\bar{u}_{R}$ is the $i$-th column of $\bar{U}$, by (10) we have

for a constant $c_{6}$ that can be made arbitrarily small by making $c>0$ an arbitarily small constant.

Now $\|R\bar{U}_{R}\|_{F}^{2}\leq(1+7\sqrt{c})k$ since event $\mathcal{E}$ occurs, and $\|R\bar{U}_{R}\|_{F}^{2}=\|R_{T}\bar{U}_{R}\|_{F}^{2}+\|R_{S}\bar{U}_{R}\|_{F}^{2}$ since the rows of $R$ are the concatenation of rows of $R_{S}$ and $R_{T}$, so combining with (11), we arrive at

for a constant $c_{7}>0$ that can be made arbitrarily small by making $c>0$ arbitrarily small.

Combining the square of (7) with (12), we thus have

where the constant $c_{8}>0$ can be made arbitrarily small by making $c>0$ arbitrarily small.

or equivalently, $\|B_{2}-B_{2}P\|_{F}^{2}\leq 3k-((2+\zeta)k-c_{8}k)-(c_{2}k)+2k(((2+\zeta)-c_{8})c_{2})^{1/2}\leq(1-\zeta)k+c_{8}k+2k(((2+\zeta)-c_{8})c_{2})^{1/2}\leq(1-\zeta)k+c_{9}k$ for a constant $c_{9}>0$ that can be made arbitrarily small by making the constant $c>0$ small enough. This intuitively says that $P$ provides a good low rank approximation for the matrix $B_{2}$. Notice that by (15),

Conditioned on event $\mathcal{E}$, $\Delta\in[k-c_{10}k,k+c_{10}k],$ where $c_{10}>0$ is a constant that can be made arbitrarily small by making $c>0$ arbitrarily small.

The lemma now follows since $\sqrt{c}$ and $c_{9}$ can be made arbitrarily small by making the constant $c>0$ small enough. ∎

For the first part of the corollary, observe that

for a constant $c_{14}>0$ that can be made arbitrarily small by making $c>0$ arbitrarily small.

where $c_{15}>0$ is a constant that can be made arbitrarily small by making the constant $c>0$ arbitrarily small.

By a union bound over the occurrence of events $\mathcal{E},\mathcal{F}$, $\mathcal{G}$, and $\mathcal{H}$, and the streaming algorithm succeeding (which occurs with probability $3/4$), it follows that Bob succeeds in solving Index with probability at least $49/51-1/4-c_{16}>2/3$, as required. This completes the proof. ∎

Related Work on Matrix Sketching and Streaming

As mentioned in the introduction, there are a variety of techniques to sketch a matrix. There are also several ways to measure error and models to consider for streaming.

In this paper we focus on the row-update streaming model where each stream elements appends a row to the input matrix $A$. A more general model the entry-update model fixes the size of $A$ at $n\times d$, but each element indicates a single matrix entry and adds (or subtracts) to its value.

The accuracy of a sketch matrix $B$ can be measured in several ways. Most commonly one considers an $n\times d$, rank $k$ matrix $\hat{A}$ that is derived from $B$ (and sometimes $A$) and measures the projection error where proj-err $=\|A-\hat{A}\|_{F}^{2}/\|A-A_{k}\|_{F}^{2}$. When $\hat{A}$ can be derived entirely from $B$, we call this a construction result, and clearly requires at least $\Omega(n+d)$ space. When the space is required to be independent of $n$, then $\hat{A}$ either implicitly depends on $A$, or requires another pass of the data. It can then be defined one of two ways: $\hat{A}=\pi_{B}^{k}(A)$ (as is considered in this paper) takes $B_{k}$, the best rank-$k$ approximation to $B$, and then projects $A$ onto $B_{k}$; $\hat{A}=\Pi_{B}^{k}(A)$ projects $A$ onto $B$, and then takes the best rank-$k$ approximation of the result. Note that $\pi_{B}^{k}(A)$ is better than $\Pi_{B}^{k}(A)$, since it knows the rank-$k$ subspace to project onto without re-examining $A$.

We also consider covariance error where covar-err $=\|A^{T}A-B^{T}B\|_{2}/\|A\|_{F}^{2}$ in this paper. One can also bound $\|A^{T}A-B^{T}B\|_{2}/\|A-A_{k}\|_{F}^{2}$, but this has an extra parameter $k$, and is less clean. This measure captures the norm of $A$ along all directions (where as non-construction, projection error only indicates how accurate the choice of subspace is), but still does not require $\Omega(n)$ space.

Sketch paradigms

Given these models, there are several types of matrix sketches. We describe them here with a bit more specificity than in the Introduction, with particular attention to those that can operate in the row-update model we focus on. Specific exemplars are described which are used in out empirical study to follow. The first approach is to sparsify the matrix , by retaining a small number of non-zero. These algorithms typically assume to know the $n\times d$ dimensions of $A$, and are thus not directly applicable in out model.

Catalog of Related Bounds.

Recall that the error can be measured in several ways.

If it is not constructive, it is denoted by a P for projection if it uses a projection $\hat{A}=\pi_{B}^{k}(A)$. Alternatively if the weaker form of $\hat{A}=\Pi_{B}^{k}(A)$, then this is denoted $\textsf{P}_{r}$ where $r$ is the rank of $B$ before the projection.

In most recent results a projection error of $1+\varepsilon$ is obtained, and is denoted $\varepsilon$R.

Although in some cases a weaker additive error of the form $\|A-\hat{A}\|^{2}_{F}\leq\|A-A_{k}\|^{2}_{F}+\varepsilon\|A\|^{2}_{F}$ is achieved and denoted $\varepsilon$A. This can sometimes also be expressed as a spectral norm of the form $\|A-\hat{A}\|^{2}_{2}\leq\|A-A_{k}\|^{2}_{2}+\varepsilon\|A\|^{2}_{F}$ (note the error term $\varepsilon\|A\|^{2}_{F}$ still has a Frobenius norm). This is denoted $\varepsilon$L2.

In a few cases the error does not follow these patterns and we specially denote it.

Algorithms are randomized unless it is specified. In all tables we state bounds for a constant probability of failure. If we want to decrease the probability of failure to some parameter $\delta$, we can generally increase the size and runtime by $O(\log(1/\delta))$.

It is worth noting that under the construction model, and allowing streaming turnstile updates to each element of the matrix, the hashing algorithm has been shown space optimal by Clarkson and Woodruff (assuming each matrix entry requires $O(\log nd)$ bits, and otherwise off by only a $O(\log nd)$ factor). It is randomized and it constructs a decomposition of a rank $k$ matrix $\hat{A}$ that satisfies $\|A-\hat{A}\|_{F}\leq(1+\varepsilon)\|A-A_{k}\|_{F}$, with probability at least $1-\delta$. This provides a relative error construction bound of size $O((k/\varepsilon)(n+d)\log(nd))$ bits. They also show an $\Omega((k/\varepsilon)(n+d))$ bits lower bound. Our paper shows that the row-update model is strictly easier with a lower upper bound.

Although not explicitly described in their paper , one can directly use their techniques and analysis to achieve a weak form of a non-construction projection bound. One maintains a matrix $B=AS$ with $m=O((k/\varepsilon)\log(1/\delta))$ columns where $S$ is a $d\times m$ matrix where each entry is chosen from $\{-1,+1\}$ at random. Then setting $\hat{A}=\pi_{B}(A)$, achieves a ${\small\textsf{proj-err}}=1+\varepsilon$, however $B$ is rank $O((k/\varepsilon)\log(1/\delta))$ and hence that is the only bound on $\hat{A}$ as well.

Experimental Results

The other three baselines are Sampling, Hashing, and Random-projection, described in Section 5.

2 Datasets

We consider the real-world dataset Birds in which each row represents an image of a bird, and each column a feature. This dataset has 11788 data points and 312 features, and we recover rank $k=100$ for projection error experiments.

3 Results

Similar results in error and runtime are shown on the real Birds data set in Figure 8.

Acknowledgments

The authors thank Qin Zhang for discussion on hardness results. We also thank Stanley Eisenstat, and Mark Tygert for their guidance regarding efficient svd rank-1 updates.

References