An error committed before the arithmetic

The problem the solver was actually given establishes that a nonlinear eigenvalue problem is approximated before it is solved, and that the approximation contributes a term no residual contains. This essay is about the decision that sets that term, which is made before the approximation is built and is not the choice of approximant.

It is the target set: the region a caller believes the eigenvalues lie in.

The decision #

The function is γ√(λ + c), with a branch point at −c. A caller does not know where the eigenvalues are — that is what they are asking — so they name a region they believe contains them, and the approximation is fitted on that region.

How far the region reaches towards the branch point is the whole difficulty. Measured with ten poles:

Five orders, from a sentence — the same size of effect the scaling that buys ten orders measures for a change of units, and from a decision of the same kind: free to make, impossible to revise, and reported by nothing. The eigenvalues have not moved. The function has not changed. The same ten poles, the same fitting procedure, the same linearisation of the same size — and the answer that comes out of the whole chain is between nine digits and three.

Why the reach costs so much #

A rational function with poles on a branch cut approximates √ geometrically, and the rate is set by the conformal distance from the target set to the singularity. A target set standing well clear sees a function that is analytic on a neighbourhood of it and converges at a rate like ρ^k for a ρ comfortably below one. A target set touching the cut sees a function with unbounded derivatives at one end, and the achievable rate degrades towards algebraic.

That is the same structure as the Zolotarev rate in where to put the poles of a rational function: a geometric rate whose ratio depends logarithmically on how hard the region is. Here the hardness is the distance to a branch point rather than the condition number of a spectrum, and the dependence has the same shape.

Walking the reach shows both halves of that: how steeply the cost falls as the set stands clear, and how the two families compare while it does.

A factor of ten in the reach buys a factor of forty-five. Across reaches of 0.02, 0.04, 0.1, 0.15 and 0.2 the eight-pole error runs 3.31·10⁻³, 1.5·10⁻³, 3.54·10⁻⁴, 1.51·10⁻⁴ and 7.42·10⁻⁵, and the degree-nine polynomial’s runs 2.16·10⁻², 1.13·10⁻², 2.95·10⁻³, 1.27·10⁻³ and 6.09·10⁻⁴ — falling by 45 and 35 times respectively, for the same eight or nine degrees of freedom.

And the rational family’s advantage widens as the set stands clear: the ratio of the two errors runs 6.5, 7.5, 8.3, 8.4 and 8.2 across the sweep. That is the opposite of the intuition that poles matter most when the function is hardest — the poles are worth least at the hardest reach, because a set touching the branch cut degrades both families towards the same algebraic behaviour and the rational one has more to lose.

The consequence is worth stating as advice rather than as a rate. A region larger than the question needs is paid for. A caller who names a generous target set out of caution — better to be safe, the eigenvalues might be anywhere — has paid orders for the caution, and paid them in a term nothing downstream will report.

The quantity in that table is not the answer #

‖g − r‖ is how well the approximation fits the function. It is not what the caller asked for, which is the eigenvalues, and the two are measured together here for the first time.

The third column falls monotonically across five and a half orders, which is the finding above and is unchanged.

The fourth column has an interior minimum. It falls to 1.10·10⁻⁷ at a reach of thirty per cent and then climbs by four orders. The target set that gives the best approximation gives the worst answer, by sixteen thousand times.

Because containing the eigenvalues is not enough #

One detail of the table is worth reading before the mechanism, because it rules out the obvious alternative explanation. The eigenvalue error is not simply tracking the approximation error with a lag or a constant: at a reach of 0.15 the two are 1.53·10⁻⁵ and 1.36·10⁻⁶, a ratio of eleven, and at 0.70 they are 3.98·10⁻⁸ and 5.64·10⁻⁴, a ratio of one part in fourteen thousand the other way. The ratio between them crosses one and keeps going, so there is no calibration between the two columns and no way to infer either from the other.

The mechanism is the margin, and it is visible in the second column.

At the tightest reach the interval’s left endpoint is 2.9390 and the smallest eigenvalue is 2.9390. The eigenvalue sits exactly on the boundary of the region the approximation was fitted on — where a rational approximant is at its worst by construction, since the error equioscillates and the endpoints are where it peaks, and where the linearisation has nothing beyond it to work with. At a reach of thirty per cent the same eigenvalue has a margin of 2.3 and comes out four orders better.

