Exact at the points that were named

The bound that is known in advance computes a reduced model with an error bound available before the model exists, and the bound turns out to be an equality. It is the best-behaved error statement on this site. It is also the one nobody uses at scale, and the reason is in the first paragraph of that essay: it costs two Lyapunov equations, each a Schur factorisation of an n × n matrix, which is O(n³).

A finite-element model of a structure has n in the hundreds of thousands. O(n³) is not expensive there; it is unavailable. Every large-scale reduction in use is the other kind, and the other kind gives up the bound entirely.

What it does instead #

Pick r points σ₁ … σᵣ in the complex plane. Build

V = span{(σᵢI − A)⁻¹B},     W = span{(σᵢI − Aᵀ)⁻¹Cᵀ},

and project: Aᵣ = WᵀAV, Bᵣ = WᵀB, Cᵣ = CV, with W rescaled so WᵀV = I. Then

Hᵣ(σᵢ) = H(σᵢ)   and   H′ᵣ(σᵢ) = H′(σᵢ)

exactly, at every one of the r points. Measured here: 4.7·10⁻¹⁶ for the values and 2.0·10⁻⁸ for the derivatives, the second being the accuracy of the finite difference used to compare them rather than of the interpolation.

Two things about that are worth separating.

It is a theorem, not an approximation. The reduced function passes through the original at the points and leaves with the same slope, in exact arithmetic, for any distinct σ. There is nothing to converge and no tolerance to set.

It costs 2r shifted solves. One per point for V and one for W. On a sparse model those are exactly the operations a sparse direct solver or a preconditioned iteration already does, so the whole method is available at a scale where nothing in the previous three essays is.

The half nobody states #

The derivative condition is the part usually left out of the sentence, and it is half the value.

r points buy 2r conditions, not r. A reduced model of order r has 2r degrees of freedom in its transfer function — r poles and r residues — so 2r conditions is exactly enough to determine it. Interpolating only the values would leave it underdetermined and something else would have to fix the rest.

The second solve per point is what buys it, and dropping it is not a saving. Build the same projection one-sided, with W = V, and the values still agree to 10⁻¹⁶ while the derivatives are out by 4.6·10⁻⁵ — three orders worse than the two-sided model’s, which is the whole of what the extra solve was for. Nothing in a statement about subspaces distinguishes the two constructions, and the measurement does.

That is the same shape as the plane survives what its vectors do not: a property that belongs to a subspace and a property that does not, told apart only by measuring the thing the subspace was supposed to guarantee.

And the exactness is exact at every order, while what happens between the points is not.

One column does not move and the other falls by four hundred thousand. At r = 2, 3, 5 and 6 the worst error at a named point is 1.01·10⁻¹⁵, 1.52·10⁻¹⁵, 8.88·10⁻¹⁶ and 1.52·10⁻¹⁵ — the rounding level, with no trend — while the worst error away from one runs 0.0607, 0.00915, 2.97·10⁻⁵ and 1.51·10⁻⁷.

So the promise is order-independent and the quality is not, which is the trade stated as two numbers rather than as a sentence. Adding a point buys nothing at all where the model was already exact — it could not — and buys about a decade and a half everywhere else. That is a good rate, and it is a rate rather than a guarantee: nothing in the construction says the next point will buy another decade and a half, and the four measured values do not fit a law tightly enough to promise one. The method’s own claim is the flat column and the useful behaviour is the falling one, and only the first of the two is something a caller can rely on before running it.

What has been given up, stated plainly #

Everything away from the points.

The error curve above falls to 10⁻¹⁶ at each σ and reaches 4·10⁻³ between and beyond them, and there is no bound on it anywhere. Not a loose bound — none. The method offers 2r exact conditions and makes no statement about any other frequency.

Set the two methods side by side:

• Balanced truncation asks the user for nothing, bounds ‖H − Hᵣ‖∞ over all frequencies before the model exists, and costs O(n³).
• Interpolation asks the user for r complex numbers with no information about where to put them, bounds nothing, and costs 2r sparse solves.

The field’s central trade is between those two sentences, and it is a trade about cost rather than about accuracy — a point the measurements in interpolating at the model’s own poles make sharper than any argument could.

Where the points have to go #

The user’s r numbers are the whole method, and there is no information available for choosing them. That is not a rhetorical flourish: at the moment the choice is made, nothing about H has been computed.

Three things are nonetheless true and worth knowing.

Points inside the range of the poles work and points beyond it do not, which is the same observation the rate the condition number predicts makes about where a Krylov method spends its effort. The model here has poles running from 9.9 to 580; the points 1.5, 5, 16, 50, 160, 500 span them. Two points both far past the last pole give two solves that are both nearly B/σ and nearly parallel, and the sixth direction does not exist — the construction refuses rather than returning a projection quietly built on rounding. That refusal is deliberate and is the sort of thing this field’s routines should do.

