Numerical Linear Algebra
where the algebra stops being the arithmetic
Eigenvalues, singular values, rank

Accurate is not a property of a method

A bidiagonal matrix whose every entry is 1 or 4096 has singular values spanning thirty decades. On it, the method recommended for small singular values loses the small one by one and a half per cent, the sweep with the theorem behind it does not converge at all, and the shift the theorem is a warning about gets every value to 5·10⁻¹⁶. Nothing there contradicts the theory.
28 min read 21 figures Identical algebra, different arithmeticWhose fault is itExact ground truthMeasured, not asserted

Worth reading first: Small compared to what · The units the matrix is measured in · The condition number is an amplifier.

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 same four routes on a bidiagonal whose entries are all 1 or 2^12Nothing about this matrix is graded: every entry is within a factor of 4096 of every other. Its singular values are not — 7 of them are clustered near 4096.9 and the last is 5.17·10⁻²⁶, a spread of 29 decades. Every ordering from the graded family turns over. One-sided Jacobi, which held everything there, returns the 7 large values exactly and σ_min wrong by a relative 0.0145. The zero-shift sweep, which has the theorem behind it, does not converge at all — 400 sweeps and still running. The shifted sweep, which the theorem is a warning about, gets every value to 5.55·10⁻¹⁶ in 16.1234567810⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹singular value, largest firstrelative errorone-sided Jacobizero-shift QRshifted QReigenvalues of BᵀBthe same four routes, reversedσ_min, exactly5.2·10⁻²⁶Jacobi's error on it0.015sweeps, zero shift400sweeps, shifted16a method is not accuratea method on a matrix is
Fig. 1 A bidiagonal whose diagonal is all ones and whose superdiagonal is all 4096. The same four routes, and a different verdict on each of them.

The second family

Take d_i = 1 for every i and e_i = 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.

How far |r_nn| sits above σ_min on Kahan's matrix, against the size and the parameter3 curves of |r_nn| ÷ σ_min against n, one per Kahan parameter. Every curve rises without turning over, reaching 7·10⁴ at n = 30, c = 0.5. Column pivoting makes no interchange at any point on any of them, so the failure is not a poor choice — there is nothing to choose.813182328110¹10²10³10⁴10⁵10⁶size of the matrix|r_nn| ÷ σ_minthe two agreec = 0.2c = 0.35c = 0.5no ceilingratio at n = 1021ratio at n = 307·10⁴interchanges, anywhere0the greedy rule never had a choiceand the gap grows with every row
Fig. 2 A close relative from the rank-revealing essays: a matrix whose smallness is a product of couplings rather than a property of any entry, and which fools a cheaper method for the same reason.
Wilkinson's 12×12 matrix and its upper triangular factorThe matrix on the left has ones on the diagonal, minus ones below it and a column of ones at the right. On the right, its U factor, whose last column doubles down the rows to sixty-four.1··········1-11·········1-1-11········1-1-1-11·······1-1-1-1-11······1-1-1-1-1-11·····1-1-1-1-1-1-11····1-1-1-1-1-1-1-11···1-1-1-1-1-1-1-1-11··1-1-1-1-1-1-1-1-1-11·1-1-1-1-1-1-1-1-1-1-111-1-1-1-1-1-1-1-1-1-1-11A1··········1·1·········2··1········4···1·······8····1······16·····1·····32······1····64·······1···128········1··256·········1·512··········11024···········2048U‖PA − LU‖/‖A‖0growth factor2048the 2ⁿ⁻¹ bound2048Row interchanges performed: 0.Partial pivoting had nothing to choose.every entry is 0, 1 or −1the bound is attained here
Fig. 3 And the classic of the type from the elimination field, where the same accumulation of ordinary-sized entries produces something enormous.

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 same four routes on a bidiagonal whose entries are all 1 or 2^4Nothing about this matrix is graded: every entry is within a factor of 16 of every other. Its singular values are not — 7 of them are clustered near 16.929 and the last is 3.711·10⁻⁹, a spread of 10 decades. Every ordering from the graded family turns over. One-sided Jacobi, which held everything there, returns the 7 large values exactly and σ_min wrong by a relative 1.6·10⁻¹³. The zero-shift sweep, which has the theorem behind it, does not converge at all — 400 sweeps and still running. The shifted sweep, which the theorem is a warning about, gets every value to 8.66·10⁻¹⁶ in 17.1234567810⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹singular value, largest firstrelative errorone-sided Jacobizero-shift QRshifted QReigenvalues of BᵀBthe same four routes, reversedσ_min, exactly3.7·10⁻⁹Jacobi's error on it1.6·10⁻¹³sweeps, zero shift400sweeps, shifted17a method is not accuratea method on a matrix is
Fig. 4 At a superdiagonal of sixteen the smallest singular value is 3.7·10⁻⁹ and Jacobi still gets it to twelve digits: the hypothesis has not failed yet.
The same four routes on a bidiagonal whose entries are all 1 or 2^8Nothing about this matrix is graded: every entry is within a factor of 256 of every other. Its singular values are not — 7 of them are clustered near 256.92 and the last is 1.388·10⁻¹⁷, a spread of 19 decades. Every ordering from the graded family turns over. One-sided Jacobi, which held everything there, returns the 7 large values exactly and σ_min wrong by a relative 1.32·10⁻⁷. The zero-shift sweep, which has the theorem behind it, does not converge at all — 400 sweeps and still running. The shifted sweep, which the theorem is a warning about, gets every value to 7.8·10⁻¹⁶ in 17.1234567810⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹singular value, largest firstrelative errorone-sided Jacobizero-shift QRshifted QReigenvalues of BᵀBthe same four routes, reversedσ_min, exactly1.4·10⁻¹⁷Jacobi's error on it1.3·10⁻⁷sweeps, zero shift400sweeps, shifted17a method is not accuratea method on a matrix is
Fig. 5 At 256 it is 1.4·10⁻¹⁷ and Jacobi is out by 1.3·10⁻⁷. Watching the slider is watching a hypothesis fail.

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.

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 “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.

