The licence is not the boundary

A library that chooses its algorithm at run time chooses it with a predicate. Something about the input is measured or estimated, a branch is taken, and the caller is told nothing about which side of the branch the answer came from, because the routine has one name and one signature either way. The length that changes the kernel measures one such branch on a dot product, where the quantity in the predicate is the length of the vector and the two implementations differ only in how they round.

A predicate like that is worth having when two things hold. It has to be computable before the work is done, since a test that needs the answer is not an admission test. And it has to be a boundary: the behaviour on one side has to differ from the behaviour on the other, at the place the predicate names rather than somewhere else.

Those two requirements are usually collapsed into one. A predicate that is computable is treated as a predicate that divides, because the arithmetic that produced it was a derivation about failure and a derivation about failure is assumed to describe where failure begins. That assumption is what gets measured here, and it is separable from the derivation being right: a bound can be correct, tight and well calibrated over twelve decades and still name a location nothing happens at.

This essay takes two predicates from the same corner of least squares, both computable in advance, both stated in the unit roundoff, and both quoted in the same tone of voice — and measures where the behaviour they are supposed to divide actually changes. One of them is a boundary. The other is a calibrated accuracy model with no edge in it at all, and the difference between the two is where the dangerous answers live.

A boundary that is where it was said to be #

The normal equations are the case where the prediction and the failure meet. Forming AᵀA squares the condition number, and below ε ≈ √u the quantity 1 + ε² rounds to 1, every entry of the cross-product matrix becomes exactly one, and the matrix a textbook would now ask for an inverse of is exactly singular. That is the road that squares the problem’s subject and its measurement; what matters here is only the shape of the agreement.

In binary32 the route stops at ε = 1.8·10⁻⁴ and √u is 2.4·10⁻⁴. The measured breakdown is 0.75 of the predicted one, and the sweep steps ε by a quarter of a decade, so those two values are one grid step apart. Nine precisions from 16 to 40 significand bits give the same reading: the failure lands at 0.59 to 0.92 of √u at every one of them, while the ε at which it happens moves from 3.2·10⁻³ to 5.6·10⁻⁷. Across twenty-four bits of significand, a factor of 4,096 in √u, the ratio does not leave a factor of 1.6.

That is what a computable boundary looks like when it is one. A caller who wants the coefficients and has a value of ε can decide in advance whether to take this route, and the decision is right at every precision, because the number in the prediction and the number in the failure are the same number.

The important part is what the prediction is about. It is not about accuracy: the normal-equations route is already far behind the QR route long before it breaks, by more than a factor of ten in the worst error across the surviving range. √u predicts the point at which the method stops producing anything, which is a discrete event with a computable location. Nothing about the digits has to be interpreted to see it, and no threshold has to be chosen to detect it.

That is why the agreement can be as sharp as it is. The event being predicted is an equality between two floating-point numbers — 1 + ε² and 1 — and it happens at the ε where ε² falls below the spacing of the numbers near one. There is no tolerance anywhere in that sentence, no choice of what counts as having failed, and no dependence on the shape of the matrix beyond the ε in it. A prediction of a discrete event can be checked by asking whether the event happened, and the only residual freedom is where the sweep put its sample points.

The other predicate, and the number that is missing from it #

The single-reduction QR is the least-squares route taken for a completely different reason. Forming G = AᵀA is a sum over the processors’ slices, so it is one reduction and everything after it is local — the fewest communication rounds any factorisation can use, which is what the message and the word counts and why the method is on the menu at all. Its licence is stated everywhere in the same form: it is usable when κ²u ≪ 1.

Read as an equality, that names a condition number. The unit roundoff of binary64 is u = 1.11·10⁻¹⁶, so κ²u = 1 at κ = 9.49·10⁷.

