Numerical Linear Algebra
where the algebra stops being the arithmetic
The eigenvalue problem that is not linear

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

Worth reading first: Why a Gramian can be truncated at all · A problem with infinitely many eigenvalues · The projection and the right angle.

A nonlinear eigenvalue problem is solved by replacing the nonlinearity with a rational function and solving the linearised problem that results. The error that substitution commits is committed before any arithmetic and is invisible to every residual the computation produces afterwards, so how good the approximation is decides how good the answer is.

The approximants this field has used so far are built by a rule. Poles clustered geometrically towards the branch point, at a ratio the author sets, because that is where the function is not analytic and where an approximation is hardest — and the contest between a good rule and a plausible one is worth five orders of magnitude.

Both rules need something nobody has: where the singularity is.

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. 1 The error against degree for support points chosen greedily by the residual, and for the same barycentric form with points spread evenly along the interval.

The algorithm

The adaptive Antoulas–Anderson algorithm keeps a rational function in barycentric form,

r(z)  =  Σⱼ wⱼ fⱼ /(z − zⱼ)   ÷   Σⱼ wⱼ /(z − zⱼ)

which interpolates f at every support point zⱼ automatically, whatever the weights are. So the weights are free to be chosen for accuracy everywhere else.

At each step it does two things. It adds as a new support point the sample where the current approximation is worst — the greedy choice, and the only place the function’s own behaviour enters the construction. Then it re-solves for the weights: a least-squares problem over the samples that are not support points, whose matrix is a Loewner matrix of divided differences and whose solution is that matrix’s smallest right singular vector.

Nothing in it is told that the function has a branch point, or where. It is handed a table of values.

What it finds

The target function is γ√(λ + c), whose branch point is at −c, on an interval that reaches to within a stated fraction of it. Running the greedy to degree twelve and sorting the chosen support points by their distance to the branch point:

reach   degree   error         nearest point   mean ratio between   span
                                to branch      consecutive          
0.05      12     8.3·10⁻¹⁴      1.60·10⁻¹          1.450             59
0.02      12     6.1·10⁻¹²      6.41·10⁻²          1.576            150
0.008     12     1.5·10⁻¹⁰      2.57·10⁻²          1.712            370

The distances to the branch point form a geometric progression, with a ratio the algorithm discovered rather than was given, and the ratio rises as the target set reaches further in — 1.45, 1.58, 1.71 — because a harder problem needs its points spread over a wider range of scales.

That is the same construction clusteredPoles is handed as a rule, arrived at from the residual. The algorithm has located a singularity it was never told about, and has resolved it at the rate the theory says is right.

Why the barycentric form and not a ratio of polynomials

The form the algorithm keeps is not incidental, and it is the reason the whole construction is numerically usable at all.

A rational function can be written as p(z)/q(z) with p and q in a monomial basis, and that representation is a disaster on a clustered set of points: the coefficient matrix is a Vandermonde matrix whose condition number grows exponentially with the degree, so recovering p and q from interpolation conditions loses every digit by degree ten or so. This site has an essay about that and the point transfers exactly.

The barycentric form has no coefficients in that sense. It carries the function’s values at the support points, which are data, and a set of weights whose scale is irrelevant — only their ratios matter, so they can be normalised to unit length and the whole representation is bounded.

Two consequences are worth having. Interpolation is exact by construction: at a support point the formula’s numerator and denominator both have a pole, they cancel, and the value returned is fⱼ to the last bit, whatever the weights are. And the conditioning of the representation does not degrade as the points cluster, which is precisely the regime the greedy drives it into.

So the choice of form and the choice of points are the same design decision seen twice. An algorithm that clusters its points has to have a representation that survives clustering, and the barycentric form is the one that does.

What it is worth

Against the same barycentric form with its support points spread evenly along the interval — the choice anybody makes who has no reason to make another:

degree    adaptive        even            advantage
  4       3.6·10⁻⁴       1.5·10⁻³            4.1×
  6       1.9·10⁻⁶       2.7·10⁻⁵           15×
  8       1.0·10⁻⁸       5.8·10⁻⁷           59×
 10       2.1·10⁻¹¹      1.3·10⁻⁸          620×

and reaching further towards the branch point widens the gap: at reach 0.02 the degree-ten advantage is 750×.

The advantage grows with the degree, which is the signature of a difference in rate rather than in constant. Both errors fall geometrically; the adaptive one falls faster, so the ratio compounds. Six degrees of separation buys two and a half orders.

Where it buys nothing

At degree three, on a target set that stops well short of the branch point, the even support is three per cent better.

