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

The conditioning that rises with the ceiling

Higher moments multiply a contour method's ceiling by K and grade its block Hankel over ρ to the 2K, so the two knobs are the same knob. One division per quadrature point separates them, and the measurement of what it is worth grows from twenty to twenty thousand.
36 min read 30 figures Where the algebra stops being the arithmeticMeasured, not assertedIdentical algebra, different arithmetic

Worth reading first: A ceiling with a knob on it · Rank is a decision · An index that is a pair.

A ceiling with a knob on it multiplies a contour eigensolver’s ceiling by the number of moments used, for no extra solves. The claim is true and it is not the whole price.

The moments are integrals of zᵖ against a fixed function, so their sizes grow like |z|ᵖ on the contour. A circle of radius four has ‖A₃‖/‖A₀‖ of order 64 before anything about the problem is taken into account. The block Hankel of K moments is therefore graded over ρ²ᴷ, its condition number rises at exactly the rate the ceiling does, and the rank decision it asks for gets harder at the rate the answer gets better.

What the higher moments cost, and the one division that removes itThe condition number of the block Hankel against K, on a contour of radius 6, integrating in z and in the scaled variable (z − c)/ρ. Unscaled, ‖Aₚ‖ grows like ρᵖ, so the block Hankel of K moments is graded over ρ²ᴷ and its condition number rises at exactly the rate the ceiling does — from 1.28 to 2.26·10⁶ across six moments. Scaled, it runs 1.28 to 1750, a factor of 1290 better at the far end. The scaled variable is on the unit circle, so every moment is the same size; the cost is one division per quadrature point and the eigenvalues come back with c + ρ× applied.123456110²10⁴10⁶moments Kκ of the block Hankelintegrated in zin (z − c)/ρone division per pointradius6unscaled, K = 62.3·10⁶scaled, K = 61750the factor between1291the ceiling rises by Kand so did the conditioning
Fig. 1 The block Hankel’s conditioning against K, integrated in z and in the scaled variable. The two separate at the rate the ceiling rises.

The measurement

Condition number of the block Hankel over the directions that carry an eigenvalue, for the delay problem with two probes:

radius, K unscaled ÷ scaled
6, K = 3 19.6
6, K = 4 53.6
15, K = 5 19,440

Two knobs, and they multiply. The gap is ρ²ᴷ, so it is set by the radius as much as by the number of moments — a contour of radius fifteen and five moments is four orders apart from its own scaled version, where a radius of three is barely apart at all.

The division

The repair is a change of variable. Integrate in

w = (z − c) / ρ

rather than in z, where c is the contour’s centre and ρ its radius. On a circle w is on the unit circle, so |w|ᵖ = 1 and every moment is the same size. The eigenvalues come back with c + ρ× applied at the end.

It costs one complex division per quadrature point — or, since the parameterisation is known, one multiplication by a unit vector already computed for the quadrature weight. In this implementation it is the difference between w = z and w = {re: cos θ, im: sin θ}, which is a line.

That the fix is a scaling should be no surprise by now. It is the units the matrix is measured in arriving in a place where the “units” are a coordinate on a contour, and it has the same character everywhere it appears: exact in both directions, approximating nothing, and deciding how many digits of the answer survive. That is the property a condition number scaling cannot move is careful to distinguish from the cases where a rescaling changes nothing — here it changes everything, because the grading is a coordinate rather than a fact about the problem.

Why it matters more than a factor of twenty sounds

The number the conditioning spoils is not an eigenvalue. It is a rank, and a rank is a decision about a gap.

Rank is a decision is the essay about that, and the point transfers exactly. A block Hankel whose singular values are graded over ρ²ᴷ has directions that are genuinely present but scaled far below the largest — and a numerical rank test cuts relative to the largest. So the eigenvalues those directions carry are not merely computed inaccurately; they are not counted, and the method reports a smaller rank than the truth.

Which is the same failure the moments were introduced to fix. The ceiling was lifted from K·ℓ and replaced by a conditioning limit, and the second limit is worse than the first in one specific way: the first is a known integer that a caller can reason about, and the second is a number that depends on the contour and moves silently.

