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

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.
26 min read 22 figures Exact ground truthTwo routes to a numberMeasured, not assertedWhere the algebra stops being the arithmetic

Worth reading first: The best approximation there is · Rank is a decision · An answer that is known.

Every routine on this site that needs a singular value calls the same function, and that function has had a comment on it since the month it was written:

One-sided Jacobi rotates pairs of columns of A until they are mutually orthogonal … It computes the small singular values to high relative accuracy, which the standard bidiagonal method does not — and since half this site is about small singular values, that is the method it should use.

Every word of that is correct. It is also, until now, the only assertion in this collection that has never been given a test it could fail. Half the site’s arguments — rank decisions, condition numbers, gaps, sketching error — rest on it, and it has been carried, all that time, as a sentence.

The reason is not carelessness. It is that measuring it requires a singular value that is known.

Relative error of each singular value of a bidiagonal matrix graded over 30 decades, by four routesThe matrix is 8×8, bidiagonal, and every entry is a power of two times a small integer, so the matrix the float routines are handed and the matrix the exact routine reasons about are the same matrix. The reference is a Sturm bisection in BigInt rationals, which never rounds. σ_min is 2.098·10⁻³⁰. One-sided Jacobi holds every value to 4.41·10⁻¹⁶ relative and the zero-shift bidiagonal sweep to 2.21·10⁻¹⁶; the eigenvalues of BᵀB are wrong by 1 — not a rounding error but the whole number — and the damage is entirely at the small end.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
Fig. 1 Four routes to one spectrum, each singular value’s own relative error, against an answer that has not been rounded.

Why the measurement was not available

A computed σ_min of 3.7·10⁻⁴¹ cannot be checked against another computed σ_min of 3.7·10⁻⁴¹ from a method with the same weakness. Comparing two routes says which disagrees; it does not say which is wrong, and where both are wrong in the same direction it says nothing at all.

This site’s usual answer to that is exact ground truth — the Hilbert matrix, whose inverse has a closed form in rationals, so an error is measured rather than estimated. That works for a solution and not for a singular value, because a singular value is a root of a polynomial and roots are not rational.

The way through is the reason every statement of this subject is about bidiagonal matrices.

The exact route

A bidiagonal matrix has 2n − 1 entries, and they are the data. BᵀB is tridiagonal with entries that are simple products of them. And the number of eigenvalues of a symmetric tridiagonal below a shift s is a count of sign changes in a recurrence with one division per row and no square roots:

q₀ = d₀² − s q_i = d_i² + e_{i−1}² − s − (d_{i−1}e_{i−1})² ⁄ q_{i−1}

The number of negative q’s is the number of eigenvalues below s. That is Sylvester’s law of inertia used as an algorithm, it is what every bisection eigensolver runs, and because it contains only addition, multiplication and division it can be run in BigInt rationals with nothing rounded anywhere.

Bisect on s and each singular value comes out to any relative accuracy asked for. The implementation here brackets each σ² between consecutive powers of two by an exponent search, then takes sixty dyadic bisections inside the bracket — which gives sixty significant bits of the value’s own size, which is what relative accuracy means.

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. 2 The habit in the form this site opened with, on a solution rather than a spectrum.

And the matrix has to be the same matrix

The entries of every test matrix here are powers of two times small integers, and that is not tidiness. If an entry were rounded on its way into the float routine, the matrix the float routine sees and the matrix the exact routine reasons about would be different matrices, and the difference between them would be indistinguishable from the error being measured.

So the construction is d_i = 2^(−g·i) and e_i = ¾·2^(−g·i) for an integer g. Exactly representable, exactly convertible, and graded over as many decades as the number of rows and the choice of g allow.

The exact route is checked before it is used, by two routes of its own. On a diagonal matrix — where the singular values are the entries — it returns 1, 9.54·10⁻⁷, 8.67·10⁻¹⁹ and 7.89·10⁻³¹, each to fourteen digits, without being told they are the answer. And the product of the values it returns is the determinant of the bidiagonal, which is the product of its diagonal: −36.1236 in logarithms, to eight digits, by a route that shares no arithmetic with the bisection.

