Sketching as a Tool for Numerical Linear Algebra

Introduction

To give the reader a flavor of results in this survey, let us first consider the classical linear regression problem. In a special case of this problem one attempts to “fit” a line through a set of given points as best as possible.

For example, the familiar Ohm’s law states that the voltage $V$ is equal to the resistance $R$ times the electrical current $I$, or $V=R\cdot I$. Suppose one is given a set of $n$ example volate-current pairs $(v_{j},i_{j})$ but does not know the underlying resistance. In this case one is attempting to find the unknown slope of a line through the origin which best fits these examples, where best fits can take on a variety of different meanings.

More formally, in the standard setting there is one measured variable $b$, in the above example this would be the voltage, and a set of $d$ predictor variables $a_{1},\ldots,a_{d}$. In the above example $d=1$ and the single predictor variable is the electrical current. Further, it is assumed that the variables are linearly related up to a noise variable, that is $b=x_{0}+a_{1}x_{1}+\cdots+a_{d}x_{d}+\gamma$, where $x_{0},x_{1},\ldots,x_{d}$ are the coefficients of a hyperplane we are trying to learn (which does not go through the origin if $x_{0}\neq 0$), and $\gamma$ is a random variable which may be adversarially chosen, or may come from a distribution which we may have limited or no information about. The $x_{i}$ are also known as the model parameters. By introducing an additional predictor variable $a_{0}$ which is fixed to $1$, we can in fact assume that the unknown hyperplane goes through the origin, that is, it is an unknown subspace of codimension $1$. We will thus assume that $b=a_{1}x_{1}+\cdots+a_{d}x_{d}+\gamma$ and ignore the affine component throughout.

Subspace Embeddings and Least Squares Regression

We are interested in fast solutions to this problem, which we present in §2.5.

Definition 2 will be used in applications throughout this book, and sometimes for convenience we will drop the word oblivious.

Returning to Definition 2, the first usage of this in the numerical linear algebra community, to the best of our knowledge, was done by Sárlos, who proposed using Fast Johnson Lindenstrauss transforms to provide subspace embeddings. We follow the exposition in Sarlós for this .

which implies all inner products are preserved up to $\varepsilon$ by rescaling $\varepsilon$ by a constant.

There are many constructions of Johnson-Lindenstrauss transforms, possibly the simplest is given by the following theorem.

We will see a proof of Theorem 4 in Lemma 18.

By an argument of , it suffices to choose $\mathcal{N}$ so that for all ${\mathbf{y}}\in\mathcal{S}$, there exists a vector ${\mathbf{w}}\in\mathcal{N}$ for which $\|{\mathbf{y}}-{\mathbf{w}}\|_{2}\leq 1/2$. We will refer to $\mathcal{N}$ as a $(1/2)$-net for $\mathcal{S}$.

To see that $\mathcal{N}$ suffices, if ${\mathbf{y}}$ is a unit vector, then we can write

where $\|{\mathbf{y}}^{i}\|\leq{1\over 2^{i}}$ and ${\mathbf{y}}^{i}$ is a scalar multiple of a vector in $\mathcal{N}$. This is because we can write ${\mathbf{y}}={\mathbf{y}}^{0}+({\mathbf{y}}-{\mathbf{y}}^{0})$ where ${\mathbf{y}}^{0}\in\mathcal{N}$ and $\|{\mathbf{y}}-{\mathbf{y}}^{0}\|\leq 1/2$ by the definition of $\mathcal{N}$. Then, ${\mathbf{y}}-{\mathbf{y}}^{0}={\mathbf{y}}^{1}+(({\mathbf{y}}-{\mathbf{y}}^{0})-{\mathbf{y}}^{1})$ where ${\mathbf{y}}^{1}\in\mathcal{N}$ and

The expansion in (3) then follows by induction. But then,

where the first equality follows by (3), the second equality follows by expanding the square, the third equality follows from (2), and the fourth equality is what we want (after rescaling $\varepsilon$ by a constant factor).

We show the existence of a small $(1/2)$-net $\mathcal{N}$ via a standard argument.

For any $0<\gamma<1$, there exists a $\gamma$-net $\mathcal{N}$ of $\mathcal{S}$ for which $|\mathcal{N}|\leq(1+4/\gamma)^{d}$.

This can be done by choosing a maximal set $\mathcal{N}^{\prime}$ of points on $\mathcal{S}^{t-1}$ so that no two points are within distance $\gamma/2$ from each other. It follows that the balls of radius $\gamma/4$ centered at these points are disjoint, but on the other hand they are all contained in the ball of radius $1+\gamma/4$ centered at the origin. The volume of the latter ball is a factor $(1+\gamma/4)^{t}/(\gamma/4)^{t}$ larger than the smaller balls, which implies $|\mathcal{N}^{\prime}|\leq(1+4/\gamma)^{t}$. See, e.g., for more details.

It follows by setting $V=\mathcal{N}$ and $f=9^{d}$ in Theorem 4, we can then apply Lemma 5 and (2) to obtain the following theorem. Note that the net size does not depend on $\varepsilon$, since we just need a $1/2$-net for the argument, even though the theorem holds for general $\varepsilon$.

