CS 5220

Sparse linear algebra

David Bindel

2024-10-31

Goal

Solve \[ Ax = b, \] where \(A\) is sparse (or data sparse).

Plan for today

  • Reminder of stationary iterations
  • Krylov idea and performance via
    • Parallelism in algorithm
    • Better convergence (preconditioning)
    • Better memory access (reordering)
  • Sparse Gaussian elimination

Reminder: Stationary Iterations

From splitting \(A = M-K\), compute: \[ x^{(k+1)} = x^{(k)} + M^{-1}(b-Ax^{(k)}). \] “Linear” rate of convergence, dependent on \(\rho(M^{-1}K)\).

Krylov Subspace Methods

What if we only know how to multiply by \(A\)? About all you can do is keep multiplying! \[ K_k(A,b) = \operatorname{span}\left\{ b, A b, A^2 b, \ldots, A^{k-1} b \right\}. \] Gives surprisingly useful information!

Example: Conjugate Gradients

If \(A\) is symmetric and positive definite, \(Ax = b\) solves a minimization: \[\begin{aligned} \phi(x) &= \frac{1}{2} x^T A x - x^T b\\ \nabla \phi(x) &= Ax - b. \end{aligned}\] Idea: Minimize \(\phi(x)\) over \(K_k(A,b)\). Basis for the method of conjugate gradients

Example: GMRES

Idea: Minimize \(\|Ax-b\|^2\) over \(K_k(A,b)\). Yields Generalized Minimum RESidual (GMRES)

Convergence of Krylov Subspace Methods

  • KSPs are not stationary (no constant fixed-point iteration)
  • Convergence is surprisingly subtle!
  • CG convergence upper bound via condition number
    • Large condition number iff \(\phi(x)\) is narrow
    • True for Poisson and company

Convergence of Krylov Subspace Methods

  • Preconditioned problem \(M^{-1} A x = M^{-1} b\)
  • Whence \(M\)?
    • From a stationary method?
    • From a simpler/coarser discretization?
    • From approximate factorization?

PCG

r = b-A*x;
p = 0; beta = 0;
z = Msolve(r);
rho = dot(r, z);
for i=1:nsteps
    p = z + beta*p;
    q = A*p;
    alpha = rho/dot(p, q);
    x += alpha*p;
    r -= alpha*q;
    if norm(r) < tol, break; end
    z = Msolve(r);
    rho_prev = rho;
    rho = dot(r, z);
    beta = rho/rho_prev;
end

PCG parallel work

  • Solve with \(M\)
  • Product with \(A\)
  • Dot products and axpys

Pushing PCG

  • Rearrange if \(M = LL^T\) is available
  • Or build around “powers kernel”
    • Old “s-step” approach of Chronopoulos and Gear
    • CA-Krylov of Hoemmen, Carson, Demmel
    • Hard to keep stable

Pushing PCG

Two real application levers:

  • Better preconditioning
  • Faster matvecs

PCG bottlenecks

Key: fast solve with \(M\), product with \(A\)

  • Some preconditioners parallelize better!
  • Balance speed with performance.
    • Speed for set up of \(M\)?
    • Speed to apply \(M\) after setup?
  • Cheaper to do two multiplies/solves at once...
    • Can’t exploit in obvious way — lose stability
    • Variants allow multiple products (CA-Krylov)
  • Lots of fiddling possible with \(M\); matvec with \(A\)?

Thinking on (basic) CG convergence

Consider 5-point stencil on an \(n \times n\) mesh.

  • Information moves one grid cell per matvec.
  • Cost per matvec is \(O(n^2)\).
  • At least \(O(n^3)\) work to get information across mesh!

Convergence by counting

  • Time to converge \(\geq\) time to move info across
  • For a 2D mesh: \(O(n)\) matvecs, \(O(n^3) = O(N^{3/2})\) cost
  • For a 3D mesh: \(O(n)\) matvecs, \(O(n^4) = O(N^{4/3})\) cost
  • “Long” meshes yield slow convergence

Convergence by counting

3D beats 2D because everything is closer!

  • Advice: sparse direct for 2D, CG for 3D.
  • Better advice: use a preconditioner!

Eigenvalue approach

