The plateau is set by the length

The factor nobody forms measured a Householder factorisation’s orthogonal factor the way libraries use it: not multiplied out into a matrix, but applied to vectors through its stored reflectors. Its departure from orthogonality — the worst $\|Q^{\mathsf T}(Qx) - x\|$ over unit vectors $x$ — was 3.1 to 4.0·10⁻¹⁵ at every block size from one to sixty-four, and one reflector at a time it grew over the first twenty or so reflectors and then stopped, instead of adding up as the reflectors’ individual errors would.

The essay explained the stop by the operator’s worst direction: each reflector moves a vector in one direction only, so the direction that is worst after the first few reflectors is touched by only a few more, and the rest leave it alone. Its last section drew the consequence for a change of size. “Twenty reflectors on 96-long vectors is where the growth stopped here. The argument above predicts the saturation count depends on the vector length and not on the number of reflectors, so that on 384-long vectors the growth would continue longer before stopping. That is one change of dimension in the same experiment.”

It is one change of dimension, and the curve does change with it — but not in the way predicted. The length moves the height of the plateau, not the place it starts.

One reflector at a time, never formed #

The factorisation is the textbook one: for each column, a reflector $I - \beta vv^{\mathsf T}$ that zeroes the entries below the diagonal, applied to the remaining columns before the next is formed. Nothing is multiplied out. The factor $Q$ is the product of the reflectors, and applying $Q$ or $Q^{\mathsf T}$ to a vector means applying the reflectors one after another, each touching the entries from its own column down. After each $j$ of them the departure is measured as before — the largest round trip through the first $j$ reflectors and back, by power iteration on the map $x \mapsto Q_j^{\mathsf T}(Q_j x) - x$, which for an operator this close to the identity converges in a few dozen steps.

The matrices have condition number $10^8$, as in the earlier essay, and five seeds at each shape: 64 columns on 96, 192 and 384 rows, and 128 columns on 384 rows. Forming nothing makes the larger shapes cheap — the earlier essay formed the factor alongside, which at 384 rows would have cost a dense product of order 384 for every reflector.

The count does not lengthen #

The figure at the top of the page is the comparison the prediction asked for, and its answer is visible without a number. All three curves rise quickly and flatten. They flatten at different heights and at about the same place.

To give “the same place” a number, the plateau is read as the median of each curve over its second half — reflectors 33 to 64 — and the count is the first reflector after which the curve has reached 90 per cent of it. On 96-long vectors that count is 11; on 192, 19; on 384, 9. At 80 per cent it is 7, 11 and 3. There is scatter in these numbers, because each curve is a median of five noisy ones, but there is no trend with the length, and the longest vectors reach their plateau soonest. The prediction’s direction — four times the length, a longer climb — is the one direction the data does not take.

The dial shows the three curves one at a time, each against its own plateau. On 384-long vectors the first reflector alone leaves a round-trip error of 1.2·10⁻¹⁵, the eighth already 5.6·10⁻¹⁵, and the curve wanders between 5.4 and 6.8·10⁻¹⁵ for the remaining fifty-six. On 96-long vectors the first reflector leaves 1.8·10⁻¹⁵, the eighth 2.5·10⁻¹⁵, and the climb to 3.6·10⁻¹⁵ is slower and lower. The long vectors do not climb for longer; they climb harder, and stop.

What the length sets #

The plateaus are 3.6·10⁻¹⁵ on 96 rows, 5.0·10⁻¹⁵ on 192 and 6.7·10⁻¹⁵ on 384. Doubling the length raises the plateau by a factor of 1.4 each time, which is $\sqrt 2$. Divided by $u\sqrt m$ — the unit roundoff, $2^{-53}$, times the square root of the length — the three plateaus are 3.28, 3.27 and 3.09, and the 128-column factor on 384 rows is 3.37.

That is the shape an analysis of one reflector predicts, and the same square root that appeared earlier in a plainer setting. The order they are added in and three walks and one bound measured how rounding errors accumulate along a long sum: the worst-case bound grows with the number of terms, the typical error with its square root, because errors of random sign walk rather than march. Applying $I - \beta vv^{\mathsf T}$ to a vector of length $m$ forms an inner product of length $m$ and an update of length $m$, and the rounding error in that inner product, for random-looking data, grows like $\sqrt m$ times the unit roundoff — the typical case, not the worst, which grows like $m$. Each reflector applied contributes an error of that size, and the departure is what the errors leave behind in the round trip. A plateau proportional to $u\sqrt m$ says that what is left behind is a fixed multiple of one reflector’s typical error, whatever the length.