κ(A) and κ of the row-equilibrated matrix, along both families, n = 8Write A = D·X with D diagonal and every row of X of unit norm. Every theorem about relative accuracy is a hypothesis on κ(X), and κ(A) appears in none of them — which is easy to read past and is the whole difference between the two families here. Along the graded family κ(A) climbs from 271 to 7.04·10⁵⁰ and κ(X) is 4.892 at every one of the six matrices — the grading is exactly what the diagonal factor absorbs. Along the uniform family κ(X) climbs to 3.87·10²⁵. That is the number that says which question has an answer, and it is not the number anybody prints.01020304050110¹⁰10²⁰10³⁰10⁴⁰10⁵⁰log₁₀ κ(A)log₁₀ of the condition numberκ(X) = κ(A)graded familyuniform familyone number decides, and it is not κκ(X), graded, at every grading4.9κ(A), graded, at the widest7·10⁵⁰κ(X), uniform, at the widest3.9·10²⁵κ(A) there8·10²⁸κ is a fact about the matrixand the hypothesis is about a factor of it
Fig. 6 κ(A) across the bottom, κ of the row-equilibrated matrix up the side. The dashed diagonal is where they would be equal.

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.

κ(A) and κ of the row-equilibrated matrix, along both families, n = 12Write A = D·X with D diagonal and every row of X of unit norm. Every theorem about relative accuracy is a hypothesis on κ(X), and κ(A) appears in none of them — which is easy to read past and is the whole difference between the two families here. Along the graded family κ(A) climbs from 4350 to 8.83·10⁴⁹ and κ(X) is 5.653 at every one of the six matrices — the grading is exactly what the diagonal factor absorbs. Along the uniform family κ(X) climbs to 1.45·10³¹. That is the number that says which question has an answer, and it is not the number anybody prints.01020304050110¹⁰10²⁰10³⁰10⁴⁰10⁵⁰log₁₀ κ(A)log₁₀ of the condition numberκ(X) = κ(A)graded familyuniform familyone number decides, and it is not κκ(X), graded, at every grading5.7κ(A), graded, at the widest8.8·10⁴⁹κ(X), uniform, at the widest1.5·10³¹κ(A) there2.4·10³¹κ is a fact about the matrixand the hypothesis is about a factor of it
Fig. 7 At n = 12 the constant is a different constant and the shape is the same shape, which is what makes it a classifier rather than a number that came out right twice.

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”.

