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

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.
31 min read 26 figures Where the algebra stops being the arithmeticMeasured, not assertedTwo routes to a number

Worth reading first: A model that is a rational function · Rank is a decision.

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.

Six is the degree, and the model reaches it from 24 numbersA reduced model of a forty-state system built out of 24 samples of its transfer function and nothing else — no state matrix, no eigenvalues, no factorisation of anything the model is of. The Loewner pencil of divided differences is formed from the samples, its singular values are read for the degree (6 here, with a gap of 2.01·10⁸ across the cliff), and the pencil is projected onto that many directions. The filled curve is the worst relative error at the samples the model was built from; the open one is the worst error at forty points it never saw. At the true degree of six the first is 6.99·10⁻¹² — interpolation is exact by construction, not by accuracy — and the second is 2.12·10⁻⁹.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
Fig. 1 A sixth-order model of a forty-state system, built from twenty-four samples and nothing else. The filled curve is the error at the samples it was built from; the open curve is the error at forty points it never saw.

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

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. 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.

Noise at 10⁻⁶ puts the degree at 10 and the model at 0.00164A reduced model of a forty-state system built out of 24 samples of its transfer function and nothing else — no state matrix, no eigenvalues, no factorisation of anything the model is of. The Loewner pencil of divided differences is formed from the samples, its singular values are read for the degree (10 here, with a gap of 6542 across the cliff), and the pencil is projected onto that many directions. The filled curve is the worst relative error at the samples the model was built from; the open one is the worst error at forty points it never saw. With relative noise of 10⁻⁶ on the samples the singular values stop falling to the rounding level, the cliff becomes a slope, and no order does better than 0.00164 — the data's own error is a floor, and taking more of the degree past 11 buys nothing.12345678910111210⁻¹³10⁻¹¹10⁻⁹10⁻⁷10⁻⁵10⁻³10⁻¹10¹order of the reduced modelworst relative error10⁻⁶, the noise in the dataa model made of measurementssamples used24degree read10gap at the cut6542best model, order11its error0.0016states in the original40filled: at its own samplesopen: everywhere else
Fig. 2 The same construction with relative noise of 10⁻⁶ on the samples, where no order does better than 1.6·10⁻³ and the interpolation error is still at the rounding level.

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. Oversampling is free and useful: extra samples do not raise the order — the rank is still the degree — and they make the rank decision better conditioned, which on noisy data is the one lever available.

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.

At other settings

