Half the conditions and a certificate

A model that cannot be run reduced a stable convection–diffusion system of thirty states to order three by matching its transfer function at three points, and got back a model with a pole at +24.1 at one placement and unstable models at four placements of nineteen. Every one of them matched the full system at its three points to the rounding level. It compared that construction with balanced truncation, which is stable by theorem, and concluded that the cheap method needs a check the expensive one does not.

Two constructions were on the table. There is a third, and it sits between them in cost and has a stability certificate of its own.

A two-sided reduction builds two bases, V from solves with the shifted state matrix and W from solves with its transpose, and projects the model as WᵀAV. The second basis is what buys the derivative conditions: exact at the points that were named measured the two-sided model matching the full system’s slope to 2·10⁻⁸ where a one-sided model was out by 4.6·10⁻⁵. A one-sided reduction skips the second basis and projects with V on both sides — a Galerkin projection, VᵀAV for an orthonormal V. It matches values and not slopes.

What it gets in exchange was never measured in these essays, because every model the reduction field had used before was symmetric, and on a symmetric model the two constructions coincide.

A pole that cannot leave a region #

The two-sided curve is the one the earlier essay drew: stable at small and large placements, a spike to +24.1 between them, four placements in the right half plane. The two one-sided curves never come near zero. The order-3 model’s rightmost pole wanders between −4.2 and −10.7 across the placements, the order-6 model’s never right of −3.2, and neither comes near the dashed line at −0.986.

That line is not a fit. It is the edge of the numerical range of the state matrix — the set of values xᵀAx over unit vectors x — and its right edge is the largest eigenvalue of the symmetric part (A + $A^{\mathsf{T}}$)/2.

The certificate takes two lines. A one-sided model’s state matrix is VᵀAV with V orthonormal. If λ is one of its eigenvalues with unit eigenvector y, then the real part of λ is the real part of y*VᵀAVy, which is (Vy)*A(Vy) with Vy also a unit vector, and so it lies inside the numerical range of A. Every pole of every Galerkin model lies left of the numerical range’s edge. If that edge is negative, every Galerkin model is stable, at every placement, at every order, whatever the shifts.

A two-sided model has no such argument. WᵀAV is not a Rayleigh quotient of A; it is a projection along one subspace onto another, and its eigenvalues can be anywhere — which is what the spike is.

Why the edge is where it is #

On this operator the symmetric part can be written down. The centred-difference convection–diffusion matrix is tridiagonal with diagonal −2ε/h² and off-diagonals ε/h² ± 1/(2h). The convection term is the antisymmetric part: it adds 1/(2h) above the diagonal and subtracts it below. So (A + $A^{\mathsf{T}}$)/2 is exactly the diffusion matrix, and its largest eigenvalue has the closed form

which at Pe = 10 and thirty states is −0.98612, the value the figure draws, to twelve digits.

That makes the certificate universal on this family and gives it a margin that can be read off without computing anything: one-sided models of a centred convection–diffusion operator are stable whatever the Péclet number, by a margin of about $\pi^2$/Pe.

Across eight Péclet numbers and nineteen placements each — 152 models of order 3 and 150 of order 6, two placements at Pe = 4 and 6 refusing to add six independent directions from shifts that close together — not one one-sided pole lies right of the edge, and not one lies in the right half plane. Two-sided models at the same placements are unstable at three to ten of nineteen, rising to eight and ten at Pe = 100 and 60.

Where the certificate and the spectrum agree #

The weakness of the certificate has a converse, and it is worth seeing before the weakness itself: at low Péclet number the certificate is nearly the whole story.

At Pe = 4 the worst one-sided pole of order 3 is −3.455 and of order 6 −3.44, against a system abscissa of −3.461: the Galerkin models keep the system’s own decay rate to the second digit at their worst placement. At Pe = 6 the three numbers are −3.13, −3.14 and −3.14. On a nearly symmetric operator the numerical range hugs the spectrum, and a Galerkin model’s poles, confined to the range, have nowhere else to go.

That is the regime the reduction field’s earlier essays lived in without naming it. On a symmetric system the numerical range is the interval between the extreme eigenvalues, the one-sided and two-sided constructions are the same construction, and the certificate is exact. What convection does is open a gap between the range and the spectrum — and the two-sided construction’s instability and the one-sided construction’s slow decay are two views of that one gap. The two-sided model, unconstrained, can put a pole anywhere the badly conditioned projection sends it; the one-sided model, constrained to the range, can put one anywhere in a range that now reaches much closer to zero than the spectrum does.

