2.1 Matrix multiplication using sampling

Suppose $A$ is an $m\times n$ matrix and $B$ is an $n\times p$ matrix, and the product $AB$ is desired. We can use sampling to get an approximate product faster than the traditional multiplication. Let $A(:,k)$ denote the $k$ th column of $A$ . $A(:,k)$ is an $m\times 1$ matrix. Let $B(k,:)$ be the $k$ th row of $B$ . $B(k,:)$ is a $1\times n$ matrix. We have

Can we estimate the sum over all $k$ by summing over a sub-sample? Let $p_{1},p_{2},\ldots ,p_{n}$ be non-negative reals summing to 1 (to be determined later). Let $z$ be a random variable that takes values in $\{1,2,\ldots ,n\}$ with Pr $(z=j)=p_{j}$ . Define an associated matrix random variable $X$ such that

Let $E(X)$ denote the entry-wise expectation, that is,

The scaling by $1/p_{k}$ makes $X$ an unbiased estimator of $AB$ . We will be interested in $E(\Vert AB-X\Vert _{F}^{2})$ , which is just the sum of the variances of all entries of $X$ ,Footnote ² since

We can ignore the $\Vert AB\Vert _{F}^{2}$ term since it does not depend on the $p_{k}$ . Now

It can be seen by calculusFootnote ³ that the minimizing $p_{k}$ must be proportional to $|A(:,k)\Vert B(k,:)|$ . In the important special case when $B=A^{T}$ , this means picking columns of $A$ with probabilities proportional to the squared length of the columns. In fact, even in the general case when $B\neq A^{T}$ , doing so simplifies the bounds, so we will use it. If $p_{k}$ is proportional to $|A(:,k)|^{2}$ , that is,

then

To reduce the variance, we can take the average of $s$ independent trials. Each trial $i$ , $i=1,2,\ldots ,s$ yields a matrix $X_{i}$ as in (2.1). We take

as our estimate of $AB$ . Since the variance of the sum of independent random variables is the sum of the variances, we find that

Let $k_{1},\ldots ,k_{s}$ be the $k$ chosen in each trial. Expanding this, we obtain

We write this as the product of an $m\times s$ matrix with an $s\times p$ matrix as follows. Let $C$ be the $m\times s$ matrix consisting of the following columns, which are scaled versions of the chosen columns of $A$ :

Note that this scaling has a nice property, which we leave to the reader to verify:

Define $R$ to be the $s\times p$ matrix with the corresponding rows of $B$ similarly scaled, namely, $R$ has rows

The reader may also verify that

From (2.2), we see that

This is represented in Figure 2.1. We summarize our discussion in Theorem 2.1.

Figure 2.1.

Approximate matrix multiplication using sampling.

When is this the error bound useful? Let us focus on the case $B=A^{T}$ so that we have just one matrix to consider. If $A$ is the identity matrix, then the guarantee is not very good. In this case, $\Vert AA^{T}\Vert _{F}^{2}=n$ , but the right-hand side of the inequality is $n^{2}/s$ . So we would need $s>n$ for the bound to be any better than approximating the product with the zero matrix.

More generally, the trivial estimate of the zero matrix for $AA^{T}$ makes an error in the Frobenius norm of $\Vert AA^{T}\Vert _{F}$ . If $\unicode[STIX]{x1D70E}_{1},\unicode[STIX]{x1D70E}_{2},\ldots$ are the singular values of $A$ , then the singular values of $AA^{T}$ are $\unicode[STIX]{x1D70E}_{1}^{2},\unicode[STIX]{x1D70E}_{2}^{2},\ldots ,$ and we have

So from the theorem we can assert

provided that

If $\operatorname{rank}(A)=r$ , then there are $r$ non-zero $\unicode[STIX]{x1D70E}_{t}$ and the best general upper bound on the ratio $(\unicode[STIX]{x1D70E}_{1}^{2}+\unicode[STIX]{x1D70E}_{2}^{2}+\cdots )^{2}/(\unicode[STIX]{x1D70E}_{1}^{4}+\unicode[STIX]{x1D70E}_{2}^{4}+\cdots )$ is $r$ , so in general, $s$ needs to be at least $r$ . If $A$ is of full rank, this means that sampling will not give us any gain over taking the whole matrix!

