The right-hand side as one more column

Two Gram–Schmidts ends its account of what a lost orthogonality costs with a sentence about least squares: the solution comes from Rx = Qᵀb, that derivation needs QQᵀ to be a projection, and with an unstable Q “the computed x can be wrong by a factor of κ more than it needs to be.” The sentence is stated there and not measured, and it is the natural reading of everything that essay does measure. Modified Gram–Schmidt loses orthogonality like κu; multiply b by a Q that is κu away from orthogonal, solve a triangular system with condition number κ, and the answer should be κ²u away from right.

It is, by one route. By another route the same routine, handed the same matrix and the same right-hand side, returns an answer more accurate than Householder’s — while its Q is exactly as far from orthogonal as before. The difference is where b enters the computation, and it is one line.

Two ways to reach Qᵀb, and one of them never multiplies by Qᵀ #

A least-squares solve by QR needs two things: the triangular factor R, and the first n entries of Qᵀb. Every Gram–Schmidt produces R directly. The question is how the second object is formed, and there are two answers that are identical in exact arithmetic.

    route 1 — through Qᵀb                 route 2 — b as one more column
    [Q, R] = mgs(A)                       [Q̃, R̃] = mgs([A  b])
    c = Qᵀ b                              R = R̃(1:n, 1:n)
    x = R \ c                             c = R̃(1:n, n+1)
                                          x = R \ c

The first route is what a reader writes from the formula. It computes each cᵢ = qᵢᵀb as a dot product of the finished basis vector with the original right-hand side.

The second route asks the orthogonalisation to treat b as an (n+1)-th column of the matrix, runs the loop one step further than it needs to, and throws away the last basis vector. The first n entries of R̃’s last column are exactly the coefficients the loop subtracted from b on its way to making it orthogonal to everything — which is Qᵀb, by the definition of the loop. Nothing is ever multiplied by Qᵀ. And in the modified loop each of those coefficients is computed against what is left of b after the previous subtractions, not against b itself.

That last clause is the whole of the difference between the two Gram–Schmidts, arriving in a new place. Two Gram–Schmidts locates the instability of the classical algorithm in one argument: projecting the original column rather than the running one. Route 1 is the classical projection applied to b by an algorithm that was modified for every other column. Route 2 is modified all the way down.

The problem is built so the answer has nowhere to hide #

A forward error needs a true answer, and a relative forward error needs a denominator that stays still. The road that squares the problem found the trap in that second requirement: on a random right-hand side the least-squares solution is dominated by the smallest singular direction, ‖x‖ grows like κ, and a relative error measured against it is quietly divided by κ. A method losing κ² digits can look like one losing κ.

So the problems here are consistent by construction. A is a 40×8 matrix built as UΣVᵀ with singular values spaced logarithmically from 1 down to 1/κ, x★ is the fixed vector (1, 1.1, …, 1.7), and b = Ax★. The least-squares solution’s norm is 3.87 at every condition number from 10² to 10¹⁴, and that is checked before the figure is drawn at all.

The reference is not x★. Rounding happened when A and b were formed, so the problem a program holds is slightly different from the one written down, and its solution is slightly different too. The reference is the exact solution of the problem as stored: AᵀA and Aᵀb formed in rational arithmetic from the doubles themselves, and solved exactly. Forming AᵀA squares the condition number, but squaring costs nothing when no digit is rounded, and every error below is the error of a route rather than of the reference. It is the same move an answer that is known makes with the Hilbert matrix’s rational inverse.

Three seeds at every κ, and the median reported.

κ² for one route and κ for the other #

Read the second column first, because it is the prediction. Modified Gram–Schmidt through Qᵀb has a fitted slope of 2.01 in κ over the first seven decades: κ²u, as the naive argument said, 0.13 at κ = 10⁸ and more than a thousand at 10¹⁰. An answer with a relative error of 1,400 is not an inaccurate answer; it is a number of the wrong size pointing in an arbitrary direction.

