Low Rank Approximation and Regression in Input Sparsity Time

Introduction

A large body of work has been devoted to the study of fast randomized approximation algorithms for problems in numerical linear algebra. Several well-studied problems in this area include least squares regression, low rank approximation, and approximate computation of leverage scores. These problems have many applications in data mining , recommendation systems , information retrieval , web search , clustering , and learning mixtures of distributions . The use of randomization and approximation allows one to solve these problems much faster than with deterministic methods.

Another problem we consider is approximating the leverage scores. Given an $n\times d$ matrix $A$ with $n\gg d$, one can write $A=U\Sigma V^{\top}$ in its singular value decomposition, where the columns of $U$ are the left singular vectors, $\Sigma$ is a diagonal matrix, and the columns of $V$ are the right singular vectors. Although $U$ has orthonormal columns, not much can be immediately said about the squared lengths $\|U_{i}\|_{2}^{2}$ of its rows. These values are known as the leverage scores, and measure the extent to which the singular vectors of $A$ are correlated with the standard basis. The leverage scores are basis-independent, since they are equal to the diagonal elements of the projection matrix onto the span of the columns of $A$; see for background on leverage scores as well as a list of applications. The leverage scores will also play a crucial role in our work, as we shall see. The goal of approximating the leverage scores is to, simultaneously for each $i\in[n]$, output a constant factor approximation to $\|U_{i}\|_{2}^{2}$. Using randomization, this can be solved in $O(nd\log n+d^{3}\log d\log n)$ time .

There are also solutions for these problems based on sampling. They either get a weaker additive error , or they get bounded relative error but are slow . Many of the latter algorithms were improved independently by Deshpande and Vempala and Sarlós , and in followup work . There are also solutions based on iterative and conjugate-gradient methods, see, e.g., , or as recent examples. These methods repeatedly compute matrix-vector products $Ax$ for various vectors $x$; in the most common setting, such products require $\Theta(\operatorname{\mathtt{nnz}}(A))$ time. Thus the work per iteration of these methods is $\Theta(\operatorname{\mathtt{nnz}}(A))$, and the number of iterations $N$ that are performed depends on the desired accuracy, spectral properties of $A$, numerical stability issues, and other concerns, and can be large. A recent survey suggests that $N$ is typically $\Theta(k)$ for Krylov methods (such as Arnoldi and Lanczos iterations) to approximate the $k$ leading singular vectors . One can also use some of these techniques together, for example by first obtaining a preconditioner using the Johnson-Lindenstrauss (JL) transform, and then running an iterative method.

We resolve the above gaps by achieving algorithms for least squares regression, low rank approximation, and approximate leverage scores, whose time complexities have a leading order term that is $O(\operatorname{\mathtt{nnz}}(A))$, sometimes up to a log factor, with constant factors that are independent of any numerical properties of $A$. Our results are as follows:

see Theorem 38. We also note improved results for constrained regression, §7.6.

2 Techniques

In fact, our subspace embedding is nothing other than the CountSketch matrix in the data stream literature , see also . This matrix was also studied by Dasgupta, Kumar, and Sarlós . Formally, $S$ has a single randomly chosen non-zero entry $S_{h(j),j}$ in each column $j$, for a random mapping $h:[n]\mapsto[t]$. With probability $1/2$, $S_{h(j),j}=1$, and with probability $1/2$, $S_{h(j),j}=-1$.

We stress that our choice of matrices $S$ does not preserve the norms of an arbitrary set of $\exp(d)$ vectors with high probability, and so the above approach cannot work for our choice of matrices $S$. We instead critically use that these $\exp(d)$ vectors all come from a $d$-dimensional subspace (namely, $C(A)$), and therefore have a very special structure. The structural fact we use is that there is a fixed set $H$ of size $d/\alpha$ which depends only on the subspace, such that for any unit vector $y\in C(A)$, $H$ contains the indices of all coordinates of $y$ larger than $\sqrt{\alpha}$ in magnitude. The key property here is that the set $H$ is independent of $y$, or in other words, only a small set of coordinates could ever be large as we range over all unit vectors in the subspace. The set $H$ selects exactly the set of large leverage scores of the columns space $C(A)$!

Given this observation, by setting $t\geq K|H|^{2}$ for a large enough constant $K$, we have that with probability $1-1/K$, there are no two distinct $j\neq j^{\prime}$ with $j,j^{\prime}\in H$ for which $h(j)=h(j^{\prime})$. That is, we avoid the birthday paradox, and the coordinates in $H$ are “perfectly hashed” with large probability. Call this event $\mathcal{E}$, which we condition on.

Given a unit vector $y$ in the subspace, we can write it as $y^{H}+y^{L}$, where $y^{H}$ consists of $y$ with the coordinates in $[n]\setminus H$ replaced with , while $y^{L}$ consists of $y$ with the coordinates in $H$ replaced with . We seek to bound

