Ax = b asks which input x the transformation A turns into b. Two lines cross once; columns mix to hit a target; and when they cannot, least squares finds the closest reachable point.
The row picture asks where the lines (or planes) of the equations meet. The column picture asks how much of each column of A to mix to build b. The first is where constraints agree; the second is what the machine can reach.
x₁·col₁ + x₂·col₂ = bRow operations turn A into an upper-triangular U; the multipliers used along the way are exactly the lower-triangular L, so A = LU. Pay O(n³) once, then every new right-hand side costs only two triangular solves, O(n²).
A = LU · A = LLᵀ (symmetric positive definite)With more data points than parameters, b is usually outside the column space, so no x is exact. Minimize the squared error instead: the normal equations AᵀA x = Aᵀb give the line that balances every residual — this is linear regression.
minimize ‖Ax − b‖² → AᵀA x = AᵀbAx = b asks which input x the machine A turns into b: square and well-behaved gives one answer found by elimination, tall and noisy gives no exact answer so least squares takes the closest, and the condition number tells you how many digits of whichever answer to believe.
By the end you will be able to solve small systems by elimination with partial pivoting, explain why LU and Cholesky amortize the work, read an overdetermined problem as a projection, derive the normal equations, use κ to predict digit loss, and run a Jacobi sweep on a system too big to factor.
A linear system is a pile of straight-line constraints, and it is also a recipe for building a target out of columns. Both readings describe the same x, and each one answers a different question.
A linear system collects m equations in n unknowns. In matrix form it collapses into one line: Ax = b. The matrix A holds the known coefficients, b holds the known outputs, and x is the vector of unknowns we want. Every equation is linear — no squares, no products of unknowns, no curves.
The row picture reads A one row at a time: each equation is a line in 2D (a plane in 3D, a hyperplane beyond that), and the solution is the point where all of them meet. The column picture reads A one column at a time: each column is a vector, and Ax is a mixture of those columns with x saying how much of each to use. The solution is the mix that lands exactly on b.
Move the six coefficients. In the row picture the solution is where the two lines cross; in the column picture it is how much of each column to mix to land exactly on b.
The row picture is two constraints; the column picture is one target built from two ingredients. Same equation, different question.
The row picture solves two equations by substitution. The column picture solves the same problem by bookkeeping on vectors. Watch them land on the same point.
Numeric check of the answer in the original equations: 2·2 + 1 = 5 ✓ and 2 − 1 = 1 ✓. The row picture and the column picture are not two methods; they are two readings of the same matrix-vector product.
Before solving anything, ask whether the columns of A point in enough independent directions to reach b. The rank is the count of those directions, and it separates the three possible outcomes.
The rank of A is the number of linearly independent rows — which equals the number of independent columns. It is computed by running elimination and counting the pivots. Rank is what survives the row operations; it measures how many genuine directions the system can move in.
For a square n × n system, rank n means the columns are independent, the determinant is non-zero, and exactly one x reaches any b. Rank less than n means the columns are dependent: some direction of b may be unreachable, and any reachable point is reachable in infinitely many ways. The system is singular.
| What elimination finds | Rank | Solutions | Picture in 2D |
|---|---|---|---|
| n pivots, no contradiction | n (full) | exactly one | lines cross at a point |
| a pivot row [0 … 0 | non-zero] | < n | none — inconsistent | parallel lines |
| a zero row [0 … 0 | 0] | < n | infinitely many | the same line twice |
Run elimination on each system and read the last row. The same 2 × 2 coefficient matrix can produce any of the three outcomes depending only on the right-hand side.
Overdetermined systems (m > n, more equations than unknowns) are the normal case in machine learning: the fourth data point rarely lands on the line through the first three. Generic data puts b outside the column space, and that is precisely when least squares takes over.
Elimination is the workhorse: two legal row moves turn any square system into an upper-triangular one whose last equation has a single unknown. Pivoting is what keeps the arithmetic honest.
Two operations leave the solution set untouched: swap two rows, and subtract a multiple of one row from another. Gaussian elimination uses them to zero out every entry below the diagonal, column by column. What remains is an upper-triangular system, and back substitution solves it from the last equation upward — one unknown at a time.
The multiplier that clears entry (i, k) is m = A[i][k] / A[k][k]. The divisor A[k][k] is the pivot. If it is zero the algorithm stops; if it is tiny, every rounding error already present in that row is multiplied by a huge number. Partial pivoting fixes both problems by swapping the row with the largest absolute entry in the pivot column into position first. The cost is O(n³) either way — about a billion multiply–subtracts for a 1000 × 1000 system.
Step through row operations until the matrix is upper triangular, then watch back substitution fill in the solution. The toggle turns partial pivoting on and off.
Try the tiny-pivot system with pivoting off: dividing by 0.001 makes the multiplier 1000 and every later digit pays for it. Pivoting swaps first and keeps every multiplier at or below 1.
Why pivoting matters — a second worked example. Take A = [[0.0001, 1], [1, 1]] and b = [1, 2]. The exact solution is x₁ = 1/0.9999 ≈ 1.0001 and x₂ ≈ 0.9999. Simulate three-significant-digit arithmetic (a teaching model, not float64):
The update R2 − m·R1 multiplies whatever rounding error lives in R1 by m. Keeping m ≤ 1 — which is exactly what swapping in the largest pivot achieves — stops that amplification. In exact arithmetic the two orders agree perfectly; in finite arithmetic they can differ in every digit.
In the unpivoted 2×2 example, the multiplier was 10,000. What actually went wrong?
Elimination secretly factors A into two triangular matrices. Once those exist, every new right-hand side costs O(n²) instead of O(n³) — and symmetric positive definite matrices get a half-price version.
During elimination, the multipliers you compute are not throwaway scrap. Collect them below the diagonal, put 1s on the diagonal, and you have a lower-triangular matrix L. The eliminated result is an upper-triangular matrix U, and together they satisfy A = LU. With row swaps the identity becomes PA = LU, where P records the permutation.
Solving with the factors splits one hard problem into two easy ones. Write Ax = b as LUx = b, then let y = Ux: first solve Ly = b top-down (forward substitution), then Ux = y bottom-up (back substitution). Each triangular solve is O(n²), so the O(n³) factorization is paid once and reused by every future b.
When A is symmetric (A = Aᵀ) and positive definite (every eigenvalue > 0), it has a matrix square root: A = LLᵀ, the Cholesky decomposition. It needs only the lower triangle, costs about n³/3 operations instead of 2n³/3, and fails loudly — a negative number under the square root — if the matrix was not positive definite after all. Covariance matrices, Gaussian-process kernel matrices, and the ridge matrix XᵀX + λI are all symmetric positive definite.
Edit the symmetric matrix. Both panels factor the same A — LU keeps both triangles, Cholesky needs only the lower one and half the work.
Half the storage, half the arithmetic: Cholesky exploits the symmetry that LU ignores. The price is the requirement: symmetric, and every eigenvalue positive.
Gaussian processes use both facts: the predictive mean needs Kα = y, solved once from the Cholesky factor, and the marginal likelihood needs log det K, which the diagonal of L gives for free. Ridge regression does the same trick for (XᵀX + λI)w = Xᵀy.
Ridge regression solves (XᵀX + λI)w = Xᵀy with Cholesky rather than LU. What makes Cholesky legal here?
Every regression problem is overdetermined: more data points than parameters, so b sits outside the column space. Least squares finds the closest reachable point — and the answer is a familiar formula.
When A has m rows and n columns with m > n, the system is overdetermined. There is no x that satisfies every equation, but there is a best x: minimize the squared error ‖Ax − b‖², the sum of squared residuals. A residual is what one equation is still off by: rᵢ = (Ax)ᵢ − bᵢ.
Geometrically, Ax can only reach points in the column space of A. The closest such point to b is its perpendicular projection, and the leftover residual is perpendicular to the column space. That perpendicularity, written in coordinates, is the pair of equations AᵀA x = Aᵀb — the normal equations. For a data matrix with one column of 1s (the intercept) and one column per feature, this is the closed-form solution to linear regression.
Pull the four points up and down. The line minimizes the summed area of the orange squares — the squared residuals — and the fit always solves the normal equations.
Those two zeros are the normal equations in disguise: at the best line, the residual vector is perpendicular to every column of A.
The normal equations are the shortest path to the answer, but they have a numerical price: κ(AᵀA) = κ(A)², so forming AᵀA squares the condition number and can turn a mildly hard problem into an impossible one. Libraries fit with QR or SVD on A itself when features are nearly collinear. For a handful of well-scaled features — the textbook case — the normal equations are perfect.
At the least-squares solution, what is true about the residual vector r = Ax − b?
The condition number measures how loudly a small error in the input echoes in the answer. It turns “the computer said so” into a specific number of trustworthy digits.
The condition number of A is κ(A) = ‖A‖ · ‖A⁻¹‖ = σmax / σmin, the ratio of the largest to the smallest singular value. It is the worst-case amplification factor for relative error: disturb b by a relative 1% and x can move by up to κ%. A matrix with κ near 1 behaves like a rotation — safe. A matrix with κ near 10¹⁶ effectively destroys its smallest direction: in float64, whose relative precision is about 10⁻¹⁶, the answer carries no information at all.
The rule of thumb is one decimal digit per power of ten: with κ ≈ 10ᵏ you lose about k of the ~15 significant digits float64 stores, leaving 15 − k. Ill-conditioning in machine learning usually means nearly collinear features: two columns that say almost the same thing make σmin tiny. Regularization is the standard fix, because adding λI lifts σmin without inflating σmax.
The two rows of A are nearly identical, so the two lines nearly coincide. Nudge b a hair and watch the intersection shoot off the chart; add λ and it comes back.
κ is the worst-case amplification of a relative error in b. Each power of ten in κ costs about one decimal digit of the answer.
Perturb the right-hand side and see how the solution responds.
Numeric check. Let κ = 10⁸ and let b be stored with standard float64 rounding, a relative error of about 10⁻¹⁶. Then the solution’s relative error can reach 10⁸ × 10⁻¹⁶ = 10⁻⁸: eight digits lost, leaving about seven good ones.
You compute a solution with a matrix whose condition number is 10¹⁰, in float64 (about 15 digits). Roughly how many digits of the answer are meaningful?
A million unknowns make O(n³) hopeless, but many huge matrices are sparse: each equation touches only a few unknowns. Iterative methods exploit that by never eliminating anything at all.
Direct methods — elimination, LU, Cholesky — change the matrix until it is triangular. That costs O(n³) time and O(n²) storage, which is fine at n = 1,000 and impossible at n = 50,000,000. Iterative methods instead produce a sequence of guesses x₀, x₁, x₂, … that (hopefully) converges to the solution, using only matrix–vector products Ax. For a sparse A, each product costs O(non-zeros), not O(n²).
The simplest example is the Jacobi method: solve the i-th equation for the i-th unknown, plugging in the previous guesses for everything else. Everything depends on one number — the spectral radius ρ of the iteration matrix, the largest absolute eigenvalue. If ρ < 1 the error shrinks by a factor of ρ each sweep; if ρ ≥ 1 it stalls or explodes. Diagonally dominant matrices (and symmetric positive definite ones) guarantee ρ < 1. Conjugate gradient is the famous upgrade: for SPD systems it finds the exact answer in at most n steps in exact arithmetic, and usually far fewer — with convergence speed set by κ once more.
Instead of factoring A, each sweep solves every equation for its own unknown using the previous guess. Strong diagonal → fast convergence; weak diagonal → divergence.
Jacobi only needs A times a vector — no factorization. That is the door into conjugate gradient and the world of million-unknown systems.
That comparison is the whole story of large-scale solvers: the best method depends on how fast the error decays, and that rate is controlled by the spectrum — which is κ again, in disguise. Better conditioning means fewer sweeps, which is a third reason regularization helps.
Answer before you look. The pivoting and Cholesky questions are the ones that separate “I can call solve” from “I know what solve does”.
0 / 5 answered · 0 correct
01What does it mean geometrically when Ax = b has no exact solution?
02Why is partial pivoting used in Gaussian elimination?
03Why is Cholesky preferred over LU for solving (XᵀX + λI) w = Xᵀy in ridge regression?
04A matrix has condition number κ = 10⁸. Using float64 (about 15 significant digits), how many digits of the solution can you trust?
05What is the main advantage of LU decomposition over Gaussian elimination when solving Ax = b for many different b vectors?
Click a card to swap the lazy description for what it actually means.
Four problems, all fully worked — two of them correct arithmetic slips in the reference notebook. Try first; the answer is one click away.
Augmented [A | b]; R2 ← R2 − 4·R1 gives [0, −3, −6 | −9]; R3 ← R3 − 7·R1 gives [0, −6, −11 | −15]; R3 ← R3 − 2·R2 gives [0, 0, 1 | 3]. Back substitution: x₃ = 3; −3x₂ − 18 = −9 so x₂ = −3; x₁ − 6 + 9 = 6 so x₁ = 3. Check: 3 − 6 + 9 = 6 ✓, 12 − 15 + 18 = 15 ✓, 21 − 24 + 30 = 27 ✓. The multipliers 4, 7, 2 become L = [[1, 0, 0], [4, 1, 0], [7, 2, 1]] with U = [[1, 2, 3], [0, −3, −6], [0, 0, 1]], and L @ U reproduces A exactly. Warning: the reference notebook prints x = (1, 1, 1), which satisfies the first two equations but gives 7 + 8 + 10 = 25, not 27, in the third — the printed answer is a typo and x = (3, −3, 3) is the solution.
A has rows [1, x], so AᵀA = [[4, 10], [10, 30]] and Aᵀb = [22, 63]. The determinant is 4·30 − 100 = 20, giving c = (22·30 − 10·63)/20 = 30/20 = 1.5 and m = (4·63 − 10·22)/20 = 32/20 = 1.6. The fit is y = 1.5 + 1.6x, with predictions 3.1, 4.7, 6.3, 7.9 and residuals (fit − y) 0.1, −0.3, 0.3, −0.1. Checks: Σr = 0 and Σx·r = 1(0.1) + 2(−0.3) + 3(0.3) + 4(−0.1) = 0, so Aᵀr = 0 — the residual is perpendicular to both columns. SSE = 0.01 + 0.09 + 0.09 + 0.01 = 0.2. The reference prints 1.7, but c = 1.5, m = 1.7 gives Aᵀr = [1.0, 3.0] (not zero) and SSE = 0.5, so 1.6 is the true minimum; the printed value is an arithmetic slip.
Two nearly parallel rows mean a tiny determinant: det A = 10⁻⁸ and κ(A) ≈ 4/det = 4 × 10⁸ — about 8.6 digits are lost, leaving roughly 6 or 7 good digits in float64. Writing ε = 10⁻⁸ and δ = 10⁻⁸, the exact solution is x = (1 − δ/ε, δ/ε) = (0, 1): the second component moved by a full unit in response to a 10⁻⁸ change in b. With λ = 0.01: AᵀA ≈ [[2, 2], [2, 2]], so AᵀA + λI = [[2.01, 2], [2, 2.01]], Aᵀb = [2, 2], det = 2.01² − 4 = 0.0401, and x = (0.02/0.0401, 0.02/0.0401) ≈ (0.499, 0.499). The regularized answer accepts a slightly larger residual in exchange for being stable and small. In regression language, λ shrinks the weights toward zero and makes XᵀX + λI safely positive definite.
The exact solution is x₁ = 1/11 ≈ 0.0909 and x₂ = 7/11 ≈ 0.6364 (check: 4/11 + 7/11 = 1 and 1/11 + 21/11 = 2). Splitting A = D + (L + U) gives the iteration matrix M = −D⁻¹(L + U) = [[0, −1/4], [−1/3, 0]], whose eigenvalues are ±√(1/12) ≈ ±0.2887, so ρ ≈ 0.2887. From x₀ = (0, 0): sweep 1 gives (0.25, 0.6667), sweep 2 gives (0.0833, 0.5833), sweep 3 gives (0.1042, 0.6389) — oscillating toward the answer. The error shrinks like ρᵏ, so k ≥ ln(10⁻⁸)/ln(0.2887) ≈ 18.42/1.243 ≈ 14.8: 15 sweeps. Conjugate gradient solves this particular SPD 2 × 2 system in at most 2 steps, which is exactly why it is preferred when n is huge.
You do not need to master these here. Each one gets a proper treatment in its own lesson; the one-line meaning is enough to keep reading. Orange dotted underlines in the prose point back to this list.
This lesson is a port of an open course. Everything here traces back to it — and the next step is running the code yourself.
Lesson text adapted from AI Engineering from Scratch (Phase 01, Lesson 17) and the Math Foundations Notebook reference build. The row/column picture lab, elimination stepper, LU-vs-Cholesky console, least-squares lab, condition-number lab and Jacobi lab are original to this page, as are the corrected worked answers to exercises 1 and 2 and the exact squared form of the ridge condition number. Every displayed number is computed in your browser from the values shown. All labs run in your browser.