Transcript Lecture 6: Matrix Computations: LU Decomposition

Numerical Simulation
CSE245 Spring 2010
Solution of Sparse Linear Systems
Thanks to: Kin Sou, Deepak Ramaswamy, Michal Rewienski,
Jacob White, Shihhsien Kuo and Karen Veroy and
especially Luca Daniel
Outline of today’s lecture
• Solution of Sparse Linear Systems
– Examples of Problems with Sparse Matrices
• Struts and joints, resistor grids, 3-D heat flow
– Tridiagonal Matrix Factorization
– General Sparse Factorization
• Fill-in and Reordering
• Graph Based Approach
– Sparse Matrix Data Structures
• Scattering
2
Sparse Matrices
Applications
Space Frame
Nodal Matrix
Space Frame
5
3
4
2
1
7
6
9
8
X
X
X



X







X
X X
X
X X
X X
X
Unknowns : Joint positions
Equations :  forces = 0
X
X
X
X
X
X
X
X
X
X
X
X
X X
X
X X
X X
X
X X
X

X= 

X
X
X
X









X

X
X
X 



2 x 2 block
3
Applications
Sparse Matrices
1
m 1
2
3
m2
m3
(m  1) (m  1)
Resistor Grid
4
m 1
m
2m
m2
Unknowns : Node Voltages
Equations :  currents = 0
4
Sparse Matrices
Nodal Formulation
Applications
Resistor Grid
Matrix non-zero locations for 100 x 10 Resistor Grid
Sparse Matrices
Nodal Formulation
Applications
Temperature in a cube
Temperature known on surface, determine interior temperature
m2  1
Circuit
Model
1
m 1
2
m2
m2  2
Outline
• Solution of Dense Linear Systems
• Solution of Sparse Linear Systems
– LU Factorization Reminder.
– Example of Problems with Sparse Matrices
• Struts and joints, resistor grids, 3-d heat flow
– Tridiagonal Matrix Factorization
– General Sparse Factorization
• Fill-in and Reordering
• Graph Based Approach
– Sparse Matrix Data Structures
• Scattering
7
Sparse Matrices
Nodal Formulation
1
2
m
Tridiagonal Example
3
X X
X X X


X X X

X X


X







Matrix Form
m 1
4
X
X
X
X
X
X










X

X X
X X 
m
How many
operations to
do LU
factorization?
Tridiagonal Example
Sparse Matrices
GE Algorithm
For i = 1 to n-1 {
For j = i+1 to n {
M ji 
M ji
M ii
“For each Row”
“For each target Row below the source”
Pivot
For k = i+1 to n { “For each Row element beyond Pivot”
M jk  M jk  M ji M ik
}
Multiplier
Order n Operations!
}
}
9
Outline
• Solution of Dense Linear Systems
• Solution of Sparse Linear Systems
– LU Factorization Reminder.
– Example of Problems with Sparse Matrices
• Struts and joints, resistor grids, 3-d heat flow
– Tridiagonal Matrix Factorization
– General Sparse Factorization
• Fill-in and Reordering
• Graph Based Approach
– Sparse Matrix Data Structures
• Scattering
10
Sparse Matrix Fill-In
Sparse Matrices
R4
Example
V3
Resistor Example
V2
V1
R1
R2
R5
R3
iS 1
Nodal Matrix
0




1 1
1
1
1




R1 R2 R4
R2
R4
1
1 1


R2
R2 R3
1
1 1


R4
R4 R5
 V1  0 
Symmetric
 V   0 
  2    Diagonally Dominant
 V3  iS 1 
11
Sparse Matrices
Matrix Non zero structure
X X X
X X 0 


 X 0 X 
Fill-In
Example
Fill Ins
X X X
X X X

0


0 X 
 X X
X= Non zero
X = fill in where 0 became non-zero
12
Sparse Matrices
Fill-In
Second Example
Fill-ins Propagate
X
X
X

0

0
X
X
X
X
X0
X
X
X
X
X
X
X0
X

X0

X0 

