Every essay — page 2
Elimination, and the swap
Gaussian elimination is the first algorithm anybody learns and the row interchange is the part nobody is given a reason for. Here it is run without one, on a matrix where that fails — quietly, returning an answer of the right shape — and the growth factor that governs the whole story is plotted against a bound it never attains.
The swap that is not optional
Run elimination without a row interchange on a matrix that needs one and nothing announces a failure. There is no division by zero, no warning, and an answer of the right shape. It is simply wrong, and how wrong depends on a number you did not look at.
The bound that is never attained
Partial pivoting's stability guarantee permits the entries to double at every step — a factor of 5.5·10¹¹ at n = 40. The measured growth on random matrices of that size is about three. The gap is eleven orders of magnitude, and the guarantee is still worth having.
The pivot that reads the units
Partial pivoting compares the entries of a column and takes the largest. Those entries carry units, so the comparison depends on them — and there is a row scaling, on the standard two-by-two that pivoting exists to fix, which makes partial pivoting perform the identical catastrophic elimination it was introduced to prevent, with no interchange at all.
A factorisation with nothing to pivot for
Cholesky's growth factor is not bounded by one. It is equal to one, at every size and every condition number, and the two-line reason is why the algorithm needs no pivoting at all — not "usually gets away without it". Its only failure is the square root of a non-positive number, which is exactly the test for definiteness, and in floating point that test moves with the precision.
When symmetry is not enough
The matrix [[0, 1], [1, 0]] is symmetric, nonsingular and perfectly conditioned, and there is no diagonal entry to pivot on. Every factorisation restricted to symmetric interchanges and one-by-one pivots fails on it, at any depth of searching, because every entry it could search is zero. The repair is to take two variables at once.
A rule that is correct and unusable
Cramer's rule gives every component of the solution in closed form, in terms of determinants, and it is a theorem. On two-by-two systems whose rows are nearly parallel it returns an answer with a backward error of 458 units of roundoff where elimination returns 1.3 — on a matrix whose condition number is 32,000 and which elimination solved perfectly.
The inverse that is never formed
x = A⁻¹b is how the solution of a linear system is written and it is not how it is computed. The usual reason given is cost — three times the arithmetic. The real reason is that one of the two routes is backward stable and the other is not, and at κ = 10¹⁴ they differ by twelve orders of magnitude in the number that says whose fault a wrong answer is.
Orthogonality, measured
Orthogonal is not an adjective, it is a number: ‖QᵀQ − I‖. Two algorithms that are the same algebra written in a different order return 10⁻¹⁵ and 1 for it on the same matrix — and the one that fails still reconstructs the matrix perfectly, which is why nothing warns you.
Orthogonal is a number
"Q is orthogonal" is a claim about a measurable quantity, ‖QᵀQ − I‖, and on the eight-by-eight Hilbert matrix two standard algorithms return 10⁻¹⁵ and 1 for it. The one that returns 1 still reconstructs the matrix perfectly, which is why nothing warns you.
Two Gram–Schmidts
One argument changes. Classical Gram–Schmidt projects the original column onto each previous direction; modified projects what is left of it. In exact arithmetic the coefficients are identical. In floating point they differ by eight orders of magnitude in the thing that matters.
A reflection cannot stop being one
Householder QR holds orthogonality at 10⁻¹⁵ whatever the condition number of the matrix, and Gram–Schmidt does not. The reason is not that it is more careful. It is that its Q is built from unit vectors, and rounding a unit vector gives a different reflection rather than a broken one.
The nearest orthogonal matrix
Every field that has to clean up a drifted rotation reaches for QR, and QR does not answer the question. The nearest orthogonal matrix is the orthogonal factor of the polar decomposition — nearer by about a tenth, and, more to the point, the same matrix whatever order the columns were written in. QR's answer changes completely.
An iteration that only multiplies
Newton's iteration for the polar factor needs an inverse every step. Newton–Schulz needs only matrix products — nothing that reads an entry, nothing that pivots — and it converges if and only if every singular value is below √3. At 1.73205 it converges and at 1.73206 it returns an orthogonal matrix that is not the answer, with a residual of 5·10⁻¹⁶ and nothing to say so.
Least squares, and the road not to take
The normal equations are taught first and used by nobody, because forming AᵀA squares the condition number and then, below a computable value of ε, breaks outright. Underneath that is a harder fact: a wide range of very different fits explain the data equally well, and no arithmetic can choose between them.
The projection and the right angle
The least-squares solution is the one whose residual is perpendicular to everything the columns can reach. That is not a mnemonic — it is an equation, Aᵀr = 0, and the computed answer satisfies it to 10⁻¹⁶.
The road that squares the problem
The normal equations are the first method every course teaches and the method no library uses. Forming AᵀA squares the condition number, and below ε = √u it does not degrade — it produces a matrix that is exactly singular, from data that was perfectly usable.
The valley with no bottom
A degree-nine fit's coefficients can be moved by a third of their own size before the residual changes in the sixth significant figure. The arithmetic did not lose those digits. The data never contained them.
When the matrix is wrong too
Every least-squares problem on this site has assumed A is exact and b is not, and moved b onto the column space of A. Where both were measured, the smallest correction that makes the system consistent moves the matrix as well — and on the problems where that answer is more accurate, it has the larger residual, by construction rather than by luck.
A correction cheaper than the problem
Sherman and Morrison's formula updates a solved system for a rank-one change to the matrix, at 4n² operations instead of (2/3)n³. It is exact algebra. On a problem whose updated matrix is the identity — condition number one, the easiest system there is — it returns a forward error of 2.5·10⁻⁴ where a direct solve returns 10⁻¹⁶.
The observation that cannot be removed
Removing a rank-one term from a Cholesky factor needs a rotation that is not orthogonal, and the number under its square root is 1 − h, where h is the leverage of the row being removed. The algorithm's breakdown condition and the statistician's warning are the same quantity, arrived at from opposite ends, and neither field states it in the other's language.
Regularisation, and the answer that is chosen
Some problems do not determine their own answer. The data is consistent with a range of solutions that differ by orders of magnitude, no arithmetic can choose between them, and something outside the data has to. That choice is the computation rather than a preliminary to it — and every published rule for making it is a heuristic scored, here, against a truth that exists only because the problem was constructed.
When the answer is a choice
A backward-stable least-squares solve of this problem returns an answer whose relative error is 5.5·10⁸. Nothing went wrong. The singular values decay exponentially with no gap anywhere in them, the data does not determine the answer, and something outside the data has to choose — which is the computation rather than a preliminary to it.
Where the answer stops being in the data
The Picard condition finds the index where a noisy right-hand side stops carrying signal, from the data alone, with no knowledge of the answer. It lands at 32 where the truncation that actually minimises the error is 28 — and at 45 where the best is 38. It overshoots at every noise level from 1% to 0.001%, and it overshoots for a reason.
Choosing without knowing
Three published rules for choosing a regularisation parameter, scored against an oracle that requires the exact answer and is therefore not a method. Generalised cross-validation lands on the oracle's λ exactly; the discrepancy principle costs 6%; the L-curve costs 129%. And told a noise level ten times too small, the discrepancy principle's error goes from 0.112 to 10,449.
The basis decides what a filter is
The vocabulary of regularisation is spectral — a method keeps a component or discards it, and the weights are a function of the singular value. Row-normalising a symmetric blur so that it preserves a constant makes it 8.6% asymmetric, and that is enough to move GMRES's weights from 7·10⁻¹⁴ off a function of σ to 4.4·10⁻².
Eigenvalues, singular values, rank
A symmetric matrix gives up its eigenvalues without complaint. An unsymmetric one can move them by the eighth root of a perturbation, so rounding error alone shifts them by a hundredth. And rank is not a property a floating-point matrix has — it is a decision about a gap, and the gap is worth printing beside it. The algorithm that libraries actually run lives here too, and it needs a shift before it is an algorithm at all.
Symmetry is worth more than precision
A symmetric matrix gives up its eigenvalues to full accuracy however ill-conditioned it is. An unsymmetric one can move them by the eighth root of a perturbation, so the rounding involved in merely storing the matrix shifts the spectrum by a hundredth.
Rank is a decision
A floating-point matrix does not have a rank. It has a spectrum of singular values, and somewhere in that spectrum is a place where the values stop being signal and start being noise. Deciding where is a judgement, and the evidence for it is a gap.