Faster least squares approximation

1. Faster least squares approximation Tamas Sarlos Yahoo! Research IPAM 23 Oct 2007 http://www.ilab.sztaki.hu/~stamas Joint work with P. Drineas, M. Mahoney, and S. Muthukrishnan.

2. Outline Randomized methods for • Least squares approximation • Lp regression • Regularized regression • Low rank matrix approximation

3. Models, curve fitting, and data analysis In MANY applications (in statistical data analysis and scientific computation), one has n observations (values of a dependent variable y measured at values of an independent variable t): Model y(t) by a linear combination of d basis functions: A is an n x d “design matrix” with elements: In matrix-vector notation:

4. Least Squares (LS) Approximation We are interested in over-constrained L2 regression problems, n >> d. Typically, there is no x such that Ax = b. Want to find the “best” x such that Ax ≈ b. Ubiquitous in applications & central to theory: Statistical interpretation: best linear unbiased estimator. Geometric interpretation: orthogonally project b onto span(A).

5. Many applications of this! • Astronomy: Predicting the orbit of the asteroid Ceres (in 1801!). Gauss (1809) -- see also Legendre (1805) and Adrain (1808). First application of “least squares optimization” and runs in O(nd2) time! • Bioinformatics: Dimension reduction for classification of gene expression microarray data. • Medicine: Inverse treatment planning and fast intensity-modulated radiation therapy. • Engineering: Finite elements methods for solving Poisson, etc. equation. • Control theory: Optimal design and control theory problems. • Economics: Restricted maximum-likelihood estimation in econometrics. • Image Analysis: Array signal and image processing. • Computer Science: Computer vision, document and information retrieval. • Internet Analysis: Filtering and de-noising of noisy internet data. • Data analysis: Fit parameters of a biological, chemical, economic, social, internet, etc. model to experimental data.

6. Projection of b on the subspace spanned by the columns of A Pseudoinverse of A Exact solution to LS Approximation Cholesky Decomposition: If A is full rank and well-conditioned, decompose ATA = RTR, where R is upper triangular, and solve the normal equations: RTRx=ATb. QR Decomposition: Slower but numerically stable, esp. if A is rank-deficient. Write A=QR, and solve Rx = QTb. Singular Value Decomposition: Most expensive, but best if A is very ill-conditioned. Write A=UVT, in which case: xOPT = A+b = V-1kUTb. Complexity is O(nd2) for all of these, but constant factors differ.

7. Reduction of LS Approximation to fast matrix multiplication Theorem: For any n-by-d matrix A, and d dimensional vector b, where n>d, we can compute x*=A+b in O(n/dM(d)) time. Here M(d) denotes the time needed to multiply two d-by-d matrices. • Current best M(d)=d, =2.376… Coppersmith–Winograd (1990) Proof: Ibarra et al. (1982) generalize the LUP decomposition for rectangular matrixes Direct proof: • Recall the normal equations ATAx = ATb • Group the rows of A into blocks of d and compute ATA as the sum of n/d square products, each d-by-d • Solve the normal equations • LS in O(nd1.376) time, but impractical. We assume O(nd2).

8. Original (expensive) sampling algorithm Drineas, Mahoney, and Muthukrishnan (SODA, 2006) • Algorithm • Fix a set of probabilities pi, i=1,…,n. • Pick the i-th row of A and the i-th element of b with probability • min {1, rpi}, • and rescale by (1/min{1,rpi})1/2. • Solve the induced problem. Theorem: Let: If the pi satisfy: for some  (0,1], then w.p. ≥ 1-, • These probabilities pi’s are statistical leverage scores! • “Expensive” to compute, O(nd2) time, these pi’s!

9. Fast LS via uniform sampling Drineas, Mahoney, Muthukrishnan, and Sarlos (2007) • Algorithm • Pre-process A and b with a “randomized Hadamard transform”. • Uniformly sample • r=O(d log(n) log(d log(n))/) • constraints. • Solve the induced problem: • Non-uniform sampling will work, but it is “slow.” • Uniform sampling is “fast,” but will it work?