Now the third column. The same algorithm, the same loop, the same matrix — and a fitted slope of 0.85. At κ = 10⁸ it is 2.7·10⁻¹⁰, which is 470 million times more accurate than route 1 on identical data. It is also fifteen times more accurate than Householder, whose column is the standard every library meets. It is ahead of Householder at all seven stops: by a factor of seven to fifteen from κ = 10⁴ to 10¹⁰, and by a few per cent at either end of the range.

And the last column is the one that makes the result strange. At κ = 10⁸ the Q that route 2’s orthogonalisation built — the same Q as route 1’s, plus one column — is 4.3·10⁻⁹ from orthogonal. The answer computed alongside it is 2.7·10⁻¹⁰ from the truth. The solution is sixteen times more accurate than its own basis is orthogonal. Any account of least squares that runs through “the answer can be no better than the Q” has to explain that number, and the account in the next two sections does.

It is worth being exact about what this figure establishes and what it does not. It establishes that ‖QᵀQ − I‖ for modified Gram–Schmidt grows like κu on square test matrices with prescribed conditioning, which is the quantity every discussion of the method starts from. It does not say anything about what that Q was for. Orthogonal is a number turned “orthogonal” from an adjective into a measurement; the measurement is real, and it is a measurement of the basis, not of any answer computed with it.

Why classical Gram–Schmidt gets nothing from the column #

The fifth column of the table has one entry per row, and that is not a typesetting economy. The two classical routes produce identical answers — not close, identical: the difference between the two computed solutions is zero at all seven condition numbers.

That is what the loop says should happen, and it is the cleanest possible control for the experiment. In the classical loop every coefficient of a column is computed against the original column. So when b is appended as an extra column, the coefficients the loop computes for it are q₁ᵀb, q₂ᵀb, …, qₙᵀb — the same dot products, of the same vectors, summed in the same order, as forming Qᵀb explicitly. The augmented route is route 1 written differently, and floating point agrees with the algebra to the last bit.

So the benefit of the column is not a benefit of appending anything. It is a benefit of computing b’s coefficients the way modified Gram–Schmidt computes every other coefficient: against the running remainder. Classical Gram–Schmidt has no running remainder to offer, and gets κ²u — 8.0 at κ = 10⁸, which is to say no digits — however b is handed to it.

Modified Gram–Schmidt is a Householder factorisation of a bigger matrix #

The explanation is a theorem of Björck and Paige, and it can be checked rather than quoted, which is the habit of this collection.

Take A and put n rows of zeros above it: an (m+n)×n matrix [0; A]. Run Householder QR on that. The result is claimed to be the modified Gram–Schmidt factorisation of A, in floating point and not only in the algebra. Specifically: the n×n R that Householder produces for the padded matrix is modified Gram–Schmidt’s R; the bottom m×n block of the Householder Q is modified Gram–Schmidt’s Q, to within the orthogonality that Q lost; and the top n×n block — which would be zero if the padding contributed nothing — is where that lost orthogonality went.

Measured over the same seven condition numbers, with the Householder diagonal signs fixed so that the comparison is not a comparison of sign conventions:

The R factors agree to 4·10⁻¹⁶ at every condition number from 10² to 10¹⁴. Not to κu, which is how closely classical and modified Gram–Schmidt’s R factors agree with each other: to the rounding of the arithmetic, independent of κ. The two algorithms are computing the same numbers.

And three quantities have the same size at every stop: the top block’s norm, the orthogonality modified Gram–Schmidt lost, and the distance between modified Gram–Schmidt’s Q and the padded Q’s bottom block — 3.2·10⁻⁹, 4.3·10⁻⁹ and 3.7·10⁻⁹ at κ = 10⁸, and never more than a factor of 2.2 apart anywhere on the axis. The padded Q is a product of reflections and is orthogonal by construction — a reflection cannot stop being one — to 1.6·10⁻¹⁵. Modified Gram–Schmidt computes the part of that orthogonal matrix that sits under A and never forms the part that sits under the zeros, and the defect everyone measures is the size of the part it did not form. The exact algebra tying the blocks together is Björck and Paige’s and is not re-derived here; what is measured is that the sizes agree at all seven condition numbers.

