A stable block is not a stable basis

Gram–Schmidt as two Gram–Schmidts measures it takes one column at a time. The version that runs on large machines takes several. A block Gram–Schmidt splits the columns into groups of p, and for each group does two things: it removes the group’s components along the blocks already finished, and it orthogonalises what is left with a small QR inside the group. The first step is where the method touches everything done so far, and it is one matrix product rather than a chain of dot products. The second step is where the block’s own columns are made perpendicular to each other, and it can be any QR at all.

That structure is why the method exists. The message and the word counts what a factorisation sends between processors, and a block method sends far fewer messages than one that visits the columns singly; the same arithmetic at a different price is the same argument about memory. It is also why the method’s stability is a question with two parts. There are two orthogonalisations stacked on top of each other, and a reader who knows that Householder cannot lose orthogonality might reasonably expect that a block method with Householder inside every block cannot lose it either.

Four variants and three matrices #

The variants differ at the two levels independently.

Between blocks, the removal can be classical — project the block against the whole finished basis at once, W ← W − Q(QᵀW) — or modified — project against one finished block at a time, updating W after each. It can be done once, or twice, which is the block form of re-orthogonalisation. Inside a block, the QR can be Householder, which is orthogonal by construction, or Cholesky QR — factor WᵀW = RᵀR and set Q = WR⁻¹ — which costs one Gram matrix and one triangular solve and is the cheapest QR there is to communicate.

Four combinations are drawn: classical once with Householder inside, modified once with Householder inside, classical twice with Householder inside, and classical twice with Cholesky QR inside. The last is close to what doing it twice recommends for a single tall block, placed inside a block method.

And the matrix is built three ways, because κ(A) is one number and the ill-conditioning of a blocked matrix lives in two places.

Across: take a generic 64×16 matrix of prescribed conditioning and replace every block of four columns by an orthonormal basis for itself. Each block now has condition number exactly one. The column space of each block is unchanged, so the angles between blocks are what they were, and κ(A) — 4.3·10⁷ when 10⁸ was asked for — lives entirely there.

Inside: build A = UΣVᵀ with V block-diagonal, so that different blocks are exactly orthogonal to each other in the algebra, and grade Σ from 1 to 1/κ within each block. Now there is nothing to remove between blocks, and κ(A) = κ lives entirely inside them.

Mixed: the generic matrix, ill-conditioning wherever a random V puts it.

Three seeds for every matrix, the median reported, and a breakdown on any seed reported as a breakdown — a method that stops on one matrix in three is not a method whose median means anything.

Across the blocks: Householder inside changes nothing #

The first figure is the case where the within-block QR has nothing to do. Every block arrives orthonormal; Householder and Cholesky QR alike are handed a matrix they return almost unchanged. Whatever happens to the basis happens in the between-block step.

Classical once, with Householder inside, is classical Gram–Schmidt. Its loss grows by about four orders for every two decades of κ, the κ² law, and it has lost orthogonality completely — ‖QᵀQ − I‖ above one — by κ = 4·10⁹. That every block was orthogonalised by the method that cannot lose orthogonality is invisible in the result, because each block was already orthonormal and what went wrong was the removal of the blocks’ components along each other.

Modified once is modified Gram–Schmidt: the κu law, 5.7·10⁻⁹ at κ = 4.3·10⁷, a factor of 700,000 better than the classical order on the same matrix. The mechanism is the one-column mechanism at block scale. The classical removal computes QᵀW against the original block, so every error already in Q enters the coefficients multiplied by the full size of W; the modified removal computes each block’s coefficients against what is left of W after the blocks before it, which for a nearly dependent block is small.

For calibration, the single-column algorithms on the same matrix: classical Gram–Schmidt returns 0.15 at κ = 10⁸ requested and modified returns 5.9·10⁻⁹. Block modified is on modified’s line almost exactly. Block classical, at 4.2·10⁻³, is thirty-six times better than single-column classical and on the same slope — fewer projections, the same law.

Both twice variants are flat. The second removal subtracts almost nothing from a block that is already nearly orthogonal to the basis, so its coefficients are small and their errors do no damage — the argument two Gram–Schmidts gives for “twice is enough”, applied to blocks. And Cholesky QR inside causes no trouble here, because the blocks it factorises are perfectly conditioned.

Inside the blocks: a second pass cannot undo a breakdown #

Move the ill-conditioning into the blocks and the picture changes at both levels.

