Notes for Randomized Numerical Linear Algebra, updated 03 Mar 2015 Starting today, we will videocapture each lecture, broadcast them live, and archive them on youtube. The link today will be https://www.youtube.com/watch?v=q6mcfCCXjUo We will post instructions for future lectures on the website. In addition to material on the class website https://who.rocq.inria.fr/Laura.Grigori/RandNLA.html we borrow from material presented in Ma221 in Fall 2014, at http://www.cs.berkeley.edu/~demmel/ma221_Fall14/Lectures/Lecture_19.txt Papers on recent communication avoiding algorithms like TSQR may be found at bebop.cs.berkeley.edu Our goals in this introductory lecture are (1) to discuss the metrics for performance, accuracy and reliability that we want to use to compare randomized and nonrandomized algorithms (both direct and iterative), which we hope will inform discussions later in the semester too; (2) recall different kinds of Least Squares (LS) problems, and for the simplest case of full-rank, overdetermined problems, the usual nonrandomized algorithms (direct: NE = Normal Equations, QR and SVD, and iterative: LSQR) and what their metrics are; (3) discuss randomized LS solvers, including Blendenpik, Clarkson-Woodruff, ... Metrics: (1) #Flops: This may just depend on the matrix dimensions, or also sparsity structure. For dense, square n x n matrices, direct linear algebra algorithms typically cost O(n^3). Given a fast matmul algorithm costing O(n^w), there are "numerically stable" O(n^(w+ eps)) algorithms, for any eps>0; see "Fast Linear Algebra is Stable" by Demmel/Dumitriu/Holtz for details. In at least one case, the nonsymmetric eigenproblem, the only algorithm attaining this uses randomization. For m x n matrices with m>n, for example the overdetermined LS problem, the complexities are O(m*n^2) and O(m*n^(w - 1 + eps)). We will mostly consider O(n^3) or O(m*n^2) algorithms, which are often (but not always) more practical. For sparse matrices, with nnz(A) = number of nonzeros in A << m*n, we seek both direct and iterative algorithms that do much less arithmetic. Suppose A has n(i) nonzeros in row i, so nnz(A) = sum_i n(i). Looking ahead, we will have nonrandomized direct algorithms for LS with m >> n that cost O(sum_i n(i)^2) flops (plus ``lower order terms''), which will be a lower bound under certain assumptions, in contrast to a randomized algorithm that costs O(sum_i n(i)) = O(nnz(A)) flops (Clarkson-Woodruff), or O(m*n) in the dense case. Iterative methods like LSQR are also used, where the inner loop is typically matrix-vector multiplication and so does 2*nnz(A) flops. The number of iterations depends on numerical properties of the matrix, like the distribution of its singular values, often summarized by the condition number kappa(A) = sigma_max(A)/sigma_min(A). Looking ahead, we will describe a nonrandomized iterative method that does O(kappa(A)) iterations, and a randomized one that does O(1) iteration w.h.p. (Blendenpik), and so again costs O(m*n) in the dense case. (2) Communication cost: Communication means moving data, either between levels of memory hierarchy or between parallel processors over a network, and costs orders of magnitude more per operation than a flop. #Flops is a good metric for small problems that fit in cache, but the motivation for randomized algorithms are very large problems where this is not the case. So our performance model for sending one "message" of n contiguous words will be alpha + n*beta = latency + n/bandwidth where a "message" could be a cache line in a shared memory system, an actual message in a distributed memory system (like Hopper and Edison), or a disk access. Our overall performance model (along critical path of an algorithm) is then gamma*#Flops + beta*#Words_sent + alpha*#Messages = gamma*F + beta*W + alpha*S where gamma << beta << alpha, with the gaps growing over time following technological trends. We can use the same model for energy, where gamma etc have units of joules instead of seconds, since energy usage is becoming at least as important metric as run-time (and they are obviously correlated). There has been a lot of recent work on communication lower bounds, starting with linear algebra, but extending to code that can be described abstractly as “for all k-tuples of integers in some some subset of Z^k, execute some code that accesses arrays whose subscripts are affine functions of those integers.” This class of algorithms clearly includes matrix multiplication C = A*B, with 3-tuples of integers (i,j,k) accessing arrays A(i,k) B(k,j) and C(i,j), but most other linear algebra algorithms (LU, QR, …), for dense or sparse matrices, and either sequential or parallel implementations. And it goes beyond classical linear algebra, including fast linear algebra (eg Strassen's method), some graph algorithms, direct n-body algorithms, tensor contractions, etc. For this class, we will just state the lower bound for (most) classical linear algebra, and use it as a goal against which to measure our algorithms. Since randomized algorithm often use classical linear algebra as building blocks, it will apply at least to components of them. To state the result, consider a 2-level memory hierarchy, called fast and slow, where fast could equal cache and slow = DRAM, or fast = DRAM and slow = disk, or fast = local memory of one parallel processor, and slow = memory of other parallel processors connected by a network. We assume that arithmetic can only be done on data in the fast memory. We want to lower bound W = number of words moved between fast and slow memory = Omega(max( #iterations / M^{1/2}, #inputs, #outputs ) where the number of iterations refers to the inner loop of the algorithm, say a multiply-add for matrix multiply or classical linear algebra more generally. In the case of a overdetermined least squares with a dense m x n input matrix, running sequentially, #iterations = Theta(m n^2), so this bound becomes W = Omega ( max( m*n^2 / M^{1/2} , m*n ) ) = Omega ( m*n ) when n <= M^{1/2} In other words, the lower bound is the size of the input, raising the question of whether we can really just access the data once, even when it is arbitrarily large, provided n <= M^{1/2}, or n^2 <= M, i.e. we can fit a square submatrix in memory. Def: A "streaming algorithm" can solve an arbitrarily large problem that is read once into fast memory (say a row at a time). Note that a streaming algorithm does not even require the data to exist entirely at one time, eg sensors could be producing data in real time that is streamed into a computer and then thrown away. If we instead assume that the data is in fact stored in slow memory, and so can be accessed more than once, we can ask for a slightly weaker goal: Def: A "weakly streaming algorithm" can solve a problem by reading it from slow to fast memory O(1) times. TSQR is a recent deterministic streaming algorithm. It targets dense matrices, but could be modified to take (some) advantage of sparsity. We will later describe randomized (weakly) streaming LS algorithms. Regarding the latency cost, given a lower bound on W, a lower bound on the number of messages is simply given by S >= W/largest_message_size. The largest message size may be determined by characteristics of the HW of SW, but is obviously bounded by M, the fast memory size to which the message is delivered. Thus S = Omega(max( #iterations / M^{3/2}, #inputs/M, #outputs/M ) Lower bounds on both W and S can be attained by deterministic algorithms, TBD for randomized ones. In the parallel case, with p processors, where fast memory is the local memory of one processor and slow memory is the memory of other processors, we may consider a load balanced algorithm so the lower bound becomes W = number of words moved between one processor and other processors = Omega( (#iterations/p) / M^{1/2} ). The smallest we can choose M, for the dense m x n case, is M = m n / p, i.e. there is just enough memory to store the problem, leading to the lower bound W = Omega ( (m n^2 / p)/(m n / p)^{1/2) = Omega ( m^{1/2} n^{3/2} / p^{1/2} ) which is attainable by recent deterministic algorithms. The latency lower bound is again a factor of M smaller: S = Omega ( (n/m)^{1/2} * p^{1/2} ) which is also attainable when m=n. But the lower bound on W remains true if M exceeds its lower bound m*n/p. For example, if we have enough memory for c>1 copies of the data, so M = c*m*n/p, then the lower bound on W decreases by a factor of c^{1/2}, and the lower bound on S decreases by c^{3/2}. This is attainable by matmul, at least when m=n, up to c = p^{1/3}, and for W by QR and LU. On the other hand, we do not expect S to decrease in all cases, because (at least for LU, and we conjecture for QR and SVD, but not NE), there is an independent lower bound W*S = Omega(n^2), so if W decreases, S must increase. We note that for "Strassen-like algorithms" (which requires some technical assumptions), the M^{1/2} in the above lower bounds is replaced by M^{w/2 - 1}, where w is the exponent of matrix multiplication, and is attainable for various nonrandomized algorithms. Finally, for iterative methods, the communication costs are typically dominated by the number of iterations, assuming the matrix is too big to fit in fast memory. To summarize goals (1) (flops) and (2) (communication): (1) Direct (nonrandomized) methods cost O(mn^2) or O(sum_i n(i)^2) flops, where n(i) is the number of nonzeros in A(i,:). Iterative (nonrandomized) methods cost O(nnz(A)) flops per iteration, and may take O(kappa(A)) iterations. Randomized goal: O(nnz(A)) flops, i.e. a constant number of flops per nonzero. (2) Direct: when n^2 <~ M = fast memory size, can access data just once Iterative: when A does not fit in fast memory, access data O(kappa(A)) times. Randomize goal: access data once, or O(1) times (3) Accuracy: For LS, we can ask either about the error in the solution x or the residual r = A*x-b. The usual nonrandomized goal is called backward stability, i.e. getting the exact solution for a slightly different input, since the input usually has uncertainty, so we ask for an answer "as accurate as the data deserves". So this metric of success is that the computed solution xhat satisfies (1) xhat = argmin || (A+E)*xhat - (b+f) || where E and f are "small" compared to A and b, resp. say measured by eps = norm(E)/norm(A) + norm(f)/norm(b) << 1 We may use eps to bound the difference between xhat and the true solution x; to first order in eps: (2) || x-xhat || / || x || <= eps*( kappa(A)/cos(theta) + kappa^2(A)*tan(theta) ) where kappa(A) is the condition number and sin(theta) = norm(r)/norm(b), i.e. theta is the angle between b and it best approximation A*x. So this error bound can be very large in two situations: A is ill-conditioned, i.e. kappa(A) is large, and theta is close to pi/2, i.e. b is nearly orthogonal to the columns of A The residuals r = A*x-b and rhat = (A+E)*xhat - (b+f) can differ by (3) || r - rhat || / || b || <= (1 + 2*kappa(A)) * eps so the residual is less sensitive. Deterministic algorithms based on QR and SVD are backward stable and so attain these bounds with eps = O(machine epsilon) = O(10^(-16)) in double precision. NE is not backward stable, and only guarantees (4) || x - xhat || / || x || <= O(machine epsilon) * kappa^2(A) even when tan(theta) is small. For example, Blendenpik (see below) is "about as stable as QR" in the words of the authors. Finally, in some applications we might only want the weaker criterion that (5) || A*xhat - b || <= (1 + eps) * || A*x-b || where xhat is the computed solution and A*xhat-b its residual, and A*x-b is the true residual, i.e. the fit is about as good in norm. Clarkson-Woodruff (see below) makes such a guarantee for their algorithm: Thm 1 (CW): With probability at least 2/3, norm(A*x-b,2) is at most (1+eps) times larger than the true minimum, with cost O(nnz(A) + "lower order terms") where the "lower order terms" grow like O(n^3*log^6(n/eps)/eps^2). This theorem illustrate a more general feature of randomized algorithms: they depend on a parameter eps, which is the desired accuracy, and the cost (eg number of random samples) grows as epsilon decreases. Furthermore, they have a probability of failure delta (delta = 1/3 above) that may also decrease as epsilon decreases. To summarize, we have presented 5 metrics of accuracy, in order from strongest to weakest ((1) implies (2) and (3), but (2) and (3) could potentially hold without (1)). The cost of a randomized algorithm will generally depend on the desired accuracy epsilon, and possibly the desired probability of failure. (It is also possible to use iterative refinement and extra precision to get even higher precision than (2) and (3), but we will not pursue this here.) Various LS problems: We will only discuss the first one, randomized algorithms for many of the others remain to be explored, afaik: Full rank, overdetermined: argmin || A*x - b || "Damped" overdetermined: argmin || [A; lambda*I] *x - [b;0] || aka ridge regression, or Tikhonov regularization with lambda*B instead Low rank, overdetermined: choose x minimizing both || A*x - b|| and || x || Weighted: argmin || W*(A*x - b) || Constrained: choose x minimizing || A*x - b || such that B*x = c, or || x || <= tau, or || B*x - c || <= tau Underdetermined versions of these Randomized LS solvers: Citations below refer to the following papers: RT = Rokhlin, Tygert MSM = Meng, Saunders, Mahoney DMMS = Drineas, Mahoney, Muthukrishnan, Sarlos DMM = Drineas, Mahoney, Muthukrishnan AMT = Avron, Maymounkov, Toledo CW = Clarkson, Woodruff There are two general approaches: (R1) Let S be a s x m random matrix, let A' = S*A and b' = S*b. Then solve the s x n LS problem argmin || A'*x' - b' ||. When s << m, and S*A is not too expensive, this can be much cheaper than the original problem. This approach includes Random row sampling as presented in DMMS ((1) below) Leverage score sampling as presented in DMMS ((2) below) CW ((7) below) (R2) Again form A'=S*A, and factor A' = Q*R; more generally we only require that the factor Q be well-conditioned, and R both square and that it be easy to solve Rx=y for x (eg R is triangular). Then solve the LS problem iteratively with R as a preconditioner; more precisely solve argmin || A*inv(R)*y - b || iteratively for y, and then let x = inv(R)*y. By construction A*inv(R) will be well-conditioned w.h.p. and so take a bounded number of iterations to converge. These iterative methods are typically mathematically equivalent to a method like CG for square linear systems applied to the normal equations (A^T*A)*x = A^T*b; A^T*A is not formed explicitly, rather we just need to multiply vectors (A^T*A)*y = A^T*(A*y) using two matrix-vector multiplications. This approach includes RT ((4) below) LSRN as presented in MSM ((5) below) Blendenpik of AMT ((6) below) Both approaches depend on the Johnson-Lindenstrauss Lemma (JL) and its successors, which say that projecting n m-vectors into a suitably randomly chosen s-dimensional subspace, i.e. multiplying them by a s x m matrix with orthonormal rows (call it S') will reasonably accurately preserve the distances between all pairs of vectors. So solving the projected least squares problem argmin || S'*(A*x -b) || should give about the same answer. S and S' are not the same; the original S' of JL, a "uniformly distributed" random orthonormal matrix, i.e. with Haar measure, may be computed by doing the QR factorization of an m x s matrix of i.i.d. N(0,1) random variables, and letting S' = Q^T. But computing S' this way is at least as expensive as solving the original LS problem directly (and more if s > n, as it must be). So later approaches to (R1) and (R2) try to use much cheaper S matrices. We now illustrate some of these approaches, from simple to more sophisticated (and without any proofs), and indicate both the papers which they appears and whether the algorithms use approach R1 or R2: (1) (DMMS / R1) Choose s rows of A and b randomly with equal probability, and call the resulting matrix A'. This is equivalent to letting S be an s x m matrix where each row is zero except for a 1 in a random column, and letting A'=S*A. Similarly, b'=S*b. Then if A is "random" and s = O(n log(m) log(n log(m))), then || A*xhat - b || <= (1 + eps) || A*x-b || with "high probability". Of course we can't assume A is random in general. (2) (DMM / R1) Suppose some column of A has just one nonzero entry, say A(i,j). Then unless we are lucky enough to sample row i, A' will have a zero column, and so be singular, making x' ill-defined. So we clearly need to assign a higher sampling probability to row i. Since A(i,j) is the only nonzero in column j, then in A=Q*R, Q(i,j+1:end) =0, since later columns of Q are all orthogonal to column j of A. Thus, recalling that Q is m x n, || Q(i,:) || = 1. This motivates the definition Def: the coherence of A is mu(A) = max_i || Q(i,:) ||^2. The coherence satisfies n/m <= mu(A) <= 1, the upper bound because Q is orthogonal, the lower bound because sum_i || Q(i,:) ||^2 = n. If mu(A) is near 1, i.e. the matrix is coherent, then some rows are more "important" than others, whereas if mu(A) is near n/m, i.e. A is incoherent, then rows are more nearly equally important. || Q(i,:) ||, also called "leverage score", measures how important row i is to sample. Let p_i = || Q(i,:) ||^2 / n be the sampling probability for row i; note that sum_{i=1 to n} p_i = 1. Thm 2: Let A' be obtained by sampling (with replacement) s rows of A where row i is chosen with probability p_i = ||Q(i,:)||^2 /n. This is equivalent to letting S be an s x m matrix with each row is zero except for a 1 in column i, with probability p_i, and letting A'=S*A. Solve argmin ||A'*x' - b'||. Then (under some additional conditions on b), we get || A*x' - b || <= (1+eps) * ||A*x-b|| with probability 1-delta provided we sample at least this many rows: s = Omega( n^2*ln(3/delta) / eps^4 ) Note that the small LS problem, argmin || A'*x' - b' ||, would cost O(sn^2) = O(n^4 * ln(3/delta) / eps^4) if done using a dense deterministic problem, which is potentially much cheaper than a deterministic O(mn^2) algorithm on the original problem, if m >> n. On the other hand, computing the p_i in the first place would cost as much as a QR decomposition, i.e. O(mn^2). (3) (AMT / R2) There is an analogous result to the above theorem, which says that if we use uniform random sampling to get A' = S*A as in (1), then for inv(R') from A' = Q'*R' to be a good preconditioner for LSQR, the number of sampled rows s has to grow at least as fast as m*mu(A), so O(m) for a coherent matrix down to O(n) for an incoherent matrix, i.e. random sampling of a coherent matrix does not work well, as before: Thm 3 (AMT): Let A' be s randomly sampled rows of A, with uniform probability. Let tau = C * sqrt(m * mu(A) * log(s) / s ), and tau/delta < 1. Then with probability at least 1 - delta, kappa(A*inv(R')) <= (1 + tau/delta)/(1 - tau/delta) We will use this to analyze Blendenpik below. (4) (RT / R1 and R2) In the above examples S was as cheap as possible to multiply by, since forming A'=S*A was equivalent to choosing randomly selected rows of A. Alternatively, we can pick a more expensive (denser) S to multiply by, where multiple nonzeros each each row of S "mix" different rows of A, and so avoid the coherence problem. There are many different examples in the literature, and we mention several of them, starting with the first (RT). Here S = G*H where G is s x m, and G = S'*F*D where S' is s x m and samples s rows uniformly at random as in the first example F is an m x m DFT, so cheap to multiply by D is an m x m diagonal matrix where D(i,i) is i.i.d. uniformly on the unit circle in the complex plane, and H is m x m and the product of m x m matrices: H = Theta * Pi * Z * Theta' * Pi' * Z' where Theta and Theta' are each the product of m-1 i.i.d. random Givens rotations [cos(t(i)), sin(t(i)) ; -sin(t(i)), cos(t(i)) ] on the coordinates (1,2), (2,3), ... (m-1,m) Pi and Pi' are independent random permutations Z and Z' are like D above The point is that multiplying a vector of A by S is relatively cheap, O(m log s) because the (sparse) DFT is the dominant cost, and that the patterns of the various factors of G and H "mix" the rows of A and avoid coherence problems. Thm 4 (RT/R1) Choose alpha > 1 and beta > 1, and the number of rows of S equal to s >= (alpha^2+1)/(alpha^2-1) * beta * (n+1)^2 Then the solution x' = argmin || A'*x' - b' || where A'=S*A and b'=S*b satisfies || A*x' - b || <= alpha * || A*x - b || with probability 1 - 1/beta. Note that A' has O(n^2) rows, so solving for x' would cost O(n^4) using a direct method, for a total cost of O(m*n*log(n) + n^4). RT also uses an iterative approach R2, forming A' = Q*R and using inv(R) as a preconditioner: Thm 5 (RT/R2) Choosing s as in Thm 4, guarantees kappa(A*inv(R)) <= 3, so an iterative algorithm converges in a finite number of steps. Unfortunately, Thm 4 says s = Theta(n^2), which is very large, again leading to a total cost of O(m*n*log(n) + n^4), but for a possibly much more accurate answer. RT's numerical experiments show that s = 4*n, which is much cheaper, is enough in practice. (5) (MSM / R2) This paper presents the LSRN algorithm. Here S is chosen with i.i.d. N(0,1) random variables, with s ~ 2*n, and the preconditioner is gotten by taking the SVD of S*A = A' = U * Sigma * V^T = U * R, and using inv(R) as a preconditioner. The paper addresses the rank deficient case too, for simplicity we consider the full rank case: Thm 6 (MSM) The reciprocal singular values of A*inv(R) have the same probability distribution as the singular values of an s x n i.i.d. N(0,1) matrix. This means that the probability of a very large or small singular value falls off exponentially fast, so convergence is fast w.h.p. Since S is a dense matrix, multiplying A' = S*A costs 2*s*nnz(A) = 4*m*n^2 in the case of dense A, so more than the usual direct algorithm in the dense case, but possibly much less in the sparse case. (6) (AMT / R2) This paper presents the Blendenpik algorithm. Here S is chosen very simply as S = F*D, where F can be any fixed s x m matrix with orthogonal rows, and D is diagonal with D(i,i) = +-1, i.i.d. with probability .5. A simple property of F guarantees that S*A will have low coherence w.h.p.: Thm 4 (AMT) Let F' be any m x m orthogonal matrix, and D = diag(+-1) where each entry is i.i.d. with probability .5. Let f = max_{ij} |F'(i,j)|^2. Then with probability at least .95, mu(F'*D*A) <= C*n*f*log(m) (here C depends on the choice .95). Since 1/m <= f <= 1, mu(F'*D*A) can potentially be as small as O(n/m * log(m)), near the lower bound n/m, which means, by Thm 3, that we only need to sample O(n*polylog(m)) rows of F*D*A to get S*D*A = A' such that inv(R') (where A'=Q'*R') is a good preconditioner with high probability. There are a number of good choices of F' including FFT, DCT, Walsh-Hadamard, etc. that minimize f and can be applied to a column of A in m*log(m) flops. Thus the overall cost of the algorithm for a dense matrix is O(m*n*log(m) + n^3) In the paper, AMT show speedups of up to 4.5x over the dense QR algorithm. (7) (CW / R1) Finally we discuss the Clarkson-Woodruff algorithm. Let eps < 1 be an error tolerance and s = O((n/eps)^2*log^6(n/eps)). Let S = F*D where F is s x m is a uniformly randomly selected column of eye(s), and D is diagonal with D(i,i) = +-1 i.i.d. with probability .5. Thus forming A'=S*A costs just O(nnz(A)). To repeat the analysis stated above: Thm 1: With probability at least 2/3, norm(A*x-b,2) is at most (1+eps) times larger than the true minimum, with cost O(nnz(A) + "lower order terms") where the "lower order terms" grow like O(n^3*log^6(n/eps)/eps^2).