One minus a leverage is a subtraction

The deletion diagnostics of a regression all divide by the same number. The deleted residual is $e_i/(1 - h_i)$. The studentised residual divides by the square root of $1 - h_i$. Cook’s distance multiplies by $h_i/(1 - h_i)^2$. The leave-one-out cross-validation score sums the squares of $e_i/(1 - h_i)$. And removing an observation from a factorisation takes the square root of $1 - h_i$ inside a hyperbolic rotation. Every one of them is most informative exactly where 1 − h is small — at the observations that hold most of a unit of influence — and every one of them is fed a 1 − h that was computed by working out h and subtracting it from one.

That subtraction is the whole subject of this essay. A leverage near one is a number near one, and a number near one stored in double precision carries its error in the last of its sixteen digits. Subtract it from one and the leading digits cancel, the error stays where it was, and what remains is a small number whose relative error is the leverage’s absolute error divided by 1 − h. Nothing is wrong with the leverage. The loss happens afterwards, in the one line of arithmetic nobody thinks of as a computation.

The essay on the influence budget closed by recommending the squared row norms of Q from a QR factorisation as the better-behaved route to the leverages. For the leverages that is right. For one minus them it is exactly as bad as any other route that subtracts, and the repair is a different sum of squares from the same factorisation.

A leverage to the last bit, and one minus it to four digits #

The measurement needs a row whose leverage is known and a reference that does not round. A forty-by-six design is built with its first row scaled until its leverage is 1 − 10⁻ᵏ, for k from 2 to 15.5, and at each k the exact value of 1 − h is computed in rational arithmetic from the stored doubles, so the reference is the true 1 − h of the matrix actually held and not of the matrix intended.

The figure at the head of the page is the relative error of three computed versions of 1 − h against that reference. The first takes the leverage as $\lVert R^{-\mathsf T} a_i \rVert^2$ from the Cholesky factor of AᵀA and subtracts it from one. The second takes it as the squared norm of row i of the thin QR factor $Q_1$ and subtracts it from one. At k = 12 they are wrong by 1.11·10⁻⁴ and 4.44·10⁻⁴ of the true value. Twelve of sixteen digits are gone.

The leverages themselves were not wrong. A relative error of 1.11·10⁻⁴ in a quantity of size 10⁻¹² is an absolute error of 1.1·10⁻¹⁶, which is one unit of roundoff: the Cholesky route computed h to the last bit a double has, and the subtraction turned that last bit into the fourth significant digit of the answer. The dotted line in the figure is 10ᵏ times the unit roundoff, and both subtracting routes lie along it at every k: the Cholesky route within a factor of four of the line from k = 2 to 15, the thin QR route within a factor of twenty. The thin route’s extra factor is its leverage carrying a few more units of roundoff than the Cholesky route’s does, which a subtraction then magnifies by the same 10ᵏ. That is the law for any route that subtracts, whatever produced the leverage: a relative error of a small multiple of u/(1 − h), where the multiple is how many units of roundoff the leverage itself carries.

The third route has no such slope. It is wrong by at most 3.03·10⁻¹⁵ of the true value at every k, including the ones where the first two have nothing left.

At the last rounding step, a zero and a negative number #

Plotted as values rather than errors, the same computation shows what the diagnostics are handed. Up to k = 14 the two subtracting routes wander around the true value by the amount the law allows. At k = 15.5, where the true 1 − h is 3.33·10⁻¹⁶, the Cholesky route returns exactly zero and the thin QR route returns −3.11·10⁻¹⁵.

Those are not small errors in the usual sense; they are values of the wrong kind. A deleted residual divided by zero is infinite. A studentised residual needs the square root of a negative number and gets NaN. Cook’s distance squares the denominator and so reports a finite, enormous and positive number for an observation whose true distance is merely large. A hyperbolic downdate asked to remove that row either divides by zero or takes the square root of a negative discriminant and fails. Four diagnostics, four different kinds of nonsense, from one subtraction that returned a number with the wrong sign.

The complementary route cannot do that, whatever its accuracy, because it is a sum of squares. It can be wrong but it cannot be zero unless every term underflows, and it cannot be negative at all. On this design at k = 15.5 it returns 3.33·10⁻¹⁶, wrong by 1.5·10⁻¹⁶ of itself — a number that sits a third of a unit of roundoff above zero, returned with every digit a double can give it.

One minus a sum of squares, or the sum of the other squares #

The route with no slope comes from the same factorisation as the second one. A full QR factorisation of the m × p design gives an m × m orthogonal Q, whose first p columns $Q_1$ span the column space and whose remaining m − p columns $Q_2$ span its orthogonal complement. Every row of an orthogonal matrix has unit length, so for row i

