The factor a sparse code keeps anyway

One minus a leverage is a subtraction established the defect and the repair together. Every deletion diagnostic divides by 1 − h; computing that as one minus a computed leverage loses digits in proportion to how close the leverage is to one, however accurate the leverage; and the squared norm of the complementary block of a full orthogonal factor is the same number with no subtraction in it, keeping every digit on a well-conditioned design.

It also named the reason the repair might not be available. A sparse QR factorisation is usually described as Q-less: the orthogonal factor is discarded because storing it would destroy the sparsity that made the factorisation affordable, and what is kept is the triangular factor alone. Without the orthogonal factor, the argument ran, the only route back to 1 − h is the subtraction.

Three ways of avoiding the subtraction were available to try. All three fail, and they fail at the same number, which is what makes the failure informative rather than a series of disappointments.

The fourth route does not avoid the subtraction. It goes and gets the orthogonal factor, from something a sparse code keeps whether it wants to or not.

Three repairs that all land on the same line #

The first is a correction step. The leverage comes from a triangular solve, that solve has a residual, and one step of iterative refinement on it is the standard repair for a factorisation routed through the normal equations. It gives a better leverage.

It gives exactly the same divisor. At a leverage of 1 − 10⁻¹² the plain subtraction is wrong by 1.11·10⁻⁴ of itself and the refined one by 1.11·10⁻⁴.

The second is to compensate the sum of squares. The leverage is a sum of squared entries, so it can be accumulated as a leading term and an error term rather than rounded into one double at every step; then one minus the leading term is exact wherever the cancellation is severe, and the digits that would have been lost are in the error term. That is the standard construction and it is the right construction for a different problem — it is what the sum that cannot be wrong is built on, where the quantity being accumulated is the answer rather than something the answer is subtracted from. Here it returns 4.23·10⁻⁵ where the plain route returns 1.11·10⁻⁴ — a factor of two and a half, at a leverage where four decades have already gone.

The third is the one that ought to work. The correction from the refinement step is a small number added to a large one, and adding it rounds it away; keeping the pair separately, and computing the leverage as the leading sum plus twice the cross term, should give a leverage accurate past a unit of roundoff and therefore a divisor to match. With the residual formed in extended precision so that the correction means something, it returns 1.46·10⁻⁵.

Better by a factor of eight, and eight is not the eleven decades that have been lost.

Why they land there #

The dashed line on the figure is a unit of roundoff divided by 1 − h, and all three routes sit on it.

That is the whole explanation, and it is a statement about arithmetic rather than about any of the three constructions. The leverage is a number near one. Whatever route produced it, it is stored as a double, and a double near one carries an absolute error of about a unit of roundoff — that is what the representation is. Subtracting it from one produces a number of size 1 − h carrying that same absolute error, so the relative error of the result is a unit of roundoff over 1 − h.

Improving the leverage’s relative accuracy does nothing, because the leverage’s relative accuracy was never the problem: at h = 1 − 10⁻¹² a leverage correct to fifteen significant figures is still a double, and its last bit is still worth 10⁴ times the answer.

So the three repairs are three ways of computing h more carefully, and h being computed carefully is not the question. The vector that hides it makes the same distinction for a sum: a better summation is not a better answer when the answer is a difference of things the summation cannot see. The question is whether the answer can be obtained without ever forming h, and none of the three does that.

There is a fourth thing worth noticing on the figure, and it is about the family rather than about the routes. The designs here are built to carry a row of a requested leverage, and a row cannot have leverage near one without the design being ill conditioned — κ runs from 9.7 at 1 − 10⁻² to 5.3·10⁷ at the far end, roughly as the reciprocal of the square root of the divisor. The two difficulties arrive together and cannot be separated by construction, which is worth stating because it means no experiment on this family can attribute the loss to the cancellation alone. What settles the attribution is that the measured loss tracks the subtraction’s law rather than anything containing κ.

What a sparse factorisation actually keeps #

The premise to attack is not the arithmetic but the storage claim, and the storage claim is imprecise in a way that matters.

A sparse QR does not store the orthogonal factor as a matrix. That is true and it is the reason the factorisation is affordable: the orthogonal factor of a sparse matrix is dense, and materialising it would cost more than the problem.