10. A structural lemma Approximate by solving where is any matrix. Let be the matrix of left singular vectors of A. assume that -fraction of mass of b lies in span(A). Lemma: Assume that:  Then, we get relative-error approximation:

11. Randomized Hadamard preprocessing Facts implicit or explicit in: Ailon & Chazelle (2006), or Ailon and Liberty (2008). Let Hn be an n-by-n deterministic Hadamard matrix, and Let Dn be an n-by-n random diagonal matrix with +1/-1 chosen u.a.r. on the diagonal. Fact 1: Multiplication by HnDn doesn’t change the solution: (since Hn and Dn are orthogonal matrices). Fact 2: Multiplication by HnDn is fast - only O(n log(r)) time, where r is the number of elements of the output vector we need to “touch”. Fact 3: Multiplication by HnDn approximately uniformizes all leverage scores:

12. A matrix perturbation theorem

13. Establishing structural conditions Lemma: |1-i2(SHUA)| ≤ 1/10, with constant probability, i.e,. singular values of SHUA are close to 1. Proof: Apply spectral norm matrix perturbation bound. Lemma: ||(SHUA)TSHb|| ≤  OPT, with constant probability, i.e., SHUA is nearly orthogonal to b. Proof: Apply Forbenius norm matrix perturbation bound. Bottom line: • Fast algorithm approximately uniformizes the leverage of each data point. • Uniform sampling is optimal up to O(1/log(n)) factor, so over-sample a bit. • Overall running time O(ndlog(…)+d3log(…)/)

14. Fast LS via sparse projection Drineas, Mahoney, Muthukrishnan, and Sarlos (2007) - sparse projection matrix from Matousek’s version of Ailon-Chazelle 2006 • Algorithm • Pre-process A and b with a randomized Hadamard transform. • Multiply preprocessed input by sparse random k x n matrix T, where • and where k=O(d/) and q=O(d log2(n)/n+d2log(n)/n) . • Solve the induced problem: • Dense projections will work, but it is “slow.” • Sparse projection is “fast,” but will it work? • -> YES! Sparsity parameter q of T related to non-uniformity of leverage scores!

15. Lp regression problems Dasgupta, Drineas, Harb, Kumar, and Mahoney (2008) We are interested in over-constrained Lp regression problems, n >> d. Lp regression problems are convex programs (or better!). There exist poly-time algorithms. We want to solve them faster!

16. Lp norms and their unit balls Recall the Lp norm for : Some inequality relationships include:

17. Solution to Lp regression Lp regression can be cast as a convex program for all . For p=1, Sum of absolute residuals approximation (minimize ||Ax-b||1): For p=∞, Chebyshev or mini-max approximation (minimize ||Ax-b||∞): For p=2, Least-squares approximation (minimize ||Ax-b||2):

18. What made the L2 result work? • The L2 sampling algorithm worked because: • For p=2, an orthogonal basis (from SVD, QR, etc.) is a “good” or “well-conditioned” basis. • (This came for free, since orthogonal bases are the obvious choice.) • Sampling w.r.t. the “good” basis allowed us to perform “subspace-preserving sampling.” • (This allowed us to preserve the rank of the matrix and distances in its span.) • Can we generalize these two ideas to p2?

19. Sufficient condition Easy Approximate minx||b-Ax||p by solving minx||S(b-Ax)||p giving x’ Assume that the matrix S is an isometry over the span of A and b: for any vector y in span(Ab) we have (1-)||y||p ||Sy||p (1+ ) ||y|| p Then, we get relative-error approximation: ||b-Ax||p  (1+ ) minx|| S(b-Ax)||p (This requires O(d/2) for L2 regression.) How to construct an Lp isometry?

20. p-well-conditioned basis (definition) Let A be an n x m matrix of rank d<<n, let p  [1,∞), and q its dual. Definition: An n x d matrix U is an (,,p)-well-conditioned basis for span(A) if: (1) |||U|||p≤, (where |||U|||p = (ij|Uij|p)1/p ) (2) for all z in Rd, ||z||q ≤  ||Uz||p. U is a p-well-conditioned basis if ,=dO(1), independent of m,n.