The error between two adjacent points falls when a point is added between them, and the far end does not improve. Past the last point there is nothing left to match, so the error there is set by the model’s behaviour at infinity rather than by how many points were bought.

Points that are close together are nearly free and nearly useless. Two σ close to each other give two nearly parallel solves, which is the beginning of the same failure the sixth point above hit outright. What that means practically is that the spread of the points is the decision and their exact placement is not — which is a relief, since nothing is known about the exact placement. It is the same relief choosing without knowing finds one field over: the parameter matters by orders and the precision with which it is chosen does not.

How far apart “not close together” is #

That relief is worth a number, because a per cent is the figure a reader would supply and it is wrong by two orders. Holding the two endpoints at 1.5 and 500 — the same span every time, so nothing about the spread changes — and sliding a pair of inner points along at a fixed separation:

separation between the inner pair    H∞ error across four placements    ratio
        0.03 decades  (7%)               9.5·10⁻⁶ to 2.9·10⁻³          311×
        0.08 decades  (20%)              9.5·10⁻⁶ to 1.6·10⁻³          169×
        0.20 decades  (a factor 1.6)     9.6·10⁻⁶ to 3.6·10⁻⁴           37×
        0.35 decades  (a factor 2.2)     9.9·10⁻⁶ to 5.6·10⁻⁵          5.6×
        0.60 decades  (a factor 4)       9.8·10⁻⁶ to 1.9·10⁻⁵          2.0×
        0.84 decades  (a factor 7)       1.09·10⁻⁵ to 1.09·10⁻⁵        1.00×

The claim holds above about a third of a decade and fails badly below it. At a separation of a factor of seven the placement is worth nothing at all — four different positions give the same error to three digits. At a separation of seven per cent it is worth a factor of three hundred, which is larger than anything else in this essay.

So the working rule is a factor of two between adjacent points, not a per cent. Below that the exact placement decides two and a half orders of the answer, which is exactly the situation the relief was supposed to rule out — and a user who reads the spread is the decision and then places two points thirty per cent apart has not followed the advice they think they have.

The first column carries a second reading that is easy to miss. The best placement at every separation gives the same error — between 9.5·10⁻⁶ and 1.1·10⁻⁵ in all six rows. Separating the points does not buy a better model; it buys protection against a bad placement. There is a good answer available at every separation and the tight ones make it hard to hit, which is why the symptom is a spread rather than a shift.

And there is a mechanism behind the number rather than a coincidence. The two solves at σ₁ and σ₂ differ by (σ₁I − A)⁻¹ − (σ₂I − A)⁻¹, which the resolvent identity writes as (σ₂ − σ₁)(σ₁I − A)⁻¹(σ₂I − A)⁻¹B — so the angle between the two directions is proportional to the separation, and the basis’s conditioning is inversely proportional to it. Below a factor of two the two directions are close enough that which of them the orthogonalisation keeps is decided by the rest of the configuration, and the model that comes out is a model of whichever direction survived. Above it the two are genuinely different directions and the projection is doing what the theorem describes.

That is the same quantity a basis that is the same subspace and not the same thing measures at one point with several moments, arriving here at several points with one moment each. Both are the resolvent’s directions running together, and both have the same repair: separate what is being asked for, or orthogonalise as the basis is built rather than after.

The subspace, and the basis that must not be the obvious one #

The interpolation conditions are conditions on the subspace V, not on the vectors spanning it. Any basis will do, in exact arithmetic.

In floating point the obvious basis — the moments (σI − A)⁻¹B, (σI − A)⁻²B, … at one point — is a Krylov sequence in disguise, and it loses independence at exactly the rate the power method converges. Its condition number reaches 7.7·10⁹ by the eighth vector while the same eight solves spent at eight different points stay at 1.0. A basis that is the same subspace and not the same thing is that measurement, and it is the reason every implementation orthogonalises as it builds rather than afterwards.

Why the reduced model is a projection and not a fit #

It is worth being clear about what kind of object has been produced, because the word interpolation invites the wrong picture.

Nothing has been fitted. There is no least-squares problem, no residual minimised, no tolerance. A subspace was built from r solves, the original model was projected onto it, and the interpolation conditions are a consequence of what that subspace contains — (σI − A)⁻¹B is in V, so the reduced resolvent reproduces the original one there.

The consequence worth carrying is that the reduced model inherits structure from the projection. Aᵣ = WᵀAV of a stable A need not be stable — a one-sided projection of a symmetric negative-definite A is, because then W = V and the projection is congruent, but a two-sided projection of a general A is not. The interpolation conditions hold either way. So a reduced model from this family can be a correct interpolant of a stable system and be unstable itself, which is a failure mode balanced truncation does not have and which no amount of adding points fixes.

