How to solve a sparse linear system with SciPy

Many scientific models produce coefficient matrices in which most entries are zero. Keeping those matrices sparse lets SciPy solve A x = b without allocating a dense array for every possible coefficient.

SciPy's spsolve() accepts a square sparse coefficient array and a vector or matrix right-hand side. The sample system uses a diagonally dominant five-row CSR array, so it has a unique direct solution while keeping only 13 values in storage.

The solution values alone do not show whether they satisfy the original equations. Multiplying the sparse array by the returned vector and measuring the norm of A @ x - b provides a fail-capable check against a numerical tolerance.

1. Create sparse_linear_system_solve.py with the imports, sparse coefficient array, and right-hand side.

   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])

   diags_array() stores the main diagonal and two neighboring diagonals directly in CSR format.

2. Append the solve and residual calculations after the right-hand side definition.

   x = spsolve(A, b)
   residual = A @ x - b
   residual_norm = np.linalg.norm(residual)
   tolerance = 1e-12

   spsolve() requires a square coefficient array with a compatible right-hand side. An exactly singular array produces MatrixRankWarning instead of a usable solution.

3. Append the sparse-format and tolerance report after the residual calculations.

   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: {residual_norm:.3e}")
   print("passes tolerance:", residual_norm < tolerance)
   print("matrix is sparse:", issparse(A))

4. Run the completed 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

   The comparison uses tolerance as its acceptance limit. A False result means the computed vector does not meet that residual limit.