However, if there is a constant $c$ and a small integer $p$ such that

then

and so $s\geq p/c^{2}$ gives us a better estimate than the zero matrix. Further increasing $s$ by a factor decreases the error by the same factor. The condition (2.5) (in words, the top $p$ singular values make up a constant fraction of the spectrum) is indeed the hypothesis of the subject of principal component analysis, and there are many situations when the data matrix does satisfy the condition and so sampling algorithms are useful.

2.2 Sketch of a large matrix

The main result of this section is that for any matrix, a sample of columns and rows, each picked according to length-squared distribution, provides a good sketch of the matrix in a formal sense that will be described briefly. Let $A$ be an $m\times n$ matrix. Pick $s$ columns of $A$ according to length-squared distribution. Let $C$ be the $m\times s$ matrix containing the picked columns scaled so as to satisfy (2.3), that is, if $A(:,k)$ is picked, it is scaled by $1/\sqrt{sp_{k}}$ . Similarly, pick $r$ rows of $A$ according to length-squared distribution on the rows of $A$ . Let $R$ be the $r\times n$ matrix of the picked rows, scaled as follows. If row $k$ of $A$ is picked, it is scaled by $1/\sqrt{rp_{k}}$ . We then have $E(R^{T}R)=A^{T}A$ . From $C$ and $R$ , we can find a matrix $U$ so that $A\approx CUR$ .

One may recall that the top $k$ singular vectors give a similar picture. But the SVD takes more time to compute, requires all of $A$ to be stored in RAM, and does not have the property that the singular vectors, the basis of the reduced space, are directly from $A$ . The last property – that the approximation involves actual rows/columns of the matrix rather than linear combinations – is called an interpolative approximation, and is useful in many contexts. Some structural results of such approximations are found in the work of Stewart (Stewart Reference Stewart1999, Stewart Reference Stewart2004, Berry, Pulatova and Stewart Reference Berry, Pulatova and Stewart2004) and Goreinov, Tyrtyshnikov and Zamarashkin (Goreinov, Tyrtyshnikov and Zamarashkin Reference Goreinov, Tyrtyshnikov and Zamarashkin1997, Goreinov and Tyrtyshnikov Reference Goreinov and Tyrtyshnikov2001).

We briefly mention two motivations for such a sketch. Suppose $A$ is the document-term matrix of a large collection of documents. We are to ‘read’ the collection at the outset and store a sketch so that later, when a query represented by a vector with one entry per term arrives, we can find its similarity to each document in the collection. Similarity is defined by the dot product. In Figure 1.1 it is clear that the matrix–vector product of a query with the right-hand side can be done in time $O(ns+sr+rm)$ , which would be linear in $n$ and $m$ if $s$ and $r$ are $O(1)$ . To bound errors for this process, we need to show that the difference between $A$ and the sketch of $A$ has small 2-norm. The fact that the sketch is an interpolative approximation means that our approximation essentially consists a subset of documents and a subset of terms, which may be thought of as a representative set of documents and terms. Moreover, if $A$ is sparse in its rows and columns – each document contains only a small fraction of the terms and each term is in only a small fraction of the documents – then this property will be preserved in $C$ and $R$ , unlike with the SVD.

A second motivation comes from recommendation systems. Here $A$ would be a customer–product matrix whose $(i,j)$ th entry is the preference of customer $i$ for product $j$ . The objective is to collect a few sample entries of $A$ and, based on these, get an approximation to $A$ so that we can make future recommendations. A few sampled rows of $A$ (all preferences of a few customers) and a few sampled columns (all customer preferences for a few products) give a good approximation to $A$ provided that the samples are drawn according to the length-squared distribution.

It now remains to describe how to find $U$ from $C$ and $R$ . There is an $n\times n$ matrix $P$ of the form $P=QR$ which acts as the identity on the space spanned by the rows of $R$ and zeros out all vectors orthogonal to this space. The matrix $Q$ is just the pseudo-inverse of $R$ .

