Doing it twice

The message and the word put three factorisations of one matrix on two axes and found they rank differently on each: 48, 4 and 4 communication rounds against 1,170, 1,170 and 2,160 words. The method with the fewest rounds — Cholesky QR, which forms G = AᵀA in a single all-reduce and factorises it locally — pays for that with the site’s oldest currency.

κ(AᵀA) = κ(A)², which is the road that squares the problem arriving in a factorisation rather than in a least-squares solve. The implied orthogonality of Cholesky QR’s Q has a fitted slope of 1.95 in κ where the Householder sweep and the reduction tree stay at 1.00, and past κ ≈ 3·10¹⁰ the Gram matrix is numerically indefinite and the factorisation stops.

The repair is one line.

The one line #

Q₁R₁ = A by Cholesky QR. Then Q₂R₂ = Q₁ by Cholesky QR again. Take Q = Q₂ and R = R₂R₁.

Fitted slopes: 1.95 for one pass, 0.93 for two. The square is gone, and what is left is κu, which is what a backward-stable factorisation costs and what the other two methods have.

The second column is not merely improved. At every condition number in the table it is at or below the Householder sweep’s — the sweep being the reference implementation of a backward-stable QR, and the method Cholesky QR was supposed to be a cheap and inferior substitute for.

Both of those hold at every width the figure will draw, and the third column is the one the field is about.

At n = 6, 8, 12 and 16 the one-pass slope reads 1.96, 1.99, 1.98 and 1.97 and the two-pass slope 0.79, 0.96, 0.95 and 0.91. Two and one, at every width — the square removed and κu left, with no width at which the repair stops working.

And the twice-run beats the Householder sweep at all four: 2.32·10⁻¹⁰ against 4.91·10⁻⁹, 3.75·10⁻⁹ against 8.52·10⁻⁹, 1.79·10⁻⁹ against 3.24·10⁻⁸, and 2.08·10⁻⁹ against 1.21·10⁻⁸ — by factors of 21, 2.3, 18 and 5.8.

The rounds are the column that says why any of this is worth doing. Across the four widths the communication rounds read 3, 6, 18, 3 at six columns; 3, 6, 24, 3 at eight; 3, 6, 36, 3 at twelve; and 3, 6, 48, 3 at sixteen. Three of the four methods are constant in the width and one of them is 3n.

So the trade the previous essay set up resolves cleanly at every width rather than at one. The Householder sweep’s round count is proportional to the number of columns; Cholesky QR’s is three whatever the shape; running it twice is six; and the reduction tree matches it at three. Doubling the work of the cheapest method in rounds still leaves it eight times cheaper than the sweep at sixteen columns, and more as the matrix widens — while returning a Q that is better than the sweep’s at every condition number and every width measured.

Why a second pass is enough, and a third is not #

The mechanism is short enough to state and is the reason this is not a trick.

The first pass returns a Q₁ whose orthogonality error is of order κ²u. That means κ(Q₁) is about 1 + κ²u rather than κ: it is nearly orthogonal, and how nearly is exactly what the first pass got wrong. The second Cholesky QR squares that condition number — and squaring a number that is nearly one leaves a number that is nearly one.

So the second pass has an easy problem precisely because the first pass had a hard one, and the composition lands at u rather than at κ⁴u or at anything else the arithmetic might have suggested. That is the same mechanism as re-orthogonalising a Gram–Schmidt pass, where twice is enough for the modified variant for exactly this reason and not for the classical one.

The two condition numbers, measured rather than asserted:

The first two columns agree to three digits over twelve decades, which is the squaring happening exactly as the algebra says — and the last column is the reason the second pass is cheap. The matrix handed to it has a condition number of one, to four decimal places, at every κ where the first pass returns anything usable at all. At κ = 10⁸, where the first pass’s Q is 0.37 from orthogonal, κ(Q₁) is 1.26.

Squaring 1.26 is 1.59. Squaring 10⁸ is 10¹⁶, which is past 1/u and is where the method stops.

A third pass has nothing left to do, and that is the refusal this essay’s library publishes: feeding the build the claim that a third pass improves the accumulated orthogonality by an order of magnitude must fail, because the second pass has already reached the floor a backward-stable factorisation has.

What it costs, counted #

One more all-reduce of n² numbers.

At p = 8 and n = 8. The sweep needs one reduction per column, so its rounds grow with n; the two Cholesky runs need one and two whatever n is.

A quarter of the sweep’s rounds for the same orthogonality is the headline, and the honest sentence beside it is that the reduction tree spends 3 rounds and 252 words and reaches the same place. On both counts this file measures, the tree is the better method.

What Cholesky QR twice has instead is its shape: two Gram matrices, two Cholesky factorisations and two triangular solves, every one of them a dense matrix–matrix operation, with no tree, no user-defined reduction operator and no dependency between processors beyond two collective calls that every library already implements. That is an argument about what a machine is fast at rather than about what it counts, and this essay does not make it — it counts.

The boundary is the first pass’s boundary #