It does store the Householder vectors. It has no choice, and an ordering that does not wait for the numbers is where prices what a sparse factorisation keeps and why. Solving a least-squares problem means applying the transpose of the orthogonal factor to the right-hand side, and the only way to do that without the matrix is to apply the reflectors one after another. A factorisation that threw its reflectors away could not solve anything; it would be a triangular matrix with no connection to the problem it came from.

So the orthogonal factor is unavailable as a matrix and available as an operator, and the complement needs only the operator.

One row of the factor is the reflectors applied, in order, to a unit vector. Each reflector is symmetric, the product of all of them is the factor, and applying them to a single basis vector costs one pass — four operations per stored entry of the factorisation, 900 of them on the design measured here, against a stored factor of 240 numbers that was never formed.

It returns the same bits #

Rebuilt from the reflectors and summed past the pth entry, the complement agrees with the one computed from a stored orthogonal factor to within 10⁻¹⁸ at every leverage measured — which is to say, bit for bit.

It has to. It is the same arithmetic in a different order: the stored factor is the reflectors applied to all m unit vectors and this is the reflectors applied to one of them, and nothing in between rounds differently.

That is the answer to the question, and the shape of it is worth keeping. The open question asked whether some other computation could recover the digits the subtraction loses. It cannot, and three attempts say why. What recovers them is the computation that was already known to work, reached by noticing that the thing it needs was never actually thrown away.

Its own limit, on a family that has one #

The reflector route is the complement, so it inherits whatever the complement’s limit is — and the leverage family above cannot show what that is, because on it the complement is right to fifteen figures everywhere.

The other family in use here puts the difficulty in the design instead. Thirty ordinary points and one far outside them, fitted by a cubic: as the far point moves out, the condition number rises from 7.2 to 6.9·10⁸ and the divisor falls from 6.4·10⁻³ to 4.2·10⁻¹⁹, without anybody asking for a leverage.

Here the reflector route is not exact. At a far point of 100 it is wrong by 7.1·10⁻¹¹, and at 1000 by 2.7·10⁻⁷. Those track the dashed line, which is the condition number times a unit of roundoff — the loss an orthogonal factorisation of an ill-conditioned matrix makes and cannot avoid.

What it is not is the subtraction’s loss. At a far point of 100 the subtraction is wrong by 5.8·10⁻⁴, at 300 by 0.53, and at 1000 it has no correct digit whatever. The gap between the two at the far end is about nine decades.

So the reflector route has a real limit and it is the right limit. A route to 1 − h cannot be better than the factorisation it comes from, and the complementary-block construction is exactly as good as the factorisation — no better, and not degraded by anything it does afterwards.

The two families have to be read together, and neither is sufficient alone. The first sets the divisor directly and drags the conditioning along with it; the second sets the conditioning and drags the divisor along. In both, the subtracting routes follow a unit of roundoff over the divisor and the reflector route follows the condition number times a unit of roundoff, and those two laws are what the pair of experiments establishes that neither could establish by itself.

What it costs, and when it is worth it #

The cost is one pass over the reflectors per row wanted. Counted, it is 316, 616, 900 and 1168 operations at two, four, six and eight columns on a forty-row design, and 516, 900 and 1284 at twenty-four, forty and fifty-six rows with six columns — linear in both, and a few per cent under 4mp because the last reflector is shorter than the first.

For a handful of diagnostics that is nothing beside a factorisation costing some multiple of mp². For all m rows it is 4m²p, which is what materialising the whole orthogonal factor costs — so the route is cheap exactly when the diagnostics are wanted for a few observations and stops being cheap when they are wanted for every one.

That is a more useful boundary than it looks, because deletion diagnostics are usually wanted for a few. The observations whose influence is in question are the ones some cheaper screen has already flagged, and two observations that hide each other is the case where the screen has to look at pairs rather than singles — a quadratic number of candidates, but only a handful of them worth an exact divisor.

Where every row’s divisor is wanted, the honest answer is the one the subtraction gives, with the caveat that it is reliable only while the leverages stay away from one. And the leverages staying away from one is not a coincidence that can be relied on: influence is decided before the data shows that a fit has exactly p units of influence to distribute among m observations, so a design that concentrates them is a design with a few rows near one and the rest near zero — which is precisely the case where the few matter.

What the two laws say to a practitioner #

The pair of experiments leaves a rule that fits in two sentences, and it is worth extracting because the temptation is to read the first family’s result as the complement is exact and the second’s as nothing is exact.

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. So the question of which to use is the question of which of those two numbers is larger, and on a design where a row is influential because it is far from the others they are not the same number at all: at the cubic’s far point of 100 the divisor is 4.2·10⁻¹³ and the condition number is 6.9·10⁵, so the subtraction owes twelve decades and the factorisation six.