We begin with some intuition. In particular, we first present a simpler idea that does not work but will then motivate the idea that does work. Write $A$ as $AI$ , where $I$ is the $n\times n$ identity matrix. Now, let us approximate the product $AI$ using the algorithm of Theorem 2.1 from the previous section, that is, by sampling $s$ columns of $A$ according to their squared lengths. Then, as in the last section, write $AI\approx CW$ , where $W$ consists of a scaled version of the $s$ rows of $I$ corresponding to the $s$ columns of $A$ that were picked. Theorem 2.1 bounds the error $\Vert A-CW\Vert _{F}^{2}$ by

But we would like the error to be a small fraction of $\Vert A\Vert _{F}^{2}$ which would require $s\geq n$ , which clearly is of no use since this would pick at least as many columns as the whole of  $A$ .

Instead, let us use the identity-like matrix $P$ instead of $I$ in the above discussion. Using the fact that $R$ is picked according to the squared length, we will show the following proposition later.

We then use Theorem 2.1 to argue that instead of doing the multiplication $AP$ , we can use the sampled columns of $A$ and the corresponding rows of $P$ . The sampled $s$ columns of $A$ form $C$ . We have to take the corresponding $s$ rows of the matrix

This is the same as taking the corresponding $s$ rows of $R^{T}$ and multiplying this by $(\sum _{t=1}^{r^{\prime }}{\unicode[STIX]{x1D70E}_{t}}^{-2\,}\mathbf{u}_{t}\mathbf{u}_{t}^{T})R.$ It is easy to check that this leads to an expression of the form $CUR$ . Moreover, by Theorem 2.1, the error is bounded by

since we will deduce the following inequality.

Combining (2.6) and Proposition 2.3, and using the fact that by the triangle inequality we have

which in turn implies that

we obtain the main result.

We see that if $s$ is fixed, the error is minimized when $r=(s/2)^{2/3}$ . Choosing $s=O(r/\unicode[STIX]{x1D700})$ and $r=1/\unicode[STIX]{x1D700}^{2}$ , the bound becomes $O(\unicode[STIX]{x1D700})\Vert A\Vert _{F}^{2}$ .

Now we prove Proposition 2.3. First,

Let us first suppose that $\mathbf{x}$ is in the row space $V$ of $R$ . We have $P\mathbf{x}=\mathbf{x}$ , so for $\mathbf{x}\in V$ we have $(A-AP)\mathbf{x}=\mathbf{0}$ . Now, since every vector can be written as the sum of a vector in $V$ plus a vector orthogonal to $V$ , this implies that the maximum must therefore occur at some $\mathbf{x}\in V^{\bot }$ . For such $\mathbf{x}$ , we have $(A-AP)\mathbf{x}=A\mathbf{x}$ . Thus, the question now becomes: For unit-length $\mathbf{x}\in V^{\bot }$ , how large can $\Vert A\mathbf{x}\Vert ^{2}$ be? To analyse this, we can write

This implies that we get $\Vert A-AP\Vert _{2}^{2}\leq \Vert A^{T}A-R^{T}R\Vert _{2}$ . So, it suffices to prove that

which follows directly from Theorem 2.1, since we can think of $R^{T}R$ as a way of estimating $A^{T}A$ by picking (according to length-squared distribution) columns of $A^{T}$ , that is, rows of $A$ . By Jensen’s inequality, this implies that

and proves Proposition 2.3.

Proposition 2.4 is easy to see. Since, by Lemma 2.2, $P$ is the identity on the space $V$ spanned by the rows of $R$ , and $P\mathbf{x}=0$ for $\mathbf{x}$ perpendicular to the rows of $R$ , we have that $\Vert P\Vert _{F}^{2}$ is the sum of its singular values squared, which is at most $r$ as claimed.

Finally, to bound the time needed to compute $U$ , the only step involved in computing $U$ is to find the SVD of $R$ . But note that $RR^{T}$ is an $r\times r$ matrix, and since $r$ is much smaller than $n,m$ , this is fast.

2.3 Low-rank approximations and the SVD