Classical and modified Gram–Schmidt: the same subtraction, in a different orderTwo panels of pseudocode differing in one argument, with the resulting pairwise dot products of the computed columns listed beneath each.for each previous column i, subtract the projection of column j onto qᵢclassicalr[i][j] = qᵢ · a[j] ↑ the ORIGINAL columnv = v − r[i][j] · qᵢthe three worst |qᵢ · qⱼ|:columns 5 and 6: 5·10⁻⁵columns 4 and 6: 7.1·10⁻⁷columns 4 and 5: 10⁻⁸modifiedr[i][j] = qᵢ · v ↑ what is LEFT of itv = v − r[i][j] · qᵢthe three worst |qᵢ · qⱼ|:columns 1 and 6: 1.8·10⁻¹⁰columns 2 and 6: 7.3·10⁻¹¹columns 3 and 6: 7.5·10⁻¹²The two R factors agree to 1.2·10⁻¹¹ relative. The two Q factors do not.the 6×6 Hilbert matrixone word, eight orders
Fig. 8 The site’s original instance of the same shape, and its cleanest: two orderings of one derivation with a difference the algebra cannot see.
Loss of orthogonality against condition number, in binary64A log–log plot of the norm of Q-transpose-Q minus the identity against condition number. Classical Gram–Schmidt rises steeply, modified Gram–Schmidt rises gently, and Householder is flat.110²10⁴10⁶10⁸10¹⁰10¹²10⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹condition number κ(A)‖QᵀQ − I‖classicalmodifiedHouseholderκ²uκu8×8, eight seeds per κ, binary64all three reconstruct A
Fig. 9 And the number that separates them once it is printed rather than assumed.

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.

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.

Singular values of a rank-4 matrix with noise of relative size 10⁻⁶Ten singular values on a logarithmic axis. The first four sit near one; the rest sit at the noise level, and the vertical distance between the two groups is the evidence for the rank.1234567891010⁻¹⁸10⁻¹⁵10⁻¹²10⁻⁹10⁻⁶10⁻³1indexsingular valuecutoff, σ₁ · 10⁻¹⁰numerical rank 10gap 8.2·10⁴an opiniontrue rank 410×10, built with 4 nonzero valuesrank is a decision
Fig. 10 The first of those, and the one this site drew earliest.
The widest box a uniqueness proof is available forVerified half-width against the distance between two roots, both axes logarithmic. The measured widths are 0.005, 0.015, 0.05 at separations of 0.01, 0.03, 0.1 — a ratio of 0.5000 at every one of them. The operator verifies exactly up to the point where the box would reach the second root, and a box containing both is refused.10⁻²10⁻¹10⁻²10⁻¹distance between the two rootswidest verified half-widththe separationverified widthuniqueness is a claim about a regionwidth ÷ separation at δ = 0.010.5width ÷ separation at δ = 0.030.5width ÷ separation at δ = 0.10.5verified up to the second rootand refused past it
Fig. 11 And the second: how wide the set of answers consistent with the data is, which is what makes the choice a choice.

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.

