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

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.
35 min read 30 figures Where the algebra stops being the arithmeticMeasured, not assertedExact ground truth

Worth reading first: A problem with infinitely many eigenvalues · Rank is a decision · A block nobody can call sparse.

A problem with infinitely many eigenvalues established three ceilings on Beyn’s contour method and lifted none of them. A probe block of ℓ columns returns at most ℓ eigenvalues; a problem of size n returns at most n; and the object that comes back distinguishes none of the three from having found everything. A caller who asks for two gets two, both genuine, with nothing to say there were twelve.

That last clause is what makes the ceilings dangerous rather than merely limiting. The repair is one line of the derivation, and it is worth stating exactly because it is the difference between a method with a ceiling and a method with a knob.

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. 1 The rank a contour method can return, against the number of moments used, with the ceiling drawn as a line it sits on.

The line

The first moment

A₀ = (1/2πi) ∮ T(z)⁻¹V dz

is n × ℓ for a probe block V of ℓ columns, so its rank is at most min(n, ℓ) however many eigenvalues are inside the contour. That is the ceiling, and it is a property of the shape of A₀ rather than of anything about the problem.

The higher moments

Aₚ = (1/2πi) ∮ zᵖ T(z)⁻¹V dz

cost one more complex multiplication each. The solves are the expense, there are points × probes of them, and they are already done at every quadrature point — so the higher moments are free in the only currency that matters.

Stack them into a block Hankel matrix of K blocks:

H₀ = [ A₀ A₁ … A_{K−1} ; A₁ A₂ … ; … ],

which is nK × ℓK. Its rank is the number of eigenvalues inside, up to K·min(n, ℓ). K moments multiply both ceilings by K, for no extra solves.

The measurement, and both halves matter

On a delay problem of size four with a probe block of two columns and a circle of radius six holding twelve eigenvalues:

moments K 1 2 3 4
ceiling K·ℓ 2 4 6 8
rank returned 2 4 6 8

The rank sits exactly on the ceiling at every K, which says the ceiling is binding — there is more inside the contour than the method can see, at every one of those settings.

Now three probes and a radius of five, holding twelve:

moments K 1 2 3 4 5
rank returned 3 6 9 12 12

It rises with K and then stops at twelve, which is the number genuinely inside.

Both halves are the finding. That it rises says the ceiling was binding rather than the count being right. That it stops says the count is real rather than the ceiling still binding. A rank that keeps rising with K is a rank that has not found everything yet — which is the diagnostic the original construction could not offer, and it is available for the cost of computing the rank at two values of K. It is the same diagnostic keeping the vectors, and losing the bound asks for in the Lanczos field: a quantity that keeps moving when it should have settled is the signal that a method has not finished rather than that it has.

What comes back, and against what

The eigenvalues themselves come from the two block Hankels: H₀ = V̂ΣŴ*, then B = V̂*H₁ŴΣ⁻¹, and the eigenvalues of B are the eigenvalues inside the contour.

B is complex. Its spectrum is read through the real embedding [[X, −Y], [Y, X]], whose eigenvalues are those of B together with their conjugates — the same device pseudo.js uses for σ_min and the same one this field uses to solve complex systems, which keeps a complex eigensolver out of the site entirely — the same economy the message and the word argues for in a different currency: one implementation used three times is one thing that can be wrong.

Checked against the Lambert-W closed form λ = μ + W_k(γe^{−μ}), one per eigenvalue of A per branch: 2.2·10⁻¹² worst — measured against a formula rather than against a better computation, which is the standing an answer that is known established and which a method with this many stages needs — for four eigenvalues inside a circle of radius four found with two probes and two moments. The whole apparatus — quadrature, moments, a block Hankel, a truncated SVD, a small complex eigenproblem read through a real embedding — measured against a formula and against nothing else.

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. 2 An infinite spectrum, a circle, and what the moments recover from inside it.

The two knobs, and why they are not interchangeable

The ceiling is K·ℓ, so K = 6 with one probe reaches the same ceiling as K = 2 with three. The costs are completely different.

A probe is a whole extra set of solves round the contour: 256 or 512 complex linear systems. A moment is one complex multiplication per quadrature point per probe, which is a few thousand flops.

So the picture is an argument for spending the budget on moments. What stops that from being the whole answer is the conditioning, which is the conditioning that rises with the ceiling and is the reason this essay stops at four rather than at forty.