That gap is not an accident of these families. A leverage approaches one when a row leaves the space the others span, and leaving that space raises the condition number as a square root rather than proportionally — which is why the two laws diverge in the direction they do, and why the route through the factor is the better one by roughly a factor of two in decades wherever the question arises at all.

The exception is a design that is ill conditioned for a reason that has nothing to do with the row in question. Scaling the columns badly does it — the scaling that buys ten orders is where a change of variable that hands the whole of such a loss back is measured, and a condition number scaling cannot move is the componentwise measure that was never inflated by it in the first place. There the factorisation’s loss is an artefact while the divisor is ordinary, and the subtraction is the better route. Equilibrating first removes the case rather than deciding it, which is the honest advice: a design whose condition number is a statement about units should be rescaled before any diagnostic is computed from it at all.

What must fail for any of this to be wrong #

Five claims and a refusal. That a correction step on the triangular solve buys under two decades in the divisor, at every leverage from 1 − 10⁻⁸ to 1 − 10⁻¹⁴. That a compensated sum of squares buys under two decades. That the same correction kept as a pair, with the residual formed in extended precision, buys under two decades. That the plain subtraction’s error sits within a factor of five of a unit of roundoff over the divisor, which is what says the loss is the cancellation rather than anything around it. And that the reflector route reproduces the stored factor’s complement to within 10⁻¹⁸ and is right to within 10⁻¹³ at every leverage.

The refusal is fed the claim that the complement needs a stored orthogonal factor, and required to fail.

The three negative claims are stated as bounds on the improvement rather than as equalities, because two of them do buy something — a factor of two and a half, a factor of eight — and a claim that they buy nothing would be false. What they must not do is buy a decade, and that is what is required of them.

What this does not settle #

Dense designs throughout. The reflector count and the cost per reflector are what a dense Householder factorisation produces, and in a sparse factorisation both depend on the pattern: the reflectors are sparse, there are more of them, and the cost of a pass is the number of their stored entries rather than 4mp. The ranking cannot change, because it does not depend on the cost, but the boundary in the previous section does.

One row at a time. Every measurement asks for the first row’s divisor. Asking for several shares no work between them in this formulation — each is its own pass — and whether a blocked application of the reflectors to several unit vectors at once recovers the usual factor is unmeasured.

And the reflectors are assumed to be the ones the factorisation produced. A code that stores them in a compressed form, or that has applied a column permutation for fill, has more bookkeeping between the stored data and the row of the orthogonal factor than this measurement models.

Still open: the sparse cost, and the other diagnostics that divide #

The pass in a sparse factorisation. The whole practical question is what one pass over the reflectors costs when they are sparse, and that is a property of the pattern rather than of the arithmetic. A measurement on a genuinely sparse design — counting the stored entries the pass touches, against the factorisation’s own entry count — would turn the boundary above from a dense estimate into a rule. It would also say whether the answer depends on the fill-reducing ordering, which would make the diagnostic’s cost a fourth thing that ordering decides.

Several rows at once. A pass per row is the formulation measured here, and the reflectors are applied to the unit vectors one at a time. Applying them to a block of unit vectors is the same operation with a matrix in place of a vector and should cost substantially less per row on a machine, for the reason every blocked kernel does. How many rows it takes before the blocked form wins is the measurement, and it is the same question a second objective that is the first one doubled asks about a hierarchical product — where the answer turned out to be that the operation count cannot see the difference at all.

The other quantities that divide by a small number. The divisor is the worst case because it approaches zero, but a deletion diagnostic also divides by an estimate of the residual variance and multiplies by a ratio of leverages. Whether any of those carries its own cancellation, and whether the reflector route supplies them too, is a question this essay has not asked and the pair-deletion case makes sharper, since a pair’s divisor is a two-by-two determinant rather than a scalar.

And the weighted case. A weighted fit has leverages that are the diagonal of a different projection, and a weight driven towards infinity drives a row towards a leverage of one — so the smallness of the divisor is chosen rather than found. Whether the reflector route survives weights spread over many decades, where the reflectors themselves are computed from a badly scaled matrix, is the direct continuation, and it is the question one minus a leverage is a subtraction left open in the same breath as this one.