Fast Regression with an $\ell_\infty$ Guarantee

2017·Arxiv

Abstract

Abstract

Sketching has emerged as a powerful technique for speeding up problems in numerical linear algebra, such as regression. In the overconstrained regression problem, one is given an matrix , as well as an , and one wants to find a vector minimize the residual error . Using the sketch and solve paradigm, one first computes for a randomly chosen matrix S, then outputs so as to minimize

The sketch-and-solve paradigm gives a bound on is well-conditioned. Our main result is that, when S is the subsampled randomized Fourier/Hadamard transform, the error behaves as if it lies in a “random” direction within this bound: for any fixed direction , we have with probability that

where are arbitrary constants. This implies is a factor smaller than . It also gives a better bound on the generalization of to new examples: if rows of A correspond to examples and columns to features, then our result gives a better bound for the error introduced by sketch-and-solve when classifying fresh examples. We show that not all oblivious subspace embeddings S satisfy these properties. In particular, we give counterexamples showing that matrices based on Count-Sketch or leverage score sampling do not satisfy these properties.

We also provide lower bounds, both on how small can be, and for our new guarantee (1), showing that the subsampled randomized Fourier/Hadamard transform is nearly optimal. Our lower bound on shows that there is an separation in the dimension of the optimal oblivious subspace embedding required for outputting an , compared to the dimension of the optimal oblivious subspace embedding required for outputting an the former problem requires dimension while the latter problem can be solved with dimension . This explains the reason known upper bounds on the dimensions of these two variants of regression have differed in prior work.

1 Introduction

Oblivious subspace embeddings (OSEs) were introduced by Sarlos [Sar06] to solve linear algebra problems more quickly than traditional methods. An OSE is a distribution of matrices with such that, for any d-dimensional subspace , with “high” probability S preserves the norm of every vector in the subspace. OSEs are a generalization of the classic JohnsonLindenstrauss lemma from vectors to subspaces. Formally, we require that with probability

Because A is a “tall” matrix with more rows than columns, the system is overdetermined and there is likely no solution to Ax = b, but regression will find the closest point to b in the space spanned by A. The classic answer to regression is to use the Moore-Penrose pseudoinverse:

is the “pseudoinverse” of A (assuming A has full column rank, which we will typically do for simplicity). This classic solution takes time, where is the matrix multiplication constant [CW90, Wil12, Gal14]: time to compute time to compute the inverse. OSEs speed up the process by replacing (2) with

for an OSE S on d + 1-dimensional spaces. This replaces the regression problem with an problem, which can be solved more quickly since lies in the d + 1-dimensional space spanned by b and the columns of A, with high probability S preserves the norm of

That is, S produces a solution which preserves the cost of the regression problem. The running time for this method depends on (1) the reduced dimension m and (2) the time it takes to multiply S by A. We can compute these for “standard” OSE types:

• If S has i.i.d. Gaussian entries, then is sufficient (and in fact, is required [NN14]). However, computing time, which is worse than solving the original regression problem (one can speed this up using fast matrix multiplication, though it is still worse than solving the original problem).

• If S is a subsampled randomized Hadamard transform (SRHT) matrix with random sign flips (see Theorem 2.4 in [Woo14] for a survey, and also see [CNW16] which gives a recent improvement) then m increases to we can compute SA using the fast Hadamard transform in O(nd log n) time. This makes the overall regression problem take

• If S is a random sparse matrix with random signs (the “Count-Sketch” matrix), then m = a decreasing function of the sparsity [CW13, MM13, NN13, BDN15, Coh16]. (The definition of a Count-Sketch matrix is, for any and the column sparsity of matrix S is s. Independently in each column s positions are chosen uniformly at random without replacement, and each chosen position is set to with probability with probability 1/2.) Sparse OSEs can benefit from the sparsity of A, allowing for a running time of nnz(A) denotes the number of non-zeros in A.

When n is large, the latter two algorithms are substantially faster than the naïve

1.1 Our Contributions

Despite the success of using subspace embeddings to speed up regression, often what practitioners are interested is not in preserving the cost of the regression problem, but rather in the generalization or prediction error provided by the vector . Ideally, we would like for any future (unseen) example with high probability.

Ultimately one may want to use to do classification, such as regularized least squares classi-fication (RLSC) [RYP03], which has been found in cases to do as well as support vector machines but is much simpler [ZP04]. In this application, given a training set of examples with multiple (non-binary) labels identified with the rows of an , one creates an each column indicating the presence or absence of one of the r possible labels in each example. One then solves the multiple response regression problem , and uses X to classify future examples. A commonly used method is for a future example a, to compute where are the columns of X. One then chooses the label is maximum.