The Fast Johnson Lindenstrauss Transform is significantly faster than the above $O({\rm nnz}({\mathbf{x}})\cdot\varepsilon^{-1}\log(f/\delta))$ time for many reasonable settings of the parameters, e.g., in a number of numerical linear algebra applications in which $1/\delta$ can be exponentially large in $d$. Indeed, the Fast Johnson Lindenstrauss Transform was first used by Sárlos to obtain the first speedups for regression and low rank matrix approximation with relative error. Sárlos used a version of the Fast Johnson Lindenstrauss Transform due to . We will use a slightly different version called the Subsampled Randomized Hadamard Transform, or SRHT for short. Later we will see a significantly faster transform for sparse matrices.

We will not present the proof of Theorem 7, instead relying upon the above intuition. The proof of Theorem 7 can be found in the references listed above.

There are distributions on matrices $S$ with the following properties:

2 Matrix multiplication

In this section we study the approximate matrix product problem.

Before giving the theorem, we need a definition.

We now show that sparse embeddings matrices satisfy the $(\varepsilon,\delta,2)$-JL-moment property. This was originally shown by Thorup and Zhang .

where the second equality uses that $h$ and $\sigma$ are independent, while the third equality uses that ${\bf E}[\sigma(j)\sigma(j^{\prime})]=1$ if $j=j^{\prime}$, and otherwise is equal to .

where the second equality uses the independence of $h$ and $\sigma$, and the third equality uses that since $\sigma$ is $4$-wise independent, in order for ${\bf E}\left[\sigma(j_{1})\sigma(j_{2})\sigma(j_{1}^{\prime})\sigma(j_{2}^{\prime})\right]$ not to vanish, it must be that either

$j_{1}=j_{2}=j_{1}^{\prime}=j_{2}^{\prime}$ or

$j_{1}=j_{2}$ and $j_{1}^{\prime}=j_{2}^{\prime}$ but $j_{1}\neq j_{1}^{\prime}$ or

$j_{1}=j_{1}^{\prime}$ and $j_{2}=j_{2}^{\prime}$ but $j_{1}\neq j_{2}$ or

$j_{1}=j_{2}^{\prime}$ and $j_{1}^{\prime}=j_{2}$ but $j_{1}\neq j_{2}$.

Note that in the last two cases, for ${\bf E}\left[\delta(h(j_{1})=i)\delta(h(j_{2})=i)\delta(h(j^{\prime}_{1})=i^{\prime})\delta(h(j^{\prime}_{2})=i^{\prime})\right]$ not to vanish, we must have $i=i^{\prime}$. The fourth equality and first inequality are based on regrouping the summations, and the sixth inequality uses that $\|{\mathbf{x}}\|_{2}=1$.

3 High probability

() Algorithm 1 outputs a subspace embedding with probability at least $1-\delta$. In expectation step 3 is only run a constant number of times.

4 Leverage scores

Proof: We will use the following matrix Chernoff bound for a sum of random matrices, which is a non-commutative Bernstein bound.

For any $j$, ${\mathbf{z}}_{j}^{T}{\mathbf{z}}_{j}/q_{j}$ is a rank-$1$ matrix with operator norm bounded by $\|{\mathbf{z}}_{j}\|_{2}^{2}/q_{j}\leq k/\beta$. Hence,

We need a version of the Johnson-Lindenstrauss lemma, as follows. We give a simple proof for completeness.

The random variable $\sum_{i=1}^{t}g_{i}^{2}$ is $\chi^{2}$ with $t$ degree of freedom. The following tail bounds are known.

(Lemma 1 of ) Let $g_{1},\ldots,g_{t}$ be i.i.d. $N(0,1)$ random variables. Then for any $x\geq 0$,

Setting $x=\varepsilon^{2}t/16$, we have that

For $t=O(\log n/\varepsilon^{2})$, the lemma follows by a union bound over $i\in[n]$.

Finally, it follows that for all $i\in[n]$,

Hence, $q_{i}\geq(1-\gamma)(1-2\gamma)p_{i}$, which for an appropriate choice of constant $\gamma\in(0,1)$, achieves $q_{i}\geq\beta p_{i}$, as desired.

5 Regression

We formally define the regression problem as follows.

Inspecting the simple proof of Theorem 21 we see that (10) in particular implies

Because of the normal equations, we may apply the Pythagorean theorem,

where the first inequality is the triangle inequality, the second inequality uses the sub-multiplicativity of the spectral norm, and the third inequality uses (12). Rearranging, we have

By the normal equations in the sketch space,

6 Machine precision regression

Here we show how to reduce the dependence on $\varepsilon$ to logarithmic in the regression application, following the approaches in .

and by choosing $\varepsilon_{0}=1/2$, say, $O(\log(1/\varepsilon))$ iterations suffice for this scheme also to attain $\varepsilon$ relative error.

7 Polynomial fitting

We now describe the problem more precisely, starting with a definition.

Vandermonde matrices of dimension $n\times q$ require only $O(n)$ implicit storage and admit $O(n~{}\log^{2}q)$ matrix-vector multiplication time (see, e.g., Theorem 2.11 of ). It is also possible to consider block-Vandermonde matrices as in ; for simplicity we will only focus on the simplest polynomial fitting problem here, in which Vandermonde matrices suffice for the discussion.

Least Absolute Deviation Regression