$$h_i = \lVert (Q_1)_{i} \rVert^2, \qquad 1 - h_i = \lVert (Q_2)_{i} \rVert^2 .$$

The second identity is the first rearranged, and in exact arithmetic the two are the same statement. In floating point they are different computations. Geometrically they are the two legs of the right angle a projection makes: the row’s unit vector split into its component in the column space and its component orthogonal to it, and 1 − h is the squared length of the second. The first route gets 1 − h by computing p squares, adding them to a number near one, and subtracting from one. The second computes m − p squares whose sum is the small number, and never forms anything near one. This is two routes to one number in its sharpest form: identical algebra, and arithmetic that differs by twelve digits.

The law does not care about the shape of the design. On sixteen observations of two columns the subtracting routes are wrong at k = 12 by 1.11·10⁻⁴ and 5.55·10⁻⁴, the same order as on forty of six, and the complement is wrong by at most 4.24·10⁻¹⁶ at any k. The subtraction’s cost is set by 1 − h and by nothing else about the matrix.

Nor is the complement expensive. A Householder QR stores the p reflectors that produce R, and row i of Q is obtained by applying those p reflectors to the i-th unit vector: the first p entries of the result are the row of $Q_1$ and the other m − p are the row of $Q_2$. So the complement’s row norm costs the same O(mp) as the thin factor’s row norm — the same reflectors applied to the same vector, with a different part of the result summed. The full m × m matrix never needs to be formed, and on a design whose rows number in the millions it must not be.

The complement is not free: it pays κ(A)·u #

The forty-by-six design was made to have a large leverage by scaling one row, and nothing else about it is difficult. That is the right control and the wrong model of how large leverages arise. In practice a leverage near one belongs to an observation far from the others along some direction of the design — a setting well outside the range of the rest — and an observation like that makes the design itself ill conditioned.

So the second experiment is a polynomial fit: thirty points scattered on [−1, 1], one more at x, and columns 1, x, …, x^d for degree d from one to three. Moving x from 2 out to 1,000 drives the far point’s 1 − h down and κ(A) up together, and the exact 1 − h is again computed in rationals from the stored design.

For a straight line nothing is at stake. At x = 1,000 the far point’s 1 − h is 9.63·10⁻⁶ and κ(A) is 183. The subtracting routes lose five digits, as the law says, and are wrong by 3.15·10⁻¹¹ and 6.07·10⁻¹¹. The complement is wrong by 1.18·10⁻¹⁴ — not the 10⁻¹⁶ it achieved on the scaled design, but eleven digits, which no diagnostic needs.

The complement’s error has grown, though, and the figure shows what it has grown with: a second dotted line at κ(A)·u, which the complement follows as the subtracting routes follow u/(1 − h). That is the reading of a backward-stable factorisation. The Householder QR computes the exact factors of a matrix within about u·‖A‖ of the one stored, and a row of $Q_2$ moves under such a perturbation by roughly κ(A) times as much, relatively. So the complement’s error scales with the conditioning of the design, as the subtraction’s scales with the smallness of 1 − h, and on the scaled design κ(A) was small and only the second mattered.

For a quadratic both laws start to bite. At x = 1,000 the far point’s 1 − h is 2.41·10⁻¹² and κ(A) is 3.23·10⁵. The subtracting routes are wrong by 1.17·10⁻⁴ and 2.56·10⁻⁴ — four digits left — and the complement by 4.38·10⁻¹¹, ten digits left. Both are near their own laws: u/(1 − h) = 4.61·10⁻⁵ for the first two and κ(A)·u = 3.59·10⁻¹¹ for the third.

Across the whole sweep of three degrees and seven positions no route is ever wrong by more than ten times its law. Over the points where each law is above rounding and below total loss, the Cholesky route sits at a median of 1.84 times u/(1 − h), ranging from a twelfth to 2.7; the thin QR route at a median of 4.08, from 1.3 to 7.7; and the complement at a median of 1.22 times κ(A)·u, from a seventh to 4.2. A rounding can land well, and none lands badly. The laws are envelopes that the measurement fills rather than exact predictions, and the distinction matters for what can be promised in advance: the digits a route keeps can be bounded before it runs and read off exactly only afterwards.

The cubic is where the two routes separate for good. At x = 100 the far point’s 1 − h is 4.2·10⁻¹³, the subtracting routes are wrong by 5.82·10⁻⁴ and 10⁻³, and the complement by 7.15·10⁻¹¹ against κ(A)·u = 7.67·10⁻¹¹ — on its law to within a tenth. At x = 300 the true 1 − h is 5.79·10⁻¹⁶, about five units of roundoff, and both subtracting routes return numbers with no correct digit, wrong by 54 and 81 per cent. The complement is wrong by 2.95·10⁻¹⁰: nine digits. At x = 1,000, where 1 − h is 4.2·10⁻¹⁹ and below the unit roundoff altogether, the complement is still wrong by only 2.66·10⁻⁷, six digits of a number no subtraction in double precision can represent at all.

