A certificate written in coordinates

Half the conditions and a certificate found a reduced model of a convection–diffusion system that could not be unstable. A one-sided projection — the state matrix restricted to a subspace and tested against the same subspace — has every pole inside the numerical range of the full state matrix, the set of values xᵀAx over unit vectors. On the centred-difference operator the right edge of that range is the largest eigenvalue of the diffusion term, about −$\pi^2$/Pe, and so every one-sided model at every placement and every Péclet number measured was stable.

The argument used a unit vector, which means it used a length, which means it used an inner product: the Euclidean one, in the coordinates the state happened to be written in. The system does not know about those coordinates. Its state is a vector of nodal values, and nothing prevents a code from storing them scaled — by a mesh weight, by a unit conversion, by a preconditioner applied to the state, or by a change of variables chosen for some other reason.

So the question is whether the certificate belongs to the system or to the way it was written down.

One system, a family of coordinates #

The rescaling is the simplest one that matters. Write the state as x = Dz with D diagonal, so that the k-th state is multiplied by γ raised to t times its index. The state matrix becomes $D^-$¹AD, the input $D^-$¹B and the output CD. The eigenvalues of $D^-$¹AD are the eigenvalues of A, and C D (sI − $D^-$¹AD)⁻¹ $D^-$¹B is C(sI − A)⁻¹B term for term: the transfer function is identical. Nothing a user of the model could measure has changed.

The tridiagonal matrix’s two off-diagonals change, though. The one below the diagonal is divided by $\gamma^t$ and the one above is multiplied by it. With γ the square root of the ratio of the two off-diagonals, the choice t = 1 makes them equal — the rescaled matrix is symmetric — and t = 2 swaps them. That square root is real because both off-diagonals are positive, which holds while the cell Péclet number is below one — the threshold the stencil that is not symmetric measured at 1.0000000000000002 — and it is 0.16 and 0.65 at the two Péclet numbers drawn here. Values of t below zero or above two make the matrix more lopsided than it was.

The symmetric part follows the off-diagonals. At t = 1 it is the whole matrix, so the numerical range shrinks to the interval between the extreme eigenvalues and its edge is the spectral abscissa itself, −3.49. At t = 0, the original coordinates, it is −0.986. At t = −0.5 it has already crossed zero, at +2.16, and at t = −3 it is +37.9. The same happens on the other side from t = 2.5.

What a positive edge lets through #

The certificate is sufficient, not necessary: a positive edge removes the guarantee and does not by itself make any model unstable. The bottom of the first figure is the measurement of what it lets through.

In the original coordinates and wherever the edge is negative, no one-sided model at any order is unstable, as the certificate requires. At t = −1, with the edge at +6.6 and the rescaling’s condition number at 112, one placement of nineteen returns an unstable model of order 1. At t = −1.5, condition number 1,200: eleven of order 1 and four of order 2. At t = −2: eleven and five. At t = −3, condition number 1.4·10⁶: eleven, six, one of order 3 and two of order 6.

The two-sided model does not notice. Its four unstable placements and its worst pole of +24.1 are the same in these coordinates as in the original ones, to the rounding level, and the figure checks that at every placement. That is not luck either: a two-sided projection builds V from solves with the rescaled matrix and W from solves with its transpose, the first basis transforms by $D^{-1}$ and the second by D, and the projected model is the same model in a different basis. The two-sided construction is invariant under a change of coordinates and has no certificate; the one-sided construction has a certificate and is not invariant. Each has what the other lacks, and for the same reason — the second basis that exact at the points that were named found buying the slope conditions is also what makes the projection independent of how the state is scaled.

Less rescaling is needed as convection grows #

At Pe = 40 the off-diagonals differ by a factor of 4.7 rather than 1.4, so each unit of t is a much larger rescaling: the symmetrising scaling alone has condition number 4.6·10⁹. The edge is −0.247 in the original coordinates and already +6.58 at t = −0.25, where the rescaling’s condition number is only 260. At that point ten placements of nineteen give an unstable model of order 1 and three give an unstable model of order 6. At t = −1, eight of order 3 are unstable.

Read against the condition number of the rescaling, the first unstable model does not order neatly — it appears at 112 at Pe = 10, 492 at Pe = 25 and 260 at Pe = 40 — but the point at which ten or more placements fail does: a condition number of 1,200 at Pe = 10, 492 at 25 and 260 at 40. The more non-normal the operator, the less rescaling it takes to make the symmetric part indefinite. The order-3 models need more rescaling — none unstable at Pe = 10 until 1.4·10⁶, up to nine at Pe = 40 — because a larger subspace averages over more of the directions a skewed inner product distorts.

