Numerical Linear Algebra
where the algebra stops being the arithmetic
Reduction, and what a model is for

Why a Gramian can be truncated at all

Every method in this field rests on one fact nobody states the reason for — the eigenvalues of a Gramian fall off a cliff. The equation defining it has a rank-one right-hand side and no low-rank structure anywhere — and the answer's decay is a rational approximation problem with a closed-form rate.
38 min read 34 figures Where the algebra stops being the arithmeticMeasured, not assertedExact ground truth

Worth reading first: The bound that is known in advance · The problem the solver was actually given · A block nobody can call sparse.

Every essay in this field so far has rested on one fact and none of them has said where it comes from: the eigenvalues of a Gramian fall off a cliff. If they did not, there would be no Hankel singular value decay, no truncation worth making, and the bound that is known in advance would report a number of the same size as ‖H‖∞ at every order.

The fact is not obvious and it deserves better than an appeal to smoothness. A Gramian is the solution of

AP + PAᵀ + BBᵀ = 0,

an n × n matrix defined by an equation whose right-hand side has rank one. There is no low-rank structure in the statement — P is generically full rank, and it is — and nothing in the equation suggests that its eigenvalues should span fifteen orders. They do, and the reason is a theorem about rational approximation with a rate that can be written down.

Why a Gramian can be truncated at all: a rational approximation problemλₖ₊₁(P)/λ₁(P) for the controllability Gramian of a 30-state model, against Z_k², where Z_k is the smallest a degree-(k,k) rational function with poles on one side can be made on the other. The Gramian's spectrum falls from 1 to 2.48·10⁻⁸ in 10 steps, and the bound falls with it. Nothing about the model enters the bound except the two ends of the spectrum of −A, 9.86 and 3834, whose ratio is 388.8 — so the decay that makes model reduction possible is predictable from two numbers before the Gramian exists. That is the fact this whole field rests on and it is a fact about rational approximation.1234567891010⁻¹⁹10⁻¹⁶10⁻¹³10⁻¹⁰10⁻⁷10⁻⁴10⁻¹kλₖ₊₁ ÷ λ₁, and the boundZ_k²the Gramianpredicted from two numbersstates30κ of the spectrum389λ11 ÷ λ₁2.5·10⁻⁸the bound there5.2·10⁻⁴the cliff everything rests onand the reason for it
Fig. 1 The cliff everything in this field rests on, and the bound that explains it.

The statement

Let a and b be the ends of the spectrum of −A, which for a stable system are both positive. Then

λₖ₊₁(P) / λ₁(P)  ≤  Z_k([a,b], [−b,−a])²,

where Z_k is a Zolotarev number: the smallest that a rational function of degree (k, k) with its poles in one interval can be made on the other, relative to its size there. It is a quantity from approximation theory with nothing to do with linear algebra, and it has a closed-form rate governed by b/a and by nothing else about the problem.

Measured on a thirty-state model whose spectrum runs from 9.9 to 3855 — a condition number of 389 — the ratio and the bound at each k:

k 2 4 6 8 10
λₖ₊₁/λ₁ 5.4·10⁻² 2.2·10⁻³ 7.0·10⁻⁵ 1.6·10⁻⁶ 2.5·10⁻⁸
Z_k² 3.8·10⁻¹ 7.4·10⁻² 1.4·10⁻² 2.7·10⁻³ 5.2·10⁻⁴

Every ratio is under its bound, and the bound falls geometrically. Nothing about the model enters it except two numbers: the ends of the spectrum of −A. Not B, not C, not where the eigenvalues sit inside the interval, not the size of n.

So the decay that makes this whole field possible is predictable before the Gramian exists, from a quantity that costs an estimate of the extreme eigenvalues of A to obtain. A code that wants to know whether a model is worth reducing has an answer for two eigenvalue estimates rather than for two Lyapunov solves.

Why a rational function has anything to do with it

The connection is not a coincidence and there is a route to it that stays concrete.

The ADI iteration for the Lyapunov equation produces, after k steps, a factor Z of size n × k with P ≈ ZZᵀ. Its error is exactly

P − ZZᵀ  =  r(A) P r(A)ᵀ,     r(z) = ∏ⱼ (z − pⱼ)/(z + pⱼ),

where the pⱼ are the shifts the user chose. So the algorithm’s accuracy after k steps is max|r| over the spectrum, squared — and the best possible choice of shifts is exactly the rational function that makes that maximum smallest, which is the Zolotarev problem.

Then the bound on the eigenvalues follows: if k steps produce a rank-k approximation with relative error Z_k², then λₖ₊₁(P)/λ₁(P) cannot exceed it, because the (k+1)th eigenvalue is what a rank-k approximation leaves behind. The decay of a Gramian and the convergence of an algorithm for it are the same statement, and the second is the concrete route to the first. It is worth putting beside an equation whose unknown is a matrix, which is where this site first solved a Lyapunov equation and where the object’s size was the whole difficulty.

The ADI error is the rational function, not merely bounded by it‖P − ZZᵀ‖/‖P‖ after k steps of low-rank ADI on a 30-state model, against max|r|² for the same shifts, with P taken from its closed form rather than from another solve. The two curves are parallel: their ratio is 0.756 at k = 2 and 0.746 at k = 10, a spread of 1.0340 while the error itself falls from 0.289 to 3.87·10⁻⁴. A bound that is a constant factor from its quantity over four orders is an identity with a projection in front of it — the constant is how much of P lies where |r| is largest, which depends on the system and not on k.234567891010⁻⁴10⁻²1ADI steps krelative error, and max|r|²max|r|²the measured errora bound that is an identitysteps drawn9ratio, worst0.77spread of the ratio1error at k = 103.9·10⁻⁴its residual1.1·10⁻⁴four orders of errorand one constant
Fig. 2 The measured error and the squared rational factor. Their ratio is one number.

An identity, not a bound

The measurement above is worth more than “the error is at most max|r|²”, because it is a ratio and the ratio does not move:

