Two Gram–Schmidts

Gram–Schmidt orthogonalisation is a construction anyone can reproduce from the idea. Take the columns of a matrix one at a time; from each, subtract its projection onto every direction already established; normalise what is left. The result is an orthonormal set spanning the same space, and the coefficients subtracted along the way assemble into the R of a QR factorisation.

Write it as code and there is one place where the description is ambiguous, and the ambiguity does not look like one.

                          classical                     modified
for j in columns:
    v = a[j]                                            v = a[j]
    for i < j:
        r[i][j] = q[i] · a[j]      ← the ORIGINAL       r[i][j] = q[i] · v      ← what is LEFT
        v = v − r[i][j] · q[i]                          v = v − r[i][j] · q[i]
    r[j][j] = ‖v‖                                       r[j][j] = ‖v‖
    q[j] = v / r[j][j]                                  q[j] = v / r[j][j]

One argument. Everything else is identical, including the operation count.

Why they are the same in exact arithmetic #

In exact arithmetic, qᵢ · aⱼ and qᵢ · v are equal, and the reason is the property being constructed.

By the time the loop reaches i, the running vector v differs from aⱼ by a combination of q₁ … qᵢ₋₁. Taking the inner product with qᵢ, every one of those terms contributes zero, because qᵢ is orthogonal to all of them. So the two expressions give the same number, and the two algorithms produce the same output.

The argument is correct and it is self-referential, which is the crack the arithmetic gets into. It assumes the previous q’s are exactly orthogonal to each other. They are not — they are the output of earlier steps of the same process — and so the terms that should contribute zero contribute a little, and the two coefficients differ.

Why the difference matters so much #

The difference between the two coefficients is small. What it multiplies is not.

In the classical version, the coefficient r[i][j] is computed against a vector that still contains the components being removed. Any error in qᵢ therefore enters r[i][j] multiplied by the full size of aⱼ. In the modified version, it enters multiplied by the size of what remains of aⱼ after previous subtractions — which for a nearly dependent column is very much smaller.

That is the whole mechanism, and it produces exactly the observed rates: the classical error is amplified by κ twice, giving κ²u, and the modified error is amplified once, giving κu.

The R factors agree, which is the surprising part #

The natural assumption is that the unstable version is unstable throughout. It is not. Both versions produce essentially the same R.

On the eight-by-eight Hilbert matrix, ‖R_classical − R_modified‖/‖R‖ is 1.2·10⁻⁶ — small, and specifically of the order κu rather than of the order of the orthogonality loss, which is 1.

So the two algorithms differ in the Q and agree in the R, and the algebra they share is visible in that agreement. This is not a case of one implementation being broken. It is a case of two correct implementations of the same derivation, differing in one respect that the derivation cannot see.

The practical consequence is unpleasant. Anything that only uses R — computing a determinant, say, or a condition estimate — is unaffected by the choice. Anything that uses Q is affected completely. So the bug can sit in a codebase for years, doing no damage, until somebody uses the other output.

Every one of these numbers has a ceiling, and they hit it in order #

Three laws have now been quoted — κ²u for the classical loss, κu for the modified loss, κu for the disagreement between the R factors — and each was quoted at a single condition number. A law read off one point is not a law; any number of other laws pass through the same point. So the three were measured across ten decades of κ, taking the median over nine random matrices of prescribed conditioning at each, 12×12, in double:

The three normalised columns are what to read. Down to κ = 10⁸ they are constants: κ²u to within a factor of 2.9, and both κu columns to within a factor of 2.3 and 3.3. Seven and eight decades is enough to call those laws confirmed rather than quoted, and the ‖ΔR‖ column in particular is now a measurement rather than one Hilbert matrix.

Then two of the three fall off a cliff, and that is the part worth carrying away.

‖QᵀQ − I‖ = 1 is the ceiling of the instrument, not a measurement of the damage. The quantity cannot exceed about √n, because that is what ‖QᵀQ − I‖ is when Q’s columns are arbitrary unit vectors, and classical Gram–Schmidt has reached it by κ = 10¹⁰. Its normalised column collapses by four orders and then by eight — not because the algorithm improved, but because the number stopped being able to grow. The 1 quoted at the top of this essay is the same 1 that would appear at κ = 10¹⁵ and at κ = 10²⁰. It says the loss ran out of room and nothing else.

