Interpolating at the model’s own poles

Exact at the points that were named ends on an uncomfortable fact: the method’s whole content is r complex numbers, and nothing is known about where to put them at the moment the choice is made.

There is one choice that is not a guess. If the interpolation points are the mirror images of the reduced model’s own poles — σᵢ = −λᵢ(Aᵣ) — then Hᵣ satisfies the first-order conditions for a best H₂ approximation of H. That is a condition on the answer rather than an input to it, so it defines a fixed point, and IRKA is the iteration that looks for one: interpolate at some points, read the poles of what comes back, mirror them, repeat.

What is checked at the end, and what is not #

The iteration converges here in thirteen steps from a logarithmically spread guess, and settles at 9.861, 40.154, 102.133, 302.827 — which is close to, and deliberately not exactly, the model’s own poles at 9.9, 39.3, 87.8, 154.6.

What is checked is not that it stopped moving. A fixed-point iteration that stops moving without satisfying its own condition has converged to nothing, and there is no convergence proof for this one to lean on. So the check is the condition itself: take the model the final points produced, compute its poles, mirror them, and compare against the points that produced it. Agreement to 2.8·10⁻¹² on this model.

That distinction is worth making loudly because the two checks look the same in an implementation and are not the same statement. move < tol is a claim about the iteration. σ = −λ(Aᵣ) is a claim about the answer, and only the second is what the theorem is about. It is the same discipline a tolerance that reads its own residual applies to an inexact Newton method: stop on the condition that matters, not on the quantity that is easy to watch.

What the condition is worth, measured #

Here is where the essay stops being about IRKA and becomes about the field.

At n = 24 both methods can be run. Balanced truncation’s two Lyapunov equations are instant, IRKA’s four solves are instant, and the two reduced models can be compared in both norms on the same system. Four systems, of different sizes, degrees and pole spreads:

system                        H₂ (BT ÷ IRKA)   H∞ (IRKA ÷ BT)
24 states, degree 8, r = 4        1.0079            1.156
24 states, degree 8, r = 3        1.0126            1.233
30 states, degree 10, r = 5       1.0042            1.123
24 states, degree 12, r = 4       1.0194            1.248

Both orderings hold on every system, so this is not a coin flip. IRKA wins H₂ every time, by 0.4 to 1.9 per cent. Balanced truncation wins H∞ every time, by 12 to 25 per cent.

Each method wins the norm it was designed for, which is the expected half. The unexpected half is the size: a method costing four sparse solves lands within half a per cent of a method costing two O(n³) factorisations, in the norm the cheap one optimises. And the expensive one’s advantage in its own norm is a fifth, not an order.

Reading that correctly #

The comparison is usually reported one way round — either as balanced truncation has a bound, so use it or as IRKA is H₂-optimal, so use it — and both readings are supported by measuring one norm and not the other. Measuring both is what makes the answer honest, and the answer is that neither method is meaningfully better than the other at accuracy.

So the field’s choice is about cost, and about what is known in advance:

• If n is small enough that O(n³) is available, balanced truncation gives an error bound before the model exists, which nothing else does.
• If it is not, IRKA gives a model within a per cent of that in H₂ and a fifth of it in H∞, for r solves and no bound at all.

That is a clean answer and it is not the one the cost difference would predict. It is worth putting beside four knobs and one floor, which found four different regularisation rules landing within a few per cent of each other for the same underlying reason: when several methods are all extracting the same information from the same problem, the differences between them are small and the differences in what they cost are not.

Why the H∞ ordering goes the other way #

It is worth one paragraph, because “each wins its own norm” is a slogan until there is a mechanism.

H₂ is an average over frequency — it is √(1/2π ∫|H − Hᵣ|² dω) — and IRKA’s conditions make the error orthogonal to the reduced model’s own tangent space, which is exactly what minimising an average means. H∞ is a maximum, and nothing in the interpolation conditions constrains a maximum — the distinction small compared to what draws between a norm and a componentwise statement, one level up: the error curve has r zeros and rises freely between them, and where it rises highest is not something the method looks at.

Balanced truncation, by contrast, drops directions ordered by a quantity that is itself a bound on the H∞ error, so its worst frequency is controlled by construction and its average is not. That is the same asymmetry the best approximation there is records for the truncated SVD, which is optimal in two norms at once and unusual for being so. The two methods are each blind to the other’s norm, and the measurement says how much that blindness costs: about a fifth, in both directions, which is smaller than either method’s advocates would guess.

What the mirrored poles mean #

The condition σ = −λ(Aᵣ) looks arbitrary written down and it is not, so it is worth unpacking once.

H₂ is a norm on the error function, and minimising it over all rational functions of degree r is a smooth optimisation over 2r real parameters. Setting the derivative to zero produces conditions on where the error function is orthogonal to the tangent space of the manifold of degree-r rational functions — and the tangent space at Hᵣ is spanned by objects whose poles are Hᵣ’s own. Working the inner products out gives Hermite interpolation of H by Hᵣ at the mirror images of Hᵣ’s poles.

So the condition is not a heuristic that happens to work. It is what stationarity looks like when it is translated into the language of interpolation. That translation is the whole reason a first-order condition can be met by an iteration that only ever builds interpolants: the optimisation problem and the interpolation problem turn out to be the same problem written in two ways.