So the rule the section above states is right about one end of the axis and wrong about the other, and the corrected version is one word longer. It is not name a tight region. It is:

Name a region that contains the eigenvalues with margin.

A region that reaches to the branch point costs orders in the approximation. A region whose boundary is an eigenvalue costs orders in the answer. Both are paid in a term nothing downstream reports, and the second is the one a caller taking this essay’s advice seriously would walk into.

Which is a caution about the objective rather than about the method #

That is worth separating from the finding, because it generalises past this approximant.

The approximation error is the quantity the fitting procedure minimises, and it is the quantity a figure can draw without solving anything. The eigenvalue error needs the whole chain — fit, linearise, solve — and needs an exact spectrum to compare against, which is available here only because the closed form is. So the cheap, drawable, optimisable quantity and the wanted one are different, they agree over most of the range, and they part company exactly where the cheap one is best.

This collection has now measured that shape three times in three fields — a column energy that is minimised while the convergence factor worsens, a step count whose residual falls while its error turns, and now a fitting error that falls while the eigenvalues get worse. In all three the optimised quantity is the tractable one and the wanted one is not, and in all three the disagreement is invisible until somebody measures both.

Where the poles go #

Given the target set, the poles are placed on the cut, and the rule is the one this collection keeps arriving at: geometric clustering towards the branch point.

Measured against evenly spaced poles, on a target set that stands off from the cut:

Clustering wins, by an amount that grows with the number of poles, which is what a rate difference looks like.

And where the rule stops #

On a target set that reaches to within a hundredth of the branch point, the same comparison:

Clustering wins nothing, and on three of the four settings loses slightly.

The reason is a real constraint rather than a numerical accident. The poles are placed geometrically from the target set’s near end down to its far end, and when the near end is a hundredth of the way to the cut, the tightest poles sit closer to the cut than any point where the function is being sampled. Two such poles produce basis functions 1/(λ − ξ₁) and 1/(λ − ξ₂) that are nearly identical over every sample point, the least-squares problem that fits the coefficients becomes ill conditioned, and the extra resolution is spent on a direction the fit cannot resolve.

So the rule is not cluster towards the singularity. It is cluster relative to the region of interest, and the two coincide only when the region stands off from the singularity. That is the same qualification the ADI shifts have — shifts outside the spectrum are wasted — arriving here in a form where it actually bites, because a target set can approach a branch point in a way a spectrum cannot approach its own end.

What is actually being fitted #

Worth one paragraph, because the fitting procedure is deliberately dull and that is a design decision.

The poles are fixed in advance, so r(λ) = c₀ + Σ αⱼ/(λ − ξⱼ) is linear in its unknown coefficients, and fitting it is one least-squares problem — a QR factorisation of a tall matrix and nothing else. No iteration, no greedy selection, no tolerance.

That matters because the approximation is supposed to be a fixed, known term in the error decomposition. An adaptive procedure that chose its own poles would have an accuracy story of its own, and the three-way separation in the previous essay would become a four-way one with a fuzzy boundary. Here ‖g − r‖ is a number that exists before anything else happens and stays fixed.

The sample points are geometric from the branch point rather than linear across the interval, for the same reason the poles are: a linear grid on an interval whose left end is a hundredth of the way to a singularity puts almost no points where the function is difficult. That is the same failure as evenly spaced ADI shifts, made one level down in the same construction, and it is easy to make twice.

The interaction nobody warns about #

Two decisions have now been named — the reach, and the pole placement — and they are not independent.

The reach sets how hard the approximation problem is. The pole placement is a response to that hardness, and its effectiveness is limited by the reach: on an easy target set clustering is worth an order and on a hard one it is worth nothing. So a caller who reaches further towards the cut both makes the problem harder and removes the tool that would have helped.

That is an unpleasant interaction and it is the practical reason to be conservative about the target set. The cost of a generous region is not the smooth degradation the first table suggests; it is a degradation that also disables the response.

The counterpart in the reduction field is milder. There the spectrum’s ends are a property of the problem rather than a choice, and geometric shifts work across the whole range — a wider spectrum costs a slower rate and nothing else. It is worth noticing which of the two situations a given construction is in — the distinction a condition number scaling cannot move draws between a difficulty that is a property and one that is a choice, because only one of them rewards care about the parameter.

