What the appended block inherits

The right-hand side as one more column found that modified Gram–Schmidt is a better least-squares solver than its own Q deserves. At κ = 10⁸ its Q is 8.3·10⁻¹⁰ from orthogonal, and multiplying b by that $Q^{\mathsf T}$ gives an answer wrong by 10⁻². Hand the same routine b as an extra column and it is right to 2.7·10⁻¹⁰, better than Householder. The reason is Björck and Paige’s: modified Gram–Schmidt is exactly Householder QR on the matrix with n rows of zeros stacked on top, and its last column of R is $Q^{\mathsf T}b$ computed the way that factorisation computes it, each coefficient against what was left of b after the previous ones. Classical Gram–Schmidt projects the original b every time, so its last column of R is the same dot products as $Q^{\mathsf T}b$ and it gains nothing, to the last bit.

A stable block is not a stable basis then moved to blocks, measured four block variants’ loss of orthogonality with the ill-conditioning placed between the blocks, inside them or everywhere, and left the connecting question open: “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.”

It has an answer, and the answer comes with a second finding the orthogonality measurements could not have produced: the variant whose Q is the most orthogonal of the four is, on one placement, the worst least-squares solver.

Four variants, two routes, one exact answer #

The matrices are the preceding essay’s: 64 rows, 16 columns, blocks of four, κ from 10² to 10¹², with the ill-conditioning either spread through a generic matrix, placed only in the angles between blocks, or placed only inside each block. The four variants are one pass of block classical Gram–Schmidt with Householder inside each block, one pass of block modified Gram–Schmidt with Householder inside, two passes of block classical with Householder inside, and two passes of block classical with Cholesky QR inside.

Each variant now also returns its R, and each is used two ways. The first multiplies b by $Q^{\mathsf T}$ and back-substitutes. The second does what a block code does when b arrives with the matrix: appends it as a block of width one and runs it through the same between-block projections as every other block — against all the finished blocks at once for a classical outer loop, against them one at a time for a modified one, and twice for the two-pass variants. Its coefficients are then the last column of R. The problems are consistent, b = Ax for a fixed x of moderate size, so that the relative error does not shrink as κ grows, and every answer is compared with the least-squares solution computed exactly in rationals from the stored doubles.

The R of a two-pass variant needs a word, because it is where the second finding lives. Each pass removes the finished blocks’ components and factorises what is left, so a block’s final columns are its original columns times the inverse of the product of the two passes’ triangular factors, and its coefficients against the finished blocks are the two passes’ coefficients carried through that product. A code that wants to solve with the factorisation has to accumulate that R, and it inherits whatever the inner factorisation got wrong in the first pass: the second pass corrects the basis, not the record of how the basis was reached. Every figure below uses R accumulated this way, and the factorisation residual ‖QR − A‖ ÷ ‖A‖ is measured for each variant alongside its orthogonality.

Block modified Gram–Schmidt inherits the rescue #

With the ill-conditioning between the blocks, block modified Gram–Schmidt through $Q^{\mathsf T}b$ is wrong by 1.3·10⁻⁵ at κ = 10⁶, 8.9·10⁻³ at 10⁸ and has no correct digit from 10¹⁰. The same Q and R with b appended as a block give 3.0·10⁻¹², 6.9·10⁻¹⁰ and 2.2·10⁻⁶ at 10¹² — against Householder’s 1.5·10⁻¹¹, 1.4·10⁻⁹ and 1.3·10⁻⁵. The appended block is not merely rescued to Householder’s level; like the column case, it lands a little below it.

Turn the dial. On the generic matrix the appended block gives 1.1·10⁻⁹ at κ = 10⁸ against $Q^{\mathsf T}b$’s 0.23; with the ill-conditioning inside the blocks, 5.6·10⁻⁹ against 7.3·10⁻². At every placement and every κ from 10⁶ up, the appended block is within ten times Householder and the $Q^{\mathsf T}b$ route is at least three orders worse. The rescue does not care which level the ill-conditioning is at, which is the first difference from the preceding essay’s orthogonality curves, where each level mattered.

The reason is the same as for columns, lifted a level. A block modified outer loop removes the finished blocks’ components from b one block at a time, updating b between them, so each block’s coefficients are computed against what was left of b after the previous blocks — the sequential structure Björck and Paige’s argument uses. Inside each finished block the coefficients are dot products against the block’s own columns all at once, which is classical within the block; but each block’s Q came from Householder and is orthogonal within itself to rounding, so projecting against it all at once commits no error that projecting one column at a time would avoid. The loss of orthogonality that ruins $Q^{\mathsf T}b$ is between blocks, and between blocks the appended block is projected sequentially.

