The scaling that buys ten orders

The previous essay ends with a problem that has lost every digit it had, under a change of variable that is exact in both directions. This one is the repair, and the repair is two lines.

γ = √(‖K‖/‖M‖),      δ = 2/(‖K‖ + γ‖C‖).

Replace λ by γμ, multiply the whole polynomial by δ, and solve the quadratic with coefficients (δγ²M, δγC, δK). Every eigenvalue of the original is γ times an eigenvalue of the scaled problem, exactly. Both parameters are computed from three norms and nothing else — no iteration, no tuning, no reference to the answer.

What it recovers, and the shape that matters #

At n = 8, over γ from 1 to 10⁸: the unscaled forward error runs 7.5·10⁻¹⁴ to 1.3·10⁻³, and the scaled one runs 1.3·10⁻¹³ to 1.6·10⁻¹³. Ten orders of magnitude, and — the half that matters more — flat.

Ten orders is the headline and flatness is the finding. Flatness is what says the repair is the right one rather than a lucky one. A device that reduced the error by a constant factor would leave a reader with a slope and a worry: it would be helping, at a rate that runs out. This leaves a horizontal line, and the reason is structural. After scaling, every stop of the sweep is the same problem in the same units, so there is nothing left for γ to do — the sweep has been collapsed to a point rather than tilted. The coefficient norms come out within a factor of a few hundred of each other at every stop, where they spanned 4·10¹⁵ before, and the backward error for the quadratic — the quantity the previous essay shows growing to 1.2·10⁻⁴ — sits at 1.5·10⁻¹⁵ throughout.

Why those two numbers #

γ makes the outer coefficients balance: after the substitution ‖γ²M‖ and ‖K‖ are equal by construction, which is the same idea as the units the matrix is measured in — a row scaling chosen so that no row is systematically larger than another — carried from a matrix to a polynomial. δ then sets the overall size, so that ‖Q̃‖ is about one and no coefficient is near the top or bottom of the format’s range.

Neither is subtle and neither is a heuristic in the sense choosing without knowing uses the word. There is no parameter being guessed against a truth that is not available: the scaling is a change of variable, it is exact, and the only thing being chosen is which exact change of variable to make. If the choice were bad the answer would still be right; it would merely be no better than before.

How wrong γ is allowed to be #

A formula invites the reading that it has to be got right, so it is worth measuring how much it does not. The badly scaled problem at γ₀ = 10⁸, with γ moved off the formula’s own value by k decades:

Six decades of γ land within three orders of the best, so the formula is one member of a wide family of adequate choices rather than a prescription. A code that estimates the three norms crudely, or reuses a γ computed for a neighbouring problem, is in no danger at all. The best γ in the sweep is not the formula’s either — it is a decade below it, 5.6·10⁻¹⁴ against 8.2·10⁻¹⁴ — which is a factor of one and a half and is the right size for the gap between a good rule and an optimum.

The second row is the one worth reading twice. The balance — the ratio of the largest scaled coefficient norm to the smallest, which is exactly the quantity γ is chosen to make small — is symmetric in k, to the digit: 10⁸ at both ±4, 10¹⁶ at both ±8. The error is not. At the same imbalance of 10¹⁶ it is 5.3·10⁻⁸ on one side and 5.6·10⁻³ on the other, five orders apart.

So the quantity the scaling is designed around does not predict the accuracy it produces. That is worth stating flatly because it undercuts an obvious-looking improvement: a code that chose γ by minimising the coefficient balance numerically, rather than by the closed form, would be minimising a function that is symmetric about a point the error is not symmetric about, and would have no way to know which side of it to err on.

Which side to err on is answerable, and the answer is that it is not a fact about the units. Run the same sweep through the trailing reduction and the basin is even-handed — 4.5·10⁻⁹ against 7.2·10⁻⁹ at ±4, 0.45 against 2.1 at ±8. The lopsidedness belongs to the leading reduction, which inverts a block holding γ²M beside a fixed identity, so moving γ up and moving it down are not the same computation however symmetric the norms look. The same measurement gives the two reductions’ other difference for free: the leading one is the more forgiving away from the optimum and the less accurate at it, 8.2·10⁻¹⁴ against 9.1·10⁻¹⁵.

What it does not do, which is where the interesting part is #

γ = √(‖K‖/‖M‖) is a single number, and it balances the polynomial as a whole. It says nothing about where the eigenvalues are, and the natural next question is whether the same three norms can say that too.

They can, half of the time, and the half they fail on is instructive.

The tropical roots are the corners of the max-plus polynomial with coefficients ‖M‖, ‖C‖ and ‖K‖: if ‖C‖² > ‖M‖‖K‖ there are two of them, at ‖K‖/‖C‖ and ‖C‖/‖M‖, and otherwise a double one at √(‖K‖/‖M‖), which is γ. They are advertised as estimates of the moduli of the two groups of eigenvalues — obtainable before anything runs, from three cheap numbers.