Noise at 10⁻¹² puts the degree at 6 and the model at 1.39·10⁻⁹A reduced model of a forty-state system built out of 24 samples of its transfer function and nothing else — no state matrix, no eigenvalues, no factorisation of anything the model is of. The Loewner pencil of divided differences is formed from the samples, its singular values are read for the degree (6 here, with a gap of 4.73·10⁶ across the cliff), and the pencil is projected onto that many directions. The filled curve is the worst relative error at the samples the model was built from; the open one is the worst error at forty points it never saw. With relative noise of 10⁻¹² on the samples the singular values stop falling to the rounding level, the cliff becomes a slope, and no order does better than 1.39·10⁻⁹ — the data's own error is a floor, and taking more of the degree past 11 buys nothing.12345678910111210⁻¹³10⁻¹¹10⁻⁹10⁻⁷10⁻⁵10⁻³10⁻¹10¹order of the reduced modelworst relative error10⁻¹², the noise in the dataa model made of measurementssamples used24degree read6gap at the cut4.7·10⁶best model, order11its error1.4·10⁻⁹states in the original40filled: at its own samplesopen: everywhere else
Fig. 3 Noise at 10⁻¹², where the cliff has narrowed and the model is still good to 10⁻⁹.
Noise at 10⁻¹⁰ puts the degree at 9 and the model at 1.63·10⁻⁷A reduced model of a forty-state system built out of 24 samples of its transfer function and nothing else — no state matrix, no eigenvalues, no factorisation of anything the model is of. The Loewner pencil of divided differences is formed from the samples, its singular values are read for the degree (9 here, with a gap of 4.41 across the cliff), and the pencil is projected onto that many directions. The filled curve is the worst relative error at the samples the model was built from; the open one is the worst error at forty points it never saw. With relative noise of 10⁻¹⁰ on the samples the singular values stop falling to the rounding level, the cliff becomes a slope, and no order does better than 1.63·10⁻⁷ — the data's own error is a floor, and taking more of the degree past 11 buys nothing.12345678910111210⁻¹³10⁻¹¹10⁻⁹10⁻⁷10⁻⁵10⁻³10⁻¹10¹order of the reduced modelworst relative error10⁻¹⁰, the noise in the dataa model made of measurementssamples used24degree read9gap at the cut4.4best model, order11its error1.6·10⁻⁷states in the original40filled: at its own samplesopen: everywhere else
Fig. 4 At 10⁻¹⁰, where the gap at the cut has fallen from 10⁸ to 4.4.
Noise at 10⁻⁸ puts the degree at 10 and the model at 1.64·10⁻⁵A reduced model of a forty-state system built out of 24 samples of its transfer function and nothing else — no state matrix, no eigenvalues, no factorisation of anything the model is of. The Loewner pencil of divided differences is formed from the samples, its singular values are read for the degree (10 here, with a gap of 6561 across the cliff), and the pencil is projected onto that many directions. The filled curve is the worst relative error at the samples the model was built from; the open one is the worst error at forty points it never saw. With relative noise of 10⁻⁸ on the samples the singular values stop falling to the rounding level, the cliff becomes a slope, and no order does better than 1.64·10⁻⁵ — the data's own error is a floor, and taking more of the degree past 11 buys nothing.12345678910111210⁻¹³10⁻¹¹10⁻⁹10⁻⁷10⁻⁵10⁻³10⁻¹10¹order of the reduced modelworst relative error10⁻⁸, the noise in the dataa model made of measurementssamples used24degree read10gap at the cut6561best model, order11its error1.6·10⁻⁵states in the original40filled: at its own samplesopen: everywhere else
Fig. 5 At 10⁻⁸, where the best model is four orders worse than the clean one.
Noise at 0.001 puts the degree at 12 and the model at 1.85A reduced model of a forty-state system built out of 24 samples of its transfer function and nothing else — no state matrix, no eigenvalues, no factorisation of anything the model is of. The Loewner pencil of divided differences is formed from the samples, its singular values are read for the degree (12 here, with a gap of ∞ across the cliff), and the pencil is projected onto that many directions. The filled curve is the worst relative error at the samples the model was built from; the open one is the worst error at forty points it never saw. With relative noise of 0.001 on the samples the singular values stop falling to the rounding level, the cliff becomes a slope, and no order does better than 1.85 — the data's own error is a floor, and taking more of the degree past 4 buys nothing.12345678910111210⁻¹³10⁻¹¹10⁻⁹10⁻⁷10⁻⁵10⁻³10⁻¹10¹order of the reduced modelworst relative error0.001, the noise in the dataa model made of measurementssamples used24degree read12gap at the cut∞best model, order4its error1.9states in the original40filled: at its own samplesopen: everywhere else
Fig. 6 And at a tenth of a per cent, where no order produces a usable model.
Two rank decisions for one integer: the McMillan degree of a 40-state modelSingular values, each divided by the largest of its own set. The Loewner matrix is built from 32 samples of H on the imaginary axis and never touches A: its 6th singular value stands 9.46·10⁴ above the next, so the rank decision is not a close call. The Hankel matrix of Markov parameters CAᵏB is built from A, B and C — strictly more information — and has no cliff, because ‖A‖ is (n+1)² and its entries span 3·10³⁷. It returns 2 where the answer is 6. The definition is not the algorithm, and the route with less information is the one that works.1357911131510⁻¹⁸10⁻¹⁵10⁻¹²10⁻⁹10⁻⁶10⁻³1indexsingular value ÷ the largestdegree 6from samples of Hfrom A, B and Cone integer, two routesstates40Loewner rank6its gap9.5·10⁴Markov rank2its dynamic range3·10³⁷the rank of a divided differenceand the rank that cannot be seen
Fig. 7 The rank that reads the degree, in the essay that measured it first.
The support points found a branch point nobody named, and cluster at it by 1.52Approximating γ√(λ + c) on an interval reaching to within 0.03 of its branch point, by a barycentric rational whose support points are chosen greedily: at each step the sample where the current approximation is worst becomes the next support point. The filled curve is the error against degree; the open one is the same barycentric form with support points spread evenly along the interval, which is the choice anybody makes who has no reason to make another. At degree 10 the adaptive rule is 702.3× better, and the advantage grows with the degree. Nothing in the algorithm is told that the function has a branch point or where it is: the points it selects lie between 0.0962 and 9.53 from it, a span of 99.06, with a mean ratio of 1.519 between consecutive distances — the geometric clustering the hand-built approximant in this field is given as a rule, arrived at from the residual.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
Fig. 8 A different construction from the same Loewner matrix, whose noise failure has the same shape.
Singular values of a rank-4 matrix with noise of relative size 10⁻⁸Ten singular values on a logarithmic axis. The first four sit near one; the rest sit at the noise level, and the vertical distance between the two groups is the evidence for the rank.1234567891010⁻¹⁸10⁻¹⁵10⁻¹²10⁻⁹10⁻⁶10⁻³1indexsingular valuecutoff, σ₁ · 10⁻¹⁰numerical rank 10gap 8.2·10⁶an opiniontrue rank 410×10, built with 4 nonzero valuesrank is a decision
Fig. 9 A threshold turning a spectrum into an integer, in the field that named it.
A randomised rank-k solve, 4 seeds a rankRelative error against the rank kept, on a logarithmic vertical axis, with a vertical bar at each rank spanning the seeds. The band is 1.68 wide at rank 8, where the method is at its worst, and 1.004 wide at rank 20, where it is at its best. The same computation on the same data returns a different answer each time.0481216202428321rank keptrelative errormedianthe answer movesspread at rank 81.7spread at rank 201best median error0.15widest where the method is worstand the bound does not say so
Fig. 10 A rank reported as a band rather than a number.
Two bases for one subspace, at an interpolation point of 20κ₂ of the moment basis (σI − A)⁻¹B, (σI − A)⁻²B, … at a single point σ = 20, against κ₂ of the same number of solves spent at points spread across the model's poles and orthogonalised as they are built. The moment basis runs 10.9 to 7.75·10⁹, multiplying by about 65.8 per vector, while the spread one stays at 1. The reason is the one behind every Krylov basis on this site: repeated application of one operator converges to a direction, so every new moment is nearly the last. The rate is set by the ratio of distances from σ to the poles, which means the user's own choice of interpolation point decides how fast their basis becomes unusable — and spending the solves at different points instead is both a better approximation and a basis that survives.2345678110²10⁴10⁶10⁸10¹⁰vectors in the basisκ₂ of the basisas derivedsolves spread outone subspace, two spanning setseight moments at one point7.7·10⁹eight points, spread1growth per vector661/u4.5·10¹⁵the same subspaceand only one of them usable
Fig. 11 Interpolation at points somebody named, which this construction replaces with points somebody measured.
Hankel singular values, and the half of them a product cannot seeσₖ divided by σ₁, for a 20-state model of McMillan degree 12. The lower curve is the square-root route — Cholesky-like factors of the two Gramians and one SVD of RᵀS — and it descends to the level at which the Gramians were computed. The upper curve eigendecomposes the product PQ, agrees for the first 6 values and then stops, flattening at 9.67·10⁻¹⁰. The dashed line is σ₁√u = 2.32·10⁻⁹, which is where a route that squares the conditioning must stop: the eigenvalues of PQ are σ², an absolute error of uσ₁² on them is a relative error of √u on σ. κ(PQ) is 1.8·10³⁶ against 1.34·10¹⁸. Every σ below the line is a term in the error bound, so this is not an academic loss.1357911131510⁻¹⁸10⁻¹⁵10⁻¹²10⁻⁹10⁻⁶10⁻³1index kσₖ ÷ σ₁σ₁√utwo factorstheir productthe normal equations, againroutes agree to k =6product floors at9.7·10⁻¹⁰σ₁√u2.3·10⁻⁹κ(P)κ(Q)1.8·10³⁶√ of it1.3·10¹⁸do not form the productthe σ below the line are the bound
Fig. 12 The singular values a reduction with an error bound is truncated at, which this construction has no access to.
Backward and forward error against the condition numberA log–log plot over twelve decades of condition number. The backward error is a flat line at ten to the minus sixteen; the forward error rises in proportion to the condition number.110²10⁴10⁶10⁸10¹⁰10¹²10¹⁴10⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹condition number κ(A)relative errorforward errorbackward errorpredicted: κ · u8×8, 20 seeds per κ; dashed is the worstthe problem worsens, not the method
Fig. 13 The two errors, of which the interpolation error here is neither.
What the solver reported, and what it achievedFinal relative residual and final relative error against the working precision, on a logarithmic vertical axis. The residual sits at the stopping tolerance at every precision — the method converged, by every test available to it. The error runs from 1.4·10⁻¹³ at 53 bits to 0.00597 at 8, tracking the unit roundoff, which is drawn beside it.614223038465410⁻¹³10⁻¹⁰10⁻⁷10⁻⁴working significand bitsrelative sizeerrorresidualuthe gap the solver cannot seeerror ÷ residual at 24 bits1.5·10⁶error ÷ residual at 16 bits2.3·10⁷error ÷ residual at 8 bits2.4·10¹⁰every convergence test passesand the answer is wrong in proportion to u
Fig. 14 A residual that says nothing about the error, which is exactly this essay’s diagnostic problem.
Conjugate gradients on an ill-posed problem at 0.10% noiseTwo curves against the step count on a logarithmic vertical axis. The relative residual falls at every one of the 120 steps without exception. The error against the true signal falls to 0.1050 at step 44 and then climbs, reaching 0.56 by the end — 5.3 times its best value.015304560759010512010⁻³10⁻²10⁻¹1steprelative sizeleast error: 44discrepancy stop: 27errorresidualthe knob is an integerleast error, at step44error there0.11error at step 1200.56the residual falls at every stepthe error turns and keeps rising
Fig. 15 An error that falls and then rises as more is asked for.
What the regularisation costs, and what 6 steps of refinement take backSolving [[H + δI, Aᵀ], [A, −δI]] instead of K gives the exact answer to a different problem, so its error is proportional to δ: measured at 1489·δ across six decades, which is a slope of one and not a trend. Refining against the unregularised matrix — the residual formed with K and the correction solved with the regularised factorisation — removes that term entirely, because the perturbation was never in the residual. It works while δ is below σ_min(K) = 6.797·10⁻⁴, marked on the axis, and stops working above it: the iteration's contraction factor is δ/σ_min and a fixed point needs that under one. So the trade-off curve every regularisation essay on this site has drawn — a term falling in δ against a term rising in it — has, here, a left branch that can simply be removed.-14-12-10-8-6-4-210⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹log₁₀ δrelative error against the exact answerδ = σ_min(K)no refinement6 stepsa left branch that can be removederror ÷ δ, unrefined1489σ_min(K)6.8·10⁻⁴refined at δ = 10⁻⁶3.6·10⁻¹⁶refined at δ = 10⁻²0.62the perturbation is known exactlybecause the code chose it
Fig. 16 A parameter deciding an answer the data does not, which the order here becomes.
The L-curve, and where four rules put λThe norm of the solution against the norm of its residual, on logarithmic axes, as λ sweeps eight decades. The curve has a corner: to the left of it the noise is being amplified and to the right the signal is being thrown away. Four points are marked — the three rules that use only the data, and the oracle, which requires the exact answer and is not a method.10⁻¹110¹10²10³10⁴10⁵‖Ax − b‖‖x‖the oraclediscrepancyL-curvegeneralisedscored against a truth none hasoracle, relative error0.14discrepancy principle, as a multiple1L-curve corner, as a multiple1.3generalised cross-validation, as a multiple1the oracle needs the exact answer and is not a methodit is the reference the others are scored on
Fig. 17 The published rules for choosing such a parameter, scored against a truth.
What a rank-10 approximation can achieve, by spectrumA semi-logarithmic plot of singular value against index for three spectra — geometric decay, algebraic decay, and flat — with the rank-ten approximation error marked on each.1112131415110⁻³10⁻²10⁻¹1index jσⱼσ11geometric, 0.85ʲalgebraic, j^−1.00flatmeasured, and equal to σ₁₁rank-10 error, geometric0.2rank-10 error, algebraic0.091rank-10 error, flat160×60, spectrum chosen rather than the entriesthe matrix decides, not the method
Fig. 18 A spectrum with no gap, which noisy Loewner singular values become.
Error of the best rank-k approximation to a 12×12 matrixApproximation error against k on a logarithmic axis, with the measured error and the next singular value drawn as separate curves lying exactly on top of one another.123456789101110⁻¹⁰10⁻⁸10⁻⁶10⁻⁴10⁻²1rank k of the approximation‖A − Aₖ‖measured, 2-normσₖ₊₁, from theorymeasured, Frobeniusthe first two agreeto 2·10⁻⁸worst |‖A−Aₖ‖₂ − σₖ₊₁| / σₖ₊₁2·10⁻⁸worst Frobenius discrepancy2·10⁻⁸κ = 10⁹; 30 random rank-3 matrices, none closerthe error is σₖ₊₁
Fig. 19 An approximation error against the rank kept.
How close Hager's estimate is to the true κ₁, over 200 seeded 8×8 matricesFive bars. The estimate is exactly the true condition number on 81% of the sample and inside ten per cent on 87%; the worst underestimate in the whole sample returns 38% of the truth. The last bar is the matrix built to defeat it, at 7.7%, well below anything the random sample reached.each bar is a percentage — of the sample, or of the true condition numberexactly right80.5%inside 10%87.0%inside a factor of 287.0%worst in the sample, ×10037.7%the constructed matrix, ×1007.7%usually exactexact share0.81worst of the sample0.38the constructed matrix0.077a routine that is right most of the timeand never wrong in the safe direction
Fig. 20 An estimate reported with its spread.
The exact solution of a 13×13 Hilbert system beside the computed oneTwo columns of numbers: the exact answer, which is the integers one to thirteen, and the answer double-precision elimination returns, with the number of correct digits beside each.H13 x = b, b formed exactly so that x = (1, 2, …, 13)12345678910111213exact1.00002.00003.00133.97865.18864.993010.46720.048921.2657-2.575319.21478.906013.5113computed6.6 correct digits4.8 correct digits3.4 correct digits2.3 correct digits1.4 correct digit0.8 correct digitno correct digitsno correct digitsno correct digitsno correct digitsno correct digits0.6 correct digit1.4 correct digitbackward error2.2·10⁻¹⁷κ = 1.7·10¹⁸The algorithm solved a neighbouring problem perfectly. That problem's answer is this one.right-hand side built in BigInt rationalsthe truth is known
Fig. 21 An answer that is known, which the modal check in this essay also has.
How much a perturbation of the right-hand side is amplified, κ = 10⁶The cumulative distribution of the amplification factor over two hundred random perturbation directions, with the condition number marked as the upper limit.110¹10²10³10⁴10⁵10⁶10⁷00.250.50.751amplification of the input perturbationfraction of directions at or belowκ = 10·10⁵worst found 7.6·10⁵6×6, 200 directionsmedian reaches 0.29 of κ
Fig. 22 The amplifier belonging to the problem.
Backward error, forward error and the condition number, with measured valuesTwo boxes at the top — the problem posed and the nearby problem the algorithm answered exactly — and two answers below them, with the distances between all four labelled by numbers from a Hilbert solve.the problem you posedA = H10b = A·(1, 2, …, 10)the problem it answered exactlyA + δA, b + δb‖δ‖ / ‖A‖ = 2.3·10⁻¹⁷the answer you wantedx = (1, 2, …, 10), exactlythe answer you gotx̂, wrong by 2.7·10⁻⁴ relativebackward error 2.3·10⁻¹⁷forward error 2.7·10⁻⁴κ = 1.6·10¹³κ · η = 3.6·10⁻⁴, and the measured forward error is 2.7·10⁻⁴.The algorithm is not at fault. The problem is.H10, LU with partial pivotingresidual and error differ
Fig. 23 The identity that has no term for an error in the data.
Forward error of a 6×6 Hilbert solve at eight precisionsA bar for each significand width from 12 to 53 bits showing the relative error in the computed solution, with the condition number times the unit roundoff marked as a prediction.κ = 1.5·10⁷ · the exact answer is (1, 2, …, 6)12 bits3.316 bits0.8620 bits4.324 bits0.09330 bits7.6·10⁻⁴36 bits6.2·10⁻⁶43 bits2·10⁻⁷53 bits4.5·10⁻¹¹dashed: κ · unit roundoffone matrix, eight arithmeticsmeasured against a known answer
Fig. 24 The precision knob, which does not move a noise floor.
Growth factor under partial pivoting to n = 40: the bound, the worst case, and realityGrowth factor against matrix size on a logarithmic vertical axis. The two-to-the-n bound rises as a straight line; Wilkinson's matrix sits exactly on it; random matrices stay near one.0816243240110²10⁴10⁶10⁸10¹⁰10¹²10¹⁴matrix size ngrowth factor max|u| / max|a|the 2ⁿ⁻¹ boundworst of 30 randommedian randomWilkinson's matrix sits on the bound30 Gaussian matrices per sizeat n = 40: bound 5.5·10¹¹, worst 4.8
Fig. 25 A bound never attained, in the field that measured it first.
The same 16 answers, backward stable under one reading and not under anotherEvery real eigenpair of an overdamped chain of 8 masses, computed once through the first companion linearisation, and its backward error measured four ways — against perturbations of all three coefficients, of K alone, and of M alone. The unrestricted number is the one every published formula computes and it sits at the rounding level across the spectrum. Restricting the perturbation charges the same residual to less data, so it can only raise the quotient, and the factor has a closed form: (|λ|²‖M‖ + |λ|‖C‖ + ‖K‖) divided by the allowed part, checked here against the two computed numbers to 10⁻⁸. Because it carries |λ|², it is not a property of the polynomial: allowing only K to move costs between 1.06 and 45.69, and allowing only M costs between 1.99 and 6367. At a tolerance of 10⁻¹⁵ 15 of these pairs are stable against the whole polynomial and 7 against K alone.10⁻²10⁻¹110¹10⁻¹⁷10⁻¹⁵10⁻¹³10⁻¹¹|λ|backward error of the computed pairone residual, four denominatorspairs measured16worst, all three1.4·10⁻¹⁵worst, K alone1.8·10⁻¹⁴worst, M alone9·10⁻¹²K-only factor, low1.1high46solid: every coefficient may movedashed: only one may
Fig. 26 A backward error under a restriction, in the previous essay of this phase.

