Faster Least Squares Approximation

Introduction

where $A^{\dagger}$ denotes the Moore-Penrose generalized inverse of the matrix $A$ . This solution vector has a very natural statistical interpretation as providing an optimal estimator among all linear unbiased estimators, and it has a very natural geometric interpretation as providing an orthogonal projection of the vector $b$ onto the span of the columns of the matrix $A$.

Recall that to minimize the quantity in eqn. (1), we can set the derivative of $\mbox{}\left\|Ax-b\right\|_{2}^{2}=(Ax-b)^{T}(Ax-b)$ with respect to $x$ equal to zero, from which it follows that the minimizing vector $x_{opt}$ is a solution of the so-called normal equations

Geometrically, this means that the residual vector $b^{\perp}=b-Ax_{opt}$ is required to be orthogonal to the column space of $A$, i.e., ${b^{\perp}}^{T}A=0$. While solving the normal equations squares the condition number of the input matrix (and thus is not recommended in practice), direct methods (such as the QR decomposition ) solve the problem of eqn. (1) in $O(nd^{2})$ time assuming that $n\geq d$. Finally, an alternative expression for the vector $x_{opt}$ of eqn. (2) emerges by leveraging the Singular Value Decomposition (SVD) of $A$. If $A=U_{A}\Sigma_{A}V_{A}^{T}$ denotes the SVD of $A$, then

It is worth noting that our second algorithm has a (marginally) less restrictive assumption on the connection between $n$ and $d$. However, the first algorithm is simpler to implement and easier to describe. Clearly, an interesting open problem is to relax the above constraints on $n$ for either of the proposed algorithms.

2 Related work

We should note several lines of related work.

First, techniques such as the “method of averages” preprocess the input into the form of eqn. (6) of Section 3 and can be used to obtain exact or approximate solutions to the least squares problem of eqn. (1) in $o(nd^{2})$ time under strong statistical assumptions on $A$ and $b$. To the best of our knowledge, however, the two algorithms we present and analyze are the first algorithms to provide nontrivial approximation guarantees for overconstrained least squares approximation problems in $o(nd^{2})$ time, while making no assumptions at all on the input data.

Second, Ibarra, Moran, and Hui provide a reduction of the least squares approximation problem to the matrix multiplication problem. In particular, they show that $MM(d)O(n/d)$ time, where $MM(d)$ is the time needed to multiply two $d\times d$ matrices, is sufficient to solve this problem. All of the running times we report in this paper assume the use of standard matrix multiplication algorithms, since $o(d^{3})$ matrix multiplication algorithms are almost never used in practice. Moreover, even with the current best value for the matrix multiplication exponent, $\omega\approx 2.376$ , our algorithms are still faster.

Third, motivated by our preliminary results as reported in and , both Rokhlin and Tygert as well as Avron, Maymounkov, and Toledo have empirically evaluated numerical implementations of variants of one of the algorithms we introduce. We describe this in more detail below in Section 1.3.

Fourth, very recently, Clarkson and Woodruff proved space lower bounds on related problems ; and Nguyen, Do, and Tran achieved a small improvement in the sampling complexity for related problems .

3 Empirical performance of our randomized algorithms

In prior work we have empirically evaluated randomized algorithms that rely on the ideas that we introduce in this paper in several large-scale data analysis tasks. Nevertheless, it is a fair question to ask whether our “random perspective” on linear algebra will work well in numerical implementations of interest in scientific computation. We address this question here. Although we do not provide an empirical evaluation in this paper, in the wake of the original Technical Report version of this paper in 2007 , two groups of researchers have demonstrated that numerical implementations of variants of the algorithms we introduce in this paper can perform very well in practice.

In 2008, Rokhlin and Tygert describe a variant of our random projection algorithm, and they demonstrate that their algorithm runs in time

Their numerical experiments on this class of matrices clearly indicate that their implementations of variants of our algorithms perform well for certain matrices as small as thousands of rows by hundreds of columns.

