The units the matrix is measured in

A linear system carries units. The first row of Ax = b might be a force balance in newtons and the second a displacement in metres; the first unknown might be a length and the third an angle. None of that is written down anywhere in the array of numbers a solver receives, and all of it is present in the numbers themselves.

Multiplying one equation through by a thousand changes nothing about the problem. The same x satisfies the scaled equation — exactly, in exact arithmetic, and in floating point too, because both sides are scaled together and the solution is never touched. A reader who has followed the condition number is an amplifier would expect the sensitivity of the answer to be unchanged as well.

It is. The condition number is not.

What the picture is of #

The matrix on the left of that figure is an ordinary one: random entries, a strengthened diagonal, a condition number of about 3.4. It is as far from a hard problem as a matrix gets.

Every point to the right of it is that same matrix with its rows multiplied by powers of ten, spread evenly over the number of decades on the axis. The right-hand side is scaled with them, so the solution is untouched. The right-hand end of the axis, at ten decades, is a matrix whose first row is in units 10¹⁰ times larger than its last.

κ_∞ at that end is 1.9·10¹⁰. It was 9.8.

The solution is the same. It is the same because the site’s habit of building problems with a known answer applies here as everywhere else — b was formed from an x of small integers, so the true answer is an input rather than an estimate. At eight decades of spread, where κ₂ has gone from 3.4 to 1.3·10⁸, the forward error of the scaled solve is 7.3·10⁻¹⁵ — not 10⁻⁸, which is what κ = 10⁸ would lead a reader to budget for.

Two definitions of one word #

The reason is that there are two ways to say how much can a small change in the data move the answer, and they differ in what “small” means.

Normwise. Perturb A by a matrix E with ‖E‖ ≤ ε‖A‖. That is one number describing the whole matrix, and κ(A) is exactly the amplification factor from it to the answer.

Componentwise. Perturb each entry a_ij by at most ε|a_ij| — every entry by a small fraction of itself. That is what rounding the data does, what storing it in a narrower format does, and what a measurement error usually looks like.

On a matrix whose rows are all the same size the two are close together. On a badly scaled matrix they are not, and the reason is worth saying plainly: a normwise perturbation of a badly scaled matrix perturbs the tiny row by the size of the big one. The row whose entries are 10⁻⁵ is being moved by ε·10⁵. Of course the answer falls apart; the perturbation obliterated an equation. It is a stupid model of what a machine does to a matrix, and κ is the exact amplification factor for it.

The componentwise number — Skeel’s, written cond(A) = ‖ |A⁻¹| |A| ‖_∞, with the absolute values taken entry by entry before the product — answers the second question, and it does not notice the units at all.

The invariance is exact #

Not approximate. Write A = D·A₀ with D diagonal and positive. Then A⁻¹ = A₀⁻¹D⁻¹, so

|A⁻¹| |A|  =  |A₀⁻¹| |D⁻¹| |D| |A₀|  =  |A₀⁻¹| |A₀|

because |D⁻¹||D| is the identity, entry by entry, before any norm is taken. There is nothing to estimate: the two diagonal factors cancel inside the absolute values, and the componentwise condition number of the scaled matrix and the unscaled one are the same number.

assertSkeelIsInvariantUnderRowScaling checks it at 10⁻⁸ relative rather than exactly, for the one reason worth stating: the two numbers are computed by inverting two different matrices, and the inversions round. The algebra is exact and the arithmetic is not, which is the sentence this site is named for, arriving in its own library.

It is also checked against a second scaling — random, twelve decades wide, unrelated to the one the test matrix was built with — because an invariance verified only on the construction that produced it is an invariance verified against itself.

Equilibration, and what a theorem does not say #

If κ is a fact about the units, the obvious repair is to fix the units: divide each row by its own largest entry, so every row has ∞-norm exactly 1. That is equilibration, and it is what LAPACK’s xGEEQU computes and xGESVX applies.

Van der Sluis’s theorem says the result is within a factor of n of the best any row scaling can achieve. That is a bound, and a bound of n is compatible with equilibration being fairly bad in practice — a factor of eight on an 8×8 is not nothing.

So it is measured. assertEquilibrationIsNearlyOptimal runs 1,500 random diagonal scalings over the same range, then a coordinate refinement from the best of them, and asks how much better than equilibration the search can get. The answer is 1.30, against a bound of 6. The theorem is loose and the practice is close to optimal, and the second sentence is the one a reader can use.

And the matrix it does nothing for #

Here is the trap the previous section sets, and it is worth walking into deliberately.

Scale the rows and the ill-conditioning goes away is false. It is false for the matrix this site has used since its first commit.

The badly scaled matrix has κ_∞ and cond seven orders of magnitude apart. The Hilbert matrix has them within three, and both enormous. That gap between the two numbers is the diagnostic: a matrix whose normwise condition number is far above its componentwise one is badly described, and a matrix where the two agree is badly conditioned. The first is repairable by a diagonal matrix and the second is not repairable at all.

An answer that is known measured the Hilbert system’s solution against the exact rational answer and found 0.049 relative error at n = 13. No scaling recovers those digits, and the reason is now expressible in one number.

What the bounds actually predict #

The two condition numbers are not decorations. Each is the constant in a perturbation bound, and the bounds can be run against a measurement.

The perturbation is entrywise and relative: each entry moved by up to ε of itself, which is what storing the matrix in a shorter format does. At ε = 10⁻⁹ the worst error over forty draws is 2.1·10⁻⁹. The componentwise bound is 4.8·10⁻⁹. The normwise bound is 0.19.

Both are bounds. Only one of them is a prediction.

