Chapter 13 Gaussian Elimination III BunchParlett diagonal pivoting

1
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

2
OutLine

• Preliminary- 1x1, 2x2 pivoting in Gaussian
  Elimination- real symmetric indefinite matrix
• Symmetric permutation
• LDL decomposition (diagonal pivoting)
• Example of complete diagonal pivoting
• Algorithm of complete diagonal pivoting

3
Pivoting in Traditional Gaussian Elimination
consider
then one step LU-factorization leads to
and
where
It is easy to show
from Gaussian Elimination
We use principal submatrix
as pivoting, called
pivoting since
or higher
pivoting, i.e.
Question can we use
4
2x2 pivoting in block Gaussian Elimination
If principal sub-matrix
is non-singular,
then one-step GE leads to
pivoting
Example 1
where
and
5
3x3 pivoting in block Gaussian Elimination
pivoting
Example 2
and
where
Question what is criterion to choose 1x1, 2x2 or
3x3 pivoting?
stability
performance
Question what is the cost to estimate the
criterion ?
6
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
7
LU-factorization for real symmetry indefinite
matrix A
factorization
factorization
where
and
Question why not LU-decomposition?
A is real symmetric, we only store lower triangle
part of A, say
8
Pointer array to store lower triangular part of
A
col-major based
row-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?
9
OutLine

• 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

10
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
1
2
3
1
3
2
interchange row 2, 3
interchange column 2, 3
2
4
5
3
5
6
3
6
5
3
5
6
2
4
5
2
5
4
11
Change diagonal element to (1,1) position
1
suppose we want to change
to position (1,1) , then consider permutation
and do symmetric permutation on A, say
9
8
6
3
1
2
3
4
3
2
1
4
6
5
2
7
7
7
2
5
6
6
5
2
interchange col 1, 3
interchange row 1, 3
3
2
1
4
3
6
8
9
8
6
3
9
10
9
7
4
4
7
9
10
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
7
2
5
6
7
2
5
6
9
3
6
8
9
3
6
1
4
7
9
10
4
7
9
10
12
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
2
3
6
1
9
6
10
4
7
9
4
7
9
does not changed
diagonal term
we have changed
does not changed
13
Change diagonal element to (1,1) position
3
6
2
2
6
2
2
6
interchange col 1, 3
6
2
interchange row 1, 3
2
6
6
9
7
4
4
7
9
9
7
4
Keep lower triangle part
8
8
6
Fill-in diagonal term
6
5
6
5
Fill-in A(3,1)
2
2
1
3
2
1
9
7
4
9
7
4
10
9
7
4
10
Question any compact form to represent above
interchange?
14
Change diagonal element to (1,1) position
4
2
6
2
6
6
2
4
7
9
9
7
4
4
7
9
case 1
4
7
9
case 2
2
6
15
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.
does not changed
keep position (k,1) and (1,k), i.e.
16
Change diagonal element to (1,1) position
6
Step 2 interchange column 1 and k below row k1,
i.e
Therefore
17
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
18
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
19
Implementation of symmetric permutation swapping
overhead 2
Step 3 interchanging column 1 and row k, i.e.
is NOT efficient since
is NOT contiguous.
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.
20
Change off-diagonal element to (2,1) position
1
suppose we want to change
to position (2,1), then consider first
permutation
which change first index from 4 to 2
2
1
4
3
1
2
3
4
1
4
3
2
4
10
9
7
7
5
2
5
6
2
7
6
interchange col 2, 4
interchange row 2, 4
3
9
8
6
3
6
8
9
3
9
8
6
5
2
7
6
4
7
9
10
4
10
9
7
and second permutation
which change second index from 3 to 1
1
4
3
2
3
4
1
2
8
9
3
6
4
10
9
7
9
10
4
7
interchange row 1, 3
9
10
4
7
interchange col 1, 3
3
9
8
6
8
9
3
6
3
4
1
2
5
5
5
2
7
6
6
7
2
6
7
2
21
Change off-diagonal element to (2,1) position
2
1
2
3
4
1
4
3
2
8
9
3
6
7
7
7
2
5
6
4
10
9
9
10
4
interchange row and col 2, 4
interchange row and col 1, 3
3
6
8
9
3
9
8
6
3
4
1
2
4
7
9
10
2
7
6
5
6
7
2
5
Observation two permutation matrices can be
computed easily since 2 ?? 4 and
1 ?? 3 are independent.
22
Change off-diagonal element to (2,1) position
3
HOW to transform
to position (2,1), under only lower triangle part
of A is stored?
Observation 1
and
does not changed
4
1
2
3
1
4
3
2
2
5
6
7
4
10
9
7
interchange row/col 2, 4
3
6
8
9
3
9
8
6
10
5
4
7
9
2
7
6
Step 1 interchange position (4,4) and (2,2),
4
1
2
3
1
2
3
4
2
10
6
7
2
5
6
7
3
6
8
9
3
6
8
9
4
7
9
5
4
7
9
10
23
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 ).
1
2
3
4
2
10
6
7
2
9
3
6
8
3
6
4
7
9
5
4
9
Step 2 interchange row 2 and 4 below col 2, i.e
2
4
9
interchange col 2, 4
4
9
interchange row 2, 4
3
6
3
6
3
9
6
4
9
2
6
2
6
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.
24
Change off-diagonal element to (2,1) position
5
2
4
3
6
3
6
4
9
2
9
Step 3 interchange column 2 (from row-3 to
row-3) and row 4 (from col-3 to col-3)
4
4
3
6
3
9
2
9
2
6
25
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.
does not changed
keep position (k,2) and (2,k), i.e.
26
Change off-diagonal element to (2,1) position
7
Step 2 interchange column 2 and k below row k1,
i.e
Therefore
27
Change off-diagonal element to (2,1) position
8
Step 3 interchange row 2 (column 1) and row k
(column1), i.e
Therefore
28
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
29
OutLine

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