Even at Pe = 4 the two-sided construction fails at three placements, and fails badly. So the gap does not have to be large for the unconstrained construction to escape; it only has to exist. A range that is exactly the spectrum — a symmetric matrix — is the only case in which the two-sided model is certified too.

The margin is not the system’s #

The certificate’s weakness is visible in the same figure, and it runs in the opposite direction from the Péclet number.

As convection grows, the spectrum of the centred operator moves left: its rightmost eigenvalue is −3.46 at Pe = 4, −3.49 at 10, −11.5 at 40 and −19 at 100. The numerical range’s edge moves the other way: −2.47, −0.986, −0.247, −0.099. At Pe = 100 the system decays at a rate of nineteen and the certificate promises only a tenth. The more non-normal the operator, the wider the gap between what the system does and what the certificate can promise about a model of it.

That is the numerical range doing what a limit the matrix never reaches and the direction the error leans found this operator’s non-normality doing everywhere else: the eigenvalues describe the long-time decay, and the symmetric part describes what can happen transiently. A certificate built from xᵀAx is a transient bound, and a reduced model can inherit the transient behaviour without inheriting the decay.

The measured poles sit well inside the certificate. The worst one-sided pole of order 3 is −4.21 at Pe = 10, −2.39 at 40, −2.15 at 100; of order 6, between −3.05 and −3.44 throughout. From Pe = 10 upward none comes within a factor of two of the edge. But none keeps the system’s own rate either: at Pe = 40 the full system decays at 11.5 and the worst reduced model at 2.4. A Galerkin model of a strongly convective system is stable and decays too slowly, and the certificate says why it cannot be guaranteed to do better.

At Pe = 40 the contrast is sharper than at 10. Six two-sided placements fail, most of the lower half of the range, and the one-sided curves run flat and stable through them — sitting between the edge and the spectrum, and closer to the edge’s side of it than the system is.

What the slope conditions were worth #

The certificate is not free. A one-sided model matches values at its interpolation points and not derivatives, so at the same order it matches half as many conditions, and it should be less accurate. It is.

At order 3 the one-sided model’s error runs from 1.4·10⁻² to 4.2·10⁻², and at nine of the ten placements where the two-sided model is stable the two-sided error is smaller, by a factor of up to five. That is the derivative conditions’ worth, and it is not small.

But the comparison at equal order is the wrong one on cost. A two-sided model of order 3 needs three solves with the shifted A and three with its transpose: six factorisations of a thirty-by-thirty tridiagonal matrix, or on a large problem six sparse factorisations. A one-sided model of order 6 needs six solves with the shifted A — the same count. Spent that way the six solves buy six value conditions instead of three values and three slopes, and a model twice the size.

At the ten placements where the two-sided model is stable, the one-sided model of order 6 is more accurate at seven: 2.9·10⁻³ against 1.2·10⁻² at the smallest placement, 5.2·10⁻³ against 6.6·10⁻³ at s = 1, 1.3·10⁻² against 1.8·10⁻² at s = 32. At the three placements where the two-sided model is unstable the one-sided one is simply the only usable model. Against balanced truncation of order 3 at 6.8·10⁻³, the one-sided order-6 model is better at eight placements of thirteen.

The pattern survives at Pe = 40, and it narrows. Every construction’s error rises — balanced truncation of order 3 by a factor of 3.7, of order 6 by forty; the one-sided order-6 model beats the stable two-sided order-3 model at six placements of eight, by margins from nearly nothing at s = 5.6 to a factor of two at s = 32; and it is never better than balanced truncation of order 6, which here reaches 4.8·10⁻³ against its best of 2.1·10⁻².

So the trade is not stability against accuracy. For a fixed number of solves the one-sided construction gives a larger model that is usually more accurate and always stable. The cost is in the model’s size — six states rather than three — which is exactly the quantity reduction exists to make small.

A certificate that can be checked before anything is built #

A stability check on a two-sided model has to happen after the model exists: compute its poles, and if one is in the right half plane, move the shifts and build again. The numerical range turns that around.

