Chapter 13 Gaussian Elimination III BunchParlett diagonal pivoting

Chapter 13 Gaussian Elimination (III) Bunch-Parlett diagonal pivoting Speaker: Lung-Sheng Chien Reference: James R. Bunch and Linda Kaufman, Some Stable Methods for Calculating Inertia and Solving Symmetric Linear Systems, Mathematics of Computation, volume 31, number 137, January 1977, page 163 -179

Out. Line • Preliminary - 1 x 1, 2 x 2 pivoting in Gaussian Elimination - real symmetric indefinite matrix • • Symmetric permutation LDL’ decomposition (diagonal pivoting) Example of complete diagonal pivoting Algorithm of complete diagonal pivoting

Pivoting in Traditional Gaussian Elimination consider then one step LU-factorization leads to and where It is easy to show We use principal submatrix Question: can we use from Gaussian Elimination as pivoting, called or higher pivoting, i. e. pivoting since

2 x 2 pivoting in block Gaussian Elimination If principal sub-matrix then one-step GE leads to Example 1: where pivoting and is non-singular,

3 x 3 pivoting in block Gaussian Elimination Example 2: pivoting and where Question: what is criterion to choose 1 x 1, 2 x 2 or 3 x 3 pivoting? stability performance Question: what is the cost to estimate the criterion ?

Real symmetry indefinite matrix A is symmetric implies and From linear algebra, if A is real symmetric, then A has real spectrum decomposition (譜分解) and or write in matrix form and without loss of generality, we can rearrange eigen-values to be increasing Definition: suppose matrix A is real symmetric, then we say “A is indefinite” if there exists an eigen-value of A less than zero. If A is nonsingular, then Thereafter we focus on LU-decomposition of symmetry indefinite matrix A

LU-factorization for real symmetry indefinite matrix A factorization where and Question: why not LU-decomposition? A is real symmetric, we only store lower triangle part of A, say

Pointer array to store lower triangular part of A row-major based col-major based Question: what is the cost to fetch one element of A ? That is, operation count for Question: can you find another representation for lower part of matrix A ? Question: if one uses row-major to store lower part of matrix A, then how to fetch a column of A efficiently?

Out. Line • Preliminary • Symmetric permutation - change diagonal element to (1, 1) position - change off-diagonal element to (2, 1) position - implementation of symmetric permutation • LDL’ decomposition (diagonal pivoting) • Example of complete diagonal pivoting • Algorithm of complete diagonal pivoting

How to do permutation on real symmetry indefinite matrix A Let Define permutation matrix interchange row 2, 3 interchange column 2, 3 Symmetric permutation: 1 2 3 2 4 5 3 5 6 1 2 3 interchange 3 5 6 row 2, 3 2 4 5 1 3 2 interchange 3 6 5 column 2, 3 2 5 4

Change diagonal element to (1, 1) position suppose we want to change [1] to position (1, 1) , then consider permutation and do symmetric permutation on A, say 1 2 3 4 2 5 6 7 3 6 8 9 4 7 9 10 3 2 1 4 interchange 6 5 2 7 col 1, 3 8 6 3 9 9 7 4 10 8 6 3 9 interchange 6 5 2 7 row 1, 3 3 2 1 4 9 7 4 10 Question: we only store lower triangle part of matrix A, above permutation does not work, how to modify it? Observation 1: 1 2 3 4 8 2 3 4 2 5 6 7 3 6 8 9 3 6 1 9 4 7 9 10

Change diagonal element to (1, 1) position [2] Observation 2: we only need to update lower triangle part of A (diagonal term is excluded ). 8 2 3 4 2 5 6 7 3 6 1 9 4 7 9 10 2 6 4 7 does not changed diagonal term we have changed does not changed 9

Change diagonal element to (1, 1) position 2 2 interchange 6 7 6 2 6 col 1, 3 6 4 [3] 9 interchange 2 6 9 7 6 2 2 row 1, 3 7 9 4 4 Keep lower triangle part 8 8 6 5 3 2 1 9 7 4 Fill-in A(3, 1) 10 6 9 Fill-in 5 2 1 7 4 6 diagonal term 10 Question: any compact form to represent above interchange? 2 9 7 4

Change diagonal element to (1, 1) position 2 6 6 4 7 case 1: case 2: 2 9 4 4 7 2 6 9 7 9 9 7 4 [4]

Change diagonal element to (1, 1) position [5] For general real symmetric matrix A with dimension n, suppose we want to change interchanges rows or columns Step 1: interchange position (1, 1) and (k, k), i. e. keep position (k, 1) and (1, k), i. e. does not changed