While least squares regression is arguably the most used form of regression in practice, it has certain non-robustness properties that make it unsuitable for some applications. For example, oftentimes the noise in a regression problem is drawn from a normal distribution, in which case least squares regression would work quite well, but if there is noise due to measurement error or a different underlying noise distribution, the least squares regression solution may overfit this noise since the cost function squares each of its summands.

where the inequality follows by Hölder’s inequality.

where $\zeta\in(0,1]$ can be thought of as a relaxation parameter which will allow for more efficient algorithms.

Using this bound together with independence of the sampled rows,

We have computed ${\bf E}[Z]$ and bounded ${\bf Var}[Z]$ as well as $\max_{i\mid p_{i}<1}Z_{i}$, and can now use strong tail bounds to bound the deviation of $Z$ from its expectation. We use the following tail inequalities.

Moreover, if $Z_{i}-{\bf E}[Z_{i}]\leq\Delta$ for all $i$, we have

It suffices to choose $\mathcal{N}$ so that for all ${\mathbf{y}}\in\mathcal{B}$, there exists a vector ${\mathbf{w}}\in\mathcal{N}$ for which $\|{\mathbf{y}}-{\mathbf{w}}\|_{1}\leq\varepsilon$. Indeed, in this case note that

If ${\mathbf{y}}-{\mathbf{w}}=0$, we are done. Otherwise, suppose $\alpha$ is such that $\alpha\|{\mathbf{y}}-{\mathbf{w}}\|_{1}=1$. Observe that $\alpha\geq 1/\varepsilon$, since, $\|{\mathbf{y}}-{\mathbf{w}}\|_{1}\leq\varepsilon$ yet $\alpha\|{\mathbf{y}}-{\mathbf{w}}\|_{1}=1$.

Then $\alpha({\mathbf{y}}-{\mathbf{w}})\in\mathcal{B}$, and we can choose a vector ${\mathbf{w}}^{2}\in\mathcal{N}$ for which $\|\alpha({\mathbf{y}}-{\mathbf{w}})-{\mathbf{w}}^{2}\|_{1}\leq\varepsilon$, or equivalently, $\|{\mathbf{y}}-{\mathbf{w}}-{\mathbf{w}}^{2}/\alpha\|_{1}\leq\varepsilon/\alpha\leq\varepsilon^{2}$. Hence,

Repeating this argument, we inductively have that

There exists an $\varepsilon$-net $\mathcal{N}$ for which $|\mathcal{N}|\leq(2/\varepsilon)^{d}$.

Then $B(\varepsilon,0)$ is a $d$-dimensional polytope with a ($d$-dimensional) volume denoted $|B(\varepsilon,0)|$. Moreover, $B(1,0)$ and $B(\varepsilon/2,0)$ are similar polytopes, namely, $B(1,0)=(2/\varepsilon)B(\varepsilon/2,0)$. As such, $|B(1,0)|=(2/\varepsilon)^{d}|B(\varepsilon/2,0)|$.

By applying (19) and a union bound over the points in $\mathcal{N}$, and rescaling $\varepsilon$ by a constant factor, we have thus shown the following theorem.

2 The Role of subspace embeddings for L1-Regression

Before discussing the existence of such embeddings, let us see how they can be used to speed up the computation of a well-conditioned basis.

3 Gaussian sketching to speed up sampling

we have that (20) holds with $\zeta=1/(4d)$.

4 Subspace embeddings using cauchy random variables

The Cauchy distribution, having density function $p(x)={1\over\pi}\cdot{1\over 1+x^{2}}$, is the unique $1$-stable distribution. That is to say, if $C_{1},\ldots,C_{M}$ are independent Cauchys, then $\sum_{i\in[M]}\gamma_{i}C_{i}$ is distributed as a Cauchy scaled by $\gamma=\sum_{i\in[M]}|\gamma_{i}|$.

The absolute value of a Cauchy distribution has density function $f(x)=2p(x)={2\over\pi}{1\over 1+x^{2}}$. The cumulative distribution function $F(z)$ of it is

Note also that since $\tan(\pi/4)=1$, we have $F(1)=1/2$, so that $1$ is the median of this distribution.

Although Cauchy random variables do not have an expectation, and have infinite variance, some control over them can be obtained by clipping them. The first use of such a truncation technique in algorithmic applications that we are aware of is due to Indyk .

Consider the event $\mathcal{E}$ that a Cauchy random variable $X$ satisfies $|X|\leq M$, for some parameter $M\geq 2$. Then there is a constant $c>0$ for which $\Pr[\mathcal{E}]\geq 1-{2\over\pi M}$ and ${\bf E}[|X|\mid\mathcal{E}]\leq c\log M,$ where $c>0$ is an absolute constant.

(see “Connection to Auerbach bases” in Section 3.1 of ) There exists a $(d,1,1)$-well-conditioned basis.

For readability, it is useful to separate out the following key lemma that is used in Theorem 36 below. This analysis largely follows that in .

Letting $\mathcal{F}=\wedge_{i=1}^{r}\mathcal{F}_{i}$, we have by another union bound that

We can perform the following manipulation (for an event $\mathcal{A}$, we use the notation $\neg\mathcal{A}$ to denote the occurrence of the complement of $\mathcal{A}$):

and $\Pr[\mathcal{F}_{i,j}]\geq 1-O(1/(C^{\prime}rd))$. Combining these two, we have