The block version is not quite the column version. Column-wise modified Gram–Schmidt on [A b] gives 2.7·10⁻¹⁰, 1.9·10⁻¹⁰ and 1.7·10⁻⁹ at κ = 10⁸ on the three placements; the block version gives 6.9·10⁻¹⁰, 1.1·10⁻⁹ and 5.6·10⁻⁹, between two and a half and six times worse. Both are inside Householder’s range at every point — Householder gives 1.4·10⁻⁹, 2.8·10⁻⁹ and 5.1·10⁻⁹ — so the difference does not change what a code should do. It is the price of projecting against four columns at once inside each block, and it is small because those four columns are orthogonal among themselves.

Block classical Gram–Schmidt gains nothing #

For block classical Gram–Schmidt the ratio is exactly one, at every κ and every placement. Its outer loop projects a new block against all the finished blocks at once, $Q(Q^{\mathsf T}W)$, so the appended b is projected as Q($Q^{\mathsf T}b$) and its column of R is $Q^{\mathsf T}b$, the same numbers in the same order. Two Gram–Schmidts located the whole difference between the classical and modified algorithms in which vector each projection is taken against, and appending b as a block is one more projection; the classical loop takes it against the original.

For block modified Gram–Schmidt the ratio grows about as κ(A) — from a few hundred at 10² to 6·10¹² at 10¹² on the generic matrix. That is the ratio of $\kappa^2$ to κ, the $Q^{\mathsf T}b$ route’s $\kappa(A)^2u$ error over the appended block’s $\kappa(A)u$. It is the column result, unchanged by blocking.

The two-pass variants, and what the appended block adds to an orthogonal Q #

The two-pass Householder variant is the preceding essay’s safe choice, orthogonal to 3·10⁻¹⁵ on every placement, and as a solver it is safe too: through $Q^{\mathsf T}b$ it is within a few times Householder at every κ and placement. The appended block still helps it a little at the top of the range. At κ = 10¹² on the generic matrix its $Q^{\mathsf T}b$ route gives 2.0·10⁻⁵ and its appended block 2.7·10⁻⁶, against Householder’s 3.2·10⁻⁵; between the blocks, 3.8·10⁻⁵ against 3.6·10⁻⁶. An orthogonal Q is not the whole of $Q^{\mathsf T}b$’s error, since forming $Q^{\mathsf T}b$ in floating point still rounds each product, and the appended block’s sequential projections round differently. The effect is a factor of up to ten and it is in the appended block’s favour everywhere it appears.

The most orthogonal Q, and the worst answer #

The two-pass variants hardly need the appended block, since their Q is orthogonal to rounding and multiplying by its transpose is safe. With Householder inside, both routes stay within a few times Householder’s error everywhere. With Cholesky QR inside, on the generic matrix and between the blocks, the same is true.

With the ill-conditioning inside the blocks it is not, and nothing in the orthogonality measurement said so.

At κ = 10⁶ the variant’s Q is orthogonal to 1.2·10⁻¹⁵ — as orthogonal as a Q can be in double precision. Its least-squares answer is wrong by 1.2·10⁻⁵, a hundred and thirty thousand times Householder’s 9.3·10⁻¹¹, and appending b as a block changes it by nothing. From κ = 10⁸ the Cholesky factorisation stops outright on at least one of three seeds, which is the breakdown the preceding essay measured.

The third line explains it. The factorisation residual, ‖QR − A‖ divided by ‖A‖, is 4.6·10⁻¹¹ at κ = 10⁶, where Householder’s is at rounding. Cholesky QR computes R as the Cholesky factor of the Gram matrix of an ill-conditioned block, and that Gram matrix has lost half its digits before the factorisation starts; the second pass reorthogonalises Q perfectly and carries R’s error forward, since it multiplies R by the second pass’s nearly-identity factor. A least-squares solve back-substitutes with R, and an R that reproduces A only to 4.6·10⁻¹¹ gives an x wrong by κ times that, about 10⁻⁵, which is what is measured.

It also failed earlier than the breakdown suggests. At κ = 10⁴, four orders of κ before Cholesky QR stops, the variant’s answer is already 3,400 times Householder’s, and at 10² it is forty-five times. The preceding essay’s orthogonality curve for this variant is flat at rounding up to the breakdown and then ends; a solve that trusted that curve would have been wrong by three orders of magnitude on a matrix whose condition number is ten thousand.

Set beside the solve figure, the orthogonality curve is right about everything it measures and silent on the one thing a least-squares code needs. Its two-pass Cholesky line is flat at 10⁻¹⁵ until κ = 10⁸ and absent after, and the only warning it gives is the absence. The solve’s error had meanwhile been growing as κ(A)²: 2.7·10⁻¹³, 2.3·10⁻⁹ and 1.2·10⁻⁵ at κ = 10², 10⁴ and 10⁶, eight and then five thousand times per hundredfold step, where Householder’s grows a hundredfold. That is the signature of a Gram matrix formed and factorised — a $\kappa^2$ that never appears in Q, because the second pass removes it from Q and leaves it in R.

