TopolearnRegression, graphics, physics and ML all reduce to a few linear-algebra calls; knowing the right one avoids slow and unstable code.
* on arrays is element-wise. For matrix multiplication use the `@` operator (np.matmul): A @ B for matrices, A @ v for matrix-vector. A 1-D vector is treated as a row or column as needed, and the result drops that dimension again.
np.linalg holds the standard toolkit, backed by optimised LAPACK/BLAS libraries:
np.linalg.solve(A, b)solvesA x = b. Prefer it tonp.linalg.inv(A) @ b; solving is faster and more numerically accurate than forming an inverse.np.linalg.lstsq(A, b)for least squares whenAis not square (e.g. fitting a line).np.linalg.norm(x)for vector length or matrix norms, withaxis=for many vectors.np.linalg.det,inv,matrix_rank,eig(general),eigh(symmetric, faster and real),svd,qr,cholesky.
Most of these functions are batched: given a stack of matrices with shape (n, k, k), @ and solve operate on each matrix at once, broadcasting over the leading axes. That avoids Python loops over many small matrices.
NumPy 2 brought helpers from the array API standard: a.mT transposes only the last two axes (unlike .T, which reverses all of them), and np.vecdot(x, y) computes dot products along the last axis of stacks of vectors.
np.dot predates @ and has odd rules for N-D inputs; use @ for matrices and np.dot or np.vecdot for vectors. Avoid the np.matrix class; it's deprecated in spirit and plain arrays do everything it did.
Singular systems raise np.linalg.LinAlgError. Nearly singular ones don't raise but give huge, unreliable results; np.linalg.cond(A) warns you.
Example
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)
print(x) # [2. 3.]
print(np.allclose(A @ x, b)) # True
# least-squares line fit y = m*t + c
t = np.array([0.0, 1.0, 2.0, 3.0])
y = np.array([1.1, 2.9, 5.2, 6.8])
X = np.column_stack([t, np.ones_like(t)])
(m, c), *_ = np.linalg.lstsq(X, y)
print(round(m, 2), round(c, 2)) # 1.94 1.09
vals, vecs = np.linalg.eigh(A) # symmetric matrix
print(vals.round(3)) # [1.382 3.618]
stack = np.stack([A, 2 * A, 3 * A]) # shape (3, 2, 2)
print((stack @ b).shape) # (3, 2) batched
print(np.linalg.norm([3.0, 4.0])) # 5.0Remember
@is matrix multiplication;*is element-wise- Use
solve(orlstsq) rather than computing an inverse np.linalgfunctions broadcast over stacks of matriceseighfor symmetric matrices; checkcondfor near-singular systems
Watch out for
- Writing
A * Bfor a matrix product and getting an element-wise result with no error when shapes happen to match.