A model with no matrices behind it

Every reduction in this field so far starts from (A, B, C) and produces a smaller (A, B, C). The state matrix is the input, the projection is computed from it, and the reduced model is a shadow of an object somebody already had.

The Loewner realisation does not start from a state matrix. It is handed values of the transfer function at a set of points — the output of an instrument, or of a simulation nobody can open — and returns a descriptor system that interpolates every one of them exactly.

The construction is two matrices of divided differences #

Split the samples into a left set {μᵢ, hᵢ} and a right set {λⱼ, gⱼ}, and form

Lᵢⱼ = (hᵢ − gⱼ)/(μᵢ − λⱼ)        Lsᵢⱼ = (μᵢhᵢ − λⱼgⱼ)/(μᵢ − λⱼ)

Then (E, A, B, C) = (−L, −Ls, h, gᵀ) is a descriptor system whose transfer function passes through every sample. That is the whole of it: no factorisation of anything, no state, no model of the physics. The data is the model.

The verification is direct. On a system whose McMillan degree is genuinely three, five or eight — a forty-state heat model whose actuator and sensor excite only the first r modes — the model of order r matches every sample to:

true degree   samples   model at r − 1   model at r      at points it never saw
    3            24        2.2·10²        1.2·10⁻¹¹          3.4·10⁻¹³
    5            24        1.2·10¹        8.7·10⁻¹⁴          2.0·10⁻¹¹
    8            24        3.1·10⁻³       5.8·10⁻¹⁴          2.4·10⁻¹⁰

Interpolation is exact by construction rather than by accuracy. At the right order the model passes through the data at the rounding level; one order lower it does not, by fourteen orders in the first row. And the last column is the part that matters: the model also matches the function at points it was never shown, to 10⁻¹⁰ or better.

The degree is a rank, and the rank has a cliff #

The order to project onto is not a free choice. It is the rank of the Loewner matrix, and on clean data that rank is not a decision at all — the singular values fall off a cliff at the McMillan degree and the gap across it is:

true degree    gap at the cut
    3            1.1·10¹⁴
    5            7.6·10¹⁴
    8            5.4·10¹²

Twelve to fourteen orders of separation. Rank is a decision is this site’s argument that a numerical rank is a threshold rather than a property — and here the threshold has fourteen empty decades to sit in, so any cut anywhere returns the same integer.

That is the same measurement the degree-detection essay makes about the pencil. What is new is the step after it: building the model, and confirming that order r interpolates while r − 1 does not.

And then the data has error in it #

Real samples do. An instrument has a tolerance, a simulation has a solver tolerance, and a measured frequency response has both.

Adding relative noise to the twenty-four samples and repeating everything:

noise      degree read   gap at the cut   best model     at order
none            6           2.0·10⁸        2.7·10⁻¹¹        12
10⁻¹⁰           9           4.4            1.6·10⁻⁷         11
10⁻⁶           10           6.5·10³        1.6·10⁻³         11
10⁻³           12           —              1.9              4

The gap is not monotone in the noise, which is the trap. Across noise of 10⁻¹², 10⁻¹⁰, 10⁻⁸ and 10⁻⁴ the gap at the cut runs 4.73·10⁶, 4.41, 6,561 and ∞ — collapsing by six orders and then recovering by three — while the degree read runs 6, 9, 10 and 12, climbing monotonically the whole way.

So a healthy-looking gap is not evidence of a correct degree. At 10⁻⁸ a caller inspecting the singular values would find a clean separation of three and a half orders — better than many honest problems offer — and read a degree of ten from a system whose degree is six. The gap recovered because the noise floor rose past the next group of singular values and made a new cliff, and a new cliff between the wrong pair of values looks exactly like the right one. The one diagnostic available is the quantity that recovers while the answer keeps degrading, which is worse than a diagnostic that simply fails.

Three things fail at once and they are worth separating, because only one of them is obvious.

The rank decision loses its cliff. At 10⁻¹⁰ noise the gap falls from 2·10⁸ to 4.4. The singular values below the degree stop falling to the rounding level and sit at the noise level instead, so there is no longer an empty region for a threshold to live in and the degree becomes a choice. At 10⁻³ there is no cut at all — every singular value is above any sensible threshold, and the routine reports the full order.

The interpolation is still exact, and that is the problem. The construction passes through the samples it was given, so it passes through the corrupted ones. A model that reproduces its data perfectly is not evidence of anything when the data is wrong, and the interpolation error — the one quantity available without knowing the answer — reports success throughout.

And the approximation is floored by the noise, at this sample count. No order does better than 1.6·10⁻⁷ at noise 10⁻¹⁰, 1.6·10⁻³ at 10⁻⁶, and 1.9 — worthless — at 10⁻³. Taking more of the degree past the best order buys nothing, and the best order is not the true degree. That floor turns out to be a property of the twenty-four samples rather than of the noise, which the section on oversampling below measures.

The projection, and why it is two subspaces rather than one #

