Eigenvalues, singular values, rank

The cheap rank and what it cannot see

Almost nobody computes singular values to decide a rank. The standard substitute is QR with column pivoting, read off the diagonal of R — and there is a triangular matrix on which the greedy rule makes no interchange at all, has no better column available at any step, and reports a matrix eight orders of magnitude further from singular than it is.

Worth reading first: Rank is a decision · The best approximation there is.

Rank is a decision settled what rank means on a machine. There is no rank; there is a list of singular values, a gap somewhere in it with luck, and a threshold somebody chose. The essay’s figure is the spectrum with the threshold drawn across it, and its argument is that the threshold is the content.

What it did not say is that almost nobody computes those singular values.

An SVD costs several times what a factorisation costs — a bidiagonalisation and then an iterative phase, against one sweep of reflections. The standard substitute is QR with column pivoting: eliminate greedily, always taking the column with the largest norm in the part still to be eliminated, and read the rank off |r_kk|. LAPACK’s xGEQP3, MATLAB’s qr(A, 0) with a permutation, scipy.linalg.qr(pivoting=True) — all the same algorithm, and all described as rank-revealing.

The diagonal of a pivoted R against the singular values, on a 50×50 Kahan matrix at c = 0.5Two curves down the index. The diagonal of R decays at one constant rate to 8.69·10⁻⁴; the singular values follow it and then fall away to 3.73·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.171319253137434910⁻¹²10⁻¹⁰10⁻⁸10⁻⁶10⁻⁴10⁻²110²indexmagnitude|r_kk|σ_kone factorisation, two verdicts‖AP − QR‖/‖A‖1.1·10⁻¹⁶|r_nn|8.7·10⁻⁴σ_min3.7·10⁻¹²column interchanges0|r_nn| is never below σ_minso the cheap verdict errs one way only
Fig. 1 The diagonal of a pivoted R against the singular values of the same 50×50 matrix. The diagonal decays smoothly at one constant rate and ends at 8.7·10⁻⁴; the singular values follow it and then fall away to 3.7·10⁻¹². Column pivoting made no interchange at all. Drag to a planted rank instead, where the two curves coincide.
A randomised rank-k solve, 4 seeds a rankRelative error against the rank kept, on a logarithmic vertical axis, with a vertical bar at each rank spanning the seeds. The band is 1.68 wide at rank 8, where the method is at its worst, and 1.004 wide at rank 20, where it is at its best. The same computation on the same data returns a different answer each time.0481216202428321rank keptrelative errormedianthe answer movesspread at rank 81.7spread at rank 201best median error0.15widest where the method is worstand the bound does not say so
Fig. 2 Where this site left the rank question: a band of thresholds, each giving a different answer, with no gap to put one in. Everything here is about which of two spectra a caller is looking at when they choose.

What the diagonal is worth, in one line

There is a guarantee, it is one line of algebra, and it holds for every triangular factor — pivoted, unpivoted, from any algorithm at all.

eᵀ_n R⁻¹ e_n = 1/r_nn, so ‖R⁻¹‖₂ ≥ 1/|r_nn|, so σ_min(R) = 1/‖R⁻¹‖₂ ≤ |r_nn|.

|r_nn| is an upper bound on σ_min. The cheap test can report a matrix as further from singular than it is, and never as closer. Its rank verdict can be too high and never too low.

That is the same asymmetry as an estimate that can be fooled, reached by a completely different mechanism, and the two together are the phase’s clearest instance of a pattern: a cheap certificate of ill-conditioning is always a witness, and a witness gives a one-sided bound.

assertTheDiagonalIsAnUpperBoundOnSigmaMin checks it across forty Gaussian matrices, three planted ranks and three Kahan matrices, with the worst overestimate in the sample at 2.7·10⁴.

And it works

The essay would be dishonest without this half. On a matrix with a genuine gap, the diagonal of R finds it.