Singular value decomposition yields the best approximation of given rank to a matrix in both spectral norm and Frobenius norm. Here, we show that a near-optimal low-rank approximation (LRA) of a matrix $A$ can be found by computing the SVD of a smaller matrix $R$ . The main result is that for any matrix $R$ with

the restriction of $A$ to the space spanned by the top few right singular vectors of $R$ is a good approximation to $A$ . This is purely based on matrix perturbation theory: there is no probability involved.

If we can find an $R$ with smaller dimensions than $A$ , with $R^{T}R\approx A^{T}A$ , then computing the SVD for $R$ would be faster. Indeed, we have already seen one way of getting such an $R$ : pick a set of $r$ rows of $A$ in $r$ i.i.d. trials according to the length-squared distribution and scale them so that $E(R^{T}R)=A^{T}A$ . We saw in Theorem 2.1 that $R^{T}R\approx A^{T}A$ .

The Hoffman–Wielandt inequality (stated below) implies that if $\Vert R^{T}R-A^{T}A\Vert _{F}$ is small, then the singular values of $A$ and $R$ are close. If the top singular vectors of $R$ and $A$ were close, the desired result would follow easily. However, as is well known, multiple singular values are points of discontinuity of singular vectors. But it is not difficult to see intuitively that LRA is not discontinuous: using a singular vector of a slightly worse singular value does not induce gross errors in LRA, and this intuition underlies the result proved formally here.

3.1 Symmetric non-negative tensors with good low-rank approximations

As we mentioned earlier, a small set of samples cannot hope to deal with every matrix. But we saw that sampling yields good results for numerically low-rank matrices (matrices where, say, $\unicode[STIX]{x1D70E}_{1}^{2}(A)+\unicode[STIX]{x1D70E}_{2}^{2}(A)+\cdots +\unicode[STIX]{x1D70E}_{p}^{2}(A)\geq c\Vert A\Vert _{F}^{2}$ for a small $p$ ). Here, we will discuss instances of tensors with a similar property.

For a matrix, the permitted scaling is to scale the rows, scale the columns or both. In this section we only deal with tensors with non-negative real entries which are also symmetric, that is, for any permutation $\unicode[STIX]{x1D70F}$ of $\{1,2,\ldots ,r\}$ , we have

We will show that two natural classes of tensors have Property X. We begin by discussing these classes for matrices. The first class of matrices are adjacency matrices of dense graphs, that is, $n\times n$ matrices with entries $\{0,1\}$ with $\unicode[STIX]{x1D6FA}(n)$ non-zeros in each row. These arise in many applications. They are easily seen to satisfy Property X with no scaling, since

Another important class of matrices represent metrics, where $A_{ij}$ is the distance between points $i$ and $j$ in a metric space, that is, distances satisfy the triangle inequality. Let $D_{i}=\sum _{j=1}^{n}A_{ij}$ be the total distance of point $i$ to other points. Then, one can prove the following result.

We call the lemma the ‘local density property’ since, for the previous example of dense graphs, each entry of the adjacency matrix is indeed at most a constant times the row and column average. Metrics also have this property for each entry.

Using this, one can show that a row and column scaled version of $A$ satisfies Property X. Let

The scaled matrix $B$ is defined by

We will not prove here that $B$ has Property X, but this will follow from the general tensor case below. In fact, we will see that for tensors, a global averaged version of the density lemma holds, and implies Property X.

Let $A$ be any symmetric non-negative tensor. Let

Let $B$ be a scaled version of $A$ :

4.1 Hoeffding–Chernoff inequality for matrix-valued random variables

Notation. For two $n\times n$ real symmetric matrices $B,C$ , we write $B\preccurlyeq C$ when $C-B$ is positive semidefinite.

4.2 Applying Hoeffding–Chernoff to length-squared sampling

Suppose we use length-squared sampling on a matrix $A$ to draw $s$ columns in $s$ i.i.d. trials. Let

be the probability of picking column $j$ . Define the real symmetric matrix-valued random variable $Y$ satisfying

We saw that $EY=AA^{T}/s$ . Then the random variable

