AI Engineering Foundations 17 | Linear Systems: Seven Methods, One Warning Line
This is lesson 17 of our walkthrough of the open-source course AI Engineering from Scratch (58k GitHub stars, MIT license). It's a Build-type lesson, 120 minutes, with the first three linear algebra lessons as prerequisites.
Every fit is doing something two thousand years old
Every time you write model.fit(X, y), the same thing happens underneath: solving a linear system Ax = b, where A is your data matrix, b is your target, and x is the weights you want. Linear regression, least squares, every neural network layer's y = Wx + b, Gaussian process kernel decompositions, Mahalanobis distances via covariance inverses — all variants of this one equation. The course opens by calling Ax = b one of the oldest problems in mathematics that still drives your neural networks: the Chinese classic *Nine Chapters on the Mathematical Art* systematically performed elimination about 1,800 years before Gauss was born.
The lesson implements every mainstream solution method from scratch, ending with a table of seven methods. They answer two questions: how to solve, and whether the solution can be trusted.
A different way to read equations: recipes
Ax = b has two readings. Row-wise, each equation is a hyperplane and the solution is the intersection. Column-wise, each column of A is an ingredient, and solving means finding the recipe. Small example: first column [2, 1], second column [1, -1], target b = [5, 1] — exactly 2 parts of column one plus 1 part of column two. The column view is more fundamental: if b lies in the span of the columns, the system has a solution; if not, find the closest point in the span — which is exactly what least squares does.
Elimination, and where it breaks
Gaussian elimination clears columns one by one into an upper triangular matrix, then back-substitutes. A 100×100 system costs about a billion floating-point operations — instant on modern hardware.
The real danger is the pivot. Elimination divides by diagonal entries; a zero pivot crashes, a tiny pivot amplifies rounding error. The course shows a classic comparison: eliminating without pivoting gives a multiplier of 1000, while swapping rows first gives 0.001. Interestingly, both versions still compute x = 1.000 in the textbook figure — in float64, the error hides beyond the eighth significant digit. The author re-ran the computation in float32: without pivoting, x₁ = 0.99993, an error of 7.3 × 10⁻⁵ — two and a half digits lost out of seven. Partial pivoting (always swap in the row with the largest absolute value in the current column) isn't ceremony; it's insurance that keeps multipliers small and error amplification minimal.
Factorize once, solve a thousand times
Re-running elimination for each new b costs O(n³) every time. LU decomposition stores the intermediate product: A = LU, with L holding the multipliers and U the upper triangle. Each new b then needs only two triangular solves, Ly = b and Ux = y, each O(n²). The course claims LU saves "1000/3" in total work for 1000 solves on the same A — a number the author could not derive. The actual savings depend on size: 16× at n = 100, 143× at n = 1000, 625× at n = 10000. The direction is right — the more solves, the more the fixed factorization cost amortizes — but "1000/3" has no derivable basis.
Two close relatives:
- QR decomposition splits A into an orthogonal matrix times an upper triangle. Orthogonal matrices preserve lengths and angles, making QR exceptionally stable — the standard for least squares.
- Cholesky decomposition serves only symmetric positive definite matrices (covariance matrices, kernel matrices, Hessians of convex functions), at half the work and half the storage. Gaussian processes and ridge regression use it daily.
When there is no exact solution
In machine learning, A is usually tall and thin — a hundred data points, three parameters — so no x satisfies all equations. The best you can do is minimize ‖Ax − b‖². The minimizer satisfies the normal equations: AᵀA x = Aᵀb, derived by setting the gradient to zero. This *is* linear regression.
The course fits a line through (1,3), (2,5), (3,6), (4,8), computing AᵀA = [[4, 10], [10, 30]], Aᵀb = [22, 63], and prints the solution x = [1.5, 1.7]. The author solved it by hand: the slope is 1.6, not 1.7 (Cramer's rule: 32/20). Verification: slope 1.6 gives residual sum 0.2 vs 0.5 for 1.7, and only 1.6 makes the residuals orthogonal to both columns — the textbook's final step is a misprint. The dangerous thing in ML is never lacking an answer; it's using a suboptimal answer as if it were optimal.
The warning line: condition number
A fast solution isn't enough — can you trust it? The condition number κ (largest singular value divided by smallest) tells you: κ ≈ 10^k means losing about k digits of precision. float64 has 16 digits, so at κ ≈ 10¹⁶ the answer is effectively all noise and the matrix is numerically singular.
Textbooks say QR is more stable than the normal equations for least squares, usually without explaining why. The mechanism is one line: the normal equations square the condition number. κ(AᵀA) = κ(A)². The author built a 50×3 matrix with two nearly identical columns: κ(X) = 37.0, κ(XᵀX) = 1371 — exactly 37². Squaring doubles the digits lost; QR operates on X directly and only endures 37.
A fallback exists: the pseudoinverse A⁺ = VΣ⁺Uᵀ via SVD. It returns the exact solution when one exists, the least-squares solution when none does, and the minimum-norm solution when there are infinitely many. NumPy's lstsq and pinv use it internally. Universal, but slow — a safety net, not a first choice.
Adding 0.01 crushes 300 million back to 0.1
Exercise 3 is the best problem of the lesson, and the author re-ran it end to end. Build a 50×2 matrix whose second column is a copy of the first plus noise at the 10⁻¹⁰ scale. Then κ(XᵀX) = 2.7 × 10¹⁶ — right on float64's death line, 16 digits theoretically wiped out.
Solving without regularization: weights of −3.39×10⁸ and +3.39×10⁸. With two nearly identical columns, the model learns "subtract one from the other," amplifying a 10⁻¹⁰ difference three hundred million times to fit noise. Mathematically valid, practically garbage.
Add 0.01 × I — ridge regression: the condition number falls from 2.7 × 10¹⁶ to 7778, and the weights become −0.082 and −0.082. Nearly identical columns get nearly identical weights — sensible. The cost: residuals rise from 42.388 to 42.428, a thousandth. Regularization never buys a better fit; it buys a solution you can believe.
Iterative methods: when the matrix is too big to factorize
For sparse systems with millions of unknowns (large PDEs, huge kernel matrices), LU and Cholesky are too expensive. Conjugate gradient (CG) skips factorization: start from zero and step along carefully chosen directions until residuals are small enough.
The course says CG finds the exact solution in at most n steps — with the parenthetical caveat *in exact arithmetic*. That parenthesis is alive. On a 100×100 symmetric positive definite matrix: κ ≈ 5 converged in 21 steps, well under 100; κ ≈ 380 needed 138 steps — over the guarantee. Floating-point error breaks the "at most n steps" promise; better conditioning means faster convergence, so regularization helps here too.
Cross-checking the numbers before closing
The lesson's capstone code aligns the course's hand-written path (elimination for least squares, hand-rolled Cholesky for ridge) with NumPy's SVD route and sklearn: least squares agrees to 8.9 × 10⁻¹⁶, ridge to 2.2 × 10⁻¹⁶ — machine precision. Two fully independent algorithmic paths landing on the same number is rare, good verification.
Exercise 1 solves [[1,2,3],[4,5,6],[7,8,10]] x = [6, 15, 27] with the course's elimination, the course's LU, and NumPy's solver: all three give [3, −3, 3], verifiable by substitution row by row. Lesson 18 covers convex optimization — from solving equations to solving them under constraints.