X0 
Fill-ins from Step 1 result in Fill-ins in step 2
13
Sparse Matrices
V3
V1
V2
0
V3
V2
V1
Fill-In
Reordering
x
x

 x
x
x
x
x

 0
x
x
x
x

x Fill-ins

x 
0

x No Fill-ins

x 
x
0
Node Reordering
- Can reduce fill-in
- Preserves properties (Symmetry, Diagonal Dominance)
- Equivalent to swapping rows and columns
14
Sparse Matrices
Fill-In
Reordering
Where can fill-in occur ?


Possible
Fill-in

x 0 x 0 x Locations


xx x0 x 
x0 xx x 

Already Factored
Multipliers








Fill-in Estimate = (Non zeros in unfactored part of Row -i)
 (Non zeros in unfactored part of Col -i)
Markowitz product

15
Sparse Matrices
Fill-In
Reordering
Markowitz Reordering
For i  1 to n
Find diagonal j  i with min Markowitz Product
Swap rows j  i and columns j  i
Factor the new row i and determine fill-ins
End
Greedy Algorithm !
16
Fill-In
Sparse Matrices
Reordering
Why only try diagonals ?
• Corresponds to node reordering in Nodal formulation
1
2
0
3
1
0
3
2
• Reduces search cost
• Preserves Matrix Properties
- Diagonal Dominance
- Symmetry
17
Fill-In
Sparse Matrices
Pattern of a Filled-in Matrix

























Very Sparse
Very Sparse
Dense

























18
Sparse Matrices
Fill-In
Unfactored Random Matrix
Symmetric
19
Sparse Matrices
Fill-In
Factored Random Matrix
Symmetric
20
Outline
• Solution of Dense Linear Systems
• Solution of Sparse Linear Systems
– LU Factorization Reminder.
– Example of Problems with Sparse Matrices
• Struts and joints, resistor grids, 3-d heat flow
– Tridiagonal Matrix Factorization
– General Sparse Factorization
• Fill-in and Reordering
• Graph Based Approach
– Sparse Matrix Data Structures
• Scattering
21
Matrix Graphs
Sparse Matrices
Construction
Structurally Symmetric Matrices and Graphs
X
X X
X
X X
X X
X
X
1
2
X
3
X X X
X
X X
4
5
• One Node Per Matrix Row
• One Edge Per Off-diagonal Pair
22
Sparse Matrices
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
Matrix Graphs
Markowitz Products
X
1
2
3
4
5
Markowitz Products = (Node Degree)2
M11  3333  9
M 22  22 22  4
M 33  3
3 33  9
M 44  22 22  4
M 55  22 2  4
(degree 1)2  9
(deg ree 2)2  4
(deg ree 3)2  9
(degree 4)2  4
(degree 5)2  4
23
Matrix Graphs
Sparse Matrices
Factorization
One Step of LU Factorization
X
X X
X
X
X X
X
X X
X
X
1
2
X
3
X X X
X
X
X
X X
4
5
• Delete the node associated with pivot row
• “Tie together” the graph edges
24
Matrix Graphs
Sparse Matrices
x
x



x

x
x
x
x
x
x
x
x
x
x
x
x
x


x

x
x

Example
1
2
3
4
5
Graph
Col Row
Markowitz products
( = Node degree)
1
3
3
 9
2
2
3
2
3
 4
 9
3
3
3
3
 9
 9
3
4
5
25
Matrix Graphs
Sparse Matrices
Example
Swap 2 with 1
x
x

x




x
x
0
X
x
0
X
x
x
x
x
x
x
x
x
x
Eliminate 1

x

x

x
x