The scaled version has neither limit until much later, which is what makes the division worth a section rather than a footnote. It is the same reading an eigenvalue count that cannot be slightly wrong gives to an integer computed from an inertia: a quantity that is an integer either arrives right or arrives obviously wrong, and one that arrives quietly small is the worse failure.

What the two versions agree about

Both return the same rank while the conditioning allows it, and that agreement is what makes the comparison a measurement rather than a preference. At radius six and K = 3 both return six, and at K = 4 both return eight — the unscaled matrix is twenty and fifty times worse conditioned and still resolves the answer, because twenty and fifty are small numbers against 1/u.

The routes part where the argument says they should: when ρ²ᴷ approaches the reciprocal of the unit roundoff. At radius fifteen and five moments the unscaled gap is 1.9·10⁴, and pushing K further on that contour is where the unscaled version stops returning the right integer while the scaled one continues.

That is the same discipline every comparison of two algebraically identical routes on this site follows, from two Gram–Schmidts onwards: the interesting statement is not that one is better but where they stop agreeing, and whether the place was predicted before it was measured.

The Vandermonde underneath

It is worth saying where the grading comes from, because it is the same object that makes the moments work at all.

The pth moment picks out the same sum over eigenvalues inside the contour that A₀ does, with each term weighted by λⱼᵖ. So the matrix relating the stacked moments to the eigenvector directions is a Vandermonde in the λⱼ, and it has full rank whenever the eigenvalues are distinct — which is precisely why K moments lift the ceiling.

A Vandermonde in points spread over a circle of radius ρ has entries growing like ρᵖ, and is notoriously ill conditioned when its points span a wide range of magnitudes. So the object that supplies the extra rank is the object that spoils the conditioning, and there is no version of the construction in which one appears without the other.

What the change of variable does is make the Vandermonde’s points lie on the unit circle instead, where they all have modulus one and the matrix is as well conditioned as a Vandermonde on distinct points can be. It does not remove the Vandermonde; it puts its points somewhere they can be told apart.

That is worth carrying past this essay, because it is the general shape of the repair whenever a graded object appears: the grading is usually a coordinate rather than a fact, and finding the coordinate in which it disappears is cheaper than working harder in the coordinate that has it.

The measurement that had to be constructed

Getting the comparison to say anything took one repair worth recording, because the first version reported the graded matrix as the better conditioned of the two.

The obvious quantity is κ = σ₁/σ_r, where r is the rank the routine decided on. Run that on both versions and at a small radius the unscaled one comes back three times better — which is not a measurement, it is an artefact: the two routines made different rank decisions, so σ_r is the rth singular value of two matrices with two different r, and the comparison is between two different quantities.

The repair is to fix r from outside both routines. Here that is available, because the number of eigenvalues inside the contour is known from the closed form, so κ is taken over exactly the directions that carry one. Then the two numbers are about the same object and the comparison means what it looks like.

The general form of that mistake is worth naming, since it is easy to make in any comparison of two routines that each decide something. A quantity indexed by a routine’s own decision is not comparable across routines. A residual over the iterations a method chose to take, an error at the rank a routine selected, a cost at the tolerance an adaptive scheme settled on — each of them looks like a number and is a number plus a decision, and comparing two of them compares the decisions as much as the numbers.

What is left binding

Three things, after the division.

The number of quadrature points, which decides the accuracy of every moment equally and is governed by the distance from the contour to the nearest eigenvalue outside it. Nothing here touches that.

The problem’s own size, since the ceiling is K·min(n, ℓ) and the n side of the minimum is not a choice. A contour holding more than K·n eigenvalues is under-resolved whatever the conditioning.

The eigenvalues’ own spread inside the contour. The Vandermonde is well conditioned when its points are spread on the unit circle and poorly conditioned when they cluster, so a contour containing a tight cluster of eigenvalues is harder than one containing the same number spread out, at every K, in both versions. That is a property of the problem rather than of the coordinate, and it is the one limit no change of variable removes.

The last of those is worth noticing because it is the same statement about the same object as the gap decides the eigenvector in the spectral field: what makes a cluster hard is that the directions being separated are nearly the same direction, and no amount of arithmetic recovers a distinction the data does not contain.

What a caller should set

Two numbers, and the advice differs from the previous essay’s because the conditioning changes the ordering.