log₁₀|det H_n| by three routes, to n = 26The exact rational determinant of the Hilbert matrix falls to 10^-387 at n = 26. The accumulated logarithm of the pivots follows it to n = 13 and then departs, reaching 10^-289 — wrong by 98 orders of magnitude and still an ordinary-looking number. The product of the pivots is a separate curve that ends at n = 26, where it underflows to exactly zero.261014182226-400-300-200-1000nlog₁₀|det H_n|exact rationalΣ log|u_ii|Π u_iithree routes, one theoremexact at the largest n-387accumulated logarithm-289decades of disagreement98smallest pivot at that n3.1·10⁻¹⁸every pivot is a normal numberat every size on this axis
Fig. 3 The second of those checks in the setting where this collection first used it: a determinant computed two ways, and where the two stop agreeing.

The four routes

One-sided Jacobi, which is what this site uses. Rotate pairs of columns until they are mutually orthogonal; the column norms are then the singular values. Every operation is a rotation, and a rotation of two columns does not mix their scales in a way that loses the smaller one.

The eigenvalues of BᵀB, which is what the definition suggests. σ² are the eigenvalues of BᵀB, so form it and diagonalise.

The implicit zero-shift QR sweep, which is Demmel and Kahan’s and is what LAPACK runs when relative accuracy is asked for. The property that makes it work is visible in its code rather than argued for: every new entry is an old entry times a cosine or a sine, so there is no subtraction anywhere in the sweep and no cancellation is possible.

And the same sweep with a Wilkinson shift, which is the ordinary method.

What they do on a graded matrix

At 29.5 decades of grading, with σ_min = 2.10·10⁻³⁰:

  • one-sided Jacobi: 4.4·10⁻¹⁶ worst relative error over the whole spectrum;
  • zero-shift QR: 2.2·10⁻¹⁶;
  • shifted QR: 8.9·10⁻¹⁶;
  • eigenvalues of BᵀB: 1.00.

That last number is not a large error. It is a relative error of one, which means the routine returned zero. The information was gone before the eigensolver saw the matrix: σ² of 10⁻⁶⁰ sits inside a sum whose largest term is 1, and 10⁻⁶⁰ added to 1 is 1.

Relative error of two algebraically identical expressions for (1 − cos x)/x², in binary64A log–log plot of relative error against x. The expression written as it reads loses accuracy as x falls and is eventually wrong in every digit; the rearranged form stays at rounding level.10⁻¹²10⁻¹⁰10⁻⁸10⁻⁶10⁻⁴10⁻²110⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹xrelative error of the computed value(1 − cos x)/x², as written2 sin²(x/2)/x²no digits left at allbinary64 throughoutone function, two spellings · zero below 1.5·10⁻⁸
Fig. 4 The mechanism at its ordinary size, from the arithmetic field: what a subtraction of two nearly equal numbers leaves, and how much of it is the answer.
Least squares by QR and by the normal equations in binary32Relative error of the computed coefficients against epsilon, on log axes. The QR route is a flat line near the bottom; the normal-equations route climbs and then stops, at the epsilon where the cross-product matrix becomes exactly singular.10⁻⁸10⁻⁷10⁻⁶10⁻⁵10⁻⁴10⁻³10⁻²10⁻¹10⁻⁹10⁻⁷10⁻⁵10⁻³10⁻¹10¹ε in Läuchli's matrix (smaller ε, larger κ)relative error in the coefficientsAᵀA exactly singularnormal equationsQRκ from 1.7·10⁸ to 17the cliff is at √u = 2.4·10⁻⁴
Fig. 5 And the same object one step earlier in the argument, from the least-squares field: forming AᵀA squares the condition number, and this is what the squaring does to the small end.

The damage is entirely at the small end, which is the shape worth carrying: the six largest values come back to 10⁻¹⁶ from every route, so a check on the largest singular value — the norm, the thing most likely to be looked at — reports that everything is fine.

And how it opens up