So an orthogonal Q is half of a factorisation. The preceding essay drew orthogonality curves, and on them this variant looked like one of the two safe ones until it broke down. As a solver it had stopped being safe three orders of κ earlier.

Orthogonality is not the axis #

Drawn together, the two quantities are unrelated in the way that matters. Block modified Gram–Schmidt’s dots run across ten decades of orthogonality, from rounding to 10⁻⁴, and sit on the line of Householder’s accuracy the whole way, because its appended block reads the right thing out of an imperfect Q. Block classical Gram–Schmidt’s dots rise with its lost orthogonality, because it has no route that avoids multiplying by its imperfect Q. And the two-pass Cholesky variant’s dots sit at the far left, orthogonal to rounding, with the three from the inside placement at forty-five, 3,400 and 130,000 times Householder’s error.

The measure that would have separated them is the factorisation residual together with the route. A block QR is fit for least squares when QR reproduces A to rounding and the right-hand side is taken through the factorisation’s own sequential projections, or when Q is orthogonal and QR reproduces A. Orthogonality alone is neither condition.

What it costs, and which to choose #

One pass of block modified Gram–Schmidt with b appended is the cheapest variant that solves every problem here as well as Householder. It costs one between-block pass and one Householder factorisation per block, where the two-pass Householder variant, which also solves everything, costs two of each. What it gives up is communication: a modified outer loop projects against the finished blocks one at a time, a synchronisation per block, where a classical loop does it in one matrix product. How wide the block should be and an orthogonalisation nobody calls one are the essays where that trade is priced; a code that is bound by synchronisations will prefer the classical loop twice, and a code that is bound by arithmetic will prefer the modified loop once.

In arithmetic, for m rows and n columns in blocks of p, one pass of the between-block projection costs about 2mn² and the Householder factorisations inside the blocks about 2mnp more; the second pass of the two-pass variants doubles the first term. At the sizes here that is roughly 33,000 against 66,000 multiply-adds for the projections, a factor that grows toward two as the blocks narrow relative to the matrix. Appending b adds one column’s worth of projection, about 2mn, to either.

The preceding essay’s rule was to decide which level the ill-conditioning is at and pay for stability there. For a least-squares solve there is a cheaper rule and a sharper warning. The cheaper rule: if b is appended as a block and the outer loop is modified, the level does not matter, because the appended block reads the answer out of the factorisation at Householder’s accuracy wherever the ill-conditioning sits. The warning: doing it twice restores the orthogonality of Cholesky QR’s Q and does not restore its R, and a solve needs R.

The placement that breaks the Cholesky variant is not an exotic one. The preceding essay pointed out that a block Krylov method builds its blocks from consecutive powers of the operator applied to a starting block, and consecutive Krylov vectors are nearly parallel by construction — which is ill-conditioning inside the blocks, arriving through the ordering of the columns. That is exactly where a cheap inner factorisation is most tempting, since the blocks are small and a Gram matrix is one matrix product, and exactly where this measurement says it turns an orthogonal basis into a wrong answer. A code that orthogonalises a Krylov block with Cholesky QR twice, and then solves a small least-squares problem with the R it accumulated, is solving with the factor that carries the $\kappa^2$.

The modified variant with the appended block has no such hazard, because its inner factorisation is Householder and its right-hand side reads the answer out sequentially. Its hazard is the one already named: synchronisations. Neither the orthogonality measurement nor this one says which of those costs a given machine will feel more.

What this does not settle #

Consistent problems only. The right-hand side as one more column found that a residual puts even the stable routes on a $\kappa^2$ law, because the least-squares problem’s own conditioning contains $\kappa^2$ times the residual; the block routes will meet the same law, and whether the appended block keeps its margin over $Q^{\mathsf T}b$ there is unmeasured.

One size, one block width, three seeds a point. The block width matters to the between-block loss — with two blocks there is one projection and classical and modified coincide — so the rescue’s size at other widths is a measurement, not an inference.

The account of why the appended block works for the modified loop is Björck and Paige’s argument lifted a level, and it is an argument. A block form of their padded-matrix equivalence, checked numerically as the column essay checked the original, would make it a measurement.

Still open: a residual, and the tree #

An inconsistent right-hand side. Adding a component of b outside the range of A brings in the $\kappa^2$ term that no route avoids. Whether the appended block and $Q^{\mathsf T}b$ then converge to the same error, as the column routes did, and at what size of residual the appended block stops being worth its synchronisations, is the measurement a least-squares code actually faces.

A tree instead of a sequence. A reduction that changes the order measured a tall-skinny QR done as a tree of small factorisations, which is how a distributed code orthogonalises a block. A tree has no sequential projection of b at all, and whether appending b to the tree’s leaves inherits Householder’s stability — each leaf’s factorisation is Householder — or loses it in the combining steps is the block question asked one level lower.