How do I solve a system of linear equations in NumPy?

The short version

Use np.linalg.solve(A, b). Do not compute the inverse with np.linalg.inv and multiply. solve is faster and numerically more accurate, and reaching for the inverse is a habit worth breaking, even though it gives the right answer on small well-behaved problems.

Solving directly

import numpy as np

A = np.array([[2.0, 1.0],
              [1.0, 3.0]])
b = np.array([5.0, 10.0])

np.linalg.solve(A, b)      # array([1., 3.])

This solves the system where 2x + y = 5 and x + 3y = 10, giving x = 1 and y = 3.

The inverse approach, and why to avoid it

np.linalg.inv(A) @ b       # also array([1., 3.])

It gives the same answer here, so what is wrong with it?

Computing an inverse does more work than solving, and it introduces additional rounding error at every step. On a well-conditioned two-by-two matrix the difference is invisible. On a large or nearly singular system, the inverted-and-multiplied result can be noticeably less accurate than the directly solved one, and in bad cases meaningfully wrong.

The rule from numerical analysis is simple and worth internalising: if you find yourself writing inv(A) @ b, write solve(A, b) instead. You almost never actually need the inverse itself.

Matrix multiplication

Use the @ operator, which is unambiguous.

A @ A        # matrix multiplication
A * A        # element-wise! a completely different thing

This is a classic error. * on NumPy arrays multiplies element by element, not as matrices. If your matrix maths gives results of the right shape but the wrong values, this is the first thing to check.

Singular matrices

If the system has no unique solution, solve raises rather than returning nonsense.

np.linalg.solve(np.array([[1.0, 2.0], [2.0, 4.0]]), b)
# LinAlgError: Singular matrix

The second row is twice the first, so the equations are not independent. The error is telling you something real about your problem, not about your code.

Nearly singular is worse than singular

A matrix that is almost singular does not raise, but the solution is extremely sensitive to tiny changes in the inputs, so the answer may be dominated by rounding error. Check the condition number.

np.linalg.cond(A)    # large means ill-conditioned, treat the result with suspicion

A condition number in the millions means you should not trust the digits you are seeing. This matters in regression, where highly correlated predictors produce exactly this situation, which is why collinearity makes coefficient estimates unstable.

Least squares for over-determined systems

If there are more equations than unknowns, there is generally no exact solution, and solve will not work. lstsq finds the best fit in the least-squares sense, which is precisely what fitting a linear regression does.

np.linalg.lstsq(A, b, rcond=None)[0]