Running it twice cannot repair a first pass that refused, and this is where the honest reporting has to be careful, because the first pass’s failure is not a threshold.

It refuses at 3·10¹⁰ and 10¹¹ and returns an answer at 3·10¹¹ and 10¹² — an answer on a worse-conditioned matrix than one it declined, which is the erratic boundary the communication field already recorded for one pass and which two passes inherit unchanged.

And where it does return past the cliff, the second pass is better than the Householder sweep: 1.2·10⁻⁶ against 3.5·10⁻⁵ at κ = 3·10¹¹. The first pass’s Q there is 1.62 from orthogonal — worse than useless — and the second pass turns that into an answer an order of magnitude better than the reference method’s.

That combination is worth sitting with. The method either declines or returns something good, and which of the two it does is not predictable from κ. A code that treats a refusal as a signal to fall back is doing the right thing; one that treats κ as a guide to whether to try is not — and κ is what an estimate would report, which makes the guide both cheap and useless here.

Every method reconstructs, including the one whose Q is worthless #

The site’s rule earns its place here more than anywhere.

At κ = 10⁸ all four methods reconstruct A to about 10⁻¹⁶:

The residual — the quantity anybody checks a factorisation with — is identical to two significant figures across a method whose Q is 37% wrong and three whose Q is right to nine digits.

That is the refutation this essay publishes, and it is not a subtle one: a residual says nothing about orthogonality. It is the same gap a small residual is not a small error measures for a solve, sharpened by the fact that here the residual is not merely uninformative but structurally guaranteed. The reason is structural rather than numerical. Q = AR⁻¹ by construction, so QR = A whatever R is, as long as R is invertible; the residual is testing the triangular solve and nothing else.

Two passes, and where the idea comes from #

The construction has a name and a literature and this essay measures it rather than deriving it, so it is worth saying what was already known and what was checked.

Known: that Cholesky QR’s orthogonality degrades like κ², that a second pass restores it while κ²u < 1, and that the composition is what the published analysis calls CholeskyQR2. None of that is a finding here.

Checked: the two slopes, on this site’s own matrices, with the site’s own implied-orthogonality measure — 1.95 and 0.93 rather than 2 and 1. The claim that the second pass reaches the Householder sweep’s accuracy rather than merely improving, which the table above makes precise and which is stronger than the usual statement. The round and word counts, from the same network model the communication field uses for the other three methods, so the four are comparable. The behaviour past the cliff, which is where the published analysis stops and where the erratic refusals live. And the residual, which is identical across a method whose Q is 37% wrong and three whose Q is right.

The last of those is the one this site adds. Every published comparison of these methods reports orthogonality; the reason the residual belongs beside it is that the residual is what a reader would otherwise check, and it says nothing at all.

What the second pass does not fix #

Three things, and the third is the one that would matter in practice.

It does not make the residual mean more. Both passes produce a Q from a triangular solve, so QR = A either way and the reconstruction is uninformative either way. A code that wants to know whether its Q is orthogonal has to measure ‖QᵀQ − I‖ — orthogonal is a number is the essay about turning the adjective into that measurement — which costs a product the factorisation was trying to avoid.

It does not extend the range. The cliff is where the first pass’s Gram matrix goes indefinite, and that is a property of κ(A)²u alone. Two passes have exactly the range one pass has.

And it does not help a matrix that is rank-deficient rather than ill-conditioned. Everything above is measured on matrices with a prescribed condition number and full numerical rank. A tall-skinny matrix with an exactly duplicated column has a singular Gram matrix, the first Cholesky refuses at that column, and there is nothing for a second pass to be run on — where the reduction tree and the sweep both return an R and let the caller decide what its smallest diagonal entry means.

That reads as the practical reason the tree stays the default in libraries, and it is worth measuring rather than stating, because it is the essay’s only argument that is not a count.

What a near-dependent column actually does #

Duplicating a column exactly reproduces the claim: Cholesky QR refuses at column 7, and the sweep and the tree return an R whose smallest diagonal entry is 3.8·10⁻¹⁶ and 8.1·10⁻¹⁷. Note already that neither is a zero — it is a rounding, and reading it as a rank deficiency is a decision rather than an observation.

Now perturb the duplicate by a relative 10⁻¹⁶. The smallest singular value is 6.0·10⁻¹⁷ of the largest, which is no less deficient than before, and Cholesky QR returns. The refusal is triggered by exactness rather than by rank deficiency, and exactness is not what floating-point data contains: a column that is a measured multiple of another, or a duplicate that has been through one arithmetic operation, is not exact.

What all four return across the near-deficient range:

perturbation    Cholesky ×1    ×2         sweep       tree
   10⁻¹⁶          1.00        0.796      0.995       0.923
   10⁻¹⁴          1.00        7.6·10⁻³   6.5·10⁻²    1.8·10⁻²
   10⁻¹²          1.00        2.1·10⁻⁴   1.4·10⁻³    6.0·10⁻⁴
   10⁻¹⁰        refused      refused     1.6·10⁻⁵    7.5·10⁻⁶
   10⁻⁸           0.945       9.9·10⁻⁹   1.5·10⁻⁷    3.4·10⁻⁸

