A pole between the samples

An error estimate made of samples took the number every Loewner code reads to choose its order, the pencil’s next singular value $\sigma_{k+1}/\sigma_1$, and found it the error of the reduced model for two orders and then optimistic by a decade an order, up to 6,680 at the sixth. What worked instead was sampling the transfer function again, at the 47 geometric midpoints between the 48 samples the model was built from: the true error was between 1.00 and 2.09 times the largest error there, at every order on every system. The essay ended by asking for something cheaper. AAA chooses its next support point where its own residual is largest, and “a Loewner model could choose its next check point the same way, where the current model and a model one order lower disagree most. The prediction is that on these four systems four such points, chosen adaptively, estimate the error within a factor of two at every order — as well as the 47 midpoints, and better than the four nearest the poles chosen by hand.”

That prediction is tested here, and with exact samples it is right by a wider margin than it asked for. Put noise on the samples and it still holds on every model where it can hold at all, and that turns out to be the smaller question. The larger one is that past a certain order a model built from noisy samples has a pole sitting between two of them, its error at that point is infinite, and no check made of a finite set of points can say so. The test that does see it was never an estimate. It is the model’s own poles, which every Loewner code already has in hand.

Four points, chosen by the next order #

The systems and samples are the earlier essay’s: three heat-equation systems showing twelve or eight modes, and a modal system with eight poles, each sampled at 48 real points spaced geometrically left of its spectrum, with the order-$k$ model the Loewner pencil projected onto its leading $k$ singular directions. The true error is the largest difference from the transfer function over 400 test points across the interval, over the largest value of $|H|$ there.

The check the prediction proposes needs one change to be cheap. A model one order lower is the comparison the prediction named, but the order-$(k+1)$ model comes out of the same pencil at the same cost, and it is the better comparison, for a reason that will matter below: the next order is more accurate than this one, while the previous order is less. So the measurement takes $m_k$ and $m_{k+1}$ from the same 48 samples, evaluates $|m_k - m_{k+1}|$ over the 400 test points, finds its local maxima, and keeps the four largest. The transfer function is sampled at those four points and nowhere else, and the largest error of $m_k$ there is the estimate.

That gives two estimates, not one. The first costs four new samples and is the error at the four chosen points. The second costs nothing at all: it is the disagreement itself, $\max |m_k - m_{k+1}|$ over the test points, read as the error of $m_k$ without any new value of the function. Both are compared against the three estimates the earlier essay already had — the next singular value, the 47 midpoints and the four midpoints nearest the poles.

The four chosen points give the error. Over the twenty orders above the arithmetic floor on the four systems, the true error over the error at the four chosen points is between 1.00 and 1.02 — not within the factor of two the prediction asked for but within two per cent, at every order, on every system. The 47 midpoints, which cost forty-three more samples, are within 2.09; the four midpoints chosen by hand near the poles within 2.72; the singular value is optimistic by up to 6,680. And the free estimate, the disagreement with no sample at all, is the error within three per cent at every order but one: between 0.973 and 1.01, with a single 2.01 at the twelve-mode system’s sixth order.

Why the next order points at the error #

The two estimates work for the same reason, and the reason says exactly where they will stop working. Write the error of the order-$k$ model as $m_k - H$, and the disagreement as $m_k - m_{k+1}$. They differ by $m_{k+1} - H$, the error of the next order. On these systems each order is about a hundred times more accurate than the one before — the eight-pole system goes $1.6 \cdot 10^{-2}$, $4.5 \cdot 10^{-4}$, $3.1 \cdot 10^{-6}$, $3.1 \cdot 10^{-8}$, $8.5 \cdot 10^{-10}$, $7.3 \cdot 10^{-12}$ — so the disagreement is the error of $m_k$ to within one per cent of itself, everywhere along the interval at once. It is not an estimate of the size of the error that happens to come out right. It is the error, as a function, with a one per cent perturbation on it.

At order three the two curves cannot be told apart at the scale of the page. The error peaks at $3 \cdot 10^{-6}$ at the end nearest the poles, dips through a few sign changes, and runs across a broad hump near $10^{-6}$ over the rest of the interval; the dashed curve lies on it all the way. Turn the dial to five and the whole picture sits three to five decades lower, with more sign changes near the poles, and the curves still coincide. So the four peaks of the disagreement are the four peaks of the error. There is nothing clever about where they sit; they sit where the error is largest, because the disagreement says where that is.

