Bracketing an error nobody can measure

The essay that showed the ADI error is the rational function rather than a bound on it put one quantity on the page five times, and every appearance of it was ‖P − ZZᵀ‖/‖P‖ — the relative distance between the low-rank factor the iteration produces and the Gramian it approximates. That number is the whole verdict on a run. It is also the one number a caller of the method never has.

It appears on those pages because P appears on them. The models the figures are drawn on are constructed with a known eigenbasis, so the Gramian is a formula rather than a solve: the right-hand side is rotated into the eigencoordinates, each entry is divided by the sum of the two eigenvalues it sits between, and the result is rotated back. A hundred-thousand-state model offers none of that. Its Gramian has 10¹⁰ entries, and not forming it is the entire reason the low-rank iteration was written. The measurement that decides whether a run worked requires the object the run exists to avoid.

So the useful question is not how accurate the iteration is. That was settled, and settled as an identity rather than an inequality. The question is what a caller holds instead, and the answer is that there are exactly two candidates, arriving at opposite ends of the run. Before a single solve there is max|r|², the squared value of the rational function the shifts amount to, computable from the two ends of the spectrum of −A and from nothing else. After the last solve there is the Lyapunov residual ‖AZZᵀ + ZZᵀAᵀ + BBᵀ‖/‖BBᵀ‖, computable from the factor alone without an n × n matrix ever existing.

What follows measures both against the truth, on models where the truth is available, and reports how far apart they sit. They do not merely both approximate the error: one is above it and one is below it at every measurement taken, so the pair is a bracket rather than two estimates.

The truth on these pages is a formula, not a better computation #

The habit that makes the comparison possible is the same one an answer that is known is built on: an error is measured against a closed form, never against a second and more careful run of the same machinery. A Gramian has a closed form whenever A’s eigendecomposition is known, because the Lyapunov equation diagonalises — in the eigenbasis every entry decouples and the solution is an entrywise division. It costs one eigendecomposition of A and two matrix products, and it is exact in the sense that matters: nothing about the iteration under test enters it.

That is why the ratio in the figure above is a measurement and not a plausibility argument. It is also why the ratio stops being available at exactly the size where anybody would want it. An eigendecomposition is O(n³) and dense, which is the cost an equation whose unknown is a matrix already found unpayable at a hundred states when the equation was written out as a linear system.

The same restriction hides a second quantity, and it is the one the whole field is built on. The decay of the Gramian’s own eigenvalues — λₖ₊₁(P)/λ₁(P), the thing that decides whether truncation is worth doing at all — is a statement about the spectrum of P, so it needs P.

Both the achieved ratio and the bound are on that picture, and only one of them survives contact with a large problem. The bound is two numbers and a formula. The ratio required forming a twenty-six by twenty-six Gramian and eigendecomposing it, and at a hundred thousand states neither step exists. Everything this field measures about accuracy is measured on models small enough not to need the field. That is not a criticism of the models; it is the reason the two substitutes below have to be priced.

One of the two substitutes is known before a single solve #

max|r|² costs an estimate of the extreme eigenvalues of A and an evaluation of a product of k Möbius factors on a logarithmic grid. It is available before the first shifted solve, which makes it the natural thing to plan against, and it is the quantity where to put the poles of a rational function minimises when it argues for geometric clustering.

Measured against the truth, across five grid sizes and every value of k the figures draw — forty-four pairs in all — the ratio of the error to the prediction never leaves the interval 0.7025 to 0.7915. That is a spread of 1.1267, over a range in which the error itself moves from 3.4·10⁻¹ to 1.8·10⁻⁵, over four orders. Written as a rule: the error is 0.746 times the prediction, with every one of the forty-four measurements inside ±6.1% of that.

The lowest values of the ratio are at the coarsest grid and the highest at the finest, but the ordering is not clean and the variation is small either way. What does move cleanly is the error itself, and it moves in the direction that decides how many shifts to buy.

Nine shifts are nine solves at either size, and they buy a relative error of 1.811·10⁻⁵ on the ten-state model and 1.676·10⁻³ on the forty-state one — ninety-three times worse for identical work. Refining the grid by four widens κ from 48.37 to 680.62, and the rational function has to be small over a wider interval with the same number of poles. A caller who fixed a shift count on a coarse mesh and kept it while refining is buying steadily less of a Gramian, and the prediction says so before any of the solves are run, which is the same accounting the accuracy worth paying for makes for a hierarchical solve.

The other is known afterwards, and needs nothing that has n² of anything #

