A condition number that is not the model's

A reduction method that begins with two Lyapunov solves ends up quoting two condition numbers, and they are quoted for one purpose and read for another.

The purpose is settled, and the product nobody had to form settles it: the Hankel singular values can be obtained by eigendecomposing PQ or by taking the singular values of RᵀS, the two are the same algebra, and one of them works at κ(P)κ(Q) while the other works at its square root. On a twenty-two-state model those are 2.72·10³⁶ and 1.65·10¹⁸, and only the second is on the near side of 1/u. That is a statement about two algorithms, it is correct, and nothing below disturbs it.

The reading is the sentence that follows it in most accounts, and it is not the same sentence. If κ(P) is enormous, the Gramian is nearly singular; a Gramian is nearly singular exactly when its eigenvalues have fallen off a cliff; and the cliff is what makes reduction possible. So κ(P) looks like a measure of reducibility with a large number attached — a model that reduces well should have a huge κ(P), and a model that reduces badly should have a modest one.

The two quantities that decide whether that is true are the numerator and the denominator of κ(P) separately. Both were measured across a sweep in which the model’s McMillan degree runs from four to fourteen at a fixed state dimension of twenty-two, so the thing that is supposed to move is moving by a factor of three and a half and nothing else is moving at all.

The numerator is the same number at every degree #

‖P‖₂ across the sweep is 0.1419552374, 0.1419605699, 0.1419623305, 0.1419627089, 0.1419627168, 0.1419627175 and 0.1419627175, at degrees four, five, six, eight, ten, twelve and fourteen. To four significant figures it is 0.1420 at every one of them. The whole movement, from one end of the sweep to the other, is 5.3·10⁻⁵ relative, and it is monotone and converging rather than wandering.

The reason is structural rather than a coincidence of this family. The model is built by driving the first r modes of a fixed twenty-two-mode operator with weights that fall geometrically, so the input vector’s component on the leading mode is the same number whatever r is: b̃₁ = 1.6297632, against a leading eigenvalue of −9.8542691. In modal coordinates the Gramian’s leading diagonal entry is −b̃₁²/(2λ₁) = 0.1347704, which is 94.9% of ‖P‖₂, and it is that quantity at every degree in the sweep. Raising the degree adds modes further down the geometric weighting, and they contribute to the Gramian’s smaller directions rather than to its largest one.

So the largest Gramian direction comes from the first excited mode, and that mode is the same throughout. Nothing a caller would call reducibility can be visible in the numerator, because the numerator has not been told about it.

The denominator is a rounding level #

The Gramian of a model of McMillan degree r has exact rank r. The input matrix is a combination of r eigenvectors, the right-hand side BBᵀ has rank one inside that r-dimensional span, and the solution of AP + PAᵀ + BBᵀ = 0 is supported on it. In exact arithmetic P has r non-zero singular values and twenty-two minus r zeros, and κ(P) is infinite.

It is not infinite when computed, because the zeros arrive as whatever the Bartels–Stewart solve left behind. Divided through by the largest, the singular values of P at degree four are

1.00  3.13·10⁻²  1.17·10⁻³  1.84·10⁻⁵  8.64·10⁻¹⁷  7.73·10⁻¹⁷  6.71·10⁻¹⁷  …  3.72·10⁻¹⁹

— four values spanning five decades, then a drop of twelve orders, then eighteen values that do not descend so much as mill about between 10⁻¹⁶ and 10⁻¹⁹. That second group is not a tail of the Gramian. It is the arithmetic, seen through a matrix that has nothing there.

σ_min(P) over the sweep is 5.27·10⁻²⁰, 5.23·10⁻²⁰, 2.51·10⁻²⁰, 1.35·10⁻²⁰, 1.13·10⁻¹⁹, 1.80·10⁻¹⁹ and 5.15·10⁻²¹. Every one of them is below ‖P‖₂u = 3.15·10⁻¹⁷, which is the level at which a backward-stable solve is entitled to leave something, and none of them is anywhere near the smallest genuine singular value — 1.84·10⁻⁵ relative at degree four, and 4.23·10⁻¹⁵ relative at degree fourteen, which is the lowest a real Gramian direction gets anywhere in the sweep.