The diagonal of a pivoted R against the singular values, on a matrix of planted rank 12Two curves down the index, lying on top of each other. Both fall by more than five orders of magnitude between index 12 and index 13, so both report rank 12 at any threshold placed in the gap.171319253137434910⁻¹³10⁻¹¹10⁻⁹10⁻⁷10⁻⁵10⁻³10⁻¹10¹indexmagnitude|r_kk|σ_kone factorisation, two verdicts‖AP − QR‖/‖A‖9.3·10⁻¹⁶|r_nn|10⁻¹²σ_min10·10⁻¹³column interchanges44|r_nn| is never below σ_minso the cheap verdict errs one way only
Fig. 3 The same comparison on a 50-column matrix built with twelve singular values at order one and the rest at 10⁻¹². The two curves are indistinguishable: both fall by more than eleven orders of magnitude between index 12 and 13, and any threshold in that gap gives rank 12 from either routine.

assertTheCheapTestFindsAPlantedGap runs planted ranks 2 through n − 2 at three shapes, and the two verdicts agree every time, through a gap in the diagonal of at least 10⁵. The generator is greedy, the greedy choice is the right one on a matrix with a real gap, and the whole thing works.

What the counterexample removes is the word guarantee.

The diagonal of a pivoted R against the singular values, on a matrix of planted rank 4Two curves down the index, lying on top of each other. Both fall by more than five orders of magnitude between index 4 and index 5, so both report rank 4 at any threshold placed in the gap.1611162126313610⁻¹³10⁻¹¹10⁻⁹10⁻⁷10⁻⁵10⁻³10⁻¹10¹indexmagnitude|r_kk|σ_kone factorisation, two verdicts‖AP − QR‖/‖A‖1.1·10⁻¹⁵|r_nn|10⁻¹²σ_min10·10⁻¹³column interchanges37|r_nn| is never below σ_minso the cheap verdict errs one way only
Fig. 4 A planted rank of four. Both curves fall eleven orders of magnitude at the same index, the greedy rule interchanges ten times, and the cheap verdict is the right one.

Kahan’s matrix

Let s² + c² = 1 and build the upper triangular matrix with sⁱ down the diagonal and −c·sⁱ everywhere above it:

K = diag(1, s, s², …) · [ 1  −c  −c  ⋯ ]
                        [    1  −c  ⋯ ]
                        [        1  ⋯ ]

Now compute the norm of column j restricted to rows k and below, which is exactly what the pivot rule looks at at step k:

‖K(k:, j)‖²  =  s^{2j} + c²(s^{2k} + s^{2k+2} + ⋯ + s^{2j−2})
             =  s^{2j} + s^{2k} − s^{2j}
             =  s^{2k}

Every trailing column norm is exactly s^{2k}, for every j ≥ k, at every step. They are all equal. The greedy rule is choosing between numbers that are the same number.

So it makes no interchange, R is K unchanged, and the diagonal ends at s^{n−1}. At n = 50 and c = 0.5 that is 8.7·10⁻⁴, and σ_min is 3.7·10⁻¹²: a factor of 2.3·10⁸.

Not a bad choice — no choice

Three things are asserted about that, because each rules out an explanation a reader will reach for.

The rule is not choosing badly. assertKahanDefeatsColumnPivoting checks that the interchange count is zero. There was no better column at any step, because there was no different column at any step.

The factorisation is not inaccurate. ‖AP − QR‖/‖A‖ is at rounding. This is not a case of cancellation destroying the diagonal; the diagonal is exactly what the matrix put there.

And there is no threshold that would have worked. Every consecutive ratio in the diagonal is the same s, to within 10⁻⁴ across the whole of it, so the diagonal has no gap in it anywhere for a threshold to sit in. The refusal a threshold offered as the repair for a diagonal with no gap in it is fed exactly that claim and fails.

The diagonal of a pivoted R against the singular values, on a 50×50 Kahan matrix at c = 0.3Two curves down the index. The diagonal of R decays at one constant rate to 0.0992; the singular values follow it and then fall away to 4.87·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.171319253137434910⁻⁷10⁻⁵10⁻³10⁻¹10¹indexmagnitude|r_kk|σ_kone factorisation, two verdicts‖AP − QR‖/‖A‖2.4·10⁻¹⁶|r_nn|0.099σ_min4.9·10⁻⁷column interchanges0|r_nn| is never below σ_minso the cheap verdict errs one way only
Fig. 5 Kahan’s matrix at c = 0.3 rather than 0.5. The diagonal still decays at one constant rate and the singular values still fall away at the end; what changes is how far.