for $C^{\prime}>0$ a sufficiently large constant. Plugging (22) into the above,

We thus have, combining (23) with Markov’s inequality,

As $C^{\prime}$ can be chosen sufficiently large, while $C$ is the fixed constant of Lemma 33, we have that

The lemma now follows by appropriately setting the constant $C$ in the lemma statement.

where the last equality used that ${\mathbf{y}}^{1},\ldots,{\mathbf{y}}^{d}$ is an Auerbach basis.

Hence, the statement of the theorem holds with probability at least $9/10$, by a union bound over the events in the dilation and contraction arguments. This concludes the proof.

5 Subspace embeddings using exponential random variables

We now describe a speedup over the previous section using exponential random variables, as in . Other speedups are possible, using , though the results in additionally also slightly improve the sampling complexity. The use of exponential random variables in is inspired by an elegant work of Andoni, Onak, and Krauthgamer on frequency moments .

An exponential distribution has support $x\in[0,\infty)$, probability density function $f(x)=e^{-x}$ and cumulative distribution function $F(x)=1-e^{-x}$. We say a random variable $X$ is exponential if $X$ is chosen from the exponential distribution. The exponential distribution has the following max-stability property.

If $U_{1},\ldots,U_{n}$ are exponentially distributed, and $\alpha_{i}>0\ (i=1,\ldots,n)$ are real numbers, then $\textstyle\max\{\alpha_{1}/U_{1},\ldots,\alpha_{n}/U_{n}\}\simeq\left(\sum_{i\in[n]}\alpha_{i}\right)\left/U\right.$, where $U$ is exponential.

The following lemma shows a relationship between the Cauchy distribution and the exponential distribution.

Let $y_{1},\ldots,y_{d}\geq 0$ be scalars. Let $U_{1},\ldots,U_{d}$ be $d$ independendent exponential random variables, and let $X=(\sum_{i\in[d]}y_{i}^{2}/U_{i}^{2})^{1/2}$. Let $C_{1},\ldots,C_{d}$ be $d$ independent Cauchy random variables, and let $Y=(\sum_{i\in[d]}y_{i}^{2}C_{i}^{2})^{1/2}$. There is a constant $\gamma>0$ for which for any $t>0$.

Proof: We would like the density function $h$ of $y_{i}^{2}C_{i}^{2}$. Letting $t=y_{i}^{2}C_{i}^{2}$, the inverse function is $C_{i}=t^{1/2}/y_{i}$. Taking the derivative, we have ${dC_{i}\over dt}={1\over 2y_{i}}t^{-1/2}$. Letting $f(t)={2\over\pi}{1\over 1+t^{2}}$ be the density function of the absolute value of a Cauchy random variable, we have by the change of variable technique,

We would also like the density function $k$ of $y_{i}^{2}E_{i}^{2}$, where $E_{i}\sim 1/U_{i}$. Letting $t=y_{i}^{2}E_{i}^{2}$, the inverse function is $E_{i}=t^{1/2}/y_{i}$. Taking the derivative, ${dE_{i}\over dt}={1\over 2y_{i}}t^{-1/2}$. Letting $g(t)=t^{-2}e^{-1/t}$ be the density function of the reciprocal of an exponential random variable, we have by the change of variable technique,

We claim that $k(t)\leq h(\gamma t)/\gamma$ for a sufficiently small constant $\gamma>0$. This is equivalent to showing that

which for $\gamma<1$, is implied by showing that

We distinguish two cases: first suppose $t\geq y_{i}^{2}$. In this case, $e^{-y_{i}/t^{1/2}}\leq 1$. Note also that $y_{i}t^{1/2}\leq t^{3/2}/y_{i}$ in this case. Hence, $\gamma^{1/2}y_{i}t^{1/2}\leq\gamma^{1/2}t^{3/2}/y_{i}$. Therefore, the above is implied by showing

which holds for a sufficiently small constant $\gamma\in(0,1)$.

Next suppose $t<y_{i}^{2}$. In this case $y_{i}t^{1/2}>t^{3/2}/y_{i}$, and it suffices to show

Using that $e^{x}\geq x^{2}/2$ for $x\geq 0$, it suffices to show

which holds for a small enough $\gamma\in(0,1)$.

where we made the change of variables $s_{i}=\kappa t_{i}$. Setting $\gamma=\kappa^{1/2}$ completes the proof.

We need a bound on $\Pr[Y\geq t]$, where $Y=(\sum_{i\in[d]}y_{i}^{2}C_{i}^{2})^{1/2}$ is as in Lemma 38.

There is a constant $c>0$ so that for any $r>0$,

Proof: For $i\in[d]$, let $\sigma_{i}\in\{-1,+1\}$ be i.i.d. random variables with $\Pr[\sigma_{i}=-1]=\Pr[\sigma_{i}=1]=1/2$. Let $Z=\sum_{i\in[d]}\sigma_{i}y_{i}C_{i}$. We will obtain tail bounds for $Z$ in two different ways, and use this to establish the lemma.

On the one hand, by the $1$-stability of the Cauchy distribution, we have that $Z\sim\|y\|_{1}C$, where $C$ is a standard Cauchy random variable. Note that this holds for any fixing of the $\sigma_{i}$. The cumulative distribution function of the Cauchy random variable is $F(z)={2\over\pi}\arctan(z).$ Hence for any $r>0$,