For this to work, we would like the inner products to be close to is the solution to is the solution to -accurate OSE on d + r dimensional spaces [Sar06], which also satisfies so-called approximate matrix multiplication with error , we get that

where is the spectral norm of , which equals the reciprocal of the smallest singular value of A. To obtain a generalization error bound for an unseen example a, one has

which could be tight if given only the guarantee in (3). However, if the difference vector were distributed in a uniformly random direction subject to (3), then one would expect an factor improvement in the bound. This is what our main theorem shows:

Theorem 1 (Main Theorem, informal)and matrix are given. Let be a subsampled randomized Hadamard transform matrix with rows for an arbitrarily small constant , and any fixed

with probability for an arbitrarily large constant C > 0. This implies that

with probability. If n > poly(d), then by first composing S with a Count-Sketch OSE with poly(d) rows, one can achieve the same guarantee.

(Here is a constant going to zero as n increases; see Theorem 10 for a formal statement of Theorem 1.)

Notice that Theorem 1 is considerably stronger than that of (4) provided by existing guarantees. Indeed, in order to achieve the guarantee (6) in Theorem 1, one would need to set existing OSEs, resulting in rows. In contrast, we achieve only rows. We can improve the bound in Theorem 1 to is a matrix of i.i.d. Gaussians; however, as noted, computing is slower in this case.

Note that Theorem 1 also makes no distributional assumptions on the data, and thus the data could be heavy-tailed or even adversarially corrupted. This implies that our bound is still useful when the rows of A are not sampled independently from a distribution with bounded variance.

The ) of Theorem 1 is achieved by applying (5) to the standard basis vectors and applying a union bound. This guarantee often has a more natural interpretation than the guarantee—if we think of the regression as attributing the observable as a sum of various factors, (6) says that the contribution of each factor is estimated well. One may also see our contribution as giving a way for estimating the pseudoinverse we get that in the sense that each entry is within additive

a lot of work on computing entries of inverses of a matrix, see, e.g., [ADL

sparse (e.g. [Lee12]). In such cases, thresholding to the top k entries will yield an guarantee a factorbetter than (3). One could ask if Theorem 1 also holds for sparse OSEs, such as the Count-Sketch. Surprisingly, we show that one cannot achieve the generalization error guarantee in Theorem 1 with high probability, say, , using such embeddings, despite the fact that such embeddings do approximate the cost of the regression problem up to a factor with high probability. This shows that the generalization error guarantee is achieved by some subspace embeddings but not all.

Theorem 2 (Not all subspace embeddings give the guarantee; informal version of Theorem 20). The Count-Sketch matrix with rows and sparsity —which is an OSE with exponentially small failure probability—with constant probability will have a result that does not satisfy the guarantee (6).

We can show that Theorem 1 holds for S based on the Count-Sketch OSE with probability. We can thus compose the Count-Sketch OSE with the SRHT matrix and obtain an time algorithm to compute achieving (6). We can also compute is a matrix of Gaussians, which is more efficient now that STA only has rows; this will reduce the number of rows to

Another common method of dimensionality reduction for linear regression is leverage score sampling [DMIMW12, LMP13, PKB14, CMM15], which subsamples the rows of A by choosing each row with probability proportional to its “leverage scores”. With rows taken, the result will satisfy the ) with probability . However, it does not give a good bound:

Theorem 3 (Leverage score sampling does not give the guarantee; informal version of Theorem 23)

small failure probability—with constant probability will have a result that does not satisfy the guarantee (6).

Finally, we show that the rows that SRHT matrices use is roughly optimal:

Theorem 4 (Lower bounds for guarantees; informal versions of of Theorem 14 and Corollary 18). Any sketching matrix distribution over matrices that satisfies either the guarantee (3) or the

Notice that our result shows the necessity of the separation between the results originally defined in Equation (3) and (4) of Theorem 12 of [Sar06]. If we want to output some vector that , then it is known that is necessary and sufficient. However, if we want to output a vector , then we show that is necessary and sufficient.

1.1.1 Comparison to Gradient Descent

While this work is primarily about sketching methods, one could instead apply iterative methods such as gradient descent, after appropriately preconditioning the matrix, see, e.g., [AMT10, ZF13, CW13]. That is, one can use an OSE with constant to construct a preconditioner for A and then run conjugate gradient using the preconditioner. This gives an overall dependence of