That turns the least-squares result into a one-line argument. Least squares by Householder on the padded augmented matrix [0 0; A b] is backward stable; that is the property Householder has for any matrix, padded or not. Applied to the 40×9 matrix [A b], the equivalence says modified Gram–Schmidt computes the same R as that factorisation, last column included, to rounding. So route 2’s c and R are a backward-stable factorisation’s c and R, and the answer is a backward-stable answer. The augmented route is Householder in disguise, and it inherits Householder’s law.

Route 1 has no such cover. It keeps modified Gram–Schmidt’s R, which is fine, and then forms Qᵀb using only the bottom block of an orthogonal matrix whose top block it never computed. The missing block is κu in size, it enters c multiplied by ‖b‖, and the triangular solve with R multiplies it by κ. That is the κ²u of the second column, derived rather than fitted.

Why route 2 is slightly more accurate than a real Householder solve is not settled by this argument, and it is not claimed anywhere. The measured advantage is a factor of seven to fifteen in the middle of the range and a few per cent at its ends; a plausible reading is that the Householder route forms Qᵀb through an m×m sequence of reflections applied to b while route 2 never touches the extra m − n directions, but that is a reading and not a measurement.

A residual puts the stable routes on the κ² law too #

A consistent system is the clean experiment and it is not the common case. A real fit has a residual, and least squares with a residual has a sensitivity of its own that no algorithm can beat: the perturbation bound carries a term proportional to κ²‖r‖/(‖A‖‖x‖). Once that term dominates, a backward-stable method is entitled to a κ² error, because the problem itself is that sensitive.

So the same four routes were measured with a component orthogonal to the column space added to b, of relative size 10⁻⁶ and then 10⁻².

At a residual of 10⁻⁶ the table is unchanged up to κ = 10⁸ — 7.1·10⁻⁹ for Householder and 4.8·10⁻⁹ for the augmented route, against 4.0·10⁻⁹ and 2.7·10⁻¹⁰ without a residual. The first route’s error at κ = 10⁸ is 0.128 in both experiments, to three figures. It has nothing to lose from the residual, because it had already reached κ²u on the consistent problem.

At a residual of 10⁻² the stable routes meet the κ² term much earlier.

Two things change, and both are the problem’s rather than the method’s. At κ = 10⁸ the stable routes now sit at 2.9·10⁻⁵ and 4.7·10⁻⁵ — the augmented route is 60 per cent worse than Householder at this stop, which is the first place in the whole sweep where it is not at least level — and the gap between them and route 1 has closed from a factor of 470 million to a factor of 2,700. And past κ = 10⁸ the least-squares solution itself moves: its norm is 6.4 at κ = 10¹⁰ where it was 3.87, and 90,000 by 10¹². The figure stops where the denominator starts moving, for the reason the road that squares the problem gave for insisting on a fixed ‖x‖.

So the honest statement has a condition in it. On a problem whose residual is small against 1/κ², modified Gram–Schmidt on [A b] is as accurate as Householder and route 1 is κ times worse than both. On a problem whose residual is large, all routes share the problem’s own κ² sensitivity, and route 1’s extra κ matters less and less. It never stops mattering in these measurements — 2,700 at the right edge of the worst figure is still more than three orders — but the dramatic ratio belongs to nearly consistent problems, which are the ones where least squares is most often trusted.

What the augmented route does not repair #

The column fixes x. It does not fix Q, and it is worth listing what still depends on Q directly, because a code that adopts the augmented solve has not thereby acquired an orthonormal basis.

A projection formed with Q. Computing the residual as (I − QQᵀ)b, or projecting another vector onto the column space, uses the κu-defective basis and inherits its defect. The residual computed as b − Ax with the augmented x is fine; the residual computed from the basis is not. The two agree only on a well-conditioned A, which is the case nobody needed to check.

Leverage and anything built from rows of Q. The diagonal of QQᵀ is the leverage of each observation, decided before the data and computed from the rows of an orthonormal basis. A basis off orthogonal by κu reports leverages whose sum is no longer exactly the number of columns.

A Krylov method. An orthogonalisation nobody calls one is the case where the basis is the product, and there the loss of orthogonality produces duplicate Ritz values and stalled convergence. No column of right-hand side rescues that, because nothing is being solved against a fixed b.