The worst relative error over the whole spectrum, against how many decades the matrix is graded over, n = 8Each point is one bidiagonal matrix and the worst relative error any of its singular values suffers, measured against a Sturm bisection in exact rationals. The route through BᵀB is as good as anything at 2.1 decades — 3.63·10⁻¹⁴ — and by 11 decades it is at 14.5, which is not an error in the answer, it is the answer. One-sided Jacobi and the zero-shift bidiagonal sweep are flat at about the unit roundoff across the whole range, which is the claim this site's own SVD has been making in a source comment since it was written and had never measured.0102030405010⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹10⁴10⁷decades of gradingworst relative errorthe answer is gonevia BᵀBone-sided Jacobizero-shift QRone axis, four routesBᵀB at the narrowest grading3.6·10⁻¹⁴and at the widest1.9·10⁷Jacobi, worst over the sweep1.5·10⁻¹⁵zero shift, worst1.1·10⁻¹⁵the definition is not a methodand squaring buries what it squares
Fig. 6 The worst relative error of each route, as the grading is opened from two decades to fifty. Two flat lines, and one that leaves the picture.

At 2.1 decades the BᵀB route is at 3.6·10⁻¹⁴ — as good as anything, which is why this distinction can be gone a career without meeting. By 10.5 decades it is at 14.5. By 19 it is at 577. Past 29 it is returning zeros.

The other two are flat at about the unit roundoff across the whole range, and the flatness is the claim: relative accuracy is not a property that degrades gracefully with the difficulty of the matrix. Either the algorithm has it or it does not, and the algorithms that have it hold it at fifty decades as firmly as at two.

Relative error of each singular value of a bidiagonal matrix graded over 2 decades, by four routesThe matrix is 8×8, bidiagonal, and every entry is a power of two times a small integer, so the matrix the float routines are handed and the matrix the exact routine reasons about are the same matrix. The reference is a Sturm bisection in BigInt rationals, which never rounds. σ_min is 0.00478. One-sided Jacobi holds every value to 1.2·10⁻¹⁵ relative and the zero-shift bidiagonal sweep to 1.09·10⁻¹⁵; the eigenvalues of BᵀB are wrong by 3.63·10⁻¹⁴ — not a rounding error but the whole number — and the damage is entirely at the small end.1234567810⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹singular value, largest firstrelative errorone-sided Jacobizero-shift QRshifted QReigenvalues of BᵀBagainst a rational bisectionσ_min, exactly0.0048worst, one-sided Jacobi1.2·10⁻¹⁵worst, zero-shift QR1.1·10⁻¹⁵worst, eigenvalues of BᵀB3.6·10⁻¹⁴a relative error is a ratioand the denominator is the answer
Fig. 7 At two decades all four routes agree and the picture is four flat lines at the roundoff.
Relative error of each singular value of a bidiagonal matrix graded over 44 decades, by four routesThe matrix is 8×8, bidiagonal, and every entry is a power of two times a small integer, so the matrix the float routines are handed and the matrix the exact routine reasons about are the same matrix. The reference is a Sturm bisection in BigInt rationals, which never rounds. σ_min is 3.726·10⁻⁴⁵. One-sided Jacobi holds every value to 2.2·10⁻¹⁶ relative and the zero-shift bidiagonal sweep to 2.19·10⁻¹⁶; the eigenvalues of BᵀB are wrong by 2.37·10⁶ — not a rounding error but the whole number — and the damage is entirely at the small end.1234567810⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹singular value, largest firstrelative errorone-sided Jacobizero-shift QRshifted QReigenvalues of BᵀBagainst a rational bisectionσ_min, exactly3.7·10⁻⁴⁵worst, one-sided Jacobi2.2·10⁻¹⁶worst, zero-shift QR2.2·10⁻¹⁶worst, eigenvalues of BᵀB2.4·10⁶a relative error is a ratioand the denominator is the answer
Fig. 8 And at forty-five, where three of them are still there.

The route nobody defends, and why it is everywhere

The BᵀB route is the one every account warns against, and it is also the one most code contains, because it is what “compute the singular values” means to anybody starting from the definition.

It appears as the covariance matrix in a statistics package. It appears as the Gram matrix in a kernel method. It appears wherever somebody needed AᵀA for another reason and then took its eigenvalues because they were there. In none of those cases is anybody thinking of themselves as choosing a singular value algorithm — the matrix was formed for a different purpose and the spectrum was taken of what was to hand.

