Numerical Linear Algebra
where the algebra stops being the arithmetic
Reduction, and what a model is for

An error estimate made of samples

A Loewner model built from samples of a transfer function comes with no error bound, and the number every code reads to choose its order — the pencil's next singular value — is the obvious stand-in. On four systems sampled at 48 points it is the error of the order-k model for the first two orders, within a factor of 2.5, and then falls away from it by about a decade an order: by the fifth it is optimistic by 109, 181 and 457, and on one system by 6,700 at the sixth. An AAA fit of the same degree to the same samples has an error within a factor of six of the Loewner model's at every order, so the shortfall is not the projection's. The model's misfit at its own samples is optimistic the same way. What works is more samples: the error at the 47 midpoints between the samples is within a factor of 2.2 of the true error at every order on every system, and the four midpoints nearest the poles alone within a factor of three.
19 min read 5 figures Measured, not assumedExact ground truthClaims that can fail

Worth reading first: A model with no matrices behind it · A model that is a rational function.

A model with no matrices behind it built a reduced model out of nothing but values of a transfer function: two matrices of divided differences, the Loewner pencil, projected onto their leading singular subspaces. At the McMillan degree the model interpolated every sample and matched the function it never saw to 10−910^{-9}. The essay ended on what the construction does not have. “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.”

That sentence is careful, and practice is not. Every code that builds a Loewner model below its degree has to choose the order, and the number it reads to choose is the pencil’s singular values: keep kk when σk+1/σ1\sigma_{k+1}/\sigma_1 has fallen below the tolerance. That reads the singular value as an estimate of the error, which the first essay declined to claim. This essay measures how good an estimate it is, finds out where it fails, and finds one that does not.

Four systems, sampled where nothing is singular

The systems are two kinds. Three come from the earlier essays’ heat equation on forty states, driven and observed so that only twelve or eight of its modes appear in the transfer function: one with twelve modes of similar weight, one with twelve whose residues fall by a factor of three a mode, and one with eight. The fourth is a modal system with eight poles between −0.5-0.5 and −9-9 and residues of both signs. Each is sampled at 48 real points spaced geometrically along an interval to the left of its whole spectrum — [−105,−2000][-10^5, -2000] for the twelve-mode systems, [−105,−1000][-10^5, -1000] for the eight-mode one, [−1000,−12][-1000, -12] for the eight-pole one — so that the function is smooth everywhere it is sampled and everywhere it is checked.

The true error of a model is read at 400 points spread over the same interval, as the largest difference from the transfer function divided by the largest value of ∣H∣|H| there. Orders one to six are built on every system. Below about 10−1310^{-13} the error is the arithmetic’s rather than the order’s, and those orders are marked and left out of every comparison: on the falling-residue system that is from the fourth order on.

The singular value is the error, for two orders

The figure at the top of the page is the eight-pole system. At order one the true error is 1.60⋅10−21.60 \cdot 10^{-2} and σ2/σ1\sigma_2/\sigma_1 is 2.19⋅10−22.19 \cdot 10^{-2}; at order two, 4.52⋅10−44.52 \cdot 10^{-4} against 3.71⋅10−43.71 \cdot 10^{-4}. A code that read the singular values would choose the right order at these tolerances. At order three the error is 3.07⋅10−63.07 \cdot 10^{-6} and the singular value 7.40⋅10−77.40 \cdot 10^{-7}, four times too small. At order five, 8.53⋅10−108.53 \cdot 10^{-10} against 1.87⋅10−121.87 \cdot 10^{-12}: the singular value says the model is accurate to twelve digits and it is accurate to nine. At order six it says fifteen and the model has eleven.

The true error of the Loewner model of each order over the pencil's next singular value, on four systemstwelve modes: order 1 1.28, order 2 1.18, order 3 1.81, order 4 12.1, order 5 109, order 6 937. twelve, falling: order 1 1.53, order 2 1.50, order 3 1.31, order 4 1.54 (floor), order 5 3.68 (floor), order 6 48.2 (floor). eight modes: order 1 1.29, order 2 1.12, order 3 2.55, order 4 18.6, order 5 181, order 6 46.4 (floor). eight poles: order 1 0.733, order 2 1.22, order 3 4.15, order 4 30.1, order 5 457, order 6 6.68e+3.123456110¹10²10³10⁴order of the modeltrue error ÷ σₖ₊₁/σ₁twelve modestwelve, fallingeight modeseight polesopen: at the arithmetic floorabout a decade an order
Fig. 1 The true error of the Loewner model of each order over the pencil’s next singular value σk+1/σ1\sigma_{k+1}/\sigma_1, on four systems. Open circles are orders at the arithmetic floor.

