Sparse Matrices & Large Systems (scipy.sparse)

In 1D, 2D, and 3D finite-difference or finite-element simulations with $N$ unknowns, each node interacts only with a few neighbors. Storing a dense $10^5 \times 10^5$ matrix requires $80\text{ GB}$ of RAM, whereas a sparse matrix uses only a few megabytes.

1. Sparse Matrix Formats in SciPy

Constructing and Inspecting Sparse Matrices

2. Solving Large Sparse Systems (scipy.sparse.linalg)

SciPy provides both direct and iterative solvers for sparse systems A x = b:

from scipy.sparse.linalg import spsolve
x = spsolve(A.tocsc(), b)

from scipy.sparse.linalg import cg
x, info = cg(A, b, rtol=1e-8)
# info == 0 indicates successful convergence

from scipy.sparse.linalg import splu
lu = splu(A.tocsc())
x_step = lu.solve(b_step)

Solving a 1D Poisson Boundary Value System

Let’s discretize $-u''(x) = 1$ with boundary conditions $u(0) = u(1) = 0$ using central differences:

3. Sparse Eigenvalue Problems (eigsh)

When finding the lowest vibrational modes or buckling loads of huge systems, you rarely need all $10^5$ eigenvalues. scipy.sparse.linalg.eigsh (Lanczos algorithm) extracts the few eigenvalues you specify:

4. Try It Yourself

Problem: Build and Solve a Sparse Diagonal System

Create a tridiagonal matrix of size $N=5$ with 4 on the main diagonal and -1 on both off-diagonals. Solve $A x = \mathbf{1}$ using sparse.linalg.spsolve.

import numpy as np
from scipy import sparse
from scipy.sparse.linalg import spsolve
 
N = 5
diags = [-np.ones(N-1), 4.0 * np.ones(N), -np.ones(N-1)]
A = sparse.diags(diags, [-1, 0, 1], format="csc")
b = np.ones(N)
 
x = spsolve(A, b)
print(f"x: {np.round(x, 4)}")