and therefore using the Taylor series for $\arctan$ for $r>1$,

On the other hand, for any fixing of $C_{1},\ldots,C_{d}$, we have

If $R\geq 0$ is a random variable with finite variance, and $0<\theta<1$, then

Applying this inequality with $R=Z^{2}$ and $\theta=1/2$, we have

Suppose, towards a contradiction, that $\Pr[Y\geq r\|y\|_{1}]\geq c/r$ for a sufficiently large constant $c>0$. By independence of the $\sigma_{i}$ and the $C_{i}$, by (26) this implies

By (25), this is a contradiction for $c>{24\over\pi}$. It follows that $\Pr[Y\geq r\|y\|_{1}]<c/r$, as desired.

Let $y_{1},\ldots,y_{d}\geq 0$ be scalars. Let $U_{1},\ldots,U_{d}$ be $d$ independendent exponential random variables, and let $X=(\sum_{i\in[d]}y_{i}^{2}/U_{i}^{2})^{1/2}$. There is a constant $c>0$ for which for any $r>0$,

Proof: The corollary follows by combining Lemma 38 with Lemma 39, and rescaling the constant $c$ from Lemma 39 by $1/\gamma$, where $\gamma$ is the constant of Lemma 38.

For the dilation, we need Khintchine’s inequality.

(). Let $Z=\sum_{i=1}^{r}\sigma_{i}z_{i}$ for i.i.d. random variables $\sigma_{i}$ uniform in $\{-1,+1\}$, and $z_{1},\ldots,z_{r}$ be scalars. There exists a constant $c>0$ for which for all $t>0$

which we denote by event $\mathcal{E}$ and condition on. Notice that the probability is taken only over the choice of the $\sigma_{i}$, and therefore conditions only the $\sigma_{i}$ random variables.

In the following, $i\in[d]$ and $j\in[r]$. Let $\mathcal{F}_{i,j}$ be the event that

Letting $p=\Pr[\mathcal{E}\wedge\mathcal{F}_{i,j}]\geq 99/100$, we have by Corollary 40,

We can perform the following manipulation:

It follows by linearity of expectation that,

Consequently, by a Markov bound, and using that $p\geq 1/2$, conditioned on $\mathcal{E}\wedge\mathcal{F}$, with probability at least $9/10$, we have the occurrence of the event $\mathcal{G}$ that

where the first inequality uses the triangle inequality, the second the occurrence of $\mathcal{E}$, and the third (28). This completes the proof.

Proof: The corollary follows by combining Theorem 30, Lemma 32 and its optimization in §3.3, and Theorem 41.

6 Application to hyperplane fitting

For the subspace approximation problem, we can write

Low Rank Approximation

We will show how to use sketching to speed up algorithms for both problems, and further variants. Our exposition is based on combinations of several works in this area by Sárlos, Clarkson, and the author .

Section Overview: In §4.1 we give an algorithm for computing a low rank approximation achieving error proportional to the Frobenius norm. In §4.2 we give a different kind of low rank approximation, called a CUR decomposition, which computes a low rank approximation also achieving Frobenius norm error but in which the column space equals the span of a small subset of columns of the input matrix, while the row space equals the span of a small subset of rows of the input matrix. A priori, it is not even clear why such a low rank approximation should exist, but we show that it not only exists, but can be computed in nearly input sparsity time. We also show that it can be computed deterministically in polynomial time. This algorithm requires several detours into a particular kind of spectral sparsification given in §4.2.1, as well as an adaptive sampling technique given in §4.2.2. Finally in §4.2.3 we show how to put the pieces together to obtain the overall algorithm for CUR factorization. One tool we need is a way to compute the best rank-$k$ approximation of the column space of a matrix when it is restricted to lie within a prescribed subspace; we defer the details of this to §4.4, where the tool is developed in the context of an application called Distributed Low Rank Approximation. In §4.3 we show how to perform low rank approximation with a stronger guarantee, namely, an error with respect to the spectral norm. While the solution quality is much better than in the case of the Frobenius norm, it is unknown how to compute this as quickly, though one can still compute it much more quickly than the SVD. In §4.4 we present the details of the Distributed Low Rank Approximation algorithm.

Rescaling $\varepsilon$ by a constant factor completes the proof.

Fortunately, we can cast this projection problem as a regression problem, and solve it approximately.

implying the first part of the theorem after rescaling $\varepsilon$ by a constant factor.

For the second part of the theorem, note that Lemma 45 gives the stronger guarantee that

2 CUR decomposition

We now outline the approach of Boutsidis and the author . A key lemma we need is the following, which is due to Boutsidis, Drineas, and Magdon-Ismail .

as needed. The lemma now follows by Markov’s bound.

To proceed, we need an algorithm in the next subsection, which uses a method of Batson, Spielman, and Srivastava refined for this application by Boutsidis, Drineas, and Magdon-Ismail .

The following theorem shows correctness of the Deterministic Dual Set Spectral Sparsification algorithm described in Algorithm 2.

We now turn to correctness. The crux of the analysis turns out to be to show there always exists an index $j$ in each iteration for which