The Lyapunov residual substitutes ZZᵀ into the equation the Gramian solves and measures what is left over, relative to the right-hand side. It is the natural stopping quantity, and it is what a tolerance that reads its own residual sets an inner accuracy from in a different field.

Written naively it forms ZZᵀ, which is n × n and defeats the purpose. It does not have to be. With G = [AZ | Z | B], which is n by 2k + m, the residual is G M Gᵀ for a fixed small M carrying two identity blocks off the diagonal and one on it, and its Frobenius norm squared is tr(MKMK) with K = GᵀG. Every object in that expression is either n by 2k + m or small and square. The residual of an n × n equation is computable in the width of the factor, which is the same move an iterate that must be made smaller makes for a low-rank iteration in the tensor field: never form the object, work in the factors.

The two routes were run against each other on all forty-four pairs. They agree to 8·10⁻⁸ or better on every model of sixteen states or more, and to 1.7·10⁻⁶ at worst, that worst case being the last step of the smallest model, where the residual is 6.89·10⁻⁶ and a trace of a product of small matrices is cancelling most of its digits. The arithmetic is worth stating in counts. At forty states with ten shifts, G holds 840 numbers and its Gram matrix 441, against 1,600 for the Gramian — a saving of nothing much. At a hundred thousand states the same two counts are 2.1·10⁶ and 10¹⁰, a factor of 4,762, and that is where the distinction stops being pedantic.

So the residual is cheap, universal and available. What it is not is the error, and this site has an essay about precisely that: a small residual is not a small error is the standing warning that substituting an answer back and finding it fits verifies the wrong thing. The measurement here is the same warning with the factor filled in.

The error is trapped between them at every size measured #

Across the forty-four pairs the residual is below the error every time, between 0.2702 and 0.6858 of it. The error is below the prediction every time, at 0.7025 to 0.7915 of it. Neither ordering has a single exception, and neither is an accident of the family: on three further models — a coarser grid, a system with a low-degree transfer function, and a modal system with twelve geometrically spaced poles — the residual stays between 0.223 and 0.645 of the error at every step from k = 2 to k = 12.

Reading the two together at a fixed number of shifts gives the width of the interval. At nine shifts:

The last row is the bracket. The error is known to lie inside a factor of 3.74 at ten states and 4.64 at forty, and the interval widens by a quarter across a fourfold refinement of the grid, in the direction that matters least: it is widest where the problem is largest. Both ends of it are computable at a scale where the quantity between them is not, which is the whole content of the result — an unmeasurable number is not an unbounded one.

The sharper estimate is the one that arrives first, and it is the one that needs calibrating #

Two numbers bracketing a third invites the question of which is the better single estimate, and the answer is lopsided. Substituting 0.746 times the prediction puts every one of the forty-four measurements inside ±6.1%. Substituting 2.32 times the residual puts them inside a factor of 1.59 — the residual-to-error ratio spans 2.54 where the error-to-prediction ratio spans 1.13. The number available before the first solve is an estimate nearly ten times tighter than the number available after the last.

That is not the ranking a reader expects, and there is a real cost hidden inside it. The 0.746 is not a constant of the construction. It is how much of P lies in the direction where |r| attains its maximum, which depends on B and on the eigenvectors of A — neither of which enters the Zolotarev argument at all. Measured on a modal system whose poles sit away from the peaks of r, the same ratio settles at 0.436 rather than 0.75, so a caller carrying the heat model’s constant to that system would overestimate its error by seventy per cent.

Obtaining the constant means running the iteration once on a problem where the error can be measured, which is a problem where P is available, which is the situation this essay opened by ruling out. So the sharper of the two estimates is sharp only after a calibration that the large problem cannot supply, and the looser one needs no calibration whatever. That is the honest division of labour between them, and it is the same division why a Gramian can be truncated at all arrived at from the other side: the interval predicts the rate and the system supplies the constant.

The lower half drifts because it is the error seen through the operator #

The upper half of the bracket is flat in k to three per cent. The lower half is not, and the way it moves is worth knowing before anybody stops a run on it. On the forty-state model the residual is 0.686 of the error at k = 2 and 0.270 at k = 10 — it falls by a factor of 2.5 while the error falls by 427. The residual gets more optimistic as the iteration converges.

The reason is an identity rather than a tendency. The Gramian satisfies AP + PAᵀ + BBᵀ = 0 exactly, so subtracting that from the residual’s definition leaves R = −(AE + EAᵀ) with E = P − ZZᵀ. The residual is not a quantity of its own: it is the error mapped through the Lyapunov operator. Computed both ways the two agree to 2.7·10⁻¹¹ relative at worst over the sizes checked.