There is also a genuine role for probes that moments cannot fill. The probe block has to have a component along every eigenvector inside the contour, and a random block of ℓ columns has one with probability one — but the size of that component decides how visible the corresponding singular value is. A single probe unlucky enough to be nearly orthogonal to one eigenvector makes that eigenvalue’s direction small in every moment, and adding moments does not help, because they all inherit the same V. Two probes make the bad case require two coincidences.

That is worth stating because it is the kind of thing a purely arithmetic comparison misses. The right answer is a small ℓ — two or three — and as many moments as the conditioning allows.

An empty contour, and a rank decision on rounding

Building this produced one defect worth recording, because it is the failure rank is a decision is about, arriving in a routine whose right answer is zero.

A contour with nothing inside it produces a moment matrix that is entirely rounding: its singular values here are 6.6·10⁻³⁴, against integrand values of 10⁻¹⁸. A rank test relative to the matrix’s own largest singular value looks at that and sees a perfectly good gap — the values are all of the same size and there is nothing to cut — so it returns full rank, confidently, on numbers that are the arithmetic error of an integral of nothing.

The fix is a scale that comes from outside the matrix. The routine tracks the size of what it integrated as it goes and cuts relative to the larger of that and the leading singular value. Then the empty contour returns zero, and the reason is stated rather than thresholded: ‖A₀‖ is 10⁻¹⁶ of the integrand’s size, which is the unit roundoff and not a rank.

The general lesson is the one that essay states and this is a clean instance of: a relative rank test has nothing to be relative to when the matrix is entirely noise, and every rank decision needs a scale that did not come from the matrix being decided about.

Why a moment lifts anything at all

The line is short and the reason it works is worth a paragraph, because “add more moments” sounds like it should not buy independent information.

T(z)⁻¹ has a simple pole at each eigenvalue λⱼ inside the contour, with residue a rank-one object built from the left and right eigenvectors. So the contour integral picks out a sum over the eigenvalues inside, and A₀ is a sum of k rank-one terms, each of which is (right eigenvector) times (left eigenvector applied to V). That sum has rank at most k — and at most ℓ, because V has ℓ columns and every term factors through it.

The pth moment picks out the same sum with each term multiplied by λⱼᵖ. Same k directions, different weights. Stack K of them and the coefficient matrix relating the stacked moments to the k eigenvector directions is a Vandermonde matrix in the λⱼ — which has full rank whenever the eigenvalues are distinct.

That is where the extra rank comes from. The directions were always there in the residues; one moment collapses them through an ℓ-column filter, and K moments look at the same collapse through K different weightings whose combination is invertible. Nothing new is computed about the problem — what changes is how much of what was already integrated survives the projection.

The Vandermonde is also the reason for the price. A Vandermonde in points spread over a circle of radius ρ has entries growing like ρᵖ, so its conditioning degrades with K exactly as fast as the rank rises, which is the subject of the next essay and is not an accident of implementation.

What the count was for

It is worth connecting back, because this field already had a way of counting eigenvalues inside a contour and it did not need any of this.

Counting what is inside a circle computes (1/2πi)∮ tr(T(z)⁻¹T′(z)) dz, an integral that is an integer, converging exponentially and reaching 1.4·10⁻¹³ at 128 quadrature points. It answers how many and says nothing about which, and it has no ceiling of any kind.

So the two constructions divide the work. The trace integral gives the count for one pass round the contour and is the right thing to run first. The moments give the eigenvalues, need enough of them to clear the ceiling, and — now that the rank rises and stops — can be checked against the count the trace integral already produced.

That check is the thing the original method most lacked. A rank of six with a count of twelve says immediately that K is too small; a rank of twelve with a count of twelve says the region has been resolved. Two integrals round the same contour, one integrand each, and the pair answers a question neither answers alone.

What the method now needs from its caller

Three inputs, and the change of emphasis is worth setting out because the moments move where the difficulty lives.

A contour. Still the hardest input, and the one nothing in this essay helps with. It has to enclose the eigenvalues wanted and to stay away from the ones not wanted, and the convergence of the quadrature is governed by the distance from the contour to the nearest eigenvalue outside it — which a problem with infinitely many eigenvalues measures: a circle passing near an eigenvalue at 1.587 reaches only 7.6·10⁻⁴ at 128 points where one comfortably clear of everything reaches 1.4·10⁻¹³.

A number of quadrature points. Decided by the same distance, and checkable: run the count at two resolutions and see whether the integer moves.