Keep the probe block small — two or three. A probe is a whole extra set of solves round the contour and a moment is a multiplication, so the arithmetic argues for spending on moments. What the conditioning adds is that a probe costs nothing in grading: the block Hankel of K moments and ℓ probes is graded over ρ²ᴷ regardless of ℓ. So probes buy ceiling at a high arithmetic price and no conditioning price, and moments buy it at no arithmetic price and a conditioning price.

Raise K until the rank stops rising, and no further. The diagnostic from the previous essay is still the right one, and the conditioning gives it a second reading: a rank that falls as K rises is the signature of grading rather than of a smaller true count, and on a scaled construction it should not happen at any K a caller would reach.

Scale the contour variable, always. It is a line, it costs a multiplication already being performed, and it removes the second failure entirely up to the point where the Vandermonde’s own clustering binds.

The ordering between the first two is what changes with the radius. On a small contour the moments are free in both currencies and there is no reason to spend on probes at all. On a large one the grading is what limits K, and past that limit the only remaining knob is the probe block — which is expensive, and is the point at which splitting the contour into several smaller ones becomes the better answer.

Splitting the contour, which is the other answer

The last paragraph named it and it deserves a section, because it is what a production implementation actually does and it changes which of the constraints binds.

A region containing forty eigenvalues can be handled as one contour with enough moments, or as four contours with ten each. The second costs four times the solves — the expensive resource — and every one of the four is a small contour with a small radius, so the grading is mild, the required K is small, and the Vandermonde’s points are spread over a region containing few eigenvalues rather than many.

The trade is therefore solves against conditioning, and it has the shape every trade in this collection has: the cheap resource runs out first. On a problem where each solve is a sparse factorisation, four contours is four times the dominant cost and is usually not worth it; on one where each solve is an iteration that converges quickly, it is nearly free and removes every conditioning question at once.

What makes the decision tractable is that both sides are measurable in advance. The number of solves is points × probes × contours and is known before anything runs. The grading is ρ²ᴷ and is known from the radius. So a caller can price both without computing an eigenvalue, which is the same arrangement the bound that is known in advance describes for a reduced model and which is rarer than it should be.

There is one genuine complication and it is worth flagging. Splitting a region requires knowing roughly where the eigenvalues are in order to split it sensibly, and an even split of a region whose eigenvalues are all in one corner produces one hard contour and three empty ones — which the empty-contour test from the previous essay at least reports rather than answering wrongly.

The two knobs, separated

The sentence this essay exists for: a construction whose accuracy knob and whose conditioning are the same quantity has one knob, and finding the coordinate in which they separate is worth more than any amount of care in turning it.

Here the knob is K, the ceiling is K·ℓ, the grading is ρ²ᴷ, and one line of change of variable makes the second independent of the first. After it the ceiling still rises with K and the conditioning does not, which is the whole difference between a parameter a caller can increase and a parameter they have to balance.

Where the argument came from

One note on provenance, since this measurement was not the one the work set out to make.

The higher moments were added to lift a ceiling that an earlier essay in this field had measured and left in place. The conditioning question arrived afterwards, from noticing that the moments’ sizes differ by ρᵖ — which is visible in the numbers as soon as they are printed, and is invisible in the derivation, where zᵖ is one symbol among several.

That is the ordinary way these things are found on this site: a quantity is printed for a different reason and its scale is wrong. The habit that produces it is printing the intermediate objects rather than only the answer, and the habit that would have missed it is checking that the answer is right — which it was, at every setting where the argument was first run.

The measurement then had to be constructed carefully, for the reason two sections above, and the first version of it said the opposite of the truth. Both halves of that are worth recording: the finding came from looking at a number nobody needed, and it came out backwards the first time it was measured.

At other settings