That is the honest end of the result and it is asserted as such: the adaptive rule is only required to win from degree five upwards, and the reason is that a shallow target set has no range of scales for the clustering to resolve. The function is smooth over the whole interval, an even spread is a perfectly good spread, and the greedy step is spending its choices on nothing.

An adaptive rule buys nothing on a problem with no scale in it, which is worth stating because the reverse is easy to assume. The gap in the table above is a property of the problem — the reach towards the branch point — and not of the algorithm’s cleverness.

The support points found a branch point nobody named, and cluster at it by 1.29Approximating γ√(λ + c) on an interval reaching to within 0.3 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 1934× 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.962 and 9.53 from it, a span of 9.906, with a mean ratio of 1.29 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 reached10error there10⁻¹⁵even support, same degree2·10⁻¹²advantage1934nearest support point0.96clustering ratio1.3filled: points chosen by the erroropen: points spread evenly
Fig. 2 A shallow target set, where the two curves nearly meet and the adaptive rule is buying almost nothing.

The greedy step is where the noise gets in

The rule is “add the sample where the error is largest”, and on data with error in it that eventually means “add the sample with the largest error in it”. The consequence is measurable and it is not a plateau.

noise on the samples   floor reached   at degree   error at degree 14
none                   7.4·10⁻¹⁵          14         7.4·10⁻¹⁵
10⁻¹⁰                  7.8·10⁻¹⁰           9         8.3·10⁻⁸
10⁻⁶                   6.6·10⁻⁶            6         1.6·10⁻³
10⁻³                   3.5·10⁻³           12         1.4·10⁻¹

With no noise the error falls monotonically for fourteen steps to the rounding level. With relative noise of 10⁻⁶ it reaches 6.6·10⁻⁶ at degree six — the noise level, as expected — and then gets worse, ending at 1.6·10⁻³ eight degrees later. Two and a half orders of degradation, bought with eight further degrees of work.

The mechanism is direct. Past the noise floor the largest residual is at whichever sample happened to be perturbed most, so the algorithm places a support point there, forces the approximant through a corrupted value exactly, and the resulting function is worse everywhere else. Every subsequent step compounds it.

That is a sharper failure than the usual one. A method that stalls at the noise level has stopped being useful; a method that degrades past it has become actively harmful, and the symptom is invisible from inside — the interpolation error at the support points is exactly zero throughout, because interpolation is what the barycentric form does.

What to do about that

The published remedy is a stopping rule: stop when the error stops falling, which requires knowing when it has stopped, which on noisy data is the same question again. In practice one of three.

Stop at an estimate of the noise. If the data’s error is known — a measurement tolerance, a solver’s own accuracy — stop when the approximation reaches it. Cheap and correct, and requires a number nobody always has.

Watch for the turn. Keep the best approximant seen rather than the last, which converts the degradation into a plateau. Costs nothing and does not prevent the wasted degrees.

Or regularise the weight solve. The least-squares step’s matrix becomes ill conditioned as the support points cluster, and truncating its small singular values damps the response to noise. That is the regularisation field’s machinery applied one level down, and it introduces a parameter of its own.

None of them is free and the first is the honest one. The general point is the one the nonlinear eigenvalue essays keep arriving at: an approximation that is better than its data is not better than its data.

The spurious pole, and why it is not in this measurement

A rational approximant built this way is known to produce Froissart doublets — a pole and a zero almost on top of each other, with a residue at the rounding level, contributing nothing to the function and appearing in the pole list.

They matter for a linearisation, because the linearised problem’s spectrum contains the approximant’s poles, and a spurious pole is an eigenvalue of the linearisation that is an eigenvalue of nothing. The eigenvalues that are answers to nothing is the essay about that, in the context of a hand-built approximant where the poles are placed by rule and the spurious ones are the ones off the target set.

This essay does not measure them, and the reason is worth stating rather than leaving as an omission: locating the poles of a barycentric rational requires solving a generalised eigenvalue problem of its own, and the measurement here is of the approximation error, which is a function of the values rather than of the pole locations. An approximant with a doublet in it can have an excellent error curve, which is exactly why the doublets are a problem for the eigensolve and not for the fit.

What the noise measurement does show is the mechanism that produces them. Past the noise floor the algorithm places support points to interpolate corrupted values, and the weights it computes to do that are the ones that produce near-cancelling pole–zero pairs. So the degradation in the table and the appearance of doublets are the same phenomenon reported two ways, and a stopping rule that prevents one prevents the other.

Why this matters for a linearisation

The reason to care about degree rather than only about error is that the degree is a size.

