The residual the solution cannot hold

The weight the factor met first left one practical rule: for a weighted fit, sort the rows by decreasing weight before a Householder factorisation, and take one minus each leverage from the complementary block $Q_2$ of the orthogonal factor rather than by subtracting a computed leverage from one. Sorted, the heavy row’s 1 − h came back right to a unit of roundoff at a weight of $4^{30}$, where it is 7.4·10⁻¹⁸ and both subtracting routes return nothing.

That rule fixes a denominator. Every deletion diagnostic divides something by 1 − h, and the something is the row’s residual. Cook’s distance for row $i$ is

$$D_i = \frac{r_i^2}{p\,s^2}\,\frac{h_i}{(1-h_i)^2},$$

where $r_i$ is the residual, $s^2$ the residual mean square and $p$ the number of columns. The studentised residual divides $r_i$ by $s\sqrt{1-h_i}$. The same essay closed by asking whether the sort that repairs the divisor also repairs what it divides — because a heavily weighted row is one the fit nearly interpolates, so its residual is as small as its 1 − h is, and a small number computed as a difference is exactly the kind the divisor had been.

It does not, and the reason is not the order. The residual as a textbook computes it is lost before any factorisation is chosen.

A deleted residual in disguise #

It helps to rewrite the diagnostic before measuring it. The deleted residual $d_i = r_i/(1-h_i)$ is the prediction error for observation $i$ from a fit that never saw it, which influence is decided before the data drew beside forty refits and found identical to the closed form. In terms of it,

$$D_i = \frac{h_i\, d_i^2}{p\, s^2}.$$

For a heavy row $h_i$ is nearly one, so Cook’s distance is the squared deleted residual over $p s^2$, and the deleted residual is an ordinary number: a straight line fitted through thirty light points predicts the thirty-first, and misses by whatever the noise and the line’s error there come to. Weighting the row by $w$ multiplies its datum and its design row by $\sqrt w$, so it multiplies the deleted residual by $\sqrt w$ too — at $4^{30}$ it is −3.8·10⁷, which is $2^{30}$ times the −0.035 it is at weight one. Nothing about it is small.

What is small is the pair it is usually computed from. The residual falls as $1/\sqrt w$ and 1 − h as $1/w$, and the quotient of the two is the ordinary number. A route that computes the pair and divides must get both small numbers to the same relative accuracy it wants in the answer. One minus a leverage is a subtraction measured what that costs on the divisor. This essay measures the numerator.

The residual by three routes #

The design is the one the weighted-divisor essay used, with a response added: thirty-one points scattered on [−1, 1], data from a straight line with noise of standard deviation 0.05, and the first point’s row and datum multiplied by $2^k$ so that its weight is $4^k$ and the scaled values are exactly the doubles stored. The exact residual is computed in rational arithmetic from those doubles, by solving the normal equations exactly — the route that is useless in floating point and perfect in rationals.

Three routes to the heavy row’s residual are compared against it.

The first is the textbook one: solve the problem by Householder QR, then form $r_i = b_i - a_i^{\mathsf T} x$.

The second is the same double solution $x$ with the subtraction and every product carried out exactly, in rationals. It has no rounding of its own; whatever error it shows belongs to $x$.

The third never forms a solution. With the full orthogonal factor $Q = [\,Q_1\;Q_2\,]$, the residual vector is the projection of $b$ onto the complement of the column space, $r = Q_2 Q_2^{\mathsf T} b$, and the heavy row’s residual is its row of $Q_2$ times the vector $Q_2^{\mathsf T} b$. It is built from the same block the divisor is built from.

With the heavy row last, all three routes lose digits as the weight grows, which is what the divisor did in that position too. At $4^{15}$ the residual is −9.3·10⁻⁶ and the textbook route is wrong by 8.5·10⁻⁸, the exact subtraction by 3.6·10⁻⁸, and the complement by 7.1·10⁻¹³. At $4^{30}$ the textbook residual has no correct digit, the exactly formed one is wrong by twenty-four times the true value, and the complement is wrong by 2.6·10⁻⁸.