The edge of the range is the largest eigenvalue of a symmetric matrix, (A + $A^{\mathsf{T}}$)/2. On a large sparse system that is a Lanczos iteration using only products with A and $A^{\mathsf{T}}$, converging quickly at an extreme eigenvalue, and it answers one question about the system rather than one question per model. If the edge is negative, every one-sided model that will ever be built from that state matrix is stable — at every shift, at every order, from any iteration that chooses the shifts. The check is paid once, before the first solve.

Compare the question it replaces. Whether the full system is stable at all needs the rightmost eigenvalue of a non-symmetric matrix, which is the harder eigenvalue problem, and on this operator the harder one in a specific way: its eigenvalues are sensitive, and an iteration aimed at the rightmost of them works against the same non-normality that makes the reduction delicate. A model that cannot be run noted that on a large system the reduced model’s poles are cheap and the full system’s are not, so a reduced pole in the right half plane is ambiguous evidence about a system of unknown stability. A negative numerical-range edge is not ambiguous: it proves the full system stable too, since the spectrum lies inside the range, and it proves every Galerkin model stable in the same breath.

The converse is where the certificate is silent. A positive edge says nothing about the system — this operator in other coordinates has one while its eigenvalues stay put — and nothing about any particular model. A code that computes a positive edge learns only that it will have to check its models one at a time, as it would have for a two-sided construction.

Where each construction belongs #

Three constructions now have measured positions on this operator, and they sort by what the reduced model is for.

Balanced truncation at two Lyapunov solves and a pair of cubic factorisations gives stability, an error bound known in advance, and the best accuracy at every order. Why a Gramian can be truncated at all is the theorem behind the first and the bound that is known in advance the measurement of the second. On a system small enough to afford the Lyapunov solves it is the answer.

Two-sided interpolation gives the most accuracy per state from a handful of solves, and two approximants and one matrix size found it within a fraction of a per cent of balanced truncation on symmetric models. On this operator it needs its poles checked, and a failed check means moving the shifts.

One-sided interpolation gives a certificate for the price of the derivative conditions. It is the construction to use when a model will be simulated without anyone checking it, when the shifts will be chosen by an iteration that cannot be supervised, or when the model will be embedded in a larger computation whose stability depends on every component’s — and it asks for the order to be doubled to recover the accuracy.

Interpolating at the model’s own poles sits beside all three: the iteration that places two-sided shifts at mirrored poles converges to a stable model, and its intermediate models need not be. A one-sided version of that iteration would have no unstable intermediates at all, at the cost of the optimality condition that makes the two-sided one worth running.

What this rests on #

One system: thirty states, a centred difference, random input and output vectors, Péclet numbers from 4 to 100. Nineteen placements for the poles and thirteen for the errors, each a geometric set of shifts spanning a factor of four. The H∞ norms are computed on a frequency grid with a golden-section refinement, not exactly. The certificate is a theorem and holds for any matrix; its closed form holds only for the centred difference, where the convection term is exactly antisymmetric. An upwind discretisation puts convection on the diagonal and changes the symmetric part, and is not measured here.

The claim that has to fail #

The tempting reading of a stability guarantee is that the reduced model decays like the system. At Pe = 10 the system’s rightmost eigenvalue is −3.49 and the numerical range’s edge — the only margin the Galerkin argument gives — is −0.986. The refusal is fed the claim that the edge is no further right than the spectrum, and fails. At Pe = 40 the gap is a factor of forty-seven.

Still open: a certificate that belongs to the coordinates #

The inner product. The certificate used the Euclidean inner product: V orthonormal, xᵀAx. Nothing about the system chose that inner product; it came with the coordinates the state happens to be written in. The transfer function does not change if the state is rescaled, but the symmetric part of the rescaled matrix does, and so does the numerical range. Whether a stable system can be written in coordinates where its Galerkin models are no longer certified — and whether they then actually fail — is measured in a certificate written in coordinates.

Spending the solves differently. The comparison here gave the one-sided construction six shifts spanning the same factor of four as the two-sided three. Six shifts spread across the whole frequency range, or three shifts each used twice in a one-sided block with the derivative direction added by hand, are other ways to spend the same solves, and whether any of them recovers the two-sided accuracy at order 3 while keeping the certificate is unmeasured.