The weight the factor met first

The factor a sparse code keeps anyway ended its account of one minus a leverage with a rule in two sentences. “A route that subtracts loses a digit for every decade the leverage is from one, whatever else it does. A route through the orthogonal factor loses a digit for every decade of the condition number, whatever else it does.” The second sentence came from one minus a leverage is a subtraction, which computed 1 − h as the squared norm of a row of the complementary block $Q_2$ of a Householder factor, and on a polynomial fit whose far point made 1 − h small found that number wrong by about κ(A)·u.

Both essays left the weighted case open, and for a reason the rule’s own wording makes pressing. A weight drives a row’s leverage towards one by choice: multiply one row of an ordinary design by $\sqrt{w}$ and 1 − h falls as 1/w while κ(A) rises as $\sqrt{w}$. If the rule holds, a weight of 4³⁰ costs the complement eight or nine digits. It costs none, on one condition, and the condition is not about the matrix.

One row, weighted by powers of four #

The design is a straight line through thirty-one points scattered on [−1, 1], with one of them multiplied by $2^k$, so that its weight is $4^k$ and the scaled row is still exactly the doubles that were stored. The exact 1 − h is computed in rational arithmetic from those doubles, as both earlier essays did. Nothing else about the design is unusual: without the weight its κ(A) is 1.8.

With the weighted row last, the figure is what the rule predicts. At a weight of $4^{15}$, 1 − h is 7.9·10⁻⁹ and κ(A) is 9.4·10³; the two subtracting routes are wrong by 4.1·10⁻⁹ and 3.8·10⁻⁸, which is u/(1 − h) as always, and the complement by 6.5·10⁻¹⁴. At $4^{30}$, 1 − h is 7.4·10⁻¹⁸, below the unit roundoff, κ(A) is 3.1·10⁸, the subtractions have nothing — one returns a negative number and the other the square root of one — and the complement is wrong by 2.2·10⁻⁹, a little under κ(A)·u. Seven digits kept of a quantity the other routes cannot see at all.

Turn the dial to the first position. The rows are the same thirty-one rows at the same weights, and only the order in which the factorisation meets them has changed. The subtracting routes do not notice: their error is the same u/(1 − h) in every frame, because a leverage computed from the Gram matrix or from $Q_1$ is order-blind. The complement is flat. At every weight from 1 to $4^{30}$ it is right to within 6.3·10⁻¹⁶, and at 4³⁰ its relative error is 2.1·10⁻¹⁶ while κ(A)·u is 3.4·10⁻⁸. The middle position looks like the last.

Against the condition number, first and elsewhere #

Drawn against the condition number, the two behaviours separate completely. On a straight line and on a cubic, with the weighted row first, the complement’s error never exceeds 1.2·10⁻¹⁵ while κ(A) climbs to 2.5·10⁹. With it in the middle or last it rises with κ(A), at a tenth to a hundredth of κ(A)·u: 2.2·10⁻⁹ for the line at κ(A) = 3.1·10⁸, and 6.8·10⁻⁹ for the cubic at 2.5·10⁹.

So on a single weighted row, κ(A)·u is an upper bound on the unsorted error and irrelevant to the sorted one. That is already enough to say the rule was stated too broadly: the same matrix, up to a permutation of its rows, is either κ(A)·u wrong or a unit of roundoff wrong, and a permutation of the rows does not change κ(A). But on one weighted row the condition number and the weight grow together, so the figure cannot yet say which of them the unsorted error is following. Two weighted rows can.

Two heavy rows, and a condition number of 3.1 #

Weight two rows instead of one. If both are heavy, they carry the whole of a straight line’s two-dimensional column space between them, and the light rows are a small perturbation of a well-conditioned two-by-two system. With weights $4^{25}$ and $4^{24}$ the design’s κ(A) is 3.1. Every account of the complement’s error in terms of conditioning says it should be right to the last bit, in any order.

Sorted, heavier first, it is exact. Reversed, lighter heavy row first, it is wrong by 4.1·10⁻¹⁶. With both after the light rows it is wrong by 2.4·10⁻¹⁰ — six digits gone, at a condition number of 3.1, where κ(A)·u is 3.4·10⁻¹⁶. The other three pairs say the same thing where κ(A) is larger: at weights $4^{30}$ and $4^{20}$, κ(A) is 1.5·10³ and κ(A)·u is 1.7·10⁻¹³; sorted the error is 6.1·10⁻¹⁶, reversed 2.0·10⁻¹⁴, and with the light rows first 6.8·10⁻⁹, forty thousand times κ(A)·u.

