Sparse Matrix Methods

1
Sparse Matrix Methods

• Day 1 Overview
• Matlab and examples
• Data structures
• Axb
• Sparse matrices and graphs
• Fill-reducing matrix permutations
• Matching and block triangular form
• Day 2 Direct methods
• Day 3 Iterative methods

D
2
(Demo in Matlab)

• Matlab sprand
• spy
• sparse, full
• matrix arithmetic and indexing
• examples of sparse matrices from different
  applications (from UF site)

3
Matlab sparse matrices Design principles

• All operations should give the same results for
  sparse and full matrices (almost all)
• Sparse matrices are never created automatically,
  but once created they propagate
• Performance is important -- but usability,
  simplicity, completeness, and robustness are more
  important
• Storage for a sparse matrix should be O(nonzeros)
• Time for a sparse operation should be
  O(flops)(as nearly as possible)

4
Data structures

• Full
• 2-dimensional array of real or complex numbers
• (nrowsncols) memory

• Sparse
• compressed column storage
• about (1.5nzs .5ncols) memory

D
5
(Demos in Matlab)

• A \ b
• Cholesky factorization and orderings

6
Solving linear equations

• x A \ b
• Is A square?
• no gt use QR to solve least squares problem
• Is A triangular or permuted triangular?
• yes gt sparse triangular solve
• Is A symmetric with positive diagonal elements?
• yes gt attempt Cholesky after symmetric minimum
  degree
• Otherwise
• gt use LU on A(, colamd(A))

7
Matrix factorizations in Matlab

• Cholesky
• R chol(A)
• simple left-looking column algorithm
• Nonsymmetric LU
• L,U,P lu(A)
• left-looking GPMOD, depth-first search,
  symmetric pruning
• Orthogonal
• Q,R qr(A)
• George-Heath algorithm row-wise Givens rotations

8
Graphs and Sparse Matrices Cholesky
factorization
Fill new nonzeros in factor
Symmetric Gaussian elimination for j 1 to n
add edges between js higher-numbered
neighbors
G(A)chordal
G(A)
9
Elimination Tree
G(A)
T(A)
Cholesky factor

• T(A) parent(j) min i gt j (i,j) in G(A)
  
• T describes dependencies among columns of factor
• Can compute T from G(A) in almost linear time
• Can compute G(A) easily from T

D
10
(Demos in Matlab)

• matrix and graph
• elimination tree
• orderings in detail

11
Fill-reducing matrix permutations

• Minimum degree
• Eliminate row/col with fewest nzs, add fill,
  repeat
• Theory can be suboptimal even on 2D model
  problem
• Practice often wins for medium-sized problems
• Nested dissection
• Find a separator, number it last, proceed
  recursively
• Theory approx optimal separators gt approx
  optimal fill and flop count
• Practice often wins for very large problems
• Banded orderings (Reverse Cuthill-McKee, Sloan, .
  . .)
• Try to keep all nonzeros close to the diagonal
• Theory, practice often wins for long, thin
  problems
• Best modern general-purpose orderings are ND/MD
  hybrids.

12
Fill-reducing permutations in Matlab

• Nonsymmetric approximate minimum degree
• p colamd(A)
• column permutation lu(A(,p)) often sparser
  than lu(A)
• also for QR factorization
• Symmetric approximate minimum degree
• p symamd(A)
• symmetric permutation chol(A(p,p)) often
  sparser than chol(A)
• Reverse Cuthill-McKee
• p symrcm(A)
• A(p,p) often has smaller bandwidth than A
• similar to Sparspak RCM

D
13
(Demos in Matlab)

• nonsymmetric LU
• dmperm, dmspy, components

14
Matching and block triangular form

• Dulmage-Mendelsohn decomposition
• Bipartite matching followed by strongly connected
  components
• Square, full rank A
• p, q, r dmperm(A)
• A(p,q) has nonzero diagonal and is in block upper
  triangular form
• also, strongly connected components of a directed
  graph
• also, connected components of an undirected graph
• Arbitrary A
• p, q, r, s dmperm(A)
• maximum-size matching in a bipartite graph
• minimum-size vertex cover in a bipartite graph
• decomposition into strong Hall blocks

15
Complexity of direct methods
n1/2
n1/3
16
The Landscape of Sparse Axb Solvers
D
17
(Demos in Matlab)

• conjugate gradients