Why half: the far point’s 1 − h is one over κ(A)² #

The two laws would say nothing useful about each other if 1 − h and κ(A) were unrelated. On these designs they are not, and the relation is a square.

From x = 30 outwards, for every degree, 1/(1 − h) lies between 3.12 and 5.06 times κ(A)². The reason is the geometry of an isolated row. The far row is long — its entries grow as x^d — and it points almost along the design’s top singular direction, so it supplies most of σ_max. The thirty ordinary rows supply σ_min. One minus the far row’s leverage is one over one plus that row’s squared length measured against the rest of the design, and that is roughly (σ_min/σ_max)², the inverse square of the condition number, up to a factor set by how the ordinary points are spread.

So on a design whose large leverage comes from an isolated observation, a route that pays u/(1 − h) pays about κ(A)²·u, and a route that pays κ(A)·u pays the square root of that. In digits, the subtraction loses about twice what the complement loses. That is the same exponent the road that squares the problem finds in the normal equations, arriving from a different direction: there the square was forming AᵀA, and here it is subtracting a number whose distance from one is already a squared quantity. The complement does not remove the conditioning of the design — nothing that computes a property of A can — but it declines to square it.

Which divisor to compute #

The practical rule follows directly. Compute 1 − h as the squared row norm of $Q_2$, not as one minus the squared row norm of $Q_1$ or of $R^{-\mathsf T} a_i$. It costs the same reflector applications, it can never return zero or a negative number, and it keeps every digit on a well-conditioned design and half of those the subtraction would lose on an ill-conditioned one. The leverage itself can still come from either route, because h is not the number the diagnostics are sensitive to.

Know how a relative error in 1 − h reaches each diagnostic. The deleted residual inherits it one for one, the studentised residual at half the relative size through the square root, and Cook’s distance at twice the relative size through the square. At h = 1 − 10⁻¹², a Cook’s distance computed from a subtracted 1 − h is wrong in its fourth digit from the second power alone. None of this shows in the ordinary residual, which is tiny at such a row for the reason a small residual is not a small error gives: the observation has bought its own fit, and every digit of evidence against it sits in the divisor.

Read the downdate the same way. A hyperbolic rotation that removes a row needs $\sqrt{1 - h}$ and forms it by subtraction inside the algorithm; its residual in the observation that cannot be removed grew with a fitted slope of 1.022 against 1/(1 − h), which is this essay’s u/(1 − h) law measured on a different computation. Factorising afresh there, and reading the complement here, are the same escape: compute the small quantity as itself.

And for a pair, read the block the same way. The masked pair of two observations that hide each other needs 1 − s, where s is the larger eigenvalue of its two-by-two leverage block. Since QQᵀ = I, the complement $I - H_{SS}$ equals $(Q_2)_S (Q_2)_S^{\mathsf T}$, so the block’s complement is a two-by-two Gram matrix of the pair’s rows of $Q_2$, and its smallest eigenvalue is 1 − s computed without subtracting s from one. A pair near a joint share of one is exactly the case in which the subtraction would matter.

The general shape is one this field keeps meeting. A quantity defined as a difference of two nearly equal numbers is catastrophically cancelled if it is computed as that difference, however accurate the two numbers are — cancellation takes the answer, not a digit, because the subtraction is itself exact and only exposes the error its operands already carried — and the escape is almost never more precision. It is an identity that expresses the small quantity directly. An answer that is known is what makes the loss visible here: without the rational reference, the subtracting routes’ 1 − h at k = 12 would look like any other number with sixteen digits printed.

Still open: the divisor without Q, and leverage under weights #

The complement needs Q, or at least the reflectors that define it. A sparse QR factorisation is usually computed Q-less — R is kept and the orthogonal factor is discarded, because storing it would destroy the sparsity that made the factorisation affordable — and without Q the only route to 1 − h left is the subtraction. Whether a correction step on the seminormal equations, or a solve against R that computes the complement’s norm some other way, can recover the lost digits without keeping Q is the measurement this essay does not make.

The other direction is weights. A weighted least-squares fit has leverages that are the diagonal of a different projection, and a weight driven towards infinity drives a row towards h = 1 — the limit in which a heavily weighted observation becomes a constraint. There the smallness of 1 − h is chosen rather than found, and the question is whether the complement’s κ(A)·u law survives weights spread over many decades, or whether the weighted design’s conditioning makes the complement pay the same square the subtraction does.