480 likes | 737 Views
Chapter 14 Gaussian Elimination (IV) Bunch-Kaufman diagonal pivoting (partial pivoting). Speaker: Lung-Sheng Chien.
E N D
Chapter 14 Gaussian Elimination (IV)Bunch-Kaufman diagonal pivoting(partial pivoting) Speaker: Lung-Sheng Chien Reference: [1] 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 [2] Cleve Ashcraft, Roger G. Grimes, and John G. Lewis, Accurate Symmetric Indefinite Linear Equation Solvers, SIAM J. MATRIX ANAL. APPL. Vol. 20, No. 2, 1998, pp. 513-561
OutLine • Review Bunch-Parlett diagonal pivoting- expensive pivot-selection strategy- pivot induces swapping row/column, not efficient • Partial pivoting: basic idea • Implementation of partial pivoting • Example
Recall Criterion for pivot strategy (complete pivoting) [1] expensive Case 1: Define permutation matrix and do symmetric permutation 1 3 2
Recall Criterion for pivot strategy (complete pivoting) [2] Then do on with pivot, where Therefore maximum of off-diagonal elements in column Observation: if we define , then
Recall Criterion for pivot strategy (complete pivoting) [3] Question: To compute is expensive, if we don’t want to compute it, holds ? and Can we have other choice such that controllability of to Observation: It suffices to change pivoting condition , then and The same as that we do by using pivoting criterion (complete pivoting) However computing is cheaper than since but
Recall Criterion for pivot strategy (complete pivoting) [4] Technical problem: computing does not attain optimal. column-major based Contiguous memory block, fast NOT contiguous memory block, slowly. How to improve this ?
Recall Criterion for pivot strategy (complete pivoting) [5] Second, if , we need to interchange row/column relates to read/write on row vector 3 1 which is NOT contiguous, slowly. 3 2 Question: why we persist on computing , why not choose first column? say and maximum of off-diagonal elements in column 1 Again, if , then we also have and
OutLine • Review Bunch-Parlett diagonal pivoting • Partial pivoting: basic idea- select pivot by at most computation of two columns- potential risk of growth rate of lower triangle matrix • Implementation of partial pivoting • Example
Partial pivoting: basic idea [1] and maximum of off-diagonal elements in column 1 , choose is determined later , then If where Pros (優點): compute and , we need to deal with the case Cons (缺點): It is difficult to satisfy holds. and carefully such that controllability of Comparison with : suppose for we CANNOT move to due to symmetric permutation, It is impossible to keep stability by just computing one column as we do in
Partial pivoting: basic idea [2] Question: how to deal with the case for <ans> from procedure of complete pivoting, we need 2x2 pivot. define permutation , change to 1 4 3 2
Partial pivoting: basic idea [3] we may say
Partial pivoting: basic idea [5] Under the case , we want to find a condition for such that is bounded. Lower bound of upper bound of define maximum of off-diagonal elements in column r 1 2 If , then holds. upper bound of Question: how to achieve lower bound of
Partial pivoting: basic idea [6] From assumption , to keep inequality, the natural choice is Since if we add assumption , then and moreover If we add third assumption 3 , then and then holds. Lower bound of under 2x2 pivoting, we need three conditions To sum up, in order to control (1) (2) (3)
Partial pivoting: basic idea [7] Question: does the three condition keep controllability of (2) (3) (1)
Partial pivoting: basic idea [8] So far, we deal with two cases in the following decision-making tree 1x1 pivot 1x1 pivot 1 ? ? 2 ? ? 3 2x2 pivot 2x2 pivot 4 Question: how to deal with the other two cases and what is the order of decision-making, briefly speaking which tree is correct, left or right?
Partial pivoting: basic idea [9] and maximum of off-diagonal elements in column 1 , we do not compute it explicitly. maximum elements of pivot 1
Partial pivoting: basic idea [10] define maximum of off-diagonal elements in column r Contiguous memory block, fast NOT contiguous memory block, slowly. Remark: when case 1 is not satisfied, we must compute for remaining decision-making. Remark: In general, may be more larger than , that is why case 2 holds 2
Partial pivoting: basic idea [11] 2 under pivot Remark: bound on lower triangle matrix is NOT good since it is proportional to however, the bound of reduced matrix is good.
Partial pivoting: basic idea [12] under pivot 3 Clearly, this is the same as case 1 if we interchange row/column define permutation , change to 1 3 2
Partial pivoting: basic idea [13] Then do on with pivot,
Partial pivoting: basic idea [14] 4 under pivot define permutation , change to 1 3 4 2
Partial pivoting: basic idea [15] Then do on with pivot, Remark: bound on lower triangle matrix is NOT good since it is proportional to however, the bound of reduced matrix is good.
Partial pivoting: basic idea [16] To sum up, maximum of off-diagonal elements in column 1 maximum of off-diagonal elements in column r pivot 1 2 pivot pivot 3 4 pivot
Partial pivoting: basic idea [17] worst case analysis : choose such that reduced matrix satisfies growth rate of pivot + pivot growth rate of pivot or equivalently growth rate of pivot growth rate of pivot Define where satisfies Remark: this optimal value is the same as we have in complete pivoting. however, the bound of lower triangle matrix is NOT good.
Complete pivoting versus partial pivoting , growth rate of reduce matrix contributes to final block diagonal matrix partial pivoting complete pivoting Remark: Bunch in [1] only deal with controllability of reduced matrix However Ashcraft in [2] point out the importance of growth rate of Reference: [1] 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 [2] Cleve Ashcraft, Roger G. Grimes, and John G. Lewis, Accurate Symmetric Indefinite Linear Equation Solvers, SIAM J. MATRIX ANAL. APPL. Vol. 20, No. 2, 1998, pp. 513-561
OutLine • Review Bunch-Parlett diagonal pivoting • Partial pivoting: basic idea • Implementation of partial pivoting • Example
Algorithm ( PAP’ = LDL’, partial pivot ) [1] Given symmetric indefinite matrix , construct initial lower triangle matrix use permutation vector to record permutation matrix and let , while we have compute update original matrix , where combines all lower triangle matrix and store in
Algorithm ( PAP’ = LDL’, partial pivot ) [2] 1 compute Case 1: pivot no interchange 2 do 1x1 pivot : then and 3 When case 1 is not satisfied, we compute maximum of off-diagonal elements in column r is used. under lower triangular part of matrix
Algorithm ( PAP’ = LDL’, partial pivot ) [3] 4 no interchange pivot Case 2: 5 do 1x1 pivot : then and 6 The same as code in case 1
Algorithm ( PAP’ = LDL’, partial pivot ) [4] Case 3: pivot do interchange define permutation to do symmetric permutation 7 To compute , we only update lower triangle of 1 1 8 2 3 3 then 2
Algorithm ( PAP’ = LDL’, partial pivot ) [5] To compute 9 We only update lower triangle matrix then 10 do 1x1 pivot : then and 11
Algorithm ( PAP’ = LDL’, partial pivot ) [6] Case 4: pivot do interchange define permutation , change to 12 To compute , we only update lower triangle of 13 do interchange row/col 1 1 2 3 4 3 4 then 2
Algorithm ( PAP’ = LDL’, partial pivot ) [7] 14 do interchange row then 15 do 2x2 pivot : then and 16 endwhile
Question: recursion implementation • Initialization- check algorithm holds for k=1 • Recursion formula- check algorithm holds for k or 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. normal abnormal Extension of algorithm
OutLine • Review Bunch-Parlett diagonal pivoting • Partial pivoting: basic idea • Implementation of partial pivoting • Example
Example (partial pivoting) [1] -6 1 invalid 6 12 3 4 12 -8 -13 invalid 2 1 3 -13 -7 -6 4 1 6 invalid 3 4 pivot 12 ,since k +1 = r , we don’t need permutation , since k +1 = r , we don’t need permutation 13 do interchange row/col 2 1 3 4
Example (partial pivoting) [2] 14 do interchange row , since k +1 = r , we don’t need permutation 15 do 2x2 pivot : -6 6 12 3 1 6 12 4 12 -8 -13 1 12 -8 1 -5.5 3 -13 -7 -0.6875 0.5938 1 2.7813 -6 4 1 6 0 -0.5 1 -5.5 8 and 16 2 0 0 0
Example (partial pivoting) [3] 6 12 1 invalid 12 -8 invalid 2 -5.5 2.7813 -5.5 8 pivot 3 7 Compute 1 6 12 12 -8 8 2 -5.5 8 -5.5 2.7813 3
Example (partial pivoting) [4] Compute 9 1 1 1 1 -0.6875 0.5938 1 0 -0.5 1 0 -0.5 1 -0.6875 0.5938 1 10 do 1x1 pivot : 6 12 1 6 12 12 -8 1 12 -8 8 1 -5.5 8 -1 -0.6875 1 -5.5 2.7813
Example (partial pivoting) [5] and 11 2 0 1 0 issue proper termination condition. To sum up 1 2 4 3 2 0 1 1 1 6 12 1 12 -8 0 -0.5 1 8 -0.6875 0.5938 -0.6875 1 -1
Example: complete pivot versus partial pivot Partial pivot -6 6 12 3 4 12 -8 -13 1 2 4 3 2 0 1 1 1 3 -13 -7 -6 4 1 6 1 6 1 12 -8 0 -0.5 1 8 Complete pivot -0.6875 0.5938 -0.6875 1 -1 2 3 4 1 2 0 1 1 1 -8 1 -13 -7 0.1327 -0.3984 1 5.8584 0.3982 -1.1681 -1.0967 1 -2.3202
Exercise How about flow-chart of left figure? Bunch-Kaufman proposed flow chart 1x1 pivot 1x1 pivot 1 ? ? 2 ? ? 3 2x2 pivot 2x2 pivot 4