That factors the ratio cleanly. The residual over the error is the operator’s amplification of the error direction, ‖AE + EAᵀ‖/‖E‖, times ‖P‖/‖BBᵀ‖. At forty states the second factor is 2.081·10⁻³ and fixed, while the first falls from 329.5 at k = 2 to 129.9 at k = 10; their product is 0.2702, which is the measured ratio. The amplification is free to sit anywhere between 2a = 19.73 and 2b = 13,428, and it sits at 130 — well up the spectrum, and drifting downwards as the shifts clear the decades one at a time. The two factors also move against each other with size: between ten states and forty the amplification at the last step rises from 50.3 to 129.9 while ‖P‖/‖BBᵀ‖ falls from 7.565·10⁻³ to 2.081·10⁻³, which is why the ratio only drifts from 0.381 to 0.270 instead of moving by orders.

This is the condition number is an amplifier arriving with a specific operator in place of a general one, and it says what to do about it: a residual read as an error should be divided by an amplification, and on this family that amplification is a two-digit number rather than a κ. What it must not be read as is an error directly, and the drift is why. A run stopped at a residual target has a forward error between one and a half and four times the number it stopped on, and which of those it is depends on how far through the run it stopped.

A bracket four wide cannot decide a contest worth two per cent #

An interval of a factor of 4.6 is enough to size a run and not enough to compare two of them, and the difference between those two uses is worth measuring rather than asserting. The comparison this field actually makes between reduction methods is the one interpolating at the model’s own poles reports, and the margins in it are small.

Under two per cent in one norm and under twenty-five in the other, and both are decided every time — the orderings never swap. A quantity known only to a factor of 4.6 decides neither. It could not even establish the sign of the H∞ margin, let alone its size, and it would report two models differing by 1.94% as indistinguishable.

The two situations are not the same question asked at different precisions. The transfer-function norms in that contest are computable on the models being compared, by sampling the difference system, and they were computed. ‖P − ZZᵀ‖ is not. So a bracket is the right instrument where nothing sharper exists, and the wrong one wherever an error is directly available — which is the reason a comparison of methods and a decision about shift count take their numbers from different places. The bracket answers how many shifts to buy, where a factor of four in the error is a shift or two in k, and the geometric fall of max|r|² makes a shift or two cheap. It does not answer which method is better.

What follows for a run that cannot see its own error #

Plan against the prediction and stop against the residual. They are available at the two moments a run needs a number, and neither requires the Gramian. The shift count is chosen from max|r|² before the first solve; the run is stopped when the residual reaches a target, which is a quantity computed from the factor’s own width.

Set the residual target below the error target, by the measured factor and not by one. On this family the error is between 1.46 and 3.70 times the residual, so a residual target equal to the accuracy wanted delivers an answer up to 3.7 times worse than asked for. Dividing the target by four covers every measurement here. That is the concrete version of what the zero that means it is finished insists on generally: a stopping quantity has to be a statement about the thing being decided, and if it is only proportional to it the constant belongs in the rule.

Calibrate the upper half once if a comparable small problem exists, and not otherwise. A coarse-grid version of the same model gives the constant, and the constant is a property of the system rather than of n — 0.7025 to 0.7915 across a fourfold refinement here. A system of a different kind does not, and 0.436 is what the wrong calibration looks like.

Report the bracket rather than either end. An error known to lie in [residual, 0.79 · max|r|²] is a more honest object than a single estimate, and it costs no extra computation because both ends were computed anyway. This is the same discipline the residual the method reports argues for when a printed residual and the true one have parted company: say which quantity is being reported.

And keep the factor factored, since the residual requires it. Recompressing Z to a smaller rank mid-run changes the residual and the error together, so the stopping rule survives the truncation — which is the arrangement the rounding that was not the problem reaches in a hierarchical arithmetic, where the object is never assembled and the question becomes when to truncate rather than whether.

The refusal this page publishes is the reason the first of those rules is worth anything. Fed the claim that the ADI error is a thousandth of its own rational factor — the shape of the belief that a bound taken over an interval must be hopelessly loose — the assertion rejects it, on a twenty-state model, at the ratio it actually measures. The prediction is not a safety margin with orders of slack in it. It is the answer times 0.75, and a method whose accuracy is knowable to six per cent before the first solve is a rare thing to have in a field where the accuracy itself is unmeasurable.