Accurate is not a property of a method

The previous essay measured a claim this site has been carrying in a source comment since its first month and found it true. One-sided Jacobi holds every singular value of a strongly graded bidiagonal to a relative 10⁻¹⁶, over fifty decades, while the route through BᵀB returns the small ones as zero.

This essay is about a second family of matrices, and on it every one of those orderings turns over.

The second family #

Take dᵢ = 1 for every i and eᵢ = 2^k for every i. Every entry of the matrix is within a factor of 4096 of every other, so nothing about it is graded in any sense that a glance would notice.

Its singular values, computed exactly:

4096.9239018723 4096.7071830588 4096.3828136592 4096.0001525879 4095.6174468246 4095.2929695291 4095.0761428197 5.1698785203·10⁻²⁶

Seven of them clustered within a twentieth of a per cent of each other, and one twenty-nine decades below. The spectrum spans what the graded family’s spanned; the entries do not.

Nothing about the construction is contrived. It is a bidiagonal with a large superdiagonal, which is what a badly scaled discretisation produces, what a companion matrix of a polynomial with large coefficients looks like, and what any process that couples neighbours much more strongly than it anchors them gives — the same shape the units the matrix is measured in produces from a change of scale rather than from a coupling.

What each route does #

One-sided Jacobi returns the seven large values to 4.4·10⁻¹⁶ — exactly, in the sense that matters — and σ_min wrong by a relative 1.45·10⁻². One and a half per cent, on the one value the method is recommended for.

The zero-shift sweep, which has the theorem behind it, does not converge. Four hundred sweeps and it is still running, with a worst relative error of 2·10⁻⁴ when it is stopped.

The shifted sweep, which the theorem is a warning about, converges in sixteen sweeps with a worst relative error of 5.6·10⁻¹⁶ across all eight values.

And BᵀB loses σ_min entirely, as it did before — the one ordering that does not reverse.

The progression across the slider is a clean ladder: Jacobi’s relative error on σ_min runs 1.6·10⁻¹³, 9.4·10⁻¹⁰, 1.3·10⁻⁷, 3.4·10⁻⁵, 1.5·10⁻² as the superdiagonal doubles three bits at a time, while its error on every larger value stays at 10⁻¹⁵ throughout.

The smallest singular value has a closed form and the slider confirms it. Across superdiagonals of 2⁴, 2⁶, 2⁸, 2¹⁰ and 2¹² it reads 3.711·10⁻⁹, 2.273·10⁻¹³, 1.388·10⁻¹⁷, 8.47·10⁻²² and 5.17·10⁻²⁶ — which is 2^(−7k) at every one, that is e^−(n−1) with e the superdiagonal and n = 8. Eighty-four binary orders of magnitude walked in five steps, with the exponent read off the drawing rather than asserted.

And κ(X) is two over it. The essay’s three values of κ(X) — 5.4·10⁸, 1.4·10¹⁷, 3.9·10²⁵ — divided into 1/σ_min at the same superdiagonals give 2.00, 1.94 and 2.02. So on this family the hypothesis of the theorem and the smallness of the answer are the same number, which is why no diagonal factoring helps and why the theorem correctly predicts nothing.

What the bound does not predict is how badly. u·κ(X) is 6·10⁻⁸ at the mildest member and 4.3·10⁹ at the fiercest — the second is vacuous, since a relative error above one means nothing at all is known. Jacobi’s measured error at those two points is 1.6·10⁻¹³ and 1.5·10⁻², so it is four hundred thousand times better than the bound at one end and two hundred billion times better at the other. The error grows as about 2^(4.6k) against the bound’s 2^(7k): the theorem is right that no relative accuracy is guaranteed and wrong by a growing margin about the size of the loss.