The measurement above prices that. At a narrow grading it costs nothing. At ten decades it costs everything below the middle of the spectrum, and the cost arrives without a warning, without a condition estimate, and with the largest singular values still correct to sixteen digits.

The one honest consolation is that it is detectable from inside. A σ_min returned as exactly zero from a matrix that is not exactly rank-deficient is a signal, and so is a σ_min that changes by orders when the matrix is scaled. Neither is checked anywhere by default.

Why it matters here specifically

A relative error is a ratio and the denominator is the answer, so the whole question is whether anybody cares about the small singular values as values rather than as a distance from zero.

On this site they are cared about constantly.

A rank decision is a decision about a gap, and a gap between σ_k and σ_{k+1} is a ratio. If the small values are known only to absolute accuracy then every gap below the noise floor is indistinguishable from every other, and the decision is being made on numbers that carry no information.

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. 9 The decision this collection has already drawn, with the noise on a slider and the gap shrinking in proportion.
The diagonal of a pivoted R against the singular values, on a 40×40 Kahan matrix at c = 0.3Two curves down the index. The diagonal of R decays at one constant rate to 0.159; the singular values follow it and then fall away to 1.08·10⁻⁵ at the last one. Column pivoting made no interchange at all, because every trailing column norm is equal at every step — the rule had nothing to choose between.1611162126313610⁻⁶10⁻⁵10⁻⁴10⁻³10⁻²10⁻¹110¹10²indexmagnitude|r_kk|σ_kone factorisation, two verdicts‖AP − QR‖/‖A‖2.3·10⁻¹⁶|r_nn|0.16σ_min1.1·10⁻⁵column interchanges0|r_nn| is never below σ_minso the cheap verdict errs one way only
Fig. 10 And the cheaper verdict that a pivoted factorisation gives instead, which errs in the flattering direction.

A condition number is a ratio of two singular values, and the one in the denominator is the small one. A κ computed from a σ_min that is accurate to 10⁻¹⁶ absolute and 10⁰ relative is a κ that is wrong by whatever factor σ_min is wrong by.

How much a perturbation of the right-hand side is amplified, κ = 10⁸The cumulative distribution of the amplification factor over two hundred random perturbation directions, with the condition number marked as the upper limit.110¹10²10³10⁴10⁵10⁶10⁷10⁸10⁹10¹⁰00.250.50.751amplification of the input perturbationfraction of directions at or belowκ = 10⁸worst found 7.6·10⁷6×6, 200 directionsmedian reaches 0.29 of κ
Fig. 11 What that number is used for once it has been computed, which is why being wrong about it is expensive.
How close Hager's estimate is to the true κ₁, over 200 seeded 10×10 matricesFive bars. The estimate is exactly the true condition number on 88% of the sample and inside ten per cent on 93%; the worst underestimate in the whole sample returns 45% of the truth. The last bar is the matrix built to defeat it, at 7.7%, well below anything the random sample reached.each bar is a percentage — of the sample, or of the true condition numberexactly right87.5%inside 10%93.0%inside a factor of 293.0%worst in the sample, ×10044.6%the constructed matrix, ×1007.7%usually exactexact share0.88worst of the sample0.45the constructed matrix0.077a routine that is right most of the timeand never wrong in the safe direction
Fig. 12 And the cheap estimate of it that most libraries print, whose under-reporting this collection has already measured.

And a low-rank approximation’s error is exactly σ_{k+1}. The Eckart–Young theorem says the best rank-k approximation is wrong by the next singular value, so a truncation is choosing an error and the error is one of these numbers.