The other systems do the same. Over the first two orders the error over σk+1/σ1\sigma_{k+1}/\sigma_1 is between 0.73 and 1.53 on all four. Then it climbs: on the twelve-mode system 1.81, 12.1, 109 and 937 at orders three to six; on the eight-mode system 2.55, 18.6 and 181 at orders three to five; on the eight-pole system 4.15, 30.1, 457 and 6,680. About a decade per order once it starts.

The climb is the difference of two rates. From the third order on, the eight-pole system’s singular values fall by 2.7 to 3.2 decades an order — σ4/σ3\sigma_4/\sigma_3 is 1.4⋅10−31.4 \cdot 10^{-3}, σ6/σ5\sigma_6/\sigma_5 is 1.8⋅10−31.8 \cdot 10^{-3} — while the true error falls by 1.6 to 2.2: from 3.07⋅10−63.07 \cdot 10^{-6} to 3.12⋅10−83.12 \cdot 10^{-8} to 8.53⋅10−108.53 \cdot 10^{-10} to 7.25⋅10−127.25 \cdot 10^{-12}. Each order the singular value gains about a decade on the error, and there is no constant that converts one into the other. That is the difference from the bound balanced truncation offers, which the bound that is known in advance found attained exactly on ordinary problems and loose by at most a factor of about two where the removals interact; a bound with a bounded slack can be trusted, and an estimate whose slack grows a decade an order cannot. The falling-residue system reaches the arithmetic floor at order four, still within a factor of 1.5, before the climb begins. The singular value is a good estimate exactly as long as nobody needs one, at the first orders where the error is large, and it is optimistic by orders of magnitude at the orders where a tolerance would stop.

A code that chose its order by σk+1/σ1<10−11\sigma_{k+1}/\sigma_1 < 10^{-11} on the eight-pole system would stop at five and deliver 8.5⋅10−108.5 \cdot 10^{-10}, eighty-five times its tolerance, with every number it computed saying it had succeeded.

Two fits, one error

There are two places the shortfall could come from. The singular values might be misreporting how well the samples can be fitted, or the Loewner model might be fitting them worse than the singular values say is possible, because projecting onto singular subspaces is one way of building a model and not necessarily the best. A second, independent fit separates them.

For every order above the arithmetic floor on four systems: the true error of the Loewner model against that of an AAA rational fit of the same degree to the same samplestwelve modes, order 1: Loewner 0.0122, AAA 0.00911; twelve modes, order 2: Loewner 5.16·10⁻⁵, AAA 7.66·10⁻⁵; twelve modes, order 3: Loewner 2.46·10⁻⁷, AAA 6.03·10⁻⁷; twelve modes, order 4: Loewner 3.2·10⁻⁹, AAA 1.04·10⁻⁸; twelve modes, order 5: Loewner 3.76·10⁻¹¹, AAA 5.28·10⁻¹¹; twelve modes, order 6: Loewner 4.25·10⁻¹³, AAA 5.99·10⁻¹³; twelve, falling, order 1: Loewner 3.72·10⁻⁶, AAA 2.77·10⁻⁶; twelve, falling, order 2: Loewner 3.1·10⁻¹⁰, AAA 3.41·10⁻¹⁰; twelve, falling, order 3: Loewner 2.56·10⁻¹³, AAA 3.53·10⁻¹³; eight modes, order 1: Loewner 0.00961, AAA 0.00749; eight modes, order 2: Loewner 3·10⁻⁵, AAA 5.03·10⁻⁵; eight modes, order 3: Loewner 1.17·10⁻⁷, AAA 2.49·10⁻⁷; eight modes, order 4: Loewner 8.66·10⁻¹⁰, AAA 2.52·10⁻⁹; eight modes, order 5: Loewner 3.52·10⁻¹², AAA 4.37·10⁻¹²; eight poles, order 1: Loewner 0.016, AAA 0.0116; eight poles, order 2: Loewner 4.52·10⁻⁴, AAA 0.0021; eight poles, order 3: Loewner 3.07·10⁻⁶, AAA 5.02·10⁻⁶; eight poles, order 4: Loewner 3.12·10⁻⁸, AAA 6.19·10⁻⁸; eight poles, order 5: Loewner 8.53·10⁻¹⁰, AAA 1.36·10⁻⁹; eight poles, order 6: Loewner 7.25·10⁻¹², AAA 1.09·10⁻¹¹. The dashed lines are a factor of six either side of equality.10⁻¹⁴10⁻¹²10⁻¹⁰10⁻⁸10⁻⁶10⁻⁴10⁻²110⁻¹⁴10⁻¹²10⁻¹⁰10⁻⁸10⁻⁶10⁻⁴10⁻²1Loewner model's errorAAA fit's errortwelve modestwelve, fallingeight modeseight polesdashed: a factor of six either waytwo fits, one error
Fig. 2 For every order above the arithmetic floor on the four systems: the true error of the Loewner model against that of an AAA rational fit of the same degree to the same 48 samples. The dashed lines are a factor of six either side of equality.