The shifted sweep does not notice any of it. Across the five superdiagonals it takes 17, 17, 17, 16 and 16 sweeps and returns worst relative errors of 8.66·10⁻¹⁶, 1.67·10⁻¹⁵, 7.8·10⁻¹⁶, 9.99·10⁻¹⁶ and 5.55·10⁻¹⁶. Eighty-four binary orders of magnitude in the answer, a hypothesis that fails harder at every step, and a method with no relative-accuracy theorem behind it that is flat in both columns. The zero-shift sweep, which has the theorem, fails to converge in four hundred sweeps at every one of the five.

Why none of this contradicts the theory #

It would be easy to read the above as a counterexample to a theorem, and it is not one. It is a counterexample to how the theorem is remembered.

The statements about relative accuracy — Demmel and Kahan’s, Demmel and Veselić’s — are not of the form, in the way a bound that is never attained is not the form its slogan takes, “algorithm M computes small singular values accurately”. They are of the form: write A = D·X with D diagonal and X of unit rows; then the computed singular values have relative error bounded by a modest multiple of u·κ(X).

κ(X), not κ(A). The hypothesis is about the matrix after the diagonal has been factored out, and the diagonal is exactly what a grading is.

That is the number to compute, and it is decisive.

Along the graded family, κ(A) climbs from 271 to 7.0·10⁵⁰ — fifty decades — and κ(X) is 4.892 at every one of the six matrices. Not approximately constant: constant, because the grading is precisely what the diagonal factor absorbs and there is nothing left over.

Along the uniform family, κ(X) climbs 5.4·10⁸, 1.4·10¹⁷, 3.9·10²⁵ as the superdiagonal grows. There is no diagonal factor to take out, so the ill-conditioning is in X where the theorem can see it, and the theorem correctly predicts that no relative accuracy is available.

So both families obey the theory exactly. What they do not obey is the sentence the theory gets compressed into.

Where each method’s mechanism runs out #

It is worth saying, for each of the three sweeps, exactly which step stops working on the second family — because “the hypothesis fails” is a statement about a theorem and the mechanisms are statements about code.

One-sided Jacobi orthogonalises pairs of columns, and its accuracy argument is that a rotation of two columns preserves the relative accuracy of both. That argument is about the columns. Here the columns are all of comparable norm — every one of them is about 4096 — and the smallness lives in their near-dependence rather than in any of their sizes. The rotation angle is computed from a difference of two nearly equal inner products, and the smallness has to survive that subtraction. It does not.

The zero-shift sweep has no subtraction in it, and its accuracy is intact. What it does not have is convergence: the whole point of a shift is to make a QR iteration converge quickly, and on a matrix whose seven large singular values are within a twentieth of a per cent of each other there is nothing for an unshifted iteration to separate. It is a correct algorithm that does not finish.

The shifted sweep subtracts, and the subtraction is d² − μ. On the graded family that would be dangerous, because μ is near the smallest value and the small entries are what the accuracy depends on. Here μ is near 4096² and the entries it is subtracted from are near 4096², so the cancellation happens in the cluster — where the values genuinely are nearly equal and there is nothing to preserve — and the small value at the bottom of the spectrum is never involved in a subtraction at all.

Three mechanisms, each correct, each failing or succeeding for a different reason, and none of them predictable from the sentence “computes small singular values accurately”.

What the classifier is worth, and what it costs #

It costs a singular value decomposition, which is more than the answer it is guarding. Nobody runs it as a pre-check.

Except that a free version separates the families anyway #

κ(X) needs a decomposition. The question it answers — is this matrix’s difficulty something a diagonal scaling can absorb — has a one-pass answer that this collection has already measured in a condition number scaling cannot move: the spread of the row norms.

   family    parameter    κ(A)      κ(X)      row spread    κ(A) ⁄ κ(X)
   graded       10       6.5·10¹⁰   4.89       3.4·10¹⁰      1.3·10¹⁰
   graded       30       6.0·10²⁹   4.89       3.2·10²⁹      1.2·10²⁹
   graded       50       7.0·10⁵⁰   4.89       3.7·10⁵⁰      1.4·10⁵⁰
   uniform       4       4.6·10⁹    5.4·10⁸    1.6·10¹       8.5
   uniform      12       8.0·10²⁸   3.9·10²⁵   4.1·10³       2.1·10³
   uniform      16       5.1·10³⁵   1.1·10³⁴   6.6·10⁴       4.9·10¹

