Linear Algebra with NumPy

NumPy's linalg module provides the standard matrix operations - products, solves, inverses, decompositions, eigenproblems - backed by optimized BLAS/LAPACK routines. The practical emphasis is on the right tool per job: dot/@ for products, linalg.solve for "Ax = b" (never the inverse), and the decompositions for everything that needs the structure of a matrix. This chapter is the daily linear algebra; the decompositions get their own chapter next.

Products

import numpy as np

v = np.array([1, 2, 3])
w = np.array([4, 5, 6])
print(v @ w)              # 32 - vector dot (scalar)
print(np.dot(v, w))       # 32

A = np.array([[1, 2], [3, 4]])
B = np.array([[5, 6], [7, 8]])
print(A @ B)              # matrix product [[19 22] [43 50]]
print(A @ v)              # matrix-vector: [ 8 14]
print(A.T @ A)            # (2,2) - the Gram matrix, always symmetric PSD

# batched: stacks of matrices multiply elementwise by stack index
X = np.random.default_rng(0).normal(size=(5, 2, 3))
W = np.random.default_rng(1).normal(size=(5, 3, 4))
print((X @ W).shape)      # (5, 2, 4)

@ handles 1-D dot, 2-D matrix product, and batched N-D products (last two axes are the matrix dimensions, leading axes are batched). A.T @ A is the workhorse expression of least squares and similarity - it is always symmetric positive semidefinite, which is why it appears everywhere in the math.

Solving Linear Systems

A = np.array([[3.0, 1.0], [1.0, 2.0]])
b = np.array([9.0, 8.0])

x = np.linalg.solve(A, b)
print(x)                          # [2. 3.]
print(A @ x - b)                  # ~[0 0] - the residual, should be ~0

# multiple right-hand sides at once: B is (2, 2), solve for both columns
B = np.array([[9.0, 1.0], [8.0, 2.0]])
X = np.linalg.solve(A, B)
print(A @ X - B)                  # ~ zero matrix

linalg.solve(A, b) is the correct way to compute A^{-1} b: it factors A once (LU) and backsolves, which is both faster and more accurate than forming the inverse and multiplying. Forming linalg.inv(A) and then doing inv(A) @ b is the classic mistake - double the work, double the roundoff. If you need the inverse as an object (it comes up), compute it, but do not use it to solve.

Norms, Rank, Condition

v = np.array([3.0, 4.0])
print(np.linalg.norm(v))          # 5.0 - the 2-norm (default)
print(np.linalg.norm(v, ord=1))   # 7.0 - L1
print(np.linalg.norm(v, ord=np.inf))  # 4.0 - max absolute value

A = np.array([[1, 2], [3, 4]])
print(np.linalg.matrix_rank(A))         # 2 - numerically, with a tolerance
print(np.linalg.cond(A))                # condition number: how ill-posed the solve is
print(np.linalg.det(A))                 # -2.0

norm defaults to the 2-norm for vectors and the 2-norm (largest singular value) for matrices. matrix_rank is a numeric rank - it counts singular values above a tolerance, which is the honest answer for floating-point data (an "exact" rank is a rational-arithmetic concept). cond is the condition number: a solve with cond ~ 1e16 in float64 has essentially no reliable digits, and that is the number to check when a system gives nonsense answers.

Eigenvalues and Eigenvectors

A = np.array([[2.0, 1.0], [1.0, 2.0]])
w, V = np.linalg.eig(A)
print(w)              # [3. 1.] - the eigenvalues
print(V)              # columns are the eigenvectors
print(A @ V[:, 0] - w[0] * V[:, 0])   # ~[0 0] - the eigen equation

# symmetric matrices: real eigenvalues, orthonormal eigenvectors
w_s, V_s = np.linalg.eigh(A)
print(V_s.T @ V_s)    # ~ identity

linalg.eig is the general (possibly complex) eigenproblem; linalg.eigh is the specialized solver for real symmetric / Hermitian matrices - faster, and it guarantees real eigenvalues with orthonormal eigenvectors. Use eigh whenever the matrix is symmetric (correlation matrices, Gram matrices, physics Hamiltonians) - it is strictly the better tool there.

QR and the Everyday Decompositions

A = np.random.default_rng(0).normal(size=(4, 3))
Q, R = np.linalg.qr(A)
print(Q.T @ Q)         # ~ identity - Q has orthonormal columns
print(Q @ R - A)       # ~ zero - the factorization holds

# least squares via the product (small case): the QR route is stable
b = np.random.default_rng(1).normal(size=4)
x, res, rank, sv = np.linalg.lstsq(A, b, rcond=None)
print(A @ x - b)       # the residual of the least-squares fit

qr factors A = Q R with orthonormal Q - the basis for stable least squares and orthogonalization. lstsq is the least-squares solver for overdetermined or rank-deficient systems, returning the solution, residual information, rank, and singular values in one call. For the full decomposition toolkit - SVD, eig, cholesky, and when each is the right one - the next chapter goes deeper.

Summary

Do linear algebra with the dedicated tools: @ for products (including batched), linalg.solve for systems - never inverse-then-multiply, norm/rank/cond for the numeric health checks, eig for general eigenproblems and eigh for the symmetric case, qr/lstsq for least squares. The condition number is the diagnostic for "why is this solve garbage", and the symmetric-specialized solvers are faster and more accurate whenever the symmetry holds. The decompositions themselves - what they compute and which one answers which question - are the next chapter.