38 · Linear algebra with np.linalg
Solve linear systems, decompose matrices, and compute norms — the foundational linalg operations behind linear regression, PCA, and every gradient-based optimiser.
`np.linalg` is NumPy's gateway to LAPACK — a battle-tested Fortran library that has been the engine of numerical linear algebra for 50 years. Three operations cover almost every ML use case: **`solve(A, b)`** for linear systems (regression normal equations, any `Ax = b`), **`eig(A)`** for eigendecomposition (covariance analysis, graph algorithms), and **`svd(A)`** for singular value decomposition (PCA, low-rank approximation, pseudo-inverse).
Without this:
Without `np.linalg`, implementing linear regression from scratch requires either a for-loop gradient descent or deriving the normal-equation closed form by hand and implementing `inv(X.T @ X) @ X.T @ y` — which is numerically fragile. `np.linalg.solve` and `np.linalg.lstsq` handle the numerical edge cases (near-singular matrices, overdetermined systems) that the manual inverse approach silently gets wrong.
Python's np.linalg module is a thin wrapper around LAPACK routines. The key functions and when to use each:
np.linalg.solve(A, b)— solvesAx = bfor squareA. Prefer this overnp.linalg.inv(A) @ b:solveuses LU factorisation internally, which is faster and avoids the numerical amplification of errors that comes with explicitly computingA⁻¹.np.linalg.lstsq(A, b)— least-squares solution whenAis not square (overdetermined systems). This is whatsklearn.LinearRegressionuses internally.np.linalg.eig(A)— returns(eigenvalues, eigenvectors). For symmetric matrices, prefernp.linalg.eighwhich exploits symmetry for speed and guarantees real output.np.linalg.svd(A)— returns(U, s, Vt)whereA = U @ np.diag(s) @ Vt. SVD exists for any matrix (unlike eigendecomposition which requires square matrices).np.linalg.norm(v)— defaults to the Euclidean (L2) norm. Passord=1for L1,ord=np.inffor max-norm.np.linalg.det(A)— determinant; non-zero means full rank.np.linalg.matrix_rank(A)— numerical rank (uses SVD internally).
Python (in browser)
`np.linalg.solve` is preferable to `inv(A) @ b` because computing the inverse amplifies floating-point errors — for a well-conditioned 2×2 system the difference is negligible, but for large or ill-conditioned systems `solve` can be orders of magnitude more accurate. The residual `norm(Ax - b)` is the standard way to verify a solution.
Python runs entirely in your browser via Pyodide (~6 MB on first Run, cached after).
Python (in browser)
The defining property of an eigenvector: `A @ v == lambda * v` — the matrix only scales the vector, never rotates it. For a covariance matrix, the eigenvectors are the principal components (directions of maximum variance) and the eigenvalues tell you how much variance each component explains. This is the mathematical core of PCA.
Python runs entirely in your browser via Pyodide (~6 MB on first Run, cached after).
Python (in browser)
`full_matrices=False` returns the 'thin' or 'economy' SVD — `U` has shape `(m, k)` instead of `(m, m)` where `k = min(m, n)`. The singular values in `s` are in descending order; truncating to the top `k` gives the best rank-`k` approximation to `A` in the Frobenius norm. PCA of a centred data matrix is exactly this: the right singular vectors `Vt` are the principal components.
Python runs entirely in your browser via Pyodide (~6 MB on first Run, cached after).
`np.linalg.lstsq` uses SVD internally to handle singular or near-singular design matrices gracefully. When you use `sklearn.LinearRegression`, it calls LAPACK's `dgelsd` — the same algorithm. Understanding this connection means you can diagnose rank-deficiency warnings ('Matrix is close to singular'), choose regularisation (`Ridge` adds `lambda * I` to make `X.T @ X` invertible), and interpret `rcond` (the numerical threshold below which singular values are treated as zero).
What does `np.linalg.norm(v)` return by default for a 1D vector `v`?
MML lesson 19 covers eigenvalues and eigenvectors — the mathematical objects that `np.linalg.eig` and `np.linalg.eigh` compute. The defining equation `Av = λv` (a matrix scales certain vectors without rotating them) is the foundation of PCA (eigenvectors of the covariance matrix), spectral graph theory, and Markov chain analysis. Lesson 22 covers SVD — the generalisation of eigendecomposition to non-square matrices that powers PCA, pseudo-inverse, and low-rank approximation.
- `np.linalg.solve(A, b)` solves `Ax = b` and is always preferred over `inv(A) @ b` — it is faster and numerically more stable.
- `np.linalg.eigh(A)` returns `(eigenvalues, eigenvectors)` for symmetric matrices. Each eigenvector satisfies `A @ v == lambda * v`.
- `np.linalg.svd(A, full_matrices=False)` returns `(U, s, Vt)`. Singular values in `s` are in descending order; truncating to the top `k` gives the best rank-`k` approximation.
- `np.linalg.lstsq(A, b)` solves the overdetermined least-squares problem — the math behind `sklearn.LinearRegression`.
- `np.linalg.norm(v)` defaults to the L2 (Euclidean) norm. Pass `ord=1` for L1 and `ord=np.inf` for max-norm.
PCA is `np.linalg.svd` (or `np.linalg.eigh`) on a centred data matrix — the top-`k` right singular vectors are the principal components. Ordinary least squares (linear regression) is `np.linalg.lstsq`. Gradient clipping is `np.linalg.norm(grad)`. Cosine similarity is `(a @ b) / (norm(a) * norm(b))`. Every recommendation-system matrix factorisation is SVD. The Gram–Schmidt process (QR factorisation, accessible via `np.linalg.qr`) underpins numerically stable implementations of both gradient descent and least squares.
If you remove it: Without `np.linalg`, implementing PCA from scratch would require coding Gram–Schmidt orthogonalisation, the QR algorithm for eigenvalues, and Golub–Reinsch for SVD — all of which took decades of research to stabilise numerically. `np.linalg` gives you production-quality implementations of all of these in one import, letting you focus on the ML problem rather than numerical analysis.