Running the same figure at zero spread is not a formality. It is the check that the argument is about the scaling and not about something else in the construction — the site’s own rule, learned the hard way in an earlier phase, that a test problem must differ in the property under test and nothing else.

Which number a library reports #

Every library reports the normwise one, and reports an estimate of it rather than the number itself. numpy.linalg.cond, MATLAB’s cond and rcond, LAPACK’s xGECON: all normwise, and an estimate that can be fooled is about how the estimate is arrived at.

Almost nothing reports the componentwise number, because computing it needs |A⁻¹| — the entries of the inverse, not merely its action — and forming an inverse costs more than the solve did. LAPACK’s expert drivers do compute a componentwise error bound by a cheaper route, and the fact that they are called expert drivers is a reasonable summary of how often anybody reaches for them.

So the number a reader is most likely to have in hand is the one that is a property of the units, and the number that would tell them what to expect is the one nobody computes.

What this does to the site’s own identity #

The spine of this site is

forward error  ⪅  condition number × backward error

and it has been used to separate two authors of a wrong answer: the algorithm, through the backward error, and the problem, through the condition number.

This essay puts a third author in the middle. The condition number is not a property of the problem alone; it is a property of the problem and of how it was written down. The identity still holds — it is an inequality about norms and nothing here disturbs it — but the term that was supposed to carry “how hard this problem is” turns out to carry “how hard this problem is, in these units”, and the second is a decision somebody made before the solver was called.

The identity has a componentwise form which does not have that problem:

forward error  ⪅  cond(A, x) × componentwise backward error

and the figure above is that inequality, measured, with the ordinary form drawn above it for scale.

Where badly scaled matrices come from #

It would be comfortable to file all of this under badly written input, and it is worth resisting, because the matrices that arrive in this state arrive there honestly.

A physical model with more than one quantity in it. A structural analysis couples displacements in metres to rotations in radians and forces in newtons; a circuit couples volts to amps to coulombs. The stiffness matrix that comes out has blocks differing by whatever the constants of the physics differ by, and nobody chose those constants.

A constrained problem. A Lagrange multiplier row is a constraint written in whatever units the constraint is natural in, next to rows written in the units of the objective. The penalty method’s whole idea is to multiply one block by a large number, which is a row scaling performed on purpose, and the large number is the parameter being swept.

Finite differences on a graded mesh. A discretisation whose cells vary in size by three orders of magnitude produces rows scaled by the cell size to some power, and the grading is there because the solution varies fast somewhere and slowly elsewhere.

And time. An implicit time step of Δt multiplies one block by Δt and leaves another alone, so the condition number a solver reports changes as the step is refined even when the underlying operator does not.

In every one of those the scaling is a consequence of the modelling rather than a mistake in typing. That is why the repair is a routine in every serious library rather than advice, and why a κ read off an unscaled assembled matrix is so often a number about the assembly.

What a solver can do about it, and when #

A library that equilibrates does three things and they are worth separating, because only the first is free.

Compute D. One pass over the matrix, n numbers out, and — in the version LAPACK uses — every multiplier rounded to a power of two, so the scaled matrix is exactly the unscaled one with exponents adjusted and the repair introduces no rounding of its own. equilibrateRowsByPowersOfTwo does that here for the same reason. A scaling that costs accuracy to apply is a scaling that has to be argued for; this one does not.

Decide whether to apply it. LAPACK’s rule is a threshold on the ratio of the largest row norm to the smallest, and the default is not to scale unless the ratio is extreme. That is a decision made for the caller, silently, and it means two runs of the same driver on two matrices can differ in whether an equilibration happened at all.

Report which condition number. Having scaled, the expert driver reports the condition number of the scaled matrix, which is the honest one for the solve it performed and is not the condition number of the matrix the caller handed it. Two libraries can therefore report condition numbers orders of magnitude apart for the same input and both be right about a different question.

None of that is hidden — it is in the documentation of every routine involved. What it is, is easy to read past, and the number that survives into a user’s notebook is a single scalar with no record of which of those three decisions produced it.

Column scaling is a different question, and a worse one #

Everything above is about scaling rows, which is scaling equations. Scaling columns is scaling unknowns — measuring a length in millimetres instead of kilometres — and it is not benign in the same way.

A row scaling leaves x alone. A column scaling changes x, because A(Dy) = b means y = D⁻¹x, and the error in y is not the error in x. So “which x is accurate” becomes a question about which units the answer is wanted in, and there is no scaling-invariant statement available: an answer accurate in metres can be inaccurate in kilometres if the two components differ enormously in size.

That is why the equilibration a library performs is a row scaling, and why the column scaling it offers has to be asked for and comes with a warning. The asymmetry is not an oversight; the two operations do different things to the object being solved for.

What is worth carrying #

A condition number is a number about a description. It is exactly the amplification factor for the normwise perturbation model, that model is a poor description of what a machine does to a badly scaled matrix, and the whole of the discrepancy in this essay follows from those two sentences.

The gap between κ and cond is the diagnostic, and it is cheap to reason about even when it is expensive to compute. If the rows of a matrix differ enormously in size, expect κ to be reporting that fact and not much else. If they do not, expect κ to be telling the truth.

And a repair that works has to be checked against a case where it does not. The Hilbert curve in the hero figure is there because the previous three paragraphs are the kind of argument that generalises itself in a reader’s memory into a rule, and the rule is wrong.

The next essay takes the componentwise number seriously as an object rather than as a foil: a condition number scaling cannot move. The one after it takes the same question into the elimination itself, where the pivot rule turns out to be reading the units too — the pivot that reads the units.