We start with a lemma which uses the Sherman-Morrison-Woodbury identity to analyze a rank-$1$ perturbation.

For the second part, we use the following well-known formula.

Letting $L^{\prime}=L+\delta_{LOW}$, we have

We also need the following lemma concerning properties of the $UP$ function.

or equivalently, $t{\mathbf{a}}^{T}{\mathbf{a}}\leq\delta_{UP}$. Hence,

Equipped with Lemma 51 and Lemma 52, we now prove the main lemma we need.

At every iteration $\tau=0,\ldots,r-1$, there exists an index $j\in\{1,2,\ldots,n\}$ for which

Indeed, if we show (33), then by averaging there must exist an index $j$ for which

We first prove the equality in (33) using the definition of $\delta_{UP}$. Observe that it holds that

We will show $\mathcal{E}\geq 0$ below. Given this, we have

where the inequality uses Lemma 51. Since $\delta_{LOW}=1$, we have

We now turn to the task of showing $\mathcal{E}\geq 0$. The Cauchy-Schwarz inequality implies that for $a_{i},b_{i}\geq 0$, one has $(\sum_{i}a_{i}b_{i})^{2}\leq(\sum_{i}a_{i}^{2}b_{i})(\sum_{i}b_{i})$, and therefore

Plugging into (34), we conclude that $\mathcal{E}\geq 0$, as desired.

By Lemma 53, the algorithm is well-defined, finding a $t\geq 0$ at each iteration (note that $t\geq 0$ since $t^{-1}\geq UP(a_{j},\delta_{UP})\geq 0$).

Finally, note that Algorithm 2 runs in $r$ steps. The vector $s$ of weights is initialized to the all-zero vector, and one of its entries is updated in each iteration. Thus, $s$ will contain at most $r$ non-zero weights upon termination. As shown above, the value $t$ chosen in each iteration is non-negative, so the weights in $s$ are non-negative.

We will also need the following corollary, which shows how to perform the dual set sparsification much more efficiently if we allow it to be randomized.

Lemma 55 gives us a way to find $4k$ columns providing an $O(1)$-approximation. We would like to refine this approximation to a $(1+\varepsilon)$-approximation using only an additional $O(k/\varepsilon)$ number of columns. To do so, we perform a type of residual sampling from this $O(1)$-approximation, as described in the next section.

2.2 Adaptive sampling

For $i=1,\ldots,n$, and some fixed constant $\alpha>0,$ let $p_{i}$ be a probability distribution such that for each $i:$

For $i=1,\ldots,n$ let $p_{i}$ be a probability distribution such that for each $i:$

To analyze the expected error of the algorithm with respect to the choices made in the sampling procedure, we have

where the second equality uses the Pythagorean theorem.

2.3 CUR wrapup

Hence, for a parameter $c_{2}>0$, if we set

3 Spectral norm error

We also collect a few facts about the singular values of a Gaussian matrix.

In order to analyze SubspacePowerMethod, we need a key lemma shown in concerning powering of a matrix.

If we raise both sides to the $1/(2(2q+1))$-th power, then this completes the proof.

We can now prove the main theorem about SubspacePowerMethod

4 Distributed low rank approximation

The main algorithm AdaptiveCompress of is given in Algorithm AdaptiveCompress below.

Graph Sparsification

We formally define the problem as follows, following the notation and outlines of . Consider an ordering on the $n$ vertices, denoted $1,\ldots,n$. We will only consider undirected graphs, though we will often talk about edges $e=\{u,v\}$ as $e=(u,v)$, where here $u$ is less than $v$ in the ordering we have placed on the edges. This will be for notational convenience only; the underlying graphs are undirected.

We call $H$ a spectral sparsifier of $G$ The usual notation for (57) is

Notice that Theorem 64 shows that if one knows the leverage scores, then by sampling $O(n\varepsilon^{-2}\log n)$ edges of $G$ and reweighting them, one obtains a spectral sparsifier of $G$. One can use algorithms for approximating the leverage scores of general matrices , though more efficient algorithms, whose overall running time is near-linear in the number of edges of $G$, are known .

A beautiful theorem of Kapralov, Lee, Musco, Musco, and Sidford is the following .

In the remainder of the section, we give an outline of the proof of Theorem 65, following the exposition given in . We restrict to unweighted graphs for the sake of presentation; the arguments generalize in a natural way to weighted graphs.

The proof of the following theorem is elementary. We believe the power in the theorem is its novel use in algorithm design.

For the second condition, for all ${\mathbf{x}}$,

Finally, for the third condition, for all ${\mathbf{x}}$,

The bounds on the eigenvalues of a Laplacian are given in (the bound on the maximum eigenvalue follows from the fact that $n$ is the maximum eigenvalue of the Laplacian of the complete graph on $n$ vertices. The bound on the minimum eigenvalue follows from Lemma 6.1 of ).

To do this, first observe that the leverage score for a potential edge $i=(u,v)$ is given by

Then, by the second property of Theorem 66,

Thus, the only task left is to implement this hierarchy of leverage score sampling using linear sketches.

For this, we need the following standard theorem from the sparse recovery literature.