What the interpolation conditions do not say #

It is worth being precise about the failure mode, because “no bound” is easy to read as “no guarantee of anything” and that is too strong.

The reduced model is a rational function of degree r that agrees with H to first order at r points. Between those points it is whatever the projection made it, and whatever the projection made it is not arbitrary — it is a rational function of degree r, so it cannot oscillate freely. What is missing is a statement about how far a degree-r rational interpolant can be from the function it interpolates, and no such statement exists without assumptions about H that nobody has.

The practical consequence shows up in the error curve. It is smooth, it has r zeros, and it rises between and beyond them in a way that looks entirely predictable once drawn. The problem is that drawing it requires H, which is the thing the reduced model exists to avoid evaluating. A user with the full model can check the reduction and does not need it; a user who needs it cannot check it.

That asymmetry is what makes the a-priori bound valuable out of proportion to its cost, and it is worth setting against the residual the method reports, which is the same problem in the iterative field: the quantity that can be computed and the quantity that matters are different quantities, and the gap between them is not bounded either.

The partial repair, and it is only partial, is to interpolate at more points than the model’s order and check the fit at the extras. That costs solves and turns the method into something closer to a fit, and it is what a careful implementation does. It still bounds nothing between the points it checked.

The cost, counted honestly #

2r shifted solves for the bases, and then the projection itself: r matrix–vector products with A per basis vector, so O(r²) products in total, and an r × r eigenproblem at the end that costs nothing.

Against balanced truncation’s two Lyapunov equations at O(n³) each, that is the entire argument. At n = 24, where every figure in this field is drawn, the comparison is meaningless — both are instant. At n = 10⁵ one of them exists and the other does not.

What makes the comparison worth drawing at n = 24 anyway is that the accuracy question can be settled at a size where both methods run. The answer, measured in interpolating at the model’s own poles, is that they land within a couple of per cent of each other in the norm the cheap one optimises — so the trade really is about cost, and the cheap method is not paying for its cheapness in the way one would expect.

The one thing that is bounded #

There is a guarantee in this method and it is not about accuracy. The interpolation conditions hold exactly, in floating point, at every point. Measured at 4.7·10⁻¹⁶ across five points, which is rounding — not “approximately”, not “to the tolerance of the solve”.

That matters more than it sounds. It means a reduced model built this way can be audited: hand it the same σ, evaluate both functions, and the agreement is a check on the whole construction — the solves, the orthogonalisation, the normalisation of W against V — in one number. A construction with a wrong sign in the biorthogonalisation still produces a plausible reduced model and fails that check immediately, which is how the implementation behind these figures was debugged.

So the method offers no bound on its answer and a perfect check on its own arithmetic, which is almost exactly the opposite arrangement to the rest of this site.

The structure the projection does and does not keep #

A reduced model is three small matrices, and what can be said about them beyond their size is worth a section, because it is where this method’s one genuine weakness sits.

Symmetry survives a one-sided projection. If A = Aᵀ and C = Bᵀ then W = V, the projection is congruent, and Aᵣ = VᵀAV is symmetric with the same sign of definiteness. That is a theorem and it carries stability with it: a symmetric negative-definite A projects to a symmetric negative-definite Aᵣ, so the reduced model is stable whatever points were chosen.

Nothing survives a two-sided projection. Aᵣ = WᵀAV with W ≠ V is a general square matrix, and there is no statement that its eigenvalues lie in the left half-plane. So a reduced model built by the method this essay is about can be unstable while interpolating a stable system exactly at r points, and no amount of adding points removes the possibility.

That is a failure with an unusual character. The model is not inaccurate — it matches H at r points to rounding, and its error curve between them is ordinary. It is unusable, because a simulation driven by it grows without bound while the system it represents does not. A caller has to check the poles of Aᵣ, and if they are in the wrong half-plane the only recourse within this method is to choose different points and try again.

Balanced truncation does not have this failure: truncating a balanced realisation of a stable system gives a stable system, by a theorem. That is a second entry on its side of the ledger, alongside the bound, and it is one the cost comparison usually leaves out.

So the honest statement of the trade is three-sided rather than two. Interpolation is cheap and exact where it is asked; balanced truncation is expensive, bounded everywhere, and cannot return something the caller must reject. Which of those matters most depends on whether the reduced model is going to be looked at or run.

The shape of the whole field, in three sentences #

A reduced model is an approximation of a function, so its error is a norm of a function and can be measured exactly rather than estimated. One method computes that error’s bound before the model exists and costs more than a large problem has; the other is exact at points the caller names, bounds nothing between them, and costs what a large problem can afford. Everything else in this field is a consequence of which of those two a problem’s size forces.