A diagonal scaling with condition number in the hundreds is not exotic. A state vector that mixes quantities in different units has one built in, a nondimensionalisation chosen per variable has one, and a code that stores a state already multiplied by a diagonal preconditioner has one. None of these changes the system and every one changes the certificate.

A rescaled projection is a weighted projection #

The two families of results — models that fail in skewed coordinates, and models that are rescued by a weight — are one result seen twice, and the algebra that joins them is short.

A one-sided projection in the rescaled coordinates builds an orthonormal basis Z for the rescaled Krylov space and forms ZᵀD⁻¹ADZ. The rescaled Krylov space is $D^{-1}$ times the original one, so Z = $D^-$¹V for some basis V of the original space, orthonormal not in the Euclidean inner product but in the weighted one: ZᵀZ = VᵀD⁻²V = I. The reduced state matrix is then VᵀD⁻²AV. A one-sided projection in rescaled coordinates is a one-sided projection in the original coordinates, in the inner product ⟨x, y⟩ = xᵀD⁻²y. The same subspace, a different test.

Every rescaling is therefore a choice of weight, and the question of which coordinates are safe is the question of which weights P = $D^{-2}$ make PA dissipative — make xᵀPAx negative for every x. The original coordinates are P = I, which works by a margin of $\pi^2$/Pe. The symmetrising coordinates are one particular diagonal P, which works with the spectrum’s margin. The lopsided coordinates are diagonal Ps that do not work at all. And the Lyapunov solution is a dense P that works in any coordinates, because it is computed from the matrix rather than inherited from the way the state was stored.

That also explains why the two-sided construction does not care. A two-sided projection already chooses its own test space — W, built from solves with the transpose — and a rescaling changes V and W in ways that cancel. It has no inner product of its own to lose, and so it has no certificate to lose either.

Where the certificate is exact #

The same rescaling can be run the other way, and at t = 1 it is the best coordinate system the certificate could ask for.

In the symmetrising coordinates the Galerkin certificate’s margin is the system’s own decay rate. The worst one-sided pole of order 3 moves from −4.21 to −3.91 and of order 6 from −3.22 to −3.49: every Galerkin model decays at least as fast as the system, which half the conditions and a certificate found impossible to guarantee in the original coordinates. And the two-sided curve is still the same four unstable placements, because nothing about it changed.

So the one-sided construction’s weakness at high Péclet number — stable, but certified only to a margin of $\pi^2$/Pe — was a property of the coordinates rather than of the construction. On a tridiagonal operator whose off-diagonals have the same sign the symmetrising scaling is known in closed form, and projecting in it recovers the spectrum’s margin. The price is that the scaling itself is badly conditioned when convection is strong, 4.6·10⁹ at Pe = 40, so the recovered margin is bought with basis vectors whose entries span nine orders of magnitude.

An inner product that does not depend on the coordinates #

There is a choice of inner product that restores the certificate in any coordinates, and it is one the reduction field already solves for. If P is symmetric positive definite and AᵀP + PA is negative definite, then xᵀPAx is negative for every nonzero x, and a projection that tests against PV rather than V — the reduced state matrix (VᵀPV)⁻¹VᵀPAV — has every pole inside the numerical range of A in the P inner product, which is in the left half plane. A stable A always has such a P: solve the Lyapunov equation AᵀP + PA = −I.

The weighted projections are stable at every rescaling measured. At t = −3, where the plain order-1 models reach +3.04 and order-2 models +19.1, the weighted ones reach −0.28 and −6.12; at t = −2, −0.81 and −5.39. In the original coordinates the weighted and plain models are close, −4.72 against −5.17 at order 1, so the weighting costs little where it is not needed.

The Lyapunov solution is exactly the object why a Gramian can be truncated at all built balanced truncation’s stability from: a positive definite solution of a Lyapunov equation whose leading block certifies the truncated system. A one-sided projection in that inner product is the same argument with the truncation replaced by a projection onto an arbitrary subspace. It costs a dense Lyapunov solve, which is most of the cost of balanced truncation, and so it is not a cheap route to a certificate. What it shows is where the certificate lives: in an inner product in which the state matrix is dissipative, and the Euclidean inner product of whatever coordinates arrived is one such product only by accident.