Also, we have $E(X^{2})\preccurlyeq E(Y^{2})$ ; this is proved along the same lines as for real random variables, where the variance of a random variable is at most the second moment. Therefore

which implies

We will also need an absolute upper bound on $\Vert X\Vert _{2}$ . For this, note that

5.1 Relative error in the residual

We have seen that if we sample a set of columns of $A$ and scale them to form an $m\times s$ matrix $C$ , the restriction of $A$ to the column space of the top $k$ right singular vectors of $C$ is a good approximation to $A$ , in the sense of Theorem 2.1. But the error is in terms of $\Vert A\Vert _{F}$ . In this section we see how to get errors of the form $O(\Vert A-A_{k}\Vert _{F})$ , where $A_{k}$ denotes the best rank- $k$ approximation to $A$ . Of course, we can get such errors by using the SVD of $A$ . The point here is that we will see that a random sample corresponding to preconditioned length-squared sampling on $A_{k}$ (instead of $A$ ) gives us a subset of columns of $A$ in whose span we find an $O(\Vert A-A_{k}\Vert _{F})$ -type error.

There are at least three approaches to getting such a relative error bound. The first is to use length-squared sampling iteratively (Deshpande et al. Reference Deshpande, Rademacher, Vempala and Wang2006). In this approach, after performing length-squared sampling on the initial matrix $A$ to obtain a sample $S_{1}$ , the columns of $A$ are sampled again, but this time with probability proportional to their squared lengths in the residual, that is, after orthogonal projection onto the span of $S_{1}$ . This process is repeated, and results in a geometrically decaying error factor with each round. Such an approach was also used by Li, Miller and Peng (Reference Li, Miller and Peng2013). More recent papers analyse a similar approach with uniform sampling (Cohen et al. Reference Cohen, Lee, Musco, Musco, Peng and Sidford2015, Cohen, Musco and Musco Reference Cohen, Musco and Musco2017).

The second approach is volume sampling, which extends length-squared sampling by picking subsets of $k$ columns jointly (Deshpande and Vempala Reference Deshpande and Vempala2006). The probability of picking a $k$ -subset is proportional to the squared volume of the $k$ -dimensional simplex induced by the columns along with the origin, a generalization of length-squared sampling. This single sample approach has expected squared error $(k+1)\Vert A-A_{k}\Vert _{F}^{2}$ , the best possible bound using $k$ columns of the matrix. This can be improved by one round to length-squared sampling in the residual to $(1+\unicode[STIX]{x1D716})$ (relative error) using $O(k/\unicode[STIX]{x1D716})$ samples.

The third approach (Drineas et al. Reference Drineas, Mahoney and Muthukrishnan2008), which we present here, uses only one round of length-squared sampling, but the lengths are according to a different matrix, which can be viewed as a preconditioned version of the original matrix, and corresponds to sampling columns according to their leverage scores. Before we go to the main proof of the preconditioned length-squared sampling method, we discuss its application to low-rank approximation in more detail.

Let $S$ be the $n\times s$ column selector (no scaling of columns) matrix, so that $C=AS$ . Let $V_{k}$ be the $n\times k$ matrix with top $k$ right singular vectors of $A$ as columns. The probability of picking column $j$ of $A$ to include in $C$ is according to the squared length of column $j$ of $A_{k}^{+}A_{k}=V_{k}V_{k}^{T}$ , where $A_{k}^{+}$ is the pseudo-inverse of $A$ . Since $V_{k}$ has orthonormal columns, this is the same as the squared length of column $j$ of $V_{k}^{T}$ . Then the approximation to $A$ in the column space of $C$ is

The error is bounded as follows:

Then

Hence $E(\Vert (A-A_{k})S\Vert _{F})\leq \Vert A-A_{k}\Vert _{F}$ . Now, $X_{3}=1/\unicode[STIX]{x1D70E}_{\min }(V_{k}^{T}S)$ (the last $V_{k}^{T}$ can be omitted since it has orthonormal rows). We can bound this as follows:

