On a real computer, linear algebra is done in finite-precision arithmetic, where the order of operations and the conditioning of a problem decide how many correct digits survive. This chapter consolidates the major matrix factorizations, introduces matrix norms and the all-important condition number, and surveys the algorithms that power , , BLAS, and LAPACK — the libraries every ML framework calls.
A field guide to the decompositions
Each factorization expresses a matrix as a product of structured pieces, trading a hard problem for easy ones (triangular or diagonal). Choosing the right one is the essence of practical numerical work.
| Factorization | Form | Use it for |
|---|---|---|
| LU (with pivoting) | \(\mP\mA=\mL\mU\) | solving square systems \(\mA\vx=\vb\); determinants |
| Cholesky | \(\mA=\mL\mL\tp\) | symmetric positive-definite systems (2\(\times\) faster); GP/covariance; sampling Gaussians |
| QR | \(\mA=\mQ\mR\) | least squares; orthonormal bases; eigenvalue iteration |
| Eigen | \(\mA=\mP\mD\mP\inv\) | square matrices; dynamics, powers, PCA via covariance |
| Spectral | \(\mA=\mQ\mLam\mQ\tp\) | symmetric matrices; PCA, kernels, Hessians |
| SVD | \(\mA=\mU\mSig\mV\tp\) | everything: any matrix, rank, low-rank, pseudoinverse |
A practical decision rule. Square and generic \(\to\) LU. Symmetric positive definite \(\to\) Cholesky (the fastest and most stable; it also fails precisely when the matrix is not positive definite, a useful test). Rectangular/least squares \(\to\) QR. Symmetric and you need eigen-structure \(\to\) spectral (eigh). Anything else, or when in doubt, or when rank/stability matters \(\to\) SVD. The SVD is the most informative and most stable, at the highest cost.
Cholesky in one line
For symmetric positive-definite \(\mA\), \(\mA=\mL\mL\tp\) with \(\mL\) lower-triangular. It uses half the work of LU, is very stable, and gives the log-determinant for free as \(2\sum_i\log L_{ii}\) — which is why it is the engine of Gaussian-process likelihoods and of drawing samples from a multivariate Gaussian (\(\vx=\vmu+\mL\vz\) with \(\vz\) standard normal).
Vector and matrix norms
Norms measure size and, crucially, how errors propagate. Beyond the vector \(p\)-norms (Chapter [ch:vectors]), the key matrix norms are:
Spectral norm \(\norm{\mA}_2=\sigma_1\), the largest singular value — the maximum factor by which \(\mA\) stretches any vector.
Frobenius norm \(\norm{\mA}_F=\sqrt{\sum_{ij}a_{ij}^2}=\sqrt{\sum_i\sigma_i^2}\) — the Euclidean norm of the flattened matrix.
Nuclear norm \(\norm{\mA}_*=\sum_i\sigma_i\) — the convex surrogate for rank, used to encourage low-rank solutions (matrix completion).
The spectral norm is an “operator” norm: it answers “what is the worst-case amplification of \(\mA\)?” The Frobenius norm treats the matrix as a long vector and is what appears in most ML loss functions (e.g. reconstruction error). The nuclear norm, being the sum of singular values, plays the role for rank that the L1 norm plays for sparsity — minimizing it tends to produce low-rank matrices.
The condition number
The condition number of an invertible matrix (in the 2-norm) is \[\cond(\mA)=\norm{\mA}_2\,\norm{\mA\inv}_2=\frac{\sigma_{\max}}{\sigma_{\min}}.\] It measures how much relative error in \(\vb\) (or in \(\mA\)) can be amplified in the solution of \(\mA\vx=\vb\). A small condition number (\(\approx1\)) is well-conditioned; a huge one is ill-conditioned.
The condition number is the most important number in numerical linear algebra. Rule of thumb: if \(\cond(\mA)\approx10^k\), you can lose about \(k\) digits of accuracy when solving with \(\mA\). In double precision (about 16 digits), a condition number of \(10^{16}\) means no correct digits. Ill-conditioning — not just singularity — is what makes problems hard, and it is why we avoid forming \(\mA\tp\mA\) (which squares the condition number) and prefer QR/SVD.
A matrix can be “technically invertible” yet useless: a tiny nonzero \(\sigma_{\min}\) gives an enormous condition number, so the computed inverse and solution are dominated by round-off. Never test singularity with the determinant; check \(\sigma_{\min}\) or the condition number. In ML, an ill-conditioned data matrix (from collinear features) both destabilizes the fit and signals that regularization is needed.
Algorithms behind the libraries
Direct solvers (LU, Cholesky, QR) compute an exact answer in \(O(n^3)\) up to round-off; ideal for dense, moderate-size systems.
Iterative solvers (conjugate gradient for symmetric PD systems; GMRES for general) only multiply by \(\mA\) and converge gradually — essential for the huge, sparse systems of scientific computing and some ML.
The power method and its refinements (Lanczos, Arnoldi, randomized SVD) find a few dominant eigen/singular vectors of enormous matrices without factoring the whole thing — exactly what PCA on big data and PageRank need.
Numerical linear algebra is the invisible foundation under every ML framework: and dispatch to BLAS/LAPACK, decades-old, exquisitely tuned Fortran/C kernels. Matrix multiply (GEMM) is so heavily optimized that the practical advice is always “express your computation as big matrix multiplications,” which is precisely how GPUs accelerate deep learning. You rarely write these algorithms, but knowing which one runs — and its conditioning behavior — is what separates robust code from mysterious numerical bugs.
Conditioning explains real training pathologies. The condition number of the loss Hessian sets gradient descent’s zig-zag and the maximum stable learning rate (Chapter [ch:eigen]); normalization layers and good initialization improve conditioning. Collinear features make \(\mX\tp\mX\) ill-conditioned, motivating Ridge regularization, which adds \(\lambda\mI\) and so raises \(\sigma_{\min}\) and lowers the condition number. Cholesky drives Gaussian processes; randomized SVD makes PCA feasible on web-scale matrices; iterative solvers appear inside second-order optimizers. Stability is not a footnote — it decides whether training converges.
import numpy as np
A = np.array([[1.0, 1.0], [1.0, 1.0001]]) # nearly singular
print(np.linalg.cond(A)) # ~ 4e4 (ill-conditioned)
s = np.linalg.svd(A, compute_uv=False)
print(s[0] / s[-1]) # same as cond(A)
# Cholesky for an SPD system, and a stable log-determinant:
M = A.T @ A + 1e-3*np.eye(2) # SPD
L = np.linalg.cholesky(M)
logdet = 2*np.sum(np.log(np.diag(L)))Chapter summary
Pick the factorization by structure: LU (square), Cholesky (SPD), QR (least squares), spectral (symmetric), SVD (anything/rank/stability).
Matrix norms: spectral (\(\sigma_1\), worst-case stretch), Frobenius (flattened L2), nuclear (\(\sum\sigma_i\), low-rank surrogate).
The condition number \(\sigma_{\max}/\sigma_{\min}\) predicts digit loss; ill-conditioning, not just singularity, is what makes problems hard.
Libraries call BLAS/LAPACK; express computation as matrix multiplies; prefer stable factorizations over explicit inverses.