Sparse and large linear systems

Discretize the Poisson equation on a 400 × 400 grid and you get 160,000 unknowns and a matrix that needs 205 GB stored dense and 10 MB stored sparse, because each row has at most five nonzero entries. Matrices like that come from every field computed on a mesh or a network: the electrostatic potential in a device, the groundwater pressure under a dam, the stresses in a bridge, a morphogen diffusing through an embryo, the Hückel Hamiltonian of a graphene flake. The question is which solver finishes.

In Python, scipy.sparse builds the matrix, from a stencil with diags and kron or from triplets of row, column, and value with coo_array, converted to CSR before any arithmetic. scipy.sparse.linalg solves it directly with spsolve and splu, iteratively with cg for symmetric positive definite matrices and gmres for the rest, preconditions with spilu, and finds a few eigenvalues with eigsh. In Julia, the standard library SparseArrays provides sparse and spdiagm, backslash calls a sparse direct solver, and Krylov.jl and IterativeSolvers.jl supply cg and gmres.

Direct solvers are the right first choice in two dimensions, but they fill in: the LU factors of the Poisson matrix on a 200 × 200 grid hold 17 times as many nonzeros as the matrix. Iterative solvers keep the memory flat and pay in iterations instead. Conjugate gradients need 304 iterations on a 100 × 100 grid, 605 on 200 × 200, and 1,156 on 400 × 400, about twice as many each time the grid spacing halves. Stopping that growth is the job of a preconditioner.

Start with storage formats and building a matrix from a stencil, then a direct solve, then conjugate gradients with and without a preconditioner. Dense matrices are linear algebra, and the meshes many of them come from are finite elements.

What belongs here

Linear algebra when the matrix is large and mostly zeros: sparse storage formats, building matrices from stencils and graphs, direct sparse solvers, iterative methods (conjugate gradients, GMRES), preconditioners, and sparse eigenvalue problems. Dense matrices and the concepts behind them belong to linear-algebra.