Solving linear systems by factoring the matrix once, inverting matrices, and fitting a curve to data that does not quite fit anything.
Chapter 28 is the numerical-methods chapter. Its centre is solving Ax = b for n equations in n unknowns, done not by the Gaussian elimination you were taught but by LU decomposition, which factors A once and then answers any number of right-hand sides cheaply. It then shows that matrix inversion and matrix multiplication are equivalent in difficulty, and closes with least-squares approximation, the workhorse of curve fitting and the ancestor of linear regression in Chapter 33.
Given a non-singular n × n matrix A and a vector b, find x with Ax = b. Two bad ideas first: Cramer’s rule is astronomically slow, and computing A⁻¹ then multiplying is both slower and numerically worse than the direct approach.
A = LU, where L is unit lower-triangular (ones on the diagonal, zeros above) and U is upper-triangular (zeros below the diagonal). Triangular systems are easy to solve, so the hard part is done once in the factorization and reused for every right-hand side.With A = LU, the equation Ax = b becomes LUx = b. Split it in two:
1. Solve Ly = b for y // forward substitution, top to bottom
2. Solve Ux = y for x // back substitution, bottom to topForward substitution works down from the first row: since L is unit lower-triangular, the first equation is y₁ = b₁, and each later row involves only already-known values.
yᵢ = bᵢ - ∑ from j=1 to i-1 of lᵢⱼ yⱼ
xᵢ = ( yᵢ - ∑ from j=i+1 to n of uᵢⱼ xⱼ ) / uᵢᵢEach substitution is Θ(n²) — a doubly nested loop. LU-SOLVE is therefore Θ(n²) given the factorization.
Ax = b for m different right-hand sides costs Θ(n³ + mn²), not Θ(mn³). In simulation and optimisation the same matrix is reused thousands of times with different data, so factoring once is the difference between feasible and hopeless.LU-DECOMPOSITION is Gaussian elimination with the multipliers recorded rather than discarded. At step k, use aₖₖ as the pivot, subtract multiples of row k from the rows below to zero out column k, and store the multiplier aᵢₖ/aₖₖ in L.
Schur complement step:
A = ┤ a₁₁ wᵀ ├ = ┤ 1 0 ├ · ┤ a₁₁ wᵀ ├
┤ v A′ ├ ┤ v/a₁₁ Iₙ₋₁ ├ ┤ 0 A′ - vwᵀ/a₁₁ ├
Then recurse on the Schur complement A′ - vwᵀ/a₁₁.Cost: Θ(n³), the same as matrix multiplication by the naive method, and CLRS notes the factorization can be done in place, overwriting A with both factors since L’s diagonal is known to be all ones.
aₖₖ is zero, the division is undefined and the algorithm stops — even when the matrix is non-singular and the system has a unique solution. And a pivot that is merely small is nearly as bad: dividing by it amplifies rounding error catastrophically, an effect called numerical instability.PA = LU. At each step, partial pivoting swaps in the row whose entry in the current column has the largest absolute value. That guarantees a non-zero pivot whenever one exists, and keeps every multiplier at most 1 in magnitude, which bounds error growth.Solving then becomes: permute b to get Pb, forward-substitute, back-substitute. Still Θ(n³) to factor and Θ(n²) to solve, and this — not plain LU — is what every numerical library actually implements.
To compute A⁻¹, solve Ax = eᵢ for each of the n standard basis vectors; the solutions are the columns of the inverse. With one LUP factorization reused n times, the cost is Θ(n³) + n · Θ(n²) = Θ(n³).
A⁻¹ and then x = A⁻¹b costs more than solving directly and is numerically worse, because forming the inverse introduces error that the multiplication then propagates. The rule in numerical computing is: never compute an inverse when a solve will do. CLRS presents inversion for its theoretical interest, not as a recommendation.A pleasing pair of theorems.
n × n matrix in I(n) time, where I(n) = Ω(n²) and satisfies a mild regularity condition, then you can multiply two n × n matrices in O(I(n)) time.The construction is elegant. Given A and B, build the 3n × 3n block matrix
D = ┤ Iₙ A 0 ├ ┤ Iₙ -A AB ├
┤ 0 Iₙ B ├ with D⁻¹ = ┤ 0 Iₙ -B ├
┤ 0 0 Iₙ ├ ┤ 0 0 Iₙ ├The product AB appears in the top-right block of the inverse. One inversion of a matrix three times the size yields the multiplication.
n × n matrices in M(n) time, you can invert a non-singular matrix in O(M(n)) time. So inversion and multiplication are asymptotically equivalent. Since Strassen multiplies in O(nlg7), matrices can be inverted in O(nlg7) too — and any future improvement to one immediately improves the other.Fitting a curve to data points that do not lie on any curve. Given m points (xᵢ, yᵢ), find a function from a chosen family — say a polynomial of degree n-1 with n < m — that comes closest.
With more equations than unknowns, Ac = y is overdetermined and generally has no exact solution. So minimise the error instead:
minimize ‖Ac - y‖² = ∑ from i=1 to m of ( f(xᵢ) - yᵢ )²c satisfiesAᵀA c = Aᵀy ⇒ c = (AᵀA)⁻¹ Aᵀ y = A⁺ yA⁺ = (AᵀA)⁻¹Aᵀ is the pseudoinverse. Note AᵀA is n × n and symmetric positive-definite when A has full column rank, so it can be factored and solved by the LU machinery above — the whole chapter connects here.
Why squared error? Because it is differentiable everywhere, which makes the minimum solvable in closed form, and because the geometry is clean: the residual Ac - y is orthogonal to the column space of A, so least squares is an orthogonal projection of y onto the space of achievable fits.
(AᵀA)⁻¹ explicitly — AᵀA squares the condition number, so QR or SVD factorizations are used instead. The theory is here; the numerics are a further course.A = LU into unit lower-triangular times upper-triangular. Θ(n³) once, then Θ(n²) per right-hand side.Ly = b, then back substitution on Ux = y.O(nlg7) transfers to inversion.‖Ac - y‖² for an overdetermined system, solved by the normal equations AᵀA c = Aᵀy.Chapter 29 generalises from equalities to inequalities. Linear programming optimises a linear objective subject to linear constraints, and it is the single most widely applicable optimisation framework there is — max flow, shortest paths, and the assignment problem are all special cases.