How to solve a sparse linear system with SciPy

Sparse linear systems appear when a model, graph, or discretized equation has many zero coefficients but still needs a numerical solution. SciPy can solve A x = b without first turning the coefficient matrix into a dense NumPy array, which keeps memory use tied to stored values instead of every possible matrix entry.

scipy.sparse.linalg.spsolve() is the direct sparse solver for a square coefficient matrix and a vector or matrix right-hand side. Passing a CSR or CSC sparse array keeps the input in a compressed format that is ready for sparse linear algebra.

A compact tridiagonal system keeps the sparse format, solution vector, and residual visible in a short terminal run. For symmetric positive-definite systems that are too large for a direct factorization, an iterative method such as scipy.sparse.linalg.cg() may fit better, but the same residual check decides whether the returned vector satisfies the equation.

1. Create a Python script named sparse_linear_system_solve.py.

   sparse_linear_system_solve.py

   import numpy as np
   from scipy.sparse import diags_array, issparse
   from scipy.sparse.linalg import spsolve
    
   n = 5
   A = diags_array(
       [-np.ones(n - 1), 4 * np.ones(n), -np.ones(n - 1)],
       offsets=[-1, 0, 1],
       shape=(n, n),
       format="csr",
   )
   b = np.array([2.0, 1.0, 0.0, 1.0, 2.0])
    
   x = spsolve(A, b)
   residual = A @ x - b
    
   print(f"matrix format: {A.format}")
   print(f"matrix shape: {A.shape}")
   print(f"stored values: {A.nnz}")
   print("solution:", np.round(x, 6))
   print(f"residual norm: {np.linalg.norm(residual):.3e}")
   print("passes tolerance:", np.linalg.norm(residual) < 1e-12)
   print("matrix is sparse:", issparse(A))

   diags_array() builds the three stored diagonals directly as a sparse array. The explicit format=“csr” keeps the coefficient matrix in compressed sparse row format before solving.

2. Run the script.

   $ python3 sparse_linear_system_solve.py
   matrix format: csr
   matrix shape: (5, 5)
   stored values: 13
   solution: [0.615385 0.461538 0.230769 0.461538 0.615385]
   residual norm: 0.000e+00
   passes tolerance: True
   matrix is sparse: True

3. Confirm that the residual norm is within the tolerance needed by your calculation.

   The residual is computed from the original sparse matrix multiplication, A @ x - b. A small residual shows that the returned vector satisfies the system more directly than inspecting the solution values alone.

4. Keep the coefficient matrix square and the right-hand side length matched to its rows.

   spsolve() is for square systems. A singular matrix or mismatched right-hand side can raise an error or warning instead of returning a usable solution.

5. Use a dense solver only when the coefficient matrix is already dense or the dense form is intentionally small.

   For a dense NumPy matrix, use scipy.linalg.solve() instead of converting a sparse workflow into dense arrays only to call spsolve().

6. Remove the demo script when it was only created for the check.

   $ rm sparse_linear_system_solve.py