An orthogonalisation nobody calls one
Every derivation of conjugate gradients ends with the same corollary, and it is always stated as a bonus rather than as a claim requiring evidence.
The residual vectors the method produces are mutually orthogonal. There are n of them at most in an n-dimensional space, so after n steps the Krylov subspace is everything, and the iterate that minimises the error over everything is the answer. Conjugate gradients terminate in at most n steps. The method is therefore a direct method that happens to be usable early.
That corollary is correct, it is provable, and it describes an algorithm nobody has ever run.
The measurement is one the site already has
This is the part worth noticing before any number appears. Conjugate gradients are not usually described as an orthogonalisation at all: they are derived by minimising a quadratic, or by constructing A-conjugate search directions, and the orthogonality of the residuals arrives as an identity along the way.
But if the residuals are supposed to be mutually orthogonal, then collecting them as the columns of a
matrix R, normalising, and computing ‖RᵀR − I‖ is exactly the measurement
orthogonal is a number is about. No new machinery is required.
orthError from the site’s matrix library, the same function that reports 1.4 for classical
Gram–Schmidt on the Hilbert matrix, applies unchanged to a basis nobody built deliberately.
That is the whole essay in one sentence: the same failure, in a method that does not look like the method it fails the same way as.
What the numbers are
On the site’s model problem — the one-dimensional discrete Laplacian at n = 30, condition number 178 — running to exactly n steps leaves the residual basis with ‖RᵀR − I‖ = 1.41.
For scale: a set of vectors that are exactly orthonormal gives 0, and a set where one pair is fully parallel gives about 1. So after thirty steps of a method whose derivation asserts mutual orthogonality, the residuals are not approximately orthogonal, not slightly degraded — they are as far from orthogonal as an arbitrary set of unit vectors would be.
Raise the condition number and it is worse. On a matrix built with κ = 10⁶ at the same size, the figure is 5.29, and the residual after all thirty steps is 5.1·10⁻⁴ — not converged, nowhere near converged, with the guarantee of exactness comprehensively spent.
Why this is the Gram–Schmidt failure and not a new one
The mechanism transfers exactly, which is what makes it worth putting the two essays beside each other.
Classical Gram–Schmidt computes a projection coefficient against the original column, and the argument that this is equivalent to using the running remainder assumes the previously computed vectors are exactly orthogonal. They are not, so the terms that should vanish contribute a little, and what they contribute is multiplied by something large. That is two Gram–Schmidts in three sentences.
Conjugate gradients do the same thing with a shorter memory. The recurrence orthogonalises the new residual against only the previous one and the one before, and it is entitled to skip the rest because the earlier residuals are already orthogonal to it. When they are not, the terms it declined to compute are not zero, and nothing in the recurrence has any way to notice.
The short recurrence is the method’s whole advantage over storing everything, and it is bought against precisely the invariant that rounding destroys.
The derivation assumed an invariant, which is the general shape
Two Gram–Schmidts states the pattern in a form that survives outside numerical analysis, and this is its cleanest additional instance:
A derivation that assumes an invariant holds exactly cannot see the cost of the invariant holding approximately.
Here the invariant is orthogonality of the residuals, the derivation is the conjugacy argument, and the cost is the entire finite-termination property. Not a degraded version of it — the property is simply absent, and its absence is not signalled anywhere in the computation.
The same shape appears in the swap that is not optional, where correctness rests on the pivot not being tiny, and in the road that squares the problem, where it rests on AᵀA being formable. Each derivation is valid; each assumption is quietly false; and in each case what follows is not small.
Where the model problem is not enough, and why that is said out loud
The first version of the site’s check for this used the two-dimensional Laplacian on a 12×12 grid — a 144-dimensional problem, run for 144 steps — and asserted that the residual would not be at rounding level. It failed, and the failure was correct.
That matrix has κ ≈ 50. Conjugate gradients reach machine precision on it around step 30, so by step 144 the method has been converged for a hundred iterations and the residual is 10⁻⁵⁴. Finite termination was never tested, because convergence got there first.
The property is only under pressure when the condition number is large enough that n steps are not enough on their own, and the check now uses a matrix built at κ = 10⁶ for exactly that reason. The substitution is recorded in the library rather than made silently, because “a different matrix was chosen” and “a matrix was chosen on which the claim is testable” describe the same action with very different standing, and only the second one is evidence.
The general form of that is worth keeping. A test that passes because the condition it tests for never arose is indistinguishable, in a green build, from a test that passes because the property holds. The only defence is to check that the test can fail, which is why every claim on this site is paired with a case it must refuse — and why the refusal here is applied to a matrix chosen to make refusal possible.
There are two orthogonalities, and both of them go
The essay has been treating conjugate gradients as having one orthogonality property. It has two, they are different properties, and separating them sharpens what is lost.
The residuals are mutually orthogonal, in the ordinary sense: rᵢᵀrⱼ = 0 for i ≠ j. This is what the figure measures and what makes the Krylov subspace argument work.
The search directions are A-conjugate, which is orthogonality in the inner product the matrix defines: pᵢᵀApⱼ = 0 for i ≠ j. This is the property the method is named for, and it is what makes the one-dimensional minimisation along each direction not spoil the minimisation already achieved along the previous ones. Without it, stepping along p₂ would undo part of the progress made along p₁ and the method would be ordinary steepest descent, which converges at a rate governed by κ rather than √κ — the difference between five hundred iterations and twenty.
Both are exact identities in exact arithmetic and both are computed by the same short recurrence, so both degrade together and for the same reason. The second is arguably the more consequential: losing residual orthogonality costs a termination guarantee nobody used, while losing conjugacy costs the convergence rate that is the method’s entire justification.
That it does not cost the convergence rate — the bound holds throughout, checked at every step — is the reconciliation the last section of this essay is about, and it is a genuinely fortunate accident rather than an obvious consequence.
What survives
It would be easy to read the above as saying conjugate gradients are broken, and they are not — they are among the most reliable methods in the subject. The precise statement is narrower and worth getting right.
Finite termination is gone. Nobody relies on it and no implementation is written to use it, so the practical loss is nil. Its real role was pedagogical, and losing it is mostly a correction to how the method is described.
The convergence bound survives. The κ bound is checked at every step in the rate the condition number predicts and it holds everywhere, with a closest approach of 43% — on runs whose residual bases have ‖RᵀR − I‖ above
- The method converges at the rate the theory gives while comprehensively violating the theory’s hypothesis.
And the delay is real. What loss of orthogonality actually costs is that the method re-explores directions it has already covered, so convergence is delayed rather than prevented. In finite precision the number of steps to a given accuracy can be several times what exact arithmetic would need, and it is still finite and still bounded.
What is asserted here
The residual basis has lost orthogonality by step n, at every position of the slider, measured
with the same orthError the factorisation essays use.
And it lost it progressively, which distinguishes a genuine drift from a single catastrophic step and is the claim the plot is actually making.
After n steps at κ = 10⁶ the residual is nowhere near rounding, at 5.1·10⁻⁴, with the residual basis at 5.29.
And the refusal: an assertion that the method terminates in n steps is fed the measured residual and must throw. It does. Without that, the entire essay would rest on a claim that had never been put in a position to fail.
Selective re-orthogonalisation, which is the practical answer
The fix in a serious code is not to abandon short recurrences — that would give up the method’s reason for existing — but to detect the loss and repair it when it matters.
The test is cheap: compare the norm of a vector before and after orthogonalisation, and if it dropped by more than about a factor of two, the subtraction cancelled and the result is untrustworthy. Re-orthogonalise only then. That is what “selective re-orthogonalisation” means, and it costs almost nothing on the steps where nothing has gone wrong.
The theorem underneath it is Kahan and Parlett’s twice is enough, and this site measures it rather than quoting it: on the eight-by-eight Hilbert matrix, one classical Gram–Schmidt pass gives ‖QᵀQ − I‖ = 1.44 and two passes give 1.09·10⁻¹². Twelve orders of magnitude for one extra pass — and still a thousand times above where Householder sits, which is the honest qualification and is in two Gram–Schmidts with the rest of the table.
For Lanczos and Arnoldi — the eigenvalue and non-symmetric-solve methods built on the same recurrence — this is not optional. Loss of orthogonality there produces ghost eigenvalues: copies of converged eigenvalues appearing repeatedly, because the basis has drifted back into a direction it already covered and the method has no way to know it has been there before.
The ghost is worth dwelling on, because it is the same phenomenon as this essay’s subject wearing a much more alarming costume. A Lanczos run asked for the ten largest eigenvalues of a matrix can return the largest one four times. Nothing has gone wrong arithmetically — each of those four values is an accurate eigenvalue of the matrix — and nothing in the output distinguishes a genuine multiplicity of four from a single eigenvalue found four times. The failure is undetectable from the answer, which is the property this whole site is organised around, and the only defence is to measure the orthogonality of the basis while the run is happening rather than to inspect the result afterwards.
Which is the practical moral of the essay. The residual is visible and orthogonality is not, so a code that watches only the residual is watching the quantity that stays healthy while the one that matters degrades. Measuring ‖RᵀR − I‖ costs a few inner products against a basis already in memory, and it is the only warning available.
The cost of doing it properly
There is an alternative to the short recurrence, and knowing what it costs is what makes the recurrence’s fragility a bargain rather than a defect.
Full orthogonalisation stores every basis vector and orthogonalises each new one against all of them. Orthogonality is then maintained to rounding level, the termination property survives, and the method is completely reliable. It is also O(k) work and O(kn) storage at step k, so the cost of the k-th step grows without limit and the memory grows with it. On a problem large enough to need an iterative method in the first place, storing two hundred vectors of length ten million is not available.
That is exactly the trade the non-symmetric case is forced into. GMRES has no short recurrence available — symmetry is what supplies one — so it must either store the whole basis or restart and give up its optimality, which is why the spectrum that predicts nothing is about a method with a memory problem as much as a convergence one.
Conjugate gradients get the short recurrence for free from symmetry, pay for it in orthogonality that rounding erodes, and the erosion turns out to cost a delay rather than a failure. Stated as a decision it is not close: three vectors and a delay, against unbounded storage and no delay.
The thing that is genuinely surprising
Not that rounding degrades an orthogonality property — that is the ordinary behaviour of the whole subject and this site has eleven essays about it.
What is surprising is that the degradation is total and the method is unaffected in the way that matters. ‖RᵀR − I‖ = 5.29 is not a small perturbation of 0. It means the residual vectors span fewer directions than they are counted as spanning, that the theoretical basis for the recurrence is gone, and that every step after the drift began is computing coefficients against an assumption that is false.
And the error still falls, still under the bound, at a rate still governed by κ.
The reconciliation is that conjugate gradients have two justifications and only one of them was destroyed. The finite-termination argument needs exact orthogonality and is gone. The optimality argument — that x_k minimises the error over the Krylov subspace built so far — degrades gracefully, because a slightly wrong subspace still contains most of what the right one contained. The method kept the property that made it useful and lost the one that made it elegant.
That is a distinction worth carrying beyond this method: when an algorithm’s guarantee fails, ask which guarantee. A proof of exactness and a proof of near-optimality fail very differently under perturbation, and the literature tends to present them together as though they stood or fell as one.