Where this goes next

10⁻¹110¹10²10³10⁻¹¹10⁻¹⁰10⁻⁹10⁻⁸10⁻⁷10⁻⁶10⁻⁵10⁻⁴10⁻³10⁻²10⁻¹frequency ω|H(iω)|one function, two routesstates24relative, ω ≤ 10010⁻¹⁴against the terms3.9·10⁻¹⁵‖H‖∞0.0924 statesone function of one variable A model that is a rational function A state matrix has a hundred thousand rows and the thing anyone wants from it is a function of one complex variable. The number that says how much of that size was ever the complexity is a rank — and the rank a derivation writes down cannot be computed, while one built from samples alone can. 2345678110²10⁴10⁶10⁸10¹⁰vectors in the basisκ₂ of the basisas derivedsolves spread outone subspace, two spanning setseight moments at one point7.7·10⁹eight points, spread1growth per vector661/u4.5·10¹⁵the same subspaceand only one of them usable A basis that is the same subspace and not the same thing The interpolation conditions are conditions on a subspace, so any basis of it will do. The one a derivation writes down reaches a condition number of 7.7·10⁹ by its eighth vector, and the rate at which it gets there is set by a number the user chose with no information. 1234567891010⁻¹⁸10⁻¹⁵10⁻¹²10⁻⁹10⁻⁶10⁻³1indexsingular valuecutoff, σ₁ · 10⁻¹⁰numerical rank 10gap 8.2·10⁶an opiniontrue rank 410×10, built with 4 nonzero valuesrank is a decision Rank is a decision A floating-point matrix does not have a rank. It has a spectrum of singular values, and somewhere in that spectrum is a place where the values stop being signal and start being noise. Deciding where is a judgement, and the evidence for it is a gap.

What links here

Computed from the collection, not written here: the essays that point at this one.

Named objects

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

Adaptive interpolationData driven realisationDescriptor systemInterpolationLoewner matrixMcMillan degreeNoise floorRank is a decisionTransfer function