In 2009, Avron, Maymounkov, Toledo introduced a randomized least-squares solver based directly on our algorithms. They call it Blendenpik, and by considering a much broader class of matrices, they demonstrate that their solver “beats LAPACK’s direct dense least-sqares solver by a large margin on essentially any dense tall matrix.” Beyond providing additional theoretical analysis, including backward error analysis bounds for our algorithm, they consider five (and numerically implement three) random projection strategies (i.e., Discrete Fourier Transform, Discrete Cosine Transform, Discrete Hartely Transform, Walsh-Hadamard Transform, and a Kac random walk), and they evaluate their algorithms on a wide range of matrices of various sizes and various “localization ” or “coherence” properties. Based on these results that empirically show the superior performance of randomized algorithms such as those we introduce and analyze in this paper on a wide class of matrices, they go so far as to “suggest that random-projection algorithms should be incorporated into future versions of LAPACK.”

4 Outline

After a brief review of relevant background in Section 2, Section 3 presents a structural result outlining conditions on preconditioner matrices that are sufficient for relative-error approximation. Then, we present our main sampling-based algorithm for approximating least squares approximation in Section 4 and in Section 5 we present a second projection-based algorithm for the same problem. Preliminary versions of parts of this paper have appeared as conference proceedings in the 17th ACM-SIAM Symposium on Discrete Algorithms and in the 47th IEEE Symposium on Foundations of Computer Science ; and the original Technical Report version of this journal paper has appeared on the arXiv . In particular, the core of our analysis in this paper was introduced in , where an expensive-to-compute probability distribution was used to construct a relative-error approximation sampling algorithm for the least squares approximation problem. Then, after the development of the Fast Johnson-Lindenstrauss transform , proved that similar ideas could be used to improve the running time of randomized algorithms for the least squares approximation problem. In this paper, we have combined these ideas, treated the two algorithms in a manner to highlight their similarities and differences, and considerably simplified the analysis.

Preliminaries

We will make frequent use of matrix and vector norms. More specifically, we let

denote the square of the Frobenius norm of $A$, and we let

2 Linear Algebra background

3 Markov’s inequality and the union bound

We will make frequent use of the following fundamental result from probability theory, known as Markov’s inequality . Let $X$ be a random variable assuming non-negative values with expectation $\mbox{}{\bf{E}}\left[X\right]$. Then, for all $t>0$,

We will also need the so-called union bound. Given a set of random events ${\cal E}_{1},{\cal E}_{2},\ldots,{\cal E}_{n}$ holding with respective probabilities $p_{1},p_{2},\ldots,p_{n}$, the probability that all events hold (i.e., the probability of the union of those events) is upper bounded by $\sum_{i=1}^{n}p_{i}$.

4 The Randomized Hadamard Transform

The Randomized Hadamard Transform was introduced in as one step in the development of a fast version of the Johnson-Lindenstrauss lemma . Recall that the (non-normalized) $n\times n$ matrix of the Hadamard transform $H_{n}$ may be defined recursively as follows:

Our algorithms as preconditioners

Both of our algorithms may be viewed as preconditioning the input matrix $A$ and the target vector $b$ with a carefully-constructed data-independent random matrix $X$. For our random sampling algorithm, we let $X=S^{T}HD$, where $S$ is a matrix that represents the sampling operation and $HD$ is the Randomized Hadamard Transform, while for our random projection algorithm, we let $X=THD$, where $T$ is a random projection matrix. Thus, we replace the least squares approximation problem of eqn. (1) with the least squares approximation problem

We explicitly compute the solution to the above problem using a traditional deterministic algorithm , e.g., by computing the vector

Alternatively, one could use standard iterative methods such as the the Conjugate Gradient Normal Residual method (CGNR, see for details), which can produce an $\epsilon$-approximation to the optimal solution of eqn. (6) in $O(\kappa(XA)rd\ln(1/\epsilon))$ time, where $\kappa(XA)$ is the condition number of $XA$ and $r$ is the number of rows of $XA$.

The two conditions that we will require of the matrix $X$ are:

for some $\epsilon\in(0,1)$. Several things should be noted about these conditions. First, although condition (9) depends on the right hand side vector $b$, Algorithms 1 and 2 will satisfy it without using any information from $b$. Second, although condition (8) only states that $\sigma_{i}^{2}(XU_{A})\geq 1/\sqrt{2}$, for all $i\in[d]$, for both of our randomized algorithms we will show that $\left|1-\sigma_{i}^{2}(XU_{A})\right|\leq 1-2^{-1/2}$, for all $i\in[d]$. Thus, one should think of $XU_{A}$ as an approximate isometry. Third, condition (9) simply states that $Xb^{\perp}=XU_{A}^{\perp}{U_{A}^{\perp}}^{T}b$ remains approximately orthogonal to $XU_{A}$. Finally, note that the following lemma is a deterministic statement, since it makes no explicit reference to either of our randomized algorithms. Failure probabilities will enter later when we show that our randomized algorithms satisfy conditions (8) and (9).

Proof: Let us first rewrite the down-scaled regression problem induced by $X$ as

Thus, by the normal equations (3), we have that

Taking the norm of both sides and observing that under condition (8) we have $\sigma_{i}((XU_{A})^{T}XU_{A})=\sigma_{i}^{2}(XU_{A})\geq 1/\sqrt{2}$, for all $i$, it follows that

To establish the first claim of the lemma, let us rewrite the norm of the residual vector as

where (19) follows since $\sigma_{min}(A)$ is the smallest singular value of $A$ and since the rank of $A$ is $d$; and (20) follows by (15) and the orthogonality of $U_{A}$. Taking the square root, the second claim of the lemma follows. $\diamond$

If we make no assumption on $b$, then (11) from Lemma 1 may provide a weak bound in terms of $\mbox{}\left\|x_{opt}\right\|_{2}$. If, on the other hand, we make the additional assumption that a constant fraction of the norm of $b$ lies in the subspace spanned by the columns of $A$, then (11) can be strengthened. Such an assumption is reasonable, since most least-squares problems are practically interesting if at least some part of $b$ lies in the subspace spanned by the columns of $A$.

Using the notation of Lemma 1 and assuming that $\mbox{}\left\|U_{A}U_{A}^{T}b\right\|_{2}\geq\gamma\mbox{}\left\|b\right\|_{2}$, for some fixed $\gamma\in(0,1]$ it follows that

Proof: Since $\mbox{}\left\|U_{A}U_{A}^{T}b\right\|_{2}\geq\gamma\mbox{}\left\|b\right\|_{2}$, it follows that

This last inequality follows from $U_{A}U_{A}^{T}b=Ax_{opt}$, which implies

By combining this with eqn. (11) of Lemma 1, the lemma follows. $\diamond$

A sampling-based randomized algorithm

In this section, we present our randomized sampling algorithm for the least squares approximation problem of eqn. (1). We also state and prove an associated quality-of-approximation theorem.

Remark: Assuming that $d\leq n\leq e^{d}$, and using $\max\{a_{1},a_{2}\}\leq a_{1}+a_{2}$, we get that

Thus, the running time of Algorithm 1 becomes

Assuming that $\frac{n}{\ln n}=\Omega(d^{2})$, the above running time reduces to

It is worth noting that improvements over the standard $O(nd^{2})$ time could be derived with weaker assumptions on $n$ and $d$. However, for the sake of clarity of presentation, we only focus on the above setting.

Remark: The assumptions in our theorem have a natural geometric interpretation.We would like to thank Ilse Ipsen for pointing out to us this geometric interpretation. In particular, they imply that our approximation becomes worse as the angle between the vector $b$ and the column space of $A$ increases. To see this, let $\mathcal{Z}=||Ax_{opt}-b||_{2}$, and note that $||b||_{2}^{2}=||U_{A}U^{T}_{A}b||_{2}^{2}+\mathcal{Z}^{2}$. Hence the assumption $||U_{A}U^{T}_{k}b||_{2}\geq\gamma||b||_{2}$ can be simply stated as

2 The effect of the Randomized Hadamard Transform

In this subsection, we state a lemma that quantifies the manner in which $HD$ approximately “uniformizes” information in the left singular subspace of the matrix $A$. We state the lemma for a general $n\times d$ orthogonal matrix $U$ such that $U^{T}U=I_{d}$, although we will be interested in the case when $n\gg d$ and $U$ consists of the top $d$ left singular vectors of the matrix $A$.