It also says where the error is largest when that is not where anyone would look. The earlier essay found the worst error in the first gap on eighteen of the twenty orders, and named the two exceptions: order three of the falling-residue system, whose worst point is at $s = -33{,}936$, and order two of the eight-mode system, at $s = -6{,}958$. The four midpoints nearest the poles miss both, because both are deep inside the interval. The four chosen points find both. On the eight-mode system the first chosen point is the worst point exactly; on the falling-residue system it is at $s = -33{,}606$, inside the same hump, where the error is within half a per cent of its peak.

The one order where the free estimate drifts is the one the arithmetic predicts. At the twelve-mode system’s sixth order the error is $4.3 \cdot 10^{-13}$ and the seventh order’s is $2.2 \cdot 10^{-13}$ — both are at the rounding, so the next order is twice as good rather than a hundred times, $m_7 - H$ is not small beside $m_6 - H$, and the disagreement is half the error rather than all of it. The four chosen points still land on the right humps and give the error within two per cent there, because they then measure the true function rather than trusting the next order. The free estimate trusts that the next order is better; the one that costs four samples checks.

Noise stops the next order from being better #

A model with no matrices behind it found that noise of $10^{-10}$ on the samples sets a floor under every quantity the construction produces: the singular values stop falling at the noise rather than at the rounding, and the gap that reads the degree off the pencil shrinks from $2 \cdot 10^{8}$ to 4.4. The earlier essay’s estimate was measured on exact samples only and said so. So the same comparison is run again on two of the systems — the twelve-mode one and the eight-pole one — with relative noise of $10^{-8}$ and $10^{-6}$ on every sample, at orders one to eight, thirty-two models in all.

The exact models fall by about two decades an order until the arithmetic stops them. The noisy ones follow them down and then stop at roughly the noise: the twelve-mode system with noise $10^{-8}$ reaches $9.8 \cdot 10^{-9}$ at order four and then reads $9.5 \cdot 10^{-8}$, $9.1 \cdot 10^{-8}$, $7.8 \cdot 10^{-8}$ and $5.2 \cdot 10^{-8}$; with noise $10^{-6}$ it reaches $8.7 \cdot 10^{-7}$ at order three and then jumps to $5.0 \cdot 10^{-4}$ at four. Past that point the next order is not a hundred times better. Sometimes it is worse. And that removes the one fact the free estimate stood on.

It also does something the first look at those dashed lines does not show, and the rings mark it. Every one of the seventeen ringed models has a real pole inside the interval the samples cover — between one and five of them, forty-three in all. None of the twenty exact models has one, at any order, on any system. With noise, every one of the four noisy cases grows its first such pole at exactly the order where the error stops falling, or one order after it, and keeps at least one at every order after that.

A pole between two samples #

A pole on the sampled interval means the model’s error there is not small or large. It is infinite: the transfer function is smooth at that point and the model goes through infinity. The 400 test points cannot see that, because the chance that one lands on the pole is zero, and what they report instead is how close the nearest one happened to fall.

The worst case among the thirty-two is the one in the figure. The order-4 model from samples with noise $10^{-6}$ has a pole at $s = -3060.9$, under one per cent from the sample at $-3032.3$, in a gap that runs to $-3295.5$. The 400 test points read $5.0 \cdot 10^{-4}$, from the one nearest the pole. The 47 midpoints read $2.0 \cdot 10^{-6}$, 252 times less, which is about what the same model reads everywhere else on the interval and about what a model at the noise level ought to read. A code checking that model at its midpoints would accept it as accurate to two parts in a million. It is accurate to a few parts in a million everywhere outside a neighbourhood about a seventh of a gap wide; inside it the error passes one part in a thousand within seven ten-thousandths of the pole’s position, one per cent within seven millionths, and then grows without limit.

The pole is the familiar object, and it has a name. The points the algorithm chose met it in AAA under noise: a Froissart doublet, a pole and a zero of the approximant almost on top of each other, with a residue at about the noise level, contributing nothing to the function and appearing in the pole list. A Loewner model built past the order its samples can support makes the same thing for the same reason. The projection has more freedom than the data has information, and it spends the spare freedom fitting the noise, which a rational function does most cheaply with a pole and a zero that nearly cancel.