A number of moments. New, and the easiest of the three, because the method now says when it is wrong. Raise K until the rank stops rising, and check the result against the count from the trace integral. Both of those are diagnostics the caller can act on rather than parameters they have to guess.

So the moments do not make the method easy. They move it from having three guesses, one of which is silent, to having two guesses and one measurement — and the silent one was the dangerous one, because a probe block too small returns a plausible short list rather than an obvious failure.

What is still not lifted

The honest accounting, since this essay is about lifting a ceiling and two of the three remain.

The distance from the contour to the nearest outside eigenvalue still sets the quadrature’s rate, and no number of moments changes it. A contour drawn close to something outside it converges slowly whatever else is done.

The rounding at the end still applies. Eight points on a circle of radius four returned a count of minus two in the original measurement — a negative number of eigenvalues from a routine that did not fail — and the moments do nothing about that, because it is a property of the quadrature rather than of the extraction. The repair there is more points, and the diagnostic is the distance from the computed value to the nearest integer.

The problem’s own size is lifted by K along with the probe count, but only to K·n. A contour holding more than K·n eigenvalues is still under-resolved, and on a small problem with a large region that is the binding constraint rather than the probe block. The measurement at four probes, where n = 4 and the ceiling stops rising at K·4, is that case arriving.

The shape of the repair

It is worth naming the shape, because it recurs and is easy to miss at close range.

The original method computed one integral and threw away everything about the integrand except its value. The repair computes the same integral against a different weight — several different weights — at no extra cost, because the expensive part is evaluating the integrand and that has already been done at every quadrature point.

That is the same move as computing several moments of a distribution from one pass over data, or several inner products from one matrix–vector product, and it has the same economics: the access is the cost and the arithmetic on what was accessed is free. This site’s cost field is built on that observation in a different currency, where the access is a memory transfer rather than a linear solve, and the same arithmetic at a different price is the essay about it.

The prompt worth carrying: on meeting a method whose expensive step produces a quantity that is then reduced to a single number, ask what else could be extracted from the same step. Here the answer was a ceiling multiplied by an integer of the caller’s choosing.

