PYTHON / NUMPY
Linear algebra with numpy.linalg
Use numpy.linalg to solve linear systems, compute determinants, norms, eigenvalues and least-squares fits, and pick between solve, lstsq and pinv.
What you will learn
- Use @ for matrix products; * stays elementwise even on 2-D arrays
- Call np.linalg.solve(A, b) instead of inv(A) @ b: less work, smaller error
- Reach for lstsq or pinv when A is non-square or rank-deficient
- np.linalg acts on the last two axes, so stacked arrays give batched results
Understanding Linear algebra with numpy.linalg
A NumPy array is not a matrix type; it is a block of numbers, and the operators keep their elementwise meaning. A * B multiplies entry by entry (with broadcasting), while A @ B is the matrix product. Everything in numpy.linalg agrees on one convention: the last two axes of an array are the rows and columns of a matrix, and any leading axes are just a stack of such matrices. That is why np.linalg.det on shape (500, 3, 3) returns 500 determinants instead of an error.
The functions themselves are thin wrappers over LAPACK, so the practical skill is choosing the right one for the shape and structure you actually have. For a square, full-rank A, np.linalg.solve factors A once with LU plus partial pivoting and back-substitutes, which costs roughly a third of the arithmetic of forming np.linalg.inv(A) and then multiplying, and it rounds fewer times. So the mental model is "solve a system", not "invert a matrix": inv is for the rare case when you need the inverse entries themselves.
When A is not square or its columns are dependent there is no unique solution, and solve will raise LinAlgError. np.linalg.lstsq minimises the 2-norm of A @ x - b, and np.linalg.pinv gives the same minimum-norm answer through the SVD. For diagnostics, ignore det: it is almost never exactly zero in floating point and it scales with the entries of A. np.linalg.matrix_rank and np.linalg.cond are the scale-aware tools, and np.linalg.norm(A @ x - b) is the cheapest sanity check that your answer is real.
import numpy as np
A = np.array([[3.0, 1.0],
[1.0, 2.0]])
b = np.array([9.0, 8.0])
x = np.linalg.solve(A, b) # solves A @ x == b
print("x =", x)
print("A @ x =", A @ x)
print("residual =", np.linalg.norm(A @ x - b))
print("det(A) =", round(float(np.linalg.det(A)), 10))
print("cond(A) =", round(float(np.linalg.cond(A)), 4))
Ainv = np.linalg.inv(A)
print("inv(A) @ b =", np.round(Ainv @ b, 12))
print("A * A (elementwise) ->", (A * A).tolist())
print("A @ A (matmul) ->", (A @ A).tolist())numpy.linalg reads the last two axes of an array as a matrix and delegates to LAPACK, so you choose a routine by the matrix's shape and structure and solve systems instead of inverting matrices.
Worked examples
Least squares for an overdetermined system
Fitting a straight line to four points, where the system has more equations than unknowns and solve cannot be used.
import numpy as np
t = np.array([0.0, 1.0, 2.0, 3.0])
y = np.array([1.0, 3.0, 5.0, 7.0]) # exactly y = 2t + 1
M = np.column_stack([t, np.ones_like(t)]) # design matrix, shape (4, 2)
coef, residuals, rank, sv = np.linalg.lstsq(M, y, rcond=None)
print("M.shape =", M.shape)
print("slope,intercept =", np.round(coef, 12))
print("rank =", rank)
print("fitted y =", np.round(M @ coef, 12))Example explained
Line 1column_stack builds a 4x2 matrix whose second column of ones supplies the intercept term.
Line 2np.linalg.solve(M, y) would raise LinAlgError here because M is not square; lstsq accepts any shape.
Line 3lstsq returns four things: the coefficients, the summed squared residual, the numerical rank, and the singular values.
Line 4rcond=None selects the modern cutoff for tiny singular values and avoids a deprecation warning on older NumPy.
Eigenvalues of a symmetric matrix with eigh
eigh exploits symmetry to return real, ascending eigenvalues and an orthonormal set of eigenvectors stored as columns.
import numpy as np
S = np.array([[2.0, 1.0],
[1.0, 2.0]])
vals, vecs = np.linalg.eigh(S)
print("eigenvalues =", np.round(vals, 12))
print("orthonormal =", np.allclose(vecs.T @ vecs, np.eye(2)))
for lam, v in zip(vals, vecs.T): # eigenvectors are COLUMNS of vecs
print(f"S @ v == {lam:.0f} * v :", np.allclose(S @ v, lam * v))
print("rebuilt S =", np.round(vecs @ np.diag(vals) @ vecs.T, 12).tolist())Example explained
Line 1eigh assumes the matrix is symmetric and reads only one triangle, so eigenvalues come back real and sorted ascending.
Line 2vecs[:, k] is the eigenvector for vals[k]; iterating over vecs.T is how you walk the columns.
Line 3vecs.T @ vecs equals the identity because the eigenvectors of a symmetric matrix are orthonormal.
Line 4vecs @ diag(vals) @ vecs.T reconstructs S, which is the spectral decomposition written as array operations.
Singular matrix: LinAlgError versus pinv
A rank-deficient system makes solve fail, while matrix_rank diagnoses it and pinv returns the minimum-norm solution.
import numpy as np
A = np.array([[1.0, 2.0],
[2.0, 4.0]]) # row 2 is exactly 2 * row 1
b = np.array([3.0, 6.0])
print("rank =", np.linalg.matrix_rank(A))
try:
np.linalg.solve(A, b)
except np.linalg.LinAlgError as err:
print("solve ->", type(err).__name__, ":", err)
x = np.linalg.pinv(A) @ b
print("pinv x =", np.round(x, 6))
print("A @ x =", np.round(A @ x, 6))Example explained
Line 1matrix_rank counts singular values above a tolerance, reporting 1 because the two rows are linearly dependent.
Line 2solve hits a zero pivot during LU factorization and raises LinAlgError('Singular matrix') instead of guessing.
Line 3pinv builds the Moore-Penrose pseudo-inverse from the SVD and returns the shortest of the infinitely many solutions.
Line 4A @ x reproduces b exactly here because b lies in the column space of A; np.linalg.lstsq(A, b, rcond=None)[0] gives the same x.
Important notes
These routines promote integer input to float64 and always return floating results; np.linalg.eig on a real non-symmetric matrix can legitimately return complex eigenvalues, so never assume a real dtype.
Since NumPy 2.0, np.linalg.solve treats a 2-D b as a matrix, not as a stack of vectors; if you mean a stack of column vectors give b an explicit trailing axis of size 1.
Common mistakes
Writing A * B when a matrix product was meant: NumPy multiplies entry by entry (or broadcasts) and returns wrong numbers with no error, so the bug shows up much later.
Computing np.linalg.inv(A) @ b instead of np.linalg.solve(A, b): more arithmetic and more rounding, and on an ill-conditioned A the explicit inverse can be nearly meaningless while solve still returns a small residual.
Testing singularity with np.linalg.det(A) == 0: floating-point determinants almost never land on exactly zero, and det scales with the entries (scaling a 3x3 matrix by 10 multiplies its det by 1000), so use matrix_rank or cond instead.
Try it yourself
Change, predict, then run
Build A = np.array([[4.0, 3.0], [6.0, 3.0]]) and b = np.array([10.0, 12.0]), solve for x with np.linalg.solve, and print both np.linalg.norm(A @ x - b) and np.linalg.cond(A). Then compute np.linalg.inv(A) @ b and check whether the two answers agree to every printed digit.
Open the Python workspaceCheck your understanding
For a square, non-singular A, both np.linalg.solve(A, b) and np.linalg.inv(A) @ b produce x. Why is solve the recommended call?
- It factors A once with LU and back-substitutes, roughly a third of the arithmetic of building a full inverse, and it avoids the extra rounding an explicit inverse introduces
- It is exact, because solve works with rational arithmetic internally while inv uses floating point
- inv only accepts symmetric positive-definite matrices, so it fails on general square matrices
- solve also accepts non-square A, which makes it the more general function
Show answer
solve does an LU factorization with partial pivoting plus two triangular substitutions (about n^3/3 operations), while inv effectively solves n systems and then needs a matrix-vector product, and every entry of A is round-tripped through the inverse, which enlarges the error. Option 4 is tempting but wrong: solve requires the last two dimensions to be square and raises LinAlgError otherwise; lstsq and pinv are the functions for non-square systems. Nothing in numpy.linalg is exact rational arithmetic, and inv works on any non-singular square matrix.