The points the algorithm chose introduced AAA: a rational approximant in barycentric form whose support points are chosen one at a time by its own residual and whose weights come from the null vector of a Loewner matrix over the other samples. It shares the samples and nothing else with the projected pencil. At every order above the floor on every system, its error is within a factor of six of the Loewner model’s, and within two on thirteen of the twenty cells: 1.04⋅10−81.04 \cdot 10^{-8} against 3.20⋅10−93.20 \cdot 10^{-9} on the twelve-mode system at order four, 1.36⋅10−91.36 \cdot 10^{-9} against 8.53⋅10−108.53 \cdot 10^{-10} on the eight-pole system at order five. Two methods that build their models differently reach the same error, and both are far above the singular value. The Loewner model is fitting about as well as these samples allow a rational function of that degree to be fitted; it is the singular value that is optimistic about what is possible.

The singular values of the Loewner matrix built from 48 samples of the twelve-mode system, of the square root at the same points, and of random numbers at the same pointstwelve-mode samples: 1, 0.0095, 4.4·10⁻⁵, 1.4·10⁻⁷, 2.6·10⁻¹⁰, 3.5·10⁻¹³, 4.5·10⁻¹⁶, 2.7·10⁻¹⁶. √(−s) at the same points: 1, 0.059, 0.0034, 1.9·10⁻⁴, 10⁻⁵, 5.3·10⁻⁷, 2.7·10⁻⁸, 1.4·10⁻⁹. random numbers: 1, 0.89, 0.74, 0.71, 0.65, 0.61, 0.45, 0.36. Each over its largest.123456789101112131410⁻¹⁷10⁻¹⁵10⁻¹³10⁻¹¹10⁻⁹10⁻⁷10⁻⁵10⁻³10⁻¹indexσᵢ ÷ σ₁twelve-mode samples√(−s) at the same pointsrandom numbersrandom data does not decaythe decay is information, not structure
Fig. 3 The singular values of the Loewner matrix built from the 48 samples of the twelve-mode system, of −s\sqrt{-s} at the same points, and of random numbers at the same points, each over its largest.

The singular values are not empty, though. Built from random numbers at the same 48 points, the Loewner matrix’s singular values barely fall: 0.89, 0.74, 0.71 and still 0.36 by the eighth. Built from −s\sqrt{-s}, which has no finite degree, they fall steadily, a little over a decade an index. Built from the twelve-mode system, they fall off a cliff to the rounding by the seventh. So their decay is information about the function and not an artefact of the divided differences; it is the wrong size, not the wrong shape. A singular value measures how far LL is from a matrix of lower rank in the Frobenius-type sense of its whole entries, while the model’s error is the worst point of a function between samples, and the two are not the same norm on the same object.

The samples flatter the model too

The quantity a code can compute without anything but its data is the model’s misfit at its own samples. Below the degree the model does not interpolate, so the misfit is not zero, and a reasonable guess is that a model that misses its samples by 10−1210^{-12} misses the function between them by about as much.

It does not. On the eight-pole system at order five the misfit at the samples is 4.7⋅10−124.7 \cdot 10^{-12} and the error between them 8.5⋅10−108.5 \cdot 10^{-10}, 182 times larger; at order six, 751 times. On the twelve-mode and eight-mode systems the ratio at order five is 40.7 and 64.2. The misfit follows the singular value down rather than the error, which is no surprise once it is said: the model is built from the samples, by a projection chosen to fit them, and a fit is always best where it was fitted. It is the residual of the model’s own construction, and a small residual is not a small error is the site’s oldest warning about reading one as the other.

The order-5 Loewner model of the eight-pole system: its error along the sampled interval, with its misfit at the 48 samples and its error at the 47 midpointsLargest error 8.53·10⁻¹⁰ at s = -12.201, next to the end of the interval nearest the poles at −12; largest misfit at the samples 4.68·10⁻¹²; largest midpoint error 5.15·10⁻¹⁰. The axis is −s, logarithmic.order 5largest error8.5·10⁻¹⁰largest misfit at samples4.7·10⁻¹²largest at midpoints5.1·10⁻¹⁰10²10³10⁻¹⁶10⁻¹⁴10⁻¹²10⁻¹⁰−s, distance left along the real axiserror, over the largest |H|errormidpointssamplesthe poles are to the left of the axis's startworst between the first samples
Fig. 4 The order-5 Loewner model of the eight-pole system: its error along the sampled interval, its misfit at the 48 samples (open) and its error at the 47 midpoints between consecutive samples (filled). The horizontal axis is −s, on a logarithmic scale; the poles lie to the left of where it starts.