Finally, we can bound $\langle Sy^{H},Sy^{L}\rangle$ as follows. Define $G\subseteq[n]\setminus H$ to be the set of coordinates $j$ for which $h(j)=h(j^{\prime})$ for a coordinate $j^{\prime}\in H$, that is, those coordinates in $[n]\setminus H$ which “collide” with an element of $H$. Then, $\langle Sy^{H},Sy^{L}\rangle=\langle Sy^{H},Sy^{L^{\prime}}\rangle$, where $y^{L^{\prime}}$ is a vector which agrees with $y^{L}$ on coordinates $j\in G$, and is on the remaining coordinates. By Cauchy-Schwarz, this is at most $\|Sy^{H}\|_{2}\cdot\|Sy^{L^{\prime}}\|_{2}$. We have already argued that $\|Sy^{H}\|_{2}=\|y^{H}\|_{2}\leq 1$ for unit vectors $y$. Moreover, we can again apply Theorem 2 of to bound $\|Sy^{L^{\prime}}\|_{2}$, since, conditioned on the coordinates of $y^{L^{\prime}}$ hashing to the set of items that the coordinates of $y^{H}$ hash to, they are otherwise random, and so we again have a mapping of our form (with a smaller $t$ and applied to a smaller $n$) applied to a vector with small infinity-norm. Therefore, $\|Sy^{L^{\prime}}\|_{2}\leq O(\varepsilon)+\|y^{L^{\prime}}\|_{2}$ with high probability. Finally, by Bernstein bounds, since the coordinates of $y^{L}$ are small and $t$ is sufficiently large, $\|y^{L^{\prime}}\|_{2}\leq\varepsilon$ with high probability. Hence, conditioned on event $\mathcal{E}$, $\|Sy\|_{2}=(1\pm\varepsilon)\|y\|_{2}$ with probability $1-\exp(-d)$, and we can complete the argument by union-bounding over a sufficiently fine net.

The first idea for bringing this down is that the analysis of can itself be tightened by using that we are applying it on vectors coming from a subspace instead of on a set of arbitrary vectors. This involves observing that in the analysis of , if on input vector $y$ and for every $i\in[t]$, $\sum_{j\mid h(j)=i}y_{j}^{2}$ is small then the remainder of the analysis of does not require that $\|y\|_{\infty}$ be small. Since our vectors come from a subspace, it suffices to show that for every $i\in[t]$, $\sum_{j\mid h(j)=i}\|U_{j}\|_{2}^{2}$ is small, where $\|U_{j}\|_{2}^{2}$ is the $j$-th leverage score of $A$. Therefore we do not need to perform this analysis for each $y$, but can condition on a single event, and this effectively allows us to increase $\alpha$ in the outline above, thereby reducing the size of $H$, and also the size of $t$ since we have $t=\Omega(|H|^{2})$. In fact, we instead follow a simpler and slightly tighter analysis of based on the Hanson-Wright inequality.

We also note that for applications such as least squares regression, it suffices to set $\varepsilon$ to be a constant in the subspace embedding, since we can use an approach in which, given constant-factor approximations to all of the leverage scores, can then achieve a $(1+\varepsilon)$-approximation to least squares regression by slightly over-sampling rows of the adjoined matrix $A\circ b$ proportional to its leverage scores, and solving the induced subproblem. This results in a better dependence on $\varepsilon$.

We can also compose our subspace embedding with a fast JL transform to further reduce $t$ to the optimal value of about $d/\varepsilon^{2}$. Since $S\cdot A$ already has small dimensions, applying a fast JL transform is now efficient.

Finally, we can use a recent result of to replace most dependencies on $d$ in our running times for regression with a dependence on the rank $r$ of $A$, which may be smaller.

3 Recent Related Work

here $\omega$ is the exponent for asymptotically fast matrix multiplication, and $\alpha>0$ is an arbitrary constant. (Some constant factors here are increasing in $\alpha$.)

Paul, Boutsidis, Magdon-Ismail, and Drineas implemented our subspace embeddings and found that in the TechTC-300 matrices, a collection of 300 sparse matrices of document-term data, with an average of 150 to 200 rows and 15,000 columns, our subspace embeddings as used for the projection step in their SVM classifier are about 20 times faster than the Fast JL Transform, while maintaining the same classification accuracy. Despite this large improvement in the time for projecting the data, further research is needed for SVM classification, as the JL Transform empirically possesses additional properties important for SVM which make it faster to classify the projected data, even though the time to project the data using our method is faster.

4 Outline

Sparse Embedding Matrices

We let ${\|A\|}_{F}$ or ${\|A\|}$ denote the Frobenius norm of matrix $A$, and ${\|A\|}_{2}$ denote the spectral norm of $A$.

$h:[n]\mapsto[t]$ is a random map so that for each $i\in[n]$, $h(i)=t^{\prime}$ for $t^{\prime}\in[t]$ with probability $1/t$.

$\Phi\in\{0,1\}^{t\times n}$ is a $t\times n$ binary matrix with $\Phi_{h(i),i}=1$, and all remaining entries .