κ(P) is therefore ‖P‖₂ divided by the smallest member of a bag of rounding errors. It is not a ratio of two things the model has. It is one thing the model has, divided by one thing the solver did, and the second is what a rounding error is: a number with no reason to take any particular value.

The sweep wanders, and its largest value is at the wrong end #

Across degrees four, five, six, eight, ten, twelve and fourteen, κ(P) is measured at

2.69·10¹⁸   2.72·10¹⁸   5.66·10¹⁸   1.05·10¹⁹   1.26·10¹⁸   7.90·10¹⁷   2.76·10¹⁹

The largest is 34.9 times the smallest, and since the numerator is fixed to five decimal places, that factor of thirty-five is entirely the factor of thirty-five between the largest and smallest rounding level: 1.80·10⁻¹⁹ against 5.15·10⁻²¹. The sequence rises, rises, rises, falls, falls, and then jumps by a factor of thirty-five in one step.

Two features of that sequence are worth separating, because only one of them is surprising.

There is no trend. A trend would be the reading the sentence at the top of this essay predicts: a model that reduces well should show a large κ(P) and one that reduces badly a small one. What the measurement shows is a sequence with no direction in it at all, whose largest value, 2.76·10¹⁹, sits at degree fourteen — the degree at which this family is least reducible. A model of McMillan degree fourteen is exactly reproduced by fourteen states and no fewer, against four at the other end of the sweep, so the ordering of reducibility across the seven is not in doubt and is not a matter of which tolerance is chosen. The reading is not merely imprecise. On this sweep it points the wrong way at the extreme.

The obvious repair — that a rounding level is at least a statistic, so its minimum should trend with how many of them there are — does not survive being counted either. The number of singular values of P below ‖P‖₂u is 18, 17, 16, 14, 12, 11 and 11 across the sweep, because raising the degree converts a rounding-level direction into a real one. A minimum taken over eighteen draws should sit lower than one taken over eleven, and the measurement runs the other way at both ends: 5.27·10⁻²⁰ out of eighteen at degree four, and 5.15·10⁻²¹ out of eleven at degree fourteen.

And the size is the same everywhere. All seven values are between 10¹⁷ and 10¹⁹, which is 1/u to within two orders. That is what a norm divided by a rounding level has to look like, and it is why the number is so easy to over-read: a κ of 10¹⁸ is genuinely enormous, it is genuinely what makes the product route fail, and it says exactly nothing about which of these seven models it came from.

The comparison the slider does not disturb is the one the figure is about. κ(P)κ(Q) is 1.11·10³⁷, 3.51·10³⁶, 3.23·10³⁷, 7.42·10³⁷, 2.72·10³⁶, 5.88·10³⁵ and 3.59·10³⁷ across the seven degrees. The smallest of those is 1.3·10²⁰ times 1/u. The product route is past the point where nothing small survives at every degree in the family, reducible or not, and the conclusion drawn from it — do not form PQ — is correct at every degree for the same reason: because it is a fact about squaring a rounding level, not about the model being squared.

One model, solved twice, has two condition numbers #

There is a sharper version of this available, because this family has a second route to P. The state matrix has a known eigenbasis, so in modal coordinates the Lyapunov equation is n² scalar divisions, P̃ᵢⱼ = −b̃ᵢb̃ⱼ/(λᵢ + λⱼ), and evaluating it is a formula rather than a solve. That is the closed form the answer is checked against in this field, and here it can be turned into an instrument.

Both routes were run at every degree. The two norms agree to ten significant figures at every one of the seven, printing the same digits. The condition numbers do not. The closed form gives 4.80·10¹⁸, 2.36·10¹⁹, 7.84·10¹⁸, 4.82·10¹⁸, 1.55·10¹⁸, 5.83·10¹⁷ and 5.55·10¹⁸, and against the solve those are ratios of 1.78, 8.71, 1.38, 0.458, 1.23, 0.737 and 0.201 — a spread of a factor of forty-three between the two extremes.

