A model that is a rational function

Every object this site has drawn so far has been a matrix, a matrix polynomial or a tensor, and the question asked of it has been solve or decompose. A linear time-invariant model is none of those. It is

ẋ = Ax + Bu,   y = Cx,

and what anybody actually wants from it is not x. It is the map from the input to the output, which in the Laplace variable is

H(s) = C(sI − A)⁻¹B.

A is n × n, and n is whatever the discretisation happened to produce — a hundred thousand for a finite-element model of a beam, or a few million for a circuit. H is a scalar function of one complex variable. The gap between those two numbers is the whole of this field, and the first thing worth establishing is that it is a real gap rather than an artefact of how the model was written down.

Two routes, so the picture means something #

The curve above is computed by factorising a complex matrix sI − A at every frequency, which is what a code does. The crosses are a closed form: the state matrix here is the second-difference operator this site has used since its foundation phase, so its eigenvalues are 4(n+1)²sin²(kπ/2(n+1)) and its eigenvectors are the discrete sine basis, both written down rather than computed. In that basis

H(s) = Σₖ (cᵀvₖ)(vₖᵀb) / (s − λₖ),

a sum of n terms with no factorisation anywhere. The two routes share the coefficient matrices and no code, and they agree to 8·10⁻¹⁵ across the range where the comparison means what it looks like. That is the standing an answer that is known insists on, and this field needs it more than most: almost every claim below is about how far an approximation is from an answer, and an approximation compared against a better approximation is not a measurement.

The qualification is worth a paragraph, because it is the first thing this model teaches. Past about ω = 10³ the two routes stop agreeing relatively — 1.4·10⁻⁹ at the top of the range — and neither of them is wrong. The residues of this system sum to cᵀb, and the sensor is not at the actuator, so they sum to zero: the leading 1/s term of the modal expansion cancels exactly, and the sum is computing a number fifteen orders below its own largest term. Measured against the size of the terms being added, which is the backward question rather than the forward one, both routes sit at 3·10⁻¹⁵ at every frequency tried. It is a small residual is not a small error read backwards: here the comparison is the thing losing digits, and the computation is fine.

The number that is the complexity #

Refine the grid and the curve above barely moves. That is not a coincidence and it is not smoothness. The modes a finer grid adds are stiff — their λ is of order n² — and an actuator with a physical extent hardly excites them, so their contribution divides by a large number twice over. The function is a property of the problem; n is a property of the discretisation.

The quantity that says how much of n was ever the complexity is the McMillan degree: the number of poles H actually has. It is a property of the input–output map and not of the coordinates, and it can be arbitrarily far below n. lowDegreeModel makes the gap explicit — the state matrix is the same second-difference operator, and the actuator and sensor are combinations of the first r sine modes only, so the remaining n − r modes are neither driven nor seen and H has exactly r poles however large n is.

That is not a trick played on the measurement. It is the ordinary situation. An actuator with a physical extent cannot excite a mode whose wavelength is shorter than it is, and it is the reason a model of a hundred thousand states can have a transfer function that six poles reproduce.

The definition, and the algorithm #

Every derivation defines the degree the same way: it is the rank of the Hankel matrix of Markov parameters CAᵏB, the Taylor coefficients of H at infinity. The definition is correct. As an algorithm it does not survive contact with floating point, and the way it fails is the kind this site keeps finding.

Two things are visible and neither is about the transfer function. The first ten parameters are exactly zero — not small, zero, with no rounding in them — because A is a three-point stencil and information takes one step per grid point to travel the ten points from the actuator to the sensor. Then they grow like ‖A‖ᵏ with ‖A‖ = (n+1)² = 625, reaching a ratio of 1.9·10³⁷ across the run.

A Hankel matrix built from those numbers has a dynamic range no rank decision can see through. Rank is a decision about where to place a gap, and a matrix whose entries span thirty-seven orders has no gap to place: the answer comes back two where it is six.

The route with less information #

The other route reads samples of H and nothing else. Pick two interleaved sets of points on the imaginary axis and form the Loewner matrix

L[i][j] = (H(μᵢ) − H(λⱼ)) / (μᵢ − λⱼ),

a divided difference of the function between two points. Its rank is the McMillan degree of whatever rational function produced the samples. A is never touched. The samples could have come from an instrument.

Measured on the same model, that rank is six, and the sixth singular value stands 10⁵ above the seventh at n = 20 and 2.9·10⁴ at n = 80. That is not a close call, and calling it a rank decision is not a hedge — it is a decision with a cliff to make it on, which is exactly what the cheap rank and what it cannot see says a rank needs.

So there are two routes to one integer. The first reads A, B and C — strictly more information — and returns the wrong answer. The second reads samples of a function and returns the right one with four orders to spare. The definition is not the algorithm, and this is the field’s first instance of a pattern that will recur: the object a derivation names and the object worth computing are different objects, and knowing more can be worse.

Why the Loewner route is well conditioned and the other is not #

Both are ranks of matrices built from the same function, so the difference has to be in what the entries are made of.