And it has no floor

A single counterexample is a curiosity. The statement worth making is that the failure is unbounded.

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 2.3·10⁸ at n = 50, 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.81624324048110²10⁴10⁶10⁸size of the matrix|r_nn| ÷ σ_minthe two agreec = 0.2c = 0.35c = 0.5no ceilingratio at n = 1021ratio at n = 502.3·10⁸interchanges, anywhere0the greedy rule never had a choiceand the gap grows with every row
Fig. 6 |r_nn| ÷ σ_min against the size of the matrix, one curve per Kahan parameter. Every curve rises without turning over — the diagonal ends at s^{n−1} and σ_min falls faster — so there is no size at which the cheap test starts working again. Drag the largest size drawn.

At c = 0.5 the ratio runs 1.2·10³, 7.0·10⁴, 4.0·10⁶ and 2.3·10⁸ at n = 20, 30, 40 and 50. Along the other axis, at n = 40, it runs 6.4·10² , 6.5·10⁴ and 4.0·10⁶ as c goes 0.2, 0.35, 0.5. Both directions are asserted to be monotone.

The verdict a caller receives

Ratios become decisions when a threshold is applied, and the decision is where the argument lands.

Put the threshold in the middle of the gap between σ_min and σ_{n−1} — the only place it could go if the singular values were known — and the singular values report rank 49 out of 50 while the diagonal of R reports 50. The two routines disagree about whether the matrix is singular, and neither has made a mistake.

assertTheTwoRoutinesDisagreeAboutTheRank states it in exactly that form: the threshold is checked to be inside the gap, and the two counts are checked to differ.

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. 7 The same sweep truncated at n = 30, where the ratios are still in the hundreds and thousands. Nothing in the shape of these curves suggests a size at which the cheap test recovers.

The second finding: a strict comparison of two equal numbers

The trailing column norms are equal in exact arithmetic. Once computed, they differ by an ulp.

So a pivot rule written with a strict > — which is the obvious way to write it, and the way it is written in every textbook — is comparing two rounding errors and calling the answer a pivot choice.

Where the same pivot rule lands on Kahan's matrix, with and without a tie tolerance, at c = 0.5Two curves of |r_nn| ÷ σ_min against n. Every trailing column norm of this matrix is exactly equal at every step, so a rule that interchanges on any strict improvement is comparing two rounding errors: it makes 0, 9, 14, 6 interchanges at n = 20, 30, 40, 50 and lands 4.5·10⁶ times nearer the truth at the largest size. It lands in the SAME place at the three largest sizes while the other rule's answer grows by orders of magnitude across them, which is what a floor set by rounding looks like rather than a rule that is finding something.18283848110²10⁴10⁶10⁸size of the matrix|r_nn| ÷ σ_minequal treated as equalstrict comparisona choice between equal numbersties found, at n = 5049interchanges, tie-tolerant0interchanges, strict6ratio between the two4.5·10⁶the improvement is real and is not a repairit is the ulps, and they are not reproducible
Fig. 8 The same rule with and without a tie tolerance. The tie-tolerant version finds a tie at every step and makes no interchange. The strict version makes 9, 14 and 6 interchanges at n = 30, 40 and 50 and lands six orders of magnitude nearer the truth. Drag the Kahan parameter.

The strict rule does better. At n = 50 it returns a ratio of 51.6 instead of 2.3·10⁸.

It is not a repair, and the count is the evidence. Fourteen interchanges at n = 40 and six at n = 50 — not monotone in the size, because it is not a function of the size. It is a function of which ulp fell where, so it does not survive a different summation order, a different blocking, or a different precision. LAPACK, which downdates its column norms rather than recomputing them, is a different summation order.

assertTheStrictRuleReordersOnRounding asserts the non-monotonicity directly, and the refusal a rounding-driven interchange count read as a pivot strategy is fed the claim that the count grows with the size.

This is the site’s absolute-tolerance rule turned around. The steady-06 phase found that an absolute epsilon inside a tick loop is a claim about the scale of the data. Here a zero tolerance inside a comparison is a claim that two computed quantities can be meaningfully ordered, and on this matrix they cannot.