The R disagreement has a ceiling too, one decade higher. At κ = 10¹² it is 1.8·10⁻⁷ against 1.7·10⁻⁷ at 10¹⁰ — it has stopped growing, and for a reason that is the mirror of the section above: two algorithms cannot disagree about columns that are numerically dependent. There is nothing left in the later columns for the ordering of the subtractions to change.

Modified Gram–Schmidt is the only one of the three still on its law at κ = 10¹². That is the sharpest statement of why it is the one to use, and it is not the statement usually made. The usual one — classical loses orthogonality, modified keeps it — is true at moderate κ and becomes vacuous at large κ, where the classical number is pinned and no comparison of magnitudes means anything. The statement that survives is about range: at every condition number tested there is still information in the modified figure, and there is none in the classical one past 10⁹.

And the order in which the ceilings arrive is what makes the section above true rather than lucky. Classical orthogonality saturates first, the R factors second, modified orthogonality not at all in this range. That ordering — and only that ordering — is what allows a code to have lost Q completely while R is still worth using, which is the situation the whole essay is about. Reverse it and the bug would announce itself in the determinant.

assertTheThreeLawsAllHaveACeiling in the matrix library measures all three columns and requires each law to be flat on its own domain and each ceiling to be present at the top. This is the same discipline as the residual the method reports: there the printed number is computed correctly and is not the quantity its name claims, and here the reported quantity is real and the thing a reader takes from it is not.

What “loss of orthogonality” costs downstream #

It is worth being concrete about the consequences, because ‖QᵀQ − I‖ = 1 sounds abstract.

Least squares by QR computes x from Rx = Qᵀb. That derivation assumes QQᵀ projects onto the column space, which requires orthogonality. With classical Gram–Schmidt on an ill-conditioned matrix, the computed x can be wrong by a factor of κ more than it needs to be — the method inherits the normal-equations behaviour it was chosen to avoid.

Krylov subspace methods build an orthonormal basis one vector at a time and are the reason Gram–Schmidt is still used at all — Householder does not fit their access pattern. Loss of orthogonality in the Arnoldi basis is the reason GMRES stalls and the reason Lanczos produces spurious duplicate eigenvalues, and re-orthogonalisation is the standard defence.

Anything that treats Q as a change of basis — an eigenvalue deflation, a subspace projection — silently uses a basis whose vectors are not independent. The failure surfaces far from its cause.

The remedies, in order of cost #

Modified Gram–Schmidt. Free. One word. Turns κ² into κ, which for κ up to about 10⁷ is enough in double precision, and it is the version anyone writing Gram–Schmidt from scratch should write.

Re-orthogonalisation. Run the inner loop twice. The second pass has almost nothing to subtract, and the classical result is that this restores orthogonality to rounding level for any Q that was not already hopeless — “twice is enough”, which is a theorem rather than a rule of thumb. It doubles the cost of that step, and it is what serious Krylov implementations do.

Householder. Roughly twice the arithmetic of Gram–Schmidt for a square matrix, and orthogonality that does not depend on the matrix at all. A reflection cannot stop being one is why. For a dense factorisation this is simply the right answer and is what every library does.

Notice what is not on the list. Precision does not fix it: the slider on the loss curve moves the classical line down without changing its slope, so a wider float buys a factor and the algorithm choice buys a rate. That is the same conclusion as the order they are added in reaches about summation, and it recurs often enough on this site to be worth treating as a general rule.

The historical note, which is unusually apt #

Modified Gram–Schmidt is older than classical Gram–Schmidt.

The algorithm now called classical appears in Gram (1883) and Schmidt (1907). The modified version appears in Laplace’s work in 1816, nearly a century earlier, and in Cauchy after him. It was rediscovered as an improvement on a method that came later.

What makes the story apt rather than merely amusing is why it happened. The modified form arises naturally when the process is done by hand, because a human computer subtracts as they go and works with the vector in front of them. The classical form arises naturally when the process is written as a formula, because a formula refers to the original column by name. The stable algorithm is the one that comes from doing the work; the unstable one comes from writing the work down.

That is not a general law, but it is a pattern worth noticing, and it appears again in the swap that is not optional, where the stable choice is the one a person eliminating by hand would make without being told.