$D$ is an $n\times n$ random diagonal matrix, with each diagonal entry independently chosen to be $+1$ or $-1$ with equal probability.

We will refer to a matrix of the form $\Phi D$ as a sparse embedding matrix.

Analysis

It will be convenient to regard the rows of $A$ and $U$ to be re-arranged so that the $u_{i}$ are in non-increasing order, so $u_{1}$ is largest; of course this order is unknown and un-used by our algorithms.

Let $T>0$ be a parameter. Throughout, we let $s\equiv\min\{i|u_{i}\leq T\}$, and $s^{\prime}\equiv\max\{i|\sum_{s\leq j\leq i}u_{j}\leq 1\}$.

We will use the notation $\llbracket P\rrbracket$, a function on event $P$, that returns 1 when $P$ holds, and 0 otherwise.

The following variation of Bernstein’s inequalitySee Wikipedia entry on Bernstein’s inequalities (probability theory). will be helpful.

For $L,T\geq 0$ and independent random variables $X_{i}\in[0,T]$ with $V\equiv\sum_{i}\operatorname{\mathbf{Var}}[X_{i}]$, if $V\leq LT^{2}/6$, then

Proof: Here Bernstein’s inequality says that for $Y_{i}\equiv X_{i}-\operatorname{\mathbf{E}}[X_{i}]$, so that $\operatorname{\mathbf{E}}[Y^{2}_{i}]=\operatorname{\mathbf{Var}}[X_{i}]=V$ and $|Y_{i}|\leq T$,

By the quadratic formula, the latter is no more than $-L$ when

which holds for $z\geq LT$ and $V\leq LT^{2}/6$.

We begin the analysis by considering $y_{s:n}$ for fixed unit vectors $y\in C(A)$. Since ${\|y\|}=1$, there must be a unit vector $x$ so that $y=Ux$, and so by Cauchy-Schwartz, ${\|y_{i}\|}^{2}\leq{\|U_{i,*}\|}^{2}{\|x\|}^{2}=u_{i}$. This implies that $\|y_{s:n}\|_{\infty}^{2}\leq u_{s}$. We extend this to all unit vectors in subsequent sections.

The following is similar to Lemma 6 of , and is a standard balls-and-bins analysis.

For $\delta_{h},T,t>0$, and $s\equiv\min\{i\mid u_{i}\leq T\}$, let $\mathcal{E}_{h}$ be the event that

where $W\equiv T\log(t/\delta_{h})+r/t$. If

then $\Pr[\mathcal{E}_{h}]\geq 1-\delta_{h}$.

Proof: We will apply Lemma 1 to prove that the bound holds for fixed $j\in[t]$ with failure probability $\delta_{h}/t$, and then apply a union bound.

Let $X_{i}$ denote the random variable $u_{i}\llbracket h(i)=j,i\geq s\rrbracket$. We have $0\leq X_{i}\leq T$, $\operatorname{\mathbf{E}}[X]=\sum_{i\geq s}u_{i}/t\leq r/t$, and $V=\sum_{i\geq s}\operatorname{\mathbf{E}}[X_{i}^{2}]=\sum_{i\geq s}u_{i}^{2}/t={\|u_{s:n}\|}^{2}/t$. Applying Lemma 1 with $L=\log(t/\delta_{h})$ gives

when ${\|u_{s:n}\|}^{2}/t\leq LT^{2}/6$, or $t\geq 6{\|u_{s:n}\|}^{2}/LT^{2}$.

Proof: We will use the following theorem, due to Hanson and Wright.

Our analysis uses some ideas from the proofs for Lemmas 7 and 8 of .

Since by assumption event $\mathcal{E}_{h}$ of Lemma 2 occurs, and for unit $y\in C(A)$, $y_{i^{\prime}}^{2}\leq u_{i^{\prime}}$ for all $i^{\prime}$, we have for $j\in[t]$ that $\sum_{i^{\prime}\in h^{-1}(j),i^{\prime}\geq s}y^{2}_{i^{\prime}}\leq W$. Hence

Putting this and (3.1) into the $Q$ of Theorem 4, we have,

2 Handling vectors with large entries

A small number of entries can be handled directly.

For given $s$, let $\mathcal{E_{B}}$ denote the event that $h(i)\neq h(i^{\prime})$ for all $i,i^{\prime}<s$. Then $\delta_{B}\equiv 1-\Pr[\mathcal{E}_{B}]\leq s^{2}/t$. Given event $\mathcal{E}_{B}$, we have that for any $y$,

Proof: Since $\Pr[h(i)=h(i^{\prime})]=1/t$, the probability that some such $i\neq i^{\prime}$ has $h(i)=h(i^{\prime})$ is at most $s^{2}/t$. The last claim follows by a union bound.

3 Handling all vectors

We have seen that $\Phi D$ preserves the norms for vectors with small entries (Lemma 3) and large entries (Lemma 5). Before proving a general bound, we need to prove a bound on the “cross terms”.

