Numerical Linear Algebra

Linear algebra that survives contact with floating point: which algorithms to use, what they cost, and how answers degrade gracefully instead of silently.

Cost and stability split the methods: direct solvers bill cubic flops once, iterative solvers bill per step, and forming the Gram matrix squares the condition number

Why Does This Exist?

Textbook linear algebra assumes exact arithmetic. Computers round every operation to sixteen decimal digits, and at scale those roundings compound, interact with algorithm choices and occasionally detonate. Numerical linear algebra is the discipline of getting trustworthy answers anyway: algorithms picked for stability (not just flop counts), error analysed backwards (the computed answer is exact for a nearby problem), and costs known before code runs.

Two themes organise the field. Direct versus iterative: factor once at $O(n^3)$ and solve exactly (up to rounding), or refine step by step at $O(n^2)$ per iteration and stop when close enough. Stability: backward-stable algorithms (Householder QR, SVD, LU with pivoting) return the exact answer to a microscopically perturbed problem, so their error is bounded by the problem's conditioning, never the algorithm's temperament. The floating-point facts live in floating point formats; the failure catalogue in numerical stability.

The one idea you need first

A flop is one floating-point add or multiply. Gaussian elimination on $n \times n$ costs $2n^3/3$ flops: at $n = 1{,}000$ that is about $6.7 \times 10^8$ , under a second on a modern core. Cost models like this decide direct versus iterative before any code runs, and computational complexity for ML continues the story onto accelerators.

Think of It Like This

Carpentry with wood that swells

A carpenter's plans assume exact lengths, but wood swells with humidity: every cut lands within a tolerance, never exactly. Fine carpentry does not demand perfect wood; it uses joints that absorb the slack (dovetails, not butt joints stressed in tension) and measures the final fit rather than every cut.

Backward-stable algorithms are dovetail joints: their rounding errors are absorbed into a microscopic perturbation of the input, so the answer is exactly right for wood cut a hair off. Unstable algorithms are tension-loaded butt joints: the same slack cracks the assembly. And the condition number is the humidity reading: damp wood defeats even good joints, which is why some problems need reformulating, not just careful code.

How It Actually Works

The working map:

Direct methods: LU with partial pivoting (general square), Cholesky (symmetric positive definite, half the flops), QR (least squares, stable), SVD (everything, slowest, most revealing). One factorisation, then cheap solves.
Iterative methods: Jacobi, Gauss-Seidel, conjugate gradients, Krylov families. Each step costs a matvec; convergence depends on conditioning and preconditioning. The only option past $n \approx 10^5$ dense or for sparse giants.
Backward error: algorithm output $\hat{x}$ satisfies $(A + \Delta A)\hat{x} = b$ with tiny $\Delta A$ . Error then splits cleanly: algorithm contributed $\Delta A$ , conditioning amplified it. Blame assigned correctly.
BLAS levels: level 1 (vector ops, memory-bound), level 2 (matvec), level 3 (matmul, compute-bound). Fast code casts work into level 3, which is why blocked algorithms win.

Worked example

Price Gaussian elimination at $n = 1{,}000$ . Elimination step $k$ updates an $(n-k) \times (n-k)$ block with one multiply-subtract per entry: $\sum_{k=1}^{n} 2(n-k)^2 \approx 2n^3/3$ flops. For $n = 1{,}000$ : $2 \times 10^9 / 3 \approx 6.7 \times 10^8$ flops. A core sustaining $10^{10}$ flops/s finishes in under a tenth of a second; at $n = 10{,}000$ the bill is a thousand times larger (cubic scaling) and an iterative method starts looking attractive. Counting first, coding second.

Code

The point is the check: solve two ways and compare against the exact answer, watching the stable route win on residuals.

1import numpy as np2
3
4A = np.array([[3.0, 1.0], [1.0, 2.0]])5b = np.array([9.0, 8.0])6x_direct = np.linalg.solve(A, b)  # LU with pivoting: backward stable7x_cramer = np.linalg.det(np.column_stack([b, A[:, 1]])) / np.linalg.det(A)8print(x_direct.tolist())  # -> [2.0, 3.0]9print(abs(x_cramer - 2.0) < 1e-12)  # -> True

Both agree here because $2 \times 2$ is kind. At $n = 500$ with a middling condition number, Cramer's determinant ratios disintegrate while solve holds: same mathematics, different joints.

Watch Out For

Solving the normal equations instead of factoring X

Forming $X^T X$ squares the condition number before the solver even starts, destroying exactly the small singular directions regularisation or truncation was meant to handle gently. Factor $X$ itself (QR for speed, SVD for insight) and solve without ever assembling the Gram matrix.

Benchmarking algorithms at toy sizes

At $n = 50$ , overhead dominates and every method reports microseconds; the $O(n^3)$ versus $O(n^2)$ -per-step trade only materialises past $n \approx 10^3$ . Benchmark at production sizes with production sparsity, or the ranking inverts on deployment day.

The Quick Version

Numerical linear algebra pairs each algorithm with its flop cost and its stability certificate.
Gaussian elimination bills $2n^3/3$ flops: 670 million at $n = 1{,}000$ , a thousand times more at $10{,}000$ .
Backward-stable methods return exact answers to nearby problems; conditioning does the rest.
Direct factors once; iterative refines per matvec and owns the sparse-giant regime.
Never form $X^T X$ when QR or SVD on $X$ answers directly.