What is asserted here #

The essay’s central pair of claims is checked on every build, and both directions can fail.

The classical version must lose orthogonality: the largest pairwise inner product among its columns must exceed the modified version’s by a factor of a hundred. If some change to the arithmetic made classical Gram–Schmidt behave, this essay would be false and the build would stop.

The R factors must agree, to within fifty times κu — an inequality scaled to the matrix rather than a round number, because the agreement is at the level κu and asserting 10⁻¹⁵ would be asserting something untrue.

Every classical column must still be a unit vector, to 10⁻¹², which is the claim that makes the failure quiet.

And Q must differ by far more than R does: the worst inner product divided by the R discrepancy must exceed 10⁵. That is the essay’s title as an assertion, and it is the one that would catch a figure whose two panels had been swapped.

On the sliding loss curve the separation is asserted across the whole κ range rather than at its right-hand end, because at low precision both methods have saturated by κ = 10¹¹ and comparing them there measures nothing. The interesting κ moves left as bits are removed, so the assertion is that the gap opens somewhere — which is true at every position of the slider, and was arrived at after the first version of the check failed at 20 bits for exactly that reason.

Doing it twice #

Re-orthogonalisation deserves more than the sentence it got above, because it is the answer in the one context where Gram–Schmidt is unavoidable.

The idea is simply to run the orthogonalisation of a new column against the existing ones a second time. The first pass removes most of the components along previous directions; the second removes what the first pass’s rounding left behind. Since the residual after the first pass is already nearly orthogonal to everything, the second pass subtracts very little, and the coefficients it computes are small — so their relative error, however large, does no damage.

The classical result is Kahan and Parlett’s, usually quoted as “twice is enough”: for any vector whose component along the existing subspace is not overwhelming, two passes give orthogonality at the level of rounding, and a third would gain nothing. The condition for it to work is a mild one and can be tested cheaply — compare the norm before and after the first pass, and re-orthogonalise only if it dropped by more than a factor of about two.

That test is what “selective re-orthogonalisation” means, and it is the reason iterative eigensolvers can afford to keep their bases orthogonal without doubling the cost of every step.

Twice is enough, measured #

The theorem has been quoted here and this site does not leave a claim in that state. On the same eight-by-eight Hilbert matrix, κ = 1.53·10¹⁰, the four routes give ‖QᵀQ − I‖ of:

Two of those confirm what has already been said: classical alone has lost orthogonality completely, and modified sits at κu, which for this matrix is 1.7·10⁻⁶.

The third is the interesting one, and it is not quite the clean win the slogan promises. A second classical pass improves matters by twelve orders of magnitude — from useless to very good — but it lands at 10⁻¹², about a thousand times above where Householder sits. “Twice is enough” is therefore enough for its own purpose, which is to make the loss of orthogonality independent of κ rather than to make it vanish: 10⁻¹² does not grow as κ² and does not grow as κ, which is the whole content of the theorem. It is not a claim that two passes reproduce Householder.

The gap is also where the theorem’s hypothesis is doing work. Kahan and Parlett’s condition asks that a column’s component along the existing subspace not be overwhelming, and on the Hilbert matrix — whose columns are as nearly dependent as any standard example — it is strained about as far as it goes. That the result still holds to 10⁻¹² on the hardest available input is a better recommendation for it than a cleaner number on an easy one would have been.

All four routes reconstruct the matrix to about 10⁻¹⁷, as usual, and as usual that is the number which settles nothing.

What this looks like at other precisions #

The slider on the loss curve is worth exploring rather than reading past, because it settles a question people frequently get backwards.

At sixteen bits, classical Gram–Schmidt has lost orthogonality entirely by κ = 10³. At fifty-three bits it holds until about 10⁵ and is gone by 10⁹. Adding thirty-seven bits of precision moved the failure point by six orders of magnitude in κ — real, and a constant.

Modified Gram–Schmidt at sixteen bits fails around κ = 10⁴; at fifty-three, around κ = 10¹². Same shift, same reason.

What does not change at any precision is the ordering or the slopes. At every position of the slider, classical is worse than modified by a factor that grows with κ, and Householder is flat.

