A model that cannot be run

Exact at the points that were named prices the two ways of building a reduced model and finds them within a fraction of a per cent of each other in the norm each was designed for. It names one advantage balanced truncation has that the cost comparison leaves out — that it cannot return an unstable model, while interpolation can — and it does not show one, because the symmetric families this field is built on never produce one.

They cannot. For a symmetric state matrix with the input and output vectors related by a transpose, the two-sided projection is a Galerkin projection onto a Krylov space, WᵀAV inherits the definiteness of A, and the reduced model is stable by construction. So the field’s own models are the wrong instrument for the question, and the sentence stayed a sentence.

Here is one.

Four of nineteen placements return a model with a pole in the right half plane, reaching +24.1 from a system whose own rightmost eigenvalue is −3.49. And the models in that band interpolate: each matches the full system’s transfer function at each of its own three points to 1.3·10⁻¹⁴.

What has and has not gone wrong #

Nothing in the construction failed. The bases were built, the projection was formed, the biorthogonality of the two bases holds to the rounding level — it is asserted inside the routine — and the interpolation conditions the method exists to satisfy are satisfied.

The reduced model is a rational function of degree three that agrees with the full system’s transfer function at three prescribed points. That is exactly what was asked for and exactly what was returned. It is also a model with a pole at +24.1, which means the differential equation it stands for has a solution growing like e²⁴ᵗ, and the system it was built from has nothing of the kind.

So the failure is not in the arithmetic and not in the algorithm. It is that the specification was incomplete, and the field’s other method satisfies the missing clause without being asked.

Why the asymmetry is necessary #

The state matrix here is a centred-difference convection–diffusion operator: tridiagonal, with off-diagonals ε/h² + 1/2h and ε/h² − 1/2h. The convection term is what makes the two unequal, and the relative asymmetry ‖A − Aᵀ‖/‖A‖ runs from 0.18 at Péclet 10 to 0.69 at 40.

That matters because the whole of the stability guarantee for a symmetric problem comes from the projection preserving a quadratic form. For symmetric A with W = V, xᵀ(VᵀAV)x = (Vx)ᵀA(Vx) < 0 for every x, so the reduced matrix is negative definite whenever A is. With W ≠ V nothing of the kind holds: WᵀAV is a projection of A onto one space along another, and a projection along a different space can move an eigenvalue anywhere.

This site already knows the family that makes such things happen. It is the one the direction the error leans and the convection field’s essays are built on: a non-normal operator, whose eigenvalues are real and negative and whose behaviour is not described by them. A projection of a non-normal operator is where the gap between the spectrum and the pseudospectrum turns into an outcome, which is the spectrum stops describing the matrix’s subject and is the honest explanation of this figure.

The band, and what decides it #

The instability is not a rule and it is not an accident: it is a band in the shift placement, and both edges of it are visible on the figure.

At small s the three interpolation points sit far below the interesting part of the frequency response; the projection is close to a moment matching at the origin, the bases are well conditioned, and the model is stable. At large s they sit far above it, the bases converge towards the leading directions, and the model is stable again. In between — where the points are near the frequencies at which the non-normality shows — the two bases are nearly orthogonal to each other, WᵀV is ill conditioned before it is normalised, and the projected matrix is a badly conditioned object whose eigenvalues have moved.

The sweep says the band moves with the Péclet number and does not close.

The count does not grow with the Péclet number and the severity does. At Pe = 4, 8, 16, 25 and 40 the placements returning an unstable model number 3, 4, 6, 4 and 6 of nineteen — between a sixth and a third, with no trend — while the worst right-half-plane pole runs 21.8, 81.6, 82.5, 614 and 106. The band’s width is roughly fixed and its depth is not, and neither is monotone.

And the interpolation is exact at every one of the five. The residual at the three matched points runs 2.94·10⁻¹⁵, 1.04·10⁻¹⁵, 6.72·10⁻¹⁵, 3.59·10⁻¹³ and 2.57·10⁻¹³ — the construction keeping its promise to twelve digits or better while returning, at Pe = 25, a model with a pole 176 times further into the right half plane than the original’s rightmost pole is into the left. A caller cannot avoid the band by a rule about where to put shifts, because its position depends on the system; and cannot detect it from the interpolation error, because there is not one.

What balanced truncation does instead #

The dashed line on the figure is balanced truncation of the same system at the same order. Its rightmost pole is −0.84, and at orders two to six the poles are −1.50, −0.84, −1.49, −2.33 and −3.27. Stable at every order, and not by luck.