How narrow each one is can be read off its residue. Near a pole $p$ the error is the residue over the distance, so with the residue divided by the function’s scale and by $|p|$ — a number $\rho$ with no units — the error passes a level $\varepsilon$ only where $|s/p - 1| < \rho/\varepsilon$. On the 43 poles measured, $\rho$ is at most $6.6 \cdot 10^{-8}$; at noise $10^{-8}$ it is between $1.2 \cdot 10^{-11}$ and $4 \cdot 10^{-10}$. Consecutive samples are 9 to 10 per cent apart. So the window in which a doublet pushes the error past ten times the noise is at most 1.3 per cent wide, about a seventh of a gap, and whether a check point lands in it is luck: a midpoint sits about four per cent from the samples on either side, and the pole can be anywhere in the gap. Over the forty-three, the closest any midpoint comes to a pole is 0.18 per cent, and even that one reads a finite number. On the model in the figure the arithmetic is exact: $\rho = 6.6 \cdot 10^{-8}$, the midpoint of that gap is 3.3 per cent from the pole, and $6.6 \cdot 10^{-8}/0.033$ is the $2.0 \cdot 10^{-6}$ the midpoints report.

Every finite check is fooled #

No finite number estimates an infinite error, so the only question for a check on such a model is whether what it reports is large enough to make a code reject it. On the seventeen it is not. The largest value any check reads on any of them — four chosen points or 47 midpoints — is $8.7 \cdot 10^{-4}$. One ten-millionth of the pole’s distance from the origin away from it, the same models’ errors are between $1.0 \cdot 10^{-3}$ and 0.66, and between 133 and 48,600 times what the checks report. A step that small is still a point on the line a code may be asked to evaluate at.

The four chosen points do one thing better: they are pulled onto the poles. On every one of the seventeen models, one of the four lies within half a per cent of a pole: the disagreement $|m_k - m_{k+1}|$ peaks there because one of the two models has the pole and the other either does not or has it somewhere else. So the chosen points know where to look. They do not report what is there, because a point half a per cent from a doublet with $\rho = 10^{-10}$ reads an error of $2 \cdot 10^{-8}$, which is the noise. On the model in the figure the first chosen point is $s = -3060.5$, four tenths of a unit from the pole, and it reads $5.0 \cdot 10^{-4}$ — 252 times what the midpoints say and the clearest warning any check gives, and still a finite number about a model whose error is not.

A check is a measurement at points, and a measurement at points bounds nothing between them. The earlier essay said as much of its midpoints — “a function with a feature narrower than the spacing between a sample and its midpoint could hide from it” — and it had in mind a feature of the transfer function. The feature that hides is the model’s, and the model makes it out of the noise.

Where the checks still hold #

Leave the seventeen out and look at the fifteen noisy models with no pole between the samples. These are the models whose error is a number, and the question of how well a check estimates it is the same question the exact samples asked.

The checks still work. The four chosen points give the error within 1.00 to 1.07 on all fifteen, and the 47 midpoints within 1.00 to 2.09 — the same bands as on exact samples, to two decimal places. Noise at these levels changes what the model is; it does not change how well a model with no pole in the way can be checked. What it does change is the free estimate, and in the direction that is safe. On ten of the fifteen the disagreement is still within four per cent of the error and on one it overstates by 1.30; on the other four it overstates by 3.4, 9.7, 14.2 and 579. The 579 is the twelve-mode system’s order-3 model with noise $10^{-6}$, whose error is $8.7 \cdot 10^{-7}$. Its disagreement is with the order-4 model, and the order-4 model is the one in the figure above: its pole is what the disagreement measures, at $5.0 \cdot 10^{-4}$ on the test points. A pessimistic estimate is a cautious one, and a code reading it would stop at the order before; but it reads the next order’s failure as this order’s, and the factor it is wrong by has no bound.

So the free estimate is good exactly while each order is a hundred times better than the last, which is the condition the noise ends, and the four samples are good wherever the model has no pole in the way. Neither says anything about a model that has one.

The model’s own poles are the check #

What says it is free. A Loewner model is a descriptor system $(E, A, B, C)$ of order $k$, its poles are the $k$ generalised eigenvalues of $sE - A$, and any code that built the model has $E$ and $A$ in hand. Computing $k$ eigenvalues of a pencil of order eight costs less than evaluating the model at the 400 test points used here to measure it. A real pole inside the interval the samples cover is not a matter of degree. It is a defect the model either has or does not, and its presence is visible without another sample of the function. Here the test was made by finding every sign change of $\det(sE - A)$ along the interval and bisecting each to the rounding, which finds every real pole there and nothing else; a code would take the eigenvalues and keep the real ones between the interval’s ends. On the twenty exact models above the arithmetic floor it finds none. On the thirty-two noisy ones it finds every one of the forty-three.