Nothing about the model changed between those two columns. The system matrix, the input, the output and the exact Gramian are identical; the arithmetic that produced the small directions is different, and κ(P) moved by up to a factor of nine. A quantity that changes by nine when only the route to it changes is a quantity about the route.

There is a second disagreement inside the same run, and it needs no alternative route at all. A model has one reducibility and two Gramians. The ratio κ(P)/κ(Q) across the sweep is 0.653, 2.10, 0.992, 1.49, 0.583, 1.063 and 21.2 — so at degree fourteen the controllability Gramian and the observability Gramian of one model, computed by one routine in one run, disagree by a factor of twenty-one about how ill-conditioned that model is. The two matrices have the same exact rank, the same driven modes and the same twenty-two-dimensional space to be rounded in; what differs is which rounding each solve happened to leave behind. A property of the model would have to be one number.

That is the same separation two condition numbers of one matrix makes for a linear system — one object, two measurements of how hard it is, and the question of which one governs decided by what the algorithm actually forms. The difference here is that neither of these two is a property of the object at all. Both are properties of a solve, and the model has no condition number in this sense to have.

A model with nothing missing reads the same number #

The sweep holds the family fixed and moves the degree, which leaves one objection open: every model in it has a Gramian of exact rank r, so perhaps κ(P) is empty only for models built that way. A model whose Gramian is not rank-deficient by construction settles that.

The heat model this family is built on top of is one. Its input excites all twenty-two modes, so P has no exact zeros at all and its singular values decay because the modes decay rather than because they are absent. Measured at twenty-two states: ‖P‖₂ = 2.00050, σ_min(P) = 2.30·10⁻¹⁸ against ‖P‖₂u = 4.44·10⁻¹⁶, a numerical rank of thirteen, and κ(P) = 8.72·10¹⁷.

That is the whole comparison in one line. The heat model needs eleven states to be reproduced to nine digits and the degree-four model needs four; the heat model’s Gramian carries thirteen resolvable directions and the degree-four model’s carries four; and the heat model reads a κ(P) three times smaller. Nothing was missing by construction in the first and eighteen directions were missing outright from the second, and their condition numbers agree to within a factor of three, in the wrong order.

The mechanism is the one already named, and it survives losing the exact rank. σ_min is still under ‖P‖₂u, so the smallest computed direction is still rounding rather than model. An exactly rank-deficient Gramian makes that vivid; a smoothly decaying one arrives at the same place anyway, because a decay of twelve orders reaches the rounding level well before it reaches the twenty-second mode.

What does report reducibility #

The quantity that answers the question κ(P) is being asked is the Hankel spectrum, and it is available from the same computation.

Over the same seven degrees, σᵣ/σ₁ — the smallest Hankel singular value the model genuinely has, relative to the largest — is 3.04·10⁻⁵, 4.58·10⁻⁷, 1.92·10⁻⁸, 2.45·10⁻¹¹, 1.09·10⁻¹⁴, 2.63·10⁻¹⁷ and 1.23·10⁻¹⁷. That falls monotonically across twelve orders of magnitude, on the same models, over the same slider, in the same run that produced the seven wandering values above.

It reads as a bound, too, which is the form a caller can act on. The smallest order whose truncation error bound is under 10⁻⁹ relative is four, five, six, seven, seven, seven and seven across the sweep — rising with the degree until the arithmetic runs out, which is the behaviour the condition number was being asked for and does not have. And it is the quantity the bound that is known in advance is a sum of, so reading the σ curve is not an extra computation but the one the method already performs.

Where both instruments stop #

The honest version of the sweep needs one more column, because the σ curve is not an unlimited instrument either and the point at which it stops is measurable.