Balanced truncation of a stable system is stable — it is a theorem, and its proof is the one why a Gramian can be truncated at all is built around: the balanced realisation has a Lyapunov solution which is diagonal and positive, and the leading block of a positive definite Lyapunov solution is a positive definite Lyapunov solution for the leading block of the system. The truncation inherits a certificate rather than being tested for one.

That is the shape of the difference, and it is worth stating in the site’s usual terms. Interpolation gives what was asked for, exactly, at r linear solves. Balanced truncation gives what was asked for approximately, at two Lyapunov equations and a pair of cubic factorisations, together with a guarantee nobody asked for and an error bound computable before the model exists.

Two approximants and one matrix size prices the first half of that trade and finds the cheap method within half a per cent of the expensive one. This figure is the other half of the price, and it is not a percentage: it is a failure mode that either happens or does not.

What the interpolation error looks like when the model is unstable #

Worth measuring rather than assuming, because a reader’s next thought is that the model must be wrong somewhere even if it is right at the three points.

It is. The interpolation conditions are three equations and the model has nine degrees of freedom, so matching at three points constrains a small part of it — everything between and beyond the points is whatever the projection produced. On a stable placement that remainder is a reasonable approximation of the response; on an unstable one it contains a pole in the right half plane, which means the model’s transfer function has a singularity where the true one is analytic.

So the honest description is not right at three points and wrong elsewhere, which is true of every interpolant. It is: the interpolant has a feature the original does not have, and that feature is not small and does not appear in the interpolation error.

Whether an error norm would have caught it #

The tempting explanation is that it could not: the H∞ norm of the difference is taken over the imaginary axis, a right-half-plane pole does not sit on it, so the error should be blind to the instability. That is the right conclusion from the wrong premise, and the measurement is more interesting than either.

At Péclet 4, 10 and 16 the norm is not blind at all: every unstable placement has a larger error than every stable one, so the ranking is exactly right. And the margin between them falls from 43% to four parts in a thousand.

At Péclet 40 it inverts. Every unstable model is more accurate in H∞ than the worst stable one — 6.1·10⁻² against 1.3·10⁻¹ — so a reader ordering the placements by error would put the models that cannot be simulated at the top of the list.

That is a worse thing to have than a blind statistic. A quantity that says nothing invites no conclusion; a quantity that ranks the failure correctly at three settings and backwards at the fourth invites the belief that a large error would have warned somebody, and there is no threshold on it that would have.

Two details worth keeping. The model with the pole at +24.1 has an H∞ error of 1.13·10⁻¹, five times the worst stable placement at that Péclet number — so perfectly ordinary overstates it in the other direction too. And the largest error on the Pe = 10 sweep, 2.72·10⁻¹, belongs to the placement whose abscissa is +2.13 rather than to the one at +24.1: the size of the right-half-plane pole does not order the error, which is what one expects if the badly conditioned projection is the cause and the pole is a symptom of it.

That is why this failure is a stability question rather than an accuracy one, and why the check is a check on the poles rather than a check on an error. It is also why the check is cheap and the alternative is not: the poles of an r × r matrix are one eigenvalue decomposition and they answer the question, while a frequency sweep is a sweep and it answers by 1.004 at one Péclet number and by the wrong sign at another.

The inversion also says something about what the H∞ error is measuring, and it is worth separating from the practical point. A reduced model’s error on the imaginary axis is dominated by wherever the two transfer functions disagree most, and on this family that is between the interpolation points rather than near them. A right-half-plane pole at +p contributes 1/(iω − p) to the model’s response, whose magnitude never exceeds 1/p — so a large instability is a small contribution to the norm, and the pole at +24.1 perturbs the response less than the pole at +2.13 does. The quantity is behaving exactly as it should; it is simply not a quantity about stability.

Which is the general form of the caveat this essay ends on. An error norm answers the question it was defined to answer, and a method that satisfies an incomplete specification produces a model the norm cannot object to. The missing clause here is one word long, and the instrument that checks it is not the instrument the field already had.

How much of the field this changes #

Not much, and saying so precisely is part of the finding.

The cost comparison stands. Two approximants and one matrix size prices interpolation at r solves against balanced truncation’s two Lyapunov equations and finds the cheap method within half a per cent in H₂. That measurement is unaffected: it was taken on symmetric systems, where the stability question does not arise, and the numbers are what they are.

The error bound stands. The 2Σσ bound the bound that is known in advance measures at a ratio of 1.0000 belongs to balanced truncation and is untouched by anything here.