k 2 4 6 8 10
error ÷ max|r|² 0.756 0.772 0.747 0.748 0.746

Constant to three digits while the error itself falls from 2.9·10⁻¹ to 3.9·10⁻⁴. A bound that stays a fixed fraction of its quantity over four orders is not a bound — it is an identity with a projection in front of it, and the 0.75 is how much of P lies in the direction where |r| attains its maximum. A property of the system, not of k.

That is the difference this site keeps drawing between a statement that happens to hold and a statement that explains. The bound that is never attained is the growth factor’s 2ⁿ⁻¹, true and useless. This one is attained to within a fixed factor at every k, which is what makes the shift choice a design problem rather than a search.

The shifts, which are the poles of the rational function

If the accuracy of k ADI steps is max|r| for a rational function whose poles are the shifts, then choosing shifts is choosing a rational approximation, and the answer is the one Zolotarev gave: poles spaced geometrically over the interval.

Geometric spacing, and the obvious alternative:

The rational function 6 ADI shifts amount to, chosen two ways|∏(λ − pⱼ)/(λ + pⱼ)| over the spectrum of −A for a 30-point discretisation, whose ends are 9.86 and 3834 and whose condition number is 388.8. The lower curve uses 6 geometrically spaced shifts and touches zero at each of them, equioscillating between with a maximum of 0.1187. The upper curve uses 6 equally spaced ones: every shift sits in the top decade of the spectrum, the small end is never approached, and the maximum is 0.8922 — a factor of 7.52 worse, which the iteration pays squared. The two runs are the same algorithm, the same number of solves, and one different line where the shifts are chosen.10¹10²10³10⁻⁹10⁻⁷10⁻⁵10⁻³10⁻¹λ, over the spectrum of −A|r(λ)|geometric maxequally spacedgeometricone line of codeshifts6κ of the spectrum389geometric max0.12equally spaced max0.89the factor between7.5the same k solvesand one choice of where
Fig. 3 The rational function six shifts amount to, chosen two ways.

The geometric choice touches zero at each of its six shifts and equioscillates between them, with a maximum of 5.2·10⁻². The evenly spaced choice puts every shift in the top decade of a spectrum spanning three, never approaches the small end, and reaches 8.5·10⁻¹ — sixteen times worse, which the iteration pays squared: the measured ADI errors after eight steps are 2.0·10⁻³ and 5.9·10⁻¹, a factor of 290.

Same algorithm, same number of solves, one different line where the shifts are chosen. It is the identical measurement the scaling that buys ten orders makes about a change of units in the polynomial field, and the identical one where to put the poles of a rational function makes about a nonlinear eigenvalue problem. Three fields, one answer, and the answer is always geometric clustering towards wherever the function is difficult.

What the bound is a bound over, and what that costs it

The Zolotarev factor is a maximum over an interval, so it is a statement about where the eigenvalues are allowed to be rather than where they are. That is what makes it computable in advance, and it is also its limitation.

Two systems spanning the same interval — the heat model, whose eigenvalues crowd towards the small end like sin², and a modal system with twelve poles spaced geometrically across it:

k 4 6 8 10
heat model 5.7·10⁻² 1.1·10⁻² 2.0·10⁻³ 3.9·10⁻⁴
modal system 1.5·10⁻² 1.3·10⁻³ 9.0·10⁻⁵ 3.2·10⁻⁶
the bound 7.4·10⁻² 1.4·10⁻² 2.7·10⁻³ 5.2·10⁻⁴

Both stay under it. The heat model tracks it at a constant 0.75 — its spectrum fills the interval, so there is always an eigenvalue where |r| is near its maximum. The modal system falls away from it, reaching 1/160 of the bound by k = 10, because its twelve poles happen to sit near the zeros of r and the maximum over the interval is attained where there is nothing.

That is the same structural fact the tropical prediction ran into in the polynomial field, — the same reading a rank that is a number of digits gives to an off-diagonal block, where what is being bounded is a worst case over a family and a particular member can be far better — and it is worth stating as a rule rather than as an instance: a quantity built out of a maximum is exact where the answer is a maximum and merely an upper bound where it is not. A system’s decay is predicted by its interval only to the extent that its spectrum fills the interval.

What follows for a large problem

The practical content is the reason the whole apparatus exists.

A Gramian of a hundred-thousand-state model cannot be formed — it has 10¹⁰ entries — and does not need to be. What is wanted is a low-rank factor Z with P ≈ ZZᵀ, and ADI produces one at k columns for k shifted solves. The bound says how many columns are needed for a given accuracy before any of them is computed: k such that Z_k² is small enough, which is a function of b/a alone.

So the sequence for a large problem is: estimate the extreme eigenvalues of A, compute the required k from the closed-form rate, choose k geometric shifts, run k solves, and hold an n × k factor instead of an n × n Gramian. Every step of that is priced in advance, which is unusual, and the whole of it is available at a scale where forming P is not. That is the same accounting the accuracy worth paying for does for a hierarchical solve: a cost stated in advance, against an accuracy stated in advance, with the trade made once rather than discovered.

It 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, hold a factor, truncate as the iteration proceeds — arriving here with a bound on the truncation that field does not have. It is also the same move the rounding that was not the problem makes about recompression: the object is kept in a factored form throughout and the question becomes when to truncate rather than whether.

What the condition number does to the rate

The bound depends on b/a and on nothing else, so the natural question is how sharply. It is worth answering with numbers because the answer decides which problems this field works on.

The rate is geometric in k with a factor that behaves like exp(−π²/log(16κ)) for large κ. That logarithm in the denominator is the important part: the rate degrades logarithmically in the condition number rather than linearly. Multiplying κ by a thousand does not multiply the number of required ADI steps by anything like a thousand; it multiplies it by roughly the ratio of the two logarithms, which is a small factor.