Convergence factor against the strength threshold, ε = 0.1A curve of convergence factor against the strength threshold on a logarithmic axis. It is flat at a high value to the left of a marked position and flat at a low value to the right of it.10⁻²10⁻¹100.20.40.60.8strength threshold θresidual reduction per cycleθ = εsemi-coarseningkept whole rowsfull coarseningone parameter, two methodsbest factor above ε0.048best factor below ε0.18the ratio across the switch3.831×31 anisotropic operatora switch, not a dial
Fig. 12 A switch of exactly that kind from the algebraic multigrid essays, where a threshold decides which of two behaviours the method has.
Two subdivision rules on roots 0.1 apart, one of which sits on the first cutSix bars. Cutting a box at its midpoint puts the root at x = 1 on a boundary of every box the search ever asks about, and the operator cannot verify a root it never sees the interior of: 0 verified after 163 evaluations. Cutting at 0.485 of the width verifies both roots in 45.both rules see the same problem and the same operatormidpoint: verified0midpoint: undecided6midpoint: evaluations163off-centre: verified2off-centre: undecided0off-centre: evaluations45where the box is cutroots verified, midpoint0roots verified, off-centre2evaluations saved118a root on a boundary is never in an interiorand 0.485 is the whole repair
Fig. 13 And another, from the branch-and-bound essay: a rule that decides rather than computes, and whose threshold is the whole content.

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.

Fitting the same degree-10 polynomial in two basesOn the left, the data and two fitted curves that lie on top of each other. On the right, the condition numbers of the two design matrices, ten orders of magnitude apart.00.10.20.30.40.50.60.70.80.910.20.40.60.81xyboth fits, drawn on top of one anothercondition number of the design matrixmonomial2.1·10⁷Chebyshev2.3largest fitted coefficientmonomial176Chebyshev0.51rms residual: 4.2·10⁻⁶ and 3.1·10⁻⁷ — the data is fitted either way.30 points, degree 10, single precisionthe basis is part of the problem
Fig. 14 Where the same squaring appears in the least-squares field, as a choice of basis rather than a choice of algorithm.

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.

Two perturbation bounds and the error that was measured, on a 12×12 matrix spread over 8 decades of unitsThree curves against the size of an entrywise relative perturbation. The normwise bound κ_∞·ε is a valid bound and sits 3.7·10⁷ times above the componentwise one cond(A, x)·ε, which is also a bound and is nearly attained by the worst of forty random perturbations at each size.10⁻¹⁴10⁻¹³10⁻¹²10⁻¹¹10⁻¹⁰10⁻⁹10⁻¹⁵10⁻¹²10⁻⁹10⁻⁶10⁻³1relative size of the entrywise perturbationrelative forward errorκ_∞ · εcond(A,x) · εmeasuredboth bounds holdκ_∞(A)3.1·10⁸cond(A, x)8.5ratio of the bounds3.7·10⁷both curves above the data are boundsand only one of them is a measurement
Fig. 15 Its nearest relative on this site, from the scaling field: a condition number that a change of units cannot move, and the diagnostic its ratio to κ provides.
Two condition numbers of one 8×8 system, as its rows are put into different unitsFour curves against the spread of the row units, in decades. κ_∞ of the scaled matrix rises from 9.83 to 223 while the componentwise condition number stays at 6.98 throughout — the same system, the same solution, and one of the two numbers is a fact about the units. Hilbert's two numbers are drawn flat beside them at 3.4·10¹⁰ and 1.2·10¹⁰: a matrix whose sensitivity no scaling repairs.012110²10⁴10⁶10⁸10¹⁰10¹²spread of the row units (decades)condition numberκ_∞(DA)cond(DA)Hilbert κ_∞Hilbert condone system, two numbersκ_∞ at no spread9.8κ_∞ at 2 decades223cond, either end7Hilbert, equilibrated1.3·10¹⁰the solution is the same at every spreadand one of these curves knows it
Fig. 16 And what a change of units does to the ordinary one, which is what makes the ordinary one an unreliable summary.
Four pivot rules on [[ε, 1], [1, 1]] with ε = 10⁻¹², before and after its first row is scaled by 1/εEight bars of forward error. On the matrix as given, no pivoting loses the answer entirely and all three pivot rules return it to rounding. Multiply the first row by 1/ε — which does not change the solution — and partial pivoting makes no interchange at all and returns the same wrong answer that no pivoting does. Scaled partial pivoting and complete pivoting are unmoved.forward error, relative to a solution of exactly (1, 1)no pivoting · as given1.3·10⁻⁴ 0 interchangesno pivoting · rows scaled1.2·10⁻⁴ 0 interchangespartial · as given2.2·10⁻¹⁶ 1 interchangepartial · rows scaled1.2·10⁻⁴ 0 interchangesscaled partial · as given2.2·10⁻¹⁶ 1 interchangescaled partial · rows scaled2.2·10⁻¹⁶ 1 interchangecomplete · as given1.1·10⁻¹⁶ 1 interchangecomplete · rows scaled2.2·10⁻¹⁶ 1 interchangethe same problem twicepartial, as given2.2·10⁻¹⁶partial, rows scaled1.2·10⁻⁴its relative residual1.2·10⁻¹⁶complete, rows scaled2.2·10⁻¹⁶the two systems have the same solutionand one pivot rule cannot see it
Fig. 17 The same idea one field over: a pivoting rule whose verdict depends on units nobody recorded.

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.