What the target set is not #

Three readings are natural and all three would make the decision sound less consequential than it is.

It is not a tolerance. Nothing about the target set says how accurate the answer should be. It says where the answer is expected to live, and the accuracy follows from that as a consequence rather than as a request. A caller who wants more accuracy on the same region adds poles; a caller who names a smaller region gets more accuracy for free.

It is not a filter. Eigenvalues outside the target set are not excluded from the computation. They are simply approximated badly, and — as the eigenvalues that are answers to nothing shows — the linearisation returns a great many numbers outside the region, most of which are not eigenvalues of anything anybody asked about. Naming a region is a statement about where the approximation is good, not about what will come back.

It is not revisable. Once the approximant exists and the n(m+1) matrix has been factorised, the answer is what it is. There is no correction step that recovers the approximation error from the computed eigenvalues, in the way iterative refinement recovers rounding: refinement works because the residual against the original problem can be computed, and here it can — which suggests a repair and is the reason the next paragraph exists.

That last point deserves the qualification, because the repair does exist and is not part of this essay’s measurement. A computed λ can be polished by Newton on the original scalar equation, using the approximant’s answer as a start, and that removes the approximation error entirely at a cost of a few evaluations of g. What it does not do is find eigenvalues the approximation moved so far that Newton converges to a different one — which is the same failure a problem with infinitely many eigenvalues measures for Newton on a determinant, where five of ten runs land beside a different eigenvalue from the one they were started near.

What a caller can actually do #

Three things, and the first is by far the most valuable.

Estimate where the eigenvalues are before naming the target set. A handful of steps of any cheap method — a few Newton iterations on det T, a coarse contour count, an eigenvalue estimate of A alone — narrows the region enormously, and every order of narrowing is an order of approximation error. The estimate does not have to be good; it has to be a bracket.

Compute ‖g − r‖ once the approximant exists. It is a maximum of a scalar function over an interval and costs a few hundred evaluations. It is the term that will dominate the answer, and it is available before any eigenvalue is computed.

Check whether clustering is helping. Fit twice, once clustered and once evenly spaced, and compare the two maxima. It costs two QR factorisations of a few hundred rows. If the ratio is near one, the target set is touching the singularity and the real fix is to narrow it rather than to add poles.

All three are cheap relative to a single eigensolve of the n(m+1) matrix at the end, which is the recurring shape of the whole subject: the decisions that dominate the answer are the ones that cost nothing, and the expensive step is the one with the least influence over the result.

The two sample grids, and a defect that was made twice #

Building this family produced one defect worth recording, because it was the same mistake at two levels of the same construction and the second instance was made after the first had been found.

The approximation is fitted by least squares over sample points on the target set, and the first version sampled them linearly across the interval. On a target set whose left end is six per cent of the way to the branch point, that puts perhaps three of two hundred and forty points in the region where the function is difficult — so the fit is optimising a quantity dominated by the easy part of the interval, and the reported maximum, measured on a different grid, is far worse than the fit’s own objective.

The symptom was that adding poles stopped helping much earlier than it should. The repair is to sample geometrically from the branch point, which puts a fixed fraction of the points in each decade of distance from it.

The same error had already been made and fixed one level up, in the pole placement, where evenly spaced poles put every pole in the easy part of the region. Finding it a second time in the sample grid is the useful part of the story: the rule is about any set of points placed across a region whose difficulty is not uniform, and a construction can contain two such sets whose spacing is decided in two different places.

approximationError measures on a geometric grid too, so the objective and the report agree. Getting those two grids to disagree is a separate way to be wrong: a fit optimised on one grid and reported on another produces a number that is neither the fit’s accuracy nor the function’s, and it is right often enough not to look broken.

What the five orders are worth in context #

It is worth putting the size of this effect beside the others this field has measured, because the ordering of the four is the field’s practical content.

All four are made before any expensive arithmetic happens, and the eigensolve — which is the whole cost — comes last and has no influence on any of them. A reader arriving from the error field will recognise the shape: the algorithm is not where the accuracy is decided, and measuring the algorithm is measuring the part that was never in doubt.

The ordering is also the order to fix them in. Units first, because ten orders; then the region, because five; then the poles and the form. It is settled rather than argued, which is the same service six routes to one spectrum performed for the two decisions it compared.