Let $U$ be an $n\times d$ orthogonal matrix and let the product $HD$ be the $n\times n$ Randomized Hadamard Transform of Section 2.4. Then, with probability at least $.95$,

Proof: We follow the proof of Lemma 2.1 in . In that lemma, the authors essentially prove that the Randomized Hadamard Transform $HD$ “spreads out” input vectors. More specifically, since the columns of the matrix $U$ (denoted by $U^{(j)}$ for all $j\in[d]$) are unit vectors, they prove that for fixed $j\in[d]$ and fixed $i\in[n]$,

From a standard union bound, this immediately implies that with probability at least $1-1/20$,

holds for all $i\in[n]$ and $j\in[d]$ . Using

for all $i\in[n]$, we conclude the proof of the lemma. $\diamond$

3 Satisfying condition (8)

We now establish the following lemma which states that all the singular values of $S^{T}HDU_{A}$ are close to one. The proof of Lemma 4 depends on a bound for approximating the product of a matrix times its transpose by sampling (and rescaling) a small number of columns of the matrix. This bound appears as Theorem 4 in the Appendix and is an improvement over prior work of ours in .

In the above, we used the fact that $U_{A}^{T}DH^{T}HDU_{A}=I_{d}$. We now can view $U_{A}^{T}DSS^{T}H^{T}HDU_{A}$ as an approximation to the product of two matrices $U_{A}^{T}DH^{T}=\left(HDU_{A}\right)^{T}$ and $HDU_{A}$ by randomly sampling and rescaling columns of $\left(HDU_{A}\right)^{T}$. Thus, we can leverage Theorem 4 from the Appendix. More specifically, consider the matrix $\left(HDU_{A}\right)^{T}$. Obviously, since $H$, $D$, and $U_{A}$ are orthogonal matrices, $\mbox{}\left\|HDU_{A}\right\|_{2}=1$ and $\mbox{}\left\|HDU_{A}\right\|_{F}=\mbox{}\left\|U_{A}\right\|_{F}=\sqrt{d}$. Let $\beta=\left(2\ln(40nd)\right)^{-1}$; since we assumed that eqn. (23) holds, we note that the columns of $\left(HDU_{A}\right)^{T}$, which correspond to the rows of $HDU_{A}$, satisfy

Thus, applying Theorem 4 with $\beta$ as above, $\epsilon=1-\left(1/\sqrt{2}\right)$, and $\delta=1/20$ implies that

holds with probability at least $1-1/20=.95$. For the above bound to hold, we need $r$ to assume the value of eqn. (26). Finally, we note that since $\mbox{}\left\|HDU_{A}\right\|_{F}^{2}=d\geq 1$, the assumption of Theorem 4 on the Frobenius norm of the input matrix is always satisfied. Combining the above with inequality (27) concludes the proof of the lemma. $\diamond$

4 Satisfying condition (9)

We next prove the following lemma, from which it will follow that condition (9) is satisfied by Algorithm 1. The proof of this lemma depends on bounds for randomized matrix multiplication algorithms that appeared in .

If eqn. (23) holds and $r\geq 40d\ln(40nd)/\epsilon$, then with probability at least .9,

Proof: Recall that $b^{\perp}=U_{A}^{\perp}{U_{A}^{\perp}}^{T}b$ and that ${\cal Z}=\mbox{}\left\|b^{\perp}\right\|_{2}$. We start by noting that since $\mbox{}\left\|U_{A}^{T}DH^{T}HDb^{\perp}\right\|_{2}^{2}=\mbox{}\left\|U_{A}^{T}b^{\perp}\right\|_{2}^{2}=0$ it follows that