12
2
12
3
4
5
3
4
5
Graph
26
Graphs
Sparse Matrices
1
m 1
2
3
m2
m3
(m  1) (m  1)
Resistor Grid Example
4
m 1
m
2m
m2
Unknowns : Node Voltages
Equations :  currents = 0
27
Sparse Matrices
Matrix Graphs
Grid Example
28
What should you pivot for?
Sparse Matrices
For growth control?
Or to avoid fill-ins?
A) LU factorization applied to a strictly diagonally
dominant matrix will never produce a zero pivot
B) The matrix entries produced by LU factorization
applied to a strictly diagonally dominant matrix will
never increase by more than a factor 2(n-1)
[which is anyway the best you can do by pivoting for growth control]
Bottom line:
if your matrix is strictly diagonally dominant
no need for numerical pivot for growth control!
29
Sparse Factorization Approach
1) Assume matrix requires NO numerical pivoting.
Diagonally dominant or symmetric positive definite.
2) Pre-Process:
•
Use Graphs to Determine Matrix Ordering
•
•
Many graph manipulation tricks used.
Form Data Structure for Storing Filled-in Matrix
•
Lots of additional non-zero added
3) Put numerical values in Data Structure and factor
Computation must be organized carefully!
Summary of Sparse Systems
• LU Factorization and Diagonal Dominance.
– Factor without numerical pivoting
• Sparse Matrices
– Struts, resistor grids, 3-d heat flow -> O(N) non-zeros
• Tridiagonal Matrix Factorization
– Factor in O(N) operations
• General Sparse Factorization
– Markowitz Reordering to minimize fill
• Graph Based Approach
– Factorization and Fill-in
– Useful for estimating Sparse GE complexity
Outline
• Solution of Sparse Linear Systems
– Example of Problems with Sparse Matrices
– Tridiagonal Matrix Factorization
– General Sparse Factorization
• Fill-in and Reordering
• Graph Based Approach
• Summary of the Algorithm
– Sparse Matrix Data Structures
• Scattering and symbolic fill in
32
Sparse Matrices
Vector of row
pointers
1
Sparse Data Structure
Arrays of Data in a Row
Val 11 Val 12
Val 1K
Col 11 Col 12
Col 1K
Matrix entries
Column index
Val 21 Val 22
Val 2L
Col 21 Col 22
Col 2L
Row pointers
Column Indices
N
Val N1 Val N2
Val Nj
Col N1 Col N2
Col Nj
Goal:
Never store a 0
Never multiply by 0
33
Sparse Matrices
Sparse Data Structure
Problem of Misses
Eliminating Source Row i from Target row j
Row i
Row j
M i ,k 1
M i ,k  7
M i ,k 15
k 1
k 7
k  15
M j ,k 1
M j ,k  4
M j , k 5
M j , k 7
M j ,k 9
M j ,k 12
M j ,k 15
k 1
k 4
k 5
k 7
k 9
k  12
k  15
Must read all the row j entries to find the 3 that match row i
34
Sparse Matrices
Sparse Data Structure
Data on Misses
Rows
Ops
Misses
Res
300
904387 248967
RAM
2806
1017289 3817587
Grid
4356
3180726 3597746
More
misses
than
ops!
Every Miss is an unneeded Memory Reference!
35
Sparse Matrices
Sparse Data Structure
Scattering for Miss Avoidance
Row j
M j ,k 1
M j ,k 1
M j ,k  4
M j , k 5
M j , k 7
M j ,k 9
M j ,k 12
M j ,k 15
k 1
k 4
k 5
k 7
k 9
k  12
k  15
M j ,k  4 M j , k 5
M j ,k  7
M j , k 9
M j ,k 12
M j ,k 15
Use target row approach – Row j is the target.
1) Read all the elements in Row j, and scatter them in an n-length vector
2) Access only the needed elements using array indexing!
36
Sparse Matrices – Another Data
Orthogonal linked list
Structure
But if fill in occurs, pointer structures change. What do we do?
Pre-compute pivoting order and sparse fill in structure symbolically
37
Summary of Sparse Systems
• LU Factorization and Diagonal Dominance.
– Factor without numerical pivoting
– hard problems (ill-conditioned)
• Sparse Matrices
– Struts, resistor grids, 3-d heat flow -> O(N) nonzeros
• Tridiagonal Matrix Factorization
– Factor in O(N) operations
• General Sparse Factorization
– Markowitz Reordering to minimize fill
• Graph Based Approach
– Factorization and Fill-in
– Useful for estimating Sparse GE complexity
• Sparse Data Structures
– Scattering and symbolic fill in
38