The error along the interval shows where. Over most of it the model is accurate to between 10−1110^{-11} and 10−1210^{-12}, dipping to zero at points where its error changes sign. At the end nearest the poles, between −s=12-s = 12 and about 20, the error rises to 8.5⋅10−108.5 \cdot 10^{-10}, with its worst point at s=−12.2s = -12.2, between the first two samples. The samples there are fitted to 10−1210^{-12} or so; the function between them, which bends hardest because it is nearest the poles, is not. The largest error is never at a spurious pole: the reduced models’ poles were computed at every order on two of the systems, and none lies inside the sampled interval. It is simply where the function is hardest and the samples, spaced evenly in the logarithm, are no closer together than anywhere else. The same holds across the whole measurement: of the twenty orders above the floor on the four systems, eighteen have their worst error in the first gap, between the sample nearest the poles and the next one. The two exceptions are order three of the falling-residue system and order two of the eight-mode one, whose worst points sit well inside the interval, at s=−33,936s = -33{,}936 and s=−6,958s = -6{,}958; the midpoint estimate holds on both, at ratios of 1.00, because a midpoint sits in every gap and not only in the first.

Midpoints are an estimate

That picture suggests the estimate. If the error lives between the samples, check it between the samples: evaluate the transfer function at the geometric midpoint of each consecutive pair, 47 more values, and compare the model there.

The true error of the Loewner model of each order over its largest error at the 47 sample midpoints, filled, and at the four midpoints nearest the poles, open, on four systemstwelve modes: order 1 1.71 and 1.71, order 2 1.11 and 1.57, order 3 1.02 and 1.02, order 4 1.25 and 1.25, order 5 1.47 and 1.47, order 6 1.66 and 1.66. twelve, falling: order 1 1.35 and 1.35, order 2 1.38 and 2.72, order 3 1.00 and 1.68. eight modes: order 1 1.64 and 1.85, order 2 1.00 and 1.48, order 3 1.05 and 1.05, order 4 1.30 and 1.30, order 5 1.50 and 1.50. eight poles: order 1 2.09 and 2.09, order 2 1.18 and 1.18, order 3 1.12 and 1.12, order 4 1.42 and 1.42, order 5 1.66 and 1.66, order 6 1.90 and 1.90. Orders at the arithmetic floor are left out.12345600.511.522.53order of the modeltrue error ÷ midpoint estimatetwelve modestwelve, fallingeight modeseight polesfilled: all 47 midpoints · open: the 4 nearest the polesa certificate made of samples
Fig. 5 The true error of the Loewner model of each order over its largest error at the 47 sample midpoints, filled, and at only the four midpoints nearest the poles, open, on four systems. Orders at the arithmetic floor are left out.

At every order above the floor on every system, the true error is between 1.00 and 2.09 times the largest midpoint error. The estimate never overstates the error and never understates it by more than a factor of 2.1, across a range of errors from 10−210^{-2} to 10−1310^{-13}, while the singular value’s understatement ran to 6,680. And most of its information is in four of the 47 points. Using only the four midpoints nearest the poles, the true error is within a factor of 2.72 at every order above the floor, and at fifteen of the twenty cells it is the same as using all 47, because the worst midpoint is one of those four.

The estimate costs four to forty-seven extra samples, and nothing else: no matrices, no Gramians, nothing the construction was invented to avoid. That is the honest replacement for the bound the first essay said the method lacks. It is not a bound — a function with a feature narrower than the spacing between a sample and its midpoint could hide from it — but it is a measurement of the right quantity, the error between the samples, where the singular value and the sample misfit are measurements of how well the samples themselves are fitted. Bracketing an error nobody can measure put two computable numbers either side of a Gramian’s error; here one computable number, bought with samples, sits within a factor of two below the model’s.

What this changes in practice

A code choosing a Loewner model’s order from data should not stop when σk+1/σ1\sigma_{k+1}/\sigma_1 passes the tolerance. On the eight-pole system at a tolerance of 10−1110^{-11} that stops at order five, eighty-five times over; on the twelve-mode system at 10−1210^{-12} it stops at five with the error at 3.8⋅10−113.8 \cdot 10^{-11}, thirty-eight times over. It should hold back some samples — a handful placed between the existing ones, densest where the function bends, which here is the end of the interval nearest the poles — and stop when the model’s error there passes the tolerance, with a safety factor of about two.