Thus, we can view $\left(S^{T}HDU_{A}\right)^{T}S^{T}HDb^{\perp}$ as approximating the product of two matrices $\left(HDU_{A}\right)^{T}$ and $HDb^{\perp}$ by randomly sampling columns from $\left(HDU_{A}\right)^{T}$ and rows/elements from $HDb^{\perp}$. Note that the sampling probabilities are uniform and do not depend on the norms of the columns of $\left(HDU_{A}\right)^{T}$ or the rows of $\mathcal{H}b^{\perp}$. However, we can still apply the results of Table 1 (second row) in page 150 of . More specifically, since we condition on eqn. (23) holding, the rows of $HDU_{A}$ (which of course correspond to columns of $\left(HDU_{A}\right)^{T}$) satisfy

for $\beta=\left(2\ln(40nd)\right)^{-1}$. Applying the result of Table 1 (second row) of we get

In the above we used $\mbox{}\left\|HDU_{A}\right\|_{F}^{2}=d$. Markov’s inequality now implies that with probability at least .9,

Setting $r\geq 20\beta^{-1}d/\epsilon$ and using the value of $\beta$ specified above concludes the proof of the lemma. $\diamond$

5 Completing the proof of Theorem 2

We now complete the proof of Theorem 2. First, let ${\cal E}_{(\ref{eqn:lem:HU_eqn2})}$ denote the event that eqn. (23) holds; clearly, $\mbox{}{\bf{Pr}}\left[{\cal E}_{(\ref{eqn:lem:HU_eqn2})}\right]\geq.95$. Second, let ${\cal E}_{\ref{lem:sample_lem20pf},\ref{lem:sample_lem40pf}|(\ref{eqn:lem:HU_eqn2})}$ denote the event that both Lemmas 4 and 5 hold conditioned on ${\cal E}_{(\ref{eqn:lem:HU_eqn2})}$ holding. Then,

In the above, $\overline{\cal E}$ denotes the complement of event ${\cal E}$. In the first inequality we used the union bound and in the second inequality we leveraged the bounds for the failure probabilities of Lemmas 4 and 5 given that eqn. (23) holds. We now let ${\cal E}$ denote the event that both Lemmas 4 and 5 hold, without any a priori conditioning on event ${\cal E}_{(\ref{eqn:lem:HU_eqn2})}$; we will bound $\mbox{}{\bf{Pr}}\left[\cal E\right]$ as follows:

In the first inequality we used the fact that all probabilities are positive. The above derivation immediately bounds the success probability of Theorem 2. Combining Lemmas 4 and 5 with the structural results of Lemma 1 and setting $r$ as in eqn. (22) concludes the proof of the accuracy guarantees of Theorem 2.

A projection-based randomized algorithm

In this section, we present a projection-based randomized algorithm for the least squares approximation problem of eqn. (1). We also state and prove an associated quality-of-approximation theorem.

In more detail, Algorithm 2 begins by preprocessing the matrix $A$ and right hand side vector $b$ with the Randomized Hadamard Transform $HD$ of Section 2.4. This algorithm explicitly computes only those rows of $HDA$ and those elements of $HDb$ that need to be accessed to perform the sparse projection. After this initial preprocessing, Algorithm 2 will perform a “sparse projection” by multiplying $HDA$ and $HDb$ by the sparse matrix $T$ (described in more detail in Section 5.2). Then, we can consider the problem

Finally, the expected running time of the algorithm is (at most)

Remark: Assuming that $d\leq n\leq e^{d}$ we get that

Thus, the expected running time of Algorithm 2 becomes

Finally, assuming $n=\Omega(d^{2})$, the above running time reduces to

It is worth noting that improvements over the standard $O(nd^{2})$ time could be derived with weaker assumptions on $n$ and $d$.

2 Sparse projection matrices

In this subsection, we state a lemma about the action of a sparse random matrix operating on a vector. Recall that given any set of $n$ points in Euclidean space, the Johnson-Lindenstrauss lemma states that those points can be mapped via a linear function to $k=O(\epsilon^{-2}\ln n)$ dimensions such that the distances between all pairs of points are preserved to within a multiplicative factor of $1\pm\epsilon$; see and references therein for details.