Where the same pivot rule lands on Kahan's matrix, with and without a tie tolerance, at c = 0.35Two curves of |r_nn| ÷ σ_min against n. Every trailing column norm of this matrix is exactly equal at every step, so a rule that interchanges on any strict improvement is comparing two rounding errors: it makes 8, 11, 7, 16 interchanges at n = 20, 30, 40, 50 and lands 8.8·10⁵ times nearer the truth at the largest size. It lands in the SAME place at the three largest sizes while the other rule's answer grows by orders of magnitude across them, which is what a floor set by rounding looks like rather than a rule that is finding something.18283848110¹10²10³10⁴10⁵10⁶10⁷size of the matrix|r_nn| ÷ σ_minequal treated as equalstrict comparisona choice between equal numbersties found, at n = 5049interchanges, tie-tolerant0interchanges, strict16ratio between the two8.8·10⁵the improvement is real and is not a repairit is the ulps, and they are not reproducible
Fig. 9 The tie comparison at c = 0.35. Eight, eleven, seven and sixteen interchanges at the four sizes — not monotone, and landing at the same ratio at the three largest, which is a floor set by rounding rather than by anything the matrix contains.

Where the rank verdict is actually used

It would be easy to read all of this as being about a diagnostic printed and forgotten. The verdict is usually an input to something.

A least-squares solve with a rank-deficient design. xGELSY runs column-pivoted QR, decides the rank from the diagonal, and solves the reduced problem. Overstate the rank and the solve inverts a direction that is not there, amplifying noise by the reciprocal of a singular value that should have been discarded. The routine’s own documentation offers xGELSD, which uses the SVD, as the alternative, and the choice between them is a choice between speed and this essay.

A rank-deficient constraint set. An optimisation code that detects redundant constraints does it by factorising the constraint Jacobian and reading off a rank. Overstating it keeps a redundant constraint, which makes the KKT system singular in a way the next factorisation has to deal with — and that is when symmetry is not enough inheriting a decision made here.

Model-order reduction and subset selection. Column-pivoted QR is the standard way to choose k columns of a matrix to keep, and the pivot order is the answer, not a diagnostic. On a matrix where the rule has nothing to choose between, the columns it keeps are the first k, which was not a choice.

In all three, the failure mode is the same shape: the rank comes out too high, so nothing is discarded, so a direction that carries no information is kept and inverted. The consequence appears later, as an answer with a large component along a direction the data never determined — which is where the answer stops being in the data arriving by a route that essay does not mention.

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 6.8·10¹⁰ at n = 64, 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.81930415263110²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 = 646.8·10¹⁰interchanges, anywhere0the greedy rule never had a choiceand the gap grows with every row
Fig. 10 Extended to n = 64. The largest ratio drawn is 10¹⁰: the smallest diagonal entry of a pivoted R reporting a matrix ten orders of magnitude further from singular than it is.

What it costs to be sure

The honest accounting, because “just use the SVD” is the obvious response and it is not free.

A Householder QR of an m×n matrix costs about 2mn² − 2n³/3. Column pivoting adds the norm bookkeeping, which is O(mn) if the norms are downdated and O(mn²) if they are recomputed — the implementation here recomputes, deliberately, so that what is measured is the pivot rule rather than one implementation’s cancellation in the downdate.

An SVD of the same matrix costs the bidiagonalisation, about 4mn² − 4n³/3, plus an iterative phase. Call it two to three times the pivoted QR in practice, and more if singular vectors are wanted.

So the substitution buys a factor of two or three. That is worth having on a large problem and it is not worth having on a small one, and the sentence a reader can use is: if the matrix is small enough that the difference would not be noticed, compute the singular values. The cases where the cheap route matters are exactly the cases where checking it is expensive, which is an uncomfortable place to leave it and is where it is.

Error of the best rank-k approximation to a 40×40 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.12345678910111213141516171819202122232425262728293031323334353637383910⁻¹⁰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. 11 The other reason to want the singular values: the error of a rank-k approximation IS σ_{k+1}, which this site checked as an equality rather than a bound. A rank verdict from a diagonal of R does not carry that number.

What is actually rank-revealing