Change diagonal element to (1, 1) position Step 2: interchange column 1 and k below row k+1, i. e Therefore [6]

Change diagonal element to (1, 1) position [7] Step 3: interchange column 1 (from row-2 to row-(k-1)) and row k (from col-2 to col-(k-1)) Therefore

Implementation of symmetric permutation: swapping overhead [1] Suppose we use col-major based pointer array to store real symmetric matrix A Step 2: interchanging column 1 and k is O. K since memory block is contiguous. Contiguous memory block

Implementation of symmetric permutation: swapping overhead Step 3: interchanging column 1 and row k, i. e. is NOT efficient since is NOT contiguous, swapping is very slow. Observation: the situation is the same even we choose row-major based pointer array as storage strategy. Solution: in order to avoid high swapping overhead, we should avoid permutation. [2]

Change off-diagonal element to (2, 1) position suppose we want to change [1] to position (2, 1), then consider first permutation which change first index from 4 to 2 1 2 3 4 2 5 6 7 3 6 8 9 4 7 9 10 4 4 2 7 6 5 interchange 4 col 2, 4 3 9 8 6 row 2, 4 3 9 8 6 10 9 7 2 7 6 5 9 3 6 10 4 7 4 which change second index from 3 to 1 4 7 interchange 9 col 1, 3 8 9 6 7 6 2 7 6 5 7 interchange 10 9 8 10 9 2 2 9 2 3 3 4 and second permutation 1 4 1 1 3 1 2 10 4 7 interchange 9 3 6 row 1, 3 3 4 1 2 2 5 6 7 2 5 8

Change off-diagonal element to (2, 1) position 1 2 3 4 2 5 6 7 interchange 4 3 6 8 9 row and col 2, 4 3 9 4 7 9 10 2 7 1 4 [2] 9 3 2 10 9 7 interchange 9 8 6 row and col 1, 3 3 4 1 2 6 5 6 7 2 5 8 3 6 10 4 7 Observation: two permutation matrices can be computed easily since 2 4 and 1 3 are independent.

Change off-diagonal element to (2, 1) position HOW to transform [3] to position (2, 1), under only lower triangle part of A is stored? Observation 1: and does not changed 1 2 3 4 2 5 6 7 interchange 4 3 6 8 9 row/col 2, 4 3 9 8 6 4 7 9 10 2 7 6 5 2 3 4 10 6 7 1 4 3 2 10 9 7 Step 1: interchange position (4, 4) and (2, 2), 1 2 3 4 1 2 5 6 7 2 3 6 8 9 4 7 9 10 4 7 9 5

Change off-diagonal element to (2, 1) position [4] Observation 2: we only need to update lower triangle part of A (diagonal term is excluded ). 2 1 2 3 4 10 6 7 2 3 6 8 9 3 4 7 9 5 4 6 9 Step 2: interchange row 2 and 4 below col 2, i. e 2 3 4 6 9 interchange 4 row 2, 4 3 2 9 6 6 interchange 4 col 2, 4 3 9 9 2 NOTE: Interchanging column 2, 4 does not change position (2, 1) and (4, 1), it suffices to interchange position (2, 1) and (4, 1) first. 6 6

Change off-diagonal element to (2, 1) position 2 3 [5] 4 6 3 9 2 9 Step 3: interchange column 2 (from row-3 to row-3) and row 4 (from col-3 to col-3) 4 3 2 4 6 3 9 2 9 6

Change off-diagonal element to (2, 1) position [6] For general real symmetric matrix A with dimension n, suppose we want to change interchanges rows or columns Step 1: interchange position (2, 2) and (k, k), i. e. keep position (k, 2) and (2, k), i. e. does not changed

Change off-diagonal element to (2, 1) position Step 2: interchange column 2 and k below row k+1, i. e Therefore [7]

Change off-diagonal element to (2, 1) position Step 3: interchange row 2 (column 1) and row k (column 1), i. e Therefore [8]

Change off-diagonal element to (2, 1) position [9] Step 4: interchange column 2 (from row-3 to row-(k-1)) and row k (from col-3 to col-(k-1)) Therefore

Out. Line • Preliminary • Symmetric permutation • LDL’ decomposition (diagonal pivoting) - growth rate of 1 x 1, 2 x 2 pivoting - criterion for pivot strategy - comparison to complete pivoting in GE • Example of complete diagonal pivoting • Algorithm of complete diagonal pivoting