Formally, let $\epsilon\in(0,1/2)$ be an error parameter, $\delta\in(0,1)$ be a failure probability, and $\alpha\in[1/\sqrt{n},1]$ be a “uniformity” parameter. In addition, let $q$ be a “sparsity” parameter defining the expected number of nonzero elements per row, and let $k$ be the number of rows in our matrix. Then, define the $k\times n$ random matrix $T$ as in Algorithm 2. Matoušek proved the following lemma, as the key step in his version of the Ailon-Chazelle result .

3 Proof of Theorem 3

In this subsection, we provide a proof of Theorem 3. Recall that by the results of Section 3.1, in order to prove Theorem 3, we must show that the matrix $THD$ constructed by Algorithm 2 satisfies conditions (8) and (9) with probability at least $.5$. The next two subsections focus on proving that these conditions hold; the last subsection discusses the running time of Algorithm 2.

In order to prove that all the singular values of $THDU_{A}$ are close to one, we start with the following lemma which provides a means to bound the spectral norm of a matrix. This lemma is an instantiation of lemmas that appeared in .

Let $M$ be a $d\times d$ symmetric matrix and define the grid

In words, $\Omega$ includes all $d$-dimensional vectors $x$ whose coordinates are integer multiples of $\left(2\sqrt{d}\right)^{-1}$ and satisfy $\mbox{}\left\|x\right\|_{2}\leq 1$. Then, the cardinality of $\Omega$ is at most $e^{4d}$. In addition, if for every $x,y\in\Omega$ we have that $\left|x^{T}My\right|\leq\epsilon^{\prime}$, then for every unit vector $x$ we have that $\left|x^{T}Mx\right|\leq 4\epsilon^{\prime}$.

We next establish Lemma 8, which states that all the singular values of $THDU_{A}$ are close to one with constant probability. The proof of this lemma depends on the bound provided by Lemma 7 and it immediately shows that condition (8) is satisfied by Algorithm 2.

Assume that Lemma 3 holds. If $q$ and $k$ satisfy:

holds for all $i\in[d]$. Here $C_{q}$ and $C_{k}$ are the unspecified constants of Lemma 6.

holds for all $i\in[d]$. Consider the grid $\Omega$ of eqn. (32) and note that there are no more than $e^{8d}$ pairs $(x,y)\in\Omega\times\Omega$, since $\left|\Omega\right|\leq e^{4d}$ by Lemma 7. Since $\mbox{}\left\|M\right\|_{2}=\sup_{\mbox{}\left\|x\right\|_{2}=1}\left|x^{T}Mx\right|$, in order to show that $\mbox{}\left\|M\right\|_{2}\leq 1-2^{-1/2}$, it suffices by Lemma 7 to show that $\left|x^{T}My\right|\leq\left(1-2^{-1/2}\right)/4$, for all $x,y\in\Omega$. To do so, first, consider a single $x,y$ pair. Let

By multiplying out the right hand side of the above equation and rearranging terms, it follows that

In order to use Lemma 6 to bound the quantities $\Delta_{1},\Delta_{2}$, and $\Delta_{3}$, we need a bound on the uniformity ratio $\mbox{}\left\|HDU_{A}x\right\|_{\infty}/\mbox{}\left\|HDU_{A}x\right\|_{2}$. To do so, note that

The above inequalities follow by $\mbox{}\left\|HDU_{A}x\right\|_{2}=\mbox{}\left\|x\right\|_{2}$ and Lemma 3. This holds for both our chosen points $x$ and $y$ and in fact for all $x\in\Omega$. Let $\epsilon_{1}=3/125$ and let $\delta=1/(60e^{8d})$ (these choices will be explained shortly). Then, it follows from Lemma 6 that by setting $\alpha=\sqrt{2d\ln(40nd)/n}$ and our choices for $k$ and $q$, each of the following three statements holds with probability at least $1-\delta$:

Thus, combining the above with eqn. (36), for this single pair of vectors $(x,y)\in\Omega\times\Omega$,