Measured on the overdamped chain, at four sizes: the large tropical root divided by the largest modulus is 0.961, 0.948, 0.944, 0.943 at n = 4, 8, 16 and 32. The small tropical root divided by the smallest modulus is 5.36, 17.08, 60.83, 229.6 over the same four sizes — a factor growing like n².

The reason is structural rather than a defect of the bound, and it is worth stating as a general sentence:

A norm is a maximum. It is exact about the end of a spectrum that is a maximum and says nothing about the end that is a minimum.

The largest eigenvalue modulus of this family is governed by ‖C‖/‖M‖, which is a ratio of norms; the smallest is governed by the smallest eigenvalue of K, which is 4sin²(π/2(n+1)) and falls like 1/n². No norm of K contains that number. Nothing about the estimate is being sloppy — it is answering the question it can answer.

This site has said something close to it once before from the other side. Small compared to what is about relative accuracy in the small eigenvalues of a symmetric matrix, and the whole difficulty there is the same: quantities that describe a matrix as a whole are quantities about its top end.

One γ, and when one is not enough #

The tropical roots are not only a prediction; they are also the argument for a second scaling. When the max-plus polynomial has two distinct corners the eigenvalues fall into two groups whose moduli are widely separated, and a single γ cannot put both of them near one. The device the literature offers is to scale twice — once with each tropical root — and to keep, from each run, the eigenvalues that group predicts.

On this family it is not needed, and the reason is worth reading because it is the same reason the prediction half-fails. The two tropical roots here are ‖K‖/‖C‖ and ‖C‖/‖M‖, and their ratio is ‖C‖²/(‖M‖‖K‖), which for the overdamped chain is about a hundred. A hundred is not a separation that hurts: the coefficients after one scaling come out within a factor of a few hundred of each other, and the backward error sits at the rounding level for every eigenvalue, small and large alike, at every stop of the sweep. Nothing in the measurement asks for a second pass.

What would ask for one is a problem whose two groups are twenty decades apart, which is exactly the shape of a palindromic quadratic — and there the repair turns out to be neither one scaling nor two, but a division. That the same symptom has three different repairs in three different families is not a defect of the theory; it is what it looks like when a phenomenon is being described from several sides.

Reading the numbers at the far end #

It is worth walking one stop of the sweep slowly, because every quantity in the argument is visible in it.

At γ = 10⁸ the coefficients are (10¹⁶M, 10⁸C, K) with norms 2.8·10¹⁶, 2.0·10⁹ and 6.8. The eigenvalues are the chain’s, divided by 10⁸, so they run from about 2·10⁻¹⁰ to 7·10⁻⁸. Now look at the scale in the backward error’s denominator: |λ|²‖M‖ is about 10⁻¹⁶ × 10¹⁶ ≈ 1, |λ|‖C‖ is about 10⁻⁸ × 10⁹ ≈ 10, and ‖K‖ is 6.8. All three terms are of order one — which is exactly right, and exactly the problem. The residual ‖Q(λ)x‖ is being compared against a quantity of order one, while the entries being cancelled to produce it are of order 10¹⁶. Sixteen digits of cancellation, in a format that has sixteen digits.

That account is about the size of the coefficients and says nothing about n, so it predicts the same recovery at every chain length. It is worth checking, because the cancellation it describes is the one thing an elimination accumulates along.

The unscaled forward error grows with the chain and the scaled one does not. Across n = 4, 6, 10 and 12 the unscaled answer is wrong by 4.33·10⁻⁴, 6.74·10⁻³, 7.41·10⁻² and 0.396 — three orders of magnitude of degradation — while its backward error stays between 7·10⁻¹⁵ and 3·10⁻¹³ and reports nothing.

The scaled column across the same four sizes is 4.77·10⁻¹⁵, 1.02·10⁻¹⁴, 2.72·10⁻¹⁴, 3.02·10⁻¹⁴: it grows by a factor of six while the unscaled one grows by nine hundred. So the recovery is not merely large at one size, it removes the size dependence — which is the stronger claim and the one the mechanism above actually implies, since a cancellation between numbers of order one has nothing to accumulate.

Scaling replaces those coefficients with three matrices of norm about one, and the same cancellation is between numbers of order one. Nothing else about the problem changed: same eigenvalues up to a factor, same conditioning, same solver, same arithmetic. That is the whole mechanism, and it is why the recovery is complete rather than partial.

The two lines are cheap, and the alternative is not #

The cost of the repair is three Frobenius norms, a square root and a division: O(n²) work in front of an O(n³) factorisation, which is free in any accounting anybody uses. It is worth stating because the natural alternative — solve, notice the answer is bad, solve again in different units — costs a factor of two and requires noticing, and the previous essay is a demonstration that nothing available at run time does the noticing.

