Small compared to what

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.

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ᵢ = dᵢ² + eᵢ₋₁² − s − (dᵢ₋₁eᵢ₋₁)² ⁄ qᵢ₋₁

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.

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ᵢ = 2^(−g·i) and eᵢ = ¾·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.

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 — the same construction two Gram–Schmidts is about, ordered so that nothing is subtracted. 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 — which is the road that squares the problem taken for a spectrum instead of a solution.

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.

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. It is also the end the best approximation there is computes a truncation error from, so a truncation is unaffected and a rank decision is not. That is a small residual is not a small error with a spectrum in place of a residual, and the same lesson about which quantity a diagnostic reaches.

And how it opens up #

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.

Flat in the grading, and not in the size #

Before the qualification, the flatness itself is worth walking rather than asserting, because it is the whole comparison.

Ten decades of grading and the three good routes do not move at all. Across 2.1, 4.2, 8.4, 10.5 and 12.6 decades one-sided Jacobi reports 1.2·10⁻¹⁵, 1.48·10⁻¹⁵, 1.01·10⁻¹⁵, 7.1·10⁻¹⁶ and 8.77·10⁻¹⁶ — no trend, and the largest value at the easiest matrix. Zero-shift QR runs 1.09·10⁻¹⁵ down to 6.5·10⁻¹⁶ and shifted QR 1.77·10⁻¹⁵ down to 6.58·10⁻¹⁶, both likewise drifting slightly better as the problem gets harder.

Read that sequence carefully, because it is not monotone and the non-monotonicity is the interesting part. A relative error of 14.5 means the computed value is fifteen times the true one; a relative error of exactly 1 means it is zero. BᵀB does not recover between 10.5 decades and 12.6 — it stops returning a wrong number and starts returning no number, which is the better of the two failures and still a failure.

The fourth route covers fourteen orders of magnitude over the same range: 3.63·10⁻¹⁴, 3.04·10⁻¹¹, 0.0205, 14.5, 1. That is the comparison in one line — three routes whose error is a property of the arithmetic and one whose error is a property of the matrix — and it is why “small compared to what” is the question rather than a quibble. A relative error of 1 is not a bad answer; it is the absence of one, arrived at by a method that never divided by anything small and never warned.

That sentence is about one axis, and it invites the reading that the constant is absolute. Holding the grading at thirty decades and varying n instead:

Both good routes rise by about a factor of four over a factor of five in n — near enough linear, which is the n-dependent constant the relative-accuracy theorems carry and which the grading axis cannot show because it is held fixed there. So the guarantee is n·u rather than u, and the flatness belongs to the grading.

That is a refinement rather than a retraction, and it is worth having in the right form: what the sweep establishes is that the difficulty of the matrix does not enter, and what it cannot establish is that the size does not. On a matrix of ten thousand rows the good routes would be at 10⁻¹² rather than 10⁻¹⁶ and the BᵀB route would still be returning zeros, so the comparison is unaffected and the absolute claim is not available.

And the shift is the wrong way round #

The fourth route is introduced as “the ordinary method” and never priced, which leaves out the more useful half of what it does. Sweep counts and rotation counts at n = 12:

The shift helps on an ungraded matrix and hurts on a graded one, and the crossover is at about five decades. At thirty the zero-shift sweep is four times cheaper in rotations and four times more accurate.

The reason is that grading is very nearly deflation. The trailing entries of a graded bidiagonal are already negligible; a zero-shift sweep deflates them at once, and a Wilkinson shift computed from those same tiny entries is a poor shift, which delays what would otherwise have happened for free.

That matters more than a factor of four, because a cost argument is the usual reason given for not running the relative-accuracy method. On this family the relative-accuracy method is also the cheaper one, and the ordinary method is dominated on both counts at every grading the essay is about. The situation is the reverse of the one symmetry is worth more than precision describes for Jacobi against tridiagonal QR, where the accurate method genuinely costs more; here the accurate choice is a one-line change that also saves work.

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 σₖ and σₖ₊₁ 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.

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.

And a low-rank approximation’s error is exactly σₖ₊₁. 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.

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.

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.

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.

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.

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.