Linearising a rational approximant of degree d on an n × n problem produces a matrix of size roughly dn, and the eigensolve costs the cube of that. So a degree-ten approximant costs eight times a degree-five one, and the choice of support points decides which is needed for a given accuracy.

From the table, reaching 10⁻⁸ needs degree eight adaptively and degree twelve or more evenly — which, cubed, is a factor of about three in the eigensolve. The support points are chosen at negligible cost, before any matrix is built, and they decide the size of the matrix that gets built.

That is the same accounting the approximant-choice essay makes about clustered against even poles, and the adaptive version is the one that does not need to be told where to cluster. On a nonlinearity given as a table of measured values — which is what a real nonlinear eigenvalue problem often is — it is the only one of the two that can be run at all.

What the two routes here are

The site’s habit is a second route, and this essay has one that is worth naming because it is not a second implementation.

The adaptive approximant and the hand-clustered one are two constructions of the same object, arrived at from opposite directions: one from a theory of where the function is hard, one from a measurement of where the error is. They agree — the geometric ratio the greedy discovers, 1.45 to 1.71, is in the range a hand-built approximant would be given — and the agreement is the check.

If they had disagreed, one of two things would have been true: either the theory’s rule is not the right one for this function, or the greedy is chasing something other than the singularity. Neither is a numerical error, and neither would have been visible in a residual. That is what makes the comparison worth making at all — it tests the modelling step, which is the one no residual on this site can see.

What the least-squares step is doing

One step in the loop deserves a closer look, because it is the only linear algebra in the algorithm and it is a familiar object wearing a new name.

Having fixed the support points, the weights are chosen to make the approximant as close as possible to the data at the remaining samples. Writing the interpolation conditions and clearing denominators turns that into a linear least-squares problem in the weights, whose matrix has entries

Aᵢⱼ  =  (fᵢ − fⱼ) / (zᵢ − zⱼ)

— a Loewner matrix of divided differences, exactly the object the model-reduction essays build a state-space realisation from. The same matrix, from the same data, used for a different purpose.

The solution wanted is the null vector, or nearest to it: the weight vector minimising ‖Aw‖ subject to ‖w‖ = 1, which is A’s smallest right singular vector. So each greedy step costs one SVD of a matrix whose height is the number of samples and whose width is the current degree — a few hundred by twelve at the end, which is nothing.

Two things follow. The rank structure of A is the degree structure of the problem: if the function being approximated were exactly rational of degree d, A would have rank d and the null vector would be exact, which is the observation the Loewner essay turns into a degree detector. And the smallest singular value is a diagnostic: it is the residual of the best fit at this degree, available without evaluating the approximant, and it falls as the greedy makes progress and stops falling when the noise floor is reached.

What this changes about how a nonlinearity is given

A closing observation about inputs, because it is the practical reason the adaptive route exists.

A hand-built approximant needs a formula. Somebody has to know that the nonlinearity is γ√(λ + c), locate the branch point at −c, and place poles clustering towards it. Every step of that requires the closed form.

Real nonlinear eigenvalue problems are often not given that way. A delay term, a boundary condition from an exterior solve, a frequency-dependent material property from a measurement — these arrive as a routine that returns a value, or as a table, and there is no formula to differentiate or singularity to locate.

The greedy needs only values. Given a table it produces an approximant whose support points cluster wherever the function is hard, without anybody knowing where that is or whether there is a branch point at all. That is a change in what the input can be, and it is a larger change than a factor of seven hundred in accuracy.

The price is the one the noise table records. A formula has no error in it; a table has whatever error the process that produced it has, and the greedy will chase that error past the point where it helps. So the two routes trade a requirement for a hazard: the rule needs knowledge nobody may have, and the adaptive route needs a stopping criterion nobody can compute from the data alone.

At other settings

