How to solve a linear system with SciPy

Linear systems appear in numerical models, simulations, and matrix calculations where a coefficient matrix and a right-hand side determine unknown values. SciPy provides dense linear algebra routines for this job, so a Python script can solve A x = b and check the returned vector against the original equation.

scipy.linalg.solve() expects a square dense coefficient matrix and a matching vector or matrix right-hand side. It calls LAPACK routines and returns the solution array without forming an explicit matrix inverse, so the residual check comes from multiplying the original matrix by the returned vector.

A compact 3 x 3 system keeps the matrix, solution, reconstructed right-hand side, residual, and condition number visible in one run. A low residual shows that the returned vector satisfies the sample equation, while a high condition number would warn that small input changes can move the solution.

1. Create a Python script named linear_system_solve.py.

   linear_system_solve.py

   import numpy as np
   from scipy.linalg import solve
    
    
   A = np.array(
       [
           [3.0, 1.0, -1.0],
           [2.0, 4.0, 1.0],
           [-1.0, 2.0, 5.0],
       ]
   )
   b = np.array([4.0, 1.0, 1.0])
    
   x = solve(A, b)
   reconstructed = A @ x
   residual = reconstructed - b
   residual_norm = np.linalg.norm(residual, ord=np.inf)
   condition_number = np.linalg.cond(A)
    
   np.set_printoptions(precision=6, suppress=True)
   print("matrix A:")
   for row in A:
       print(row)
   print("right hand side:", b)
   print("solution x:", x)
   print("A @ x:", reconstructed)
   print("residual:", residual)
   print(f"residual_inf_norm: {residual_norm:.2e}")
   print(f"condition_number: {condition_number:.2f}")
   print("verified:", bool(residual_norm < 1e-10))

   solve() is for square systems. Use scipy.linalg.lstsq() for least-squares problems where the matrix is rectangular.

2. Run the script.

   $ python linear_system_solve.py
   matrix A:
   [ 3.  1. -1.]
   [2. 4. 1.]
   [-1.  2.  5.]
   right hand side: [4. 1. 1.]
   solution x: [ 2. -1.  1.]
   A @ x: [4. 1. 1.]
   residual: [ 0.  0. -0.]
   residual_inf_norm: 6.66e-16
   condition_number: 5.15
   verified: True

3. Confirm that the reconstructed right-hand side matches the original values.

   The A @ x line should match b. A displayed value such as -0. is floating-point roundoff, not a separate negative result.

4. Check the residual norm before using the solution in later calculations.

   The residual_inf_norm value is near machine precision for this matrix, and verified: True confirms it is below the 1e-10 tolerance used by the script.

5. Review the condition number when the input matrix comes from real data.

   A singular matrix raises a linear algebra error, and an ill-conditioned matrix can produce a solution that changes sharply when the input values change.

6. Use the sparse solver when most matrix entries are zero.

   scipy.sparse.linalg.spsolve() keeps sparse matrices compressed instead of converting them to dense arrays.
   Related: How to solve a sparse linear system with SciPy

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

   $ rm linear_system_solve.py