The reversed order is worth reading carefully, because it is the one a user is likeliest to produce by accident. The two heavy rows are both above the light ones, but in the wrong order among themselves. The loss is real and it is small — between one and three orders of magnitude where the light-first order loses five or seven — because the rows that pay for it are the second heavy row, which is heavy enough to absorb most of it. The rule that removes all of it is the full sort, by decreasing weight; the rule that removes nearly all of it is “heavy before light”.

What the unsorted error is following #

With the condition number ruled out, the natural candidate is the one a constraint is a weight at infinity found for the solution of the same kind of problem. There, stacking heavily weighted constraint rows below the data rows cost ten digits of the answer at a weight of 10¹⁴, and the mechanism was that Householder QR is backward stable normwise: the perturbation it commits is of size u times the norm of the whole matrix, which is the heavy rows’ scale, and put the light rows after the heavy ones’ reflector has been applied and the perturbation lands on the heavy rows, which can afford it; put them before and it lands on the light rows, which cannot.

The leverage of a heavy row is decided by the light rows: its 1 − h is roughly the light rows’ share of that direction, divided by the weight. So a relative perturbation of size $u\sqrt{w}$ on the light rows is a relative perturbation of about that size on 1 − h. Measured over thirty-two unsorted cases — one weighted row in the middle or last on a line and a cubic, and the four pairs with their light rows first — the complement’s error sits between 8.0·10⁻⁴ and 0.17 times u times the heavy row’s scale, and never above it. The open dots, κ(A)·u for the same cases, agree on the single rows, where the scale and the condition number are the same number, and miss the pairs by up to six orders.

The rule, restated: the complement’s error is u times the scale of the rows the factorisation meets after the ones it is protecting. It is κ(A)·u only when the ill-conditioning is the scale, and even then only when the scale arrives late.

Two neighbouring results say why no condition number could have predicted this. A condition number scaling cannot move measured Skeel’s componentwise condition number, which is invariant under any scaling of the rows because the diagonal factors cancel entry by entry. A weight is a row scaling, so by that measure the weighted design is exactly as well conditioned as the unweighted one — which is the right answer for the sorted factorisation, whose error does not move with the weight at all. The normwise κ(A) sees the weight and predicts a loss; the componentwise one does not see it and predicts none. Neither sees the order, and the order is what decides.

And a reflection cannot stop being one is the reason the loss hides so well. Householder’s computed $Q$ is orthogonal to about 10⁻¹⁵ whatever the matrix, in either order, because each reflector is built from a unit vector and a rounded unit vector is a slightly different reflection rather than a broken one. So the heavy row’s full row of $Q$ has unit norm to working accuracy in both orders. What differs is how that unit norm is shared between $Q_1$ and $Q_2$, and the complement is the small share. Its entries are of size one over the heavy row’s scale, and an absolute error of u in each of them is a relative error of u times the scale — the law the unsorted cases follow.

Why first is free #

The measurement says the first position costs nothing, and the reason can be read off the first reflector, which is where the heavy row’s scale enters the factorisation.

A Householder reflector for the first column is built from that column and its norm. When the heavy row is first, its entry is the pivot: the reflector’s vector is the column with the pivot entry pushed further from zero by the column’s norm, a sum of two numbers of the same sign and so free of cancellation. The heavy row’s own row of the reflector is then, to within rounding in each entry separately, the column divided by its norm — small numbers for the light rows, each computed as a product and each accurate to u relative to itself. Every later reflector starts below the first row and never touches it again except by carrying those small entries through transformations built from the light rows alone, which are ordinary. The complement of the heavy row is assembled from quantities each of which was computed relatively accurately. Nothing was subtracted.

When the heavy row is last, the pivot is a light entry and the heavy entry is somewhere below it. The reflector’s row for the heavy observation now has a diagonal entry of one minus twice its component squared over the vector’s squared norm — and that component is nearly the whole of the vector. The diagonal entry is one minus a number just under one: the same subtraction the complement was meant to avoid, moved inside the first reflector. Its rounding is u in absolute terms and u times the ratio of heavy to light in relative terms, which is the law the unsorted cases follow.

That is an argument, not a measurement: the diagonal entry of the reflector is not computed exactly here and its error is not isolated. What the measurements do show is that nothing else in the factorisation distinguishes the two orders — same rows, same reflector count, same operation count, and a factor of up to ten million between the answers.

The far point is not a weight #

The polynomial fit that established the κ(A)·u law is the counter-case that keeps this from being a slogan. Its far point has a small 1 − h because its row is large in the high-degree columns and ordinary in the constant one — it is not a scaled copy of an ordinary row, and no reordering turns it into one.

