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 drag says the band moves with the Péclet number and does not close: three placements of nineteen at Pe = 4, six at 16 and at 40. So a caller cannot avoid it by a rule about where to put shifts, because where the band is depends on the system.

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, does not appear in the interpolation error, and is not what an error norm would report either — the H∞ norm of the difference is computed over the imaginary axis, and a right-half-plane pole does not sit on it. A model with a pole at +24.1 can have a perfectly ordinary-looking frequency response.

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.

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.

At other settings #