Transcript No Slide Title

Direct Methods for Sparse Linear
Systems
Lecture 4
Alessandra Nardi
Thanks to Prof. Jacob White, Suvranu De, Deepak Ramaswamy,
Michal Rewienski, and Karen Veroy
Last lecture review
• Solution of system of linear equations
• Existence and uniqueness review
• Gaussian elimination basics
– GE basics
– LU factorization
– Pivoting
Outline
• Error Mechanisms
• Sparse Matrices
– Why are they nice
– How do we store them
– How can we exploit and preserve sparsity
Error Mechanisms
• Round-off error
– Pivoting helps
• Ill conditioning (almost singular)
– Bad luck: property of the matrix
– Pivoting does not help
• Numerical Stability of Method
Ill-Conditioning : Norms
• Norms useful to discuss error in numerical
problems
• Norm
 : V  
(1)
x  0 if x  0, x  V
(2)
x   x if   , x  V
(3)
x  y  x  y if x, y  V
Ill-Conditioning : Vector Norms
L2 (Euclidean) norm :
x
2
1
n


i 1
L1 norm :
x
Unit circle

xi
i 1
2
1
1
n

x
2
xi
1
x
1
1
L norm :
x

 max xi
i
Unit square
x

1
Ill-Conditioning : Matrix Norms
A  max
Vector induced norm :
x
Ax
x
 max Ax
x 1
Induced norm of A is the maximum “magnification” of
A
A
A
1

2
 max
j
n

i 1
 max
i
xby A
Aij = max abs column sum
n

j 1
Aij = max abs row sum
= (largest eigenvalue of ATA)1/2
Ill-Conditioning : Matrix Norms
• More properties on the matrix norm:
I 1
AB  A B
• Condition Number:
k ( A)  A1 A
-It can be shown that: k ( A)  1
-Large k(A) means matrix is almost singular (ill-conditioned)
Ill-Conditioning: Perturbation Analysis
What happens if we perturb b? b  b  b  x  x  x
M ( x  x )  b  b
Mx  Mx  b  b
Mx  b
x  M 1b
x  M 1 b
x  M
x
x
 M
1
1
b
M
b
b
 M
b
b
1

b
b
x
x
M x
 k (M )
b
b
 k(M) large is bad
Ill-Conditioning: Perturbation Analysis
What happens if we perturb M? M  M  M  x  x  x
( M  M )( x  x)  b
x
M
1
 M M
x  x
M

x
M
 k (M )
x  x
M
 k(M) large is bad
Bottom line:
If matrix is ill-conditioned, round-off puts you in troubles
Ill-Conditioning: Perturbation Analysis
Geometric Approach is more intuitive
M  [ M1 M 2 ], Solving M x  b is finding x1 M1  x2 M 2  b
x2
M2
|| M 2 ||
x2
M1
|| M 1 ||
x1
Vectors are orthogonal
M2
|| M 2 ||
M1
|| M 1 ||
x1
Vectors are nearly aligned
When vectors are nearly aligned,
Hard to decide how much of M 1versus how much of M 2
Numerical Stability
• Rounding errors may accumulate and
propagate unstably in a bad algorithm.
• Can be proven that for Gaussian elimination
the accumulated error is bounded
Summary on Error Mechanisms for GE
• Rounding: due to machine finite precision we have an
error in the solution even if the algorithm is perfect
– Pivoting helps to reduce it
• Matrix conditioning
– If matrix is “good”, then complete pivoting solves any
round-off problem
– If matrix is “bad” (almost singular), then there is nothing to
do
• Numerical stability
– How rounding errors accumulate
– GE is stable
LU – Computational Complexity
• Computational Complexity
– O(n3), where M: n x n
• We cannot afford this complexity
• Exploit natural sparsity that occurs in circuits
equations
– Sparsity: many zero elements
– Matrix is sparse when it is advantageous to exploit
its sparsity
• Exploiting sparsity: O(n1.1) to O(n1.5)
LU – Goals of exploiting sparsity
(1) Avoid storing zero entries
–
–
Memory usage reduction
Decomposition is faster since you do need to access
them (but more complicated data structure)
(2) Avoid trivial operations
–
–
Multiplication by zero
Addition with zero
(3) Avoid losing sparsity
Sparse Matrices – Resistor Line
1
2
m
3
X X
X X X


X X X

X X


X







m 1
4
X
X
X
X
X
X










X

X X
X X 
m
Tridiagonal Case
GE Algorithm – Tridiagonal Example
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!
Sparse Matrices – Fill-in – Example 1
R4
V3
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
  e1   0 
 e    0 
 2   
 e3  is1 
Symmetric
Diagonally Dominant
Sparse Matrices – Fill-in – Example 1
Matrix Non zero structure
Matrix after one LU step
X X X
X X 0 


 X 0 X 
X X X
X X X

0


