Numerical Linear Algebra
where the algebra stops being the arithmetic
Orthogonality, measured

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.
16 min read 6 figures Identical algebra, different arithmeticMeasured, not assertedWhere the algebra stops being the arithmeticPrecision as the knob

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.

Classical and modified Gram–Schmidt: the same subtraction, in a different orderTwo panels of pseudocode differing in one argument, with the resulting pairwise dot products of the computed columns listed beneath each.for each previous column i, subtract the projection of column j onto qᵢclassicalr[i][j] = qᵢ · a[j] ↑ the ORIGINAL columnv = v − r[i][j] · qᵢthe three worst |qᵢ · qⱼ|:columns 7 and 8: 1columns 6 and 8: 0.13columns 6 and 7: 0.13modifiedr[i][j] = qᵢ · v ↑ what is LEFT of itv = v − r[i][j] · qᵢthe three worst |qᵢ · qⱼ|:columns 1 and 8: 4.4·10⁻⁷columns 2 and 8: 2.7·10⁻⁷columns 3 and 8: 2.4·10⁻⁸The two R factors agree to 1.2·10⁻⁶ relative. The two Q factors do not.the 8×8 Hilbert matrixone word, eight orders
Fig. 1 The two versions side by side, with the pairwise inner products of the columns they produce. On the eight-by-eight Hilbert matrix, the classical version leaves pairs of “orthonormal” columns whose inner product is of order one; the modified version leaves them at 10⁻⁸. The two R factors agree to within κu.

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.

Loss of orthogonality against condition number, in binary64A log–log plot of the norm of Q-transpose-Q minus the identity against condition number. Classical Gram–Schmidt rises steeply, modified Gram–Schmidt rises gently, and Householder is flat.110²10⁴10⁶10⁸10¹⁰10¹²10⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹condition number κ(A)‖QᵀQ − I‖classicalmodifiedHouseholderκ²uκu8×8, eight seeds per κ, binary64all three reconstruct A
Fig. 2 The two rates, measured over eleven decades of condition number with eight seeds at each point, and the two predictions drawn as dashed lines. Classical follows κ²u; modified follows κu; Householder does not follow anything, because it does not lose orthogonality at all.

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.

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.

QᵀQ from classical Gram–Schmidt and from Householder on the 8×8 Hilbert matrixTwo eight-by-eight tables of QᵀQ. The upper one has ones on the diagonal and entries as large as one off it; the lower one is the identity to three decimal places everywhere.A = H8 · κ = 1.5·10¹⁰ · both factorisations reconstruct A to 6·10⁻¹⁷the diagonal is 1 in both — every column is a unit vector either way1.000000000001.000000000001.000000000001.000000000001.00000.002-0.002000001.0000.125-0.13300000.0020.1251.000-1.0000000-0.002-0.133-1.0001.000classical Gram–Schmidt1.000000000001.000000000001.000000000001.000000000001.000000000001.000000000001.000000000001.000Householderclassical ‖QᵀQ − I‖1.4Householder ‖QᵀQ − I‖1.4·10⁻¹⁵largest off-diagonal 1 against 3.1·10⁻¹⁶length is not angle
Fig. 3 The two Q’s, judged. QᵀQ printed for the classical factorisation and for a Householder one on the same matrix. Both have an exact 1 down the diagonal — the columns are unit vectors either way — and only one of them is the identity anywhere else.

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.

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. So the question of precision against algorithm has the same answer everywhere on the range, and the algorithm is free.