There is a third option that a reader might reach for: do the arithmetic in higher precision. It works, and it is worth pricing against this. Doubling the significand buys sixteen more digits and costs a factor of anywhere from eight to a hundred in time depending on the hardware, which is the measurement buying the accuracy back makes in the opposite direction. Scaling buys ten of those digits for three norms. Where the two meet is a real engineering question and it has an obvious answer.

The ordering, which is the practical finding #

Two knobs have now been identified: which units and which linearisation. A reader who has to fix one first is entitled to be told which.

The measurement says the units, by a long way. At γ = 10⁶, the six ways of solving one problem — three linearisations, each reduced two ways — give forward errors from 5.5·10⁻⁶ to 2.2·10⁻⁴, a spread of 41. Scaling the same problem takes every one of them to about 10⁻¹³. So the choice of linearisation is worth a factor of tens and the choice of units is worth ten orders, and the sensible order to fix them in is not a matter of taste.

That ordering is the sort of thing this site’s combination field exists to find — four knobs and one floor is the same shape of result for a different set of parameters — and it arrives here for free, because both knobs were measured against the same closed form.

The scaling is not optional in a narrow format #

There is a second reason to scale, and in a narrow arithmetic it is the only one that matters: at some point γ²M stops being representable.

The boundary is √(largest finite value ÷ largest entry of M), which is a closed form computed from the format’s own parameters. In binary64 it is 1.3·10¹⁵⁴ and nobody meets it by accident. In binary32 it is 1.8·10¹⁹, and in fp16 it is 256.

Past that boundary the entries are infinities, the linearisation contains them, and there is no answer rather than a poor one. The scaled coefficients are representable at every stop, because δγ²M has norm about one whatever γ was. So in a format with a short exponent the two lines are not an accuracy device — they are what makes the problem exist, which is the other half of a format being the half that decides.

The trap inside that is worth its own sentence. The norm the scaling is computed from can overflow before the scaling that would prevent the overflow has been computed. A naive Frobenius norm sums squares, so it fails at √(largest finite value) even when the answer is comfortably representable — measured in binary32, the naive route returns an infinity from 10¹⁹ upwards while the true norms run to 5.8·10²² and are perfectly representable. float.js’s scaled norm, which is the repair a norm that overflows before it is a norm is about, returns them correctly at every stop. The same trap, one level up, and with a consequence rather than an inconvenience: the quantity that overflows is γ’s own input.

What a code should actually do #

Collecting the measurements into a procedure, because the field is at the point where one can be written down.

Compute the three norms first, with a scaled norm rather than a naive one. They are the diagnostic and they are γ’s input, and the naive route fails before the problem does.

Scale, always. There is no case in which it hurts: when the coefficients are already balanced γ comes out near one and δ near a constant, and the substitution is the identity to within a rounding. A caller who scales unconditionally pays three norms and can stop thinking about units.

Then choose the linearisation, and it matters less. The first companion form reduced by inverting the leading coefficient is the default for a reason; the trailing reduction is worth having when the eigenvalues near zero are the ones wanted, because it maps them to the large end where a Schur factorisation resolves them best.

And report η for the quadratic rather than the residual for the matrix. One matrix–vector product with each coefficient, three norms already computed. It is the only number in the whole computation that answers the caller’s question.

That last item is the one this field would change about existing software if it could change one thing. Everything above it is a repair; that one is a measurement, and the measurements are what this site is for.

What is left after the scaling #

The honest boundary of this essay is that scaling fixes the part of the problem that is about units and none of the part that is about conditioning. If a quadratic has a genuinely ill conditioned eigenvalue — one whose κ is 10¹⁰ because Q′(λ) is nearly singular along the eigenvector, not because the coefficients are badly balanced — then γ and δ change nothing about it. The condition number that the previous essay measures at 4.98 is invariant under the substitution in both directions: that is what made the loss attributable to the method, and it is equally what makes the repair unable to help when the loss is not.

Which is the same division of labour this site has drawn from its first essay onwards, arriving in a field where the third party — the reformulation — is new. Scaling removes the reformulation’s contribution exactly. The problem’s own difficulty is still the problem’s.

Why this is a change of variable and not a preconditioner #

One last distinction, because the words are close enough to cause trouble.

A preconditioner replaces a problem by a different problem whose answer is related to the original’s, and the relation costs something to undo: changing the condition number on purpose solves M⁻¹Ax = M⁻¹b and the answer x is the same, but the residual is not, the iteration count is not, and the theory has to be redone for the pair rather than for A.

This is not that. The scaled polynomial’s eigenvalues are the original’s divided by γ, exactly, and the eigenvectors are unchanged — not approximately, not up to a rounding, but as an algebraic identity in one line. Nothing has to be undone except a multiplication, and no property has to be re-derived. That is what makes it worth doing unconditionally: a preconditioner is a bet, and a change of variable is not.