LDL’ decomposition: 1 x 1, 2 x 2 pivoting • diagonal pivoting method with complete pivoting: Bunch-Parlett, “Direct methods fro solving symmetric indefinite systems of linear equations, ” SIAM J. Numer. Anal. , v. 8, 1971, pp. 639 -655 • diagonal pivoting method with partial pivoting: Bunch-Kaufman, “Some Stable Methods for Calculating Inertia and Solving Symmetric Linear Systems, ” Mathematics of Computation, volume 31, number 137, January 1977, page 163 -179

Growth rate in LDL’ decomposition [1] is of order maximal element of matrix A Define and constant maximal diagonal element of matrix A Case 1: s = 1 with where with and with

Growth rate in LDL’ decomposition [2] Therefore such that Observation: in PA=LU, we choose implies Therefore Observation: large element growth will not occur for a is large relative to pivot if

Growth rate in LDL’ decomposition [3] Case 2: s = 2 where and with and implies with

Growth rate in LDL’ decomposition for Therefore [4]

Growth rate in LDL’ decomposition Therefore [5]

Criterion for pivot strategy Fix a number [1] ( later on, we will determine optimal value of ) maximal element of matrix A maximal diagonal element of matrix A Case 1: suppose we interchange 1 st and k-th rows and columns by introducing permutation matrix Then do on and do symmetric permutation with pivot, written as

Criterion for pivot strategy Recall for pivot, [2] , we have growth rate and Now for and Case 2: suppose , we interchange q-th and 1 st rows and columns, and then r-th and 2 nd rows and columns by introducing permutation matrix and do symmetric permutation Question: we must change q 1 first, then change r 2, why?

Criterion for pivot strategy [3] transforms 1 2 (boundary case) 3 but

Criterion for pivot strategy [4] pivot Then do Recall for on pivot, with pivot, written as , we have growth rate and

Criterion for pivot strategy [5] Now for To sum up Case 1: and Case 2: and Question: How to choose value of

Criterion for pivot strategy worst case analysis : choose growth rate of pivot + or equivalently growth rate of such that pivot growth rate of pivot Exercise: verify pivot growth rate of Define where [6] satisfies and pivot

Diagonal pivoting versus complete pivoting of Gaussian Elimination [1] 1 We must search for the largest element in each reduced matrix, this is a complete pivoting strategy analogous to Gaussian Elimination with complete pivoting maximal element of matrix A maximal diagonal element of matrix A complete pivoting of GE (Gaussian Elimination) such that (1) find a (2) swap row and row (2) swap column then and column with

Diagonal pivoting versus complete pivoting of Gaussian Elimination [2] 2 growth rate diagonal pivoting complete pivoting in GE Remark: Bunch [1] proves that the element growth in the diagonal pivoting method with complete pivoting is bounded by in comparison with for Gaussian Elimination with complete pivoting, where [1] J. R. Bunch, “Analysis of the diagonal pivoting method”, SIAM J. Numer. Anal. , v. 8, 1971, pp. 656 -680

Out. Line • • • Preliminary Symmetric permutation LDL’ decomposition (diagonal pivoting) Example of complete diagonal pivoting Algorithm of complete diagonal pivoting

Example (complete pivoting) -6 6 12 12 -8 -13 4 3 -13 -7 1 -6 4 3 1 6 6 we choose swap row/column and -8 12 6 3 -13 -7 -13 3 -7 4 -6 1 -6 4 1 [1] 6 6 pivot by permutation matrix -8 -13 12 4 -13 -7 3 1 12 3 6 -6 4 1 -6 6

Example (complete pivoting) -8 -13 12 4 -13 -7 3 1 12 3 6 -6 0. 3982 -1. 1681 4 1 -6 6 0. 1327 -0. 3894 1 1 [2] -8 -13 -7 4. 7257 -6. 4248 5. 8584 Recursively do the same procedure for we choose swap row/column by permutation matrix pivot such that

Example (complete pivoting) Do symmetry permutation for [3] and -8 -13 -13 -7 4. 7257 -6. 4248 5. 8584 -6. 4248 4. 7257 1 1 1 0. 3982 -1. 1681 0. 1327 -0. 3894 we do 1 1 pivot on 1 0. 1327 -0. 3894 0. 3982 -1. 1681 5. 8584 -6. 4248 4. 7257 1 1

Example (complete pivoting) 5. 8584 -6. 4248 1 -6. 4248 4. 7257 -1. 0967 [4] 5. 8584 1 1 -2. 3202 -1. 0967 1 Or write in original matrix form -8 -13 -7 1 1 5. 8584 -6. 4248 1 -6. 4248 4. 7257 -1. 0967 -8 -13 -7 5. 8584 1 -2. 3202