Define the condition number for \(\kappa(L)\) s.p.d: \[\kappa(L) = \frac{\lambda_{\max}(L)}{\lambda_{\min}(L)}\] Describes how elongated the level surfaces of \(\phi\) are.

Eigenvalue approach

  • For Poisson, \(\kappa(L) = O(h^{-2})\)
  • Steps to halve error: \(O(\sqrt{\kappa}) = O(h^{-1})\).

Similar back-of-the-envelope estimates for some other PDEs. But these are not always that useful... can be pessimistic if there are only a few extreme eigenvalues.

Frequency-domain approach

FFT of e_0 FFT of e_{10}

Error \(e_k\) after \(k\) steps of CG gets smoother!

Preconditioning Poisson

  • CG already handles high-frequency error
  • Want something to deal with lower frequency!
  • Jacobi useless
    • Doesn’t even change Krylov subspace!

Preconditioning Poisson

Better idea: block Jacobi?

  • Q: How should things split up?
  • A: Minimize blocks across domain.
  • Compatible with minimizing communication!

Multiplicative Schwartz

Generalizes block Gauss-Seidel

Restrictive Additive Schwartz (RAS)

  • Get ghost cell data (green)
  • Solve everything local (including neighbor data)
  • Update local values for next step (local)
  • Default strategy in PETSc

Multilevel Ideas

  • RAS moves info one processor per step
  • For scalability, still need to get around this!
  • Basic idea: use multiple grids
    • Fine grid gives lots of work, kills high-freq error
    • Coarse grid cheaply gets info across mesh, kills low freq

Tuning matmul

Can also tune matrix multiply

  • Represented implicitly (regular grids)
    • Example: Optimizing stencil operations (Datta)
  • Or explicitly (e.g. compressed sparse column)
    • Sparse matrix blocking and reordering
    • Packages: Sparsity (Im), OSKI (Vuduc)
    • Available as PETSc extension

Or further rearrange algorithm (Hoemmen, Demmel).

Reminder: Compressed sparse row

for (int i = 0; i < n; ++i) {
  y[i] = 0;
  for (int jj = ptr[i]; jj < ptr[i+1]; ++jj)
    y[i] += A[jj]*x[col[jj]];
}

Where is the problem for memory access?

Memory traffic in CSR multiply

Memory access patterns:

  • Elements of \(y\) accessed sequentially
  • Elements of \(A\) accessed sequentially
  • Access to \(x\) are all over!

Can help by switching to block CSR. Switching to single precision, short indices can help memory traffic, too!

Parallelizing matvec

  • Each processor gets a piece
  • Many partitioning strategies
  • Idea: re-order so one of these strategies is “good”

Reordering for matvec

SpMV performance goals:

  • Balance load?
  • Balance storage?
  • Minimize communication?
  • Good cache re-use?

Reordering also comes up for GE!

Reminder: Sparsity and reordering

Permute unknowns for better SpMV or

  • Stability of Gauss elimination,
  • Fill reduction in Gaussian elimination,
  • Improved performance of preconditioners...

Reminder: Sparsity and partitioning

Want to partition sparse graphs so that

  • Subgraphs are same size (load balance)
  • Cut size is minimal (minimize communication)

Matrices that are “almost” diagonal are good?

Reordering for bandedness

Natural order RCM order

Reverse Cuthill-McKee

  • Select “peripheral” vertex \(v\)
  • Order according to breadth first search from \(v\)
  • Reverse ordering

From iterative to direct

  • RCM ordering is great for SpMV
  • But isn’t narrow banding good for solvers, too?
    • LU takes \(O(nb^2)\) where \(b\) is bandwidth.
    • Great if there’s an ordering where \(b\) is small!

Skylines and profiles

  • Profile solvers generalize band solvers
  • Skyline storage for lower triangle: for each row \(i\),
    • Start and end of storage for nonzeros in row.
    • Contiguous nonzero list up to main diagonal.
  • In each column, first nonzero defines a profile.
  • All fill-in confined to profile.
  • RCM is again a good ordering.

Beyond bandedness

  • Minimum bandwidth for 2D model problem? 3D?
  • Skyline only gets us so much farther

Beyond bandedness