30
LDL decomposition 1x1, 2x2 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

31
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
with
and
32
Growth rate in LDL decomposition
2
Therefore
such that
Observation in PALU, we choose
implies
Therefore
Observation large element growth will not occur
for a
pivot if
is large relative to
33
Growth rate in LDL decomposition
3
Case 2 s 2
where
and
with
and
with
and
with
implies
34
Growth rate in LDL decomposition
4
for
Therefore
35
Growth rate in LDL decomposition
5
Therefore
36
Criterion for pivot strategy 1
Fix a number
( later on, we will determine optimal value of
)
maximal element of matrix A
maximal diagonal element of matrix A
Case 1
suppose
we interchange 1st and k-th rows and columns by
introducing
permutation matrix
and do symmetric permutation
Then do
on
with
pivot, written as
37
Criterion for pivot strategy 2
Recall for
pivot,
, we have growth rate
and
Now for
and
Case 2
, we interchange q-th and 1st rows and columns,
suppose
and then r-th and 2nd rows and columns by
introducing permutation matrix
and do symmetric permutation
Question we must change q??1 first, then change
r ??2, why?
38
Criterion for pivot strategy 3
transforms
1
2
(boundary case)
3
but
39
Criterion for pivot strategy 4
pivot
Then do
on
with
pivot, written as
Recall for
pivot,
, we have growth rate
and
40
Criterion for pivot strategy 5
Now for
To sum up
Case 1
and
Case 2
and
Question How to choose value of
41
Criterion for pivot strategy 6
worst case analysis choose
such that
growth rate of
pivot
pivot
growth rate of
pivot
or equivalently
growth rate of
pivot
growth rate of
pivot
Define
where
satisfies
and
Exercise verify
42
Diagonal pivoting versus complete pivoting of
Gaussian Elimination 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
1
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
and column
then
with
43
Diagonal pivoting versus complete pivoting of
Gaussian Elimination 2
growth rate
2
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
44
OutLine

• Preliminary
• Symmetric permutation
• LDL decomposition (diagonal pivoting)
• Example of complete diagonal pivoting
• Algorithm of complete diagonal pivoting

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



-8



1


-13
-7
3
12
-8

12
6

3
12
3
6
-6
3
-13
-7

-13
-7

4
6
4
1
-6
6
-6
4
1
6
-6
1
46
Example (complete pivoting) 2
-8
-13
12
4
1
-8
-13

-13
-7
3
1
1
-13
-7

12
3
6
-6
0.3982
-1.1681
1


4.7257
-6.4248
4
1
-6
6
0.1327
-0.3894
1


-6.4248
5.8584
Recursively do the same procedure for
we choose
pivot
such that
swap row/column
by permutation matrix
47
Example (complete pivoting) 3
Do symmetry permutation for
and
-8
-13

-8
-13

-13
-7

