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

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.
36 min read 31 figures Where the algebra stops being the arithmeticMeasured, not assertedWhose fault is it

Worth reading first: The problem the solver was actually given · Why a Gramian can be truncated at all · A problem with infinitely many eigenvalues.

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 error that happens before any arithmetic: |g − r| on the target setγ√(λ + c) minus its approximant, plotted against the distance from the branch point at −0.4. The target set runs from -0.2 to 8.99, and its left end is 0.2 from the branch point — a choice made before anything is computed. The rational approximant with 8 poles clustered on the cut reaches 8.46·10⁻⁴; the polynomial of degree 9, which linearises to a matrix of the same size, reaches 0.00676. Both errors are committed at this stage. No residual computed later contains either of them, because a solver's residual is taken against whichever of these two functions it was handed.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
Fig. 1 The approximation error over the target set, for two approximants of the same linearisation size. The slider moves the left end towards the branch point.

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:

target set reaches to ‖g − r‖
the smallest eigenvalue 6.8·10⁻⁹
30% of the way to the cut 1.4·10⁻⁶
6% of the way 1.3·10⁻⁴
1% of the way 1.6·10⁻³

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.

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.

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:

poles 4 6 8 10
even ÷ clustered 2.2 3.5 5.6 9.0

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:

poles 4 6 8 10
even ÷ clustered 0.72 0.50 0.63 0.92

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.

decision worth
the units the polynomial is written in ten orders
the target set for the approximation five orders
the placement of the poles one order, where there is room
the choice of linearisation a factor of forty

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.

At other settings