The row spread bounds κ(A)/κ(X) in every case and is within a factor of three of it on the graded family — because there the whole of the conditioning is the grading, which is what κ(X) = 4.89 at every decade says.

And the two families are forty orders apart on it. The graded family’s spread starts at 3.4·10¹⁰ and rises; the uniform family’s ends at 6.6·10⁴. No threshold placed anywhere in between could be a delicate choice, and there is nothing to tune.

So the pre-check nobody runs is one pass over the matrix rather than a decomposition. It does not give κ(X), and it does not have to: the question is which family a matrix is in, and the answer is readable from the row norms before anything is factorised.

That is the same diagnostic, doing the same job, in a third field — the scaling essay uses it to predict how much equilibration will buy, and here it predicts whether relative accuracy is available. Both are asking how much of this matrix’s difficulty is a choice of units, and the row norms answer it for free in both.

There is a caveat that keeps it honest, and it is the last row of the table. At a superdiagonal of 2¹⁶ the uniform family’s spread is 6.6·10⁴ while κ(A)/κ(X) is only 49 — the bound is loose by three orders there, because the row norms differ by 6.6·10⁴ while the row directions are what κ(X) measures and they have stopped being informative. So the proxy is a bound and not an estimate, and its looseness runs in the safe direction: it can say equilibration might buy a lot when it will buy little, and it cannot say equilibration will buy little when it would buy a lot.

For the classifier’s purpose that is exactly the right direction to be wrong in. A large row spread means look further; a small one means there is no grading here, and the small singular values are not determined by the entries. The second is the verdict that matters, it is the one the essay’s uniform family is about, and it is the one the free quantity gets right every time.

What it is for is knowing which question has an answer. Two matrices with the same κ, the same σ_min and the same size can be in different families, and the difference decides whether “the small singular values, to relative accuracy” is a request that any algorithm can satisfy or a request that no algorithm can. That is not a question about which method to use. It is a question about whether the quantity being asked for is determined by the data at all.

Stated that way it joins a group this collection has been building for some time.

Rank is a decision, because a floating-point matrix does not have one. A regularisation parameter is a choice, because the data does not determine the answer. And now: relative accuracy in the small singular values is available or not, as a property of the matrix, and no amount of algorithmic care creates it where it is absent.

And a smaller lesson about defaults #

The switch a library makes between the two sweeps is worth looking at, because it is what a sensible answer to this looks like in practice.

LAPACK’s bidiagonal singular value routine does not always use the zero shift and does not always use a shift. It estimates whether the matrix’s smallest and largest values are far enough apart that a shift would be dangerous, and picks. That estimate is a crude relative of the classifier above, made from quantities the routine already has.

Two things follow that this collection would say about any such switch.

A switch made on a measurement is not a compromise. It is a routine that knows which of two regimes it is in and acts accordingly, which is strictly better than either fixed choice — as the two families here demonstrate, since each fixed choice is catastrophic on one of them.

And the criterion is where the interesting failure lives. The sweeps are both correct. What can be wrong is the decision about which to use, and that decision is a threshold — which is the subject the last two essays of this collection’s current run are about.

The reversal is not symmetric #

One route does not reverse, and it is worth being explicit about which.

The eigenvalues of BᵀB lose the smallest singular value on both families. On the graded one they return zero; on the uniform one they return zero as well. There is no matrix in this essay or the previous one on which forming BᵀB is the right thing to do, and there is no family constructed here on which it recovers.

