An estimate that can be fooled

The condition number is the most quoted quantity in this subject and almost nobody computes it.

κ(A) = ‖A‖·‖A⁻¹‖ needs the inverse. Forming an inverse costs three matrix multiplications’ worth of work on top of the factorisation, which is more than the solve the caller actually wanted. So no library does it. What numpy.linalg.cond returns on its 1-norm branch, what MATLAB’s rcond returns, what LAPACK’s xGECON returns, and what every “check the condition number” instruction in every user guide is pointing at, is an estimate: ‖A‖₁, which is exact and free, multiplied by an estimate of ‖A⁻¹‖₁ obtained from four or five solves with factors that already exist.

This essay is that estimator, the guarantee it comes with, and the matrix it is wrong about.

What Hager’s algorithm does #

‖B‖₁ is the largest absolute column sum, so it is maxⱼ ‖B eⱼ‖₁ over the n unit vectors. Trying all n costs n products, which is the whole matrix, which defeats the purpose.

Hager’s observation is that x ↦ ‖Bx‖₁ is convex, and a convex function on the unit ball of the 1-norm attains its maximum at a vertex — and the vertices are ±eⱼ. So the maximum is reachable, and a gradient method that walks from vertex to vertex will reach a local maximum in a handful of steps.

The walk is four lines:

1. Start at x = e/n, the centre.
2. Compute y = Bx and take ξ = sign(y).
3. Compute z = Bᵀξ. If maxₖ |zₖ| ≤ zᵀx, stop — no vertex looks better than the current one.
4. Otherwise move to x = eⱼ for j = argmax |zₖ|, and repeat.

Each step is one product with B and one with Bᵀ. In the real routine those are triangular solves with the LU factors and their transposes, so a product with A⁻¹ costs O(n²) against the O(n³) already spent. Four or five steps, and the estimate is the largest ‖Bx‖₁ seen along the way.

And the size of the matrix is the size of the lie, which the slider says in three columns at once.

At n = 4, 8, 12, 16 and 24 the estimate reads 4.00, 8.00, 12.00, 16.00 and 24.00 — it returns the dimension of the matrix, exactly, at every size — while the true 1-norm reads 10.64, 48.64, 114.00, 206.72 and 474.24. The ratio is 0.3759, 0.1645, 0.1053, 0.0774 and 0.0506, and multiplied by n those are 1.50, 1.32, 1.26, 1.24 and 1.21.

So the estimate is about 1.2/n of the truth, and the underestimate has no bound. A twentieth at twenty-four, a fiftieth at sixty, a hundredth at a hundred and twenty. Whatever tolerance a caller sets on rcond, there is a matrix of this shape large enough to slip under it.

The last two columns of that readout are the ones that make it silent. The walk visits exactly one column and spends exactly five products with the matrix at every size on the slider. Its own stopping test fires: no vertex looks better than the one it is standing on, which is true, because the construction makes the neighbouring gradient point back at the column it started from. So the routine does not run short of budget, does not hit an iteration cap, and does not report anything unusual. It converges, in the sense its own code means by the word, and returns n.

The guarantee, and its direction #

‖Bx‖₁ for a specific x with ‖x‖₁ = 1 is a lower bound on ‖B‖₁. The estimate is a maximum over the x’s visited, so it is a maximum of lower bounds, so it is a lower bound.

The estimator can never return more than the truth.

That is the dangerous direction, and it is worth being explicit about because “it is a lower bound” sounds like a safety property. An estimate that was sometimes too large would occasionally warn about a matrix that was fine — an annoyance. This one is sometimes too small, and too small means reporting that a matrix is better conditioned than it is. The failure mode is a clean bill of health.

assertTheEstimateIsAlwaysALowerBound checks the one-sidedness over four unrelated families — Gaussian matrices, inverses of Gaussian matrices, the badly scaled matrix from the units the matrix is measured in, and Hilbert — because a structural property checked on one family is a property of that family.

And it is usually exact #

This is the half that explains why the routine is in every library, and leaving it out would make the essay an attack on a tool that works.

Over 300 seeded 8×8 matrices, measured through the route a library takes — estimating ‖A⁻¹‖₁ from solves rather than running the walk on A itself — four out of five come out exactly right, and the worst underestimate in the whole sample returns 38% of the truth. The constructed matrix above sits well below anything random sampling reached.

Two details of that measurement are worth stating because getting either wrong would misreport the routine by a factor of two.

It is measured on A⁻¹, not on A. The walk run directly on a Gaussian matrix is exact 36% of the time; run on its inverse — which is what the library does — it is exact 82%. The reason is that the columns of an inverse have far more spread in their 1-norms than a Gaussian’s do, so there is a clear winner for the walk to walk to. Quoting the first number as though it were the second would understate the routine by half.

