Direct Methods for Sparse Linear Systems

1. Direct Methods for Sparse Linear Systems Lecture 4 Alessandra Nardi Thanks to Prof. Jacob White, Suvranu De, Deepak Ramaswamy, Michal Rewienski, and Karen Veroy

2. Last lecture review • Solution of system of linear equations • Existence and uniqueness review • Gaussian elimination basics • GE basics • LU factorization • Pivoting

3. Outline • Error Mechanisms • Sparse Matrices • Why are they nice • How do we store them • How can we exploit and preserve sparsity

4. 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

5. Ill-Conditioning : Norms • Norms useful to discuss error in numerical problems • Norm

6. Ill-Conditioning : Vector Norms L2 (Euclidean) norm : Unit circle L1 norm : 1 1 L norm : Unit square

7. Vector induced norm : Ill-Conditioning : Matrix Norms Induced norm of A is the maximum “magnification” of by = max abs column sum = max abs row sum = (largest eigenvalue of ATA)1/2

8. Ill-Conditioning : Matrix Norms • More properties on the matrix norm: • Condition Number: • It can be shown that: • Large k(A) means matrix is almost singular (ill-conditioned)

9. Ill-Conditioning: Perturbation Analysis What happens if we perturb b?  k(M) large is bad

10. Ill-Conditioning: Perturbation Analysis What happens if we perturb M?  k(M) large is bad Bottom line: If matrix is ill-conditioned, round-off puts you in troubles

11. Vectors are orthogonal Vectors are nearly aligned Ill-Conditioning: Perturbation Analysis Geometric Approach is more intuitive When vectors are nearly aligned, Hard to decide how much of versus how much of

12. 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

13. 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

14. 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)

15. LU – Goals of exploiting sparsity • 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

16. Sparse Matrices – Resistor Line Tridiagonal Case m

17. Pivot Multiplier GE Algorithm – Tridiagonal Example For i = 1 to n-1 {“For each Row” For j = i+1 to n {“For each target Row below the source” For k = i+1 to n {“For each Row element beyond Pivot” } } } Order N Operations!

18. Sparse Matrices – Fill-in – Example 1 Nodal Matrix 0 Symmetric Diagonally Dominant

19. X X Sparse Matrices – Fill-in – Example 1 Matrix Non zero structure Matrix after one LU step X X X X X= Non zero

20. X X X X X Sparse Matrices – Fill-in – Example 2 Fill-ins Propagate X X X X X Fill-ins from Step 1 result in Fill-ins in step 2

21. Fill-ins 0 No Fill-ins 0 Sparse Matrices – Fill-in & Reordering Node Reordering Can Reduce Fill-in - Preserves Properties (Symmetry, Diagonal Dominance) - Equivalent to swapping rows and columns

22. 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

23. 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 Where can fill-in occur ? Already Factored Possible Fill-in Locations Multipliers

24. Sparse Matrices – Fill-in & Reordering Markowitz Reordering (Diagonal Pivoting) Greedy Algorithm (but close to optimal) !

25. Sparse Matrices – Fill-in & Reordering Why only try diagonals ? • Corresponds to node reordering in Nodal formulation 1 2 3 1 0 0 3 2 • Reduces search cost • Preserves Matrix Properties - Diagonal Dominance - Symmetry

26. Sparse Matrices – Fill-in & Reordering Pattern of a Filled-in Matrix Very Sparse Very Sparse Dense

27. Sparse Matrices – Fill-in & Reordering Unfactored random Matrix

28. Sparse Matrices – Fill-in & Reordering Factored random Matrix

29. 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

30. Sparse Matrices – Data Structure 1 Orthogonal linked list

31. Sparse Matrices – Data Structure 2 Vector of row pointers Arrays of Data in a Row Matrix entries Val 11 Val 12 Val 1K 1 Column index Col 11 Col 12 Col 1K Val 21 Val 22 Val 2L Col 21 Col 22 Col 2L Val N1 Val N2 Val Nj N Col N1 Col N2 Col Nj

32. Sparse Matrices – Data StructureProblem of Misses Eliminating Source Row i from Target row j Row i Row j 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!

33. Sparse Matrices – Data Structure Scattering for Miss Avoidance Row j 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!

34. Sparse Matrices – Graph Approach Structurally Symmetric Matrices and Graphs 1 X X X X X X X 2 X X X X 3 X X X 4 X X X 5 • One Node Per Matrix Row • One Edge Per Off-diagonal Pair

35. 1 2 X X X X 3 X X X 4 X X X X 5 X X X X X X Sparse Matrices – Graph ApproachMarkowitz Products Markowitz Products = (Node Degree)2

36. Sparse Matrices – Graph ApproachFactorization One Step of LU Factorization 1 X X X X X X X 2 X X X X X 3 X X X 4 X X X X X 5 • Delete the node associated with pivot row • “Tie together” the graph edges

37. 1 2 3 4 5 1 2 3 4 5 Sparse Matrices – Graph ApproachExample Graph Markowitz products ( = Node degree)

38. 1 3 4 5 Sparse Matrices – Graph ApproachExample Swap 2 with 1 Graph

39. 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