That is why the whole field survives discretisation refinement. The second-difference operator has κ growing like n², so a grid ten times finer has a hundred times the condition number — and needs about a third more ADI steps, not a hundred times as many. On the models here, κ runs from 32 at n = 10 to 674 at n = 40, and the number of steps to reach 10⁻⁶ goes from six to about nine.

It also says which problems are hard. A system whose spectrum spans twenty decades — a multiphysics model with a stiff chemical timescale beside a slow thermal one — has log(16κ) around fifty, and the required k is several times larger. That is a real cost and it is a predictable one, which is the whole point of having a rate rather than an observation.

There is a limit to how far this reasoning goes, and it is worth marking. The bound is about a real interval, and it applies as stated when the spectrum of A is real. A system with lightly damped oscillatory modes has complex poles close to the imaginary axis, the relevant approximation problem is over a region in the plane rather than an interval, and both the rate and the shift selection are harder. Everything in this essay is a statement about the real case, which is the case of a diffusion, and the oscillatory case is where the field’s remaining difficulty lives.

The two things the shifts have to be

There is a distinction buried in the shift discussion that is easy to slide past, and it separates two different requirements that happen to have the same answer here.

The shifts have to make max|r| small over the spectrum, which is the approximation requirement and is what the Zolotarev problem answers. And each shift has to give a solvable and well-conditioned shifted system, because the ADI step is a solve with A − pI: a shift landing on an eigenvalue makes that system singular.

Geometric spacing satisfies both at once, and it is worth noticing that this is not automatic. The shifts sit inside the interval the eigenvalues occupy, so a shift can in principle land near one — and on a spectrum with n eigenvalues in an interval and k ≪ n shifts, the chance of that is small but is not zero. A careful implementation checks the residual of each ADI solve rather than assuming it.

The sign of the shift is the other place this construction can go quietly wrong, and it is where the implementation behind these figures did. The shifts are stated here as positive numbers on the spectrum of −A, because that is the parameterisation the Zolotarev factor is written in; the ADI shift is their negative, so the matrix solved with is A − pI and not A + pI. With p taken positive the second is singular exactly where the shift is a good one, which produces an iteration that diverges fastest when its shifts are best — a failure mode with the right shape to be mistaken for a theoretical problem.

The search, and whether geometric is optimal

Geometric shifts are not exactly the Zolotarev optimum, which involves elliptic functions. The difference is worth measuring rather than citing, and the measurement has to be constructed carefully: a descent started at the answer that finds no improvement is indistinguishable from a descent that does not work.

So the search starts from the evenly spaced shifts, which are orders worse. It walks them from a factor of 8.9·10⁻¹ to 1.9·10⁻¹ — a real improvement, so the descent works — and lands beside the geometric value of 1.2·10⁻¹, within a factor of 1.6 of it.

Two things follow. The geometric choice is near-optimal and a code without an elliptic-function library loses very little by using it. And the objective has a wide, flat optimum: a set of shifts within a factor of two of the best is not hard to find, which is consistent with everything else in this field where the spread of a set of points is the decision and the exact placement is not.

At other settings