-13
-7



4.7257
-6.4248


5.8584
-6.4248


-6.4248
5.8584


-6.4248
4.7257
1
1
1
1
0.3982
-1.1681
1
0.1327
-0.3894
1
0.1327
-0.3894
1
0.3982
-1.1681
1
5.8584
-6.4248
we do
pivot on
-6.4248
4.7257
48
Example (complete pivoting) 4
5.8584
-6.4248
1
5.8584
1
-1.0967
-6.4248
4.7257
-1.0967
1
-2.3202
1
Or write in original matrix form
-8
-13

1

-8
-13

-13
-7

1
-13
-7



5.8584
-6.4248
1


5.8584


-6.4248
4.7257
-1.0967
1


-2.3202
49
Example (complete pivoting) 5
1


-8
-13


1

-13
-7

0.1327
-0.3984
1


5.8584
0.3982
-1.1681
-1.0967
1


-2.3202
In practice, we have two key issues
col-major based
we only store lower part of matrix
1
6
6
6
-8
-7
12
-8
12
-13
1
3
-13
-7
3
4
-6
4
1
6
-6
50
Example (complete pivoting) 6
we store decomposition
and
into original
2
1
-8

1
-13
-7
0.1327
-0.3984
1


5.8584
0.3982
-1.1681
-1.0967
1


-2.3202
-8
-13
-7
0.1327
-0.3984
5.8584
0.3982
-1.1681
-1.0967
-2.3202
Question How can you judge correct decomposition
from
51
Example (complete pivoting) 7
Case 1 four 1x1 pivot
1
-8
13
1
-7
0.1327
-0.3984
1


5.8584
0.3982
-1.1681
-1.0967
1


-2.3202
Case 2 two 2x2 pivot
1
-8

1
-13
-7
0.1327
-0.3984
1


5.8584
0.3982
-1.1681
1


-1.0967
-2.3202
52
Example (complete pivoting) 8
Case 3 1x1 pivot 2x2 pivot 1x1 pivot
-8
1

-7
-13
1

-0.3984
5.8584
0.1327
1


-2.3202
0.3982
-1.1681
-1.0967
1
Case 4 2x2 pivot 1x1 pivot 1x1 pivot
1
-8

1
-13
-7
0.1327
-0.3984
1


5.8584
0.3982
-1.1681
-1.0967
1


-2.3202
Solution we need an array to record pivot
sequence
1x1 pivot
pivot
2
0
1
1
2x2 pivot
53
OutLine

• Preliminary
• Symmetric permutation
• LDL decomposition (diagonal pivoting)
• Example of complete diagonal pivoting
• Algorithm of complete diagonal pivoting

54
Algorithm ( PAP LDL ) 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
55
Algorithm ( PAP LDL ) 2
compute
and
1
Case 1
define permutation
to do symmetric permutation
2
To compute
, we only update lower triangle of
1
1
3
2
3
3
then
2
56
Algorithm ( PAP LDL ) 3
To compute
We only update lower triangle matrix
4
then
do 1x1 pivot
5
then
and
6
57
Algorithm ( PAP LDL ) 4
Case 2
define permutation
to interchange q-th and k-th rows and columns,
and then r-th and (k1)-th rows and columns
7
and then
To compute
, we only update lower triangle of
(1) do interchange row/col
1
1
8
2
3
3
then
2
58
Algorithm ( PAP LDL ) 5
(2) do interchange row/col
9
1
2
3
4
1
3
4
then
where
2
59
Algorithm ( PAP LDL ) 6
To compute
10
(1) do interchange row
11
(2) do interchange row
then
60
Algorithm ( PAP LDL ) 7
do 2x2 pivot
12
then
and
13
means
is 2x2 pivot
endwhile
61
Question recursion implementation

• Initialization- check algorithm holds for k1
• Recursion formula- check algorithm holds for k
  or k1, if k-1 is true
• Termination condition- check algorithm holds for
  kn-1
• Breakdown of algorithm- check which condition
  PAPLDL fails
• No extension algorithm works only for square,
  symmetric indefinite matrix.

normal
abnormal
Extension of algorithm
62
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
6
12
3
-6
12
-8
-13
4
2
3
4
1
return
1
3
-13
-7
2
0
1
1
-6
4
1
6
1
-8

1
-13
-7
0.1327
-0.3984
1


5.8584
0.3982
-1.1681
-1.0967
1


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