That asymmetry is the reason the previous essay’s conclusion survives this one. Two of the four routes have a hypothesis, and this essay is about the hypothesis being checkable rather than about the routes being unreliable. The third has no hypothesis and no case in which it is preferable: the squaring destroys the small end whatever else is true, and the only thing that varies is how much grading it takes before the destruction is visible.

So the practical summary is not “it depends”. It is: never form BᵀB if the small end matters; otherwise compute κ of the row-equilibrated matrix and let it choose between the other two.

Two matrices that look identical #

The pair worth remembering from this essay is not either family; it is a comparison that a caller could be handed.

Matrix one: 8×8 bidiagonal, entries between 2⁻⁹⁸ and 1, κ = 6·10²⁹, σ_min = 2.1·10⁻³⁰.

Matrix two: 8×8 bidiagonal, entries 1 and 4096, κ = 8·10²⁸, σ_min = 5.2·10⁻²⁶.

Same size, same shape, condition numbers within a decade of each other, smallest singular values within four. Every summary statistic anybody prints about a matrix says these two are the same kind of problem.

On the first, one-sided Jacobi gets σ_min to sixteen digits. On the second it gets it to two.

The number that separates them is κ of the row-equilibrated matrix: 4.89 against 3.9·10²⁵. It is computed by nobody, printed by nothing, and it is the only quantity in the paragraph above that predicts what will happen.

What the exact route made possible #

None of this could have been measured without an answer that had not been rounded, and it is worth saying what the exact route bought, because it was not merely a tighter tolerance.

Every number above is a relative error against a value of size 10⁻²⁶. Two float routes disagreeing about such a value cannot be adjudicated: each is entitled to be wrong by its own size. What the rational bisection provides is not a better estimate — it is a different kind of object, an answer whose relative accuracy was chosen in advance and is limited only by how many bisection steps were taken.

That is what allowed a 1.45 per cent error to be identified as an error rather than as a disagreement, and it is the strongest case this site has made for keeping exact arithmetic beside the floating-point work. The exact route runs on 8×8 matrices and would not run on 800×800 ones. What it establishes is a fact about the methods, and a method does not change size.

What this changes about the rest of the site #

The claim being checked is load-bearing here in a way it is not on most sites, so it is worth stating what the check licenses and what it does not.

Every rank decision, condition number and gap on this site is computed by one-sided Jacobi, and every one of them is drawn on a matrix built with a prescribed spectrum or a prescribed condition number — which means built as UΣVᵀ with the grading in Σ, which is the first family. The comment on the routine was right about the matrices this site draws, and the essays that rest on it rest on solid ground.

What the second family says is that the licence does not extend to a matrix somebody else hands over. A rank decision on a discretised operator with a large off-diagonal coupling, or on a companion matrix, or on any of the objects whose smallness is a product rather than a size, is a rank decision made on singular values that may be wrong in the second digit — and there is nothing in the spectrum itself that says which case it is.

The honest form of that is a note rather than a repair: the site’s own routine is the right default and its guarantee is conditional, and the condition is computable. It is now written down where it can be found rather than in a comment on a function.

The refusal #

The assertion is fed an infinite entry handed to the exact converter.

The converter’s job is to turn a double into a rational with nothing lost, and it can, because every finite double is a dyadic rational. An infinity is not, and the tempting behaviour — substitute the largest finite value, or clamp — produces a reference that is exact about a different matrix.

That is the failure mode of exactness used carelessly and it is the one this collection keeps guarding: an answer that is perfectly correct about an object nobody asked about, carrying no evidence of the substitution that made it so. A relative error measured against such a reference is a number about something else, and it looks exactly like a number about this.

What is next #

Three essays about what a small quantity means and two about whether it can be computed at all. The next pair changes subject once more, to the sentence this whole site is built on — that a good algorithm returns the exact answer to a nearby problem — and asks a question nobody has asked here in a hundred and eighteen essays: nearby problem of what kind?