At other settings

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 1 column and a circle of radius 5 holding 12 eigenvalues. The dashed line is the ceiling K·min(n, ℓ) = 1K, and the rank sits exactly on it while it binds. The solves are the expense and there are 256 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 itprobes1eigenvalues inside12rank at K = 11rank at K = 66solves, at every K256one ceiling per probeand K of them per moment
Fig. 3 One probe, where every moment adds exactly one to the ceiling.
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. 4 Three, where the ceiling reaches eighteen by K = 6.
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 4 columns and a circle of radius 5 holding 12 eigenvalues. The dashed line is the ceiling K·min(n, ℓ) = 4K, 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 1024 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.12345604812162024moments Krank of the block Hankel12 eigenvalues insidethe ceiling K·ℓthe ranka ceiling with a knob on itprobes4eigenvalues inside12rank at K = 14rank at K = 612solves, at every K1024one ceiling per probeand K of them per moment
Fig. 5 And four, where the problem’s own size begins to bind instead.
2 eigenvalues inside a circle of radius 2, and 1 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 1 contour moment against a probe block of 2 recover, agreeing to 2.44·10⁻¹⁵. With one moment the method could have returned at most 2 of them and would have said nothing about the rest.-2-1012-2-1012real partimaginary partan infinite spectrum, finitely askedinside the contour2probes2moments used1values returned4worst against the closed form2.4·10⁻¹⁵infinitely many eigenvaluesand a question with an answer
Fig. 6 A small circle, holding few eigenvalues and needing one moment.
3 eigenvalues inside a circle of radius 3, 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 7.11·10⁻¹⁵. With one moment the method could have returned at most 2 of them and would have said nothing about the rest.-3-2-10123-3-2-10123real partimaginary partan infinite spectrum, finitely askedinside the contour3probes2moments used2values returned6worst against the closed form7.1·10⁻¹⁵infinitely many eigenvaluesand a question with an answer
Fig. 7 A larger one.
12 eigenvalues inside a circle of radius 5, 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 5.02·10⁻¹⁴. With one moment the method could have returned at most 2 of them and would have said nothing about the rest.-6-4-20246-6-4-20246real partimaginary partan infinite spectrum, finitely askedinside the contour12probes2moments used6values returned24worst against the closed form5·10⁻¹⁴infinitely many eigenvaluesand a question with an answer
Fig. 8 Five, where the number inside has outrun a single moment several times over.
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 six, which is where the conditioning starts to be the binding constraint.
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. 10 What the lift costs, and the change of variable that removes it.
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. 11 The other integral, which answers how many and has no ceiling.
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. 12 The ceiling as it was measured before the moments lifted it.
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. 13 The infinite spectrum all of this is asked about.
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. 14 The decision the block Hankel’s rank is, and the gap it needs.
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. 15 Where a rank decision stops being one.
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. 16 Deciding a number is zero, which is what an empty contour needs.
How close Hager's estimate is to the true κ₁, over 200 seeded 10×10 matricesFive bars. The estimate is exactly the true condition number on 88% of the sample and inside ten per cent on 93%; the worst underestimate in the whole sample returns 45% of the truth. The last bar is the matrix built to defeat it, at 7.7%, well below anything the random sample reached.each bar is a percentage — of the sample, or of the true condition numberexactly right87.5%inside 10%93.0%inside a factor of 293.0%worst in the sample, ×10044.6%the constructed matrix, ×1007.7%usually exactexact share0.88worst of the sample0.45the constructed matrix0.077a routine that is right most of the timeand never wrong in the safe direction
Fig. 17 An estimate whose spread decides whether a threshold is safe.
Two trace estimators against their budget, on a 80×80 matrix whose spectrum decays at 0.7Two curves of relative error against the number of products with A, both axes logarithmic, as medians over 32 seeds. Hutchinson's fitted exponent is -0.31 and Hutch++'s is -5.07. The deflation changes the exponent rather than the constant, which is what makes it worth two thirds of the budget.10¹10²10⁻⁷10⁻⁶10⁻⁵10⁻⁴10⁻³10⁻²10⁻¹1products with Arelative errorHutchinsonHutch++measured at equal costfitted rate, Hutchinson-0.31fitted rate, Hutch++-5.1error at 96, Hutchinson0.051error at 96, Hutch++2.1·10⁻⁶both axes count products with Aso the sketch is paid for in the picture
Fig. 18 Reading a matrix through probes rather than through its entries.
160 single-probe estimates of one 10×10 trace, from the two standard probe distributionsTwo clouds of points, each one estimate of the trace from one random vector. Both are centred on the true value of 29.29 — the estimator has no bias at all — and the dashed bands are the standard deviation each distribution is known in closed form to have: 11.27 for the ±1 probe and 17.61 for the normal one. Measured over these draws they come out 10.64 and 20.47.0275481108135-129.514160.028290.5423121.056151.57drawzᵀAz from one probenormal±1two routes to one spreadthe trace29±1 spread, predicted11±1 spread, measured11normal ÷ ±11.6no bias in either cloudand one of them is narrower for a reason
Fig. 19 And the spread across probe draws, which is why ℓ is not one.
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. 20 The other route in this field, which returns everything and needs a filter.
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. 21 And the residual that does the filtering there.
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. 22 The approximation this route does not make.
σ_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. 23 The resolvent, which is what is being integrated.
Where the missing eigenvalues went: an 6 × 6 quadratic with 1 of its masses removedA chain of 6 masses with the first 1 of them set to zero, so M is singular of rank 5. det Q(λ) is a polynomial of degree at most 12; interpolated exactly in BigInt rationals at 13 nodes it has degree 11, so 1 of the 12 eigenvalues are at infinity — the same object a descriptor pencil has, arriving here because a degree of freedom with no inertia is an algebraic constraint. The float route counts the singular values of M judged to be zero and reaches 1, backed by a gap of ∞ between consecutive singular values. One integer, two routes, and only the second of them is a decision.finite eigenvalues (degree of det Q)11at infinity (2n − degree)1at infinity, by the rank of M12n, if M were nonsingular12a degree, not a decisiondegree of det Q11at infinity1by the rank of M1singular-value gap∞the count is a degreeand the other route is a judgement
Fig. 24 A count that is a degree rather than an integral.
The error of a rank-2 Tucker representation built from a sketch, against the columns of oversamplingThe flat line is the deterministic decomposition — d matrix SVDs — at 0.0671, and nothing random can go below it, since every curve here is a projection onto a subspace of the same size. The upper two are medians over five seeds: a dense Gaussian sketch, which the randomised field's bounds cover, and a Khatri–Rao sketch, whose columns are outer products of small random vectors and which no bound in that field applies to. With no oversampling the dense sketch is 2.78 times the deterministic answer; at twelve extra columns it is 1.20. The structured one costs 39.4 per cent more at four extra columns and 6.6 per cent at twelve — measured, because there is nothing else to say about it.-113579111310⁻¹1extra columns in the sketchrelative errora Khatri–Rao sketch: no bound covers ita dense Gaussian sketchthe decomposition, which nothing beatsrank 2, five seedsdeterministic0.067dense, p = 00.19dense, p = 120.081structured, p = 120.086structured ⁄ dense1.1a bound and a measurementand only one of them is available
Fig. 25 Oversampling a random probe, in the randomised field.
What extra columns buy a construction that never sees the matrix, over five seedsThe band is the range across five seeds and the line is their mean; the flat lower line is the best representation of rank 8 there is, which a decomposition of every block would find. With no oversampling the construction is 11.6 times behind it and the spread across seeds is 44 per cent of the mean. Eight extra columns bring it to 3.0× and the spread to 12 per cent. Both halves matter and only one of them is visible in a single run: what oversampling buys is a better answer and a more predictable one, and a figure drawn from one seed would have shown the first and hidden the second. No draw is ever better than the flat line, which is the check that the comparison is honest — a sample cannot beat the decomposition it is approximating.048121610⁻⁸10⁻⁷10⁻⁶10⁻⁵extra columns in the sample, p‖A − A_H‖ ⁄ ‖A‖the best rank-8 representationfive draws, mean and rangea band, not a lineexcess at p = 012excess at p = 83spread at p = 00.44spread at p = 80.12the optimum of this rank7.8·10⁻⁸one seed shows the meanand five show the risk
Fig. 26 And what it buys there.
The randomised SVD against the optimum it cannot beatA semi-logarithmic plot of approximation error against target rank. A shaded band shows the spread across seeds, a solid line the optimal error from the exact singular values, and a dashed line the published probabilistic bound well above both.04812162010⁻¹10⁻⁰.⁵1target rank k‖A − A_k‖₂published boundrandomisedσₖ₊₁, optimalhow far apart the three areworst seed spread1.6bound / median at k = 125.9median / optimum at k = 121.960×60, 6 seeds, oversampling p = 5band is best to worst
Fig. 27 A random block that has to have a component along what is wanted.
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. 28 The hero again, beside the count it should be checked against.
What 3 shifts do to each direction of the starting vectorAmplification of each eigen-component against its eigenvalue, on a logarithmic vertical axis. The points are measured — the component of the vector along each eigenvector, before and after — and the curve is |Π(λ − θ)| with the roots at the 3 discarded Ritz values. They agree to 3.9·10⁻¹¹. The 3 wanted directions are amplified by at least 1 and everything else by at most 0.802.1357910⁻⁵10⁻⁴10⁻³10⁻²10⁻¹1eigenvalue λamplificationkeptdiscardeda filter on the starting vectormeasured against the polynomial3.9·10⁻¹¹worst kept direction1best discarded direction0.8the roots are the discarded Ritz valuesand the vertical lines are where they sit
Fig. 29 A subspace method that also has to decide how many directions to keep.
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. 30 And a block size that has to exceed a multiplicity to see it.