How far a perturbation of size ε moves an eigenvalue, 8×8A log–log plot of eigenvalue movement against perturbation size. The symmetric case lies on a line of slope one; the non-symmetric case lies on a line of slope one eighth, and at a perturbation of ten to the minus sixteen it has already moved by a hundredth.10⁻¹⁶10⁻¹³10⁻¹⁰10⁻⁷10⁻⁴10⁻¹10⁻¹⁶10⁻¹³10⁻¹⁰10⁻⁷10⁻⁴10⁻¹size of the perturbation ‖δA‖how far the eigenvalues moveJordan block, ε^(1/8)symmetric, ≤ ‖δA‖rounding error alone moves it to 10⁻²six seeds per symmetric point; Jordan is closed formsymmetry beats precision
Fig. 4 Where the choice has consequences beyond the factorisation. An eigenvalue computation on a basis that is not orthogonal is a computation on a perturbed problem, and for a non-symmetric matrix a perturbation of 10⁻⁸ can move the spectrum by 10⁻¹. The loss of orthogonality is the perturbation.

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 least-squares solution as a projection, with the right angle measuredThe column space drawn edge-on as a plane, the data vector above it, and the perpendicular dropped to the plane, with the residual marked at a right angle to it.everything Ax can reachb = (1.1, 0.4, 1.5)Ax, the closest reachable pointr = b − Ax‖Aᵀr‖ / (‖A‖‖r‖)1.7·10⁻¹⁶‖b‖² − ‖Ax‖² − ‖r‖²1.3·10⁻¹⁵‖r‖1.3200 random nearby points of the plane were tried; none is closer.a 3×2 system, Householder QRperpendicularity is checked
Fig. 5 The purpose that needs Q. The least-squares solution is defined by a perpendicularity condition, and the QR route enforces it through Q. The badge here reports how well the condition is met — 10⁻¹⁶, with a Householder factorisation. With a classical Gram–Schmidt Q on an ill-conditioned matrix, the same figure would report a number many orders of magnitude larger.

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.

Growth factor under partial pivoting: the bound, the worst case, and realityGrowth factor against matrix size on a logarithmic vertical axis. The two-to-the-n bound rises as a straight line; Wilkinson's matrix sits exactly on it; random matrices stay near one.0816243240110²10⁴10⁶10⁸10¹⁰10¹²10¹⁴matrix size ngrowth factor max|u| / max|a|the 2ⁿ⁻¹ boundworst of 30 randommedian randomWilkinson's matrix sits on the bound30 Gaussian matrices per sizeat n = 40: bound 5.5·10¹¹, worst 4.8
Fig. 6 The same question in the elimination field. There, a stability property rests on a measured quantity staying near three when its bound permits 5.5·10¹¹; here, it rests on which of two vectors a coefficient is computed against. Both are decided empirically, both matter, and neither is visible in a derivation.

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.

Where this goes next

A = H8 · κ = 1.5·10¹⁰ · both factorisations reconstruct A to 6·10⁻¹⁷the diagonal is 1 in both — every column is a unit vector either way1.000000000001.000000000001.000000000001.000000000001.00000.002-0.002000001.0000.125-0.13300000.0020.1251.000-1.0000000-0.002-0.133-1.0001.000classical Gram–Schmidt1.000000000001.000000000001.000000000001.000000000001.000000000001.000000000001.000000000001.000Householderclassical ‖QᵀQ − I‖1.4Householder ‖QᵀQ − I‖1.4·10⁻¹⁵largest off-diagonal 1 against 3.1·10⁻¹⁶length is not angle 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. the mirrorx, length 4.000Hx = (-4.000, 0)v = x − αe₁safe sign: α = −‖x‖, so v is formed from a sum and nothing cancelsunsafe sign: α = +‖x‖ gives ‖v‖ only 35.1% of ‖x‖ + ‖x‖ — the digits go‖HᵀH − I‖5·10⁻¹⁶‖Hx‖ − ‖x‖8.9·10⁻¹⁶second component2.2·10⁻¹⁶built from a unit vectororthogonality is structural 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. [ ε 1 ; 1 1 ] x = [ 1 ; 2 ], exact answer (1.000000, 1.000000)with partial pivoting1101U after elimination1.0000001.000000computed xbackward error 0forward error 0without10⁻¹⁷10-1·10¹⁷U after elimination0.0000001.000000computed xbackward error 0.25forward error 0.71no error is raisedgrowth 10¹⁷ 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.