And the two shares are different questions. 82% exactly right and 88% within ten per cent: the gap between those is the cases where the walk found a good column and not the best one.

The matrix it is wrong about #

Now the other half. The walk stops when no other column looks better from where it is standing, and “looks better” means |sign(Beⱼ)ᵀ B eₖ| — the inner product of column k with a vector of ±1s.

A column can hide from that test by having a large 1-norm and a zero inner product with the sign vector. For a sign vector of all ones, that means a column whose entries alternate and cancel.

So the construction is three sentences:

• Column 1 is all ones. Its 1-norm is n. It is where the walk stops.
• Column 2 alternates ±t. Its 1-norm is n·t — the true answer — and its entries sum to exactly zero, so the dual vector never points at it.
• Every other column is 0.9 throughout, plus a small bump on its own diagonal so the matrix is not singular. That keeps their sums below column 1’s, so the walk prefers column 1, and keeps every row sum positive, so the opening probe’s signs are all +1.

The last clause is the only constraint on t: the row sums are 1 ± t + 0.9(n − 2), so t can be anything below 1 + 0.9(n − 2).

t therefore grows linearly with n, and the estimate is exactly 1/t of the truth.

Three things that are not the explanation #

It is not an iteration limit. assertItStopsBecauseItIsSatisfied runs the walk at a cap of forty instead of five and gets the same number after the same two iterations. The stopping test fires. The walk is standing on a local maximum of a convex function on a polytope, which is a perfectly respectable place to stop, and the global maximum is one vertex away and invisible from there.

It is not LAPACK’s extra probe. The real routine adds one more estimate from a vector with alternating signs and increasing magnitude — xᵢ = (−1)^(i+1)(1 + (i−1)/(n−1)), scaled by 2/(3n) — and it is there precisely to catch cancellations the unit-vector probes miss. On this matrix it returns less than the walk did and is discarded. The assertion says so.

And it is not an exotic matrix. Its entries take four distinct values: 1, ±t, 0.9, and 0.9 plus a bump. Written out at n = 8 it fits in eight lines. The refusal a counterexample dismissed as a contrived matrix is fed the claim that it has many distinct entries, and fails.

This is not the first estimator with this problem #

Hager’s is the second one LINPACK and LAPACK have shipped, and the first had the same shape of hole, found the same way.

The Cline–Moler–Stewart–Wilkinson estimator solves Aᵀy = d for a right-hand side d whose signs are chosen greedily to make y large, then solves Az = y, and reports ‖z‖/‖y‖. It is a witness search too: it constructs a vector it hopes lies along a small singular direction. It is cheap, it is usually good, and counterexamples were published for it within a few years.

That history is the reason the theorem here is stated the way it is. Higham’s 1988 result is not “Hager’s estimator can be bad”; it is that for any algorithm which examines a fixed number of vectors, there are matrices it is arbitrarily wrong about — and the proof is a counting argument rather than an attack on a particular rule. A method that looks at k vectors has k pieces of information about n columns, and for k < n an adversary has room to hide one.

So the construction in this essay is an instance of something that could not have been avoided by choosing a better rule. What a better rule buys is that the adversary has to work harder, and the adversary is not usually present.

What the block version changes, and what it does not #

LAPACK’s shipped routine is dlacn2, which is Hager’s walk with a block of t probe vectors instead of one — Higham and Tisseur’s block generalisation. Instead of standing on one vertex it stands on t of them, keeps a record of which columns it has already visited so it cannot cycle, and stops when a whole block fails to improve.

That is strictly better and it is better in the way the counting argument predicts: it examines t times as many vectors, so the matrix that defeats it has to hide a column from t sign vectors at once rather than from one. The default t is 2.

It does not change the direction of the error, because nothing can: every vector examined produces a lower bound and the estimate is their maximum. And it does not close the hole, because the counting argument does not care how large a fixed t is — it cares that t is fixed and n is not.

What it changes is the size of matrix on which the construction here starts to bite, and that is a real improvement and is not the same as a guarantee.

What it costs, which is the reason for all of this #

The accounting is worth writing out, because the estimator exists for exactly one reason and it is not subtle.

A factorisation of an n×n matrix costs about 2n³/3 operations. Forming the inverse from it costs about 4n³/3 more — twice the factorisation, and three times it in total. At n = 1000 that is the difference between one solve and four, for a diagnostic.

The estimator costs five triangular solves, at 2n² each: about 10n². At n = 1000 that is 10⁷ operations against 1.3·10⁹ for the inverse, a factor of 130, and the factor grows linearly with n.

So there is no version of this where the honest number is computed by default. The choice is between an estimate and nothing, and an estimate is better than nothing in the same careful sense as everything else in this essay: a small one is informative and a large one is not.

What a reader can actually do #

Four things, in increasing order of cost.