Why a Gramian can be truncated at all: a rational approximation problemλₖ₊₁(P)/λ₁(P) for the controllability Gramian of a 10-state model, against Z_k², where Z_k is the smallest a degree-(k,k) rational function with poles on one side can be made on the other. The Gramian's spectrum falls from 1 to 8.7·10⁻¹⁴ in 9 steps, and the bound falls with it. Nothing about the model enters the bound except the two ends of the spectrum of −A, 9.8 and 474.2, whose ratio is 48.37 — so the decay that makes model reduction possible is predictable from two numbers before the Gramian exists. That is the fact this whole field rests on and it is a fact about rational approximation.12345678910⁻¹⁹10⁻¹⁶10⁻¹³10⁻¹⁰10⁻⁷10⁻⁴10⁻¹kλₖ₊₁ ÷ λ₁, and the boundZ_k²the Gramianpredicted from two numbersstates10κ of the spectrum48λ10 ÷ λ₁8.7·10⁻¹⁴the bound there2.6·10⁻⁵the cliff everything rests onand the reason for it
Fig. 4 Ten states, where the spectrum is narrow and the decay is fastest.
Why a Gramian can be truncated at all: a rational approximation problemλₖ₊₁(P)/λ₁(P) for the controllability Gramian of a 16-state model, against Z_k², where Z_k is the smallest a degree-(k,k) rational function with poles on one side can be made on the other. The Gramian's spectrum falls from 1 to 1.24·10⁻¹¹ in 10 steps, and the bound falls with it. Nothing about the model enters the bound except the two ends of the spectrum of −A, 9.84 and 1146, whose ratio is 116.5 — so the decay that makes model reduction possible is predictable from two numbers before the Gramian exists. That is the fact this whole field rests on and it is a fact about rational approximation.1234567891010⁻¹⁹10⁻¹⁶10⁻¹³10⁻¹⁰10⁻⁷10⁻⁴10⁻¹kλₖ₊₁ ÷ λ₁, and the boundZ_k²the Gramianpredicted from two numbersstates16κ of the spectrum116λ11 ÷ λ₁1.2·10⁻¹¹the bound there6.7·10⁻⁵the cliff everything rests onand the reason for it
Fig. 5 Sixteen, and κ has grown with n².
Why a Gramian can be truncated at all: a rational approximation problemλₖ₊₁(P)/λ₁(P) for the controllability Gramian of a 22-state model, against Z_k², where Z_k is the smallest a degree-(k,k) rational function with poles on one side can be made on the other. The Gramian's spectrum falls from 1 to 1.19·10⁻⁹ in 10 steps, and the bound falls with it. Nothing about the model enters the bound except the two ends of the spectrum of −A, 9.85 and 2106, whose ratio is 213.7 — so the decay that makes model reduction possible is predictable from two numbers before the Gramian exists. That is the fact this whole field rests on and it is a fact about rational approximation.1234567891010⁻¹⁹10⁻¹⁶10⁻¹³10⁻¹⁰10⁻⁷10⁻⁴10⁻¹kλₖ₊₁ ÷ λ₁, and the boundZ_k²the Gramianpredicted from two numbersstates22κ of the spectrum214λ11 ÷ λ₁1.2·10⁻⁹the bound there2.1·10⁻⁴the cliff everything rests onand the reason for it
Fig. 6 Twenty-two, where the bound has loosened and the decay has slowed with it.
Why a Gramian can be truncated at all: a rational approximation problemλₖ₊₁(P)/λ₁(P) for the controllability Gramian of a 40-state model, against Z_k², where Z_k is the smallest a degree-(k,k) rational function with poles on one side can be made on the other. The Gramian's spectrum falls from 1 to 1.11·10⁻⁷ in 10 steps, and the bound falls with it. Nothing about the model enters the bound except the two ends of the spectrum of −A, 9.86 and 6714, whose ratio is 680.6 — so the decay that makes model reduction possible is predictable from two numbers before the Gramian exists. That is the fact this whole field rests on and it is a fact about rational approximation.1234567891010⁻¹⁹10⁻¹⁶10⁻¹³10⁻¹⁰10⁻⁷10⁻⁴10⁻¹kλₖ₊₁ ÷ λ₁, and the boundZ_k²the Gramianpredicted from two numbersstates40κ of the spectrum681λ11 ÷ λ₁1.1·10⁻⁷the bound there0.001the cliff everything rests onand the reason for it
Fig. 7 Forty, and both curves have moved together.
The ADI error is the rational function, not merely bounded by it‖P − ZZᵀ‖/‖P‖ after k steps of low-rank ADI on a 10-state model, against max|r|² for the same shifts, with P taken from its closed form rather than from another solve. The two curves are parallel: their ratio is 0.732 at k = 2 and 0.703 at k = 9, a spread of 1.0415 while the error itself falls from 0.119 to 1.81·10⁻⁵. A bound that is a constant factor from its quantity over four orders is an identity with a projection in front of it — the constant is how much of P lies where |r| is largest, which depends on the system and not on k.2345678910⁻⁵10⁻³10⁻¹ADI steps krelative error, and max|r|²max|r|²the measured errora bound that is an identitysteps drawn8ratio, worst0.73spread of the ratio1error at k = 91.8·10⁻⁵its residual6.9·10⁻⁶four orders of errorand one constant
Fig. 8 The identity on a small model, where the constant is different and still constant.
The ADI error is the rational function, not merely bounded by it‖P − ZZᵀ‖/‖P‖ after k steps of low-rank ADI on a 22-state model, against max|r|² for the same shifts, with P taken from its closed form rather than from another solve. The two curves are parallel: their ratio is 0.776 at k = 2 and 0.769 at k = 10, a spread of 1.0296 while the error itself falls from 0.248 to 1.61·10⁻⁴. A bound that is a constant factor from its quantity over four orders is an identity with a projection in front of it — the constant is how much of P lies where |r| is largest, which depends on the system and not on k.234567891010⁻⁴10⁻²1ADI steps krelative error, and max|r|²max|r|²the measured errora bound that is an identitysteps drawn9ratio, worst0.79spread of the ratio1error at k = 101.6·10⁻⁴its residual4.6·10⁻⁵four orders of errorand one constant
Fig. 9 And on a larger one, where the same k buys less.
The ADI error is the rational function, not merely bounded by it‖P − ZZᵀ‖/‖P‖ after k steps of low-rank ADI on a 40-state model, against max|r|² for the same shifts, with P taken from its closed form rather than from another solve. The two curves are parallel: their ratio is 0.766 at k = 2 and 0.752 at k = 10, a spread of 1.0346 while the error itself falls from 0.336 to 7.87·10⁻⁴. A bound that is a constant factor from its quantity over four orders is an identity with a projection in front of it — the constant is how much of P lies where |r| is largest, which depends on the system and not on k.234567891010⁻⁴10⁻²1ADI steps krelative error, and max|r|²max|r|²the measured errora bound that is an identitysteps drawn9ratio, worst0.78spread of the ratio1error at k = 107.9·10⁻⁴its residual2.1·10⁻⁴four orders of errorand one constant
Fig. 10 Forty states, and the ratio has still not moved with k.
The rational function 2 ADI shifts amount to, chosen two ways|∏(λ − pⱼ)/(λ + pⱼ)| over the spectrum of −A for a 30-point discretisation, whose ends are 9.86 and 3834 and whose condition number is 388.8. The lower curve uses 2 geometrically spaced shifts and touches zero at each of them, equioscillating between with a maximum of 0.6181. The upper curve uses 2 equally spaced ones: every shift sits in the top decade of the spectrum, the small end is never approached, and the maximum is 0.9731 — a factor of 1.57 worse, which the iteration pays squared. The two runs are the same algorithm, the same number of solves, and one different line where the shifts are chosen.10¹10²10³10⁻⁹10⁻⁷10⁻⁵10⁻³10⁻¹λ, over the spectrum of −A|r(λ)|geometric maxequally spacedgeometricone line of codeshifts2κ of the spectrum389geometric max0.62equally spaced max0.97the factor between1.6the same k solvesand one choice of where
Fig. 11 Two shifts, where the equioscillation has one interior maximum.
The rational function 4 ADI shifts amount to, chosen two ways|∏(λ − pⱼ)/(λ + pⱼ)| over the spectrum of −A for a 30-point discretisation, whose ends are 9.86 and 3834 and whose condition number is 388.8. The lower curve uses 4 geometrically spaced shifts and touches zero at each of them, equioscillating between with a maximum of 0.2711. The upper curve uses 4 equally spaced ones: every shift sits in the top decade of the spectrum, the small end is never approached, and the maximum is 0.9341 — a factor of 3.45 worse, which the iteration pays squared. The two runs are the same algorithm, the same number of solves, and one different line where the shifts are chosen.10¹10²10³10⁻⁹10⁻⁷10⁻⁵10⁻³10⁻¹λ, over the spectrum of −A|r(λ)|geometric maxequally spacedgeometricone line of codeshifts4κ of the spectrum389geometric max0.27equally spaced max0.93the factor between3.4the same k solvesand one choice of where
Fig. 12 Four, and the geometric curve has four zeros.
The rational function 8 ADI shifts amount to, chosen two ways|∏(λ − pⱼ)/(λ + pⱼ)| over the spectrum of −A for a 30-point discretisation, whose ends are 9.86 and 3834 and whose condition number is 388.8. The lower curve uses 8 geometrically spaced shifts and touches zero at each of them, equioscillating between with a maximum of 0.05198. The upper curve uses 8 equally spaced ones: every shift sits in the top decade of the spectrum, the small end is never approached, and the maximum is 0.8497 — a factor of 16.3 worse, which the iteration pays squared. The two runs are the same algorithm, the same number of solves, and one different line where the shifts are chosen.10¹10²10³10⁻⁹10⁻⁷10⁻⁵10⁻³10⁻¹λ, over the spectrum of −A|r(λ)|geometric maxequally spacedgeometricone line of codeshifts8κ of the spectrum389geometric max0.052equally spaced max0.85the factor between16the same k solvesand one choice of where
Fig. 13 Eight, where the evenly spaced curve has barely moved.
The rational function 12 ADI shifts amount to, chosen two ways|∏(λ − pⱼ)/(λ + pⱼ)| over the spectrum of −A for a 30-point discretisation, whose ends are 9.86 and 3834 and whose condition number is 388.8. The lower curve uses 12 geometrically spaced shifts and touches zero at each of them, equioscillating between with a maximum of 0.009967. The upper curve uses 12 equally spaced ones: every shift sits in the top decade of the spectrum, the small end is never approached, and the maximum is 0.766 — a factor of 76.9 worse, which the iteration pays squared. The two runs are the same algorithm, the same number of solves, and one different line where the shifts are chosen.10¹10²10³10⁻⁹10⁻⁷10⁻⁵10⁻³10⁻¹λ, over the spectrum of −A|r(λ)|geometric maxequally spacedgeometricone line of codeshifts12κ of the spectrum389geometric max0.01equally spaced max0.77the factor between77the same k solvesand one choice of where
Fig. 14 And twelve, where the two are five orders apart at the small end.
Hankel singular values, and the half of them a product cannot seeσₖ divided by σ₁, for a 20-state model of McMillan degree 12. The lower curve is the square-root route — Cholesky-like factors of the two Gramians and one SVD of RᵀS — and it descends to the level at which the Gramians were computed. The upper curve eigendecomposes the product PQ, agrees for the first 6 values and then stops, flattening at 9.67·10⁻¹⁰. The dashed line is σ₁√u = 2.32·10⁻⁹, which is where a route that squares the conditioning must stop: the eigenvalues of PQ are σ², an absolute error of uσ₁² on them is a relative error of √u on σ. κ(PQ) is 1.8·10³⁶ against 1.34·10¹⁸. Every σ below the line is a term in the error bound, so this is not an academic loss.1357911131510⁻¹⁸10⁻¹⁵10⁻¹²10⁻⁹10⁻⁶10⁻³1index kσₖ ÷ σ₁σ₁√utwo factorstheir productthe normal equations, againroutes agree to k =6product floors at9.7·10⁻¹⁰σ₁√u2.3·10⁻⁹κ(P)κ(Q)1.8·10³⁶√ of it1.3·10¹⁸do not form the productthe σ below the line are the bound
Fig. 15 The numbers this decay eventually produces.
The error a reduced model has, and the bound computed before it existed‖H − Hᵣ‖∞ for balanced truncation of a 20-state model, against 2Σ_{k>r}σₖ — a number available from two Lyapunov solves before any reduced model is formed — and against σᵣ₊₁, which is a lower bound no reduced model of that order can beat. The measured error and the upper bound are the same curve: the ratio is 1.0000 to 1.0000 across the sweep, spread 1.0000. The bound is attained rather than approached, because on a spectrum falling this fast the tail sum is its own first term and removing one state costs exactly twice it. The measured error sits at 2.03 times the lower bound, so the whole question of what order to truncate at is answered by a curve nobody has to compute a model to see.123456710⁻¹¹10⁻⁹10⁻⁷10⁻⁵10⁻³10⁻¹order r kepterror and bounds2Σσ, and the errorσᵣ₊₁a bound that is an equalityorders7bound ÷ error, worst1spread over the sweep1error ÷ σᵣ₊₁, worst2.1computed before the modeland attained by it
Fig. 16 And the bound they are a sum of.
Why the product route stops: two condition numbers and their productThe two Gramians of a 22-state model of McMillan degree 10, and what each route to the Hankel singular values is charged. κ(P) = 1.26·10¹⁸ and κ(Q) = 2.16·10¹⁸. The square-root route works with RᵀS, whose condition number is their geometric mean, 1.65·10¹⁸; the route that eigendecomposes PQ works at 2.72·10³⁶, which is past 1/u = 4.5·10¹⁵ — the point at which nothing small survives at all. The Gramians themselves are right: this one agrees with its closed form to 2.15·10⁻¹³. What is lost is lost in the last step, to a product nobody had to form.condition numbers, on a logarithmic scaleκ(P)1.26·10¹⁸κ(Q)2.16·10¹⁸√(κ(P)κ(Q)) — the SVD route1.65·10¹⁸κ(P)κ(Q) — the product route2.72·10³⁶1/u4.5·10¹⁵the mean, or the productand only one of them fits
Fig. 17 The conditioning that is the same fact seen from another angle.
The two 3×3 matrices of a Sylvester equation, and the 9×9 matrix it meansA and B are 3×3 with 5 nonzeros each. The coefficient matrix of the map X ↦ AX + XB is 9×9 with 21 nonzeros, drawn at the same scale so the ratio is visible rather than quoted. Bartels and Stewart's algorithm never forms it: two Schur reductions and a back-substitution give the same X to 5.5·10⁻¹⁷, and that X satisfies AX + XB = C to 8.8·10⁻¹⁷.A3×3B3×3I ⊗ A + Bᵀ ⊗ I9×9one equation, two objectsentries in A and B18entries in the coefficient matrix81two routes, relative gap5.5·10⁻¹⁷‖AX + XB − C‖/‖C‖8.8·10⁻¹⁷the small squares are the problemand the large one is the notation
Fig. 18 The equation whose solution decays.
sep(A, B) and the smallest |λᵢ(A) + μⱼ(B)| as the departure from normality grows, n = 6A and B are upper triangular with diagonals 1, 2, …, 6 and μ above them, so every eigenvalue sum is at least 2 at every point on the axis and the flat line is exact rather than nearly flat. sep — the smallest ‖AX + XB‖_F over ‖X‖_F = 1, and the reciprocal of the amplification a perturbation of C receives — starts equal to the gap at μ = 0, where the matrices are normal, and falls to 9.47·10⁻⁴ at μ = 8. The number a reader is invited to consult is the one that does not move.01234567810⁻⁵10⁻⁴10⁻³10⁻²10⁻¹110¹μ, the entry above the diagonalsep, and the eigenvalue gapmin |λᵢ + μⱼ|sep(A, B)the spectra never moveeigenvalue gap, throughout2sep at μ = 02sep at μ = 89.5·10⁻⁴amplification there1056solvability is the eigenvaluesand conditioning is not
Fig. 19 And the quantity that governs how well it can be solved.
Multiplications in a 2-dimensional model solve, through the eigenbasis against a dense factorisationA Kronecker sum's eigenvectors are the Kronecker products of its factors' eigenvectors, so the change of basis that diagonalises an operator with 1,024 rows is 2 changes of basis along 2 indices. The whole solve is a transform, 1,024 divisions and a transform back: 1.31·10⁵ multiplications at n = 32, against a dense factorisation's 7.16·10⁸, a factor of 5461. The fitted exponent is 3.00 against 3d = 6. The only decompositions taken are of the 2 one-dimensional factors, which on the model problem are the same matrix — so there is exactly one, of size 32. Every solve reproduces its right-hand side to 1.44·10⁻¹³.10¹10²10⁴10⁶10⁸n, points along one axismultiplicationsa dense factorisationthrough the eigenbasis2 decompositions of an n × nunknowns1024multiplications1.3·10⁵dense factorisation7.2·10⁸fitted exponent3‖Ax − b‖ ⁄ ‖b‖1.4·10⁻¹³nothing of size n^dis ever factorised
Fig. 20 The route to it nobody takes, which is n² × n².
What two lines of scaling are worth: the same quadratic in nine systems of unitsThe forward error against the closed form for an overdamped chain of 8 masses, before and after Fan–Lin–Van Dooren scaling. The unscaled curve runs 7.49·10⁻¹⁴ to 0.001258 — every digit gone by the far end — and the scaled one runs 1.26·10⁻¹³ to 8.23·10⁻¹⁴, flat to within a factor of 2.01. The scaling is γ = √(‖K‖/‖M‖) and δ = 2/(‖K‖ + γ‖C‖), computed from three norms and nothing else, and the map back is λ = γμ with no rounding in the statement. Flatness is the half that matters: after scaling every stop of the sweep IS the same problem, so there is nothing left for the change of units to do.0246810⁻¹⁶10⁻¹³10⁻¹⁰10⁻⁷10⁻⁴10⁻¹log₁₀ γ, the change of unitsforward erroras writtenafter scalingone change of variableunscaled, worst0.0013scaled, worst1.7·10⁻¹³orders recovered10scaled coefficient spread4.5the answer was never the problemthe units were
Fig. 21 A choice made before the computation, worth ten orders, in another field.
A prediction from three numbers, exact at one end of the spectrum and wrong by n² at the otherThe tropical roots of the max-plus quadratic with coefficients ‖M‖, ‖C‖ and ‖K‖, divided by the extreme moduli they are supposed to estimate, for chains of 4 to 32 masses with C = 6M + 0.5K. The large root is within 5.7 per cent of the largest modulus at every size. The small root is out by 5.36, 17.1, 60.8 and 229.6 — a factor growing like the square of the size. The reason is structural: a norm is a maximum, the small end of this spectrum is set by the SMALLEST eigenvalue of K, which is 4sin²(π/2(n+1)), and no norm of K contains that number. A quantity built out of maxima is exact where the answer is a maximum and silent where it is a minimum.11.31.610⁻¹110¹10²10³log₁₀ ntropical root ÷ actual modulusexactsmallest moduluslargest modulusthree norms, two endslarge root, n = 40.96large root, n = 320.94small root, n = 45.4small root, n = 32230a maximum predicts a maximumand says nothing about a minimum
Fig. 22 A prediction built out of maxima, exact at one end and not the other.
Growth factor under partial pivoting to n = 12: the bound, the worst case, and realityGrowth factor against matrix size on a logarithmic vertical axis. The two-to-the-n bound rises as a straight line; Wilkinson's matrix sits exactly on it; random matrices stay near one.02468101214110¹10²10³10⁴10⁵10⁶matrix size ngrowth factor max|u| / max|a|the 2ⁿ⁻¹ boundworst of 30 randommedian randomWilkinson's matrix sits on the bound30 Gaussian matrices per sizeat n = 12: bound 2048, worst 3.6
Fig. 23 A bound that is never attained, for contrast.
What a rank-10 approximation can achieve, by spectrumA semi-logarithmic plot of singular value against index for three spectra — geometric decay, algebraic decay, and flat — with the rank-ten approximation error marked on each.1112131415110⁻³10⁻²10⁻¹1index jσⱼσ11geometric, 0.85ʲalgebraic, j^−0.00flatmeasured, and equal to σ₁₁rank-10 error, geometric0.2rank-10 error, algebraic1rank-10 error, flat160×60, spectrum chosen rather than the entriesthe matrix decides, not the method
Fig. 24 A spectrum with no decay, where none of this applies.
Error of the best rank-k approximation to a 12×12 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.123456789101110⁻¹⁰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. 25 What a decay is worth once it exists.
Conjugate gradients on the three-dimensional model problem with every iterate cut to rank 1The falling curves are relative residuals: the lower one is the same method with no budget, which reaches 1.61·10⁻¹⁴ in 40 steps, and the upper one is the budgeted run, which stops at 0.119. The two step curves near the top are ranks, on their own scale: the un-truncated step asks for 5 at every step from the third onwards and the budget allows 1. The truncation is therefore not an occasional tidy-up — it is happening at every step, and the distance between the two residual curves is what it costs. The solution of this problem is itself a train of rank five, so a budget of five or more removes nothing and the two curves coincide.081624324010⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²conjugate gradient steprelative residualabove, on their own scale: rank asked 5, rank kept 1solid: the budgeted residual · light: no budgetbudget 1rank asked for5rank kept1residual, budgeted0.12residual, unbudgeted1.6·10⁻¹⁴numbers stored30the step asks for moreat every step
Fig. 26 Truncating inside an iteration rather than after it.
The residual a rank budget leaves, against the budget, on two right-hand sides of the same operatorBoth runs use the same method on the same matrix for 70 steps and differ only in b. With a constant right-hand side the floor falls 0.118, 0.0107, 3.61·10⁻⁴, 5.4·10⁻⁶, 1.07·10⁻¹³, 1.07·10⁻¹³, 1.07·10⁻¹³ — about two decades per column — and disappears at rank 5, because the solution *is* a train of that rank and the truncation has stopped removing anything. With independent normal entries the same ladder gives 0.969, 0.925, 0.874, 0.801, 0.735, 0.656, 0.47: 8 columns of a 10-column answer buy a factor of 2.06. So a floor is not a property of the truncation. It is the distance from the answer to the set the truncation projects onto, and the same budget is free on one problem and useless on the next.012345678910⁻¹³10⁻¹⁰10⁻⁷10⁻⁴10⁻¹rank budgetresidual the run stalls atb with no structure at allb constant: the answer is a traintwo ladders, one truncationstructured, rank 10.12structured, rank 45.4·10⁻⁶structured, rank 81.1·10⁻¹³unstructured, rank 10.97unstructured, rank 80.47the floor is not the truncation'sit is the answer's
Fig. 27 The floor a truncated iteration eventually reaches.
Numbers stored per unknown, against the size of the matrix, at 10⁻¹⁰The dense matrix stores n numbers per unknown, which is the straight line through the origin and doubles whenever n does. The hierarchical representations store 57, 87, 118, 151 and 56, 80, 108, 140, adding about 31 per doubling — a constant per doubling is a logarithm, and it is the whole claim. At n = 512 that is 30 and 27 per cent of the dense matrix, and the share falls at every size. The exponent between consecutive sizes is 1.60, 1.45, 1.35, falling towards one and never arriving, which is what a logarithm looks like from inside.56789100108216324432log₂ nnumbers stored per unknowndenseweak: every off-diagonal blockstrong: only the admissible onesa constant per doublingstrong, n = 5127.7·10⁴weak, n = 5127.2·10⁴dense, n = 5122.6·10⁵per doubling31‖A − A_H‖ ⁄ ‖A‖4.9·10⁻¹²the dense line doublesand the other two add a constant
Fig. 28 What holding a factor instead of a matrix buys.
The recovered signal at three λ, 0.10% noiseThree recovered signals against the truth. At the best λ the two bumps are recovered and the step is rounded off. Fifty times smaller and the answer is dominated by amplified noise; fifty times larger and the step has gone entirely. The step is the feature that costs the most to represent and it is the first thing any filter loses.0816243240485600.511.5positionvaluethe truthbest λtoo littletoo mucheither side of the bestbest λ, relative error0.11fifty times smaller0.73fifty times larger0.15the step is the first thing every filter losesand it does not come back
Fig. 29 A different rational filter, in the regularisation field.
Two filters on one sum, λ = 0.001The weight each term of the solution is given, against its index. Truncation is a step: one for the first 32 terms and zero after. Tikhonov is σ²/(σ² + λ²), which falls smoothly through the same place. The unregularised solution is the constant one, which is why it divides noise by a σ of 1.7·10⁻¹³.081624324048566400.250.50.751index kfilter factor fₖno regularisation: fₖ = 1truncationTikhonovthe same sum, three weightsTikhonov, relative error0.16truncation, relative error0.13no filter at all5.5·10⁸both filters are one expression with a different weightfₖ = 1 is the catastrophe
Fig. 30 And the factors it applies.
Two routes to e^A as a 8×8 bidiagonal matrix walks up to a defective oneA_δ is upper bidiagonal with diagonal λ, λ+δ, …, so its eigenvalues are exact and its eigenvectors have a closed form. As δ falls the eigendecomposition route's relative error rises to 2.92·10⁶⁵ while scaling and squaring's falls to 4.09·10⁻¹², which is the distance ‖A_δ − A₀‖ = 7·10⁻¹² to the matrix whose exponential is known in closed form. The two agree at the right of the axis, where the matrix really is far from the limit; the divergence is κ(V), which reaches 3.3·10⁸².10⁻¹²10⁻¹⁰10⁻⁸10⁻⁶10⁻⁴10⁻²10⁻¹⁴10⁻¹10¹²10²⁵10³⁸10⁵¹10⁶⁴δ, the gap between consecutive eigenvaluesrelative error in e^Aan answer with no correct digitsV f(Λ) V⁻¹scaling and squaring‖A_δ − A₀‖exact eigenvalues throughoutκ(V) at the smallest δ3.3·10⁸²eigen route2.9·10⁶⁵scaling and squaring4.1·10⁻¹²distance to the limit7·10⁻¹²the eigenvalues are the diagonaland they are exact at every stop
Fig. 31 Rational approximation of a matrix function, one field over.
Relative error of each singular value of a bidiagonal matrix graded over 11 decades, by four routesThe matrix is 8×8, bidiagonal, and every entry is a power of two times a small integer, so the matrix the float routines are handed and the matrix the exact routine reasons about are the same matrix. The reference is a Sturm bisection in BigInt rationals, which never rounds. σ_min is 1.934·10⁻¹¹. One-sided Jacobi holds every value to 7.1·10⁻¹⁶ relative and the zero-shift bidiagonal sweep to 4.41·10⁻¹⁶; the eigenvalues of BᵀB are wrong by 14.5 — not a rounding error but the whole number — and the damage is entirely at the small end.1234567810⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹singular value, largest firstrelative errorone-sided Jacobizero-shift QRshifted QReigenvalues of BᵀBagainst a rational bisectionσ_min, exactly1.9·10⁻¹¹worst, one-sided Jacobi7.1·10⁻¹⁶worst, zero-shift QR4.4·10⁻¹⁶worst, eigenvalues of BᵀB15a relative error is a ratioand the denominator is the answer
Fig. 32 A spectrum spanning ten decades, where the rate is slow.
The transfer function of a 24-state model, by a solve and by a formula|H(iω)| for the second-difference operator on 24 points with one actuator and one sensor. The curve is C(iωI − A)⁻¹B, one factorisation of a complex 24 × 24 matrix per frequency, and the crosses are Σ(cᵀvₖ)(vₖᵀb)/(iω − λₖ) with the eigenvalues and eigenvectors written down in closed form. The two share the coefficient matrices and no code. Below ω = 100 they agree to 1.05·10⁻¹⁴ relatively; above it the modal sum is computing a number far below its own largest term, because the residues sum to cᵀb and the sensor is not at the actuator, so the relative comparison stops meaning what it looks like. Measured against the size of the terms being added, the two agree to 3.85·10⁻¹⁵ at every frequency. A state dimension of 24 and a function of one variable: the gap between those is what the field is about.10⁻¹110¹10²10³10⁻¹¹10⁻¹⁰10⁻⁹10⁻⁸10⁻⁷10⁻⁶10⁻⁵10⁻⁴10⁻³10⁻²10⁻¹frequency ω|H(iω)|one function, two routesstates24relative, ω ≤ 10010⁻¹⁴against the terms3.9·10⁻¹⁵‖H‖∞0.0924 statesone function of one variable
Fig. 33 The model all of this is about.
Why a Gramian can be truncated at all: a rational approximation problemλₖ₊₁(P)/λ₁(P) for the controllability Gramian of a 30-state model, against Z_k², where Z_k is the smallest a degree-(k,k) rational function with poles on one side can be made on the other. The Gramian's spectrum falls from 1 to 2.48·10⁻⁸ in 10 steps, and the bound falls with it. Nothing about the model enters the bound except the two ends of the spectrum of −A, 9.86 and 3834, whose ratio is 388.8 — so the decay that makes model reduction possible is predictable from two numbers before the Gramian exists. That is the fact this whole field rests on and it is a fact about rational approximation.1234567891010⁻¹⁹10⁻¹⁶10⁻¹³10⁻¹⁰10⁻⁷10⁻⁴10⁻¹kλₖ₊₁ ÷ λ₁, and the boundZ_k²the Gramianpredicted from two numbersstates30κ of the spectrum389λ11 ÷ λ₁2.5·10⁻⁸the bound there5.2·10⁻⁴the cliff everything rests onand the reason for it
Fig. 34 The hero again, beside the shifts that achieve it.

