Conditioning, stability, and cost
Every result so far has been exact: a matrix either has an inverse or it does not, a system either has a solution or it does not. Real computation does not get that luxury. Entries are stored to a finite number of digits, the right-hand side is measured rather than dictated, and the answer is only as trustworthy as the ratio of the largest stretch to the smallest. This part is about that ratio. It is called the condition number, and it decides both how wrong a solution can be and which decomposition you should pay for.
The question
When the data wobbles, how much does the answer wobble?
Solving A x = b is a map in two directions. Forward, A takes an input x and produces an output b; if A squashes some direction almost to nothing, then different inputs can produce nearly the same output, and the map is hard to read backwards. Backward, we invert: x = A−1b. The same squashing now becomes a huge stretching, so a small error in b arrives at x magnified. Conditioning is the size of that magnification.
This is not a defect of any particular algorithm. Even if you could solve the system exactly, an uncertainty of one part in a thousand in b would still come out as an uncertainty in x that the matrix controls, not the method. What an algorithm can and cannot change is the extra error it introduces on top of that floor. A stable algorithm leaves the floor alone; an unstable route, or a careless rearrangement like the normal equations, raises the floor by orders of magnitude. The two ideas — how badly conditioned the problem is, and how faithfully a method solves it — are the whole of numerical linear algebra, and keeping them apart is the point of this part.
A small change with a large echo
Watch one perturbation get amplified
Here is the effect in a single picture. The faint ellipse is the image of the unit circle under A−1: the set of all answers x = A−1b as b runs around the unit circle. Its long semi-axis is 1/σmin and its short semi-axis is 1/σmax, so the more lopsided the ellipse, the more the inverse can amplify. The two dashed spokes are those axes. Drag the conditioning slider toward zero and watch the ellipse stretch into a needle: the matrix is collapsing one direction in the output, and the inverse is doing the opposite.
The two solid arrows make the warning concrete. Start with a target b, solve for the exact x, then nudge b by a small amount Δb and solve again. The move from x to x + Δx is the dashed arrow, and the readout compares its relative size to the relative size of the nudge. The perturbation is deliberately aimed along the direction the matrix squeezes hardest, which is where the amplification is worst. Edit the entries of A if you like; the picture recomputes the ellipse and the answers from scratch.
A small change Δb (aimed along the weakest direction) produces the dashed jump Δx. The ellipse is the unit circle seen through A−1.
The rule the demo is checking is the standard error bound. If Δb is the only uncertainty, then up to a term that vanishes when Δb is truly small,
Push the conditioning slider to its extreme and the observed amplification in the readout approaches the condition number. That is the sense in which κ(A) is the worst case: for most perturbations you get less, but for one special direction Δb you get exactly the bound. There is a mirror bound for a perturbation of the matrix itself, and it scales the same way.
The condition number
One ratio, read straight off the SVD
The ellipses of Part 16 already contain the answer. Every matrix sends the unit circle to an ellipse whose semi-axes are the singular values σ1 ≥ σ2 ≥ … ≥ σn ≥ 0, and those numbers describe the largest and smallest stretches the matrix can apply. The condition number is their ratio:
Read it as a statement about directions. The matrix can stretch some input direction by σmax and another by only σmin, so the inverse must divide by σmin and therefore multiplies by at least 1/σmin. A condition number of one is the best possible: rotations and reflections change nothing about lengths, so their ellipses are circles and the inverse is as tame as the map. The identity has κ = 1. As σmin → 0 the ellipse flattens and κ → ∞ exactly as A becomes singular.
Three properties make κ the right summary. It is scale invariant, because doubling every entry of A doubles the numerator and the denominator together. It is the exact worst-case amplification above, so it is a fact about the matrix and not about any solver. And it is the number of digits you lose. In floating point with about sixteen decimal digits, a solve whose condition number is 10k can surrender up to k of them, leaving roughly 16 − k trustworthy digits in x. A perfectly stable algorithm still pays that toll; a careless one can pay twice.
That distinction is worth naming precisely, because the two failure modes get blamed on each other. Backward error asks how much you would have to change the input to make your computed answer exact for the changed input. A stable algorithm keeps that change at the level of rounding, regardless of κ. Forward error asks how far your answer is from the true one, and it is bounded by the backward error times the condition number. So the condition number is a property of the problem, stability is a property of the algorithm, and the forward error is the product. QR and LU are backward stable; forming AᵀA is not, and the next demo shows what it costs.
Squared by the normal equations
The one rearrangement you should never make
When A x = b has no exact solution, Part 10 turned the problem into a projection and wrote down the normal equations AᵀA x = Aᵀb. Algebraically they are perfect. Numerically they are a trap, and the reason is one line from the SVD: squaring the matrix squares the singular values, so the condition number squares too.
If A has κ = 104 — a mild loss of about four digits on its own — then AᵀA has κ = 108 before you have solved anything at all. The QR route of Part 9 never forms that product. It factors A = Q R with Q orthogonal, solves R x = Qᵀb, and keeps the condition number at κ(A) instead of κ(A)2.
The chart below checks the claim numerically. It builds matrices whose conditioning exponent p means κ(A) = 10p, then measures the relative error of two solves: solve(A, b) against the normal equations solve(AᵀA, Aᵀb). On log–log axes a straight line means a power law, and the slope is the exponent. The stable line climbs with slope one; the normal-equations line climbs with slope two, then leaves the top of the plot. The vertical axis is the true computed error, so its position is not fixed — but its slope is, and the slope is the lesson.
Relative error of x against the conditioning exponent p, where κ(A) = 10p. Slope one is the stable solve; slope two is the normal equations.
Two practical corollaries follow. First, if you must use the normal equations because AᵀA is all you can afford to form — which happens constantly in large sparse problems — then at least solve the resulting symmetric positive-definite system with Cholesky and never invert it. You have already paid the condition-number tax; do not pay an instability tax as well. Second, the same squaring is why ridge regularisation works: adding λ2I lifts every squared singular value by λ2, capping the condition number at σmax2/λ2 and buying back the digits at the cost of a little bias, as Part 18 showed.
What each decomposition costs
Accuracy is not free, and neither is robustness
Conditioning tells you which algorithms are trustworthy. Cost tells you which ones you can afford. For an n × n dense matrix the leading term of every classical method is proportional to n3, and the differences live entirely in the constant. LU elimination is about 2n3/3 operations; if the matrix is symmetric positive definite, Cholesky halves that to n3/3. A Householder QR is about 2n3. A symmetric eigendecomposition runs to roughly 9n3 once eigenvectors are included, and a full SVD to something like 12n3, which is why you compute one only when the singular values themselves are the point.
Those constants are why the choice of method is a choice about structure, not just speed. Drag the slider to change n and read the counts off the chart. Because everything scales as n3, the lines are parallel on log–log axes and their vertical gaps are the ratios of the constants: Cholesky is roughly thirty-six times cheaper than a full SVD, QR about three times more expensive than LU. Structure is the only thing that breaks the pattern; a banded or sparse matrix can drop below n3 entirely, which is exactly why the applications below care.
Leading-order operation counts as functions of n = 2p, on log–log axes. The marker follows the slider.
| Decomposition | Use it when | Leading cost | Watch out for |
|---|---|---|---|
| Cholesky | A is symmetric positive definite — including a well-conditioned AᵀA | n3/3 | Fails at once if A is not positive definite; squares κ when applied to a Gram matrix. |
| LU (with pivoting) | A general square nonsingular system, or repeated solves with the same A | 2n3/3 | Partial pivoting is essential; near-singular pivots are the visible symptom of conditioning. |
| QR | Rectangular or least-squares problems; any time you would otherwise form AᵀA | 2n3 | Keeps the condition number of A, not its square — the safe default for least squares. |
| Symmetric eigen | Quadratic forms, covariance matrices, PCA, vibration modes | ~9n3 | Only for symmetric matrices; eigenvectors cost more than eigenvalues alone. |
| SVD | Conditioning, rank, pseudoinverse, low-rank approximation — when the singular values are wanted | ~12n3 | The most robust and the most expensive; use one-sided iteration near singularity, never eig(AᵀA). |
Costs are leading-order flop counts for dense square matrices and the constants vary with the implementation; the ordering, not the third digit, is the thing to remember.
Where this shows up
One ratio, two worlds
Why bundle adjustment avoids the normal equations
A structure-from-motion problem has one camera block that is nearly singular whenever a camera has little parallax, so the reduced system JᵀJ is exactly the squared-condition-number trap. Practical solvers either factor the sparse Jacobian with QR or run Cholesky on the reduced system with careful damping, and they always keep a trust region so the step stays finite. The pose-graph part of the optimization guide works through the same reduced system and its gauge freedom.
Mixed precision and ill-conditioned layers
Low-precision matrix multiplies are only safe while the matrices stay well-conditioned; when a weight matrix or an attention block is nearly low-rank, the rounding that mixed precision tolerates quietly becomes a large forward error, and the loss curve spikes. That is why large training runs keep sensitive reductions in higher precision and monitor scale factors, and why quantisation of an ill-conditioned layer is where accuracy is lost first — see training at scale and quantisation.
Further reading
- Grant Sanderson, "Singular value decomposition", Essence of Linear Algebra, 3Blue1Brown — the ellipse whose axes become σmax and σmin, and therefore κ(A).
- Gilbert Strang, 18.06 Linear Algebra, MIT OpenCourseWare — Lecture 12 on graphs and networks, and Lectures 29–30 on the SVD and its numerical behaviour.
- Lloyd N. Trefethen and David Bau III, Numerical Linear Algebra, Lectures 12 and 15 — condition numbers and the conditioning of least-squares problems, including the normal-equations warning.
- Cleve Moler, Numerical Computing with MATLAB, chapter 2 — a short, practical account of roundoff, pivoting and why the condition number predicts lost digits.