Turn the dial to the first position. The complement’s line falls flat, within a few units of roundoff at every weight — no worse than 3.4·10⁻¹⁴ on the line or the cubic, and 9.1·10⁻¹⁶ at $4^{30}$, where the residual is −2.8·10⁻¹⁰. The two routes through the solution do not move. At $4^{15}$ the textbook residual is wrong by 7.0·10⁻⁷ and the exact subtraction by 3.6·10⁻⁷ — slightly worse than with the row last — and at $4^{24}$ they are wrong by 2.9·10⁻² and 1.6·10⁻². The sort bought them nothing.

What the solution cannot resolve #

The second route is the one to look at, because it removes every excuse. Its subtraction is exact, its products are exact, and it is still wrong by a percent at $4^{24}$. The error is in $x$, and it does not need $x$ to be a bad solution.

A least-squares solution is a vector of doubles, and each of its components carries a rounding of at least half a unit in its last place, however it was computed. The heavy row’s fitted value $a_i^{\mathsf T} x$ is of the size of its datum, $\sqrt w$ times an ordinary number, so a rounding of $x$ in its last place moves the fitted value by about $u\sqrt w$ times an ordinary number. The residual is that fitted value’s distance from the datum, and it is an ordinary number divided by $\sqrt w$. So the relative error in the residual from any double solution — the exact one rounded to doubles included — is about $u\,|b_i|/|r_i|$, which grows as $w$.

That is the loss any subtraction of two numbers of the datum’s size must suffer, and it is the whole of the error. Across thirty-six cases — a line and a cubic, the heavy row first, in the middle and last, and weights from $4^9$ to $4^{24}$ — the textbook residual’s error sits between 0.15 and 23 times $u\,|b_i|/|r_i|$, a band that neither rises nor falls with the weight and that the order does not move. There is no mechanism in that figure except the cancellation.

It is worth stating the consequence plainly, because it is not the usual story about residuals. A small residual is not a small error argued that a residual can be tiny while the solution is wrong. Here the solution is as right as a double can be and the residual computed from it is wrong, because the thing being asked of the solution is to hold a misfit smaller than its own last digit. At $4^{27}$ the heavy row’s residual is 2.3·10⁻⁹ and its datum about $10^8$, a ratio beyond the sixteen digits a double has. No solver, no refinement and no extended-precision residual computation can recover from $x$ a quantity $x$ does not contain.

The divisor had the same law and the same cure. One minus a computed leverage loses $u/(1-h)$, and $1/(1-h)$ grows as $w$; the complement computes it as a sum of the other squares and loses nothing. The residual formed from the solution loses $u\,|b_i|/|r_i|$, which also grows as $w$; the complement computes it as a projection and, with the row first, loses nothing. The two quantities are one problem, and it has one answer: take both of them from $Q_2$.

The quotient, three ways #

With the pieces measured, the deleted residual — the number the diagnostic actually needs — can be compared by the routes a code would use. The textbook quotient divides the textbook residual by one minus the squared row norm of the thin factor $Q_1$. The complement’s quotient divides the complement’s residual by the complement’s 1 − h. And the definition refits the other thirty rows and predicts the heavy one.

The textbook quotient divides a number wrong by about $u\,w$ by a number wrong by about $u\,w$, and gets the sum. At $4^{21}$, where the deleted residual is −7.5·10⁴, it is wrong by 6.2·10⁻³; at $4^{27}$ it has no digit; at $4^{30}$ its divisor has become zero or negative and it returns an infinity. The complement’s quotient is right to 6.2·10⁻¹⁵ at $4^{21}$ and 5.8·10⁻¹⁶ at $4^{30}$.

The refit is right to 8.6·10⁻¹⁵ at every weight and, not shown, in every order. That is worth a sentence, because it says where the difficulty lives. The deleted residual is a well-conditioned quantity: the thirty light rows determine a line, the line predicts the heavy point, and the prediction’s error is $\sqrt w$ times a difference of two ordinary numbers of about the same size. Nothing is small, so nothing cancels. The difficulty is entirely in the shortcut, which reaches an ordinary number as the ratio of two tiny ones. Refitting costs a second factorisation per row examined — for the one flagged heavy row a code has reason to examine, that is affordable; for all thirty-one it is thirty-one factorisations where the complement costs one.