For $W$ as in Lemma 2, suppose the event $\mathcal{E}_{h}$ and $\mathcal{E}_{B}$ hold. Then for unit vector $y\in C(A)$, with failure probability at most $\delta_{C}$,

Proof: With the event $\mathcal{E}_{B}$, for each $i\geq s$ there is at most one $i^{\prime}<s$ with $h(i)=h(i^{\prime})$; let $z_{i}\equiv y_{i^{\prime}}D_{i^{\prime}i^{\prime}}$, and $z_{i}\equiv 0$ otherwise. We have for integer $p\geq 1$ using Khintchine’s inequality

where $C_{p}\leq\Gamma(p+1/2)^{1/p}=O(p)$, and the last inequality uses the assumption that $\mathcal{E}_{h}$ holds, and $\sum_{i^{\prime}<s}y_{i^{\prime}}^{2}\leq 1$. Putting $p=\log(1/\delta_{C})$ and applying the Markov inequality, we have

Therefore, with failure probability at most $\delta_{C}$, we have

Suppose the events $\mathcal{E}_{h}$ and $\mathcal{E}_{B}$ hold, and $W$ is as in Lemma 2. Then for $\delta_{y}>0$ there is an absolute constant $K_{y}$ such that, if $W\leq K_{y}\epsilon^{2}/\log(1/\delta_{y})$, then for unit vector $y\in C(A)$, with failure probability $\delta_{y}$, ${\|\Phi Dy\|}_{2}=(1\pm\varepsilon){\|y\|}_{2}$, when $\delta_{y}\leq 1/2$.

Proof: Assuming $\mathcal{E}_{h}$ and $\mathcal{E}_{B}$, we apply Lemmas 5, 3, and 6 , and have with failure probability at most $\delta_{L}+\delta_{C}$,

for the given $W$, putting $\delta_{L}=\delta_{C}=\delta_{y}/2$ and assuming $\delta_{y}\leq 1/2$. Thus $K_{y}\leq 1/9(K_{L}+K_{C})^{2}$ suffices.

$|E|\leq e^{cr}$ for $c=(\frac{1}{\gamma}+2)$.

For any $r\times r$ matrix $J$, if for every $u,v\in E$ we have $|u^{\top}Jv|\leq\varepsilon$, then for every unit vector $w$, we have $|w^{\top}Jw|\leq\frac{\varepsilon}{(1-\gamma)^{2}}$.

The following is our main theorem in this section.

There is $t=O((r/\epsilon)^{4}\log^{2}(r/\epsilon))$ such that with probability at least $9/10$, $\Phi D$ is a subspace embedding matrix for $A$; that is, for all $y\in C(A)$, ${\|\Phi Dy\|}_{2}=(1\pm\varepsilon){\|y\|}_{2}$. The embedding $\Phi D$ can be applied in $O(\operatorname{\mathtt{nnz}}(A))$ time. For $s=\min\{i^{\prime}\mid u_{i^{\prime}}\leq T\}$, where $T$ is a parameter in $\Omega(\epsilon^{2}/r\log(r/\epsilon))$, it suffices if $t\geq\max\{s^{2}/30,r/T\}$.

Proof: For suitable $t$, $T$, and $s$, with failure probability at most $\delta_{h}+\delta_{B}$, events $\mathcal{E}_{h}$ and $\mathcal{E}_{B}$ both hold. Conditioned on this, and assuming $W$ is sufficiently small as in Lemma 7, we have with failure probability $\delta_{y}$ for any fixed $y\in C(A)$ that $\|\Phi Dy\|_{2}=(1\pm\varepsilon)\|y\|_{2}$. Hence by Lemma 8, with failure probability $\delta_{h}+\delta_{B}+\delta_{y}K_{sub}^{r}$, $\|\Phi Dy\|_{2}=(1\pm 6\varepsilon)\|y\|_{2}$ for all $y\in C(A)$. We need $\delta_{h}+\delta_{B}+\delta_{y}K_{sub}^{r}\leq 1/10$, and the parameter conditions of Lemmas 2, Lemma 3, and Lemma 7 holding. Listing these conditions:

$\delta_{h}+\delta_{B}+\delta_{y}K_{sub}^{r}\leq 1/10$, where $\delta_{B}$ can be set to be $s^{2}/t$;

$t\geq 6{\|u_{s:n}\|}^{2}/\log(t/\delta_{h})T^{2}$;

$W=T\log(t/\delta_{h})+r/t\leq K_{y}\epsilon^{2}/\log(1/\delta_{y})$.