The Gramians agree with their closed form to between 2.14·10⁻¹³ and 2.15·10⁻¹³ relative at every degree — which is very good and is not sixteen digits. Below that level a Hankel singular value is not a Hankel singular value. The numerical rank of P at a 10⁻¹³ threshold is 4, 5, 6, 8, 9, 10 and 10 against degrees 4, 5, 6, 8, 10, 12 and 14, so a model whose McMillan degree is fourteen has ten directions in double precision and the other four are not there to be kept. That is where the last two entries of the σᵣ/σ₁ list come from: 2.63·10⁻¹⁷ and 1.23·10⁻¹⁷ are below the level at which the Gramians were computed, and they are the same kind of number as σ_min(P) — a rounding, printed to three figures.

So the sweep has three regimes and only the middle one is a measurement of the model. Below degree ten every Hankel singular value the model has is resolvable, and the σ curve reports its degree exactly. At twelve and fourteen the model’s own decay has fallen through the accuracy of the solve, and the instrument reports ten in both cases. Deciding where that happens is a decision backed by a gap rather than a property read off the matrix, and the gap here is the twelve orders between 4.23·10⁻¹⁵ and 7.79·10⁻¹⁷ in the degree-fourteen Gramian’s own spectrum.

The refusal this essay is published with is that same fact at the far end. A twenty-state model of McMillan degree four, asked for a reduction of order nineteen, has σ₁₉ = 0 exactly, and the assertion that the truncation sits above the numerical rank of the Gramians is fed that case and rejects it. The states are not there. No condition number, however large, is a warning that they might be — a zero has to be decided to have arrived, and the decision is made against a gap in the σ, not against κ.

What follows for anything that reads a Gramian #

Do not report κ(P) as a property of the model. It is a property of the solve that produced P, it moves by a factor of forty-three between two routes to the identical matrix, and its size — 10¹⁷ to 10¹⁹ throughout — is set by u and by ‖P‖₂ rather than by anything the model does. The same caution an amplifier earns applies here in reverse: κ is a measurement of a problem when the problem is what is being perturbed, and a numerically singular Gramian is not that problem.

Do report κ(P)κ(Q), for the one thing it says. It decides between two routes to the σ, it is past 1/u at every degree measured here, and that is the sentence it supports. It supports no sentence about which of the seven models it came from.

Read reducibility off the σ, and say which order the reading is good to. σᵣ/σ₁ falls twelve orders across this sweep where κ(P) wanders by a factor of thirty-five in no direction, and it stops being a reading at the accuracy of the Lyapunov solves. Both halves are needed: a σ curve quoted without the level it flattens at invites exactly the over-reading this essay is about, one index further down.

And do not confuse this with a Gramian’s decay having no cause. It has a very definite one, and why a Gramian can be truncated at all derives the rate in closed form from the spectral interval. The decay is a property of the model; the condition number of the computed Gramian is not a way of seeing it. The first is the cliff, the second is the distance from the top of the cliff to the noise on the floor, and only the first is about the geology.

A model can be barely reducible and still numerically singular. That is the combination degree fourteen is, and the two facts are independent: it needs fourteen states for an exact reproduction and seven for nine digits, and its Gramian has a numerical rank of ten out of twenty-two. A condition number cannot hold both, which is why it holds neither, and it is the same failure of one scalar to carry a spectrum that a determinant commits at the other end of the site.

And the smallest singular value needs a scale before it is small. σ_min(P) = 5.15·10⁻²¹ is a number with no meaning until it is put beside ‖P‖₂u = 3.15·10⁻¹⁷, at which point it says the solve was clean rather than that the matrix was singular — which is what small compared to what is about, and κ(P) is precisely the ratio that supplies the wrong scale.

The general shape is one this collection meets often enough to name. A quantity is defined correctly, computed correctly, and quoted correctly for one purpose; a second purpose is inferred from its size rather than from its definition; and the inference survives because the number is large in both readings and nobody sweeps the parameter that would separate them. Sweeping it costs seven Lyapunov solves here, and the sweep is what the seven wandering values are.