Minus 2.3 standard deviations, from two numbers near zero #

The studentised residual is the diagnostic a reader of a regression table actually looks at, and it shows the problem in its starkest form. It is $r_i/(s\sqrt{1-h_i})$, which for the heavy row settles at −2.3 once the weight passes about $4^6$ and stays there: the observation sits a little over two standard deviations from where the other thirty put it, whatever weight it is given, because the weight scales the numerator and the square root of the denominator by the same amount. A reader sees an unremarkable number. It is the quotient of a residual of 2.8·10⁻¹⁰ and a square root of 7.4·10⁻¹⁸.

Computed the textbook way — residual from the solution, 1 − h by subtraction, $s$ from the same residuals — it is wrong by 2.0·10⁻⁹ at $4^9$, 1.3·10⁻⁴ at $4^{18}$, 2.8·10⁻² at $4^{24}$, and from $4^{27}$ it is not a number at all, because the subtraction under the square root has gone negative. With the heavy row first and every piece from the complement, it is right to 1.8·10⁻¹⁵ at $4^{30}$. With the row last and every piece from the complement it keeps eight or nine digits, wrong by 1.6·10⁻⁹. A table that prints −2.3 in all three cases gives no sign which one it is printing.

What a late heavy row does to the others #

The weighted-divisor essay closed by noting that it had asked about the heavy row only, and that whether the light rows’ diagnostics suffer from a late heavy row was unasked. A diagnostic reads every row’s residual — through $s^2$ if nothing else — so the question matters here.

They do suffer. With the heavy row last, the worst of the thirty light residuals is wrong by 5.8·10⁻⁶ at $4^{30}$, by either route. With the heavy row first it is wrong by 7.2·10⁻¹³, which is the ordinary cancellation of a light residual of about 0.05 against data of about one. So a late heavy row does not merely cost its own diagnostic; it degrades the solution, and a constraint is a weight at infinity measured exactly that on the solution itself. The light rows’ residuals inherit it because they are formed from that solution or from the same perturbed factor.

And yet $s^2$ formed from the damaged textbook residuals is right to 6.3·10⁻¹⁵. That is not luck. The residual of a least-squares problem is orthogonal to the column space, so if the solution moves by $\delta$ the squared residual norm becomes $\lVert r\rVert^2 + \lVert A\delta\rVert^2$: the first-order term vanishes. The sum of squares is stationary at the answer, and an error in $x$ that moves every light residual by a part in $10^6$ moves their sum of squares by a part in $10^{12}$ or less. It is the right angle of the projection and the right angle, put to a use that essay did not need, and it is why the residual norm is a poor detector of a wrong solution and a good estimate of the noise.

The complement’s $s^2$ has no such protection. With the heavy row last it is wrong by 5.8·10⁻⁸, which is twice the 2.9·10⁻⁸ error in the norm of $Q_2^{\mathsf T} b$. The unsorted complement perturbs the whole residual vector’s scale, not its direction, and a sum of squares feels every bit of a scale.

Why an unsorted Cook’s distance is better than its parts #

That observation has a curious consequence for the unsorted complement, and it is the one place where computing the pieces badly gives a whole that is better than the pieces.

Cook’s distance is a ratio of squared residuals: $d_i^2$ over $s^2$. If both are computed from a residual vector whose scale is wrong by a factor $1 + \epsilon$, both acquire a factor $(1+\epsilon)^2$ and the ratio does not. On the straight line with the heavy row last, the deleted residual is wrong by 2.9·10⁻⁸, $s^2$ by 5.8·10⁻⁸, and Cook’s distance by 9.3·10⁻¹⁰ — thirty times better than either. Most of the unsorted complement’s error on this design is a common scale, and the diagnostic divides it out.