Before asking what happens there, it is worth checking whether the quantity in the predicate predicts anything, and it does. On 256×8 matrices built with a prescribed condition number, ‖QᵀQ − I‖ of the implied Q = AR⁻¹ — the factor a caller actually meets, since a solve through R is exactly this Q applied to a right-hand side — sits at a fixed fraction of κ²u across the whole usable range. The sweep steps κ by half a decade from 10² to 10⁸; six of its thirteen readings are these:

Over thirteen condition numbers from 10² to 10⁸ the ratio stays between 0.12 and 0.88 — a spread of 7.2, while the quantity being predicted moves through twelve decades. As a model of how much orthogonality is lost, κ²u is excellent, and that is a stronger result than the licence needs.

It is also the whole of what the licence says. The derivation behind κ²u ≪ 1 bounds the loss of orthogonality of the implied Q; it is a statement about how many digits survive, and the measurement above is that it is a good one. It contains no statement about the Cholesky step completing, and no statement about anything discrete at all — there is no equality between two floats in it, nothing that either happens or does not, and therefore nothing for a boundary to be the location of.

κ²u is a prediction of the error, and ≪ is not a number. A caller who wants the implied Q orthogonal to 10⁻⁸ is out of licence somewhere between κ = 10⁴ and κ = 3.2·10⁴, three and a half to four decades below 9.49·10⁷; a caller who will accept 10⁻⁵ has until about 10⁶. The equality point is where the model predicts a loss of one, which is total — orthogonal is a number and one is the value it takes when the columns have stopped being a basis. So κ = 9.49·10⁷ is not the edge of the method’s usefulness. It is somewhere well past it, and how far past depends on a tolerance the licence does not mention.

Past the licence, the method goes on answering #

At κ = 10⁸, just past the equality point, the method returns. What it returns has ‖QᵀQ − I‖ = 0.37, so the answer is not a factorisation in any sense a caller could use, and it comes back with no diagnostic attached: a triangular R of the right shape, a residual as good as the other methods’, and a basis that is not one.

Two and a half decades further on, at 3·10¹⁰, the Cholesky step meets a pivot of −6.5·10⁻¹⁷ and stops. The Gram matrix has ceased to be numerically positive definite, which is the same verdict a definiteness test returns and the same mechanism, arriving here as the outcome of a network decision rather than as a question about a matrix. There is nothing to hand back, and that is the safe failure: a caller who gets it knows.

And at 10¹², a decade and a half further still, the method returns again — a factor whose implied Q has ‖QᵀQ − I‖ = 1.49, which is further from orthogonal than the one at 10⁸ and further than the one at 10¹⁰. So the method’s behaviour past its licence is not degradation. It is a mixture of two outcomes, one safe and one not, and neither of them is what a licence stated as a single inequality would lead a reader to expect. The column sweep and the reduction tree return a usable factor of every one of these four matrices, with ‖QᵀQ − I‖ between 8.5·10⁻⁹ and 1.4·10⁻⁴ — so nothing about the matrices themselves prevents a factorisation, which is what a reduction tree is for.

Where it declines is a fact about the list it was asked #

Move the sample. The same experiment at 10⁸, 2·10¹⁰, 4·10¹⁰ and 10¹² reports the first refusal at 2·10¹⁰, on a pivot of −2.7·10⁻¹⁶, and reports a second at 4·10¹⁰. The previous reading — one refusal, at 3·10¹⁰, on a pivot of −6.5·10⁻¹⁷ — is nowhere in it, and neither is any of the numbers that reading was quoted with.

Nothing has changed about the method, the matrices are built the same way, and the two lists overlap at both ends. What changed is which condition numbers were asked. A report of the form this method declines at 3·10¹⁰ is therefore a statement about a list of four numbers, and a caller reading it as a property of the method has been handed the wrong kind of fact. The pivot magnitude moves with it: −6.5·10⁻¹⁷ and −2.7·10⁻¹⁶ are both nonpositive and both tiny, and neither is close enough to the other to be a measurement of anything.