The same holds for the other reductions in this field that are built from values. Exact at the points that were named made the point for moment matching: a model that is exact at chosen points bounds nothing anywhere else. The cure is the same, and it is cheap. And where to put the poles of a rational function suggests where to spend it: geometric clustering toward the hard end, which is where every worst error measured here sat. A point that was not used to build the model is the only kind that can say how the model does where it was not built.

Why a ratio of two is enough

A factor of 2.1 is the whole of the midpoint estimate’s slack, and it has a reason that can be read off the error curve. The error between two samples rises from near the samples’ own misfit to a peak somewhere in the gap and falls again, and a geometric midpoint sits near that peak when the gap is small compared with the function’s scale of variation, as it is everywhere here: 48 samples over 1.7 to 2 decades put consecutive samples 9 to 10 per cent apart. The estimate misses the peak by however far the peak is from the midpoint, which for a smooth bump in a gap is a modest factor and nothing like the orders of magnitude separating the samples’ misfit from the peak. The same reasoning says when it will fail: when the gap is not small compared with the function’s variation, near a resonance on the imaginary axis or with too few samples, the peak need not be near any midpoint.

What four systems do not show

Four smooth functions, sampled on the real axis to the left of their spectra. A real system is sampled on the imaginary axis, where the transfer function’s peaks sit close to the samples and the error between them could be concentrated much more sharply than at the end of a real interval; there the midpoints might need to be denser than one per gap. No noise: the first essay found that noise of 10−1010^{-10} sets a floor on every quantity, and a midpoint estimate built from noisy midpoints would stop at that floor too. One construction, the projected pencil with its alternating split of the samples; a split into two contiguous halves gives a different pencil and could change where the error sits. And orders up to six, where the arithmetic floor ends every comparison: a system with a slower-falling error would test the estimate over more orders.

Still open: the imaginary axis, and a model that chooses its own checks

The imaginary axis. Sampled at frequencies s=iωs = i\omega, a lightly damped system has resonant peaks between the samples. The prediction with a sign is that on the eight-pole system with its poles moved off the real axis to damping ratios of 0.05, the midpoint estimate at one midpoint per gap understates the true error by more than a factor of ten at some order, and that two points per gap bring it back within a factor of three.

Checks the model chooses. AAA chooses its next support point where its residual is largest; 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.

Where this goes next

12345678910111210⁻¹³10⁻¹¹10⁻⁹10⁻⁷10⁻⁵10⁻³10⁻¹10¹order of the reduced modelworst relative errora model made of measurementssamples used24degree read6gap at the cut2·10⁸best model, order12its error2.7·10⁻¹¹states in the original40filled: at its own samplesopen: everywhere else A model with no matrices behind it Twenty-four numbers — values of a transfer function at twenty-four points — produce a sixth-order model of a forty-state system that passes through every sample to 10⁻¹² and matches the function it never saw to 10⁻⁹. Noise of 10⁻¹⁰ on those numbers takes the rank decision's gap from 10⁸ to 4. 12345678910111210⁻¹⁵10⁻¹³10⁻¹¹10⁻⁹10⁻⁷10⁻⁵10⁻³10⁻¹degree of the approximantworst relative error on the target setchosen by the residualdegree reached12error there10⁻¹²even support, same degree1.1·10⁻⁷advantage702nearest support point0.096clustering ratio1.5filled: points chosen by the erroropen: points spread evenly The points the algorithm chose A rational approximant whose support points are picked by its own residual clusters geometrically at a branch point nobody named — recovering by measurement the rule a hand-built approximant is given. At degree ten it is seven hundred and fifty times more accurate than the same form with its points spread evenly. 234567891010⁻⁴10⁻²1ADI steps krelative error, and max|r|²max|r|²the measured errora bound that is an identitysteps drawn9ratio, worst0.77spread of the ratio1error at k = 103.9·10⁻⁴its residual1.1·10⁻⁴four orders of errorand one constant Bracketing an error nobody can measure The error of a low-rank Gramian factor is the one quantity a caller cannot compute, because computing it needs the Gramian the factor exists to avoid forming. Two numbers that can be computed sit either side of it — a rational factor known before the run, and a residual known after — and they stay a factor of four apart across a fourfold change of size.

Shares its objects with

Essays that name at least two of the same things, and that neither author linked.

Named objects

A flat tag is an object no other essay names yet.

Data-driven realisationError estimateInterpolationLoewner matrixMcMillan degreeRational approximationSingular value decompositionTransfer function