First or last, the far point loses digits along κ(A)·u. At x = 1,000, where κ(A) is 6.9·10⁸, it is wrong by 2.7·10⁻⁷ placed first and 2.7·10⁻⁸ placed last, against κ(A)·u = 7.7·10⁻⁸. Putting it first is, if anything, slightly worse. A weighted row at the same condition number, placed first, is right to a unit of roundoff.

So the earlier law stands for the case it was measured on. A design that is ill conditioned because one of its observations sits far from the rest in a direction the others do not span pays κ(A)·u through the orthogonal factor, in any order. A design that is ill conditioned because one of its observations was given a large weight pays nothing, if the factorisation meets that observation first, and pays by its scale otherwise. The first is geometry and the second is bookkeeping.

The control that was sorted by accident #

This also explains a number in the first of these essays that was reported as a baseline and deserves another look. Its first experiment made a leverage near one by scaling one row of a random forty-by-six design, and found the complement right to the last bit at every h up to 1 − 10⁻¹⁵.5, where κ(A) is 3·10⁷. It then observed that “in practice a leverage near one belongs to an observation far from the others,” moved to the polynomial family, and found the complement at κ(A)·u. The contrast was read as: a clean design where nothing is at stake, and a realistic one where the condition number bites.

The scaled row in that first design was row zero. It was sorted by construction. Moved to the last position, the same design’s complement is wrong by 3.1·10⁻¹¹ at h = 1 − 10⁻¹² and 2.7·10⁻¹⁰ at 1 − 10⁻¹⁵. The first experiment was not a clean design: it was a weighted design that happened to be in the one order that makes a weight free, and its perfect score was a statement about the order.

What to compute, and in what order #

The practical rule is short, and it adds one step to the rule the preceding two essays left.

For a weighted fit, sort the rows by decreasing weight before the Householder factorisation, and compute every 1 − h from the complementary block — through the reflectors a sparse code keeps if $Q$ was never formed. The sort is a permutation of the equations, costs nothing but the comparison, changes neither the solution nor the leverages in exact arithmetic, and turns the complement’s error from u times the largest weight’s square root into u. The subtracting routes gain nothing from it; the far point of an unweighted design gains nothing from it; everything that is a weight gains everything.

This is the same sort a constraint is a weight at infinity asked for to get the solution right, and it is worth noticing that the two needs are the same need. A code that sorts for the solution and computes its diagnostics from the same factor gets both. A code that sorts for the solution and then computes its diagnostics by subtraction, from the Gram matrix or from $Q_1$, gets the solution to the last digit and the deletion diagnostic of the heaviest observation to no digit at all — which is exactly the observation, as two observations that hide each other argued for pairs, whose deletion matters most.

It is also the answer to the question the scaling that buys ten orders raises from the other side. Rescaling columns changes κ(A) and not the problem; reordering rows changes neither and still decides the error. A condition number is a property of a matrix up to column scaling and row permutation, and the complement’s error, under weights, is a property of the particular permutation.

What this does not settle #

Straight lines and cubics through thirty-one points, one or two weighted rows, weights that are powers of four so that the scaled rows are exact. A weight that is not a power of four rounds the scaled row, and that rounding is a perturbation of the heavy row of size u relative to itself — harmless in principle, since the heavy row can absorb it, but not measured.

Many weighted rows at many scales — a fit whose weights spread over ten decades across a thousand observations — is the realistic case and is not measured. The rule “sort by decreasing weight” is the natural generalisation and the reversed pair suggests the partial sorts that a real code produces cost one to three orders rather than five to seven, but that is an extrapolation from two rows.

Householder QR only. Givens rotations applied row by row, and the modified Gram–Schmidt route through $Q$, have their own orderings and their own ways of committing a perturbation, and neither is measured. Column pivoting reorders the unknowns rather than the equations and is independent of all of this.

Still open: many weights, and the diagnostics that multiply #

A thousand rows over ten decades. The measurement that would turn the rule into a statement about real weighted fits is a design whose weights are spread log-uniformly over many decades, factorised in the stored order, in a partial sort of the kind a streaming code produces, and fully sorted — with the complement’s error for every row, not only the heaviest. Whether the light rows’ own 1 − h values, which are not small, suffer at all from a late heavy row is also unasked.

The other factors in a deletion diagnostic. Cook’s distance multiplies by h/(1 − h)² and divides by an estimate of the residual variance; the studentised residual divides by the square root of 1 − h times that estimate. Under weights the residual variance is dominated by the light rows and the heavy row’s residual is nearly zero, so the numerator of the diagnostic has a cancellation of its own. Whether the sort that fixes the divisor fixes the numerator, or whether the heavy row’s residual needs the same care the divisor did, is the measurement that would finish the diagnostic rather than its denominator.