A Markov parameter is a coefficient of the expansion of H at infinity, and that expansion has a radius of convergence set by the largest pole. Its coefficients therefore grow geometrically at a rate given by ‖A‖ — which is a property of the grid, not of the function, since a finer discretisation of the same operator has a larger norm. Every one of the six directions the Hankel matrix ought to carry is present in it; they are just scaled apart by thirty orders, and a numerical rank is blind to a direction thirty orders below the largest.

A Loewner entry is a divided difference of H between two points where H is of order one. Nothing in it grows. The matrix has six directions of size 10⁻¹ to 10⁻¹⁰ and the rest at rounding, which is precisely the shape a rank decision can read.

This is the same argument as the road that squares the problem: two constructions that are algebraically equivalent, one of which puts a large number where the other does not, and the arithmetic is decided entirely by that. It is worth noticing that here the losing construction does not square anything — it exponentiates, which is worse, and the essay’s lesson generalises past its own algebra.

What a rank of six is worth knowing #

A degree is not by itself a reduced model. What it is, is the answer to the question this field opens on: how large would a model have to be to reproduce this function? Six, on a system with twenty-four states or eighty; and the number is available from samples, before any reduction is attempted, without the state matrix.

That has three consequences worth carrying into the rest of the field, and it is the same shape as an eigenvalue count that cannot be slightly wrong: an integer, computed by a route that cannot be slightly wrong about it, standing in for a measurement that can.

A reduced model of order below the degree cannot be exact, and one at the degree can. So the whole question of how wrong is a reduction only becomes interesting below r = 6, and above it the answer is “not at all, up to rounding”. The essays that follow are about the region below.

The degree bounds nothing about how easy the reduction is. A function with six poles spread over four decades and a function with six poles in a cluster have the same degree and very different behaviour under approximation. What decides that is the decay of a different sequence of numbers entirely — the same distinction the best approximation there is draws between a rank and the singular values below it — and the bound that is known in advance is the essay about them.

The exact degree is almost never the number wanted. Real models have no exact degree: their Loewner singular values decay rather than falling off a cliff, and where the cut falls is a decision with the same shape as every other rank decision on this site. The clean cliff drawn above is a property of a model built to have one, and it is drawn that way so that the comparison between the two routes is not itself a judgement call.

The sampled route, and what it costs #

Nothing above is free. The Loewner matrix needs 2p samples of H, and each sample is a factorisation of a complex n × n matrix — so on a model with a hundred thousand states, thirty-two samples is thirty-two large sparse solves. That is not nothing, and it is exactly the same operation every method in this field is built out of: a shifted solve, which is what restarting is a filter spends its budget on for a different reason entirely.

What it is not is a function of the degree. The number of samples decides how large a degree the rank can resolve, and the samples themselves cost what they cost whatever comes back. So a code that wants to know whether a model is worth reducing pays thirty-two solves to find out, and if the answer is six it has also — as exact at the points that were named will show — nearly paid for the reduced model itself.

Where the frequencies have to be #

One decision is buried in the Loewner construction and it is the same decision this whole phase turns on: where the samples are taken.

The rank result is exact for any distinct points, in exact arithmetic. In floating point the sampling has to cover the region where H varies, or the divided differences are all nearly equal and the matrix is numerically rank one. Sampled on the real axis over the same range, this model’s Loewner singular values decay by a factor of about fifty per index rather than dropping off a cliff, and the sixth-to-seventh gap falls from 10⁵ to about 20.

The reason is worth stating because it recurs: on the imaginary axis the samples sit between the poles and each one sees a different combination of them; on the positive real axis they all sit on the same side of every pole and the resulting matrix is Cauchy-like and graded. Same theorem, same number of samples, and one is a measurement while the other is a judgement call.

The three counts, and which of them is a property of the problem #

It is worth listing them together, because a reader arriving from another field will expect the first to matter and it is the only one that does not.

• n, the state dimension. A property of the discretisation. Doubling the grid doubles it and changes the function in the fifth digit.
• The McMillan degree. A property of the input–output map. It is what a reduced model has to carry, and it is invariant under any change of coordinates.
• The number of samples. A property of the question being asked. It sets how large a degree can be resolved and nothing else.

Every remaining essay in this field is about a fourth number that is none of the three: the order r at which somebody chooses to truncate, which is smaller than the degree and is a decision about accuracy rather than about the model. What makes that decision possible to make well is that, uniquely in this field, there is a bound on the error before the reduced model exists.

One more thing the degree is not #

A reader arriving from the sparsity field will want to know whether a low McMillan degree means the state matrix has low rank, and it does not. A is the second-difference operator: full rank, well conditioned relative to nothing, and entirely ordinary. What is low rank is the Hankel operator that maps past inputs to future outputs, which is an object built from A, B and C together and is nowhere in the model as written.

That is worth one sentence of caution because it is the natural analogy and it points the wrong way. Every other kind of smallness on this site — sparsity, off-diagonal rank, a tensor format — is a property of a matrix’s entries or its blocks, visible by looking at it. This one is a property of a triple, invisible in any of the three, and computable only by asking the triple a question.

At other settings #