Example (complete pivoting) 1 1 0. 1327 -0. 3984 0. 3982 -1. 1681 -1. 0967 -8 -13 -7 1 5. 8584 1 -2. 3202 In practice, we have two key issues 1 [5] col-major based we only store lower part of matrix 6 12 -8 3 -13 -7 -6 4 1 6 6 -8 -7 12 -13 1 3 4 -6 6

Example (complete pivoting) 2 we store decomposition and [6] into original 1 -8 1 -13 0. 1327 -0. 3984 0. 3982 -1. 1681 -1. 0967 -7 1 5. 8584 1 -2. 3202 -8 -13 -7 0. 1327 -0. 3984 0. 3982 -1. 1681 -1. 0967 -2. 3202 5. 8584 Question: How can you judge correct decomposition from

Example (complete pivoting) [7] Case 1: four 1 x 1 pivot 1 -8 13 1 -7 0. 1327 -0. 3984 0. 3982 -1. 1681 -1. 0967 1 5. 8584 1 -2. 3202 Case 2: two 2 x 2 pivot 1 -8 1 0. 1327 -0. 3984 0. 3982 -1. 1681 -13 1 -7 5. 8584 1 -1. 0967 -2. 3202

Example (complete pivoting) [8] Case 3: 1 x 1 pivot + 2 x 2 pivot + 1 x 1 pivot -8 1 -13 -7 1 0. 1327 0. 3982 -0. 3984 1 -1. 1681 -1. 0967 5. 8584 -2. 3202 1 Case 4: 2 x 2 pivot + 1 x 1 pivot 1 -8 1 -13 0. 1327 -0. 3984 1 0. 3982 -1. 1681 -1. 0967 5. 8584 1 -2. 3202 Solution: we need an array to record pivot sequence 1 x 1 pivot 2 0 -7 1 1 2 x 2 pivot

Out. Line • • • Preliminary Symmetric permutation LDL’ decomposition (diagonal pivoting) Example of complete diagonal pivoting Algorithm of complete diagonal pivoting

Algorithm ( PAP’ = LDL’ ) Given symmetric indefinite matrix use permutation vector let , [1] , construct initial lower triangle matrix to record permutation matrix and while we have compute update original matrix , where combines all lower triangle matrix and store in

Algorithm ( PAP’ = LDL’ ) compute 1 [2] and Case 1: define permutation to do symmetric permutation 2 To compute , we only update lower triangle of 1 3 2 1 3 3 then 2

Algorithm ( PAP’ = LDL’ ) To compute 4 We only update lower triangle matrix then 5 do 1 x 1 pivot : then 6 and [3]

Algorithm ( PAP’ = LDL’ ) [4] Case 2: define permutation to interchange q-th and k-th rows and columns, and then r-th and (k+1)-th rows and columns and then 7 To compute , we only update lower triangle of (1) do interchange row/col 1 8 2 1 3 3 then 2

Algorithm ( PAP’ = LDL’ ) 9 [5] (2) do interchange row/col 1 2 3 4 then where 2 1

Algorithm ( PAP’ = LDL’ ) To compute 10 (1) do interchange row 11 (2) do interchange row then [6]

Algorithm ( PAP’ = LDL’ ) 12 [7] do 2 x 2 pivot : then 13 and means endwhile is 2 x 2 pivot

Question: recursion implementation normal abnormal Extension of algorithm • Initialization - check algorithm holds for k=1 • Recursion formula - check algorithm holds for k+1, if k-1 is true • Termination condition - check algorithm holds for k=n-1 • Breakdown of algorithm - check which condition PAP’=LDL’ fails • No extension: algorithm works only for square, symmetric indefinite matrix.

MATLAB implementation [1] Given a (full) symmetric indefinite matrix , compute factorization return four quantities Remark: try to neglect upper triangle part of in MATLAB implementation Example : 12 12 -8 -13 4 3 -13 -7 1 -6 4 3 -6 6 1 2 3 4 1 2 0 1 1 return : 6 1 -8 1 -13 0. 1327 -0. 3984 1 0. 3982 -1. 1681 -1. 0967 -7 5. 8584 1 -2. 3202

MATLAB implementation [2] When factorization is complete, 1 define 2 define , you need to provide linear solver by forward substitution , then by diagonal block inversion , then Example: 3 define , then by backward substitution However you cannot transpose 4 , you must scan explicitly in MATLAB once to determine