The main drawback of this approach is that one loses the ability to save on storage space or number of passes when A appears in a stream, or to save on communication or rounds when A is distributed. Given increasingly large data sets, such scenarios are now quite common, see, e.g., [CW09] for regression algorithms in the data stream model. In situations where the entries of A appear sequentially, for example, a row at a time, one does not need to store the full A but only the

Also, iterative methods can be less efficient when solving multiple response regression, where one wants to minimize . This is the case when is constant and t is large, which can occur in some applications (though there are also other applications for which is very small). For example, conjugate gradient with a preconditioner will take time while using an OSE directly will take only time (since one effectively replaces n with O (d) after computing ), separating t from d. Multiple response regression, arises, for example, in the RLSC application above.

1.1.2 Proof Techniques

Theorem 1. As noted in Theorem 2, there are some OSEs for which our generalization error bound does not hold. This hints that our analysis is non-standard and cannot use generic properties of OSEs as a black box. Indeed, in our analysis, we have to consider matrix products of the form for our random sketching matrix S and a fixed matrix U, where k is a positive integer. We stress that it is the same matrix S appearing multiple times in this expression, which considerably complicates the analysis, and does not allow us to appeal to standard results on approximate matrix product (see, e.g., [Woo14] for a survey). The key idea is to recursively reduce using a property of S. We use properties that only hold for specifics OSEs S: first, that each column of S is unit vector; and second, that for all pairs , the inner product between is at mostwith probability

Theorems 20 and 23. To show that Count-Sketch does not give the guarantee, we construct a matrix A and vector b as in Figure 1, which has optimal solution with all coordinates . We then show, for our setting of parameters, that there likely exists an index ing the following property: the jth column of S has disjoint support from the kth column of S for all except for a single k > d, for which share exactly one common entry in their support. In such cases we can compute explicitly, getting . By choosing suitable parameters in our construction, this gives that . The lower bound for leverage score sampling follows a similar construction.

Theorem 14 and Corollary 18. The lower bound proof for the guarantee uses Yao’s minimax principle. We are allowed to fix an sketching matrix S and design a distribution over [A b]. We first write the sketching matrix in its singular value decomposition (SVD). We choose the d + 1 columns of the adjoined matrix [A, b] to be random orthonormal vectors. Consider an orthonormal matrix R which contains the columns of V as its first m columns, and is completed on its remaining columns to an arbitrary orthonormal basis. Then . Notice that is equal in distribution to [A, b], since R is fixed and [A, b] is a random matrix with d + 1 orthonormal columns. Therefore, is equal in distribution to corresponds to the first m rows of an uniformly random matrix with orthonormal columns.

A key idea is that if , then by a result of Jiang [Jsubmatrix of a random orthonormal matrix has o(1) total variation distance to a of i.i.d. N(0, 1/n) random variables, and so any events that would have occurred had G and h been independent i.i.d. Gaussians, occur with the same probability for our distribution up to an factor, so we can assume G and h are independent i.i.d. Gaussians in the analysis.

The optimal solution in the sketch space equals , and by using that SA has the form , one can manipulate to be of the form , where the SVD of We can upper bound , since it is just the maximum singular value of a Gaussian matrix, where r is the rank of S, which allows us to lower bound Then, since h is i.i.d. Gaussian, this quantity concentrates to for a vector h of i.i.d. N(0, 1/n) random variables. Finally, we can lower bound by the Pythagorean theorem, and now we have that is the identity, and so this expression is just equal to the rank of , which we prove is at least d. Noting that for our instance, putting these bounds together gives . The last ingredient is a way to ensure that the rank of S is at least d. Here we choose another distribution on inputs A and b for which it is trivial to show the rank of S is at least d with large probability. We require S be good on the mixture. Since S is fixed and good on the mixture, it is good for both distributions individually, which implies we can assume S has rank d in our analysis of the first distribution above.

1.2 Notation

For a positive integer, let [n] = {1, 2, . . . , n}. For a vector . For a matrix to be the spectral norm of to be the Frobenius norm of to denote the Moore-Penrose pseudoinverse of is its SVD (where ), is given by

In addition to notation, for two functions f, g, we use the shorthand indicate that ) for an absolute constant for constants c, C.

Figure 1: Our construction of A and b for the proof that Count-Sketch does not obey the guarantee.