holds with probability at least $1-3\delta$. Next, recall that there are no more than $e^{8d}$ pairs of vectors $(x,y)\in\Omega\times\Omega$, and we need eqn. (37) to hold for all of them. Since we set $\delta=1/(60e^{8d})$ then it follows by a union bound that eqn. (37) holds for all pairs of vectors $(x,y)\in\Omega\times\Omega$ with probability at least .95. Additionally, let us set $\epsilon_{1}=3/125$, which implies that $\left|x^{T}My\right|\leq 9/125\leq\left(1-2^{-1/2}\right)/4$ thus concluding the proof of the lemma.

Finally, we discuss the values of the parameters $q$ and $k$. Since $\delta=1/(60e^{8d})$, $\epsilon_{1}=3/125$, and $\alpha=\sqrt{2d\ln(40nd)/n}$, the appropriate values for $q$ and $k$ emerge after elementary manipulations from Lemma 6. $\diamond$

3.2 Satisfying condition (9)

In order to prove that condition (9) is satisfied, we start with Lemma 9. In words, this lemma states that given vectors $x$ and $y$ we can use the random sparse projection matrix $T$ to approximate $\left|x^{T}y\right|$ by $\left|x^{T}T^{T}Ty\right|$, provided that $\mbox{}\left\|x\right\|_{\infty}$ (or $\mbox{}\left\|y\right\|_{\infty}$, but not necessarily both) is bounded. The proof of this lemma is elementary but tedious and is deferred to Section 6.2 of the Appendix.

The following lemma proves that condition (9) is satisfied by Algorithm 2. The proof of this lemma depends on the bound provided by Lemma 9. Recall that $b^{\perp}=U_{A}^{\perp}{U_{A}^{\perp}}^{T}b$ and thus $\mbox{}\left\|b^{\perp}\right\|_{2}=\mbox{}\left\|U_{A}^{\perp}{U_{A}^{\perp}}^{T}b\right\|_{2}={\cal Z}$.

Assume that eqn. (23) holds. If $k\geq 60d/\epsilon$ and $q\geq 2n^{-1}\ln(40nd)$, then, with probability at least .9,

Proof: We first note that since $U_{A}^{T}b^{\perp}=0$, it follows that ${U_{A}^{(j)}}^{T}b^{\perp}={U_{A}^{(j)}}^{T}DH^{T}HDb^{\perp}=0$, for all $j\in[d]$. Thus, we have that

We now bound the expectation of the left hand side of eqn. (38) by using Lemma 9 to bound each term on the right hand side of eqn. (38). Using eqn. (24) of Lemma 3 we get that

holds for all $j\in[d]$. By our choice of the sparsity parameter $q$ the conditions of Lemma 9 are satisfied. It follows from Lemma 9 that

The last line follows since $\mbox{}\left\|\left(HDU_{A}\right)^{(j)}\right\|_{2}=1$, for all $j\in[d]$. Using Markov’s inequality, we get that with probability at least $.9$,

The proof of the lemma is concluded by using the assumed value of $k$. $\diamond$

3.3 Proving Theorem 3

By our choices of $k$ and $q$ as in eqns. (31) and (30), it follows that both conditions (8) and (9) are satisfied. Combining with Lemma 1 we immediately get the accuracy guarantees of Theorem 3. The failure probability of Algorithm 2 can be bounded using an argument similar to the one used in Section 4.5.

References

Appendix

for all $i\in[n]$ for some constant $\beta\in(0,1]$. Let $\epsilon\in(0,1)$ be an accuracy parameter and assume $\mbox{}\left\|A\right\|_{F}^{2}\geq 1/24$. If

then, with probability at least $1-\delta$,

Proof: Consider the Exactly$(c)$ algorithm. Then

for $i\in[n]$. The matrix $C=AS$ has columns $\frac{1}{\sqrt{c}}y^{1},\frac{1}{\sqrt{c}}y^{2},\ldots,\frac{1}{\sqrt{c}}y^{c}$, where $y^{1},y^{2},\ldots,y^{c}$ are $c$ independent copies of $y$. Using this notation, it follows that

We can now apply Lemma 1, p. 3 of . Notice that from eqn. (41) and our assumption on the spectral norm of $A$, we immediately get that