Divided by $u\sqrt m$, the three curves fall on top of each other to within their scatter: between 1.6 and 2.6 after eight reflectors, between 2.5 and 3.8 from thirty-two on, at every length. So the plateau is a fixed multiple of a single reflector’s typical rounding error, about three, and the number of reflectors it takes to accumulate to that multiple is the same — a dozen or two — at every length. The length is in the unit of the curve and nowhere else.

Except at the very start, and the exception is worth stating because it narrows what “one reflector’s error” means. The round trip through the first reflector alone leaves 1.76·10⁻¹⁵ on 96 rows, 1.78·10⁻¹⁵ on 192 and 1.20·10⁻¹⁵ on 384: no growth with the length at all. Divided by $u\sqrt m$ those are 1.62, 1.16 and 0.55, and the scaled curves only come together from about the fourth reflector, where they read 1.77, 1.49 and 1.82. So the plateau is 2.0 times the first round trip on the shortest vectors and 5.6 times it on the longest — the long vectors start from the same place and climb further. The square root is not in any single round trip; it is in what the round trips leave behind once several have been stacked, which is the accumulation read below under why the curve stops, and it is why the fit should be read as a law for the plateau and not for each reflector.

The reflectors’ own count #

The prediction’s other half was that the number of reflectors does not set the saturation, and that half holds. On 384 rows, a factor of 128 reflectors settles at 7.3·10⁻¹⁵ where one of 64 settles at 6.7·10⁻¹⁵ — nine per cent more for twice the reflectors, within the scatter of five seeds — and reaches 90 per cent of its plateau after nine reflectors, the same as the 64-reflector factor. The last hundred reflectors of the longer factor add nothing to the departure. If each added an independent error the departure would grow like the square root of the count, $\sqrt 2$ for twice as many; it grows by a tenth of that.

Why the curve stops #

The earlier essay’s account of the stop was that each reflector touches the worst direction only through its own component along it. That account is right as far as it goes, and the data adds what it leaves out: the stop happens at a count that does not depend on the length, and at a height that does, which says the mechanism is about how errors of one reflector are undone by the next rather than about how many directions there are to fill.

A reflector is its own inverse. In the round trip $Q^{\mathsf T}Qx$, each reflector is applied once on the way and once on the way back, and in exact arithmetic the two applications cancel pairwise from the middle out. In floating point each pair leaves a residue of about one reflector’s rounding error, and the residues of different pairs, made by different vectors at different stages, are largely independent. If they simply added in quadrature the departure would grow like the square root of the number of reflectors without end. It does not, because each residue is a small perturbation of the vector that the remaining reflectors then act on, and a reflector applied to a perturbed vector carries the perturbation through unchanged except in its own direction; the round trip’s later pairs cancel what they do to it. The perturbation that survives is the part each reflector puts in directions the later ones do not reach, and once the first dozen reflectors have spread their residues over most of the directions that matter, later residues land on top of earlier ones rather than beside them.

That reading is an interpretation of the measured curves, not a derivation; but it predicts the two things the measurement shows — a count set by how fast the residues fill the reachable directions, which is a property of the reflectors and not of the length, and a height set by how big one residue is, which is the length’s.

Why a reflector holds orthogonality at all #

The whole measurement rests on a property the first essays on Householder established and is easy to forget once it is taken for granted. A reflection cannot stop being one found that Householder’s factor stays orthogonal to rounding level whatever the matrix’s condition number, because each reflector is built from a single vector and is, as computed, the exact reflection in a slightly different vector — an orthogonal operator, rounded or not. One number that has to be right found the one place that guarantee can fail: the scalar $\beta = 2/v^{\mathsf T}v$, which must match the vector it scales, or the operator is no longer a reflection.

So the departure measured here is not a reflector failing to be orthogonal. Each reflector, as stored, is orthogonal to within a rounding of its own norm. The departure is the arithmetic of applying it: the inner product $v^{\mathsf T}x$ and the update $x - \beta(v^{\mathsf T}x)v$, rounded at every entry of a vector of length $m$. That is why the length is what sets the plateau. A taller matrix has longer reflectors, and every application of a longer reflector rounds more numbers.

What this means for a caller #

