Sparse Matrix Methods
Sparse matrices store and operate on only their nonzero entries, making it possible to solve systems with millions of unknowns.
Exploiting mostly-zero matrices
Discretizing a PDE couples each unknown to only a few neighbors, so the resulting matrix is sparse: the overwhelming majority of its entries are zero. Storing and computing with only the nonzeros turns intractable dense problems, which cost O(n^2) memory and O(n^3) work, into feasible ones.
Storage formats
- COO (coordinate): lists of row, column, value, simple to build.
- CSR/CSC (compressed sparse row/column): efficient for matrix-vector products and solves.
- Diagonal and banded formats: for structured stencils where nonzeros lie on a few diagonals.
Sparse direct solvers and fill-in
Sparse LU or Cholesky factorization can solve sparse systems directly, but factorization creates fill-in: new nonzeros where the original had zeros. Reordering the rows and columns, using algorithms such as approximate minimum degree or nested dissection, drastically reduces fill and is essential to sparse direct performance.
import numpy as np
from scipy.sparse import diags
from scipy.sparse.linalg import spsolve
n = 1000
A = diags([-1, 2, -1], [-1, 0, 1], shape=(n, n)).tocsr()
b = np.ones(n)
x = spsolve(A, b) # sparse direct solveWhen to go iterative
For very large three-dimensional problems, fill-in makes even reordered direct solvers too costly, and iterative Krylov methods with preconditioning become the only practical option. They need only sparse matrix-vector products, whose cost is proportional to the number of nonzeros, and can run matrix-free.
Sparse storage and preconditioned iterative solvers make the multi-million-unknown field and transport systems of breeder Hyperion simulations tractable on available hardware.