The support points found a branch point nobody named, and cluster at it by 1.35Approximating γ√(λ + c) on an interval reaching to within 0.15 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 566.4× 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.481 and 9.53 from it, a span of 19.81, with a mean ratio of 1.348 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 reached11error there5.4·10⁻¹⁵even support, same degree4.5·10⁻¹¹advantage566nearest support point0.48clustering ratio1.3filled: points chosen by the erroropen: points spread evenly
Fig. 3 Reaching to fifteen per cent of the branch point, where the advantage is starting to open.
The support points found a branch point nobody named, and cluster at it by 1.45Approximating γ√(λ + c) on an interval reaching to within 0.05 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 616.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.1603 and 9.53 from it, a span of 59.44, with a mean ratio of 1.45 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 there8.3·10⁻¹⁴even support, same degree1.3·10⁻⁸advantage616nearest support point0.16clustering ratio1.4filled: points chosen by the erroropen: points spread evenly
Fig. 4 Five per cent, where the degree-ten advantage is a factor of six hundred.
The support points found a branch point nobody named, and cluster at it by 1.58Approximating γ√(λ + c) on an interval reaching to within 0.02 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 747.2× 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.06413 and 9.53 from it, a span of 148.6, with a mean ratio of 1.576 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 there6.1·10⁻¹²even support, same degree4.7·10⁻⁷advantage747nearest support point0.064clustering ratio1.6filled: points chosen by the erroropen: points spread evenly
Fig. 5 Two per cent, where it is seven hundred and fifty.
The support points found a branch point nobody named, and cluster at it by 1.71Approximating γ√(λ + c) on an interval reaching to within 0.008 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 812.2× 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.02565 and 9.53 from it, a span of 371.5, with a mean ratio of 1.712 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 there1.5·10⁻¹⁰even support, same degree7·10⁻⁶advantage812nearest support point0.026clustering ratio1.7filled: points chosen by the erroropen: points spread evenly
Fig. 6 And eight parts in a thousand, where the chosen points span a factor of three hundred and seventy.
Two approximants at equal linearisation size, at a reach of 0.06‖g − r‖ against the size of the matrix the eigensolver is then handed: n(m+1) for a rational approximant with m poles and nd for a polynomial of degree d, so d = m + 1 is the fair comparison. At this target set the ratio of the polynomial's error to the rational's runs 1.2, 2.9, 8.0, 23.7, 73.2 across the five sizes. On a target set that stands clear of the branch point the two are the same measurement and the ordering flips with the noise; on one that reaches it the rational pulls away, because a polynomial converges algebraically against an algebraic singularity. Neither method is better — the question is different, and which question it is was decided by whoever said where to look.263850627410⁻⁵10⁻³10⁻¹rows in the linearisation‖g − r‖ on the target setpolynomialrationalequal cost, two basesreach0.06largest size drawn78rational there2·10⁻⁵polynomial there0.0015the ratio73no winner on an easy targetand two orders on a hard one
Fig. 7 The hand-built rule this algorithm rediscovers, in the essay that introduced it.
Three residuals, and the one a solver returns is the one about nothingFor a rational approximant with m poles, on a target set reaching to 0.3: ‖T̃(λ)x‖ — the residual against the problem the eigensolver was handed — stays at 1.2·10⁻¹⁵ to 1.5·10⁻¹² — it does not fall with m, and past a point it slowly rises, because each added pole makes the fit's own basis worse conditioned; ‖T(λ)x‖, which costs one further evaluation of γ√(λ + c), falls with the approximation; and the forward error against the closed form falls with it, staying 20.5 to 24.8 times larger. The free number rises by 1306 across the sweep while the answer improves by 885. So the free number says nothing and the nearly-free number says almost everything. Evaluate the residual against the function you asked about, not against the one you handed over is the whole practical content of this field.23456710⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²poles in the approximantresiduals and errorforward error‖T(λ)x‖‖T̃(λ)x‖one extra evaluationagainst the approximant1.5·10⁻¹²against the problem asked1.9·10⁻⁶forward error4.8·10⁻⁵‖g − r‖ there8.5·10⁻⁵the free residual is flatand the answer is not
Fig. 8 The error committed before any arithmetic, beside the two the computation reports.
The 16 eigenvalues of an 8 × 8 quadratic eigenvalue problem, computed and in closed formλ²M + λC + K for a chain of 8 masses with C = 0.3M + 0.1K. The crosses are the closed form — one scalar quadratic per eigenvalue of K, whose roots are known exactly — and the discs are the eigenvalues a real Schur factorisation returns from the 16 × 16 first companion linearisation. There are 16 of them for a matrix with 8 rows, of which 16 are complex and arrive in conjugate pairs, so the eigenvectors cannot be independent: 9 vectors in 8 dimensions never are. The worst disagreement between the two routes is 2.12·10⁻¹⁵, and the routes share nothing but the three coefficient matrices.00.0939569-2-1.24498-0.4899610.2650581.020081.7751real partimaginary parttwo routes, one spectrumeigenvalues16rows8complex16against the closed form2.1·10⁻¹⁵n rows and 2n eigenvaluesso the eigenvectors are not a basis
Fig. 9 The eigenvalue problem a linearised approximant is solved as.
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. 10 A different construction from the same Loewner matrix, three essays along.
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. 11 And the same construction on noisy data, where the failure has the same shape as this one.
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. 12 A parameter deciding an answer the data does not, which a stopping degree here also is.
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. 13 The published rules for choosing such a parameter, scored against a truth.
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. 14 An error that falls and then rises, which is exactly the shape of the noisy runs here.
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. 15 A threshold deciding an integer, which a stopping degree is a version of.
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. 16 An estimate reported with its spread.
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. 17 A residual that says nothing about the error, which the interpolation error here does not.
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. 18 An answer that is known, against which this approximation is scored.
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. 19 The amplifier belonging to the problem, which the linearisation inherits.
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. 20 The two errors, of which the approximation error is a third.
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. 21 The identity that has no term for an approximation committed before it.
Relative error of two algebraically identical expressions for (1 − cos x)/x², in binary32A log–log plot of relative error against x. The expression written as it reads loses accuracy as x falls and is eventually wrong in every digit; the rearranged form stays at rounding level.10⁻⁸10⁻⁶10⁻⁴10⁻²110⁻⁹10⁻⁶10⁻³1xrelative error of the computed value(1 − cos x)/x², as written2 sin²(x/2)/x²no digits left at allbinary32 throughoutone function, two spellings · zero below 3.5·10⁻⁴
Fig. 22 Cancellation, which the barycentric form is chosen to avoid.
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. 23 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. 24 A bound never attained, in the field that measured it first.
The accuracy asked for, against the accuracy obtained, for two kernels on one partitionε is applied to each block against that block's own largest singular value; the error is then reported against the whole matrix's norm. The two are not the same number and the dashed diagonal is where they would be. For 1/r the obtained error is 29 times smaller than the tolerance asked for; for log r on the identical partition it is 13 times smaller. Both curves are straight and parallel to the diagonal, so the knob does what a knob should — a decade in buys a decade out — and neither of them sits on it. The factor is how much of the matrix's mass lives off the diagonal, which is a property of the kernel; it is not the number of blocks, which is 112 and would put the curves on the other side of the diagonal.10⁻¹³10⁻¹⁰10⁻⁷10⁻⁴10⁻¹10⁻¹⁵10⁻¹²10⁻⁹10⁻⁶10⁻³ε asked for, per block‖A − A_H‖ ⁄ ‖A‖ obtainedasked = obtainedlog r1 ⁄ rtwo numbers, not one1/r, obtained at 10⁻⁸3.4·10⁻¹⁰log r, obtained at 10⁻⁸7.8·10⁻¹⁰1/r, obtained ⁄ asked0.034log r, obtained ⁄ asked0.078blocks in the partition112the tolerance is per blockand the error is per matrix
Fig. 25 A tolerance asked for against one obtained.
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 next essay of this field.