Definition 5 (Subspace Embedding)-subspace embedding for the column space of an is a matrix S for which for all

Definition 6 (Approximate Matrix Product)be a given approximation parameter. Given matrices A and B, where A and B each have n rows, the goal is to output a matrix C so that . Typically C has the form , for a random matrix S with a small number of rows. In particular, this guarantee holds for the subsampled randomized Hadamard transform

2 Warmup: Gaussians OSEs

We first show that if S is a Gaussian random matrix, then it satisfies the generalization guarantee. This follows from the rotational invariance of the Gaussian distribution.

has full column rank. If the entries of , then for any vectors , we have, with probability

Because SA has full column rank with probability . Therefore

Thus it suffices to only consider vectors , or equivalently . In such cases, SU will be independent of Sb, which will give the result. The proof is in Appendix A.

3 SRHT Matrices

We first provide the definition of the subsampled randomized Hadamard transform(SRHT): let diagonal matrix with i.i.d. diagonal entries , for which in uniform on . The matrix is the Hadamard matrix of size , and we assume n is a power of coordinates of an n dimensional vector uniformly at random.

For other subspace embeddings, we no longer have that SU and Sb are independent. To analyze them, we start with a claim that allows us to relate the inverse of a matrix to a power series.

, and define

Suppose SA has linearly independent columns and

Proof.

We then bound the kth term of this sum:

be the subsampled randomized Hadamard transform, and let a be a unit vector. Then with probability

Hence, for , this is at most with probability at least

We defer the proof of this lemma to the next section, and now show how the lemma lets us prove that SRHT matrices satisfy the generalization bound with high probability:

has full column rank with subsampled randomized Hadamard transform with

vectors

with probability

Proof. Define

e.g., Theorem 2.4 of [Woo14] and references therein), so we condition on. Hence for . In particular, and we can apply Claim 8.

we may assume without loss of generality.

with probability. Suppose this event occurs, bounding the k = 0 term of (7). Hence it suffices to show that the terms of (7) are bounded by By approximate matrix product (Definition 6), we also have with probability that

Combining with we have for any k that

Since this decays exponentially in k at a rate of , the sum of all terms greater than k is bounded by the kth term. As long as

On the other hand, by Lemma 9, increasing factor, we have for all k that

with probability at least , as long as . Since the can be expanded as a sum of terms of this form, we get that

with probability at least , as long as for a sufficiently large constant C. Combining with (9), the result holds as long as

for any

Combining Different Matrices. In some cases it can make sense to combine different matrices that satisfy the generalization bound.

be drawn from distributions of matrices that are -approximate OSEs and satisfy the generalization bound (6). Then RS satisfies the generalization bound with a constant factor loss in failure probability and approximation factor.

We defer the details to Appendix B.

4 Proof of Lemma 9

Proof. Each column of the subsampled randomized Hadamard transform has the same distribution as is a random sign. It also has

with probability . See, e.g., [LDFU13]. By expanding the following product into a sum, and rearranging terms, we obtain

Note that is independent of . We observe that in the above expression if , then the sum over these indices equals in this case for all c. Moreover, the sum over all indices conditioned on equal to 0. Indeed, in this case, the expression can be factored into the form random variable

Let W be a matrix with . We need Khintchine’s inequality:

Fact 12 (Khintchine’s Inequality)be i.i.d. sign random variables, and let be real numbers. Then there are constants

We note that Khintchine’s inequality sometimes refers to bounds on the moment of though the above inequality follows readily by applying a Markov bound to the high moments.

We apply Fact 12 to each column of W, so that if -th column, we have by a union bound that with probability simultaneously for all columns i. It follows that with the same probability, We condition on this event in the remainder.

Thus, it remains to bound . By squaring and using that otherwise, we have,

for any Having computed the expectation of , we now would like to show concentration. Consider a specific

Thus, we can apply Khintchine’s inequality recursively over all the from which it follows that with probability , for each such ). We thus have with this probability, that completing the proof.

5 Lower bound for ℓ2 and ℓ∞ guarantee

We prove a lower bound for the guarantee, which immediately implies a lower bound for the guarantee.

Definition 13. Given a matrix We say that an algorithm A(A, b, S) that outputs a vector “succeeds” if the following property holds:

is a distribution over with the property that for any

Proof. The proof uses Yao’s minimax principle. Let D be an arbitrary distribution over then E

averaging argument implies the existence of a fixed matrix

Thus, we must construct a distribution