The classical and modified orders no longer differ. They cannot: in the algebra there is nothing between the blocks to remove, and the order in which nothing is removed does not matter. The two columns agree to within a factor of 2.4 at every stop.

And yet both lose orthogonality like κu, 9.8·10⁻³ at κ = 10¹⁴, although nothing was there to remove and every block was orthogonalised by Householder. The loss is manufactured by the two steps together. The between-block projection W − Q(QᵀW) removes nothing but leaves rounding behind, of size about u‖W‖, pointing along the finished basis. The Householder factorisation inside the block is then asked for W = QR, and Q = WR⁻¹ with R the block’s own triangular factor: the rounding is carried into the new basis vectors multiplied by ‖R⁻¹‖, which is the block’s own condition number. So a block with κ = 10¹⁴ inside it turns a residue of 10⁻¹⁶ into a component of 10⁻² along blocks it was supposed to be perpendicular to. The inner QR is perfectly orthogonal within the block, and it is the thing that amplifies the outer step’s rounding.

The second pass repairs exactly that, which is why classical twice with Householder is flat at 2.7·10⁻¹⁵ to 3.4·10⁻¹⁵ here too. The second removal acts on that block’s Q, whose columns are orthonormal and perfectly conditioned, and removes the κu components without anything to amplify them.

Cholesky QR inside cannot be repaired by a second pass, because it does not return a bad basis for the second pass to improve; it returns no basis. Cholesky QR factorises the block’s Gram matrix WᵀW, whose condition number is κ(W)², and once that reaches about 1/u the Gram matrix is no longer reliably positive definite in floating point and the factorisation can stop on a non-positive pivot. A 4×4 block graded to κ = 10⁸ has a Gram matrix at κ = 10¹⁶, on the edge: one seed of three stops there, in its fourth block, and all three stop by 10¹⁰. The breakdown is on the first pass, and “twice” is a promise about what the second pass does to a first pass that finished.

That is the result doing it twice found for Cholesky QR on its own, now located. Cholesky QR twice is excellent when the thing it factorises is not too ill-conditioned; a block method decides what it factorises, and when the ill-conditioning is inside the blocks, that is exactly what Cholesky QR is handed.

The generic matrix behaves like the across case #

Neither construction is what a real matrix looks like, so the third is there to say which one a generic matrix resembles. At blocks of four it resembles across: classical once at 5.0·10⁻³ against across’s 4.2·10⁻³, modified once at κu, and Cholesky QR inside never stopping in any of the three seeds up to κ = 10¹⁴. A random orthogonal V spreads the small singular values over all sixteen columns, so each block of four sees a slice of the spectrum and is far better conditioned than the matrix, while the angles between blocks carry the rest.

That is a statement about this construction and not a theorem. A matrix whose columns arrive in an order that groups nearly dependent columns together — a polynomial basis in increasing degree, a Krylov sequence — is closer to inside than a random V is, and the placement a real code meets is decided by the order its columns were written in. The measurement says what the two extremes do; it does not say where an arbitrary application sits between them.

How many blocks the κ² law needs #

The block size is not a detail, and the across construction shows why in three figures.

With blocks of two, classical once is at 2.9·10⁻² where blocks of four gave 4.2·10⁻³, and modified once is at 7.7·10⁻⁹ against 5.7·10⁻⁹. The classical figure moves towards single-column classical Gram–Schmidt as the blocks shrink, as it must: blocks of one are classical Gram–Schmidt.

With blocks of eight there are two blocks, and the classical and modified removals are the same computation: the second block is projected against the first, once, and there is no second finished block for the orders to disagree about. The two curves coincide to every printed digit, and they sit on the κu law, not κ². The κ² law of classical Gram–Schmidt is a law about repeated projection — it needs a second finished block whose errors enter the coefficients multiplied by the full size of the block being projected — and a factorisation with two blocks never reaches it.

That is the same observation the single-column algorithm makes when n = 2, and it is worth saying at block scale because a block size is usually chosen for the machine: a block size is a property of the machine finds the best one from the size of fast memory. On this measurement the same choice also decides which law classical removal follows, and a code that doubles its block size to fit a larger cache has changed its stability as a side effect.

The inside construction at blocks of eight moves the Cholesky breakdown the other way.