What changes is the decision the comparison feeds. A reader who takes the cost comparison as a recommendation is choosing a method whose failure mode this essay measures, and the comparison did not price it. The repair is one clause: interpolation is the cheaper method and requires a check the other does not, and the check is a line of code — which is why the honest conclusion is a caveat rather than a reversal.

And one thing gets worse rather than better. On a large system the reduced model’s poles are cheap to compute and the full system’s are not, so the check confirms that the model is stable without confirming that it should have been. A reduced model with a right-half-plane pole from a system whose stability is unknown is ambiguous evidence, and the field has no cheap instrument for resolving it.

What a practitioner does about it #

Four things, and the ordering is the useful part.

Check. The eigenvalues of a 3 × 3 or 10 × 10 reduced state matrix cost nothing. A method that can return an unstable model and is used without checking is a method used carelessly, and the check is one line.

Move the shifts and try again. The band is a band. Four placements of nineteen fail here, so a second attempt at a different scale usually succeeds — which is what makes this a nuisance rather than a catastrophe, and also what makes it easy to never notice.

Or use IRKA, and check anyway. The iteration that places its shifts at the mirror images of its own poles converges, when it converges, to a model satisfying a necessary condition for optimality — interpolating at the model’s own poles is that essay. A fixed point of it has poles in the left half plane by construction of the mirroring, so the converged model is stable; the intermediate ones need not be, and a run stopped early can return one of them.

Or pay for the certificate. If the reduced model is going to be simulated, or put inside a controller, or handed to somebody who will not check, the method whose stability is a theorem is worth two Lyapunov solves. The reason this figure exists is that the cost comparison in this field was quoting a price without that clause in it.

Why this is the field’s own version of a familiar shape #

The machine field arrives at nearly the same sentence from the opposite direction, which is worth recording.

There, a computation satisfies every guarantee anybody stated — every run backward stable, every answer within its bound — and two of them disagree, because the specification had no clause about agreement between runs. Here, a construction satisfies every condition anybody stated — the interpolation holds, the biorthogonality holds — and the result is unusable, because the specification had no clause about stability.

In both cases the missing clause is a property nobody thought to ask for, and in both cases the more expensive alternative supplies it without being asked: order-independent summation there, balanced truncation here. And in both cases the price of the expensive option had been quoted with the extra property left out of the comparison, which is exactly the arithmetic this site keeps finding.

The general form, since this collection has now met it twice: when a cheap method and an expensive one are compared on the quantity the cheap one optimises, the cheap one wins, and the comparison is not the decision.

What the site’s own gate would have caught #

A closing note about the shape of the omission, since it is one this collection makes a habit of recording.

Every generator on this site asserts something about its own arguments before it draws. The reduction field’s generators assert interpolation errors, bound ratios, Gramian agreement with a closed form — every quantity the essays argue about. None of them asserted that a reduced model was stable, because on the symmetric families used there, stability was never in question.

That is exactly the pattern the fleet’s own rules warn about: an assertion that has never rejected anything proves nothing, and an assertion nobody wrote because the case never arose is worse. The routine that builds a projection now returns its poles’ rightmost real part, and the figure asserts that the sweep contains both stable and unstable models — which would fail if either the construction stopped producing the failure or started producing nothing else.

The system, and why it is the right instrument #

A note on the construction, because a demonstration that a method fails is only worth having if the case is not contrived.

The state matrix is a centred-difference discretisation of −εu″ + u′ on the unit interval: the convection–diffusion operator this site’s convection essays are built on, at thirty states. The actuator and the sensor are random vectors from a seeded stream, so the transfer function has no structure a projection could accidentally respect.

Three things make it the right instrument rather than a construction aimed at the answer.

It is stable, and that is checked rather than assumed. The rightmost eigenvalue is −3.49 at Péclet 10 and negative at every Péclet number in the sweep, checked whenever the figure is drawn. A demonstration that a projection destroys stability is worth nothing if the system did not have any.

It is a real discretisation. Nothing about it was chosen to fail; it is the operator a convection-dominated flow gives, at a Péclet number the centred difference can carry, and the same operator the site has drawn a dozen times for other reasons.

And the failure is not everywhere. Fifteen of nineteen placements return a stable model, which is the reason this is a hazard rather than a defect — and the reason the assertion published with it has to be a refusal of the claim that every placement fails.

One line #

A reduced model that matches a stable system’s transfer function exactly at every point it was asked about can have a pole in the right half plane, and the method whose stability is a theorem costs two Lyapunov solves more.