Where this goes next

-2-101234-5-3-1135real partimaginary part4 insidea countable spectruminside the contour4drawn12existing∞worst branch residual1.6·10⁻¹⁵there is no last eigenvalueso the question has to change A problem with infinitely many eigenvalues Let the matrix depend on λ through something that is not a polynomial and three things stop being true at once. There is no linearisation, there is no characteristic polynomial, and "compute the spectrum" is not a request that can be granted — the only finite question is how many eigenvalues are inside this circle. 11.31.61.92.210⁻¹⁶10⁻¹³10⁻¹⁰10⁻⁷10⁻⁴10⁻¹10²log₁₀ quadrature pointsdistance from the true counthalf an eigenvaluean integer, eventuallytrue count2at 4 points2at 128 points2finest error1.4·10⁻¹³the integral is an integerand a rounding hides how far it was Counting what is inside a circle A trace of a matrix nobody wants to form, integrated around a contour, gives an integer — how many eigenvalues are inside. It converges exponentially, it is estimated with random probes, and the probe block is a ceiling that the answer does not mention. 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.

What links here

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

Reads more easily once this is understood

Essays that name this one as worth reading first.

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.

Complex arithmeticContour eigensolverContour integralExact ground truthMomentNonlinear eigenvalue problemNumerical rankQuadratureRank