Not all of it, and not on every design. On the cubic with the heavy row last the same four numbers at $4^{30}$ are 1.9·10⁻⁷, 9.5·10⁻⁸, 4.8·10⁻⁸ and 2.9·10⁻⁷: the heavy row’s own residual carries an error of its own beyond the scale, and Cook’s distance is no better than the worst of its parts. So the cancellation is an observation about one design rather than a guarantee, and the practical rule does not change: sort, and the question never arises.

Cook’s distance at 3.4·10¹⁷ #

Put together, the diagnostic a code reports is either right to the last bit or wrong in every digit, depending on two choices that no textbook formula distinguishes.

Cook’s distance for the heavy row grows as the weight does, from 0.037 at weight one to 3.4·10¹⁷ at $4^{30}$ — an observation carrying that much weight is by construction the most influential in the fit, and the diagnostic says so. With the row first, the complement’s value keeps fourteen or fifteen digits across the whole range. The textbook value keeps fourteen at weight one, eleven at $4^6$, eight at $4^{12}$, four at $4^{18}$, one at $4^{24}$ and none from $4^{27}$. It loses one digit for every factor of ten in the weight, which is the rate of the residual and the rate of the divisor, because the two losses are the same loss.

What to compute #

The rule the weighted-divisor essay stated needs one more clause, and the clause is about where the residual comes from.

Sort the rows by decreasing weight. Factorise by Householder. Take 1 − h and the residual from the complementary block — $\lVert Q_2(i,:)\rVert^2$ and $Q_2(i,:)\,Q_2^{\mathsf T} b$ — and divide the one by the other for the deleted residual. For a code that does not form $Q$, both are available by applying the stored reflectors, which the factor a sparse code keeps anyway showed costs one pass: the row of $Q$ from the reflectors applied to a unit vector, and $Q^{\mathsf T} b$ from the same reflectors applied to the data, which a least-squares solve has already computed. Take $s^2$ from either route; with the rows sorted both are right, and unsorted the textbook one is the safer of the two.

Do not form the heavy row’s residual from the solution. It is the step every statistics text writes first, and it is the one that cannot be repaired afterwards: sorting does not help, refinement does not help, and exact arithmetic does not help, because the information is not in $x$.

If a factor is not available at all — a fit produced by some other method and handed on as a solution — refit without the row. It costs a factorisation per observation examined and is right in every order at every weight, because the deleted residual is well conditioned and only the shortcut to it is not.

What this does not settle #

One or two weighted rows, weights that are powers of four so that the scaled rows are exact, noise that is Gaussian and modest. A heavy observation whose own datum is badly wrong — the case Cook’s distance exists to find — has a larger deleted residual, which helps the textbook route by shrinking $|b_i|/|r_i|$, and has not been measured.

Householder QR only, as before. The complement’s accuracy with the heavy row first is argued, in the weighted-divisor essay, from the first reflector pivoting on it; nothing here isolates the reflector’s own rounding, and a Givens or modified Gram–Schmidt route would have its own ordering rule.

The cancellation in the unsorted Cook’s distance is measured on two designs and explained on one. Why the cubic’s heavy row carries an error beyond the common scale, when the line’s does not, is not explained by anything here.

Still open: many weights at once, and a leverage the complement cannot see #

A thousand rows over ten decades. Still the realistic case, and still unmeasured. The numerator adds a reason to expect the partial sorts a streaming code produces to cost more than the divisor alone suggested: every row factorised after a heavier one sees its residual scale perturbed, so a diagnostic over all rows reads a mixture of scales. Whether the common-scale cancellation survives a mixture, or breaks as it did on the cubic, is the direct question.

The masking pair under weights. Two observations that hide each other found that a pair of influential points needs the two-by-two block of the hat matrix, not two diagonal entries. The complement gives that block as two rows of $Q_2$ and their inner product, which is not a sum of squares and can cancel. Whether the off-diagonal entry keeps its digits with both heavy rows first, as the diagonal ones do, is the measurement that would carry this rule from single deletions to joint ones.