Where this goes next

23456710⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²poles in the approximantresiduals and errorforward error‖T(λ)x‖‖T̃(λ)x‖one extra evaluationagainst the approximant1.5·10⁻¹²against the problem asked1.9·10⁻⁶forward error4.8·10⁻⁵‖g − r‖ there8.5·10⁻⁵the free residual is flatand the answer is not The problem the solver was actually given A linearisation is exact — it has the polynomial's eigenvalues, with their multiplicities, and the whole loss is arithmetic. A nonlinear eigenvalue problem does not offer that. Every algorithm replaces the function first, and the term that replacement contributes is committed before any number is rounded and appears in no residual. 110¹10⁻⁷10⁻⁵10⁻³distance from the branch point, λ + c|g − r|the eigenvaluesdegree 98 poleschosen, not computedreach0.06left end, from the cut0.2rational, worst8.5·10⁻⁴polynomial, worst0.0068linearisation size, both54committed before the solveand invisible to it An error committed before the arithmetic Before a nonlinear eigenvalue problem is solved, somebody says where they think the eigenvalues are. That sentence sets the accuracy of everything that follows by five orders, costs nothing to say, and cannot be revised once the approximation built on it is in hand. 123456024681012moments Krank of the block Hankel12 eigenvalues insidethe ceiling K·ℓthe ranka ceiling with a knob on itprobes2eigenvalues inside12rank at K = 12rank at K = 612solves, at every K512one ceiling per probeand K of them per moment A ceiling with a knob on it A contour method returns at most as many eigenvalues as its probe block has columns, and the object that comes back does not distinguish that from having found everything. One line of the derivation multiplies the ceiling by a number the caller chooses, and it costs no extra solves at all.

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 interpolationApproximation before linearisationBarycentric formBranch pointData driven realisationGreedy algorithmLoewner matrixRational approximationSupport points