Graph Algorithms in Scientific Computing: Past, Present, and Future

1
Graph Algorithms in Scientific Computing
Past, Present, and Future

• John R. Gilbert
• MIT and UC Santa Barbara
• August 8, 2002

2
Outline

• Past (one extended example)
• Present (two quick examples)
• Future (6 or 8 examples)

3
A few themes

• Paths
• Locality
• Eigenvectors
• Huge data sets
• Multiple scales

4
PAST
5
Graphs and sparse Gaussian elimination (1961-)
Fill new nonzeros in factor
Cholesky factorization for j 1 to n add
edges between js higher-numbered neighbors
G(A)chordal
G(A)
6
Fill-reducing matrix permutations

• Theory approx optimal separators gt approx
  optimal fill and flop count
• Orderings nested dissection, minimum degree,
  hybrids
• Graph partitioning spectral, geometric,
  multilevel

7
Directed graph
A
G(A)

• A is square, unsymmetric, nonzero diagonal
• Edges from rows to columns
• Symmetric permutations PAPT

8
Symbolic Gaussian elimination
A
G (A)

• Add fill edge a -gt b if there is a path from a to
  b through lower-numbered vertices.

9
Strongly connected components
G(A)
PAPT

• Symmetric permutation to block triangular form
• Solve Axb by block back substitution
• Irreducible (strong Hall) diagonal blocks
• Row and column partitions are independent of
  choice of nonzero diagonal
• Find P in linear time by depth-first search

10
Sparse-sparse triangular solve
G(LT)
L
x
b

• Symbolic
• Predict structure of x by depth-first search from
  nonzeros of b
• Numeric
• Compute values of x in topological order

Time O(flops)
11
Column intersection graph
A
G?(A)
ATA

• G?(A) G(ATA) if no cancellation (otherwise
  ?)
• Permuting the rows of A does not change G?(A)

12
Filled column intersection graph
A
chol(ATA)

• G?(A) symbolic Cholesky factor of ATA
• In PALU, G(U) ? G?(A) and G(L) ? G?(A)
• Tighter bound on L from symbolic QR
• Bounds are best possible if A is strong Hall

13
Column elimination tree
T?(A)
A
chol(ATA)

• Elimination tree of ATA (if no cancellation)
• Depth-first spanning tree of G?(A)
• Represents column dependencies in various
  factorizations

14
That story continues . . .

• Matlab 4.0, 1992
• Left-looking column-by-column LU factorization
• Depth-first search to predict structure of each
  column
• Slow search limited speed
• BLAS-1 limited cache reuse
• . . .
• SuperLU supernodal BLAS-2.5 LU
• UMFPACK multifrontal BLAS-3 LU Matlab 6.5,
  2002
• Ordering for nonsymmetric LU is still not well
  understood

15
PRESENT
16
Support graph preconditioning

• Define a preconditioner B for matrix A
• Explicitly compute the factorization B LU
• Choose nonzero structure of B to make factoring
  cheap (using combinatorial tools from direct
  methods)
• Prove bounds on condition number using both
  algebraic and combinatorial tools
• Bruce Hendricksons talk

17
Support graph preconditioner example
Vaidya
G(A)
G(B)

• A is symmetric positive definite with negative
  off-diagonal nzs
• B is a maximum-weight spanning tree for A (with
  diagonal modified to preserve row sums)
• Preconditioning costs O(n) time per iteration
• Eigenvalue bounds from graph embedding product
  of congestion and dilation
• Condition number at most n2 independent of
  coefficients
• Many improvements exist

18
Support graph preconditioner example
Vaidya
G(A)
G(B)

• Improve congestion and dilation by adding a few
  strategically chosen edges to B
• Cost of factorsolve is O(n1.75), or O(n1.2) if A
  is planar

19
Link analysis of the world-wide web

• Web page vertex
• Link directed edge
• Link matrix Aij 1 if page i links to page j

20
Web graph PageRank (Google) Brin,
Page
An important page is one that many important
pages point to.

• Markov process follow a random link most of the
  time otherwise, go to any page at random.
• Importance stationary distribution of Markov
  process.
• Transition matrix is pA (1-p)ones(size(A)),
  scaled so each column sums to 1.
• Importance of page i is the i-th entry in the
  principal eigenvector of the transition matrix.