The exact solution of a 13×13 Hilbert system beside the computed oneTwo columns of numbers: the exact answer, which is the integers one to thirteen, and the answer double-precision elimination returns, with the number of correct digits beside each.H13 x = b, b formed exactly so that x = (1, 2, …, 13)12345678910111213exact1.00002.00003.00133.97865.18864.993010.46720.048921.2657-2.575319.21478.906013.5113computed6.6 correct digits4.8 correct digits3.4 correct digits2.3 correct digits1.4 correct digit0.8 correct digitno correct digitsno correct digitsno correct digitsno correct digitsno correct digits0.6 correct digit1.4 correct digitbackward error2.2·10⁻¹⁷κ = 1.7·10¹⁸The algorithm solved a neighbouring problem perfectly. That problem's answer is this one.right-hand side built in BigInt rationalsthe truth is known
Fig. 18 The original of the habit, on a solution rather than a spectrum, and the reason this collection built an exact arithmetic in its first month.

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.

The same four routes on a bidiagonal whose entries are all 1 or 2^6Nothing about this matrix is graded: every entry is within a factor of 64 of every other. Its singular values are not — 7 of them are clustered near 64.925 and the last is 2.273·10⁻¹³, a spread of 14 decades. Every ordering from the graded family turns over. One-sided Jacobi, which held everything there, returns the 7 large values exactly and σ_min wrong by a relative 9.44·10⁻¹⁰. The zero-shift sweep, which has the theorem behind it, does not converge at all — 400 sweeps and still running. The shifted sweep, which the theorem is a warning about, gets every value to 1.67·10⁻¹⁵ in 17.1234567810⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹singular value, largest firstrelative errorone-sided Jacobizero-shift QRshifted QReigenvalues of BᵀBthe same four routes, reversedσ_min, exactly2.3·10⁻¹³Jacobi's error on it9.4·10⁻¹⁰sweeps, zero shift400sweeps, shifted17a method is not accuratea method on a matrix is
Fig. 19 At a superdiagonal of sixty-four the error on σ_min is 9·10⁻¹⁰: the ladder is continuous and the hypothesis fails gradually.
κ(A) and κ of the row-equilibrated matrix, along both families, n = 6Write A = D·X with D diagonal and every row of X of unit norm. Every theorem about relative accuracy is a hypothesis on κ(X), and κ(A) appears in none of them — which is easy to read past and is the whole difference between the two families here. Along the graded family κ(A) climbs from 66.7 to 8.7·10⁴⁹ and κ(X) is 4.273 at every one of the six matrices — the grading is exactly what the diagonal factor absorbs. Along the uniform family κ(X) climbs to 2.31·10¹⁸. That is the number that says which question has an answer, and it is not the number anybody prints.01020304050110¹⁰10²⁰10³⁰10⁴⁰10⁵⁰log₁₀ κ(A)log₁₀ of the condition numberκ(X) = κ(A)graded familyuniform familyone number decides, and it is not κκ(X), graded, at every grading4.3κ(A), graded, at the widest8.7·10⁴⁹κ(X), uniform, at the widest2.3·10¹⁸κ(A) there4.7·10²¹κ is a fact about the matrixand the hypothesis is about a factor of it
Fig. 20 And the classifier at six rows, where the constant is a different constant and the shape is the same.
The same four routes on a bidiagonal whose entries are all 1 or 2^10Nothing about this matrix is graded: every entry is within a factor of 1024 of every other. Its singular values are not — 7 of them are clustered near 1024.9 and the last is 8.47·10⁻²², a spread of 24 decades. Every ordering from the graded family turns over. One-sided Jacobi, which held everything there, returns the 7 large values exactly and σ_min wrong by a relative 3.39·10⁻⁵. The zero-shift sweep, which has the theorem behind it, does not converge at all — 400 sweeps and still running. The shifted sweep, which the theorem is a warning about, gets every value to 9.99·10⁻¹⁶ in 16.1234567810⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹singular value, largest firstrelative errorone-sided Jacobizero-shift QRshifted QReigenvalues of BᵀBthe same four routes, reversedσ_min, exactly8.5·10⁻²²Jacobi's error on it3.4·10⁻⁵sweeps, zero shift400sweeps, shifted16a method is not accuratea method on a matrix is
Fig. 21 At 1024 the error on the smallest singular value is 3·10⁻⁵, one rung below the picture at the top of this essay.

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?

