How to solve linear equations with NumPy

Linear systems appear whenever several unknown values must satisfy equations at the same time. NumPy handles that calculation as a matrix solve, so a coefficient array and right-hand-side vector can produce the unknown vector without manually expanding elimination steps.

The np.linalg.solve() function solves A @ x = b for square, full-rank coefficient matrices. The coefficient array carries one equation per row, while the right-hand-side vector holds the target values for those rows.

Reconstructing A @ x after the solve keeps the result tied to the original equations. np.allclose() allows normal floating-point roundoff, and a residual value near zero shows that the computed vector satisfies the supplied system.

1. Create a Python script that solves a square coefficient matrix against a right-hand-side vector.

   linear-equation-solve.py

   import numpy as np
    
   coefficients = np.array(
       (
           (3.0, 1.0, 2.0),
           (2.0, 6.0, 1.0),
           (1.0, -1.0, 4.0),
       )
   )
   rhs = np.array((11.0, 1.0, 15.0))
    
   solution = np.linalg.solve(coefficients, rhs)
   reconstructed = coefficients @ solution
   residual = reconstructed - rhs
    
   print("solution:", np.round(solution, 6))
   print("reconstructed rhs:", np.round(reconstructed, 6))
   print("residual norm:", round(float(np.linalg.norm(residual)), 12))
   print("matches rhs:", np.allclose(reconstructed, rhs))

   Each row in coefficients represents one equation. The vector rhs must have one value for each row.

2. Run the script and compare the reconstructed right-hand side.

   $ python linear-equation-solve.py
   solution: [ 2. -1.  3.]
   reconstructed rhs: [11.  1. 15.]
   residual norm: 0.0
   matches rhs: True

   matches rhs: True confirms that coefficients @ solution matches rhs within the np.allclose() tolerance.

3. Confirm that a singular matrix fails before using np.linalg.solve() for dependent equations.

   $ python - <<'PY'
   import numpy as np
   
   try:
       np.linalg.solve(np.array(((1.0, 2.0), (2.0, 4.0))), np.array((3.0, 6.0)))
   except np.linalg.LinAlgError as error:
       print(type(error).__name__ + ":", error)
   PY
   LinAlgError: Singular matrix

   np.linalg.solve() requires a square, full-rank coefficient matrix. Use np.linalg.lstsq() for rectangular systems or least-squares fits.