What the higher moments cost, and the one division that removes itThe condition number of the block Hankel against K, on a contour of radius 3, integrating in z and in the scaled variable (z − c)/ρ. Unscaled, ‖Aₚ‖ grows like ρᵖ, so the block Hankel of K moments is graded over ρ²ᴷ and its condition number rises at exactly the rate the ceiling does — from 5.67 to 3071 across six moments. Scaled, it runs 5.67 to 12.2, a factor of 251 better at the far end. The scaled variable is on the unit circle, so every moment is the same size; the cost is one division per quadrature point and the eigenvalues come back with c + ρ× applied.123456110²10⁴moments Kκ of the block Hankelintegrated in zin (z − c)/ρone division per pointradius3unscaled, K = 63071scaled, K = 612the factor between251the ceiling rises by Kand so did the conditioning
Fig. 2 A small contour, where the grading is mild and the two versions nearly agree.
What the higher moments cost, and the one division that removes itThe condition number of the block Hankel against K, on a contour of radius 5, integrating in z and in the scaled variable (z − c)/ρ. Unscaled, ‖Aₚ‖ grows like ρᵖ, so the block Hankel of K moments is graded over ρ²ᴷ and its condition number rises at exactly the rate the ceiling does — from 1.28 to 2.26·10⁶ across six moments. Scaled, it runs 1.28 to 695, a factor of 3250 better at the far end. The scaled variable is on the unit circle, so every moment is the same size; the cost is one division per quadrature point and the eigenvalues come back with c + ρ× applied.123456110²10⁴10⁶moments Kκ of the block Hankelintegrated in zin (z − c)/ρone division per pointradius5unscaled, K = 62.3·10⁶scaled, K = 6695the factor between3251the ceiling rises by Kand so did the conditioning
Fig. 3 Five, where the gap has become a factor of ten.
What the higher moments cost, and the one division that removes itThe condition number of the block Hankel against K, on a contour of radius 10, integrating in z and in the scaled variable (z − c)/ρ. Unscaled, ‖Aₚ‖ grows like ρᵖ, so the block Hankel of K moments is graded over ρ²ᴷ and its condition number rises at exactly the rate the ceiling does — from 1.28 to 2.26·10⁶ across six moments. Scaled, it runs 1.28 to 3.37·10⁴, a factor of 67.1 better at the far end. The scaled variable is on the unit circle, so every moment is the same size; the cost is one division per quadrature point and the eigenvalues come back with c + ρ× applied.123456110²10⁴10⁶moments Kκ of the block Hankelintegrated in zin (z − c)/ρone division per pointradius10unscaled, K = 62.3·10⁶scaled, K = 63.4·10⁴the factor between67the ceiling rises by Kand so did the conditioning
Fig. 4 Ten, and it grows with every added moment.
What the higher moments cost, and the one division that removes itThe condition number of the block Hankel against K, on a contour of radius 15, integrating in z and in the scaled variable (z − c)/ρ. Unscaled, ‖Aₚ‖ grows like ρᵖ, so the block Hankel of K moments is graded over ρ²ᴷ and its condition number rises at exactly the rate the ceiling does — from 1.22 to 7.37·10⁸ across six moments. Scaled, it runs 1.22 to 3952, a factor of 1.9·10⁵ better at the far end. The scaled variable is on the unit circle, so every moment is the same size; the cost is one division per quadrature point and the eigenvalues come back with c + ρ× applied.123456110²10⁴10⁶10⁸moments Kκ of the block Hankelintegrated in zin (z − c)/ρone division per pointradius15unscaled, K = 67.4·10⁸scaled, K = 63952the factor between1.9·10⁵the ceiling rises by Kand so did the conditioning
Fig. 5 And fifteen, where five moments are four orders apart.
K moments multiply a contour method's ceiling by K, for no extra solvesThe rank of the block Hankel of the first K contour moments, for a delay eigenvalue problem of size 4 with a probe block of 2 columns and a circle of radius 5 holding 12 eigenvalues. The dashed line is the ceiling K·min(n, ℓ) = 2K, and the rank sits exactly on it while it binds, then stops at 12 — the number genuinely inside. The solves are the expense and there are 512 of them whatever K is: a higher moment is one more multiplication at each quadrature point. So the ceiling the previous work on this problem measured and did not lift is lifted by a knob, and what it costs is not arithmetic.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
Fig. 6 The ceiling the moments lift, and the rank that sits on it.
K moments multiply a contour method's ceiling by K, for no extra solvesThe rank of the block Hankel of the first K contour moments, for a delay eigenvalue problem of size 4 with a probe block of 3 columns and a circle of radius 5 holding 12 eigenvalues. The dashed line is the ceiling K·min(n, ℓ) = 3K, and the rank sits exactly on it while it binds, then stops at 12 — the number genuinely inside. The solves are the expense and there are 768 of them whatever K is: a higher moment is one more multiplication at each quadrature point. So the ceiling the previous work on this problem measured and did not lift is lifted by a knob, and what it costs is not arithmetic.1234560369121518moments Krank of the block Hankel12 eigenvalues insidethe ceiling K·ℓthe ranka ceiling with a knob on itprobes3eigenvalues inside12rank at K = 13rank at K = 612solves, at every K768one ceiling per probeand K of them per moment
Fig. 7 With three probes, where the ceiling rises faster and the conditioning with it.
4 eigenvalues inside a circle of radius 4, and 2 moments to reach themA delay eigenvalue problem T(λ) = A − λI + e^{−λ}I, which has infinitely many eigenvalues and no linearisation, so the only finite question is what lies inside a given contour. The crosses are the Lambert-W closed form λ = μ + W_k(γe^{−μ}), one per eigenvalue of A per branch; the discs are what 2 contour moments against a probe block of 2 recover, agreeing to 2.16·10⁻¹². With one moment the method could have returned at most 2 of them and would have said nothing about the rest.-4-3-2-101234-4-3-2-101234real partimaginary partan infinite spectrum, finitely askedinside the contour4probes2moments used2values returned8worst against the closed form2.2·10⁻¹²infinitely many eigenvaluesand a question with an answer
Fig. 8 The eigenvalues the whole apparatus returns.
12 eigenvalues inside a circle of radius 6, and 6 moments to reach themA delay eigenvalue problem T(λ) = A − λI + e^{−λ}I, which has infinitely many eigenvalues and no linearisation, so the only finite question is what lies inside a given contour. The crosses are the Lambert-W closed form λ = μ + W_k(γe^{−μ}), one per eigenvalue of A per branch; the discs are what 6 contour moments against a probe block of 2 recover, agreeing to 4.22·10⁻¹³. With one moment the method could have returned at most 2 of them and would have said nothing about the rest.-7-5-3-11357-7-5-3-11357real partimaginary partan infinite spectrum, finitely askedinside the contour12probes2moments used6values returned24worst against the closed form4.2·10⁻¹³infinitely many eigenvaluesand a question with an answer
Fig. 9 And on a larger contour, where more moments are needed.
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. 10 The count that should be checked against the rank.
How many eigenvalues a contour method can return: the probe block is a ceiling and does not say soBeyn's method probes the contour with a random 4 × ℓ block V and reads the eigenvalues out of the rank of ∮T(z)⁻¹V dz. That rank is at most ℓ whatever is inside, so a probe block narrower than the number of eigenvalues present returns ℓ of them — with small residuals, and with nothing in the returned object to say that there were more. Measured with 2 eigenvalues inside |z| = 2: the rank is 1, 2, 2, 2, 2 at ℓ = 1, 2, 3, 4, 5, rising with the probe count and then stopping at 2. A caller who asked for two gets two. The counting integral is the check that says so, and it costs one more pass round the same contour.1 probe12 probes23 probes24 probes25 probes22 insidea ceiling nothing announceseigenvalues inside2rank at 1 probe1rank at 5 probes2the binding ceiling2the residuals are smalland half the answer is missing
Fig. 11 The ceiling as it stood before the moments.
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 820decisivetrue rank 410×10, built with 4 nonzero valuesrank is a decision
Fig. 12 The decision the conditioning spoils.
A randomised rank-k solve, 4 seeds a rankRelative error against the rank kept, on a logarithmic vertical axis, with a vertical bar at each rank spanning the seeds. The band is 1.68 wide at rank 8, where the method is at its worst, and 1.004 wide at rank 20, where it is at its best. The same computation on the same data returns a different answer each time.0481216202428321rank keptrelative errormedianthe answer movesspread at rank 81.7spread at rank 201best median error0.15widest where the method is worstand the bound does not say so
Fig. 13 Where that decision stops being one.
The worst relative error over the whole spectrum, against how many decades the matrix is graded over, n = 12Each point is one bidiagonal matrix and the worst relative error any of its singular values suffers, measured against a Sturm bisection in exact rationals. The route through BᵀB is as good as anything at 3.3 decades — 9.14·10⁻¹³ — and by 10 decades it is at 2.02, which is not an error in the answer, it is the answer. One-sided Jacobi and the zero-shift bidiagonal sweep are flat at about the unit roundoff across the whole range, which is the claim this site's own SVD has been making in a source comment since it was written and had never measured.0102030405010⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹10⁴10⁷decades of gradingworst relative errorthe answer is gonevia BᵀBone-sided Jacobizero-shift QRone axis, four routesBᵀB at the narrowest grading9.1·10⁻¹³and at the widest1Jacobi, worst over the sweep1.5·10⁻¹⁵zero shift, worst1.8·10⁻¹⁵the definition is not a methodand squaring buries what it squares
Fig. 14 Grading, measured directly, in the field that names it.
Two condition numbers of one 8×8 system, as its rows are put into different unitsFour curves against the spread of the row units, in decades. κ_∞ of the scaled matrix rises from 9.83 to 1.9·10¹⁰ while the componentwise condition number stays at 6.98 throughout — the same system, the same solution, and one of the two numbers is a fact about the units. Hilbert's two numbers are drawn flat beside them at 3.4·10¹⁰ and 1.2·10¹⁰: a matrix whose sensitivity no scaling repairs.0246810110²10⁴10⁶10⁸10¹⁰10¹²spread of the row units (decades)condition numberκ_∞(DA)cond(DA)Hilbert κ_∞Hilbert condone system, two numbersκ_∞ at no spread9.8κ_∞ at 10 decades1.9·10¹⁰cond, either end7Hilbert, equilibrated1.3·10¹⁰the solution is the same at every spreadand one of these curves knows it
Fig. 15 And the change of units that removes it there.
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. 16 The same repair in the polynomial field, worth ten orders.
Four pivot rules on [[ε, 1], [1, 1]] with ε = 10⁻¹², before and after its first row is scaled by 1/εEight bars of forward error. On the matrix as given, no pivoting loses the answer entirely and all three pivot rules return it to rounding. Multiply the first row by 1/ε — which does not change the solution — and partial pivoting makes no interchange at all and returns the same wrong answer that no pivoting does. Scaled partial pivoting and complete pivoting are unmoved.forward error, relative to a solution of exactly (1, 1)no pivoting · as given1.3·10⁻⁴ 0 interchangesno pivoting · rows scaled1.2·10⁻⁴ 0 interchangespartial · as given2.2·10⁻¹⁶ 1 interchangepartial · rows scaled1.2·10⁻⁴ 0 interchangesscaled partial · as given2.2·10⁻¹⁶ 1 interchangescaled partial · rows scaled2.2·10⁻¹⁶ 1 interchangecomplete · as given1.1·10⁻¹⁶ 1 interchangecomplete · rows scaled2.2·10⁻¹⁶ 1 interchangethe same problem twicepartial, as given2.2·10⁻¹⁶partial, rows scaled1.2·10⁻⁴its relative residual1.2·10⁻¹⁶complete, rows scaled2.2·10⁻¹⁶the two systems have the same solutionand one pivot rule cannot see it
Fig. 17 A decision that looks like an algorithm and is a scaling.
QᵀQ from classical Gram–Schmidt and from Householder on the 10×10 Hilbert matrixTwo eight-by-eight tables of QᵀQ. The upper one has ones on the diagonal and entries as large as one off it; the lower one is the identity to three decimal places everywhere.A = H10 · κ = 1.6·10¹³ · both factorisations reconstruct A to 6.4·10⁻¹⁷the diagonal is 1 in both — every column is a unit vector either way1.00000000000001.00000000000001.00000000000001.00000000000001.0000-0.001-0.001-0.001-0.001000001.000-0.053-0.042-0.035-0.0280000-0.001-0.0531.0001.0001.0001.0000000-0.001-0.0421.0001.0001.0001.0000000-0.001-0.0351.0001.0001.0001.0000000-0.001-0.0281.0001.0001.0001.000classical Gram–Schmidt1.00000000000001.00000000000001.00000000000001.00000000000001.00000000000001.00000000000001.00000000000001.00000000000001.00000000000001.000Householderclassical ‖QᵀQ − I‖3.5Householder ‖QᵀQ − I‖1.9·10⁻¹⁵largest off-diagonal 1 against 4·10⁻¹⁶length is not angle
Fig. 18 Two routes that agree and then stop, in the field that first drew the pair.
Hankel singular values, and the half of them a product cannot seeσₖ divided by σ₁, for a 20-state model of McMillan degree 12. The lower curve is the square-root route — Cholesky-like factors of the two Gramians and one SVD of RᵀS — and it descends to the level at which the Gramians were computed. The upper curve eigendecomposes the product PQ, agrees for the first 6 values and then stops, flattening at 9.67·10⁻¹⁰. The dashed line is σ₁√u = 2.32·10⁻⁹, which is where a route that squares the conditioning must stop: the eigenvalues of PQ are σ², an absolute error of uσ₁² on them is a relative error of √u on σ. κ(PQ) is 1.8·10³⁶ against 1.34·10¹⁸. Every σ below the line is a term in the error bound, so this is not an academic loss.1357911131510⁻¹⁸10⁻¹⁵10⁻¹²10⁻⁹10⁻⁶10⁻³1index kσₖ ÷ σ₁σ₁√utwo factorstheir productthe normal equations, againroutes agree to k =6product floors at9.7·10⁻¹⁰σ₁√u2.3·10⁻⁹κ(P)κ(Q)1.8·10³⁶√ of it1.3·10¹⁸do not form the productthe σ below the line are the bound
Fig. 19 Another pair that agrees and then stops, at a floor that was predicted.
Why the product route stops: two condition numbers and their productThe two Gramians of a 22-state model of McMillan degree 10, and what each route to the Hankel singular values is charged. κ(P) = 1.26·10¹⁸ and κ(Q) = 2.16·10¹⁸. The square-root route works with RᵀS, whose condition number is their geometric mean, 1.65·10¹⁸; the route that eigendecomposes PQ works at 2.72·10³⁶, which is past 1/u = 4.5·10¹⁵ — the point at which nothing small survives at all. The Gramians themselves are right: this one agrees with its closed form to 2.15·10⁻¹³. What is lost is lost in the last step, to a product nobody had to form.condition numbers, on a logarithmic scaleκ(P)1.26·10¹⁸κ(Q)2.16·10¹⁸√(κ(P)κ(Q)) — the SVD route1.65·10¹⁸κ(P)κ(Q) — the product route2.72·10³⁶1/u4.5·10¹⁵the mean, or the productand only one of them fits
Fig. 20 And the conditioning that explains it.
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. 21 A Vandermonde whose conditioning decides whether a fit is determined.
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. 22 And the same basis deciding whether roots can be recovered from coefficients.
What a perturbation of 10⁻⁶ moves, against the gap it is applied acrossThree curves against the gap between two eigenvalues, both axes logarithmic. Two are flat; the third rises steeply as the gap closes and then levels off.10⁻⁹10⁻⁷10⁻⁵10⁻³10⁻⁸10⁻⁵10⁻²gap between the two eigenvalueshow far it movedthe eigenvectorsthe eigenvaluestheir plane‖E‖ / gapone perturbation, three answerseigenvalue shift, spread over the sweep1plane angle, spread over the sweep1eigenvector angle, spread1.6·10⁵the dashed line is Davis–Kahan's ‖E‖/gaptwo of the three never noticed
Fig. 23 A cluster that no arithmetic separates.
The six largest Ritz values, against a doubled eigenvalue at 10Distance from 10 for each of the six largest Ritz values, on a logarithmic vertical axis. The single-vector run of 12 steps has one value at 10 to 4.6·10⁻⁸ and its second is 1.730 away — it has found the eigenvalue once. The block of two has two values at 10, to 1.1·10⁻⁹ and 4.7·10⁻⁹.12345610⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹Ritz value, largest firstdistance from 10one vectora block of twoa space, not a ratecopies found, one vector1copies found, block of two2Krylov dimension, one vector12the second copy is not in the spaceat any number of steps
Fig. 24 And a block that has to be wide enough to see one.
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. 25 The infinite spectrum all of this is asked about.
σ_min(zI − A) over the complex plane, for a bidiagonal 6×6 matrix with every eigenvalue at 0.8 and 2 above the diagonalA square of the complex plane shaded by how small σ_min(zI − A) is. Every eigenvalue is at 0.8, marked by the crosshair; the 10⁻³ level reaches out to 1.36, well outside the unit circle drawn through the picture. A perturbation of the matrix at the level of double-precision rounding can put an eigenvalue anywhere in that region.darker is smaller: 10⁻¹⁰, 10⁻⁶, 10⁻³, 10⁻¹one spectrum, two matricesspectral radius0.8reach of the 10⁻³ level1.4eigenvalues, all at0.8the circle is |z| = 1and every eigenvalue is well inside it
Fig. 26 The resolvent being integrated.
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. 27 The other route in this field, which has a different sorting problem.
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. 28 And the residual that sorts it.
Four candidate singularity tests on three matrices, at n = 120A table of four quantities on three matrices. 0.1·I at n = 120 has a condition number of exactly 1 and a determinant of 10⁻¹²⁰; a matrix with κ = 10¹⁰ normalised to |det| = 1 has a determinant of one. The determinant and its nth root give the wrong verdict on both, the smallest singular value gives the wrong verdict on the normalised matrix because it carries the units, and only σ_min/σ_max is right on all three. A test is read as saying *singular* when its number is below 10⁻⁶; green cells agree with the true verdict and red ones do not.does this matrix look nearly singular?green: the test agrees with the truth · red: it does not · the bar under each number is its magnitude, over sixty-two decades|det A||det A|^(1/n)σ_min1/κ = σ_min/σ_max0.1·I at n = 120perfectly conditioned10⁻¹²⁰0.10.11κ = 10¹⁰, |det| = 1nearly singular1110·10⁻⁶10·10⁻¹¹Hilbert at n = 8nearly singular2.7·10⁻³³8.5·10⁻⁵1.1·10⁻¹⁰6.6·10⁻¹¹the two counterexamplesκ of the scaled identity1its determinant10⁻¹²⁰κ of the normalised matrix10¹⁰its determinant1det(cA) = cⁿ det(A)so a determinant carries the units n times over
Fig. 29 Deciding a number is zero, which an empty contour needs.
What the higher moments cost, and the one division that removes itThe condition number of the block Hankel against K, on a contour of radius 6, integrating in z and in the scaled variable (z − c)/ρ. Unscaled, ‖Aₚ‖ grows like ρᵖ, so the block Hankel of K moments is graded over ρ²ᴷ and its condition number rises at exactly the rate the ceiling does — from 1.28 to 2.26·10⁶ across six moments. Scaled, it runs 1.28 to 1750, a factor of 1290 better at the far end. The scaled variable is on the unit circle, so every moment is the same size; the cost is one division per quadrature point and the eigenvalues come back with c + ρ× applied.123456110²10⁴10⁶moments Kκ of the block Hankelintegrated in zin (z − c)/ρone division per pointradius6unscaled, K = 62.3·10⁶scaled, K = 61750the factor between1291the ceiling rises by Kand so did the conditioning
Fig. 30 The hero again, beside the ceiling it prices.