What each repair costs #

Three weights have now been measured to restore the certificate, and they are priced very differently.

The symmetrising scaling is a vector computed from the ratios of the off-diagonals. It costs nothing, gives the spectrum’s own margin, and exists only for operators with that structure — a tridiagonal matrix with off-diagonals of one sign, or more generally a matrix similar to a symmetric one by a diagonal. Its conditioning grows geometrically with the number of states and the convection, so on a fine grid at high Péclet number the basis it produces carries entries that differ by more than a double can hold, and the margin it promises cannot be computed in the arithmetic it is promised in.

The Lyapunov solution is a dense matrix equation, and on a large system it is exactly the cost the interpolatory constructions were adopted to avoid. Two approximants and one matrix size priced that difference as the reason to prefer interpolation at all. A code that pays for a Lyapunov solve to certify a Galerkin model has paid most of the price of balanced truncation and received a weaker model.

The Euclidean inner product of well-chosen coordinates costs nothing and certifies with a smaller margin, and the measurement here says what well-chosen means: coordinates in which the symmetric part of the state matrix is negative definite. For a discretised physical operator those are usually the coordinates of its energy norm, which a code knows before it builds anything. The failure measured here is what happens when the state is stored in some other coordinates and the projection is done in them anyway.

Two limits of the measurement #

A positive edge in the other direction did not produce a failure. For t = 2.5 and 3 at Péclet 10 the edge is +2.16 and +6.59, the same values as at −0.5 and −1, and no one-sided model at any order was unstable there. The skew in that direction makes the upper off-diagonal the large one; the rational Krylov bases built from the input vector apparently do not align with the directions in which xᵀAx is positive. This measurement does not explain the asymmetry, and it is a reminder that the certificate is a sufficient condition: removing it permits failure and does not guarantee it.

The Lyapunov solver has limits of its own. At Péclet 25 in strongly rescaled coordinates, and above about 50 in any, the real Schur form of the state matrix acquires complex pairs that the Bartels–Stewart Lyapunov solver used here declines, so the weighted projection is drawn only at Péclet 10. The argument holds wherever P exists; the measurement is at one Péclet number.

What a code that reduces models should take from this #

A one-sided projection is certified only in the inner product it is computed in. A code that builds a Galerkin model and relies on its stability should either compute the numerical range’s edge in the coordinates it projects in — a symmetric eigenvalue problem, and a cheap one — or project in an inner product that is known to make the state matrix dissipative: the energy inner product of the underlying physics, a symmetrising scaling when one exists, or a Lyapunov solution when neither is known.

Nothing about this is special to reduction. A model that cannot be run traced the two-sided failure to the same non-normality that the direction the error leans found in iterative solves and a limit the matrix never reaches found in powers of the matrix; the numerical range is the quantity that governs all three transiently, and its dependence on the inner product is why the same operator can look tame in one basis and violent in another. Which basis a code works in is a decision about stability whether or not the solver knows it has made one.

The claim that has to fail #

The claim is the natural reading of the previous measurement: that one-sided projection of a stable system is stable, full stop. At Péclet number 10, in coordinates rescaled by a diagonal matrix of condition number 1.3·10⁴, the system’s eigenvalues and transfer function are unchanged and one-sided models of order 1 are unstable at eleven placements of nineteen. The refusal is fed the claim that none is and fails.

Still open: an inner product the physics already has, and a cheaper dissipative weight #

The energy inner product. A convection–diffusion equation has a natural norm — the $L^2$ norm of the continuous field — in which the diffusion term is dissipative and the convection term conserves energy. On a uniform grid that norm is the Euclidean one scaled by the mesh width, which is why the original coordinates happened to be safe. On a graded mesh the energy norm weights each node by its cell size, and projecting in it rather than in the Euclidean norm of the nodal vector is the natural repair. How much grading it takes before the Euclidean projection fails, and whether the mass-weighted projection then keeps the certificate at no extra cost, is unmeasured.

A diagonal weight found by an iteration. The Lyapunov solution is dense and expensive. A diagonal positive P with AᵀP + PA negative definite exists for this operator — the symmetrising scaling supplies one — and for many operators with a sign structure like it, can be searched for by a few steps of a simple iteration where no closed form is known, and would give a certificate for the cost of a vector. Whether one exists at high Péclet number, and how its margin compares with the symmetrising scaling’s, is the natural next measurement.