Strong rank-revealing factorisations exist — Gu and Eisenstat’s is the standard one — and they work by doing what column pivoting does and then checking: after the greedy pass, look for a column swap that would increase |det R₁₁| by more than a factor f, and perform it if there is one. Iterate until none exists.

That gives a genuine bound, σ_k / (|r_kk|) ≤ p(n, k) with p polynomial, and it costs the greedy pass plus however many corrective swaps are found. On most matrices there are none and the cost is the greedy pass; on Kahan’s matrix there are many.

The distinction is exactly the one this essay is about. Column pivoting is a heuristic that usually finds the gap; Gu–Eisenstat is an algorithm that verifies it found one. The word “rank-revealing” is attached to both in common usage, and only the second has a theorem.

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. 12 Where this started: the spectrum with a threshold across it, and the decision that is the content. Everything here is about which of two ways of estimating that spectrum a caller is actually looking at.
How close Hager's estimate is to the true κ₁, over 200 seeded 12×12 matricesFive bars. The estimate is exactly the true condition number on 84% of the sample and inside ten per cent on 90%; the worst underestimate in the whole sample returns 36% 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 right83.5%inside 10%89.5%inside a factor of 289.5%worst in the sample, ×10035.7%the constructed matrix, ×1007.7%usually exactexact share0.83worst of the sample0.36the constructed matrix0.077a routine that is right most of the timeand never wrong in the safe direction
Fig. 13 The other cheap surrogate with a one-sided error. Both work by exhibiting a witness, both prove a lower bound on the trouble, and both fail by finding nothing and reporting that as nothing being there.

What LAPACK’s downdate changes

One implementation detail is worth separating from the argument, because it is the most common objection to this essay’s second finding.

xGEQP3 does not recompute the trailing column norms at every step. It downdates them: having eliminated one entry from a column, it subtracts that entry’s square from the stored norm rather than re-adding the rest. That is O(n) a step instead of O(mn), and it is why the routine is affordable.

It is also a subtraction of two nearly equal quantities, which on a matrix whose columns shrink steadily loses relative accuracy exactly where the accuracy matters. LAPACK guards against it by recomputing whenever the downdated norm falls below a fraction of the original — the standard trick, and one more absolute-versus-relative tolerance of the kind this site has now met four times.

So the ulps this essay’s second finding is about are, in the shipped routine, a different set of ulps: they come from the downdate rather than from the sum. That changes which interchanges get made and does not change the conclusion, which is that the comparison being made is between numbers that are equal, and its outcome is therefore a property of the arithmetic rather than of the matrix.

qrColumnPivoted here recomputes, deliberately, so that what is measured is the pivot rule and not one implementation’s cancellation in its bookkeeping.

The other direction of the bound

|r_nn| ≥ σ_min is the half this essay is about. There is a bound the other way, and knowing what it costs is what makes the first one worth stating.

For column-pivoted QR the classical result is σ_k(A) ≤ |r_kk|·2^{k} at worst — an exponential factor, and Kahan’s matrix is the construction showing it is not merely an artefact of the proof. So the two-sided statement is: |r_nn| is above σ_min and within an exponential factor of it, which brackets nothing useful.

Gu and Eisenstat’s strong rank-revealing factorisation replaces the exponential with a polynomial, at the cost of the corrective swaps described above. That is the whole difference between a heuristic and an algorithm here: not the greedy pass, which both perform, but whether anything checks the result.

What is worth carrying

The diagonal of a pivoted R bounds σ_min from above, always, and the bound can be off by any factor. The direction is structural and the size is not bounded by anything.

Kahan’s matrix does not defeat the rule by being adversarially pivoted. It defeats it by making every choice identical, so the greedy rule has nothing to be greedy about, and no gap appears in the diagonal for a threshold to find.

A strict comparison of two quantities that are equal in exact arithmetic is a comparison of their rounding errors, and an improvement obtained that way is not reproducible across implementations even when it is real.

And “rank-revealing” names two different things. One of them has a theorem.

What links here

Computed from the collection, not written here: the essays that point at this one.

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.

Column pivotingCondition estimationCounterexampleHouseholder reflectionKahan's matrixLower boundNumerical rankQR factorisationRank-revealing QRSingular values