But more general solvers have similar structure

  • Ordering (minimize fill)
  • Symbolic factorization (where will fill be?)
  • Numerical factorization (pivoting?)
  • … and triangular solves

Troublesome Trees

One step of Gaussian elimination completely fills this matrix!

Terrific Trees

Full Gaussian elimination generates no fill in this matrix!

Quiz

How many fill elements are there for elimination on \[A = \begin{bmatrix} x & x & 0 & 0 & x & 0 & 0 & x \\ x & x & x & 0 & 0 & 0 & 0 & 0 \\ 0 & x & x & x & 0 & 0 & 0 & 0 \\ 0 & 0 & x & x & x & 0 & 0 & 0 \\ x & 0 & 0 & x & x & x & 0 & 0 \\ 0 & 0 & 0 & 0 & x & x & x & 0 \\ 0 & 0 & 0 & 0 & 0 & x & x & x \\ x & 0 & 0 & 0 & 0 & 0 & x & x \end{bmatrix}\]

Graphic Elimination

Consider first steps of GE

A(2:end,1)     = A(2:end,1)/A(1,1);
A(2:end,2:end) = A(2:end,2:end)-...
                 A(2:end,1)*A(1,2:end);

Nonzero in the outer product at \((i,j)\) if A(i,1) and A(j,1) both nonzero — that is, if \(i\) and \(j\) are both connected to 1.

General: Eliminate variable, connect remaining neighbors.

Terrific Trees Redux

Order leaves to root \(\implies\)
on eliminating \(i\), parent of \(i\) is only remaining neighbor.

Nested Dissection

  • Idea: Think of block tree structures.
  • Eliminate block trees from bottom up.
  • Can recursively partition at leaves.

Nested Dissection

  • Rough cost estimate: how much just to factor dense Schur complements associated with separators?
  • Notice graph partitioning appears again!
    • And again we want small separators!

Nested Dissection

Model problem: Laplacian with 5 point stencil (for 2D)

  • ND gives optimal complexity in exact arithmetic
    (George 73, Hoffman/Martin/Rose)
  • 2D: \(O(N \log N)\) memory, \(O(N^{3/2})\) flops
  • 3D: \(O(N^{4/3})\) memory, \(O(N^2)\) flops

Minimum Degree

  • Locally greedy strategy
    • Want to minimize upper bound on fill-in
    • Fill \(\leq\) (degree in remaining graph)\(^2\)
  • At each step
    • Eliminate vertex with smallest degree
    • Update degrees of neighbors
  • Problem: Expensive to implement!
    • But better varients via quotient graphs
    • Variants often used in practice

Elimination Tree

  • Variables (columns) are nodes in trees
  • \(j\) a descendant of \(k\) if eliminating \(j\) updates \(k\)
  • Can eliminate disjoint subtrees in parallel!

Cache locality

Basic idea: exploit “supernodal” (dense) structures in factor

  • e.g. arising from elimination of separator Schur complements in ND
  • Other alternatives exist (multifrontal solvers)

Pivoting

Pivoting is painful, particularly in distributed memory!

  • Cholesky — no need to pivot!
  • Threshold pivoting — pivot when things look dangerous
  • Static pivoting — try to decide up front

What if things go wrong with threshold/static pivoting?
Common theme: Clean up sloppy solves with good residuals

Direct to iterative

Can improve solution by iterative refinement: \[\begin{aligned} PAQ &\approx LU \\ x_0 &\approx Q U^{-1} L^{-1} Pb \\ r_0 &= b-Ax_0 \\ x_1 &\approx x_0 + Q U^{-1} L^{-1} P r_0 \end{aligned} \] Looks like approximate Newton on \(F(x) = Ax-b = 0\).
This is just a stationary iterative method!
Nonstationary methods work, too.

Variations on a theme

If we’re willing to sacrifice some on factorization,

  • Single precision factor + double precision refinement?
  • Sloppy factorizations (marginal stability) + refinement?
  • Modify \(m\) small pivots as they’re encountered (low rank updates), fix with \(m\) steps of a Krylov solver?

Parting advice

  • Sparse direct for 2D problems
  • Gets more expensive for 3D problems
  • Approximate direct solves make good preconditioners