The points the algorithm chose

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

The clustering and the advantage rise together as the target set approaches the singularity. Across reaches of 0.15, 0.05, 0.02 and 0.008 the clustering ratio runs 1.348, 1.45, 1.576 and 1.712 and the advantage over even spacing runs 566, 616, 747 and 812 times — while the absolute error gets steadily worse, 5.41·10⁻¹⁵, 8.28·10⁻¹⁴, 6.07·10⁻¹² and 1.54·10⁻¹⁰.

So the adaptive choice is worth most exactly where the problem is hardest, which is the opposite of how most of this collection’s advantages behave — the tree preconditioner, the Zolotarev bound and the admissibility rule all pay best on the easy cases. The reason is visible in the clustering column: as the set reaches towards the singularity there is more for the algorithm to decide, and an even spacing has no way to spend its points where the function is changing. At reach 0.15 the function is nearly benign and the algorithm’s freedom is worth 566×; at 0.008 it is worth 812×, and the degree needed has moved by one.

That the degree stays at 11 or 12 across a factor of twenty in the reach is worth noting beside it. The algorithm does not answer a harder problem by using more points — it answers it by putting the same twelve somewhere else.

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

Locating the poles of a barycentric rational in general requires a generalised eigenvalue problem of its own, which is why they are usually left unmeasured beside an approximation error. On a real interval they are cheaper than that: D(x) = Σ wⱼ/(x − zⱼ) changes sign between consecutive support points, so a bisection on those sign changes finds every real pole, and the residue at a pole p is N(p)/D′(p) — two more sums.

Poles inside the sampled interval, and their residues relative to the function’s scale:

With clean data there are no interior poles at all. The algorithm converges at degree 14 and every pole it places lies outside the target set — which is where a branch cut’s poles belong, and is the approximant doing exactly what it should.

With any noise the interior poles appear, and their residues are a fixed fraction of it: 7.9·10⁻³ of the noise, then 1.0·10⁻³, 9.1·10⁻³, 9.2·10⁻³, 6.7·10⁻³ — constant to within a factor of nine across four decades. So the degradation in the table above and the appearance of doublets are indeed the same phenomenon, and the measurement says how tightly: they switch on together at the same noise level, and one is a constant multiple of the other.

That also shows why a residue threshold is the wrong doublet test on its own. At noise 10⁻¹² the interior residues are at 10⁻¹⁵ and look like rounding; at 10⁻⁴ they are at 10⁻⁷ and do not. The same poles are equally spurious in both cases — they are interpolating noise — and what identifies them is that they are inside the set the function is being approximated on, not that their residues are small. A stopping rule keyed to the degree, which is what the table above measures, prevents both.

One more column in that table is worth reading, because it is the reason the degree matters more than the error does. The degree column is 14 in the clean case and 24 — the cap — in every noisy one. The algorithm is not slowly degrading as the noise rises; it stops converging at all the moment the noise exceeds the accuracy it could otherwise reach, and then runs to whatever limit it is given. A code that sets maxDegree to 40 and reads the error will get a larger linearisation and a spurious pole for every extra support point, at an approximation error that has stopped improving.

So the practical stopping rule is not “stop when the error stops falling”, which is a comparison of two small numbers on a noisy curve. It is “stop when the degree stops being able to buy error”, which the table makes unambiguous: the clean run halts on its own, and no noisy run ever does.

assertTheInteriorPolesArriveWithTheNoise measures the pole count and the residues, and requires the clean case to have none, every noisy case to have some, the degree to run to the cap in every noisy case, and the residues to be proportional to the noise across decades.

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.