Logic I, Fall 2009, Session 14 Handout

18.409 The Behavior of Algorithms in Practice 28 Feb. 2002 Lecture 6 1 Outline • When Gaussian elimination with partial pivoting fails. • Analysis of complete pivoting. • Speeding up the solution of linear systems. 2 2.1 When partial pivoting fails An example In lecture 3, we saw an example of a linear system that Matlab fails to solve. Here is a similar system: ⎛ 1 ⎜ ⎜−1 ⎜ ⎜ ⎜−1 ⎜ ⎜−1 ⎝ −1 ⎞ ⎛ ⎞ 1 ⎟ ⎜ ⎟ ⎟ ⎜ ⎟ 1 0 0 1⎟ ⎜1⎟ ⎟ ⎜ ⎟ −1 1 0 1⎟ x = ⎜1⎟ . ⎟ ⎜ ⎟ ⎜1⎟ −1 −1 1 1⎟ ⎠ ⎝ ⎠ −1 −1 −1 1 1 0 0 0 1 We will solve this system using Gaussian elimination with partial pivoting. First, we permute the rows so that the (1, 1)-entry has the largest magnitude in the first column. Actually, we do not need to permute; our pivot is already in the right place. Eliminating the nonzero entries in the first column, ⎛ 1 0 0 0 ⎜ ⎜0 1 0 0 ⎜ ⎜ ⎜0 −1 1 0 ⎜ ⎜0 −1 −1 1 ⎝ 0 −1 −1 −1 1 we get ⎞ 1 ⎟ 2⎟ ⎟ ⎟ 2⎟ . ⎟ 2⎟ ⎠ 2 The rest of the reduction proceeds similarly, obtaining ⎛ ⎞ 1 0 0 0 1 ⎜ ⎟ ⎜0 1 0 0 2 ⎟ ⎜ ⎟ ⎜ ⎟ ⎜0 0 1 0 4 ⎟ . ⎜ ⎟ ⎜0 0 0 1 8 ⎟ ⎝ ⎠ 0 0 0 0 16 The exponentially large magnitude of the lower-right entry is an artifact of the reduction. Because of rounding, the lower-right entry of the original matrix may not be recoverable from the reduced matrix. (Technically, Matlab solves this system, because every entry in the last column is exactly a power of two. See the file kahan2.m for a very similar system that Matlab cannot solve.) 2.2 Smoothed analysis What does smoothed analysis tell us about the example? On one hand, we have already seen that if an entire matrix has noise in it, then Gaussian elimination requires relatively few bits of precision, even without pivoting. But on the other hand, if we perturb only the nonzero entries in the above system, then the error still creeps in. Hence, our current application of smoothed analysis to Gaussian elimination fails for this example. Applying smoothed analysis to Gaussian elimination with complete pivoting may resolve this issue. 3 Complete pivoting Gaussian elimination with complete pivoting permutes rows and columns to ensure that the pivot is the largest magnitude entry in the entire submatrix that remains to be row reduced. It is easy to check that complete pivoting guarantees �L�∞ < 1. The algorithm also gives Theorem 1 (Wilkinson). 1 �U �∞ ≤ n 4 ((log n)+1) . �A�∞ Therefore, complete pivoting requires no more than log2 n bits of precision for best possible accuracy. 2 However, people do not use complete pivoting in practice. Partial pivoting rarely fails, and complete pivoting requires twice as many floating point operations to find the largest magnitude entry (the pivot). We should also note that it may be possible to improve upon Wilkinson’s theorem. Nearly always, the growth factor is actually less than n. Perhaps by assuming that there is noise in the input, we can find a better bound. Or perhaps we can calculate a better bound directly. Proof of Wilkinson. Let A(r) denote the lower-right r-by-r submatrix that is obtained after n − r − 1 eliminations. Assume that complete pivoting has already been performed, so that the largest magnitude entry in A(r) is the upper-left entry. Let pr denote this entry, i.e. the pivot. Then �A�∞ = pn and �A(r) �∞ = pr , so we need to bound |pr |/|pn |. (Note that maxr �A(r) � ∞ �A�∞ �U �∞ �A�∞ ≤ .) Wilkinson’s theorem will be a consequence of the next theorem. 2 Theorem 2. 1 1/2 1/3 1/(n−1) + 1 log n |p1 | ) 2 . ≤ 2 2 log2 (2·3 ·4 ·····n |pn | Proof. The result follows from two facts: (r) det(A )= r � pi , (1) i=1 √ det(A(r) ) ≤ ( r |pr |)r . (2) (Inequality (2) is Hadamard’s inequality.) Immediately, we have r � √ pi ≤ ( r |pr |)r i=1 If qr = log |pr |, then �r i=1 qi ≤ r 2 log r + rqr , so r−1 � i=1 qi ≤ r log r + (r − 1)qr , 2 (3) which implies r−1 � i=1 � 1 � q log r qr 1 qi r + = log r r−1 + . ≤ r(r − 1) 2(r − 1) r 2 r 3 (4) Now, rewriting (1), we obtain �n i=1 qi n−1 = log(det(A)) . n−1 (5) Summing inequality (4) over r = 2, . . . , n − 1 with (5), we get n � � n−1 � � qr qn 1 log(det(A)) qr−1 1/2 1/3 1/(n−2) + ≤ log 2 · 3 · 4 · · · · · (n − 1) + + , 2 r−1 n−1 r n−1 r=2 r=2 so � � log(det(A)) qn 1 1/2 1/3 1/(n−2) q1 + ≤ log 2 · 3 · 4 · · · · · (n − 1) + . 2 n−1 n−1 Using (1) and (3), we can bound the last term in the RHS: n log(det(A)) 1 � n = qi ≤ log n + qn . n−1 n−1 2(n − 1) i=1 Thus, q1 + � � qn 1 n ≤ log 2 · 31/2 · 41/3 · · · · · (n − 1)1/(n−2) + log n + qn , 2 n−1 2(n − 1) so qn (n − 2) n−1 � � 1 1 n−1 ≤ log 2 · 31/2 · 41/3 · · · · · (n − 1)1/(n−2) + log n1/(n−1) + log n, 2 2 2(n − 1) q1 − qn ≤ q1 − and finally, q1 − qn ≤ 4 � � 1 1 log 2 · 31/2 · 41/3 · · · · · (n − 1)1/(n−2) · n1/(n−1) + log n. 2 2 2 Speeding up the solution of linear systems If a matrix A = (aij ) is banded within b diagonals of the main diagonal, i.e. aij �= 0 implies |i − j | ≤ b − 1, then Gaussian elimination with partial pivoting can reduce the matrix in time O(nb2 ) and space O(nb). (You can check that the upper-triangular matrix U in the factorization is banded by 2b (approximately).) Given an arbitrary matrix, we would like to permute the rows and columns so that the matrix is banded. If the matrix is symmetric, then this problem becomes a graph problem. Define the graph of a symmetric n-by-n A = (aij ) to be an undirected graph on n vertices in which there is an edge (i, j) if and only if aij �= 0. (The diagonal of A is irrelevant.) 4 4.1 4.1.1 Breadth-first search Cuthill-McKee ’69 One of the first algorithms to transform matrices into banded matrices was developed by Cuthill and McKee in 1969. Given an arbitrary vertex v, the algorithm is as folows: 1. Perform a breadth-first search on the graph of A rooted at v. Label each vertex with its distance from v, so that the vertices are partitioned into distinct “levels.” 2. Order the vertices by level, i.e. find a permutation σ such that if σ(i) < σ(j), then the level of i is no higher than the level of j. 3. Permute the rows and columns of A according to the permutation σ. Note that if the BFS partitions the vertices so that each level contains relatively few vertices, then the resulting matrix is banded into a relatively small region. To be precise, if Li is the set of vertices in level i, then the bandwidth of the permuted matrix is ≤ 3 · maxi |Li |. 4.1.2 Turner ’86 In 1986, Turner analyzed Cuthill-McKee, finding an improvement in the process. His im­ provement was 1. Start with a random root v. 2. Run Cuthill-McKee. 3. Let w be the last vertex (corresponding to row and column n). 4. Run Cuthill-McKee from w. Turner’s analysis was probablistic. For a fixed b, he constructed a random (0, 1)-matrix in which each entry within b diagonals from the main diagonal was 1 with probability p and all other entries were 0. He found that for p ≥ log n n , the bandwidth of the matrix was almost certainly b. Cuthill-McKee produced a permuted matrix with bandwith about 3b. In contrast, Turner’s improvement achieved a bandwith of b + O(log n). 5 4.1.3 Feige-Krauthgamer ’98 Feige and Krauthgamer performed a more general analysis. Starting with an arbitrary (0, 1)-matrix of bandwidth b, they flipped entries in the band with probability p. (The case in which the original matrix is the zero matrix is Turner’s analysis.) Feige and Krauthgamer’s analysis deals with a certain class of perturbed matrices, similar to smoothed complexity analysis. 4.2 Graph Separators For a graph G, an α-vertex separator is a subset C ⊆ V such that V \ C is the disjoint union of sets A and B with no edges between A and B, and |A|, |B | ≤ α|V |. If we can find an α-vertex separator, then we can make progress toward rearranging the matrix into banded form, since there are no edges between the vertices in A and the vertices in B. Theorem 3 (Lipton-Tarjan). Every n-node planar graph has a √ satisfying |C | ≤ 8n. 2 3 -vertex separator C Definition 4. If S is a family of graphs closed under subgraph, then S satisfies an f (n)separator theorem if there exists an α < 1 such that all graphs in S have α-separators of size ≤ f (|V |). Theorem 5 (Lipton-Rose-Tarjan). If S has cnγ -separators, then if G is the non-zero structure of a linear system, and G ∈ S, then we can solve G in space and time given by the following table. γ= γ> 4.3 1 2 1 2 space = fill time O(n log n) O(n3/2 ) n2γ n3γ Non-planar graphs What if the graph of a matrix is not planar? There are a number of heuristic contenders. One of the first was developed by Kernighan and Cuthill. A newer, promising approach uses the eigenvector structure of the Laplacian matrix of the graph of A. Each diagonal entry of a Laplacian matrix of a graph is the degree of the corresponding vertex. The off-diagonal entries are negated edge weights. 6 The Laplacian matrix is positive semidefinite, with 0 as an eigenvalue with multiplicity 1 (if the graph is connected). The eigenvector corresponding to the second smallest eigenvalue can be used to order the vertices. The ordering can be used to cut the graph efficiently. Details will be provided in the next lecture. 7