Error of the best rank-k approximation to a 12×12 matrixApproximation error against k on a logarithmic axis, with the measured error and the next singular value drawn as separate curves lying exactly on top of one another.123456789101110⁻¹⁰10⁻⁸10⁻⁶10⁻⁴10⁻²1rank k of the approximation‖A − Aₖ‖measured, 2-normσₖ₊₁, from theorymeasured, Frobeniusthe first two agreeto 2·10⁻⁸worst |‖A−Aₖ‖₂ − σₖ₊₁| / σₖ₊₁2·10⁻⁸worst Frobenius discrepancy2·10⁻⁸κ = 10⁹; 30 random rank-3 matrices, none closerthe error is σₖ₊₁
Fig. 13 The theorem drawn, from the essay that measured it: the approximation error and the singular value it is equal to.
What a rank-10 approximation can achieve, by spectrumA semi-logarithmic plot of singular value against index for three spectra — geometric decay, algebraic decay, and flat — with the rank-ten approximation error marked on each.1112131415110⁻³10⁻²10⁻¹1index jσⱼσ11geometric, 0.85ʲalgebraic, j^−0.70flatmeasured, and equal to σ₁₁rank-10 error, geometric0.2rank-10 error, algebraic0.19rank-10 error, flat160×60, spectrum chosen rather than the entriesthe matrix decides, not the method
Fig. 14 And the spectra where truncation is worth anything at all.

The cost of the good routes

Neither of the two accurate routes is free, and it is worth saying what each costs before recommending either.

One-sided Jacobi is the slow one. It sweeps over all n(n−1)/2 column pairs and repeats until the columns are orthogonal, which on a general matrix takes several sweeps; the standard bidiagonal method is a reduction plus a QR iteration and is several times faster at any size worth timing. The trade is the plainest in this essay: an algorithm that is a small multiple slower and gets an answer the fast one does not have.

It also parallelises unusually well, because independent column pairs can be rotated at once — which is why it keeps reappearing whenever the hardware changes shape, and why the fastest known high-accuracy singular value codes are Jacobi variants rather than in spite of being them.

The zero-shift sweep is the slow one in a different way. A shift is what makes a QR iteration converge quickly, and removing it means more sweeps for the same deflation. On the graded family here that cost does not appear — the zero-shift sweep converges in three sweeps against the shifted one’s seven, because on a strongly graded matrix the shift is tiny and buys nothing while the extra rotations cost — but on an ordinary matrix it is real, and it is why LAPACK switches between the two on a criterion rather than always taking the safe one.

That criterion is worth noting for what the next essay is about: the switch is made on a measurement of the matrix, not on a preference, and the measurement is close to the classifier the next essay ends on.

Words moved against the block size, n = 96, M = 100A U-shaped curve on logarithmic axes. Too small a block refactors the panel too often; too large a block does not fit in the fast memory and the tiled update evicts what it is about to read. The minimum is at b = 8, which is √M − 2. The three-line count says √(M/3) = 5, which has the right scaling and the wrong constant.110¹10²10⁵block size bwords movedthe count: √(M/3) = 5measured best: b = 8words movedat the best block1.2·10⁵at b = 15.5·10⁵at b = 243.2·10⁵derived from M with no measurement, and scannedthe two agree
Fig. 15 The reason a slower algorithm can be the faster one, from the cost field: what a computation actually spends its time on stopped being the arithmetic some decades ago.
Four factorisations of one 512×12 matrix, κ = 7151A ranked bar chart of the departure from orthogonality, on a logarithmic scale. Classical Gram–Schmidt loses 4.6·10⁻¹⁰, modified 1.5·10⁻¹², Householder 2·10⁻¹⁴. The reduction tree, which touches the matrix once instead of 12 times, is of the same size as Householder's — and its error does not grow with the depth of the tree.classical Gram–Schmidt4.62·10⁻¹⁰modified Gram–Schmidt1.49·10⁻¹²Householder, one sweep2.03·10⁻¹⁴reduction tree, 16 leaves1.48·10⁻¹⁵departure from orthogonality, logarithmicthe tree, at four depths‖AᵀA − RᵀR‖/‖AᵀA‖, depth 13.4·10⁻¹⁵‖AᵀA − RᵀR‖/‖AᵀA‖, depth 24.3·10⁻¹⁵‖AᵀA − RᵀR‖/‖AᵀA‖, depth 31.7·10⁻¹⁵‖AᵀA − RᵀR‖/‖AᵀA‖, depth 41.5·10⁻¹⁵the same algebra, four timestwo of them are products of reflections
Fig. 16 And the shape that argument takes when the parallelism is what is being bought.

What the two good routes have in common