cannot hold for any which does not satisfy The proof can be split into three parts. First, we prove a useful property. Second, we prove a lower bound for the case . Third, we show why is necessary.

(I) We show that [SA, Sb] are independent Gaussian, if both [A, b] and S are orthonormal matrices. We can rewrite SA in the following sense,

where S is the complement of the orthonormal basis identity matrix, and the left submatrix of . Thus, using [J] as long as (because of the total variation distance between [SA, Sb] and a random Gaussian matrix is small, i.e.,

where each entry of H is i.i.d. Gaussian N(0, 1/n). (II) Here we prove the theorem in the case when (we will prove this is necessary in part III. Writing in its SVD, we have

where . By a similar argument in Equation (12), as long as we have that G

also can be approximated by a Gaussian matrix, where each entry is sampled from i.i.d. N(0, 1/n).

Similarly, also can be approximated by a Gaussian matrix, where each entry is

sampled from i.i.d. N(0, 1/n).

again, we have

where the the first equality follows by Equation (13), the second equality follows by the SVD of G, the third and fifth equality follow by properties of the pseudo-inversehas orthonormal rows and is a diagonal matrix, and the fourth equality follows since is an orthonormal basis.

Because each entry of is sampled from an i.i.d. Gaussian N(0, 1), using the result of [Ver10] we can give an upper bound for the maximum singular value of with probability at least .99. Thus,

Because h is a random Gaussian vector which is independent of

where each entry of h is sampled from i.i.d. Gaussian N(0, 1/n). Then, using the Pythagorean Theorem,

from ; with probability 1/2, A, b are sampled from In distribution is a random

orthonormal basis and d is always orthogonal to A. In distribution identity matrix in the top-d rows and 0s elsewhere, while b is a random unit vector. Then, for any (A, b) sampled from needs to work with probability at least 9/10. Also for any (A, b) sampled from S needs to work with probability at least 9/10. The latter two statements follow since overall S succeeds on with probability at least 19/20. Consider the case where A, b are sampled from distribution

Then consider which is the optimal solution to

where S can be decomposed into two matrices Plugging into the original regression problem, , which is at most is a submatrix of S, the rank of S is also d.

It remains to define several tools which are used in the main proof of the lower bound.

Claim 15. For any matrix , if each entry of a vector is chosen from an i.i.d Gaussian

Proof.

Let random variables. The random variables t degree of freedom. Furthermore, the following tail bounds are known.

Fact 16 (Lemma 1 of [LM00])random variables. Then for any

and

Definition 17. Given a matrix We say that an algorithm B(A, b, S) that outputs a vector “succeeds” if the following property holds:

is a distribution over with the property that for any

References

[ADLPatrick Amestoy, Iain S. Duff, Jean-Yves L’Excellent, Yves Robert, François-Henry Rouet, and Bora Uçar. On computing inverse entries of a sparse matrix in an out-of-core environment. SIAM J. Scientific Computing, 34(4), 2012.

[AMT10] Haim Avron, Petar Maymounkov, and Sivan Toledo. Blendenpik: Supercharging lapack’s least-squares solver. SIAM J. Scientific Computing, 32(3):1217–1236, 2010.

[BDN15] Jean Bourgain, Sjoerd Dirksen, and Jelani Nelson. Toward a unified theory of sparse dimensionality reduction in euclidean space. In STOC, 2015.

[CMM15] Michael B Cohen, Cameron Musco, and Christopher Musco. Ridge leverage scores for low-rank approximation. arXiv preprint arXiv:1511.07263, 2015.

[CNW16] Michael B Cohen, Jelani Nelson, and David P. Woodruff. Optimal approximate matrix product in terms of stable rank. In Proceedings of the 43rd International Colloquium on Automata, Languages and Programming (ICALP 2016), Rome, Italy, July 12-15, arXiv preprint arXiv:1507.02268, 2016.

[Coh16] Michael B. Cohen. Nearly tight oblivious subspace embeddings by trace inequalities. In Proceedings of the Twenty-Seventh Annual ACM-SIAM Symposium on Discrete Algorithms, SODA 2016, Arlington, VA, USA, January 10-12, 2016, pages 278–287, 2016.

[CW90] Don Coppersmith and Shmuel Winograd. Matrix multiplication via arithmetic progressions. J. Symb. Comput., 9(3):251–280, 1990.

[CW09] Kenneth L. Clarkson and David P. Woodruff. Numerical linear algebra in the streaming model. In Proceedings of the forty-first annual ACM symposium on Theory of computing, pages 205–214. ACM, 2009.

[CW13] Kenneth L Clarkson and David P. Woodruff. Low rank approximation and regression in input sparsity time. In Proceedings of the forty-fifth annual ACM symposium on Theory of computing, pages 81–90. ACM, 2013.

[DMIMW12] Petros Drineas, Malik Magdon-Ismail, Michael W. Mahoney, and David P. Woodruff. Fast approximation of matrix coherence and statistical leverage. Journal of Machine Learning Research, 13:3475–3506, 2012.

[DMMS11] Petros Drineas, Michael W Mahoney, S Muthukrishnan, and Tamás Sarlós. Faster least squares approximation. Numerische Mathematik, 117(2):219–249, 2011.

[Gal14] François Le Gall. Powers of tensors and fast matrix multiplication. In International Symposium on Symbolic and Algebraic Computation, ISSAC ’14, Kobe, Japan, July 23-25, 2014, pages 296–303, 2014.

[JTiefeng Jiang et al. How many entries of a typical orthogonal matrix can be approximated by independent normals? The Annals of Probability, 34(4):1497–1529, 2006.

[LAKD08] Song Li, Shaikh S. Ahmed, Gerhard Klimeck, and Eric Darve. Computing entries of the inverse of a sparse matrix using the FIND algorithm. J. Comput. Physics, 227(22):9408–9427, 2008.

[LDFU13] Yichao Lu, Paramveer Dhillon, Dean Foster, and Lyle Ungar. Faster ridge regression via the subsampled randomized hadamard transform. In Proceedings of the Neural Information Processing Systems (NIPS) Conference, 2013.

[Lee12] Jeff Leek. Prediction: the lasso vs. just using the top 10 predictors. http://simplystatistics.org/2012/02/23/ prediction-the-lasso-vs-just-using-the-top-10/, 2012.

[LM00] Beatrice Laurent and Pascal Massart. Adaptive estimation of a quadratic functional by model selection. Annals of Statistics, pages 1302–1338, 2000.

[LMP13] Mu Li, Gary L Miller, and Richard Peng. Iterative row sampling. In Foundations of Computer Science (FOCS), 2013 IEEE 54th Annual Symposium on, pages 127–136. IEEE, 2013.

[MM13] Xiangrui Meng and Michael W Mahoney. Low-distortion subspace embeddings in input-sparsity time and applications to robust linear regression. In Proceedings of the forty-fifth annual ACM symposium on Theory of computing, pages 91–100. ACM, 2013.

[NN13] Jelani Nelson and Huy L Nguyên. Osnap: Faster numerical linear algebra algorithms via sparser subspace embeddings. In Foundations of Computer Science (FOCS), 2013 IEEE 54th Annual Symposium on, pages 117–126. IEEE, 2013.

[NN14] Jelani Nelson and Huy L. Nguyên. Lower bounds for oblivious subspace embeddings. In Automata, Languages, and Programming - 41st International Colloquium, ICALP 2014, Copenhagen, Denmark, July 8-11, 2014, Proceedings, Part I, pages 883–894, 2014.

[PKB14] Dimitris Papailiopoulos, Anastasios Kyrillidis, and Christos Boutsidis. Provable deterministic leverage score sampling. In Proceedings of the 20th ACM SIGKDD international conference on Knowledge discovery and data mining, pages 997–1006. ACM, 2014.

[RYP03] Ryan Rifkin, Gene Yeo, and Tomaso Poggio. Regularized least-squares classification. Nato Science Series Sub Series III Computer and Systems Sciences, 190:131–154, 2003.

[Sar06] Tamás Sarlós. Improved approximation algorithms for large matrices via random projections. In Symposium on, pages 143–152. IEEE, 2006.

[Ver10] Roman Vershynin. Introduction to the non-asymptotic analysis of random matrices. arXiv preprint arXiv:1011.3027, 2010.

[Wil12] Virginia Vassilevska Williams. Multiplying matrices faster than coppersmithwinograd. In Proceedings of the 44th Symposium on Theory of Computing Conference, STOC 2012, New York, NY, USA, May 19 - 22, 2012, pages 887–898, 2012.

[Woo14] David P. Woodruff. Sketching as a tool for numerical linear algebra. Foundations and Trends in Theoretical Computer Science, 10(1-2):1–157, 2014.

[ZF13] Anastasios Zouzias and Nikolaos M. Freris. Randomized extended kaczmarz for solving least squares. SIAM J. Matrix Analysis Applications, 34(2):773–793, 2013.

[ZP04] Peng Zhang and Jing Peng. Svm vs regularized least squares classification. In ICPR (1), pages 176–179, 2004.

Appendix A Proof for Gaussian case

Lemma 19. If the entries of

for any vectors a, b with probability

Proof. With probability 1, the matrix SA has linearly independent columns, and so

It is well-known (stated, for example, explicitly in Theorem 2.3 of [Woo14]) that with probability , the singular values of . We condition on this event. It follows that

where the first equality uses that V is a rotation, the first inequality follows by sub-multiplicativity, and the second inequality uses that the singular values of SU are in the range with probability

The main observation is that since is statistically independent from Sb. Hence, Sb is distributed as , conditioned on the vector . It follows that conditioned on the value of is distributed as

and so using (14) , with probability

B Combining Diﬀerent Matrices

In some cases it can make sense to combine different matrices that satisfy the generalization bound.

be drawn from distributions of matrices that are -approximate OSEs and satisfy the generalization bound (6). Then RS satisfies the generalization bound with a constant factor loss in failure probability and approximation factor.

Proof. For any vectors we want to show

As before, it suffices to consider the case. We have with probability

suppose this happens. We also have by the properties of R, applied to SA and Sb, that

Because S is an OSE, we have . Therefore

We describe a few of the applications of combining sketches.

B.1 Removing dependence on n via Count-Sketch

One of the limitations of the previous section is that the choice of k depends on n. To prove that theorem, we have to assume that log d > log log n. Here, we show an approach to remove that assumption.

The main idea is instead of applying matrix directly, we pick two matrices is FastJL matrix and C is Count-Sketch matrix with s = 1. We first compute , then compute . The benefit of these operations is S only needs to multiply with a matrix (CA) that has poly(d) rows, thus the assumption we need is log d > log log(poly(d)) which is always true. The reason for choosing C as a Count-Sketch matrix with (2) The running time is

B.2 Combining Gaussians and SRHT

By combining Gaussians with SRHT matrices, we can embed into the optimal dimension with fast embedding time.

B.3 Combining all three

By taking Gaussians times SRHT times Count-Sketch, we can embed into the optimal dimension embedding time.

C Count-Sketch does not obey the ℓ∞ guarantee

Here we demonstrate an A and a b such that Count-Sketch will not satisfy the guarantee with constant probability, so such matrices cannot satisfy the generalization guarantee (6) with high probability.

be drawn as a Count-Sketch matrix with s nonzeros per column. There exists a matrix such that, if , then the “true” solution and the approximation have large distance with constant probability:

Plugging in we find that

even though such a matrix is an OSE with probability exponential in s. Therefore there exists a constant c for which this matrix does not satisfy the generalization guarantee (6) with probability.

Proof. We choose the matrix A to be the identity on its top

set the value of the first d coordinates of vector and set the value to be for the next coordinates, with the remaining entries all zero. Note that and

Let denote the kth column vector of matrix . We define two events, Event I, ; Event II, such that , and all other . Using Claim 21, with probability at least .99 there exists a k for which both events hold.

Given the constructions of A and b described early, it is obvious that

Conditioned on event I and II are holding, then denote . Consider the terms involving in the quadratic form

it can be written as . Hence the optimal , which is different from the desired . Plugging in our requirement of

where the last inequality follows by . Thus, we get the result.

, with probability at least .99 there exists a for which both event I and II hold.

Figure 2: Count Sketch matrix Event I, for any . Event II, there exists a unique intersect at exactly one location(row index).

, then for any be an indicator that the entries of column i are disjoint from all , so by Markov’s inequality, with probability .99, we have .99d columns having this property (indeed, the expected value of is at most ). Define Event E to be that .99d columns of first d columns have the property that the entries of that column are disjoint from all the other columns. Let S be the set of these .99d columns. Let N be the union of supports of columns in S.

Each column non-zero entries. Define event F( which is similar as event E) to be that columns of the next columns have the property that the entries of that column are disjoint from all the other columns. By the same argument, since with probability .99, we have columns in being disjoint from other columns in . Condition on event F holding. Let L be the multiset union of supports of all columns in be the union of supports of all columns in , that is, the set union rather than the multiset union. Note that because of columns are disjoint from each other.

The intersection size x of N and M is hyper-geometrically distributed with expectation

By a lower tail bound for the hypergeometric distribution

where

where the last inequality follows by setting . Thus, we get with probability .99, the intersection size is at least

Now let W be the distinct elements in L\M, so necessarily . By an upper tail bound for the hypergeometric distribution, the intersection size y of W and N satisfies

where , we again get

If |W| = 0, then y = 0. Otherwise, we can set so that this probability is less than .01, and we get with probability .99, the intersection size y is at most Note that we have that are bounded by suffices to ensure y is at most , and earlier that x is at least

The probability one of the |S| blocks in N has two or more intersections with M is less than

times the probability two random distinct items in the intersection land in the block. This

probability is

So the expected number of such blocks is which is less than , which we have. So, there are at least x/2 blocks which have intersection size exactly 1 with N. Note that the number of intersections of the |S| blocks with W is at most y, which is at most , and therefore there exists a block, that is, a column among the first d columns, which intersects M in exactly one position and does not intersect W. This is our desired column. Thus, we complete the proof.

D Leverage score sampling does not obey the ℓ∞ guarantee

Not only does Count-Sketch fail, but so does leverage score sampling, which is a technique that takes a subsample of rows of A with rescaling. In this section we show an A and a b such that leverage score sampling will not satisfy the guarantee. We start with a formal definition of leverage scores.

Definition 22 (Leverage Scores). Given an arbitrary denote the matrix consisting of the d left singular vectors of denote the th row of the matrix is a row vector. Then the leverage scores of the rows of A are given by

The leverage score sampling matrix can be thought of as a square diagonal matrix with diagonal entries chosen from some distribution. If , it means we do not choose the i-th row of matrix , it means we choose that row of the matrix A and also rescale that row. We show that the leverage score sampling matrix cannot achieve guarantee, nor can it achieve our notion of generalization error.

be a leverage score sampling matrix with m nonzeros on the diagonal. There exists a matrix and a vector such that, if , then the “true” solution and the approximation have large distance with constant probability:

Therefore there exists a constant c for which this matrix does not satisfy the generalization guarantee (6) with probability.

Proof. We choose the matrix A to be the identity on its top d rows, and L scaled identity matrices

for the next dL rows, where (to normalize each column of A), which implies . Choose some . Set the value of the first d coordinates of vector b to be and set the value to be for the next coordinates, with the remaining entries all zero. Note that

the following term,

where the optimal should be set to 1/d. The other involves the term:

where the optimal should be set to . Because we are able to choose

Second, we compute . With high probability, there exists a j satisfying Equation (16), but after applying leverage score sampling, the middle term of Equation (16) is removed. Let denote the leverage score of each of the top denote the leverage score of each of the next Ld rows of A. We need to discuss the cases separately. If m > d, then the following term involves

to be

Third, we need to compute . It is easy to see that is an orthonormal matrix. The upper bound for , and the lower bound is also a constant, which can be proved in the following way:

E Bounding E[∥Z∥2F]

Before getting into the proof details, we define the key property of S being used in the rest of the proofs.

Definition 24 (All Inner Product Small(AIPS) Property)

we say that S satisfies the “AIPS” property.

is a subsampled Hadamard transform matrix, then the AIPS property holds with probability at least

Proof. From the structure of , we have with probability . Applying a union bound over pairs, we obtain that

The main idea for bounding is to rewrite it as AIPS holds AIPS does not hold ]. Because Pr[ AIPS does not hold] is at most 1/poly(n), the first term dominates the second term, which means we only need to pay attention to the first term. We repeatedly apply this idea until all the S are removed.

We start by boundinbg by squaring and using that otherwise. Then, we obtain,

where the first equality is from our conditioning, and the second equality is the definition of Hence,

where the first equality follows from the property just mentioned, and the inequality follows by including back the tuples of indices for which , using that each summand is non-negative.

Our next step will be to bound the term

where the first equality is the definition of B, the second equality follows by separating out the index , the third equality uses that , that is, the squared norm of the , the inequality follows since all rows of a projection matrix norm at most 1, and the final equality uses that no longer appears in the expression.

We now merge our bounds for the terms A and B in the following way:

where the first two inequalities and first equality are by definition of A and B above. The first inequality follows by induction, since at this point we have replaced , and can repeat the argument, incurring another multiplicative factor of O(d(log n)/r)+1. Repeating the induction in this way we arrive at the last inequality. Finally, the last inequality follows by plugging in the definition of , using that

where the inequality follows since each row of has norm at most 1, and a is a unit vector. The final result is that

designed for accessibility and to further open science