Look at the matrix. If the rows differ enormously in size, κ is reporting the units and the units the matrix is measured in is the essay for that. This costs one pass and is the most common case by a distance.

Read the estimate as one-sided. A reported κ of 10¹² means the matrix is at least that bad. A reported κ of 10² means nothing was found, which is not the same as nothing being there.

Ask for the componentwise error bound. LAPACK’s expert drivers return one, computed by Skeel’s refinement argument rather than by forming an inverse, and it answers the question a caller usually meant to ask.

Or check the answer. A residual costs one product. It does not bound the forward error — the whole of a small residual is not a small error is that it does not — but combined with an estimate that is a lower bound on κ it gives a one-sided statement in the useful direction: large residual times large κ-estimate is definitely bad news, and a solver that reports both has told a reader more than one that reports either.

What a user reads #

The number a library hands back is ‖A‖₁ · est, and everything above was about the second factor. So the last step is to put the construction where a caller meets it — as the inverse of the matrix being solved, since that is what the estimator is applied to.

At n = 16 the true κ₁ is 3.5·10³ and the reported one is 2.7·10². A factor of 12.9, which is exactly t.

A matrix whose true condition number is 10⁴ is reported at 10³ — the gap the condition number is an amplifier prices in digits. Both of those numbers are inside every rule of thumb a reader has been given and past none of them, and the difference between them is three digits of the answer.

Hilbert, by contrast — the matrix an answer that is known is built on — is estimated exactly: 3.3872790725·10¹⁰ by both routes at n = 8, to every digit. The routine is not broken. It has a hole in it, and the hole is shaped like a matrix nobody writes down by accident and somebody could write down on purpose.

The shape this shares with a cheaper certificate #

The estimator’s error is one-sided, and the direction it errs in is the one that flatters the matrix. That is not the only cheap surrogate on this site with that property.

The cheap rank and what it cannot see is about column-pivoted QR standing in for an SVD, and |rₙₙ| ≥ σₘᵢₙ holds for every triangular factor — one line of algebra, no pivoting needed — so that verdict can only report a matrix as further from singular than it is, and never as closer.

Two independent mechanisms, two different quantities, and the same asymmetry. It is worth asking whether that is a coincidence, and the answer is that it is not: both surrogates work by finding a witness. A witness to a large ‖A⁻¹‖ is a vector; a witness to a large ‖R⁻¹‖ is the last diagonal entry. A search that finds a witness proves a lower bound on the badness, and failing to find one proves nothing at all.

The tempting conclusion is that every cheap certificate of ill-conditioning is a lower bound. That is wider than the reason supports, and a third surrogate settles it. Forty random 16×16 matrices:

The first two behave exactly as described. The estimator’s maximum is 1.000 — it is exact or low and never high — and |rₙₙ| never falls below σ_min, so it reports the matrix as further from singular by a median factor of 4.8 and by 171 at worst.

Gershgorin goes the other way, at every one of the forty. Its radius is two to two-and-four-fifths times the true spectral radius and never below it. It is cheap — one pass over the entries, cheaper than either of the others — and it errs entirely towards caution, in the way a bound that is proved requires of anything called a bound.

So cheapness is not what decides the direction. The form of the argument is. Gershgorin is a containment argument: it proves every eigenvalue lies in a union of discs, and a proof of containment can only be conservative, because widening it keeps it true. A witness search is the dual — it proves badness by exhibiting an instance, and failing to find one keeps nothing true.

The useful form of the rule is therefore the one the mechanism gives rather than the one two examples suggested: a surrogate that proves badness under-states it, and a surrogate that proves goodness over-states it. Both are safe, and they are safe about different things. A reader needs to know which kind is in front of them before deciding whether the number they have been handed is a reassurance or a warning — the same reading rank is a decision asks for of a verdict — and the estimator’s failure mode in this essay is exactly what happens when a lower bound on badness is read as an upper bound on it.

What is worth carrying #

The number a library prints is not the condition number. It is a lower bound on it, obtained by maximising a convex function over a handful of vertices, and it is right far more often than a bound suggests.

When it is wrong, it is wrong in the direction that says the matrix is fine. There is no version of this estimator that errs the other way, because the whole method is the exhibition of a witness.

And it can be wrong by any factor. Not through bad luck and not through an iteration limit — through stopping, correctly, at a local maximum, on a matrix with four distinct entries in it.

The practical reading is not “distrust rcond”. It is that a small rcond is evidence and a large one is the absence of evidence, and those are different things.

A cheaper certificate that is one-sided the same way #

The estimator on this page errs in the direction that flatters the matrix because it works by exhibiting a witness. So does the stopping test of a Krylov method for a matrix function, which reports on the subspace rather than on the answer — and so, in a different way, does a determinant, which is not a certificate of anything.