At more samples than the degree the pencil (L, Ls) is singular — its rank is the degree, and a pencil with a rank deficiency has no well-defined transfer function. The construction therefore projects, and the projection is not the usual one.

Two subspaces are needed, not one. Y is taken from the left singular vectors of the wide matrix [L Ls], and X from the right singular vectors of the tall matrix [L; Ls], each truncated to k columns. The reduced pencil is then (−YᵀLX, −YᵀLsX) with B = YᵀV and C = WX.

The asymmetry is the point. A Galerkin projection uses one subspace on both sides and is what balanced truncation does; a Petrov–Galerkin projection uses two, and it is the two-sided structure that makes the interpolation exact rather than approximate. The same distinction is what two-sided interpolatory projection is about, and it comes with the same hazard: a two-sided projection does not inherit stability, so a reduced model built this way can have a pole in the right half-plane even when the system it came from does not.

That hazard is measured in a model that cannot be run for the interpolatory case, and it applies here unchanged. Nothing in this construction checks that the model it returns is stable, and on data from a non-symmetric system there is no reason it should be.

Why the failure has this particular shape #

The three failures above are one failure seen three ways, and the mechanism is the same one the adaptive approximant has.

A construction that interpolates exactly has no freedom to disagree with its data. Whatever error the data carries is reproduced in the model, at full size, and the model’s own error against the truth is therefore bounded below by the data’s — which is the floor in the last column. Nothing about the linear algebra is at fault; the model is the exact answer to the question it was asked.

What makes this worse than an ordinary noise floor is that the diagnostic available from inside the computation is the one that does not see it. The interpolation error is computable without knowing the true system, and it says 10⁻¹⁴ throughout. The approximation error is what matters and is not computable without the answer.

So a Loewner realisation from noisy data reports perfect success and delivers whatever the noise allows, and the only way to know which is to hold samples back — fit on some, test on the rest, which is cross-validation and is not part of the construction as usually stated.

What the second route is #

The site’s habit is a second route, and there is one here that is unusually strong because it does not share the construction at all.

The reduced model is built from samples. The system it approximates is available as a state-space model, so its transfer function can be computed at any point by a complex solve of sI − A. The comparison in the last column of the first table is between two objects whose only common ancestor is the physical system: a model built from twenty-four numbers, and a forty-state matrix nobody handed the construction.

They agree to 2.1·10⁻⁹ on the low-degree family and to 1.3·10⁻¹³ on a modal system whose transfer function has a closed form. That is not a check of an implementation against itself; it is a check that a construction consuming only data has recovered an object it was never shown.

What this makes possible #

The reason the field bothers is that (A, B, C) is often not available.

A measured system. A frequency response from a network analyser is a table of complex numbers. There is no state matrix to project, and every method in this field before this one requires one.

A simulation nobody owns. A commercial solver, a legacy code, a coupled multiphysics run. Each evaluation is expensive and opaque, and the interface is “give me a frequency, get a number”.

A system too large to hold. A discretisation with millions of states has a transfer function that can be evaluated by one sparse solve per point, and that is the only operation it supports at a sensible cost.

In all three, the number of evaluations is the budget, and the Loewner construction spends it once — 2r samples for a model of order r, with no iteration and no restarts. That is a different economy from the interpolatory methods that choose their points iteratively, which need a good initial guess and converge, and it is why the two coexist.

What the samples cost, and where to put them #

The construction is linear in the data and it is not indifferent to where the data is.

The measurements here use points spread geometrically along the negative real axis, well away from the spectrum, which keeps every quantity real and makes the arithmetic checkable against the site’s own routines. Real practice samples the imaginary axis in conjugate pairs, which is the same construction seen through a unitary change of basis and produces real matrices by the same argument.

What matters in either case is coverage. The model interpolates at its samples exactly and approximates between them, so a region with no samples is a region the model is guessing in. The error curve in the figure is measured at forty points between the samples for exactly that reason: an error measured only at the samples is measuring the construction’s definition rather than its quality.

The number of samples must also be at least twice the degree, and the routine refuses an odd count because the pencil cannot be split.

What oversampling does, which is two opposite things #

The natural claim about extra samples is that they are free and useful — that they do not raise the order, since the rank is still the degree, and that they make the rank decision better conditioned, which on noisy data would be the one lever available. Sweeping the sample count at three noise levels on the degree-six system says both halves are true of clean data and the second is backwards on noisy data, which is the only case the advice is for:

noise      samples   gap at the cut   degree read   error at order 6
none          16       4.3·10¹⁵            6           6.9·10⁻¹²
none          32       2.9·10¹¹            6           6.4·10⁻¹¹
none          64       1.0·10¹¹            6           6.8·10⁻¹²
10⁻⁸          16       3.0·10⁸             6           3.5·10⁻⁵
10⁻⁸          32       2.1·10¹              7           5.4·10⁻⁶
10⁻⁸          64       1.8                 12           4.4·10⁻⁷
10⁻⁶          16       3.0·10⁶             6           3.5·10⁻³
10⁻⁶          32       2.2                11           5.4·10⁻⁴
10⁻⁶          64       1.7                20           4.4·10⁻⁵