It is worth naming, because it is the only part of this that generalises.

Neither of them ever forms a quantity in which a small singular value has to survive beside a large one.

One-sided Jacobi works on the columns of A directly. The rotation that orthogonalises columns p and q mixes those two columns and nothing else, so a column whose norm is 10⁻³⁰ is never added to a column whose norm is 1 unless the rotation angle says to, and the angle is computed from the ratio rather than the difference.

The zero-shift sweep works on the bidiagonal entries directly. Every entry it produces is a previous entry multiplied by a sine or a cosine — never a difference — so the relative accuracy of every entry is preserved through the sweep, and the relative accuracy of the entries is what the relative accuracy of the singular values rests on.

The BᵀB route forms sums of squares. The shifted sweep forms d² − μ. Both contain a subtraction of quantities that can be nearly equal, and a subtraction of nearly equal quantities is where digits go.

Relative error of three summation algorithms in binary32A log–log plot of relative error against the number of terms for naive, pairwise and compensated summation, each measured against the exactly rounded sum.10¹10²10³10⁴10⁵10⁶10⁻⁹10⁻⁷10⁻⁵10⁻³10⁻¹number of terms addedrelative error against the exact sumin orderin a treecompensatedbinary32 · terms are 1/icompensated: 3·10⁻⁸
Fig. 17 The general statement of it in the field that owns it: the order of a sum decides its value, because addition is where information about small terms is lost.

One more consequence, and it is about scaling

If the small singular values are computed to relative accuracy then a diagonal scaling of the matrix does not damage them, because a scaling multiplies each row by a number and relative accuracy is what survives multiplication.

That connects this essay to the scaling field’s two, and it sharpens both. The condition number of a matrix is a choice of units; the componentwise condition number is not, because the diagonal factors cancel entry by entry. What is measured here is the algorithmic half of the same statement: an algorithm has relative accuracy exactly when it is insensitive to the row scaling that a change of units performs.

So “computes small singular values to relative accuracy” and “is unaffected by a change of units” are the same property said twice. One-sided Jacobi has both. The route through BᵀB has neither — scale one row by 10¹⁵ and the σ that route returns for the rest change, which is a fact about the algorithm and not about the matrix.

That gives a second, cheaper test for the property, which does not need an exact reference at all: scale and re-run. A method with relative accuracy returns the same answers to within a rounding; one without does not. It is the two-routes habit again, with the second route being the first one on a matrix that has been multiplied by a diagonal.

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 1.9·10¹⁰ 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.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
Fig. 18 The units doing that work, from the scaling field: the same matrix and the same answer, with a condition number that moves by decades.
Two perturbation bounds and the error that was measured, on a 8×8 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 4·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)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
Fig. 19 And the condition number that does not move, which is what the property looks like from the perturbation side rather than the algorithm side.

The refusal

The assertion is fed a graded family asked for at less than a bit of grading per row.

It is the case where the whole experiment is about nothing: a “graded” matrix whose entries all sit within a factor of two of each other is not graded, every route gets every value right, and a figure drawn on it shows four flat lines and demonstrates that the distinction does not exist.

Refusing it rather than drawing it is the same discipline as refusing a κ sweep that never leaves one decade. A comparison needs a range over which the thing being compared actually differs, and a comparison drawn where it does not is a comparison that has quietly become a control.