with probability at least $1-\left(2c\right)^{2}\exp\left(-\frac{c\epsilon^{2}}{16M^{2}+8M^{2}\epsilon}\right)$. Let $\delta$ be the failure probability of Theorem 4; we seek an appropriate value of $c$ in order to guarantee $\left(2c\right)^{2}\exp\left(-\frac{c\epsilon^{2}}{16M^{2}+8M^{2}\epsilon}\right)\leq\delta$. Equivalently, we need to satisfy

Recall that $\epsilon<1$, and combine eqns. (42) and (39) to get $M^{2}\leq\mbox{}\left\|A\right\|_{F}^{2}/\beta$. Combining with the above equation, it suffices to choose a value of $c$ such that

We now use the fact that for any $\eta\geq 4$, if $x\geq 2\eta\ln\eta$ then $\frac{x}{\ln x}\geq\eta$. Let $x=2c/\sqrt{\delta}$, let $\eta=96\mbox{}\left\|A\right\|_{F}^{2}/\left(\beta\epsilon^{2}\sqrt{\delta}\right)$, and note that $\eta\geq 4$ if $\mbox{}\left\|A\right\|_{F}^{2}\geq 1/24$, since $\beta$, $\epsilon$, and $\delta$ are at most one. Thus, it suffices to set

which concludes the proof of the theorem. $\diamond$

2 The proof of Lemma 9

Let $t_{(i)}$ be the $i$-th row of $T$ as a row vector, for $i\in[k]$, in which case

Rather than computing $\mbox{}{\bf{E}}\left[\Delta^{2}\right]$ directly, we will instead use that $\mbox{}{\bf{E}}\left[\Delta^{2}\right]=\left(\mbox{}{\bf{E}}\left[\Delta\right]\right)^{2}+\mbox{}{\bf{Var}}\left[\Delta\right]$. We first claim that $\mbox{}{\bf{E}}\left[\Delta\right]=0$. By linearity of expectation,

We first analyze $t_{(i)}=t$ for some fixed $i$ (w.l.o.g. $i=1$). Let $t_{i}$ denote the $i$-th element of the vector $t$ and recall that $\mbox{}{\bf{E}}\left[t_{i}\right]=0$, $\mbox{}{\bf{E}}\left[t_{i}t_{j}\right]=0$ for $i\neq j$, and also that $\mbox{}{\bf{E}}\left[t_{i}^{2}\right]=1/k$. Thus,

By combining the above with eqn. (44), it follows that $\mbox{}{\bf{E}}\left[\Delta\right]=0$, and thus that $\mbox{}{\bf{E}}\left[\Delta^{2}\right]=\mbox{}{\bf{Var}}\left[\Delta\right]$. In order to provide a bound for $\mbox{}{\bf{Var}}\left[\Delta\right]$, note that

Eqn. (45) follows since the $k$ random variables $x^{T}t_{(i)}^{T}t_{(i)}y-\frac{1}{k}x^{T}y$ are independent (since the elements of $T$ are independent) and eqn. (46) follows since $\frac{1}{k}x^{T}y$ is constant. In order to bound eqn. (46), we first analyze $t_{(i)}=t$ for some $i$ (w.l.o.g. $i=1$). Then,

We will bound the $\mbox{}{\bf{E}}\left[(x^{T}t^{T}ty)^{2}\right]$ term directly:

Notice that if any of the four indices $i_{1},i_{2},j_{1},j_{2}$ appears only once, then the expectation $\mbox{}{\bf{E}}\left[t_{i_{1}}t_{i_{2}}t_{j_{1}}t_{j_{2}}\right]$ corresponding to those indices equals zero. This expectation is non-zero if the four indices are paired in couples or if all four are equal. That is, non-zero expectation happens if

In the above we used $(x^{T}y)^{2}\leq\mbox{}\left\|x\right\|_{2}^{2}\mbox{}\left\|y\right\|_{2}^{2}$. Since we assumed that $\mbox{}\left\|x\right\|_{\infty}\leq\alpha$, the second term on the right hand side of eqn. (49) is bounded by $\frac{\alpha^{2}}{kq}\mbox{}\left\|y\right\|_{2}^{2}$ and the lemma follows since we have assumed that $q\geq\alpha^{2}$.