Notes for Ma221 Lecture 6, Feb 15, 2022 Goals for today: Gaussian elimination for special structures, that let you save flops (time) and space if you recognize them. Symmetric positive definite matrices (Cholesky) save half flops, space vs GEPP Symmetric matrices save half flops, space vs GEPP Band matrices cost goes from O(n^3) to O(bw^2*n), space from n^2 to O(bw*n) Sparse matrices cost, space can drop a lot, depends a lot on sparsity pattern complicated algorithms, many software libraries available "Structured matrices" dense, but depend on O(n) parameters, e.g. Vandermonde V(i,j) = x(i)^(j-1), arising in polynomial interpolation Toeplitz T(i,j) = t(i-j), so constant along diagonals, arising in time series there are O(n^2) or faster algorithms for many of these Can imagine as many special structures as individual problems, above are most common SYMMETRIC (HERMITIAN) POSITIVE DEFINITE Def: A real and spd iff A=A^T and x^T*A*x > 0 for all x neq 0; A complex and Hpd iff A = A^H and x^H*A*x > 0 for all x neq 0; Lemma (just for the real case) 1: X nonsing => A spd iff X^T*A*X spd Pf: A spd and x neq 0 => X*x neq 0 => (X*x)^T*A*X*x = x^T*X^T*A*X*x neq 0 other direction analogous. 2: A spd and H=A(j:k,j:k) => H spd (H called a "principal submatrix") Pf: A spd and y neq 0 => x = [0;y;0] neq 0 where y in rows j:k => x^T*A*x = y^T*H*y neq 0 3: A spd iff A=A^T and its eigenvalues lambda_i > 0 Pf: A spd => A*Q = Q*Lambda, Q orthogonal eigenvector matrix and Lambda diagonal matrix of eigenvalues => Q^T*A*Q = Lambda so by 1 above, A is spd iff Lambda is spd, i.e. x^T*Lambda*x = sum_i x(i)^2*lambda(i) > 0 for all nonzero x, which is clearly true iff all lambda(i) > 0. 4: A spd => A(i,i)>0 and max_ij |A(i,j)| = max_i |A(i,i)| (largest entry on diagonal) Pf: e(i)^T*A*e(i) = A(i,i) > 0. Suppose A(i,j) with i neq j were the largest entry in absolute value. Let x be a vector of zeros except for x(i)=1 and x(j)=-sign(A(i,j)). Then x^T*A*x = A(i,i) + A(j,j) - 2*abs(A(i,j)) <= 0, a contradiction. 5: A spd iff A=L*L^T where L lower triangular with positive diagonal entries Pf: We use induction on n, showing Schur complement spd: Let A = [ A11 A12 ] = [ sqrt(A11) 0 ] * [ sqrt(A11) A12/sqrt(A11) ] [ A12' A22 ] [ A12'/sqrt(A11) I ] [ 0 S ] where A11 is 1 x 1, A12' is A12^T (for short), and A22 = S + A12'*A12/A11 so S = A11 - A12'*A12/A11 = S^T, and A = [ sqrt(A11) 0 ] * [ 1 0 ] * [ sqrt(A11) A12/sqrt(A11) ] [ A12'/sqrt(A11) I ] [ 0 S ] [ 0 I ] = X^T * [ 1 0 ] * X [ 0 S ] So [ 1 0 ] is spd, and so S is spd. [ 0 S ] By induction S = Ls * Ls^T where Ls lower triangular with positive diagonal, and A = [ sqrt(A11) 0 ] * [ 1 0 ] * [ sqrt(A11) A12/sqrt(A11) ] [ A12'/sqrt(A11) I ] [ 0 Ls*Ls^T ] [ 0 I ] = [ sqrt(A11) 0 ] * [ sqrt(A11) A12/sqrt(A11) ] [ A12/sqrt(A11) Ls ] [ 0 Ls^T ] = L*L^T Def: Suppose A spd, then A = L*L^T with L lower triangular with positive diagonal is called the Cholesky factorization Relationship of Cholesky to LU decomposition without pivoting: Lemma: If A spd and A=L*U with L unit lower triangular, let D be diagonal with D(i,i) = sqrt(U(i,i)) > 0. Then A = (L*D)*(inv(D)*U) = Ls*Ls^T is the Cholesky factorization. Proof: Homework! This lemma suggests that any algorithm for LU can be modified to do Cholesky, and in fact use half the work and space since L and U are simply related. Here is the simplest version: Cholesky algorithm to compute A = L*L^T just follows induction proof: for j = 1:n L(j,j) = sqrt( A(j,j) - sum_{i=1}^{j-1} L(j,i)^2 ) L(j+1:n,j) = (A(j+1:n,j) - L(j+1:n,1:j-1)*L(j,1:j-1)^T)/L(j,j) All other ideas for LU (blocking, recursion, communication lower bounds etc) apply. Error analysis: Same approach as to LU yields analogous bound: (A+E)(x+dx) = b where |E| <= 3*n*macheps*|L|*|L^T| But since we do not pivot, we need to use a different approach to get a bound: (|L|*|L|^T)_{ij} = sum_k |L_{ik}|*|L_{jk}| <= norm(L(i,:))*norm(L(j,:)) ... Cauchy-Schwartz inequality = sqrt(A_{ii})*sqrt(A_{jj}) <= max entry in A SYMMETRIC INDEFINITE It is also possible to save half the space and flops by just preserving symmetry, without assuming positive definiteness. The traditional approach is called Bunch-Kaufman factorization: A = P*L*D*L^T*P^T where P is a permutation, L is unit lower triangular, and D has 1-by-1 and 2-by-2 diagonal blocks. More complicated pivoting (with 2x2 blocks) is needed to preserve symmetry: consider A = [[0 1];[1 0]]. Instead of just searching a column for a pivot, once needs to search along a row too. It is more complicated to apply ideas from LU and Cholesky to this factorization (blocking, recursion, communication lower bounds): see LAPACK routine ssytrf.f. A more numerically stable version is called rook-pivoting, which may search more than one row and column of A, but has better numerical stability guarantees, see LAPACK routine ssytrf_rk.f. And there is another approach (called Aasen factorization): A = P*L*T*L^T*P^T where now T is a symmetric tridiagonal matrix, and so costs just O(n) to solve T*x=b. See "Implementing a Blocked Aasen's Algorithm with A Dynamic Scheduler on Multicore Architectures" at bebop.cs.berkeley.edu. Aasen is more amenable to optimizations to minimize communication. (We pause the recorded lecture here.) BAND MATRICES These are the simplest sparse matrices, and have solvers in LAPACK. Def lbw(A) = lower bandwidth of A, ubw(A) = upper bandwidth of A Case without pivoting (eg Cholesky): ubw(U)=ubw(A), lbw(L)=lbw(A) (picture) Cost = 2*n*ubw*lbw + n*lbw = O(n) for narrow bandwidths Case with pivoting: ubw(U)=ubw(A)+lbw(A) "lbw(L)" = lbw(A) (doesn't look banded, but uses same amount of storage) picture, Matlab spyplots: n=20,lbw=4,ubw=3; A = tril(triu(randn(n,n),-lbw),ubw); Ad = A + 100*eye(n); % no pivoting required [Ld,Ud,Pd]=lu(Ad); figure(1), subplot(221),spy(Ad),title('Ad'),subplot(222),spy(Ld),title('Ld'), subplot(223),spy(Ud),title('Ud'),subplot(224),spy(Pd),title('Pd') figure(2), [L,U,P]=lu(A); subplot(221),spy(A),title('A'),subplot(222),spy(L),title('L'), subplot(223),spy(U),title('U'),subplot(224),spy(P),title('P') Band matrices often arise from discretizing differential equations, or other physical simulations, where each unknown only depends on its neighbors, so A(i,j) is nonzero only when i and j are close, i.e. A(i,j) near diagonal. Ex: Sturm-Liouville problem for -y''(x) + q(x)*y(x) = r(x) on an interval [0,1], y(0)=alpha, y(1)=beta, q(x) >= qbar > 0, discretize at x(i) = ih, h=1/(N+1), unknowns y(1) = y(x(1)),...,y(N) = y(x(N)). Approximate y''(i) = (y(i+1) - 2*y(i) + y(i-1))/h^2 so Letting q(i) = q(x(i)) and r(i) = r(x(i)) we get linear system -(y(i+1) - 2*y(i) + y(i-1))/h^2 + q(i)*y(i) = r(i) or A*y = v where v = r + [alpha/h^2 0 ... 0 beta/h^2]' A = diag(2/h^2 + q(i)) ... diagonal + diag(-1/h^2 , 1 ) + diag(-1/h^2 , -1 ) ... so A symmetric We can prove A is spd, so we can use Cholesky, by using Gershgorin's Theorem: All eigenvalues lambda of A lie in n circles in the complex plane, where circle i has center A(i,i) and radius sum_{j=1 to n except i} |A(j,i)| proof: if A*x = lambda*x , suppose |x(i)| = largest entry in absolute value of x. Thus (A(i,i)-lambda) = sum_{j=1 to n except i} A(i,j)*(x(j)/x(i)) and so |A(i,i)-lambda| <= sum_{j=1 to n except i} |A(i,j)|*|x(j)/x(i)| <= sum_{j=1 to n except i} |A(i,j)|*1 as desired Apply Gershgorin to A=A^T: all eigenvalues in circles centered at 2/h^2 + q(i) >= 2/h^2 + qbar, with radius 2/h^2, so must be real and positive. Ex: Now consider Poisson's equation in 2 dimensions: d^2 u(x,y)/dx^2 + d^2 u(x,y)/dy^2 + q(x,y)*u(x,y) = r(x,y) in a square: 0 <= x <= 1 and 0 <= y <= 1 with boundary conditions: u(x,y) given on the boundary of square. Typical example: Given temperature on boundary of a square of uniform material, what is the temperature everywhere in the middle? Here q=r=0. We discretize the square with a nxn mesh: h=1/(n+1) x(i) = i*h for i=1:n, so x(0)=0, x(n+1) = 1 y(j) = j*h for j=1:n, so y(0) = 0, y(n+1) = 1 u(i,j) = u(x(i),y(j)), same for q(i,j) and r(i,j) Again approximate d^2 u(x,y)/dx^2 ~ [ u(x-h,y) - 2*u(x,y) + u(x+h,y) ] / h^2 d^2 u(x,y)/dy^2 ~ [ u(x,y-h) - 2*u(x,y) + u(x,y+h) ] / h^2 Add these, evaluate at (x(i),y(j)), so d^2 u(i,j)/dx^2 + d^2 u(i,j)/dy^2 ~ [ u(i-1,j) + u(i+1,j) + u(i,j-1) + u(i,j+1) - 4*u(i,j)]/h^2 (Illustrate for n=4, 6 using A=TwoDL(n)) Using "natural order" (rowwise or columnwise ordering of mesh points (i,j)), we get an spd band matrix of dimension n^2 and bandwidth n. (Later we will have explicit expressions for all eigenvalues and eigenvectors of this matrix, using the FFT). A is banded but at with most 5 nonzeros per row, so should really think of it as sparse... (We pause the recorded lecture here.) GENERAL SPARSE MATRICES Here is a small matlab example to show important of ordering on "arrowhead matrix" A = eye(8);A(1,2:8)=.1;A(2:8,1)=.1; [L,U]=lu(A);A,L,U, figure(1), clf, spy(A,'k'), pause spy(L,'b'), hold on, spy(U,'r'), spy(A,'k') Result is dense! So performing LU and storing L and U would cost as many flops and words as a dense matrix, O(n^3) and O(n^2). Now let's try LU on P*A*P^T, where P is a permutation that reverses the order of rows and columns Arev=A(8:-1:1,8:-1:1); [Lrev,Urev]=lu(Arev); Arev, Lrev, Urev figure(2), clf, spy(Arev,'k'), pause spy(Lrev,'b'), hold on, spy(Urev,'r'), spy(Arev,'k') Result is very sparse! And if you rewrite LU to take advantage of this, to neither store nor compute entries that are certain to be zero, the cost in flops and words would drop to O(n), much smaller than before. Matlab has many built-in sparse matrix operations (LU, Cholesky, etc), that take advantage of sparsity like this; do "help sparse" to learn more. Let's see what happens to the 2D Poisson equation: A = TwoDL(6); R = chol(-A); L = R'; ... R = upper triangular Cholesky factor figure(1), clf, spy(L,'b'),hold on, spy(A,'k') So the band fills in entirely, costing O(N*bw^2) = N*sqrt(N)^2 = N^2 = n^4 for an nxn mesh with N = n^2 unknowns. We'll see that we can do much better, O(n^3) instead, for a clever choice of permutation P in P*A*P^T. EX: sparse Cholesky on finite element structure (Figures 2.9-2.11 in text) MechStructureMesh1.eps - basic mesh, 483 nodes, one unknown for each node, a(i,j) nonzero if edge connects nodes i and j Def: A weighted, undirected graph G is a collection V,E,W of sets V = vertices (aka nodes), E = edges connected pairs of vertices (u,v) (note: (u,u) is allowed; undirected means (u,v) is the same as (v,u) W = weight (number) associated with each edge There is an obvious one-to-one relationship between a weighted undirected graph G and a symmetric matrix, and we will use the notations interchangeably: V = one row and column per vertex, E = locations of nonzeros (note: (u,u) on diagonal) W = values of nonzeros Examples: MechStructureMesh2.eps - to make a matrix, need to number nodes in some way MechStructureANatural.eps - color coded matrix to match structure => 3971 nonzeros => 3971/483^2 = 1.7% filled MechStructureNatural.eps - original matrix and Cholesky factor nnz = 11533, up from 3971/2, 3e5 flops (vs (1/3)n^3 = 3.7e7) Note the "skyline" structure of matrix, which some sparse Cholesky routines exploit. Red entries are new nonzeros created by factorization, called "Fill-in". MechStructureRCM.eps - original matrix with reverse Cuthill-McKee ordering tends to make bandwidth narrow, nnz = 9073 (better), 1.8e5 flops (better) MechStructureMinDeg.eps - another popular order, nnz = 8440, 2e5 flops MechStructureRand.eps - a bad ordering: nnz = 37906, 5.6e6 flops Discretization of air around wing tells similar story, called airfoil in Matlab (figure 2.12 in text): load('airfoil'), for k=1:length(i), plot([x(i(k)),x(j(k))],[y(i(k)),y(j(k))],'k'), hold on, end (We pause the recorded lecture here.) Challenges to factoring sparse matrices (LU and Cholesky): Storing and operating on only nonzeros in an efficient way: Simplest imaginable data structure: COO = coordinate format: list of nonzero entries and their locations (a(i,j),i,j) in some order (example) A = [ 2 0 7 0 5 ] , [ 0 1 4 0 3 ] [ 0 0 8 0 0 ] COO = ((2,1,1),(7,1,3),(5,1,5),(1,2,2),(4,2,3),(3,2,5),(8,3,3)) We can do better: CSR = Compressed Sparse Row: val = array of nonzero entries in each row, from row 1 to row n, from left to right col_index = columns index of each nonzero (val(i) is in col(a) of sparse matrix) rowbegin = pointer to start of each row: entries of row i lies in entries rowbegin(i) through rowbegin(i+1)-1 of val() and col() may be 2/3 as much memory as COO (example) For same A as above: val = (2 7 5 1 4 3 8), col_index = (1 3 5 2 3 5 3), rowbegin = (1 4 7 8) many other formats: defer details until discussion of iterative methods Picking the "best order" of rows and columns: We want to minimize time, memory, get the right answer (stability). Order effects fill-in, which effects time and memory (less is better!) so we want to pick order to minimize fill-in. Order effects backward stability for LU (and symmetric indefinite LDL^T), but not Cholesky (recall that the stability proof for Cholesky did not depend on order). Cholesky is easiest case: can pick order just to minimize fill in. What is the best order for Cholesky? Thm: Picking the optimal order (out of all n! possibilities) is NP-hard, i.e. any algorithm will have cost exponential in n, overwhelming everything else. So we need to use heuristics: RCM = Reverse Cuthill-McKee = Breadth-First Search order, reversed Consider matrix as Graph Def: A path in a graph from vertex v1 to vk is a set of edges (v1,v2), (v2,v3), ... , (v_k-1,vk) connected end-to-end Def: The distance between any two vertices is the length of the shortest path (fewest edges) connecting them RCM: (1) pick any vertex, call it the root (2) compute distance of all other vertices from root, and label them by distance, (so distance=0 is the root, distance=1 are the vertices with an edge to the root, etc) Cost = O(#nonzeros in matrix), using Breadth-First-Search, cheap! (3) Order vertices in reverse order by distance (farthest first) Fact: Vertices at distance k can only be directly connected to vertices at distance k-1 or k+1 (otherwise distance would be wrong!) Explains why RCM ordering tends to make matrix entries close to diagonal: matrix becomes block tridiagonal, with vertices at distance k making up diagonal block k. So it is like a band matrix, and limits fill-in to the band. Note: called symrcm in Matlab MD = Minimum Degree Def: the degree of a vertex in the graph is the number of edges for which it is an endpoint Fact: if vertex corresponding to row/column i is used as a pivot, and has degree d, then we need to do d^2 muls and adds, and so can fill in at most d^2 entries (if they were originally zero, and ignoring symmetry) MD: at each step pick pivot to have minimum degree (note: as factorization proceeds, need to update matrix and its graph, so degree can change; use latest degree) Question: How does graph change after one elimination step? Variants called amd, symamd, colamd in matlab Nested Dissection General idea: "Bisect" vertices of graph into 3 sets V = V1 u V2 u Vs (s stands for "separator") such that (1) about the same number of vertices in V1 and V2 (2) Vs much smaller than V1 and V2 (3) there are no edges directly connecting vertices in V1 and V2 Many algorithms for this (see graph partitioning lecture from CS267) Now reorder vertices with V1 first, then V2, then Vs: What does matrix look like? A = [ A1 0 A1s ] so after Cholesky L = [ L1 0 0 ] [ 0 A2 A2s ] [ 0 L2 0 ] [ A1s' A2s' As ] [ L1s L2s Ls ] in particular, A1 = L1*L1^T and A2 = L2*L2^T are independent Cholesky factorizations, and (2,1) of L block stays zero. So how do we do Cholesky on A1 and A2? Use same idea recursively. (see CS267, Spring 14, lecture 14, slide #9) Illustrate on matrix from 2D mesh (draw ordering on mesh) (MatLab: m=3 (n=2^m-1), NestedDissection) Thm (George, Hoffman/Martin/Rose, Gilbert/Tarjan, 1970s & 1980s): Any ordering for Cholesky on a 2D n x n mesh (so the matrix is n^2 x n^2) has to do at least Omega(n^3) flops, and this is attained by nested dissection ordering. (If it were dense, it would cost O(n^6), a lot more). This extends to "planar graphs", i.e. graph you can draw on paper without crossing any edges. Thm (Ballard, D., Holtz Schwarz, 2009) Lower bound on #words moved for sparse LU and Cholesky is Omega(#flops/sqrt(M)), where #flops is number actually performed. So the lower bound on #words_moved for Cholesky on nxn mesh is Omega(n^3/sqrt(M)). Thm (David, D., Grigori, Peyronnet 2010). Possible to achieve lower bound for nxn mesh of Omega(n^3/sqrt(M)) words moved by Cholesky with nested dissection ordering, if properly implemented Contrast: band solver would cost O(dimension*bw^2) =n^(2+2)=n^4 flops (see examples from CS267 Lecture 14, www.cs.berkeley.edu/~demmel/cs267_Spr14/lecture14_partition_jwd14.ppt) It's a good idea for 3D meshes too: n x n x n mesh (so matrix is n^3 by n^3) costs n^6 to factor using Cholesky with nested dissection ordering (versus n^9 if it were dense, or n^(3+2+2)=n^7 if banded) Lots of clever algorithms and software available for partitioning graphs: See Lecture 14 in CS267 from Spring 14 Packages include METIS, ParMETIS, Scotch and Zoltan (We pause the recorded lecture here.) To summarize, the overall sparse Cholesky algorithm is as follows: Choose ordering (RCM, MD, ND) Build data structure to hold A and (to be computed) entries of L There is some optional material at the end of the typed notes that goes into more detail about how to do this efficiently. Perform factorization What about nonsymmetric matrices, how do we pivot for sparsity and stability? Threshold pivoting: among pivot choices in column within a factor of 2 or 3 of the largest, pick the one with the least fill-in (so may need to change data structure on the fly) Static Pivoting (used in SuperLU_DIST) (1) Permute and scale matrix to try to make diagonal as large as possible: Thm: For any A there is a permutation P and diagonal matrices D1 and D2 such that B = D1*A*P*D2 has B(i,i) = 1 and all other |B(i,j)| <= 1 (2) Reorder using similar ideas as for Cholesky: this leaves all the diagonal entries on the diagonal. Then build the data structures for L and U, before doing LU, rather than figuring them out on the fly, which would be slower. (3) During factorization, if a prospective pivot is too small, just make it big.(This happens rarely in practice.) As a result we get the LU factorization of A + a low rank change (the rank = the number of modified pivots), which we can fix with the Sherman-Morrison-Woodbury formula or an iterative method like GMRES (to be discussed in Chap 6). This brief discussion suggests that there is a large body of software that has been developed for solving sparse Ax=b, and which one is faster may depend on your matrix (and computer). There is a survey of available sparse matrix software (link on class web page); "Updated Survey of sparse direct linear equation solvers" http://crd.lbl.gov/~xiaoye/SuperLU/SparseDirectSurvey.pdf Finally, we give a brief tour of structured matrices, which may be dense but depend only on O(n) parameters. The number of structures that can arise, and be used to get faster algorithms, is large, eg depending on the underlying physics of the problem being modeled. We will do one set of common examples that share a structure: Vandermonde: V(i,j) = x(i)^(j-1) - arises in polynomial interpolation Cauchy: C(i,j) = 1/(x(i)+y(j)) - arises in rational interpolation Toeplitz: T(i,j) = x(i-j), constant along diagonals, arises in convolution Hankel: H(i,j) = x(i+j), constant along "anti-diagonals", similar to Toeplitz Eg: Solving V*z = b means finding z(i) so that sum_{j=1 to n} x(i)^(j-1)*z(j) = b(i) i.e. finding polynomial coefficients z(j) so to match values b(i) at points x(i), i.e. polynomial interpolation; we know how to do this in O(n^2) time using, eg using "Newton interpolation." Another trick works for V^T*z = b . Eg: Multiplying by a Toeplitz matrix is same as convolution, which we know how to do fast by FFT. Eg: Solving with a Cauchy, also interpretable as interpolation Common structure of X: A*X + X*B = low rank for some simple, nonsingular A and B Def: The rank is called the "displacement rank" Ex: Let D = diag(x(i)), then D*V = "V shifted left" V * P = V* [ 0 0 0 ... 0 1 ] = "V shifted left", [ 1 0 0 ... 0 0 ] [ 0 1 0 ... 0 0 ] ... [ 0 0 0 ... 1 0 ] difference D*V - V*P = nonzero in last column only = rank 1 Ex: T Toeplitz implies that P*T - T*P = T shifted down - T shifted left = nonzero only in first row, last col = rank 2 Ex: C Cauchy implies that diag(x_i)*C + C*diag(y_i) = all ones =rank 1 Theorem (Kailath et al): there are O(n^2) solvers if displacement rank low (stability can be a problem). (We pause the recorded lecture here.) The following optional material addresses one more level of detail, how to build the data structure for the Cholesky factor L more cheaply than just performing the factorization to see what L's nonzero entries are. We will use graph theory, and assume that there is no cancellation, i.e. that once an entry of A gets filled in, it never cancels out later (which would be highly unlikely anyway). To motivate the result, let us ask how L(i,j), i>j, could come to be nonzero. Recall that at elimination step k the current entry of the Schur Complement stored at A(i,j) gets updated to be A(i,j) - L(i,k)*(L^T)(k,j) = A(i,j) - L(i,k)*L(j,k) (1) A(i,j) could be initially nonzero. So in the graph G corresponding to A, there is an edge (i,j), so a path of length 1 from i to j. (2) A(i,j) might be zero, but it may get filled in because there is some pivot step k, k < i and k < j, that fills it in, because A(i,k) and so L(i,k) is nonzero, and A(j,k) and so L(j,k) is nonzero. Thus there are edges (i,k)(k,j) in G, a path of length 2 from i to j where k < i and k < j (3) Next suppose A(i,j) and A(i,k) are zero, but L(i,k) is nonzero because it has been previously filled in by some A(i,k) = A(i,k) - L(i,m)*L(k,m) where m < i and m < k where A(i,m) and so L(i,m) is nonzero and A(k,m) and so L(k,m) is nonzero. Thus there are edges (i,m),(m,k),(k,j) in G, a path of length 3 from i to j with m and k < i and < j. More generally, we get Thm (George & Liu, 1981): L(i,j) is nonzero (assuming no exact cancellation) if there is a path from i to j in the graph G of A passing through vertices k satisfying k < i and k < j. Thm (Gilbert & Peierls) Cost of determining structure of L no more than actually doing the factorization.