Reading one κ across the whole range is what makes the ordering concrete. Fix κ = 10¹¹ and add precision: the classical process loses 5.5, 4.5, 4.4, 3.5, 3.2 and 2.4 as the significand goes from twenty bits to fifty-three. A loss of one means no orthogonality at all, so every one of those numbers is a total failure, and thirty-three additional bits have moved it from thoroughly failed to thoroughly failed.

Householder over the same range goes 1.6·10⁻⁵, 2.3·10⁻⁷, 4.4·10⁻⁹, 6·10⁻¹¹, 8.3·10⁻¹³, 1.8·10⁻¹⁵ — ten orders of magnitude for the same thirty-three bits.

So precision and algorithm are not substitutes, and the measurement says so in the strongest available form. Every bit added is delivered in full to the orthogonal method and discarded entirely by the classical one, because Householder’s loss is bounded by a multiple of the unit roundoff while classical’s is bounded by a multiple of κu — and at κ = 10¹¹ that product is already past one at fifty-three bits. Buying a wider format to rescue classical Gram–Schmidt is therefore not merely poor value; it purchases nothing at all in the regime where the choice matters. The algorithm is free and the precision is not, which is the whole of the practical advice.

Why the R factor is not the whole story #

It is worth returning to the observation that both methods produce essentially the same R, because it explains why this bug is so durable.

A QR factorisation is used for several distinct purposes, and only some of them touch Q.

Solving least squares uses both: Rx = Qᵀb. Affected. Computing a determinant uses only the diagonal of R. Unaffected. Estimating a condition number uses only R. Unaffected. Producing an orthonormal basis — for a projection, a deflation, a subspace iteration — is entirely about Q. Affected completely.

So a codebase can use classical Gram–Schmidt for years, in a context where only R matters, and be perfectly correct. The failure arrives when somebody uses the same routine for a different purpose, and it arrives without any change to the routine.

That is the argument for measuring the property rather than trusting the code path: the code was not wrong, the use was new, and only a number attached to the output would have said so.

The generalisation worth taking away #

The pattern here is not confined to Gram–Schmidt, and it is worth stating in a form that survives outside numerical analysis.

A derivation that assumes an invariant holds exactly cannot see the cost of the invariant holding approximately. The proof that the two Gram–Schmidts agree uses orthogonality of the previous columns; in the computation, those columns are only nearly orthogonal, and the proof’s “zero” becomes a small number multiplied by something large.

That is a general failure mode of correctness arguments in the presence of approximation, and it recurs across the site. Elimination without a swap is correct under the assumption that the pivot is not tiny. The normal equations are correct under the assumption that AᵀA can be formed and inverted. Each derivation is valid; each assumption is quietly false; and in each case the resulting error is not small.

The remedy is the same in all of them: identify what the derivation assumed, measure it, and print the measurement beside the answer.

For the wider argument see the exact answer to a nearby problem, which gives the vocabulary, and the swap that is not optional, which is the same pattern in the first algorithm anyone learns.

Where the classical form still appears #

Given all of the above, it is fair to ask why the unstable version is written at all, and the answer is not only inattention.

Classical Gram–Schmidt has one genuine advantage: its inner loop computes all the projection coefficients against the same vector, so they are independent and can be computed as a single matrix–vector product. Modified Gram–Schmidt’s coefficients depend on each other in sequence, and cannot. On parallel hardware that difference is real — the classical form is one BLAS-2 call where the modified form is a chain of BLAS-1 calls, and the gap can be large.

Which is why the algorithm that actually gets used in high-performance Krylov codes is neither: it is classical Gram–Schmidt with re-orthogonalisation, which keeps the parallel structure and pays for it with a second pass. Two BLAS-2 calls beat a serial chain, and the orthogonality comes out at rounding level.

So the honest summary is not “classical is wrong”. It is that classical alone is wrong, classical twice is right and fast, and modified alone is right and slow. The version to avoid is the one that looks cheapest, which is the usual shape of this subject’s traps and the reason measured, not asserted is a thread rather than a slogan.

A third answer, which is neither of these #

Both algorithms on this page produce a Q whose columns span the same nested subspaces as A’s, in order. That property is what makes them differ from each other in the way this page measures — and it is also what makes both of them the wrong answer to a question that has no order in it.