Eight blocks and sixty-four reflections found the formed factor’s departure growing with the number of reflectors, and the factor nobody forms found the applied one flattening. The plateau’s dependence on the length, measured here, finishes the account a caller needs. A least-squares solve with a tall matrix applies $Q^{\mathsf T}$ to its right-hand side once; the orthogonality it relies on is lost to about $3u\sqrt m$, set by the matrix’s height and nothing else. Doubling the number of unknowns, at a fixed height, costs nothing in orthogonality. Doubling the height costs a factor of $\sqrt 2$, which at any height a double-precision code meets is still far below the conditioning errors of the problem itself.

Orthogonal is a number began this field by insisting that orthogonality be measured rather than assumed. The measurement here gives it a formula for the applied Householder factor that a backward-error analysis would not: the standard bound grows like $mn\,u$ — the product of the dimensions — where the measured loss grows like $\sqrt m\,u$ and is independent of $n$. The bound is a worst case over every possible accumulation of errors; the data says the accumulation stops after a dozen reflectors, and that the worst direction is not the worst case.

The same distinction matters for anyone choosing between Householder and Gram–Schmidt for a tall matrix. Two Gram–Schmidts showed the modified algorithm’s loss of orthogonality growing with the condition number; Householder’s does not, and what it does grow with — the square root of the column length — is mild enough that a matrix a million rows tall would lose about three thousand times $u$ on that law, still twelve digits of orthogonality.

Blocks change the order, not the count #

Libraries do not apply reflectors one at a time; they group them into blocks and apply each block as $I - YTY^{\mathsf T}$ with two matrix products. The earlier essay found that this changes nothing in the applied departure — 3.1 to 4.0·10⁻¹⁵ at every block size on 96 rows — and the measurement here explains why it should not. A block of $b$ reflectors applied through $Y$ and $T$ forms $b$ inner products of length $m$ at once, then a small triangular product, then $b$ updates of length $m$: the same number of length-$m$ roundings as the reflectors applied one at a time, grouped differently. The residue each leaves has the same size, $u\sqrt m$ times a constant, and the plateau is set by that size and by the dozen or so of them it takes to fill the directions that matter. On 384 rows the prediction this makes is a plateau near 6.7·10⁻¹⁵ at every block size, as at 96 rows the plateau was the same at every block size; that is a statement of the measurement’s logic, not a measurement, and it is the cheapest of the predictions here to check, since the blocked application at 384 rows costs no more than the one already run.

The formed factor is a different object, and the earlier essays found its departure growing with the number of reflectors. Forming $Q$ multiplies the reflectors into an $m \times m$ matrix, and every entry of the product is a sum over all of them; the cancellation the round trip enjoys, each reflector meeting itself on the way back, never happens. So a caller who forms $Q$ pays for the number of reflectors, and one who applies it pays for the length. On a tall, thin matrix — many rows, few columns — the two differ little. On a square one the applied factor’s advantage is the whole of the difference between $\sqrt m$ and the formed factor’s growth in $n$.

What five matrices and three lengths do not show #

The lengths stop at 384, and the $\sqrt m$ fit rests on three of them a factor of two apart; a length of ten thousand might show the worst-case growth beginning to matter, though nothing in the three fitted constants, 3.28, 3.27 and 3.09, drifts upward. The matrices are dense random ones of a single condition number; a structured matrix whose columns share a common direction could make the residues correlated, and correlated residues would add rather than land on top of each other. The plateau is a median over the second half of each curve, and the reach counts read off medians of five are noisy enough that 9, 11 and 19 should be read as “about a dozen, with no trend” rather than as three different numbers. And all of it is in IEEE double precision, with the default rounding to nearest and no fused multiply-adds anywhere in the arithmetic; the $\sqrt m$ law is a statement about rounding errors of random sign, and in a lower precision, or with rounding that is not to nearest, the sign randomness is the first thing to test.

Still open: rounding that does not average, and a column of shared direction #

Rounding toward zero. Round-to-nearest gives errors of random sign, and $\sqrt m$ is what random signs add up to over an inner product of length $m$. Directed rounding — every operation rounded toward zero, as some low-precision hardware does — gives errors of one sign. The prediction with a sign is that under rounding toward zero the plateau grows like $m$ rather than $\sqrt m$, by a factor near 2 rather than 1.4 for each doubling of the length, while the count to reach it stays near a dozen.

Columns that share a direction. If every column of the matrix has a large component along one fixed vector, every reflector’s residue lies partly along it, and the residues stop landing in independent directions. The prediction is that on matrices whose columns are a shared unit vector plus a small random part, the plateau grows with the number of reflectors — like its square root — instead of stopping, so that the stop measured here is a property of matrices without a dominant common direction.