Several right-hand sides. The augmented trick costs one extra column per right-hand side, and for k of them the natural economy — factorise once, apply Qᵀ to each b — is route 1, with route 1’s error. Appending all k columns at once works, at the cost of orthogonalising k extra columns against the n basis vectors and against each other.

The practical rule therefore has two halves. If only x is wanted, modified Gram–Schmidt on [A b] is a backward-stable least-squares solver and needs no reorthogonalisation. If Q is wanted, the orthogonality has to be paid for — by a second pass, as doing it twice does for Cholesky QR, or by Householder.

What the old sentence got right #

The sentence this essay started from is not wrong. It is a sentence about route 1, stated as a sentence about Gram–Schmidt, and the measurement puts a number on each half.

For route 1 it understates. “Wrong by a factor of κ more than it needs to be” is exactly right as a ratio against Householder: 0.13 against 4.0·10⁻⁹ at κ = 10⁸ is a factor of 3.2·10⁷, near κ. But the phrase “can be” suggests a worst case, and the second column of the table is not a worst case — it is the median of three random problems, and it is on its κ²u law at every one of them.

For the algorithm it is wrong in the way the exact answer to a nearby problem teaches a reader to look for: the error bound quoted belongs to one composition of the parts, and the parts compose another way with a different bound. The loss of orthogonality is a real property of the Q; the least-squares error is a property of how b met the loop. Two quantities that the derivation treats as one — “the Q” and “what the Q computes” — separate by a factor of sixteen in one direction and by 470 million in the other.

And it is a third instance of a pattern that keeps recurring. In the road that squares the problem the seminormal equations inherit κ² from composing two triangular solves, and one correction step removes it. In a small residual is not a small error the quantity that looks like the answer’s quality is not. Here the quantity that looks like the method’s quality — ‖QᵀQ − I‖, the number two Gram–Schmidts plots — is not the quantity the most common use of the method pays for.

Three numbers checked at every stop #

The figure’s own checks are statements about the claim, made at every stop rather than at one. At every condition number drawn, the augmented route must be within a factor of ten of Householder, and the two classical routes must differ by less than 10⁻¹⁴. Somewhere past κ = 10⁶ the Qᵀb route must be at least a thousand times worse than the augmented one. The residual placements reuse the same checks, which is why the one stop where the augmented route is sixty per cent behind Householder is allowed and a stop where it fell a factor of ten behind would not be.

Behind the figures, the measurement requires the two slopes — between 1.7 and 2.3 for route 1 and between 0.7 and 1.3 for route 2 — over the decades where neither has saturated, and requires that the reference solution’s norm is flat to one per cent, so that no figure’s comparison can be quietly divided by a moving denominator. The padded comparison requires that the R factors agree to 5·10⁻¹⁵, that the padded Q is orthogonal, and that its top block and modified Gram–Schmidt’s lost orthogonality are within a factor of thirty of each other at every κ. The factor of thirty is loose on purpose: the theorem ties the two to first order, and the measurement finds them within 2.2.

The refusal is the claim in the title of the old sentence. A claim that the augmented error is at least a tenth of κ times the orthogonality defect is fed the κ = 10⁸ problem and must fail.

Where Gram–Schmidt goes next #

Everything above orthogonalises one column at a time, and almost no large computation does. A communication-avoiding or cache-aware code works in blocks of columns: it removes a block’s components along the blocks already done, and orthogonalises what is left with a smaller QR inside the block — Householder, or Cholesky QR, or a tree of small factorisations. That is two orthogonalisations stacked, and each has its own stability.

The figure says the stable inner factorisation does not make the block method stable, and a stable block is not a stable basis takes that apart by placing the ill-conditioning between the blocks, inside them, and both.

Two further questions sit behind it and are not measured here. The first is whether the augmented column survives blocking: a block method solving least squares would append b to the last block, and whether the Björck–Paige equivalence has a block form that covers it is open. The second is Gram–Schmidt in an inner product other than the Euclidean one, xᵀMy, which generalised eigenvalue problems require; the published bounds put a factor of κ(M) into the loss of orthogonality, and whether “twice is enough” survives an ill-conditioned M is a measurement not yet made.