A block of eight graded from 1 to 10⁻⁸ has the same condition number as a block of four graded over the same range, and its Gram matrix the same 10¹⁶ — yet on all three seeds its Cholesky factorisation completes at κ = 10⁸ and stops only at 10¹⁰. Why is not settled here. The pivots of a Cholesky factorisation depend on how the small singular values are arranged against the columns and not only on the largest and smallest of them, and eight graded values are arranged differently from four. What the measurement does settle is the reading of the rule: κ(W)² reaching 1/u is where a breakdown becomes possible, not where it must happen — which is also why one seed of three broke at 10⁸ with blocks of four and the other two did not.

Which level each failure belongs to #

Put the three placements beside each other and each variant has one condition number it cannot see and one it cannot survive.

Classical once, Householder inside: survives ill-conditioning inside blocks at κu, loses κ² to ill-conditioning between blocks once there are at least three of them.

Modified once, Householder inside: κu from either placement. The between-block order protects against the across case; nothing at that level protects against the inside case, where the loss is the inner factorisation amplifying the outer step’s rounding.

Classical twice, Cholesky QR inside: indifferent to ill-conditioning between blocks, at 1·10⁻¹⁵, and liable to stop on ill-conditioning inside them once κ(W)² reaches about 1/u.

Classical twice, Householder inside: 2.5·10⁻¹⁵ to 4.2·10⁻¹⁵ at every stop of all nine combinations of placement and block size measured, seven condition numbers each.

The last is the only variant that needs no knowledge of where the ill-conditioning sits, and it is also the most expensive: two between-block removals per block and a Householder factorisation inside each, which the message and the word prices as the variant with the most communication. The cheap variants are cheap because each trusts one level. That trust is a claim about the matrix, and a single κ(A) cannot confirm it — 10⁸ was the requested κ in both the across and inside tables, and the variant flat on one broke down on the other.

What the measurements here require #

The checks are written about the claim at each placement rather than about one matrix. Classical twice with Householder inside must be below 10⁻¹³ at every stop, whatever the placement and block size. Inside, the classical and modified orders must agree within a factor of ten at every κ. With two blocks they must agree to six digits, because they are the same computation. With three or more blocks across, classical once must be at least a thousand times worse than modified once at κ = 10⁸. And Cholesky QR inside must not stop on the across or mixed placements, which is the claim that its breakdown is a property of the block rather than of the matrix.

The measurement behind them requires the κ² slope for classical once on the across and mixed placements, the breakdown of Cholesky QR between 10⁸ and 10¹⁰ on the inside placement and its survival below the first breakdown, and the κu loss of both once-variants inside at 10¹⁴. Its refusals feed the claim this essay was written against — Householder in every block, therefore stable — a matrix that breaks it, and require a second pass offered as the repair for a Cholesky breakdown to fail as well.

Where this leaves the practical choice #

A block Gram–Schmidt in a Krylov solver — an orthogonalisation nobody calls one is the single-vector version, and how wide the block should be the question of what a block can reach — faces the mixed case with a twist: its columns are a sequence, and consecutive Krylov vectors are nearly parallel by design, which is the inside placement’s structure arriving through the ordering. That is the setting in which the inner factorisation is the one at risk, and the measurement here says why a Cholesky QR inside the block would fail there specifically rather than everywhere.

A tall-skinny factorisation for a least-squares solve faces something closer to the across case, and there the cheap inner QR is safe. The rule the three tables support is therefore not “always pay for Householder twice”. It is: decide which level the ill-conditioning is at, and pay for stability at that level — a second between-block pass when the blocks are nearly dependent on each other, a robust inner factorisation when the columns within a block are — and pay for both when the question cannot be answered.

Where block Gram–Schmidt goes next #

Three measurements are natural from here and none is made.

Low-synchronisation variants. Block methods have been proposed that replace the second between-block pass with a cheaper correction built from the Pythagorean identity for the block’s norm, aiming at one all-reduce per block. Whether that correction keeps the across case at rounding, and what it does to the inside case, is the first question a practitioner with a communication budget would ask.

Randomised Gram–Schmidt. Orthogonalising against a sketch of the basis rather than the basis itself changes the inner products the loop computes, and with them which errors enter the coefficients. The sketching essays measure what a sketch preserves; whether it preserves the between-block stability is open.

The augmented column in blocks. The right-hand side as one more column found that modified Gram–Schmidt on [A b] is a backward-stable least-squares solver. A block solver would append b to its last block, and whether the equivalence behind that result has a block form — and which of the four variants here inherits it — is the question that connects this essay back to the one before it.