Linear algebra with numpy.linalg
Matrix products, solving systems without inverting them, and the decompositions that answer real questions reliably.
matmul, dot and the asterisk
The single most expensive misunderstanding in NumPy: * multiplies element by element, while @ performs a matrix product. They are different operations with different costs and different shapes.
A = np.array([[3.0, 1.0], [1.0, 2.0]])
B = np.array([[1.0, 0.0], [0.0, 1.0]])
A * B # element-wise, shape (2, 2)
A @ B # matrix product - use this
np.matmul(A, B)
A.dot(B) # same result for 2-D inputs
v = np.array([1.0, 2.0])
A @ v # matrix times vector -> shape (2,)
np.inner(v, v) # 5.0
np.outer(v, v) # shape (2, 2) outer product| Operation | Shapes | Note |
|---|---|---|
A * B | Broadcast to a common shape | Element-wise; a common source of silent bugs |
A @ B | (m, k) @ (k, n) | Preferred syntax; batches over leading axes |
A.dot(B) | (m, k) . (k, n) | Same for 2-D, different rules for higher dimensions |
np.einsum | Any | Explicit indices when the intent is unusual |
Solve, do not invert
A = np.array([[3.0, 1.0], [1.0, 2.0]])
b = np.array([9.0, 8.0])
x = np.linalg.solve(A, b) # the right way
np.allclose(A @ x, b) # True
np.linalg.inv(A) @ b # same answer: slower and less accurate
# many right-hand sides at once: b has shape (2, k), solved in one call
B = np.column_stack([b, np.array([1.0, 0.0])])
np.linalg.solve(A, B)
# overdetermined (more rows than columns): least squares
M = np.array([[1.0, 0.0], [0.0, 1.0], [1.0, 1.0]])
coef, residuals, rank, sv = np.linalg.lstsq(M, np.array([1.0, 2.0, 3.0]), rcond=None)solve uses a factorisation (typically LU) instead of forming an inverse, so it does roughly half the work and keeps more precision. Reach for inv only when the inverse itself is the answer you need.
Decompositions and norms
S = np.array([[2.0, 1.0], [1.0, 2.0]])
vals, vecs = np.linalg.eig(S) # general eigenvalues
vals, vecs = np.linalg.eigh(S) # symmetric/Hermitian: faster and real
U, s, Vt = np.linalg.svd(M, full_matrices=False) # singular values
L = np.linalg.cholesky(S) # S = L @ L.T for positive-definite S
np.linalg.norm(v) # Euclidean length
np.linalg.norm(M, axis=0) # per-column norms
np.linalg.matrix_rank(M)
np.linalg.cond(A) # large condition number means unstable
np.linalg.det(A)⚠️
A large condition number (
cond) means the answer is hypersensitive to tiny input changes. When you see one, prefer lstsq or add regularisation rather than trusting solve.FAQ
When should I use lstsq instead of solve?
Whenever the system has no exact solution — more equations than unknowns, or a rank-deficient matrix.
lstsq returns the closest fit and its residuals instead of raising.Why is my symmetric matrix giving complex eigenvalues?
Tiny asymmetries from floating-point arithmetic. Use
np.linalg.eigh, which assumes symmetry, returns real values and is considerably faster.Related
Performance: einsum, ufunc tricks and avoiding copies Views, copies and memory layout
Last refreshed 2026-09-18.