Sparse Matrices
Matrices dominated by zeros, stored and solved with methods that touch only the nonzeros.
What sparsity means
A matrix is sparse when most of its entries are zero. Sparsity is not a curiosity but the norm for large problems: discretizing a physical domain couples each grid point only to its neighbors, so the resulting matrix has a handful of nonzeros per row regardless of size. Storing and computing only the nonzeros is what makes million-variable problems tractable.
Storage formats
- COO (coordinate): lists of row index, column index, and value, easy to build
- CSR (compressed sparse row): fast row access and matrix-vector products
- CSC (compressed sparse column): fast column access, common in solvers
- diagonal and banded formats for structured patterns
Why direct solvers struggle
Factorizing a sparse matrix with LU or Cholesky can create nonzeros where the original had zeros, a phenomenon called fill-in. Fill-in can turn a lean sparse matrix into a dense factor, exhausting memory. Reordering the rows and columns to minimize fill, using algorithms like approximate minimum degree or nested dissection, is essential for sparse direct methods.
Iterative solvers
For the largest systems, iterative methods such as conjugate gradient (for symmetric positive-definite matrices) and GMRES (for general matrices) avoid factorization entirely. They need only the ability to multiply the matrix by a vector, an operation whose cost scales with the number of nonzeros, and they converge faster with a good preconditioner.
import numpy as np
from scipy.sparse import csr_matrix
from scipy.sparse.linalg import spsolve
A = csr_matrix(np.array([[2.0, 0.0, 1.0], [0.0, 3.0, 0.0], [1.0, 0.0, 2.0]]))
b = np.array([3.0, 6.0, 3.0])
print(spsolve(A, b))Finite-element and finite-difference models of magnets, structures, and plasma fields produce enormous sparse matrices; exploiting their sparsity is the difference between a simulation that runs and one that cannot fit in memory.