Noting that $V_{k}V_{k}^{T}=I_{k}$ , we can view this as selecting columns (as in $S$ ) according to length-squared sampling on $V_{k}$ to multiply $V_{k}^{T}$ by $V_{k}$ . By the guarantee for approximate matrix multiplication using length-squared sampling, the error is at most $(2+\unicode[STIX]{x1D700})\Vert A-A_{k}\Vert _{F}$ . This can be improved to $(1+\unicode[STIX]{x1D716})\Vert A-A_{k}\Vert _{F}$ .

5.2 Leverage score sampling

Here we present the main guarantee for sampling according to the leverage scores, that is, the preconditioned squared lengths.

6.1 Very sparse subspace embeddings

Suppose $A$ is an $n\times d$ given data matrix. For this section we assume that $n$ is much larger than $d$ , so $A$ is a tall skinny matrix. For positive reals $a,b$ , we say $a\cong _{\unicode[STIX]{x1D700}}b$ if $a\in ((1-\unicode[STIX]{x1D700})b,(1+\unicode[STIX]{x1D700})b)$ .

We aim to bound $t$ by a function of only $d$ and $\unicode[STIX]{x1D700}$ . The point of subspace embedding is as follows: we can compute with $SA$ instead of $A$ , and if $t\ll n$ , the matrix $SA$ , which is $t\times d$ , is much smaller than $A$ . For example, if we want to find the top singular value (vector) of $A$ within relative error $\unicode[STIX]{x1D700}$ , we can work on $SA$ instead of $A$ . The next theorem is from Clarkson and Woodruff (Reference Clarkson and Woodruff2013).

Subspace embeddings as presented here are due to Clarkson and Woodruff (Reference Clarkson and Woodruff2013) using the work of Dasgupta, Kumar and Sarlós (Reference Dasgupta, Kumar and Sarlós2010). For further improvements and applications, see Woodruff (Reference Woodruff2014).

6.2 Elementary subspace embedding via Johnson–Lindenstrauss

First, we prove a simpler theorem using the Johnson–Lindenstrauss random projection theorem. It says that we can project high-dimensional vectors to a lower-dimensional space and still preserve lengths approximately.

The projection $f:\mathbb{R}^{n}\rightarrow \mathbb{R}^{k}$ that we will examine (in fact, many related projections are known to work as well) is as follows. Pick $k$ vectors $\mathbf{u}_{1},\mathbf{u}_{2},\ldots ,\mathbf{u}_{k}$ in $\mathbb{R}^{n}$ , by choosing all the $nk$ coordinates independently, each from the Gaussian distribution

Then, for any vector $\mathbf{v}$ in $\mathbb{R}^{n}$ , we define the projection $f(\mathbf{v})$ by

We will show that with high probability $\Vert f(\mathbf{v})\Vert \approx \sqrt{k}\Vert \mathbf{v}\Vert$ , so if we have to find, say, $\Vert \mathbf{v}\Vert$ , it suffices to find $\Vert f(\mathbf{v})\Vert$ (since the factor of $\sqrt{k}$ is known). The original proof was to project $\mathbf{v}$ onto a random $k$ -dimensional subspace of $\mathbb{R}^{n}$ . The proof is more complicated for that situation, since projecting onto a random subspace is not equivalent to picking $k$ vectors independently at random and taking dot products. The bound below for a Gaussian random projection is simpler to prove, and still has the same set of applications (Indyk and Motwani Reference Indyk and Motwani1998, Arriaga and Vempala Reference Arriaga and Vempala1999, Dasgupta and Gupta Reference Dasgupta and Gupta2003, Arriaga and Vempala Reference Arriaga and Vempala2006, Vempala Reference Vempala2004).

The theorem deals with one vector. However, subspace embeddings have to work for infinitely many vectors. We cannot do a simple union bound because the number of vectors is infinite. Instead, we use what is called an $\unicode[STIX]{x1D700}$ -net.

The next lemma bounds the size of the $\unicode[STIX]{x1D716}$ -net for a sphere.