The worst relative error over the whole spectrum, against how many decades the matrix is graded over, n = 6Each point is one bidiagonal matrix and the worst relative error any of its singular values suffers, measured against a Sturm bisection in exact rationals. The route through BᵀB is as good as anything at 1.5 decades — 6.79·10⁻¹⁵ — and by 11 decades it is at 55.4, which is not an error in the answer, it is the answer. One-sided Jacobi and the zero-shift bidiagonal sweep are flat at about the unit roundoff across the whole range, which is the claim this site's own SVD has been making in a source comment since it was written and had never measured.0102030405010⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹10⁴10⁷decades of gradingworst relative errorthe answer is gonevia BᵀBone-sided Jacobizero-shift QRone axis, four routesBᵀB at the narrowest grading6.8·10⁻¹⁵and at the widest1Jacobi, worst over the sweep1.1·10⁻¹⁵zero shift, worst10⁻¹⁵the definition is not a methodand squaring buries what it squares
Fig. 20 At six rows the same two flat lines and the same departure, with the crossing a little earlier because each row carries more of the grading.
Relative error of each singular value of a bidiagonal matrix graded over 11 decades, by four routesThe matrix is 8×8, bidiagonal, and every entry is a power of two times a small integer, so the matrix the float routines are handed and the matrix the exact routine reasons about are the same matrix. The reference is a Sturm bisection in BigInt rationals, which never rounds. σ_min is 1.934·10⁻¹¹. One-sided Jacobi holds every value to 7.1·10⁻¹⁶ relative and the zero-shift bidiagonal sweep to 4.41·10⁻¹⁶; the eigenvalues of BᵀB are wrong by 14.5 — not a rounding error but the whole number — and the damage is entirely at the small end.1234567810⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹singular value, largest firstrelative errorone-sided Jacobizero-shift QRshifted QReigenvalues of BᵀBagainst a rational bisectionσ_min, exactly1.9·10⁻¹¹worst, one-sided Jacobi7.1·10⁻¹⁶worst, zero-shift QR4.4·10⁻¹⁶worst, eigenvalues of BᵀB15a relative error is a ratioand the denominator is the answer
Fig. 21 And at ten decades, which is where the route through BᵀB has already lost the small end.
The worst relative error over the whole spectrum, against how many decades the matrix is graded over, n = 12Each point is one bidiagonal matrix and the worst relative error any of its singular values suffers, measured against a Sturm bisection in exact rationals. The route through BᵀB is as good as anything at 3.3 decades — 9.14·10⁻¹³ — and by 10 decades it is at 2.02, which is not an error in the answer, it is the answer. One-sided Jacobi and the zero-shift bidiagonal sweep are flat at about the unit roundoff across the whole range, which is the claim this site's own SVD has been making in a source comment since it was written and had never measured.0102030405010⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹10⁴10⁷decades of gradingworst relative errorthe answer is gonevia BᵀBone-sided Jacobizero-shift QRone axis, four routesBᵀB at the narrowest grading9.1·10⁻¹³and at the widest1Jacobi, worst over the sweep1.5·10⁻¹⁵zero shift, worst1.8·10⁻¹⁵the definition is not a methodand squaring buries what it squares
Fig. 22 And at twelve rows, where each row carries less of the grading and the crossing moves right.

What is next

Everything above says the comment on this site’s own routine is true. The next essay is about a second family of matrices, on which every ordering in it reverses — where one-sided Jacobi loses the small singular value, the sweep with the theorem behind it does not converge at all, and the shifted sweep that the theorem is a warning about gets everything. Nothing there contradicts the theory. What it contradicts is the way the theory is usually remembered.

Where this goes next

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 Rank is a decision A floating-point matrix does not have a rank. It has a spectrum of singular values, and somewhere in that spectrum is a place where the values stop being signal and start being noise. Deciding where is a judgement, and the evidence for it is a gap. 10⁻⁸10⁻⁷10⁻⁶10⁻⁵10⁻⁴10⁻³10⁻²10⁻¹10⁻⁹10⁻⁷10⁻⁵10⁻³10⁻¹10¹ε in Läuchli's matrix (smaller ε, larger κ)relative error in the coefficientsAᵀA exactly singularnormal equationsQRκ from 1.7·10⁸ to 17the cliff is at √u = 2.4·10⁻⁴ The road that squares the problem The normal equations are the first method every course teaches and the method no library uses. Forming AᵀA squares the condition number, and below ε = √u it does not degrade — it produces a matrix that is exactly singular, from data that was perfectly usable. 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 An answer that is known Almost every demonstration of numerical error estimates the error by computing the same thing more carefully. The Hilbert matrix does not need that: its inverse is a closed form in integers, so the true answer is available exactly and the error is measured rather than approximated.

Reads more easily once this is understood

Essays that name this one as worth reading first.

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 matrixCancellationCondition squaringExact arithmeticGraded matrixJacobi's eigenvalue methodRelative accuracySingular valuesSturm sequence