Several standard consequences of this theorem, as observed in , can be derived by setting $\eta={\varepsilon\over C\log n}$ for a constant $C>0$, which is the setting of $\eta$ we use throughout. Of particular interest is that for $0<\varepsilon<1/2$, from $w_{i}$ one can determine if ${\mathbf{x}}_{i}\geq{1\over C\log n}\|{\mathbf{x}}\|_{2}$ or ${\mathbf{x}}_{i}<{1\over 2C\log n}\|{\mathbf{x}}\|_{2}$ given that it satisfies one of these two conditions. We omit the proof of this fact which can be readily verified from the statement of Theorem 67, as shown in .

Sketching Lower Bounds for Linear Algebra

While sketching, and in particular subspace embeddings, have been used for a wide variety of applications, there are certain limitations. In this section we explain some of them.

While our focus in this section is on lower bounds, we mention that for integers $p\geq 1$, there is the following simple algorithm for estimating Schatten norms which has a good running time but requires multiple passes over the data. This is given in .

Proof: Let $r=C/\varepsilon^{2}$ for a positive constant $C>0$. Suppose ${\mathbf{g}}^{1},\ldots,{\mathbf{g}}^{r}$ are independent $N(0,1)^{d}$ vectors, that is, they are independent vectors of i.i.d. normal random variables with mean and variance $1$.

where we use that ${\bf E}[({\mathbf{h}}^{i})_{i}^{2}]=1$ for all $i$. We also have that

This shows correctness. The running time follows from our bound on $r$ and the number $s$ of passes.

2 Sketching the operator norm

The idea is to use the min-max principle for singular values.

Proof: The min-max principle for singular values says that

The lemma now follows from the min-max principle for singular values, since every vector in the range has its norm preserved up to a factor of $1+\varepsilon$, and so this also holds for any $i$-dimensional subspace of the range, for any $i$.

$\mu$ is the distribution on $n\times n$ matrices with i.i.d. $N(0,1)$ entries.

It follows with probability $1-O(1/n)$, by the triangle inequality

In our proof we need the following tail bound due to Latała Suppose that $g_{i_{1}},\dots,g_{i_{d}}$ are i.i.d. $N(0,1)$ random variables. The following result, due to Latała , bounds the tails of Gaussian chaoses $\sum a_{i_{1}}\cdots a_{i_{d}}g_{i_{1}}\cdots g_{i_{d}}$. The proof of Latała’s tail bound was later simplified by Lehec .

where $C_{d},c_{d}>0$ are constants depending only on $d$.

To prove the main theorem of this section, we need a few facts about distances between distributions.

In general $D_{\phi}(\mu||\nu)$ is not a distance because it is not symmetric.

The total variation distance between $\mu$ and $\nu$, denoted by $d_{TV}(\mu,\nu)$, is defined as $D_{\phi}(\mu||\nu)$ for $\phi(x)=|x-1|$. It can be verified that this is indeed a distance. It is well known that if $d_{TV}(\mu,\nu)\leq c<1$, then the probability that any, possibly randomized algorithm, can distinguish the two distributions is at most $(1+c)/2$.

The $\chi^{2}$-divergence between $\mu$ and $\nu$, denoted by $\chi^{2}(\mu||\nu)$, is defined as $D_{\phi}(\mu||\nu)$ for $\phi(x)=(x-1)^{2}$ or $\phi(x)=x^{2}-1$. It can be verified that these two choices of $\phi$ give exactly the same value of $D_{\phi}(\mu||\nu)$.

We can upper bound total variation distance in terms of the $\chi^{2}$-divergence using the next proposition.

$d_{TV}(\mu,\nu)\leq\sqrt{\chi^{2}(\mu||\nu)}$.

The next proposition gives a convenient upper bound on the $\chi^{2}$-divergence between a Gaussian distribution and a mixture of Gaussian distributions.

We can now prove the main impossibility theorem about sketching the operator norm up to a constant factor.

Partition of size 1. The only possible partition is $\{1,2,3,4\}$. We have

Partitions of size 3. Up to symmetry, there are only two partitions to consider: $\{1,2\},\{3\},\{4\}$, and $\{1,3\},\{2\},\{4\}$. We first consider the partition $\{1,2\},\{3\},\{4\}$. We have

where the first inequality follows from Cauchy-Schwarz. We now consider the partition $\{1,3\},\{2\},\{4\}$. We have

where the second equality follows by Cauchy-Schwarz.

Partition of size 4. The only partition is $\{1\},\{2\},\{3\},\{4\}$. Using that for integers $a,b$, $a\cdot b\leq(a^{2}+b^{2})/2$, we have

Latała’s inequality (Theorem 72) states that for $t>0$,

The above holds with no conditions imposed on ${\mathbf{u}},{\mathbf{v}},{\mathbf{u}}^{\prime},{\mathbf{v}}^{\prime}$. For convenience, we let

Note that conditioned on $\mathcal{E}_{{\mathbf{u}},{\mathbf{v}}}$ and $\mathcal{E}_{{\mathbf{u}}^{\prime},{\mathbf{v}}^{\prime}}$,