Where this goes next

1234567810⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹singular value, largest firstrelative errorone-sided Jacobizero-shift QRshifted QReigenvalues of BᵀBagainst a rational bisectionσ_min, exactly2.1·10⁻³⁰worst, one-sided Jacobi4.4·10⁻¹⁶worst, zero-shift QR2.2·10⁻¹⁶worst, eigenvalues of BᵀB1a relative error is a ratioand the denominator is the answer Small compared to what This site's own singular value routine has carried a sentence since the month it was written — that one-sided Jacobi computes the small singular values to high relative accuracy and the standard method does not. It has never been measured here, because measuring it needs a σ that is known rather than computed. A bidiagonal matrix and a Sturm count in exact rationals supply one. 10⁻¹⁴10⁻¹³10⁻¹²10⁻¹¹10⁻¹⁰10⁻⁹10⁻¹⁵10⁻¹²10⁻⁹10⁻⁶10⁻³1relative size of the entrywise perturbationrelative forward errorκ_∞ · εcond(A,x) · εmeasuredboth bounds holdκ_∞(A)1.9·10⁸cond(A, x)4.8ratio of the bounds4·10⁷both curves above the data are boundsand only one of them is a measurement A condition number scaling cannot move Skeel's componentwise condition number is invariant under any row scaling — exactly, before any norm is taken, because two diagonal factors cancel entry by entry. It is never larger than the normwise one and can be arbitrarily smaller, and the ratio between them is a diagnostic for which kind of ill-conditioning a matrix has. 0246810110²10⁴10⁶10⁸10¹⁰10¹²spread of the row units (decades)condition numberκ_∞(DA)cond(DA)Hilbert κ_∞Hilbert condone system, two numbersκ_∞ at no spread9.8κ_∞ at 10 decades1.9·10¹⁰cond, either end7Hilbert, equilibrated1.3·10¹⁰the solution is the same at every spreadand one of these curves knows it The units the matrix is measured in One linear system, written twice. The rows of the second are the rows of the first in different units, the solution is identical to the last bit, and the condition number has moved by eight orders of magnitude. One of those two numbers is a fact about the problem and the other is a fact about the notation.

Shares its objects with

Essays that name at least two of the same things, and that neither author linked.

Named objects

A flat tag is an object no other essay names yet.

Bidiagonal matrixCondition numberConvergence rateEquilibrationExact arithmeticGraded matrixJacobi's eigenvalue methodRelative accuracySingular values