The habit this site runs on says to give a claim a test it could fail, and the test here is simply to sample more finely. Six condition numbers to a decade — one, 1.5, 2.2, 3.2, 4.6 and 6.8 times each power of ten, from 10⁷ to 10¹³ — is thirty-seven matrices at 256×8. The method returns a factorisation at nineteen of them and refuses at eighteen. The first refusal is at 2.2·10⁸, on a pivot of −1.5·10⁻¹⁶ met at the seventh column of eight: a factor of 2.3 past the licence rather than two decades past it, and nearly two orders of magnitude below the earliest refusal either four-point list reported.

The set it refuses is not an interval #

That fine sweep says something the four-point lists can only hint at. Of the nineteen condition numbers at which the method returns a factorisation, eleven are above the first refusal, and they run from 3.2·10⁸ to 10¹³. The last of them is five decades past the licence and four and a half past the first refusal, and the factor it hands back has ‖QᵀQ − I‖ = 1.69.

So there is no second boundary to put beside the first. The condition numbers the method declines form a scattered set that begins near 2·10⁸ and thickens as κ grows without ever closing, and the condition numbers it answers at are interleaved through it all the way to 10¹³. Which outcome a particular matrix gets is decided by the sign of a quantity that has already lost every digit it had: the pivots the refusals stop on range from −1.8·10⁻¹⁸ to −4.2·10⁻¹⁵ across the sweep, and the column at which the factorisation stops moves between the fifth, sixth and seventh of eight. Over the same thirty-seven matrices the column sweep’s implied Q never gets further from orthogonal than 9.4·10⁻⁴ and the reduction tree’s than 8.4·10⁻⁴, so every one of these matrices has a factorisation and only one of the three methods cannot find it.

The loss of orthogonality is not monotone either, which is the same fact read on the other side. At 10⁸ the implied Q is 0.373 away from orthogonal and at 1.5·10⁸ it is 0.105 — a harder matrix, answered three times better — and then 2.2·10⁸ is the first refusal and 3.2·10⁸ returns 0.443. Once κ²u is of order one the model that fits so well over twelve decades has nothing left to predict, because it is predicting a quantity that has reached its own ceiling.

And which condition numbers it refuses moves with the shape of the matrix #

Narrow the matrix from eight columns to six and every number above changes. At κ = 10⁸ the implied Q is 0.17 away from orthogonal instead of 0.37, which is the licence being obeyed rather than broken: fewer columns is a shorter Cholesky and less accumulation. But the outcome at 10¹⁰ flips from a returned factor with ‖QᵀQ − I‖ = 1.27 to a refusal on a pivot of −3.6·10⁻¹⁵, and the outcome at 10¹² flips the same way.

The fine sweep measures how much of this there is. Run all thirty-seven condition numbers at six columns and at eight, and the two disagree about the outcome at eleven of them — 30 per cent of the range. Four of the eleven are cases where six columns answer and eight refuse, the remaining seven are the other way round, and they are spread across the whole seven decades rather than clustered at one end. The first refusal itself moves only from 2.2·10⁸ to 4.6·10⁸, two steps of the grid, so the onset barely moves; what moves is everything after it. And what six columns return in place of the refusal is not better: at 3.2·10⁸, where eight columns hand back a factor 0.44 away from orthogonal, six columns hand back one that is 5.15 away.

That is the property that makes the second boundary unusable even as an estimate. A licence has to be a function of what a caller can measure before the work is done, and κ is such a quantity — the condition number is an amplifier and it belongs to the problem. The set of condition numbers this method declines is not a function of κ alone. It is a function of κ, of the column count, and of the rounding of one subtraction near the bottom of a Cholesky, and a caller can compute the first two.

The interval between the two boundaries #