We put $\delta_{y}=K_{sub}^{-r}/30$, $\delta_{h}=1/30$, and require $t\geq s^{2}/30$. For the last condition it suffices that $T=O(\epsilon^{2}/r\log(t))$, and $t=\Omega(r^{2}/\epsilon^{2})$. The last condition implies the fourth condition for small enough constant $\varepsilon$. Also, since ${\|u_{s:n}\|}^{2}=\sum_{i\geq s}u_{i}^{2}\leq\sum_{i\geq s}u_{i}T\leq rT$, the bound for $T$ implies that $t=O((r/\epsilon)^{2}\log(t))$ suffices for Condition 3. Thus when the leverage scores are such that $s$ is small, $t$ can be $O((r/\epsilon)^{2}\log(r/\epsilon))$. Since $\sum_{i}u_{i}=r$, $s\leq r/T$ suffices, and so $t=O((r/T)^{2})=O((r/\epsilon)^{4}\log^{2}(r/\epsilon))$ suffices for the conditions of the theorem.

Partitioning Leverage Scores

Let $q\equiv\log_{2}1/T=O(\log(r/\varepsilon))$. We partition the leverage scores $u_{i}$ with $u_{i}\geq T$ into groups $G_{j}$, $j\in[q]$, where

Let $\beta_{j}\equiv 2^{-j}$, and $n_{j}\equiv|G_{j}|$. Since $\sum_{i=1}^{n}u_{i}=r$, we have for all $j$ that $n_{j}\leq r/\beta_{j}$.

We may also use $G_{j}$ to refer to the collection of rows of $U$ with leverage scores in $G_{j}$.

For given hash function $h$ and corresponding $\Phi$, let $G^{\prime}_{j}\subset G_{j}$ denote the collision indices of $G_{j}$, those $i\in G_{j}$ such that $h(i)=h(i^{\prime})$ for some $i^{\prime}\in G_{j}$. Let $k_{j}\equiv|G^{\prime}_{j}|$.

First, we bound the spectral norm of a submatrix of the orthogonal basis $U$ of $C(A)$, where the submatrix comprises rows of $G^{\prime}_{j}$.

We want to bound the spectral norm of the matrix $\hat{B}$ whose rows comprise those rows of $U$ in the collision set $G^{\prime}_{j}$. We let $t=\Theta(r^{2}q^{6}/\epsilon^{2})$ be the number of hash buckets. The expected number of collisions in the $t$ buckets is $\operatorname{\mathbf{E}}[|G^{\prime}_{j}|]=\frac{{n_{j}\choose 2}}{t}\leq\frac{n_{j}^{2}}{2t}.$ Let $\mathcal{D}_{j}$ be the event that the number $k_{j}\equiv|G^{\prime}_{j}|$ of such collisions in the $t$ buckets is at most $n_{j}^{2}q^{2}/t$. Let $\mathcal{D}=\cap_{j=1}^{q}\mathcal{D}_{j}$. By a Markov and a union bound, $\Pr[\mathcal{D}]\geq 1-1/(2q)$. We will assume that $\mathcal{D}$ occurs.

While each row in $B$ has some independent probability of participating in a collision, we first analyze a sampling scheme with replacement.

Our analysis will use a special case of the version of matrix Bernstein inequalities described by Recht.

Also $M\equiv{\|H_{m}\|}_{2}\leq 2\beta_{j}+\frac{1}{n_{j}}$.

Applying the above fact with these bounds for $\rho_{m}^{2}$ and $M$, we have

With probability $1-o(1)$, for all leverage score groups $G_{j}$, and for $U$ an orthonormal basis of $C(A)$, the submatrix $\hat{B}_{j}$ of $U$ consisting of rows in $G^{\prime}_{j}$, that is, those in $G_{j}$ that collide in a hash bucket with another row in $G_{j}$ under $\Phi$, has squared spectral norm $O(q(\beta_{j}+1/n_{j}+\sqrt{r/t}+n_{j}/t))$.

Proof: Fix a $j\in[q]$. If $n_{j}\equiv|G_{j}|\leq\sqrt{t}/q^{2}$, then with probability $1-o(1/q)$, the items in $G_{j}$ are perfectly hashed into the $t$ bins. So with probability $1-o(1)$, for all $j\in[q]$, if $n_{j}\leq\sqrt{t}/q^{2}$, then there are no collisions. Condition on this event.

Now consider a $j\in[q]$ for which $n_{j}\geq\sqrt{t}/q$. Then

When sampling with replacement, the expected number of distinct items is

2 Within-Group Errors

In this subsection, we show that for all $y\in L_{j}$, the error in estimating ${\|y\|}^{2}$ using $y^{\top}D\Phi^{\top}\Phi Dy$ is at most $O(\varepsilon)$.

For $y\in L_{j}$, the error in estimating ${\|y\|}^{2}$ by using $y^{\top}D\Phi^{\top}\Phi Dy$ contributed by collisions among coordinates $y_{i}$ for $i\in G_{j}$ is

and we need a bound on this quantity that holds with high probability.

By a standard balls-and-bins analysis, every bucket has $O(\log t)=O(q)$ collisions, with high probability, since $n_{j}\leq r/T\leq O(r^{2}/\varepsilon^{2})=O(t)$; we assume this event.

