Why a Gramian can be truncated at all

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.

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ₖ([a,b], [−b,−a])²,

where Zₖ 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:

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ₖ², 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.

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:

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.

Both halves of that last sentence are checkable, and the second is the one with a number missing:

   model         k=2     k=4     k=6     k=8    k=10    k=12   spread
   heat n=20    0.744   0.754   0.736   0.736  0.735   0.735    1.03
   heat n=30    0.756   0.772   0.747   0.748  0.746   0.746    1.03
   heat n=40    0.766   0.778   0.756   0.753  0.752   0.752    1.03
   lowDeg r=6   0.798   0.799   0.784   0.784  0.783   0.783    1.02
   modal geo12  0.579   0.461   0.437   0.437  0.436   0.432    1.34

Constant in k to three per cent on four of the five, across a range of k over which the error itself falls by orders — so not of k holds. And the value runs from 0.43 to 0.80 across the five, so a property of the system holds too, and the 0.75 is this model’s number rather than a constant of the construction. The quantity spans nearly a factor of two, and a reader carrying 0.75 to another system would be out by up to that.

The exception is the row that belongs there. The modal system is the one whose spectrum does not fill the interval — it is the same system the table two sections down shows falling away from the Zolotarev bound — and it is the one whose ratio is not constant at small k: 0.579 at k = 2, settling to 0.43 by k = 6. Its poles sit away from where |r| peaks, so which direction of P the maximum picks out depends on where the shifts have put the peak, and that stops moving only once there are enough shifts for the equioscillation to be established.

So the identity is an identity in k on any system whose spectrum fills the interval, and becomes one after half a dozen shifts on a system whose does not — which is the same condition, arrived at from a third direction, that the two sections below are about.

What a per-system constant is worth #

A ratio that is constant in k and varies by a factor of two across systems is an awkward object, and it is worth saying which of its two properties is usable.

The constancy is the usable half. Measure the ratio once, at any k, and the ADI error at every other k is max|r|² times that number — so a single run of the iteration calibrates a prediction for all the others, and the prediction costs a Zolotarev factor rather than a Lyapunov solve. That is a genuinely practical thing: choosing how many ADI steps to take is choosing k against a target accuracy, and with the constant in hand the choice is a closed-form inversion rather than a sweep.

The variation across systems is what stops it being a formula. There is no way to know 0.43 from 0.78 without running the iteration once, because the number is how much of P points where |r| peaks and that is a property of B and of the eigenvectors of A — neither of which appears anywhere in the Zolotarev argument. The bound is built out of two numbers and the identity needs a third that the bound was constructed to avoid needing.

Which is the honest summary of what an interval-based theory can and cannot do, and it is the same summary the section on the modal system arrives at: the interval predicts the rate and the system supplies the constant. One is free and the other is one iteration.

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 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:

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.

The size of the system moves the looseness too, and in the direction that is easy to get backwards.

At ten states the bound is loose by a factor of 3·10⁸. At sixteen it is loose by 5·10⁶.

The bound tracks κ and the truth does not. Between ten states and sixteen the bound loosens by a factor of 2.6 while the actual ratio rises by a factor of 143, so the slack falls from 3·10⁸ to 5·10⁶. The bound is becoming less wrong as the problem gets harder, which is the opposite of how a bound usually behaves and the reason it is worth reading across the slider rather than at one n.

So the bound tightens as the problem grows, by four and a half orders of magnitude between ten states and forty. Both quantities rise with n — the achieved ratio from 8.7·10⁻¹⁴ to 1.11·10⁻⁷ and the bound from 2.58·10⁻⁵ to 1.05·10⁻³ — but the achieved ratio rises far faster, because κ rises with it (48 to 681) and the interval the maximum is taken over widens. That is the useful direction: the bound is least informative exactly where the truncation is easiest and most informative where it is about to become a real decision. A reader who calibrated their trust in it on a small system would be trusting it least where it works best.

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ₖ² 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.