No method returns an orthogonal Q, including the two that never refuse: the sweep’s is 0.995 from orthogonal at 10⁻¹⁶ and the tree’s 0.923. The tree is not letting a caller decide anything the sweep is not; both are returning a Q that is wrong and an R diagonal that has to be thresholded.

And two passes give a better Q than the sweep at every level where they return — by five to fifteen times, which is the same relationship the well-conditioned table shows and which the deficiency does not disturb.

So the closing argument does not survive its own measurement. What survives is narrower and is still worth having: Cholesky QR is the only one of the four that ever says no. On the exact duplicate it refuses loudly where the other three return a plausible answer quietly, and a loud refusal is worth something even when it fires on a case that never arises — because the alternative, on the row above it, is three methods returning a Q that is 92% wrong with nothing in their output to say so. The refusal is erratic in the perturbation just as it is in the condition number, firing at 10⁻¹⁰ and not at 10⁻⁸, so it cannot be relied on; it can only be believed when it happens.

Where this sits among the site’s second passes #

Running something twice to buy back accuracy is a shape this site has met three times now, and the three are worth laying beside each other because they are not the same trade.

Iterative refinement solves the system again against the residual of the first solve, and buys back the digits a low-precision factorisation lost — but only while κu < 1, and the mixed-precision field measured the threshold at exactly that boundary.

Reorthogonalising a Gram–Schmidt step projects a second time against the vectors already orthogonalised against, and takes ‖QᵀQ − I‖ from 10⁻² to 10⁻¹⁶ — the twice is enough, measured table the standard pass added to the orthogonality field.

And this. A second Cholesky QR, which is the same idea one level up: the object being re-processed is a whole factorisation rather than a vector or a right-hand side.

All three work for one reason. The first pass produces something whose error is a known function of the input’s conditioning, and the second pass sees an input whose conditioning is that error rather than the original problem’s. Composing a squaring with a nearly-one leaves a nearly-one — which is why the answer is always twice, in all three cases, and never three times.

When to use which #

The measurements support a short table rather than a recommendation, and the entries are about what a machine is slow at rather than about the matrices.

If the reduction is cheap — a small machine, a fast network, a factorisation that is not on the critical path — the Householder sweep is the reference and there is no reason to leave it. It costs n reductions and nothing else on this page beats it for simplicity.

If the reduction is expensive and κ is under about 10⁷, two Cholesky passes give the sweep’s orthogonality for two collective calls, in an arithmetic that is entirely matrix–matrix.

If κ is unknown, the refusal is the feature. The first pass declines rather than lying — usually — and a code that falls back on a refusal is correct; a code that trusts κ to predict the refusal is not, because the boundary is erratic and it returned an answer at 3·10¹¹ having refused at 10¹¹.

And if the reduction is expensive and the matrix might be badly conditioned, the reduction tree does everything this method does with fewer rounds and fewer words, and it never refuses. That is the conclusion this essay’s own counts point to, and the argument against it is one about implementation shape rather than about arithmetic — which is the honest place to leave it.

What the residual is measuring instead #

One more observation about the residual table, because it explains why every number in that column is the same.

Q is formed as AR⁻¹ in all four methods — that is what impliedOrthogonality does, and it is what a caller who solves a least-squares problem through R is implicitly doing. So QR is A R⁻¹ R, and the residual ‖A − QR‖ is measuring the accuracy of a triangular solve followed by a triangular multiply, which is a backward-stable pair whatever R happens to be.

That makes the residual a test of the triangular arithmetic and nothing else. It would catch an R with a garbage entry, a wrong dimension, a transposed factor. It cannot catch an R that is a perfectly good triangular factor of the wrong Gram matrix, which is exactly what Cholesky QR returns past its cliff.

The measurement that does catch it is ‖QᵀQ − I‖, and it costs a product the whole family of methods was designed to avoid — n² inner products of length m, which is the same order as the factorisation itself. That is the reason the check is not routinely made, and the reason this site prints it on every figure that draws a factorisation.

What the drag does #

The slider is the number of columns, and it moves one thing and leaves another alone.

What moves is the sweep’s cost: one reduction per column, so at n = 16 the Householder sweep spends 48 rounds where the two Cholesky passes still spend 6. The advantage the second pass has in synchronisations grows linearly with the width.

What does not move is either slope. κ² is a property of forming AᵀA and κ is a property of not having to, and neither depends on the shape of the matrix — which is why the figure asserts the two slopes at every position of the slider rather than at one.

The slider starts at six columns rather than four, and the reason is the cliff. At four columns the Gram matrix of a κ = 10⁸ matrix is already numerically indefinite, the first pass refuses, and the curve cannot be drawn through a refusal — which is the same boundary the table above is about, arriving in the figure’s own assertion.

A second pass that is not available #

Re-orthogonalising repairs a loss of orthogonality because there is something to re-project against. Removing a rank-one term has no such remedy: a hyperbolic sweep is not orthogonal, and its coefficients have no bound.