• But, the matrix is 2,000,000,000 by 2,000,000,000.

21
Web graph Hubs and authorities
Kleinberg
Authorities trekbikes.com shimano.com campagnol
o.com
Hubs bikereviews.com phillybikeclub.org yahoo.c
om/cycling
A good hub cites good authorities A good
authority is cited by good hubs

• Each page has hub score xi and authority score yi
• Repeat y ? Ax x ? ATy normalize
• Converges to principal eigenvectors of AAT and
  ATA(left and right singular vectors of A)

22
FUTURE
23
Combinatorial issues in numerical linear algebra

• Approximation theory for nonsymmetric LU ordering
  
• Preconditioning
• Complexity of matrix multiplication

24
Biology as an information science

• The rate-limiting step of genomics is
  computational science. - Eric
  Lander
• Sequence assembly, gene identification, alignment
• Large-scale data mining
• Protein modeling discrete and continuous

25
Linear and sublinear algorithms for huge problems

• Can we understand anything interesting about
  our data when we do not even have time to read
  all of it? - Ronitt
  Rubinfeld

26
Fast Monte Carlo algorithms for finding low-rank
approximations Frieze, Kannan,
Vempala

• Describe a rank-k matrix B0 that is within e of
  the best rank-k approximation to A
• A B0F ? minB A - BF e AF
• Correct with probability at least 1 d
• Time polynomial in k, 1/e, log(1/d) independent
  of size of A
• Idea using a clever distribution, sample an
  O(k-by-k) submatrix of A and compute its SVD
• (Need to be able to sample A with the right
  distribution)

27
Approximating the MST weight in sublinear time
Chazelle, Rubinfeld, Trevisan

• Key subroutine estimate number of connected
  components of a graph, in time depending on
  expected relative error but not on size of
  graph
• Idea for each vertex v define
• f(v) 1/(size of component containing v)
• Then Sv f(v) number of connected components
• Estimate Sv f(v) by breadth-first search from a
  few vertices

28
Approximate number of connected components

• for each of r randomly chosen vertices u1 ... ur
  do
• breadth-first search to distance 1 from ui
• repeat until bi is set
• if coin flip yields heads or BFS has reached at
  least w vertices
• then set bi 0 and break repeat
• continue breadth-first search to double
  vertices reached
• if BFS has reached all vertices in component
  containing ui
• then set bi 2(flips) / (vtxs reached) and
  break repeat
• end repeat
• end for
• return (n/r) sum (bi)

29
Managing data and computation on the web

• Mapping, balancing, scheduling for parallel
  computing
• Data placement, compression, communication
  optimization
• Sophisticated cacheing (Akamai)
• Internet backplane (UTK)
• Active buffering, compressed migration (UIUC)
• Shannon may become nearly as important as
  Wilkinson to scientific computing. -
  Rob Schreiber

30
Modeling and distributed control
Imagining a molecular-scale crossbar switch
Heath et al., UCLA

• Multiresolution modeling for nanomaterials and
  nanosystems
• Distributed control for systems of micromachines

31
Active surface airjet paper mover
Berlin, Biegelsen et al., PARC
PC-hosted DSP control _at_ 1 kHz
12 x 12 board
sensors 32,000 gray-level pixels in 25 linear
arrays
576 valves (144 per direction)
32
A hard question

• How will combinatorial methods be used by people
  who dont understand them in detail?

33
Matrix division in Matlab

• x A \ b
• Works for either full or sparse A
• 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))

34
Matching and depth-first search in Matlab

• dmperm 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

35
A few themes

• Paths
• Locality
• Eigenvectors
• Huge data sets
• Multiple scales
• Usability

36
Morals

• Things are clearer if you look at them from two
  directions
• Combinatorial algorithms are pervasive in
  scientific computing and will become more so
• What are the implications for teaching?
• What are the implications for software
  development?

37
Thanks

• Patrick Amestoy, Erik Boman, Iain Duff, Mary Ann
  Branch Freeman, Bruce Hendrickson, Esmond Ng,
  Alex Pothen, Padma Raghavan, Ronitt Rubinfeld,
  Rob Schreiber, Sivan Toledo, Paul Van Dooren,
  Steve Vavasis, Santosh Vempala, ...