The error that happens before any arithmetic: |g − r| on the target setγ√(λ + c) minus its approximant, plotted against the distance from the branch point at −0.4. The target set runs from 2.94 to 8.99, and its left end is 3.34 from the branch point — a choice made before anything is computed. The rational approximant with 8 poles clustered on the cut reaches 4.3·10⁻⁸; the polynomial of degree 9, which linearises to a matrix of the same size, reaches 1.91·10⁻⁷. Both errors are committed at this stage. No residual computed later contains either of them, because a solver's residual is taken against whichever of these two functions it was handed.10¹10⁻¹²10⁻¹⁰10⁻⁸10⁻⁶distance from the branch point, λ + c|g − r|the eigenvaluesdegree 98 poleschosen, not computedreach1left end, from the cut3.3rational, worst4.3·10⁻⁸polynomial, worst1.9·10⁻⁷linearisation size, both54committed before the solveand invisible to it
Fig. 2 A target set that stops at the smallest eigenvalue, where both approximants reach nine digits.
The error that happens before any arithmetic: |g − r| on the target setγ√(λ + c) minus its approximant, plotted against the distance from the branch point at −0.4. The target set runs from 0.602 to 8.99, and its left end is 1 from the branch point — a choice made before anything is computed. The rational approximant with 8 poles clustered on the cut reaches 2.22·10⁻⁵; the polynomial of degree 9, which linearises to a matrix of the same size, reaches 1.71·10⁻⁴. Both errors are committed at this stage. No residual computed later contains either of them, because a solver's residual is taken against whichever of these two functions it was handed.110¹10⁻⁸10⁻⁶10⁻⁴distance from the branch point, λ + c|g − r|the eigenvaluesdegree 98 poleschosen, not computedreach0.3left end, from the cut1rational, worst2.2·10⁻⁵polynomial, worst1.7·10⁻⁴linearisation size, both54committed before the solveand invisible to it
Fig. 3 Thirty per cent of the way to the cut.
The error that happens before any arithmetic: |g − r| on the target setγ√(λ + c) minus its approximant, plotted against the distance from the branch point at −0.4. The target set runs from -0.333 to 8.99, and its left end is 0.0668 from the branch point — a choice made before anything is computed. The rational approximant with 8 poles clustered on the cut reaches 0.00331; the polynomial of degree 9, which linearises to a matrix of the same size, reaches 0.0216. Both errors are committed at this stage. No residual computed later contains either of them, because a solver's residual is taken against whichever of these two functions it was handed.10⁻¹110¹10⁻⁷10⁻⁵10⁻³10⁻¹distance from the branch point, λ + c|g − r|the eigenvaluesdegree 98 poleschosen, not computedreach0.02left end, from the cut0.067rational, worst0.0033polynomial, worst0.022linearisation size, both54committed before the solveand invisible to it
Fig. 4 Two per cent, where the clustering has run out of room.
The error that happens before any arithmetic: |g − r| on the target setγ√(λ + c) minus its approximant, plotted against the distance from the branch point at −0.4. The target set runs from -0.373 to 8.99, and its left end is 0.0267 from the branch point — a choice made before anything is computed. The rational approximant with 8 poles clustered on the cut reaches 0.00729; the polynomial of degree 9, which linearises to a matrix of the same size, reaches 0.0377. Both errors are committed at this stage. No residual computed later contains either of them, because a solver's residual is taken against whichever of these two functions it was handed.10⁻¹110¹10⁻⁶10⁻⁴10⁻²distance from the branch point, λ + c|g − r|the eigenvaluesdegree 98 poleschosen, not computedreach0.008left end, from the cut0.027rational, worst0.0073polynomial, worst0.038linearisation size, both54committed before the solveand invisible to it
Fig. 5 And under one per cent, where the fit is fighting its own basis.
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. 6 What the term costs downstream, against the residuals that cannot see 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 1: ‖T̃(λ)x‖ — the residual against the problem the eigensolver was handed — stays at 1.8·10⁻¹⁵ to 2·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 19.5 to 21.2 times larger. The free number rises by 1.14·10⁷ across the sweep while the answer improves by 2.79·10⁴. 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 approximant2·10⁻⁸against the problem asked2.5·10⁻⁸forward error5.4·10⁻⁷‖g − r‖ there3.4·10⁻⁷the free residual is flatand the answer is not
Fig. 7 And on the easy target set, where it stops being the dominant term.
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. 8 The two approximants at equal linearisation size.
36 eigenvalues come back, 6 of them mean somethingThe real eigenvalues of the n(m+1) = 36 that a rational approximant with 5 poles produces, on the real line, with the branch point at −0.4 marked. The 6 discs to the right are the answers, agreeing with the closed form to 7.74·10⁻⁴. The 6 crosses to the left are exact eigenvalues of the approximant lying where γ√(λ + c) is not a real number at all — so the residual that says nothing about accuracy is decisive as a filter: against T̃ all 12 look perfect, and against T the spurious ones do not have a residual. The remaining 24 are complex and are discarded the same way.-10123456789-101λbranch point at −0.4the closed formwhich of these is an answereigenvalues returned36wanted6past the branch point6complex24worst against the closed form7.7·10⁻⁴all of them exactfor a problem nobody asked
Fig. 9 What the linearisation returns besides the answer.
The rational function 6 ADI shifts amount to, chosen two ways|∏(λ − pⱼ)/(λ + pⱼ)| over the spectrum of −A for a 30-point discretisation, whose ends are 9.86 and 3834 and whose condition number is 388.8. The lower curve uses 6 geometrically spaced shifts and touches zero at each of them, equioscillating between with a maximum of 0.1187. The upper curve uses 6 equally spaced ones: every shift sits in the top decade of the spectrum, the small end is never approached, and the maximum is 0.8922 — a factor of 7.52 worse, which the iteration pays squared. The two runs are the same algorithm, the same number of solves, and one different line where the shifts are chosen.10¹10²10³10⁻⁹10⁻⁷10⁻⁵10⁻³10⁻¹λ, over the spectrum of −A|r(λ)|geometric maxequally spacedgeometricone line of codeshifts6κ of the spectrum389geometric max0.12equally spaced max0.89the factor between7.5the same k solvesand one choice of where
Fig. 10 The same pole-placement question in the reduction field.
The rational function 12 ADI shifts amount to, chosen two ways|∏(λ − pⱼ)/(λ + pⱼ)| over the spectrum of −A for a 30-point discretisation, whose ends are 9.86 and 3834 and whose condition number is 388.8. The lower curve uses 12 geometrically spaced shifts and touches zero at each of them, equioscillating between with a maximum of 0.009967. The upper curve uses 12 equally spaced ones: every shift sits in the top decade of the spectrum, the small end is never approached, and the maximum is 0.766 — a factor of 76.9 worse, which the iteration pays squared. The two runs are the same algorithm, the same number of solves, and one different line where the shifts are chosen.10¹10²10³10⁻⁹10⁻⁷10⁻⁵10⁻³10⁻¹λ, over the spectrum of −A|r(λ)|geometric maxequally spacedgeometricone line of codeshifts12κ of the spectrum389geometric max0.01equally spaced max0.77the factor between77the same k solvesand one choice of where
Fig. 11 And with twice the poles, where the gap has widened.
Why a Gramian can be truncated at all: a rational approximation problemλₖ₊₁(P)/λ₁(P) for the controllability Gramian of a 30-state model, against Z_k², where Z_k is the smallest a degree-(k,k) rational function with poles on one side can be made on the other. The Gramian's spectrum falls from 1 to 2.48·10⁻⁸ in 10 steps, and the bound falls with it. Nothing about the model enters the bound except the two ends of the spectrum of −A, 9.86 and 3834, whose ratio is 388.8 — so the decay that makes model reduction possible is predictable from two numbers before the Gramian exists. That is the fact this whole field rests on and it is a fact about rational approximation.1234567891010⁻¹⁹10⁻¹⁶10⁻¹³10⁻¹⁰10⁻⁷10⁻⁴10⁻¹kλₖ₊₁ ÷ λ₁, and the boundZ_k²the Gramianpredicted from two numbersstates30κ of the spectrum389λ11 ÷ λ₁2.5·10⁻⁸the bound there5.2·10⁻⁴the cliff everything rests onand the reason for it
Fig. 12 Where that question came from there.
The ADI error is the rational function, not merely bounded by it‖P − ZZᵀ‖/‖P‖ after k steps of low-rank ADI on a 30-state model, against max|r|² for the same shifts, with P taken from its closed form rather than from another solve. The two curves are parallel: their ratio is 0.756 at k = 2 and 0.746 at k = 10, a spread of 1.0340 while the error itself falls from 0.289 to 3.87·10⁻⁴. A bound that is a constant factor from its quantity over four orders is an identity with a projection in front of it — the constant is how much of P lies where |r| is largest, which depends on the system and not on k.234567891010⁻⁴10⁻²1ADI steps krelative error, and max|r|²max|r|²the measured errora bound that is an identitysteps drawn9ratio, worst0.77spread of the ratio1error at k = 103.9·10⁻⁴its residual1.1·10⁻⁴four orders of errorand one constant
Fig. 13 And what the answer is worth on the iteration.
What two lines of scaling are worth: the same quadratic in nine systems of unitsThe forward error against the closed form for an overdamped chain of 8 masses, before and after Fan–Lin–Van Dooren scaling. The unscaled curve runs 7.49·10⁻¹⁴ to 0.001258 — every digit gone by the far end — and the scaled one runs 1.26·10⁻¹³ to 8.23·10⁻¹⁴, flat to within a factor of 2.01. The scaling is γ = √(‖K‖/‖M‖) and δ = 2/(‖K‖ + γ‖C‖), computed from three norms and nothing else, and the map back is λ = γμ with no rounding in the statement. Flatness is the half that matters: after scaling every stop of the sweep IS the same problem, so there is nothing left for the change of units to do.0246810⁻¹⁶10⁻¹³10⁻¹⁰10⁻⁷10⁻⁴10⁻¹log₁₀ γ, the change of unitsforward erroras writtenafter scalingone change of variableunscaled, worst0.0013scaled, worst1.7·10⁻¹³orders recovered10scaled coefficient spread4.5the answer was never the problemthe units were
Fig. 14 A free decision worth ten orders, in this field.
Where binary32 stops being able to hold the linearisation, and where the scaled problem does notForming the first companion linearisation needs γ²M to be representable. The upper row is whether every entry of (γ²M, γC, K) is finite in binary32, at twelve changes of units; the lower row is the same question after Fan–Lin–Van Dooren scaling. The boundary is √(largest finite value ÷ largest entry of M) = 1.8447·10¹⁹, computed from the format's parameters and nothing measured — and the sweep agrees with it at every stop: the last change of units that survives is 10¹⁹ and the first that does not is 10²⁰. The scaled row never fails, because δγ²M has norm about one whatever γ was. In a narrow format the scaling is not an accuracy device: it is what makes the problem exist.02468101201change of units, by exponent10⁰10¹10²10³10⁴10⁶10¹⁰10¹⁶10¹⁹10²⁰10⁴⁰10¹⁵⁰as writtenafter scalinga range questionlargest finite value3.4·10³⁸predicted boundary γ1.8·10¹⁹last γ that forms10¹⁹stops where scaling fails0not a poor answerno answer at all
Fig. 15 Where the same decision stops being about accuracy.
Fitting the same degree-10 polynomial in two basesOn the left, the data and two fitted curves that lie on top of each other. On the right, the condition numbers of the two design matrices, ten orders of magnitude apart.00.10.20.30.40.50.60.70.80.910.20.40.60.81xyboth fits, drawn on top of one anothercondition number of the design matrixmonomial2.1·10⁷Chebyshev2.3largest fitted coefficientmonomial176Chebyshev0.51rms residual: 4.2·10⁻⁶ and 3.1·10⁻⁷ — the data is fitted either way.30 points, degree 10, single precisionthe basis is part of the problem
Fig. 16 A basis whose conditioning decides whether a fit exists.
How far the coefficients can move without changing the fit, degree 10Relative increase in the residual against relative change in the coefficients, along the least determined direction. The residual does not move measurably until the coefficients have changed by more than a factor of one.10⁻⁴10⁻³10⁻²10⁻¹110¹10²10³10⁻¹⁶10⁻¹³10⁻¹⁰10⁻⁷10⁻⁴10⁻¹relative change in the coefficients, along the worst directionrelative increase in the residualcoefficients doubled245% change, fit unmoved in the sixth digit308×: the third digit movesκ(A) = 2.2·10⁷. Exact arithmetic would pick one point on this floor. It would not raise it.24 points, degree 10, monomial basisthe data leaves them free
Fig. 17 And a fit with many equally good answers.
The roots of ∏(x − i) up to degree 10, computed from the coefficients that define themA polynomial whose roots are the integers 1 to n, expanded exactly in BigInt and then handed to the route every library takes: the eigenvalues of the companion matrix. The measured relative error runs 1.711·10⁻¹³ at degree 6 to 1.127·10⁻¹⁰ at degree 10 — the third digit of an answer whose truth is an integer. The other curve is the root condition number times the unit roundoff, which is a prediction made from the coefficients alone, and it tracks the measurement to within two orders at every degree. Every coefficient here is still exactly representable, which is what makes the failure a conditioning one rather than a storage one.5791110⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²degreerelative error of the worst rootmeasuredκ × ua polynomial is its coefficientsworst root, measured1.1·10⁻¹⁰predicted, κ × u1.4·10⁻¹⁰root condition number6.1·10⁵largest coefficient1.3·10⁷the roots are integersand the coefficients are not the roots
Fig. 18 Coefficients that stop determining roots.
A spectrum with no last eigenvalue: 12 of the infinitely many solutions of a delay problemT(λ) = A − λI + 0.5e^−λI for the model matrix at n = 4. The problem decouples in A's eigenvector basis, so each of its 4 eigenvalues μ contributes the scalar equation μ − λ + γe^−λ = 0, whose solutions are λ = μ + W_k(γe^−μ) — one for every branch of the Lambert W function, and therefore countably infinitely many, running off to the left along two curves. The marks are the closed form, each checked by substitution to 1.2·10⁻¹⁵. There is no computation that returns this set: the only well-posed question is which eigenvalues lie inside a region, and the circle of radius 4 contains 4 of them.-3-2-101234-5-3-1135real partimaginary part4 insidea countable spectruminside the contour4drawn12existing∞worst branch residual1.2·10⁻¹⁵there is no last eigenvalueso the question has to change
Fig. 19 The other nonlinear problem in this field, whose function has no branch point.
Counting the eigenvalues inside a circle: an integral that is an integer once it has convergedThe number of eigenvalues of the delay problem inside |z| = 1.5, computed as (1/2πi)∮tr(T(z)⁻¹T′(z))dz by the trapezoidal rule. The true count is 1, from the closed form. The error runs 0.473, 0.39, 0.289, 0.142, 0.0268, 7.58·10⁻⁴ at 4 to 128 points — a straight line on this axis is a constant factor a doubling, and this falls faster than that, which is what exponential convergence looks like when the integrand is analytic on the contour. The answer is then ROUNDED, and 0 of the 6 quadratures round to the wrong integer while looking exactly as confident as the ones that do not. The only evidence available is the distance from the nearest integer, which is the quantity plotted.11.31.61.92.210⁻¹⁶10⁻¹³10⁻¹⁰10⁻⁷10⁻⁴10⁻¹10²log₁₀ quadrature pointsdistance from the true counthalf an eigenvaluean integer, eventuallytrue count1at 4 points1at 128 points1finest error7.6·10⁻⁴the integral is an integerand a rounding hides how far it was
Fig. 20 And the question that can be asked of it without approximating anything.
One quadratic eigenvalue problem in 9 systems of units: what the solver reports and what the answer is worthAn overdamped chain of 6 masses, with λ replaced by γμ so that the coefficients become (γ²M, γC, K). That substitution is exact in both directions and divides the spectrum by γ exactly, so the closed form is still available and every error here is measured against it. The backward error of the eigenpair for the LINEARISED MATRIX — the residual a solver's own error analysis is about — is 7.1·10⁻¹⁶ at γ = 1 and 2.57·10⁻¹² at γ = 108 — it moves by a factor of 7022 while the other two move by 7.37·10¹¹. The backward error for the QUADRATIC, which is what the person who posed the problem is entitled to, grows by 7.03·10¹¹ across the same sweep, and the forward error follows it: 9.14·10⁻¹⁵ to 0.00674. Nothing went wrong with the solver at any stop.0246810⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹log₁₀ γ, the change of unitsrelative errorforward errorη, the quadraticη, the linearisationagainst a closed formη(linearisation), worst2.6·10⁻¹²η(quadratic), worst4.6·10⁻⁴forward error, worst0.0067coefficient spread4.2·10¹⁵the solver is right at every stopabout a problem nobody asked
Fig. 21 The backward error the approximation term sits in front of.
The condition number the problem has, and the one the solver's error analysis is written againstThe same eigenvalue of the same overdamped chain of 10 masses, in seven systems of units. Its condition number as an eigenvalue of the QUADRATIC — Tisseur's, with the three coefficient norms in the numerator and yᵀQ′(λ)x in the denominator — is 5.554 at γ = 1 and 5.554 at γ = 10⁶, a spread of 1 over six decades: it cannot move, because a change of units is not a change of problem. Its condition number as an eigenvalue of the LINEARISED MATRIX runs 55.22 to 8.685·10¹¹, a factor of 1.57·10¹⁰. The forward error follows the second one, and the first one is the honest description of the problem — so the substitution has manufactured an ill conditioning that belongs to the algorithm rather than to the question.012345610⁻¹10²10⁵10⁸10¹¹log₁₀ γ, the change of unitscondition numberthe linearisationthe quadraticone problem, two amplifiersκ(quadratic), first5.6κ(quadratic), last5.6κ(linearisation), last8.7·10¹¹how far the first moved1the problem is as well conditioned as everand the method is not
Fig. 22 And the condition number that does not multiply it.
Two routes to e^A as a 8×8 bidiagonal matrix walks up to a defective oneA_δ is upper bidiagonal with diagonal λ, λ+δ, …, so its eigenvalues are exact and its eigenvectors have a closed form. As δ falls the eigendecomposition route's relative error rises to 2.92·10⁶⁵ while scaling and squaring's falls to 4.09·10⁻¹², which is the distance ‖A_δ − A₀‖ = 7·10⁻¹² to the matrix whose exponential is known in closed form. The two agree at the right of the axis, where the matrix really is far from the limit; the divergence is κ(V), which reaches 3.3·10⁸².10⁻¹²10⁻¹⁰10⁻⁸10⁻⁶10⁻⁴10⁻²10⁻¹⁴10⁻¹10¹²10²⁵10³⁸10⁵¹10⁶⁴δ, the gap between consecutive eigenvaluesrelative error in e^Aan answer with no correct digitsV f(Λ) V⁻¹scaling and squaring‖A_δ − A₀‖exact eigenvalues throughoutκ(V) at the smallest δ3.3·10⁸²eigen route2.9·10⁶⁵scaling and squaring4.1·10⁻¹²distance to the limit7·10⁻¹²the eigenvalues are the diagonaland they are exact at every stop
Fig. 23 Approximating a function of a matrix, one field over.
Two filters on one sum, λ = 0.001The weight each term of the solution is given, against its index. Truncation is a step: one for the first 32 terms and zero after. Tikhonov is σ²/(σ² + λ²), which falls smoothly through the same place. The unregularised solution is the constant one, which is why it divides noise by a σ of 1.7·10⁻¹³.081624324048566400.250.50.751index kfilter factor fₖno regularisation: fₖ = 1truncationTikhonovthe same sum, three weightsTikhonov, relative error0.16truncation, relative error0.13no filter at all5.5·10⁸both filters are one expression with a different weightfₖ = 1 is the catastrophe
Fig. 24 A rational filter chosen for a different purpose.
Singular values and |uₖᵀb|, with and without 0.10% noiseThree curves against the index on a logarithmic vertical axis. The singular values fall exponentially to the level of rounding. With an exact right-hand side the coefficients fall faster and every term of the solution stays bounded. With noise they flatten at index 32, and from there on each term is noise divided by a σ of 6.8·10⁻⁴.081624324048566410⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹index kmagnitudethe floor: k = 32best truncation: k = 28σₖ|uₖᵀb| exact|uₖᵀb| with noisetwo different indicesthe crossing, from the data alone32the truncation that is actually best28relative error there0.11the exact coefficients never flattenthe noisy ones stop at ‖e‖/√n
Fig. 25 Where a problem stops determining its own answer.
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. 26 Choosing a parameter with no information, in the regularisation field.
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. 27 Points placed badly in a different construction.
A reduced model of order 4, and the 4 places it is exact|H(s) − Hᵣ(s)| ÷ |H(s)| along the real axis, for a rational-Krylov reduction of a 24-state model at the interpolation points 1.5, 5, 16, 50. At each of them the curve falls to 5.3·10⁻¹⁶ — the reduced function passes through the original, and its derivative does too, because the projection is two-sided. Between and beyond them it reaches 8.37·10⁻⁴, and there is no bound on it: the method buys 8 exact conditions for 4 solves and offers nothing anywhere else. That is the trade against balanced truncation, which asks for nothing and bounds everything at a cost of two Lyapunov solves.10⁻¹110¹10²10³10⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵s, on the real axis|H − Hᵣ| ÷ |H|exact where askedpoints4conditions bought8worst at a point5.3·10⁻¹⁶worst away from one8.4·10⁻⁴4 points, 8 conditionsand no bound in between
Fig. 28 And placed well.
Relative error of each singular value of a bidiagonal matrix graded over 11 decades, by four routesThe matrix is 8×8, bidiagonal, and every entry is a power of two times a small integer, so the matrix the float routines are handed and the matrix the exact routine reasons about are the same matrix. The reference is a Sturm bisection in BigInt rationals, which never rounds. σ_min is 1.934·10⁻¹¹. One-sided Jacobi holds every value to 7.1·10⁻¹⁶ relative and the zero-shift bidiagonal sweep to 4.41·10⁻¹⁶; the eigenvalues of BᵀB are wrong by 14.5 — not a rounding error but the whole number — and the damage is entirely at the small end.1234567810⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹singular value, largest firstrelative errorone-sided Jacobizero-shift QRshifted QReigenvalues of BᵀBagainst a rational bisectionσ_min, exactly1.9·10⁻¹¹worst, one-sided Jacobi7.1·10⁻¹⁶worst, zero-shift QR4.4·10⁻¹⁶worst, eigenvalues of BᵀB15a relative error is a ratioand the denominator is the answer
Fig. 29 A spectrum spanning ten decades, which is the reduction field’s version of a hard region.
Where a quadratic stops being hyperbolic, located by a Cholesky and by a sineThe overdamping margin min over modes of (βκ)² − 4κ, for a chain of 10 masses, against β. It reaches zero at β* = 1/sin(π/2(n+1)) = 7.02667418333, which is the closed form. Bisecting on a completely different question — does a Cholesky of −Q(μ) complete for some μ — gives 7.02667418333, agreeing to 14 digits. Neither route computes an eigenvalue. The marks below the axis are the largest imaginary part in the computed spectrum, which is zero to the rounding level above β* and not below it, so a third route agrees with the other two about where the boundary is.56.171117.342228.513349.6844510.855610⁻¹110¹10²stiffness damping βoverdamping marginβ* = 7.02667two routes to a boundaryclosed form β*7by certificate7difference1.9·10⁻¹³bisection steps44a factorisation that completesand a sine, agreeing to twelve digits
Fig. 30 A boundary located by bisection on a certificate, which needs no target set at all.
The error that happens before any arithmetic: |g − r| on the target setγ√(λ + c) minus its approximant, plotted against the distance from the branch point at −0.4. The target set runs from -0.2 to 8.99, and its left end is 0.2 from the branch point — a choice made before anything is computed. The rational approximant with 8 poles clustered on the cut reaches 8.46·10⁻⁴; the polynomial of degree 9, which linearises to a matrix of the same size, reaches 0.00676. Both errors are committed at this stage. No residual computed later contains either of them, because a solver's residual is taken against whichever of these two functions it was handed.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
Fig. 31 The hero again, beside the residuals that do not contain it.

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. 10¹10²10³10⁻⁹10⁻⁷10⁻⁵10⁻³10⁻¹λ, over the spectrum of −A|r(λ)|geometric maxequally spacedgeometricone line of codeshifts6κ of the spectrum389geometric max0.12equally spaced max0.89the factor between7.5the same k solvesand one choice of where Where to put the poles of a rational function Three times in one field the same question arrives from different directions — ADI shifts, rational approximation of a square root, the decay of a Gramian — and it has one answer. Cluster them geometrically towards wherever the function is difficult, and the alternative that looks reasonable costs orders. 0246810⁻¹⁶10⁻¹³10⁻¹⁰10⁻⁷10⁻⁴10⁻¹log₁₀ γ, the change of unitsforward erroras writtenafter scalingone change of variableunscaled, worst0.0013scaled, worst1.7·10⁻¹³orders recovered10scaled coefficient spread4.5the answer was never the problemthe units were The scaling that buys ten orders Two lines computed from three norms, a change of variable that is exact in both directions, and the whole of the loss the previous essay measured comes back — flat, at every stop, because after scaling every stop is the same problem.

What links here

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

Shares its objects with

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

Named objects

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

Approximation before linearisationBranch pointCondition numberLeast squaresNonlinear eigenvalue problemRational approximationShift selectionTarget setZolotarev number