Where this goes next

123456710⁻¹¹10⁻⁹10⁻⁷10⁻⁵10⁻³10⁻¹order r kepterror and bounds2Σσ, and the errorσᵣ₊₁a bound that is an equalityorders7bound ÷ error, worst1spread over the sweep1error ÷ σᵣ₊₁, worst2.1computed before the modeland attained by it The bound that is known in advance Almost every error on this site is measured after the fact. Balanced truncation has one that is computable before the reduced model exists, in a norm of a function rather than of a residual — and on ordinary problems it is not an upper bound that is loose. It is attained. 081624324010⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²conjugate gradient steprelative residualabove, on their own scale: rank asked 5, rank kept 3solid: the budgeted residual · light: no budgetbudget 3rank asked for5rank kept3residual, budgeted3.8·10⁻⁴residual, unbudgeted1.6·10⁻¹⁴numbers stored150the step asks for moreat every step An iterate that must be made smaller Applying a Kronecker-sum operator to a low-rank iterate multiplies its ranks by d and adding two of them adds their ranks, so a solver in a compressed format cannot keep what it produces. Every step is followed by a truncation — and whether that truncation is a floor on the residual depends on the right-hand side rather than on the truncation. 02040608010010⁻¹¹10⁻¹⁰10⁻⁹10⁻⁸truncations performed by the factorisationrelative errorthe representation's own error‖A − LLᵀ‖ ⁄ ‖A‖no decomposition without its residualresidual, 0 truncations2.1·10⁻¹⁰residual, 98 truncations1.1·10⁻⁹representation, deepest1.4·10⁻⁹residual ⁄ representation0.81levels, deepest5a hundred approximate stepsand an exact-looking factorisation The rounding that was not the problem A rank-k block plus a rank-k block is a rank-2k block, exactly, so every arithmetic in this format truncates after every addition. A Cholesky performed inside it does ninety-eight of those and its residual is 1.14·10⁻⁹ against a representation error of 1.40·10⁻⁹ — the roundings cost nothing measurable.

What links here

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

Reads more easily once this is understood

Essays that name this one as worth reading first.

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.

ADI iterationCondition numberEquioscillationGramianHankel singular valuesLow-rank approximationLyapunov equationRational approximationSingular valuesZolotarev number