The rank decision gets steadily worse. At 10⁻⁶ noise the gap collapses from three million at sixteen samples to 1.7 at sixty-four, and a fixed threshold reads twenty where the degree is six. The reason is countable: every extra sample adds another noise-driven direction to the Loewner matrix, and the largest of m noise singular values grows faster than the six signal ones do, so the cliff fills in from below. Oversampling is not a lever on the rank decision; it is the thing that destroys it.

The model at a known order gets steadily better. The last column falls by a factor of eighty for four times the samples, at both noise levels. So the sentence earlier in this essay — no order does better than 1.6·10⁻³ at noise 10⁻⁶ — is a statement about twenty-four samples and not about the noise. Averaging works. What it cannot do is work through the rank decision, because the rank decision is a maximum over directions rather than a mean over samples, and a maximum does not average.

The procedure that falls out of the two columns is not the one either half suggests alone: read the order from a small sample set, where the cliff is still there, and then fit at that order with every sample the budget allows. Sixteen samples give the degree correctly at every noise level measured; sixty-four give a model eighty times better at that degree and an order nobody should believe. Doing both is one extra Loewner pencil on sixteen points, which costs nothing beside the sampling.

There is a general shape in that split and it is worth naming, because this collection keeps arriving at it from different directions. A quantity that is a maximum does not average and a quantity that is a mean does. The model’s error is an approximation error — a least-squares-like object that improves as the data accumulates, because the noise enters it through inner products that cancel. The rank decision is a threshold applied to the largest of a set of singular values, and adding data adds candidates to that set, so the largest can only rise. Every extra sample is evidence for the model and a fresh chance for the decision to go wrong.

The same asymmetry is why an estimate that can be fooled is a lower bound with a probability rather than a bound: it too is a maximum over sampled directions, and more sampling raises it without ever making it certain. What is different here is that both quantities come out of the same computation, so a caller reading only one of them gets a consistent and wrong picture in either direction — a degree of twenty with a model that is excellent at six, or a model at six with no reason to have chosen six.

Where this leaves the field #

The reduction field’s arc has been about what an approximation costs when it is committed before the arithmetic. Balanced truncation commits it with a bound known in advance. Interpolatory projection commits it at points somebody named, and can return a model that cannot be run. This one commits it at points somebody measured, and its error is bounded below by the measurement’s own.

That last is the honest summary of a data-driven method and it is not a criticism. A model built from measurements is as good as the measurements, exactly, and the construction adds nothing to the error except at the rank decision. Every other method in this field requires a state matrix that a measurement cannot supply, so the comparison is not between this and a better method — it is between this and nothing.

What the whole phase has in common #

This essay closes a run of three that were written to clear gaps recorded rather than to open new ground, and the three turn out to share a shape worth naming.

The support points chosen by a residual build an approximant from values of a function, and degrade past the noise in those values. A restricted backward error asks how far an answer is from exact under an assumption about which data is trustworthy. And this one builds a whole model out of measurements and inherits their error as a floor.

All three are about where the uncertainty in a computation actually lives, and in all three it is not in the arithmetic. The linear algebra in each is backward stable and well conditioned; the thing that decides the answer is a property of the data or a decision about it. That is the same conclusion the graph field’s ranking essay reaches by a completely different route, and it is unusual enough on this site to be worth stating twice.

It does not weaken the site’s premise; it sharpens it. The identity forward ⪅ condition × backward is a statement about a computation performed on data that is taken as given, and every essay here that finds the arithmetic innocent has found a term the identity does not carry — the machine field found one, the approximation-before-linearisation essays found another, and these three found a third. The useful discipline is the same in all of them: measure the computation’s contribution, state it, and then measure what else there is.

What is not measured here #

Two things this construction is usually asked about are absent from the measurements above, and both absences are deliberate.

There is no error bound. Balanced truncation comes with one — twice the sum of the discarded Hankel singular values, known before the reduction is performed. The Loewner construction has nothing comparable: the singular values of the pencil say what the degree is, and they do not bound the error of a model built below it. The last column of the first table is a measurement rather than a certificate, and on a system whose degree is genuinely infinite there would be no clean order to stop at.

And the model is not balanced. Its state has no interpretation, its Gramians are whatever the projection produced, and quantities like the Hankel singular values are not available from it without solving two Lyapunov equations — which requires the state matrix the construction was invented to avoid needing.

Both absences are the price of the input being weaker. A method that is handed only values cannot report a quantity defined in terms of a state-space realisation, and a bound that is derived from Gramians is not available to a method that never forms one. What is available is the interpolation, which is exact, and the degree, which is a rank — and this essay’s whole business is measuring what happens to the second when the data is not clean.