The squared Euclidean norm of the vector of all $y_{i}$ that appear in the summands, that is, with $i\in G^{\prime}_{j}$, is at most $\beta_{j}+1/n_{j}+\sqrt{r/t}+n_{j}/t)$ by Lemma 13. Thus the squared Euclidean norm of the vector comprising all summands in (2) is at most

By Khintchine’s inequality, for $p\geq 1$,

and therefore $|\kappa_{j}|^{2}$ is less than the last quantity, with failure probability at most $4^{-p}$.

Putting $p=k^{\prime}_{j}\equiv\min\{r,k_{j}\}$, with failure probability at most $4^{-k^{\prime}_{j}}$, for any fixed vector $y\in L_{j}$, the squared error in estimating ${\|y\|}^{2}$ using the sketch of $y$ is at most $O(k^{\prime}_{j}(q^{2}\beta_{j}(\beta_{j}+1/n_{j}+\sqrt{r/t}+n_{j}/t))$. Assuming the event $\cal D$ from the section above, we have $k^{\prime}_{j}\leq\min\{r,q\cdot n_{j}^{2}q/t\}$. We have, using $\beta_{j}n_{j}\leq r$,

and $r\cdot q^{2}\beta_{j}n_{j}/t\leq q^{2}r^{2}/t$, and finally

using $\beta_{j}n_{j}\leq r$. Putting these bounds on the terms together, the squared error is $O(q^{4}r^{2}/t)$, or $\epsilon^{2}/q^{2}$, for $t=\Omega(q^{6}r^{2}/\epsilon^{2})$, so that the error is $O(\epsilon/q)$.

Since the dimension of $L_{j}$ is bounded by $k^{\prime}_{j}$, it follows from the net argument of Lemma 8 that for all $y\in L_{j}$, ${\|Sy\|}^{2}={\|y\|}^{2}\pm O(\epsilon/q)$, and so the total error for unit $y\in C(A)$ is $O(\epsilon)$.

There is an absolute constant $C^{\prime}>0$ for which for any parameters $\delta_{1}\in(0,1)$, $P\geq 1$, and for sparse embedding dimension $t=O(P(r/\varepsilon)^{2}\log^{6}(r/\varepsilon))$, for all unit $y\in C(A)$, $\sum_{j\in[q]}{\|Sy^{j}\|}=1\pm C^{\prime}\epsilon/P\delta_{1}$, with failure probability at most $\delta_{1}+O(1/\log r)$, where $y^{j}$ denotes the member of $L_{j}$ derived from $y$.

3 Handling the Cross Terms

To complete the optimization, we must also handle the error due to “cross terms”.

Let $\delta_{1}\in(0,1)$ be an arbitrary parameter. For $j\neq j^{\prime}\in\{1,\ldots,q\}$, let the event $\mathcal{E}_{j,j^{\prime}}$ be that the number of bins containing both an item in $G_{j}$ and in $G_{j^{\prime}}$ is at most $\frac{n_{j}n_{j^{\prime}}q^{2}}{t\delta_{1}}.$ Let $\mathcal{E}=\cap_{j,j^{\prime}}\mathcal{E}_{j,j^{\prime}}$, the event that no pair of groups has too many inter-group collisions.

Proof: Fix a $j\neq j^{\prime}\in\{1,\ldots,q\}$. Then the expected number of bins containing an item in both $G_{j}$ and in $G_{j^{\prime}}$ is at most $t\cdot\frac{n_{j}}{t}\cdot\frac{n_{j^{\prime}}}{t}=\frac{n_{j}n_{j^{\prime}}}{t},$ and so by a Markov bound the number of bins containing an item in both $G_{j}$ and $G_{j^{\prime}}$ is at most $\frac{n_{j}n_{j^{\prime}}q^{2}}{t\delta_{1}}$ with probability at least $1-\delta_{1}/q^{2}$. The lemma follows by a union bound over the ${q\choose 2}$ choices of $j,j^{\prime}$.

In the remainder of the analysis, we set $t=P(r/\varepsilon)^{2}q^{6}$ for a parameter $P\geq 1$. Let $\mathcal{F}$ be the event that no bin contains more than $Cq$ elements of $\cup_{i=1}^{q}G_{j}$, where $C>0$ is an absolute constant.

Proof: Observe that $|\cup_{i=1}^{q}G_{j}|=\sum_{i=1}^{q}n_{j}\leq r\sum_{i=1}^{q}2^{j}\leq 2r^{2}/\varepsilon^{2}.$ By standard balls and bins analysis with the given $t$, with $P\geq 1$, with probability at least $1-1/r$ no bin contains more than $Cq$ elements, for a constant $C>0$.

Condition on events $\mathcal{E}$ and $\mathcal{F}$ occurring. Consider any unit vector $y=Ax$ in the column space of $A$. Consider any $j\neq j^{\prime}\in[q]$. Define the vector $y^{j}$: $y^{j}_{i}=y_{i}$ for $i\in G_{j}$, and $y^{j}_{i}=0$ otherwise. Then,

Proof: Since $\mathcal{E}$ occurs, the number of bins containing both an item in $G_{j}$ and $G_{j^{\prime}}$ is at most $n_{j}n_{j^{\prime}}q^{2}/(t\delta_{1})$. Call this set of bins $S$. Moreover, since $\mathcal{F}$ occurs, for each bin $i\in S$, there are at most $C\log r$ elements from $G_{j}$ in the bin and at most $C\log r$ elements from $G_{j^{\prime}}$ in the bin. Hence, for any $S=\Phi\cdot D$, we have, using $n_{j}\beta_{j}\leq r$ for all $j$,

The following is our main theorem concerning cross-terms in this section.

There is an absolute constant $C^{\prime}>0$ for which for any parameters $\delta_{1}\in(0,1)$, $P\geq 1$, and for sparse embedding dimension $t=O(P(r/\varepsilon)^{2}\log^{6}r)$, the event

occurs with failure probability at most $\delta_{1}+\frac{1}{r}$, where $y^{j},y^{j^{\prime}}$ are as defined in Lemma 17.

Proof: The theorem follows at once by combining Lemma 15, Lemma 16, and Lemma 17.

4 Putting it together

Putting the bounds for within-group and cross-term errors together, and replacing the use of Lemma 5 in the proof of Theorem 11, we have the following theorem.

There is an absolute constant $C^{\prime}>0$ for which for any parameters $\delta_{1}\in(0,1)$, $P\geq 1$, and for sparse embedding dimension $t=O(P(r/\varepsilon)^{2}\log^{6}(r/\varepsilon))$, for all unit $y\in C(A)$, ${\|Sy\|}=1\pm C^{\prime}\epsilon/P\delta_{1}$, with failure probability at most $\delta_{1}+O(1/\log r)$.

Generalized Sparse Embedding Matrices

Proof: We use Theorem 20 together with Lemma 8; for the latter, we need that for any fixed $y\in L$, $\|By\|_{2}^{2}=(1\pm\varepsilon/6)\|y\|_{2}^{2}$ with probability at least $1-\delta_{sub}$. By Theorem 20, we have this for $\delta_{sub}=(\delta\varepsilon/r)^{C_{KN}}$ for an arbitrarily large constant $C_{KN}>0$. Hence, by Lemma 8, there is a constant $K_{sub}>0$ so that with probability at least $1-(K_{sub})^{O(\log(r/\varepsilon\delta))}(\delta\varepsilon/r)^{C_{KN}}=1-(\delta\varepsilon/r)^{C_{subKN}}$, for all $y\in L$, $\|By\|_{2}^{2}=(1\pm\varepsilon)\|y\|_{2}^{2}$. Here we use that $C_{KN}>0$ can be made arbitrarily large, independent of $K_{sub}$.

2 The construction

Let $a=\Theta(\varepsilon^{-1}\log(r/\varepsilon\delta))$ and $v=\Theta(\varepsilon^{-1})$, be such that Theorem 20 and Corollary 21 apply with parameters $a$ and $v$, for a sufficiently large constant $C_{subKN}>0$. Further, let

where $C_{t}>0$ is a sufficiently large absolute constant, and let $t\equiv avq$.

Let $h:[n]\rightarrow[q]$ be a random hash function. For $i=1,2,\ldots,q$, define $a_{i}=|h^{-1}(i)|$. Note that $\sum_{i=1}^{q}a_{i}=n$.

We choose independent matrices $B^{(1)},\ldots,B^{(q)}$, with each $B^{(i)}$ as in Theorem 20 with parameters $a$ and $v$. Here $B^{(i)}$ is a $va\times a_{i}$ matrix. Finally, let $P$ be an $n\times n$ permutation matrix which, when applied to a matrix $A$, maps the rows of $A$ in the set $h^{-1}(1)$ to the set of rows $\{1,2,\ldots,a_{1}\}$, maps the rows of $A$ in the set $h^{-1}(2)$ to the set of rows $\{a_{1}+1,\ldots,a_{1}+a_{2}\}$, and for a general $i\in[q]$, maps the set of rows of $A$ in the set $h^{-1}(i)$ to the set of rows $\{a_{1}+a_{2}+\cdots+a_{i-1}+1,\ldots,a_{1}+a_{2}+\cdots+a_{i}\}$.

The map $S$ is defined to be the product of a block-diagonal matrix and the matrix $P$:

$S\cdot A$ can be computed in $O(\operatorname{\mathtt{nnz}}(A)(\log(r/\varepsilon\delta))/\varepsilon)$ time.

Proof: As $P$ is a permutation matrix, $P\cdot A$ can be computed in $O(\operatorname{\mathtt{nnz}}(A))$ time and has the same number of non-zero entries of $A$. For each non-zero entry of $P\cdot A$, we multiply it by $B^{(i)}$ for some $i$, which takes $O(a)=O(\log(r/\varepsilon\delta)/\varepsilon)$ time. Hence, the total time to compute $S\cdot A$ is $O(\operatorname{\mathtt{nnz}}(A)(\log(r/\varepsilon\delta))/\varepsilon)$.

3 Analysis

where $C_{T}$ is a sufficiently large absolute constant.

Let $s\equiv\min\{i^{\prime}\mid u_{i}\leq T\}$, and for $y^{\prime}\in C(A)$ of at most unit norm, let $y\equiv y^{\prime}_{s:n}$. Since $y_{i}^{2}\leq u_{i}$, this implies that $\|y\|_{\infty}^{2}\leq T$. Since $P$ is a permutation matrix, we have $\|Py\|_{\infty}^{2}\leq T$.

$N_{j}$ is a random sparse embedding matrix with $qv=t/a$ rows and $n$ columns.

Proof: $N_{j}$ has a single non-zero entry in each column, and the value of this non-zero entry is random in $\{+1,-1\}$. Hence, it remains to show that the distribution of locations of the non-zero entries of $N_{j}$ is the same as that in a sparse embedding matrix. This follows from the distribution of the values $a_{1},\ldots,a_{q}$, and the definition of $P$.

Let $\delta\in(0,1)$. For $j=1,\ldots,a$, let $\mathcal{E}_{h}^{j}$ be the event $\mathcal{E}_{h}$ of Lemma 2, applied to matrix $N_{j}$, with $\delta_{h}\equiv\delta/a$, and $W\equiv T\log(qv/\delta_{h})+r/qv\leq 2r/C_{T}q$. Suppose $\cap_{j\in[a]}\mathcal{E}_{h}^{j}$ holds. This event has probability at least $1-\delta$. Then there is an absolute constant $K_{L}$ such that with failure probability at most $\delta_{L}$,

Proof: We apply Lemma 3 with $N_{j}$ the sparse embedding matrix $\Phi D$, and $qv$, the number of rows of $N_{j}$, taking on the role of $t$ in Lemma 2, so that the parameter $W=T\log(qv/\delta_{h})+r/qv$ as in the lemma statement. (And since $t=avq$, $qv/\delta_{h}=t/\delta$, so $W=r/C_{T}q+r/qv\leq 2r/C_{T}q$.) Since ${\|u_{s:n}\|}^{2}\leq rT$, it suffices for Lemma 2 if $qv$ is at least $2rT/T^{2}\log(t/\delta_{h})=2C_{T}q$, or $v\geq 2C_{T}$.

With $\delta_{h}=\delta/a$, by a union bound $\cap_{j\in[a]}E_{h}^{j}$ occurs with failure probability $\delta$, as claimed.

We have, for given $N_{j}$, that with failure probability $\delta_{L}/a$, $|{\|N_{j}y_{s:n}\|}^{2}-{\|y_{s:n}\|}^{2}|\leq K_{L}\sqrt{W\log(a/\delta_{L})}$. Applying a union bound, and using

3.2 Vectors with large entries

Again, let $s\equiv\min\{i^{\prime}\mid u_{i^{\prime}}\leq T\}$. Since $\sum_{i}u_{i}=r$, we have

The following is a standard non-weighted balls-and-bins analysis.

Suppose the previously defined constant $C_{t}>0$ is sufficiently large. Let $\mathcal{E}_{nw}$ be the event that $|h^{-1}(i)\cap[s]|\leq C_{t}\log(r/\varepsilon\delta)$, for all $i\in[q]$. Then $\Pr[\mathcal{E}_{nw}]\geq 1-\delta/r$.

Hence, by a Chernoff bound, for a constant $C_{t}>0$,

The lemma now follows by a union bound over all $i\in[q]$.

Assume that $\mathcal{E}_{nw}$ holds. Let $\mathcal{E}_{s}$ be the event that for all $y\in C(A)$, $\|Sy_{1:(s-1)}\|^{2}=(1\pm\varepsilon/2)\|y_{1:(s-1)}\|^{2}$. Then $\Pr[\mathcal{E}_{s}]\geq 1-\delta/r$.

Proof: For $i=1,2,\ldots,q$, let $L^{i}$ be the at most $C_{t}\log(r/\varepsilon\delta)$-dimensional subspace which is the restriction of the column space $C(A)$ to coordinates $j$ with $h(j)=i$ and $j<s$. By Corollary 21, for any fixed $i$, with probability at least $1-(\delta\varepsilon/r)^{C_{subKN}}$, for all $y\in L^{i}$, $\|Sy\|^{2}=(1\pm\varepsilon)\|y\|^{2}$. By a union bound and sufficiently large $C_{subKN}>0$, this holds for all $i\in[q]$ with probability at least $1-q(\delta\varepsilon/r)^{C_{subKN}}>1-\delta/r$. This condition implies $\mathcal{E}_{s}$, since $y_{1:(s-1)}$ can be expressed as $\sum_{i\in[q]}y^{(i)}$, where each $y^{(i)}\in L^{i}$, and letting $\hat{B}^{(i)}$ denote the $va$ rows of $S$ corresponding to entries from $B^{(i)}$,