We now claim that $tD={25t\over n}\leq cf(t)/2$ for all $\sqrt{k}<t<16k$, provided $k=o(n^{2})$. First, note that since $t<16k$, ${t\over\sqrt{k}}=O(\sqrt{t})$. Also, since $t>\sqrt{k}$, ${t\over\sqrt{k}}\leq{t^{2}\over k}$. Hence, $f(t)=\Theta(\min(t/\sqrt{k},t^{2/3}/n^{1/3}))$. Since $k=o(n^{2})$, if $t/\sqrt{k}$ achieves the minimum, then it is larger than ${25t\over n}$. On the other hand, if $t^{2/3}/n^{1/3}$ achieves the minimum, then $cf(t)/2\geq{25t\over n}$ whenever $t=o(n^{2})$, which since $t<16k=o(n^{2})$, always holds.

3 Streaming lower bounds

In this section, we explain some basic communication complexity, and how it can be used to prove bit lower bounds for the space complexity of linear algebra problems in the popular streaming model of computation. We refer the reader to Muthukrishnan’s survey for a comprehensive overview on the streaming model. We state the definition of the model that we need as follows. These results are by Clarkson and the author , and we follow the exposition in that paper.

4 Communication complexity

For lower bounds in the turnstile model, we use a few definitions and basic results from two-party communication complexity, described below. We refer the reader to the book by Kushilevitz and Nisasn for more information . We will call the two parties Alice and Bob.

For a function $f:\mathcal{X}\times\mathcal{Y}\rightarrow\{0,1\}$, we use $R_{\delta}^{1-way}(f)$ to denote the randomized communication complexity with two-sided error at most $\delta$ in which only a single message is sent from Alice to Bob. Here, only a single message $M(X)$ is sent from Alice to Bob, where $M$ is Alice’s message function of her input $X$ and her random coins. Bob computes $f(M(X),Y)$, where $f$ is a possibly randomized function of $M(X)$ and Bob’s input $Y$. For every input pair $(X,Y)$, Bob should output a correct answer with probability at least $1-\delta$, where the probability is taken over the joint space of Alice and Bob’s random coin tosses. If this holds, we say the protocol is correct. The communication complexity $R_{\delta}^{1-way}(f)$ is then the minimum over correct protocols, of the maximum length of Alice’s message $M(X)$, over all inputs and all settings to the random coins.

We also use $R_{\mu,\delta}^{1-way}(f)$ to denote the minimum communication of a protocol, in which a single message from Alice to Bob is sent, for solving $f$ with probability at least $1-\delta$, where the probability now is taken over both the coin tosses of the protocol and an input distribution $\mu$.

In the augmented indexing problem, which we call $AIND$, Alice is given ${\mathbf{x}}\in\{0,1\}^{n}$, while Bob is given both an $i\in[n]$ together with ${\mathbf{x}}_{i+1},{\mathbf{x}}_{i+2},\ldots,{\mathbf{x}}_{n}$. Bob should output ${\mathbf{x}}_{i}$.

() $R_{1/3}^{1-way}(AIND)=\Omega(n)$ and also $R_{\mu,1/3}^{1-way}(AIND)=\Omega(n)$, where $\mu$ is uniform on $\{0,1\}^{n}\times[n]$.

We start with an example showing how to use Theorem 74 for proving space lower bounds for the Matrix Product problem, which is the same as given in Definition 11. Here we also include in the definition the notions relevant for the streaming model.

Proof: Throughout we shall assume that $1/\varepsilon$ is an integer, and that $c$ is an even integer. These conditions can be removed with minor modifications. Let $Alg$ be a $1$-pass algorithm which solves Matrix Product with probability at least $4/5$. Let $r=\log_{10}(cn)/(8\varepsilon^{2})$. We use $Alg$ to solve instances of $AIND$ on strings of size $cr/2$. It will follow by Theorem 74 that the space complexity of $Alg$ must be $\Omega(cr)=\Omega(c\log(cn))/\varepsilon^{2}$.

It follows that the parties can solve the AIND problem with probability at least $4/5-25/198>2/3$. The theorem now follows by Theorem 74.

4.2 Regression and low rank approximation

One can similarly use communication complexity to obtain lower bounds in the streaming model for Regression and Low Rank Approximation. The results are again obtained by reduction from Theorem 74. They are a bit more involved than those for matrix product, and so we only state several of the known theorems regarding these lower bounds. We begin with the formal statement of the problems, which are the same as defined earlier, specialized here to the streaming setting.

() Let $\varepsilon>0$ and $k\geq 1$ be arbitrary. Then,

5 Subspace embeddings

6 Adaptive algorithms

In this section we would like to point out a word of caution of using a sketch for multiple, adaptively chosen tasks.

with probability at least $1-\delta$ for some parameter $\delta>0$. In this section we will focus on the case in which $\delta<n^{-c}$ for every constant $c>0$, that is, $\delta$ shrinks faster than any inverse polynomial in $n$

Further, the algorithm runs in $O(k^{3})$ time and the first $r-1$ queries can be chosen non-adaptively (so the algorithm makes a single adaptive query, namely, ${\mathbf{x}}^{r}$).

Proof: The algorithm first queries the sketch on the vectors

We state some of the intuition behind the proof of Theorem 83 below.

Formally, Theorem 83 makes use of a conditional expectation lemma showing that there exists a choice of $\tau$ for which

Open Problems

We have attempted to cover a number of examples where sketching techniques can be used to speed up numerical linear algebra applications. We could not cover everything and have of course missed out on some great material. We encourage the reader to look at other surveys in this area, such as the one by Mahoney , for a treatment of some of the topics that we missed.

References