Now we may apply Theorem 6.3 as follows. Suppose $V$ is a $d$ -dimensional subspace of $\mathbb{R}^{n}$ , for example, $V=\{A\mathbf{x}\}$ for an $n\times d$ matrix $A$ . The set of unit vectors in $V$ has an $\unicode[STIX]{x1D700}$ -net $N$ of size at most $\text{e}^{cd\ln (1/\unicode[STIX]{x1D700})}$ . Choose $k\geq cd\ln (1/\unicode[STIX]{x1D700})/\unicode[STIX]{x1D700}^{2}$ so that, for a single $\mathbf{x}\in V$ , the probability that $\Vert M\mathbf{x}\Vert \not \cong _{\unicode[STIX]{x1D700}}\sqrt{k}\Vert \mathbf{x}\Vert$ is at most $1/|N|^{2}$ . Then, just by union bounds, we get

Now $f(\mathbf{v})$ can be written as a matrix product, $f(\mathbf{v})=M\mathbf{v}$ , where $M$ is a $k\times n$ matrix with $\mathbf{u}_{1},\mathbf{u}_{2},\ldots ,\mathbf{u}_{k}$ as its rows. We claim that the above suffices to prove that $\Vert M\mathbf{x}\Vert \cong _{3\unicode[STIX]{x1D700}}\sqrt{k}\Vert \mathbf{x}\Vert$ for every $x\in V$ .

Here is a first attempt to prove the claim. For any $\mathbf{x}\in V$ , there is some $\mathbf{y}\in N$ with $\Vert \mathbf{x}-\mathbf{y}\Vert \leq \unicode[STIX]{x1D700}$ . So, $\Vert M\mathbf{x}\Vert =\Vert M\mathbf{y}\Vert +\unicode[STIX]{x1D6FF}$ , where $|\unicode[STIX]{x1D6FF}|\leq \Vert M(\mathbf{x}-\mathbf{y})\Vert$ . But bounding $\Vert M(\mathbf{x}-\mathbf{y})\Vert$ needs a good bound on $\Vert M\Vert _{2}$ which we do not yet have. Instead, we will see that by repeated use of the $\unicode[STIX]{x1D700}$ -net, we can write $\mathbf{x}$ as a convergent series of linear combinations of elements in $N$ .

First, there exists a $\mathbf{y}^{(1)}\in N$ with $\Vert \mathbf{x}-\mathbf{y}^{(1)}\Vert =\unicode[STIX]{x1D703}_{1}\leq \unicode[STIX]{x1D700}$ . Now there exists a $\mathbf{y}^{(2)}\in N$ with

Continuing, we get that $\mathbf{x}=\mathbf{y}^{(1)}+\unicode[STIX]{x1D703}_{1}\mathbf{y}^{(2)}+\unicode[STIX]{x1D703}_{1}\unicode[STIX]{x1D703}_{2}\mathbf{y}^{(3)}+\cdots ,$ where $\unicode[STIX]{x1D703}_{i}\in [0,\unicode[STIX]{x1D700}]$ . Therefore,

A similar argument also proves an upper bound of $1+3\unicode[STIX]{x1D700}$ on $\Vert M\mathbf{x}\Vert$ .

Now we return to the situation where $A$ is an $n\times d$ matrix of rank $r$ . The set $\{A\mathbf{x}:\mathbf{x}\in \mathbb{R}^{d}\}$ is an $r$ -dimensional subspace of $\mathbb{R}^{n}$ . Hence, by the above, as long as $k\geq \operatorname{poly}(r/\unicode[STIX]{x1D700})$ , we deduce that $M$ is a subspace embedding for $A$ . But $M$ is dense, and computing $MA$ is costly (and in general, superlinear).

6.3 Sparse embeddings

Now consider the sparse matrix $S$ defined in Theorem 6.2. Let $V$ be the subspace spanned by the columns of $A$ . Unfortunately, we cannot carry out the same sort of argument as above, because it is not true that for every (individual) $\mathbf{x}\in V$ we have

