Solving Linear Systems
Solve Ax=b using np.linalg.solve and compare to the naive inverse-multiplication approach for numerical stability.
What Is a Linear System?
A system of linear equations can be written in matrix form as Ax = b, where A is the coefficient matrix, x is the vector of unknowns, and b is the right-hand side vector. For example, two equations with two unknowns can be expressed as a 2×2 matrix equation. Solving the system means finding x such that A multiplied by x equals b. NumPy's np.linalg.solve(A, b) handles this efficiently and accurately.
import numpy as np
# System: 2x + y = 5
# x + 3y = 10
A = np.array([[2.0, 1.0],
[1.0, 3.0]])
b = np.array([5.0, 10.0])
x = np.linalg.solve(A, b)
print('Solution x:', x)
print('Verify A @ x == b:', np.allclose(A @ x, b))Why Not Use the Inverse?
Mathematically, the solution to Ax = b is x = A⁻¹b. But computing np.linalg.inv(A) @ b is slower and less accurate than calling np.linalg.solve(A, b) directly. solve() uses LU decomposition internally, which requires fewer operations and accumulates less floating-point error. The rule of thumb: if you are solving for x in Ax=b, always use solve, never inv.
import numpy as np
A = np.random.rand(100, 100)
b = np.random.rand(100)
# Avoid this pattern
x_inv = np.linalg.inv(A) @ b
# Prefer this
x_solve = np.linalg.solve(A, b)
print('Max difference:', np.max(np.abs(x_inv - x_solve)))
# solve is faster AND more accurate