0 X 
 X X
X= Non zero
Sparse Matrices – Fill-in – Example 2
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
Sparse Matrices – Fill-in & Reordering
V3
V1
V2
0
V3
V2
V1
0
x
x

 x
x
x
x

 0
x
x
x
x
x
x

x Fill-ins

x 
0

x No Fill-ins

x 
Node Reordering Can Reduce Fill-in
- Preserves Properties (Symmetry, Diagonal Dominance)
- Equivalent to swapping rows and columns
Exploiting and maintaining sparsity
• Criteria for exploiting sparsity:
– Minimum number of ops
– Minimum number of fill-ins
• Pivoting to maintain sparsity: NP-complete problem
 heuristics are used
– Markowitz, Berry, Hsieh and Ghausi, Nakhla and Singhal
and Vlach
– Choice: Markowitz
• 5% more fill-ins but
• Faster
• Pivoting for accuracy may conflict with pivoting for
sparsity
Sparse Matrices – Fill-in & Reordering
Where can fill-in occur ?
Already Factored
Multipliers








x
x
x
x
x
x

 Possible Fill-in

x  Locations

x
x

Fill-in Estimate = (Non zeros in unfactored part of Row -1)
 (Non zeros in unfactored part of Col -1)
Markowitz product

Sparse Matrices – Fill-in & Reordering
Markowitz Reordering (Diagonal Pivoting)
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 (but close to optimal) !
Sparse Matrices – Fill-in & Reordering
Why only try diagonals ?
• Corresponds to node reordering in Nodal formulation
1
2
0
3
1
0
3
• Reduces search cost
• Preserves Matrix Properties
- Diagonal Dominance
- Symmetry
2
Sparse Matrices – Fill-in & Reordering
Pattern of a Filled-in Matrix

























Very Sparse
Very Sparse
Dense

























Sparse Matrices – Fill-in & Reordering
Unfactored random Matrix
Sparse Matrices – Fill-in & Reordering
Factored random Matrix
Sparse Matrices – Data Structure
• Several ways of storing a sparse matrix in a
compact form
• Trade-off
– Storage amount
– Cost of data accessing and update procedures
• Efficient data structure: linked list
Sparse Matrices – Data Structure 1
Orthogonal linked list
Sparse Matrices – Data Structure 2
Vector of row
pointers
1
N
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
Val N1 Val N2
Val Nj
Col N1 Col N2
Col Nj
Sparse Matrices – Data Structure
Problem of Misses
Eliminating Source Row i from Target row j
Row i
Row j
M i ,i 1
i 1
M i ,i  7
i7
M i ,i 15
i  15
M i ,i 1
i 1
M i ,i  4
i4
M i ,i  5
i5
M i ,i  7
i7
M i ,i  9
i9
M i ,i 12
i  12
M i ,i 15
i  15
Must read all the row j entries to find the 3 that match row i
Every Miss is an unneeded memory reference (expensive!!!)
Could have more misses than ops!
Sparse Matrices – Data Structure
Scattering for Miss Avoidance
Row j
M i ,i 1
M i ,i 1
i 1
M i ,i  4
i4
M i ,i  4 M i ,i  5
M i ,i  5
i5
M i ,i  7
M i ,i  7
i7
M i ,i  9
M i ,i  9
i9
M i ,i 12
i  12
M i ,i 12
M i ,i 15
i  15
M i ,i 15
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!
Sparse Matrices – Graph Approach
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
Sparse Matrices – Graph Approach
Markowitz Products
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
1
X
2
Markowitz Products =
3
4
5
(Node Degree)2
2
(degree
1)
9
M 11 3  3  9
2
(deg
ree
2)
4
M 22 2  2  4
2
(deg
ree
3)
9
M 33 3  3  9
2
(degree
4)
4
M 44 2  2  4
2
(degree
5)
4
M  22  4
Sparse Matrices – Graph Approach
Factorization
One Step of LU Factorization
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
Sparse Matrices – Graph Approach
Example
x
x



x

x
x
x
x
x
x
x
x
x
x
x
x
x


x

x
x

Markowitz products
( = Node degree)
1
2
3
4
Graph
1
2
3
4
5
Col Row
3
3
 9
2
2
 4
3
3
 9
3
3
3
3
 9
 9
5
Sparse Matrices – Graph Approach
Example
Swap 2 with 1
x
x

x




x
x
x
x
x
x
x
x
x
x
x
x
x
x

x

x

x
x

1
3
Graph
4
5
Summary
• Gaussian Elimination Error Mechanisms
– Ill-conditioning
– Numerical Stability
• Gaussian Elimination for Sparse Matrices
– Improved computational cost: factor in O(N1.5)
operations (dense is O(N3) )
– Example: Tridiagonal Matrix Factorization O(N)
– Data structure
– Markowitz Reordering to minimize fill-ins
– Graph Based Approach