Put the two together. The licence is at κ = 9.49·10⁷ and is computable from the arithmetic alone. The first condition number at which the method declines is 2.2·10⁸ for eight columns and 4.6·10⁸ for six, and after that the outcome is scattered as far as anyone cares to look. Between the two lies a factor of 2.3 in κ in which a run passes its own admission test, is refused nothing, and returns a triangular factor whose implied Q is 0.37 away from orthogonal at 10⁸ and 0.105 at 1.5·10⁸ — and beyond it lie five more decades in which the same thing happens intermittently, eleven times in thirty-seven samples.

Nothing a caller normally prints will separate those runs from good ones. The backward error will not: R is the triangular factor of this matrix to rounding by every route, and the Gram residual ‖AᵀA − RᵀR‖/‖AᵀA‖ is the smallest of the three methods at κ = 10³. That is a small residual is not a small error in its sharpest available form, because here the small residual is not merely uninformative but is the best of the three on the method that is wrong. The quantity that does separate them is ‖QᵀQ − I‖, which costs an n×n product of the factor a caller already has, and which two Gram–Schmidts is the site’s other measurement of.

And the licence’s own input is not exact either. κ is estimated in practice rather than computed, and an estimate that can be fooled is what that estimate is: a condition estimator can be wrong by orders of magnitude on a matrix built to fool it. A predicate whose named boundary sits three decades above the last condition number a caller could have used and a factor of 2.3 below the first one the method declines, evaluated on a number that may itself be a decade out, is not an admission test. It is an assertion with the shape of one.

What a caller does with a predicate that is only half a boundary #

Take the licence as an error model and not as a permission. κ²u predicts ‖QᵀQ − I‖ to within a factor of 7.2 over twelve decades, which is a genuinely useful thing to have. Used that way it needs a tolerance supplied by the caller, and the answer to how conditioned may this matrix be is then a number rather than an inequality with a symbol in it.

Measure the thing the model predicts, rather than inferring it. ‖QᵀQ − I‖ on the returned factor settles in one product what no amount of reasoning about κ will settle, and it is the only quantity in this comparison that distinguishes the three factorisations at all.

Run it twice if the single reduction is worth keeping. Factoring A = Q₁R₁ and then Q₁ = Q₂R₂ takes ‖QᵀQ − I‖ from 0.37 to 3.8·10⁻⁹ at κ = 10⁸ and from 1.27 to 3.7·10⁻⁷ at 10¹⁰, which is the reduction tree’s own orthogonality on the same two matrices, one and two decades past the licence. It is doing it twice applied to a factorisation and is the same identity as buying the accuracy back with a factorisation in place of a precision. The repair costs a second reduction, which is the whole of the objection to it, and it cannot repair a refusal: where the first pass declines there is no Q₁ to factor again.

And treat the refusal as the good outcome. A run that stops on a nonpositive pivot has told the truth about a matrix. A run that returns from the same region has not, and the two are separated by the sign of a number near 10⁻¹⁶. A caller who catches the refusal and falls back to another factorisation has built the branch the licence was supposed to be, out of the one signal in the system that is never wrong when it fires — at the cost of the branch firing at only some of the condition numbers it should.

What is not measured here #

One construction, one arithmetic, one seed. The matrices are built as UΣVᵀ with logarithmically spaced singular values, so κ is prescribed rather than discovered, and a different distribution of singular values would put the refusals in different places — which is the point being made rather than a caveat against it, but it means the specific condition numbers listed here are not a table to look anything up in.

The sweep also stops at 10¹³, and a factor still comes back there at both column counts — 1.49 with six columns and 1.69 with eight. Nothing measured here establishes an upper end to the scattered set, and the claim being made is about the absence of a boundary rather than about the location of one.

And the two counts of columns measured here are eight and six, which is a narrow range chosen to keep the matrices comparable. What the disagreement between them shows is that the column count is one of the arguments the refusal set depends on, not how the set behaves as a function of it. That would be a different measurement, over more shapes, and it would not change the reading: an argument a licence does not mention is an argument a licence cannot account for, whether its effect is orderly or not.