Counter-examples are ‘lumpy’ $\mathbf{x}$ with some ‘big’ coordinates: for example, if $x_{1}=x_{2}=1/\sqrt{2}$ , with probability $\operatorname{poly}(1/t)$ , we could have $S_{11},S_{12}\not =0$ , whence $\Vert S\mathbf{x}\Vert \not \cong _{\unicode[STIX]{x1D700}}|\mathbf{x}|$ . On the other hand, ‘non-lumpy’ $\mathbf{x}$ – with, say, all $\Vert x_{i}\Vert \leq c\Vert \mathbf{x}\Vert /\sqrt{n}$ – can be shown to satisfy $\text{Pr}(\Vert S\mathbf{x}\Vert \not \cong _{\unicode[STIX]{x1D700}}\Vert \mathbf{x}\Vert )\leq \exp (-cd)$ . We do not do this since it is not required for our proof.

Here is a simple and crucial point.

The lemma asserts that only coordinates in $H$ (for heavy) can be bigger than $\sqrt{\unicode[STIX]{x1D6FC}}\Vert A\mathbf{x}\Vert$ in absolute value. Note that $H$ is independent of $\mathbf{x}$ .

Let $[n]\setminus H=L$ ( $L$ for light). We can write any $\mathbf{y}$ in the column space of $A$ as the sum of two vectors $\mathbf{y}^{H},\mathbf{y}^{L}$ , where $y_{i}^{H}=y_{i}$ for $i\in H$ and $y_{i}^{H}=0$ for $i\in L$ , and similarly $\mathbf{y}^{L}$ is $\mathbf{y}$ restricted to $L$ . Then

We deal with each of these three terms. For each $j\in \{1,2,\ldots ,n\}$ , we picked an $i\in \{1,2,\ldots ,t\}$ at random and set $S_{i,j}=\pm 1$ . It will be convenient to call the $i$ we picked the ‘hash of $j$ ’. Since $n\gg t$ , $j\in \{1,2,\ldots ,n\}$ and $i\in \{1,2,\ldots ,t\}$ , it resembles a random hash function.

For bounding the other two terms of (6.3), the following theorem due to Dasgupta et al. (Reference Dasgupta, Kumar and Sarlós2010) (which we do not prove here) is useful. Call a vector $\mathbf{y}=A\mathbf{x}$ non-lumpy if no coordinate of $\mathbf{y}$ is too big in absolute value. The theorem asserts that with high probability, for all non-lumpy vectors simultaneously, we have bounds on the difference $\Vert S\mathbf{y}\Vert -\Vert \mathbf{y}\Vert$ . Since we have already taken care of the heavy coordinates, this will help us take care of the rest.

Theorem 6.9 can be used to imply immediately that

This takes care of the second term on the right-hand side of (6.3).

The third term requires a further argument. Let

Note that $L^{\prime }$ depends only upon the choices made for $S$ , that is, it does not depend on any particular $y$ . We make three assertions, all holding with high probability:

Statement (6.5) is easy to check. For (6.6), note that for a $j\in L^{\prime }$ conditioned on its hash being one of the $|H|$ hashes of $j^{\prime }\in H$ , it is uniform random. So, $S$ restricted to columns in $L^{\prime }$ is an $|H|\times |L^{\prime }|$ matrix constructed the same way as the whole of $S$ : pick a random row for each column and set the entry in that row and column to $\pm 1$ with probability $1/2$ each. So we may apply Theorem 6.9 to this situation to get (6.6). If the actual $|H|$ is much smaller than the upper bound in (6.2), we can just augment it to get to size $cd\ln d/\unicode[STIX]{x1D700}^{2}$ . Statement (6.7) will follow using a version of Bernstein’s inequality (often called Freedman’s inequality) from probability theory, as follows.

To prove (6.7), consider a single $\mathbf{y}$ first. Since for $j\in L$ , $y_{j}^{L^{\prime }}=y_{j}\mathbf{1}(j\in L^{\prime })$ (where $\mathbf{1}$ is an indicator random variable), we obtain

using (6.2) and (6.4). Thus

Hence Lemma 6.10 with $a=\unicode[STIX]{x1D700}^{2}$ implies that

This probability is small enough that we can take a union bound over an $\unicode[STIX]{x1D700}$ -net and then extend to all vectors again by expressing the vector as a convergent series of net vectors.