That makes a stopping rule. Raise the order while the model has no pole inside the sampled interval; at each order sample the function at the four points where the next order disagrees most; stop when the error there passes the tolerance, or at the last order before a pole appears. Where to put the poles of a rational function found that a rational approximant’s accuracy is bought by where its poles sit, outside the region it must be good on; a pole inside that region buys nothing and costs everything. On the four noisy cases that last order is four, four, three and four, and its error on the test points is $9.8 \cdot 10^{-9}$, $5.4 \cdot 10^{-9}$, $8.7 \cdot 10^{-7}$ and $7.1 \cdot 10^{-6}$. The best any order reaches is the same on the first three and $2.9 \cdot 10^{-6}$ on the last, the eight-pole system with noise $10^{-6}$, whose order four is already a little worse than its order three without yet having a pole. So the rule gives up at most a factor of 2.45, and on three of four cases nothing. It never returns a model with infinite error. Reading $\sigma_{k+1}/\sigma_1$ against the tolerance cannot say that much: with noise the singular values stop falling at the noise, which leaves a code choosing between a model that is too low and one with a pole in it.

Interpolating at the model’s own poles trusted a reduced model’s poles enough to choose the next interpolation points from them; here they are a test the model can fail. The eigenvalues that are answers to nothing met the same kind of object in a linearisation: exactly right about the approximation, corresponding to nothing in the problem. The arithmetic marks none of them; only knowing where the function is analytic — here, the whole sampled interval — calls them spurious.

What the chosen points are worth now #

Four samples instead of forty-seven is the trade the prediction proposed, and on every model where a check means anything it is a good one: four points chosen by the next order give the error within seven per cent where the 47 midpoints give it within a factor of 2.1, and they find the worst point when it is not where a hand-chosen four would look. Exact at the points that were named made the general case against trusting a model at the points it was built from; the chosen points are the opposite habit, points the model did not use, placed where it is least sure of itself.

But the order of the two tests matters. A code that checks the error first will, on the model in the figure, read $2.0 \cdot 10^{-6}$ at the midpoints or $5.0 \cdot 10^{-4}$ at the chosen points, and may accept either. A code that checks the poles first never gets that far: the pole test is cheaper, exact rather than estimated, and disqualifies the model before an estimate can mislead. A small residual is not a small error is the general warning this is one case of. Every check above is a residual of the model against the function at some points, and a model can be wrong between its residuals without bound.

The bound that is known in advance is what balanced truncation offers instead, and it does not have this problem: a truncated balanced realisation of a stable system keeps its poles in the left half-plane by construction, and its error is bounded on the whole imaginary axis before it is built. A Loewner model is a fit, and a fit to noisy data has no such guarantee. The poles at least say where it has gone wrong.

What thirty-two noisy models do not show #

Two systems, two noise levels, one kind of noise: independent Gaussian perturbations, relative to each sample’s magnitude, with one seed. Noise that is correlated from sample to sample — a smooth miscalibration rather than a scatter — might be fitted by the projection without any doublet, and a model that fits a smooth error looks right to every test made here, its own poles included. Real samples only, on an interval left of the spectrum where the transfer function is real and smooth. Sampled on the imaginary axis the model is complex, a spurious pole need not sit on the sampled segment to spoil the model near it, and the test becomes a distance from the segment rather than a sign change. A pole a short way off the axis gives a finite error, which the checks would see, and the boundary between a pole the checks see and one only the eigenvalues see is not measured here. Eight orders, which reach five poles on the worst case; at higher orders with more samples, a model could plausibly carry a doublet outside the interval that still lies close enough to spoil it, and the rule above does not look there.

Still open: the imaginary axis, and what a doublet says about the noise #

The imaginary axis. The prediction carried over from the earlier essay still stands untested, because it needs the complex form of the construction: sampled at frequencies $s = i\omega$, a lightly damped system has resonant peaks between the samples. With the eight-pole system’s poles moved off the real axis to damping ratios of 0.05, the prediction with a sign is that the 47 midpoints understate the true error by more than a factor of ten at some order, that two points per gap bring it back within three, and — this essay’s addition — that four points chosen by the next order’s disagreement give the error within a factor of 1.1 at every order where the model has no pole within one sample spacing of the sampled segment.

The noise, read from the doublets. The largest residue among a model’s spurious poles, in the units used above, is between 0.010 and 0.040 times the noise at $10^{-8}$ and between 0.005 and 0.066 times it at $10^{-6}$. If that ratio is a property of the construction rather than of these two levels, the doublets measure the noise without replicate samples. The prediction is that at relative noise from $10^{-10}$ to $10^{-4}$, on both systems, the largest doublet residue over the noise level stays between 0.003 and 0.1 at every order that has one — so that the first doublet to appear estimates the noise within a factor of thirty, and the order it appears at estimates the order the samples can support.