21. p-well-conditioned basis (existence) Let A be an n x m matrix of rank d<<n, let p  [1,∞), and q its dual. Theorem: There exists an (,,p)-well-conditioned basis U for span(A) s.t.: if p < 2, then  = d1/p+1/2 and  = 1, if p = 2, then  = d1/2and  = 1, if p > 2, then  = d1/p+1/2and  = d1/q-1/2. U can be computed in O(nmd+nd5log n) time (or just O(nmd) if p = 2).

22. p-well-conditioned basis (construction) • Algorithm: • Let A=QR be any QR decomposition of A. • (Stop if p=2.) • Define the norm on Rd by ||z||Q,p ||Qz||p. • Let C be the unit ball of the norm ||•||Q,p. • Let the d x d matrix F define the Lowner-John ellipsoid of C. • Decompose F=GTG, • where G is full rank and upper triangular. • Return U = QG-1 • as the p-well-conditioned basis.

23. Subspace-preserving sampling Let A be an n x m matrix of rank d<<n, let p  [1,∞). Let U be an (,,p)-well-conditioned basis for span(A), Theorem: Randomly sample rows of A according to the probability distribution: where: Then, with probability 1- , the following holds for all x in Rm:

24. Regularized regression • If you don’t trust your data (entirely) • Minimize ||b-Ax||p+||x||r, where r arbitrary • Lasso is especially popular: minx||b-Ax||2+||x||1 • Sampled Lp regression works for the regularized version too • We want more: imagine 

25. Ridge regression • Tikhonov regularization, back to L2 • Minimize ||b-Ax||22+2||x||22 • Equivalent to augmenting n-by-d A with dand b with 0 • SVD of A UDiag(sqrt(i))VT • Let D = i sqrt(i) << d • For (1+)relative error: • precondition with randomized FFT • sample r=O(D log(n) log(D log(n))/) rows uniformly.

26. SVD of a matrix Any m x n matrix A can be decomposed as: U (V): orthogonal matrix containing the left (right) singular vectors of A. S: diagonal matrix containing the singular values of A, ordered non-increasingly. : rank of A, the number of non-zero singular values. Exact computation of the SVD takes O(min{mn2 , m2n}) time.

27. SVD and low-rank approximations Truncate the SVD by keeping k ≤  terms: • Uk (Vk): orthogonal matrix containing the top k left (right) singular vectors of A. • k: diagonal matrix containing the top k singular values of A. • Ak is the “best” matrix among all rank-k matrices wrt. to the spectral and Frobenius norms • This optimality property of very useful in, e.g., Principal Component Analysis (PCA), LSI, etc.

28. Approximate SVD • Starting with the seminal work of Frieze et al (1998) • and Papadimitriou et al (1998) till 2005 several results • of the form • ||A-Ak’||F ≤ ||A-Ak||F + t||A||F • ||A||F might be a significantly larger than ||A-Ak||F • Four independent results [H-P, DV, MRT, S 06] on • ||A-Ak’||F ≤ (1+)||A-Ak||F • This year Shyanalkumar and Deshpande …

29. Relative error SVD via random projections • Theorem Let S be n by r=O(k/) random matrix with • 1 entries. Project A to the subspace spanned by the rows of SA and then compute Ak’ as • the rank-k SVD of the projection. Then • ||A-Ak’||F ≤ (1+)||A-Ak||F with prob. at least 1/3 • By keeping an orthonormal basis of SA total time is O(Mr+(n+m)r^2) with M non-zeroes • We need to make only 2 passes over the data • Can use fast random projections and avoid the SA projection, spectral norm (talks by R and T)

30. SVD proof sketch Can show: ||A-Ak’||F2 ≤ ||A-Ak||F2 + ||Ak-Ak(SAk)+(SA)||F2 For all j columns of A consider the regression A(j)≈Akx Exact: Ak(j), error is ||A(j)-Ak(j)||2 Approximate: Ak-Ak(SAk)+SA(j)

31. Conclusion Fast methods based on randomized preprocessing and sampling or sparse projections for • Least squares approximation • Regularized regression • Low rank matrix approximation

32. Thank you!