What that does not give is a guarantee of reaching it. Stationarity is a local condition, the manifold is not convex, and a fixed point of this iteration is a local optimum at best. The measurements above compare against balanced truncation rather than against the global optimum, because the global optimum is not available to compute.

The iteration’s own failure modes #

Three, and all three are visible rather than silent.

It need not converge, which is the shape the same zero, and nothing was found gives to a different iteration in the Krylov field: a method whose stopping test can be met without the thing it tests for having happened. There is no proof, and on problems with clustered poles it can cycle. The routine reports whether it reached a fixed point rather than assuming it; a model returned from a run that did not converge is an interpolant at wherever it stopped, which is a perfectly valid reduced model with no optimality claim attached — the same distinction an answer that changes with the seed draws between a computed object and the guarantee it was supposed to carry.

Its poles can be complex. A complex interpolation point needs its conjugate present in the basis to keep the reduced model real, which is bookkeeping this implementation states rather than performs — it keeps the real part and reports how many complex poles it saw. On the models drawn here that count is zero, because the state matrix is symmetric and its reduced projections are too.

The reduced model can be unstable. A two-sided projection of a stable A need not be stable, and IRKA mirrors the poles of whatever it produced — so an unstable reduced model produces interpolation points in the wrong half-plane and the iteration walks away. This is the one failure this method has that balanced truncation cannot: balanced truncation of a stable system is stable, by a theorem.

What it costs per step #

r solves for V, r for W, and an r × r eigenproblem, repeated until the points stop moving. Thirteen steps at r = 4 is 104 sparse solves.

Against balanced truncation’s two O(n³) factorisations that is still the cheap option at any large n, but it is not the four solves the previous essay priced. The optimality is bought with a factor of the iteration count, and the iteration count is not known in advance and has no bound.

That changes the practical reading slightly. A single interpolation at sensibly spread points costs r solves and gets within a few per cent of optimal on these models; IRKA costs 13r and gets to optimal. Whether the last few per cent is worth thirteen times the solves is exactly the kind of question a parameter chosen on a smaller problem answers by measuring rather than by arguing, and the answer here is usually no — the same shape what a rebuild is worth prices for a factorisation that could be refreshed and mostly should not be.

The starting guess, and whether it matters #

An iteration with no convergence proof invites the question of how much the starting points decide, and the answer here is: the number of steps, and not the destination.

Started from points spread logarithmically over the frequency range the iteration takes thirteen steps at r = 4. Started from a tighter cluster it takes more and arrives at the same fixed point to the digits printed. That is the behaviour a local optimum with a wide basin produces, and it is consistent with what is generally reported for this iteration.

It is not a guarantee, and the honest statement is the narrow one: on the models drawn here, every start tried reached the same fixed point, and the routine reports the step count so that a run which took unusually long is visible rather than hidden inside a returned model. A method whose only guarantee is a fixed-point condition owes its caller the information that it reached one, and the number of times it went round.

The one thing the start does decide is whether a complex pole appears mid-iteration. A badly placed cluster can produce a reduced model with a conjugate pair, at which point the implementation’s restriction to real shifts bites and it keeps the real part. On these symmetric models that never happens; on a general one it is the first thing to check when the iteration behaves oddly.

What the comparison needed in order to be fair #

The four rows of numbers above are the point of the essay, and they took more care to produce than they take to read. Three things had to be held fixed for them to mean anything.

The same order. A reduced model of order four and one of order five are not comparable, and the temptation is real because IRKA’s natural output is whatever order its fixed point wanted. Both methods were asked for exactly r here.

The same measurement. Both errors are computed by forming the difference system H − Hᵣ and measuring it, rather than by sampling |H − Hᵣ| on a grid. A sampled maximum is a lower bound and would flatter whichever method happens to put its worst frequency between two samples. The H₂ norm comes from a Lyapunov solve on the difference system, which is a different routine from either reduction and therefore cannot favour one.

The same models, more than one of them. A single system would leave the ordering a coincidence. Four systems of different sizes, degrees and pole spreads, with both orderings holding on all four, is what turns “IRKA was better here” into “IRKA is better in this norm”.

The third of those is the one this site keeps having to insist on. One matrix is an anecdote, and the habit of checking a claim about typical behaviour across several problems is the same one randomisation does not create structure applies to seeds. Here there is no randomness — every model is deterministic — and the argument transfers anyway, because what is being ruled out is not noise but a coincidence of construction.

What could still be wrong with the comparison is worth stating too. All four models come from the same operator, so they share a spectral shape; a system with clustered poles or a non-symmetric state matrix might order the two methods differently, and nothing here rules that out. The claim these numbers support is narrow: on models of this kind, at these orders, each method wins its own norm by a margin small enough that cost decides. That is the claim, and it is the one the field’s practice already assumes without having measured.

What to take from four rows of numbers #

Not the ranking. The ranking is the expected result and it is worth about a fifth in one norm and half a per cent in the other.

What to take is that a method with a proved bound and a cubic cost, and a method with no bound at all and a linear one, produce reduced models of essentially the same quality on the problems where both can be run. That is not an argument for ignoring bounds — the bound is what makes the expensive method usable without checking, and checking is what a large problem cannot do. It is an argument for being clear about what a bound buys, which here is confidence rather than accuracy.

At other settings #