Where this goes next

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. 1234567891010⁻¹⁸10⁻¹⁵10⁻¹²10⁻⁹10⁻⁶10⁻³1indexsingular valuecutoff, σ₁ · 10⁻¹⁰numerical rank 10gap 8.2·10⁶an opiniontrue rank 410×10, built with 4 nonzero valuesrank is a decision Rank is a decision A floating-point matrix does not have a rank. It has a spectrum of singular values, and somewhere in that spectrum is a place where the values stop being signal and start being noise. Deciding where is a judgement, and the evidence for it is a gap. 0246810110²10⁴10⁶10⁸10¹⁰10¹²spread of the row units (decades)condition numberκ_∞(DA)cond(DA)Hilbert κ_∞Hilbert condone system, two numbersκ_∞ at no spread9.8κ_∞ at 10 decades1.9·10¹⁰cond, either end7Hilbert, equilibrated1.3·10¹⁰the solution is the same at every spreadand one of these curves knows it The units the matrix is measured in One linear system, written twice. The rows of the second are the rows of the first in different units, the solution is identical to the last bit, and the condition number has moved by eight orders of magnitude. One of those two numbers is a fact about the problem and the other is a fact about the notation.

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.

Condition numberContour eigensolverContour integralMomentNonlinear eigenvalue problemNumerical rankQuadratureRankScaling