Numerical Linear Algebra
where the algebra stops being the arithmetic
Eigenvalues, singular values, rank

A problem with no answer

If two matrices share a null vector then det(A − λB) is identically zero and every λ is an eigenvalue, which means none of them is. Perturb such a pencil by a ten-billionth and a solver returns six numbers with residuals below 10⁻⁹. Change the seed and it returns six different numbers, spread over forty-four, with residuals just as small.
28 min read 23 figures Where the algebra stops being the arithmeticWhose fault is itExact ground truthAssertions that reject

Worth reading first: Two matrices and one problem · The condition number is an amplifier · The exact answer to a nearby problem.

A pencil has n eigenvalues because det(A − λB) is a polynomial of degree n. The previous essay was about the degree dropping when B is singular, and about the roots that go to infinity when it does.

There is one more thing that polynomial can do, and it is worse than dropping degree. It can be identically zero.

That happens as soon as A and B have a common null vector x: then (A − λB)x = 0 for every λ at once, so A − λB is singular for every λ, so the determinant is zero for every λ. Every complex number is an eigenvalue. The set of eigenvalues is the whole plane, which is another way of saying the question has no content.

Such a pencil is called singular — a word doing a different job from the one it does about a matrix, and one worth keeping separate: a pencil with a singular B is a perfectly good problem, and a singular pencil is not a problem at all.

Eigenvalues of a singular 6×6 pencil after a perturbation of 10^-8, over 8 seedsA and B share a null vector, so det(A − λB) is identically zero and every λ is an eigenvalue: the question has no answer. Perturbing by 10^-8 makes the pencil regular, and a solver then returns 6 eigenvalues with residuals no larger than 1.8·10⁻⁹ — small by any standard a caller would apply, and 3.7·10¹⁰ times smaller than the disagreement between seeds. Each row here is one seed. No two rows agree, the values are spread over 66.4, and nothing in the output distinguishes this from a well-posed problem. The only thing that says so is the exact characteristic polynomial, which has no coefficients at all.-40-27-14-11225380123456789computed eigenvalueseedevery mark has a residual below 10⁻⁸a small residual, and no answercoefficients of det(A − λB)0worst residual over all seeds1.8·10⁻⁹spread of the answers66seeds drawn8the residual is small at every markand none of the marks means anything
Fig. 1 One singular pencil, perturbed by 10⁻⁸, eight times with eight seeds. Each row is what a solver returned. Every mark has a residual below 10⁻⁹.

Where it comes from, which is a modelling error with no symptom

Nobody sets out to build one. Three ordinary things produce them.

A redundant equation. Two rows of a model that say the same thing, because a constraint was written twice or because a conservation law was imposed both globally and cell by cell. The two rows are dependent in both matrices at once, which is exactly the condition.

A variable nothing depends on. A state that appears in no equation, usually because a term was dropped during simplification or a coupling was set to zero to test something and never restored. Its column is zero in both matrices.

Or a symmetry that has not been quotiented out. A structure with a free rigid-body motion, a circuit with no ground node, a Laplacian with no boundary condition. The constant vector is in the null space of both operators, and it is there for a physical reason.

The last is the common one, and it is uncomfortable because the model is right. A floating structure really does have six rigid-body modes; the model really is invariant under them; and the pencil really is singular until somebody pins it down. What is wrong is asking an eigensolver about it.

Two penalties on the same problem, against the offset in the signalBest relative error for each penalty at four offsets. With no offset the two are within 1% of each other. At an offset of 10 the derivative penalty is 1.95 times better, because a constant lies in its null space and costs it nothing, while the norm penalty pays for the whole offset at every λ.best relative error at each offset‖x‖, offset 00.1906‖L₁x‖, offset 00.1880‖x‖, offset 20.0552‖L₁x‖, offset 20.0516‖x‖, offset 50.0310‖L₁x‖, offset 50.0241‖x‖, offset 100.0248‖L₁x‖, offset 100.0127what the null space buysadvantage at offset 01advantage at offset 21.1advantage at offset 51.3advantage at offset 101.9the norm penalty pays for a constantthe derivative penalty does not
Fig. 2 The standard remedy from the constrained-problem side: push the redundancy into the objective and watch what it does to the null space.
The arrowhead matrix, eliminated from each endThree sparsity plots. The first shows an arrowhead matrix with a dense first row and column. The second shows its Cholesky factor, completely dense. The third shows the factor obtained after moving the dense row to the end, which has no fill at all.the matrix43 entriestip eliminated first253 entriestip eliminated last43 entries‖A − LLᵀ‖/‖A‖, tip first1.4·10⁻¹⁶‖A − LLᵀ‖/‖A‖, tip last0dense factor is n(n+1)/2 = 253 · sparse factor is 2n − 1 = 43one row swapped to the endnothing numerical chose between them
Fig. 3 And the shape a constraint takes when it is kept rather than penalised, from the sparsity field.

What the exact route says, which is nothing at all

The polynomial is computed here in BigInt rationals, by evaluating the determinant at n + 1 integer points and interpolating, with nothing rounded anywhere.

For the pencil in the hero it comes back with zero coefficients. Not small coefficients. Not coefficients below a threshold somebody chose. The list is empty, because every value the determinant was evaluated at was exactly zero and the interpolation of a function that is zero at n + 1 points is the zero polynomial.

That is a fact rather than a judgement, and it is the only unambiguous statement anybody is going to get about this pencil. Everything else in this essay is what happens when it is not available.

The singular values of B, and the number of infinite eigenvalues the exact polynomial says there areThe pencil has 2 algebraic constraints, so 2 of its 6 eigenvalues are infinite. Two routes say so and they are not the same kind of statement. The exact route computes det(A − λB) in BigInt rationals and finds it has degree 4, so n − 4 = 2 eigenvalues are at infinity — an integer, arrived at without rounding. The float route counts the singular values of B below a cut, and the values run 4.27, 2.6, 0.918, 0.416, 1.13·10⁻¹⁶, 4.1·10⁻¹⁵⁷: the gap is a factor of 2.75·10¹⁴⁰ and it falls in the right place, so the decision is comfortable here. It is still a decision, and it is the same decision this site's rank essay is about.σ14.27σ22.6σ30.918σ40.416σ51.13·10⁻¹⁶σ610⁻¹⁸n = 6, and det(A − λB) has degree 4an integer, and a judgementdegree of det(A − λB), exactly4infinite eigenvalues, from the degree2singular values below the cut2largest gap in the spectrum2.7·10¹⁴⁰a degree cannot be nearly threeand a singular value can be nearly zero
Fig. 4 The comparison, from the previous essay: a regular pencil, where the degree is n minus the number of constraints and both routes agree about it.
log₁₀|det H_n| by three routes, to n = 22The exact rational determinant of the Hilbert matrix falls to 10^-274 at n = 22. The accumulated logarithm of the pivots follows it to n = 13 and then departs, reaching 10^-225 — wrong by 49 orders of magnitude and still an ordinary-looking number. The product of the pivots is a separate curve that ends at n = 22, where it underflows to exactly zero.25811141720-300-250-200-150-100-500nlog₁₀|det H_n|exact rationalΣ log|u_ii|Π u_iithree routes, one theoremexact at the largest n-274accumulated logarithm-225decades of disagreement49smallest pivot at that n4·10⁻¹⁷every pivot is a normal numberat every size on this axis
Fig. 5 And what an exact route is worth on a determinant generally, which this collection measured before it needed it here.

What a solver says, which is six numbers

Perturb A and B by ε times a random matrix. The pencil is now regular with probability one — the common null vector is gone, det(A − λB) has degree n again — and every routine that is asked returns n eigenvalues.

At ε = 10⁻⁸ the residuals of those eigenvalues are at most 1.8·10⁻⁹.

That number is the whole problem. It is small by any standard a caller applies. It is what a backward-stable solver produces on a well-posed problem. There is nothing about it that says the computation is meaningless, and every check a careful user would run comes back clean.

Then change the seed.

The eight rows in the hero are eight seeds. No two agree. The values are spread over forty-four — which is larger than any of them is — and the spread is 2.4·10¹⁰ times the residual.

Eigenvalues of a singular 6×6 pencil after a perturbation of 10^-2, over 8 seedsA and B share a null vector, so det(A − λB) is identically zero and every λ is an eigenvalue: the question has no answer. Perturbing by 10^-2 makes the pencil regular, and a solver then returns 6 eigenvalues with residuals no larger than 9.63·10⁻⁴ — small by any standard a caller would apply, and 6·10⁴ times smaller than the disagreement between seeds. Each row here is one seed. No two rows agree, the values are spread over 58.1, and nothing in the output distinguishes this from a well-posed problem. The only thing that says so is the exact characteristic polynomial, which has no coefficients at all.-39-26-1301326390123456789computed eigenvalueseedevery mark has a residual below 10⁻⁸a small residual, and no answercoefficients of det(A − λB)0worst residual over all seeds9.6·10⁻⁴spread of the answers58seeds drawn8the residual is small at every markand none of the marks means anything
Fig. 6 At a perturbation of 10⁻², where the residuals are larger and the scatter is not.
Eigenvalues of a singular 6×6 pencil after a perturbation of 10^-14, over 8 seedsA and B share a null vector, so det(A − λB) is identically zero and every λ is an eigenvalue: the question has no answer. Perturbing by 10^-14 makes the pencil regular, and a solver then returns 6 eigenvalues with residuals no larger than 2.11·10⁻¹⁵ — small by any standard a caller would apply, and 1.9·10²³ times smaller than the disagreement between seeds. Each row here is one seed. No two rows agree, the values are spread over 4.1·10⁸, and nothing in the output distinguishes this from a well-posed problem. The only thing that says so is the exact characteristic polynomial, which has no coefficients at all.-40-27-14-11225380123456789computed eigenvalueseedevery mark has a residual below 10⁻⁸a small residual, and no answercoefficients of det(A − λB)0worst residual over all seeds2.1·10⁻¹⁵spread of the answers4.1·10⁸seeds drawn8the residual is small at every markand none of the marks means anything
Fig. 7 And at 10⁻¹⁴, where the residuals are at the roundoff and the scatter is exactly what it was.

The scatter does not converge

That last pair is the difference between this and an ill-conditioned problem, and it is worth being precise about.

An ill-conditioned regular pencil has eigenvalues that move a great deal under a perturbation, and they move towards something: halve ε and the answers get closer to a limit, at whatever rate the conditioning allows. The answer exists and the computation is having difficulty reaching it.

Here there is no limit. Across six decades of ε the spread does not shrink. The answers are rearranged by the seed, at every size of perturbation, because there is no answer for them to be approximations to — the perturbation is not revealing a hidden true spectrum, it is creating one, and it creates a different one each time.

This is what an ill-posed problem looks like from inside a solver, and this collection has met the shape before in the regularisation field: a problem whose data does not determine its answer, where something outside the data has to choose. The difference is that a regularised problem knows it is choosing. Here nothing knows.

The widest box a uniqueness proof is available forVerified half-width against the distance between two roots, both axes logarithmic. The measured widths are 0.005, 0.015, 0.05 at separations of 0.01, 0.03, 0.1 — a ratio of 0.5000 at every one of them. The operator verifies exactly up to the point where the box would reach the second root, and a box containing both is refused.10⁻²10⁻¹10⁻²10⁻¹distance between the two rootswidest verified half-widththe separationverified widthuniqueness is a claim about a regionwidth ÷ separation at δ = 0.010.5width ÷ separation at δ = 0.030.5width ÷ separation at δ = 0.10.5verified up to the second rootand refused past it
Fig. 8 The regularisation field’s version: how wide the set of answers consistent with the data is, drawn rather than described.
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. 9 And the valley with no bottom — a fit whose coefficients can move by a third of their own size without the residual moving in the sixth digit.

The three things a caller checks, and what each of them says here

Put the ordinary diligence beside the ordinary output and the gap is the point of the essay.

Did it converge? Yes. The eigensolver reached its tolerance in the usual number of iterations, with no warning and no fallback path taken.

Is the residual small? Yes — below 1.8·10⁻⁹, which is what a backward-stable computation on a perturbed matrix produces and is several orders better than most people demand.

Are the eigenvalues plausible? Yes. They are ordinary real numbers spread over a range a modeller would look at without concern. Nothing is at 10³⁰⁰, nothing is NaN, nothing is complex when it should not be.

Three checks, three passes, and the answer is noise. The only check that fails is the one nobody runs: run it again with a different seed. That is why the last section of this essay is about a habit rather than about an algorithm.

Why the residual cannot see it

The residual asks whether the computed pair (λ, x) nearly satisfies Ax = λBx on the perturbed matrices, and it does, because the perturbed pencil is a genuine regular pencil with genuine eigenvalues and the solver found them accurately.

Every word of that is true and none of it is about the pencil that was handed over.

This is the sharpest instance this collection has of a statement it has been making since its second essay: a small residual is a statement about the problem the algorithm solved, not about the problem was posed, and when the two differ by a perturbation of 10⁻⁸ that is normally an excellent bargain. Here it is not a bargain at all, because the map from the problem to the answer is not continuous, so being close to the problem buys nothing about being close to the answer.

Backward and forward error against the condition numberA log–log plot over twelve decades of condition number. The backward error is a flat line at ten to the minus sixteen; the forward error rises in proportion to the condition number.110²10⁴10⁶10⁸10¹⁰10¹²10¹⁴10⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹condition number κ(A)relative errorforward errorbackward errorpredicted: κ · u50×50, 20 seeds per κ; dashed is the worstthe problem worsens, not the method
Fig. 10 The identity that normally rescues this, and the reason it does not here: forward error is condition number times backward error, and a condition number that is infinite makes the product uninformative whatever the second factor is.
Distance from the truth and residual, for both methods, with 60% of the noise in the matrixFour bars, medians over 40 seeds. The upper pair is how far each answer is from the coefficients the problem was built from; the lower pair is ‖Ax − b‖ on the problem as given. Ordinary least squares minimises the lower quantity by definition, so its bar is the shorter of the two whatever happens above — and at this share it is the less accurate answer.the upper pair is distance from the truth; the lower pair is ‖Ax − b‖least squares · error0.05224total least squares · error0.03071least squares · ‖Ax − b‖3.398total least squares · ‖Ax − b‖3.467two orderingserror ratio (ls ÷ tls)1.7residual ratio (tls ÷ ls)1seeds40no vector makes the residual smallernot even the one the problem was built from
Fig. 11 And a case from the least-squares field where the residual orders two methods backwards, which is the same failure at a smaller size.

What can be detected, and what it costs

The exact polynomial is not available on real data, so the question is what a float computation can say. Three things, in increasing order of cost and usefulness.

The cheap check: a common null vector. Compute the smallest singular value of A and of B, and if both are near zero, check whether their null directions agree. That catches the construction here and every case where the redundancy is a single vector. It costs two singular value decompositions and it misses everything more subtle than one shared direction.

The honest check: the smallest singular value of the pencil over a set of λ. A pencil is singular exactly when σ_min(A − λB) is zero for every λ. Sampling that function over a grid and finding it uniformly tiny is strong evidence; finding it small at isolated points is an ordinary spectrum. It is the same idea as a pseudospectrum, which this collection has already drawn, and it costs one singular value decomposition per grid point.

σ_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. 12 The one-matrix version of exactly that function: σ_min of A − zI over the plane, which is what a pseudospectrum is.
Where 24 perturbations of size 10⁻¹² put the eigenvalues of two 6×6 matrices with the same spectrumTwo clouds of eigenvalues in the complex plane. Both matrices have every eigenvalue at 0.8 exactly; both were perturbed by the same 24 random matrices of norm 10⁻¹². The normal matrix's eigenvalues stay within 6.51·10⁻¹³ of where they were — the size of the perturbation — and the bidiagonal's spread out to 0.0151, a factor of 2.3·10¹⁰ further on the same data.100.1750.35real partimaginary part0.8bidiagonalnormalone spectrum, two matricesthe perturbation10⁻¹²normal, furthest moved6.5·10⁻¹³bidiagonal, furthest0.015ratio2.3·10¹⁰the two matrices have identical eigenvaluesand one of them says so under perturbation
Fig. 13 And what it is diagnosing there — a spectrum that moves under a perturbation far more than the eigenvalues alone suggest it should.

And the right check, which nobody runs: the staircase. There is a canonical form for pencils — the Kronecker canonical form — that separates the regular part from the singular part exactly, and the staircase algorithms that compute it are the correct answer to this problem. They are also delicate, they need a sequence of rank decisions each of which can be wrong, and they are not what a caller gets by asking for eigenvalues.

The practical position is therefore uncomfortable and worth stating plainly: the standard tool returns an answer, the answer is meaningless, the diagnosis is available, and it is not run by default.

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 Why the diagnosis is delicate: every step of a staircase is one of these, and a wrong one changes the structure rather than a digit.
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 And the honest form of such a decision, which is a band rather than a verdict.

The perturbation is doing the choosing

There is a sharper way to say what the perturbation does, and it makes the failure sound less like bad luck and more like what it is.

A singular pencil is a point on a surface in the space of all pencils, and the surface has codimension greater than one — moving off it in any direction produces a regular pencil, and the eigenvalues that come out depend on the direction rather than on the distance. Halving ε moves half as far off the surface along the same ray and returns the same spectrum; changing the seed changes the ray and returns a different one.

So the answers are a function of the perturbation’s direction, which is exactly the thing a rounding error has no opinion about. Every arithmetic path through the computation is a different direction, which is why the same solver on the same data can give different answers on two machines, in two compilers, or with two orderings of the same sum — and the differences are the size of the answers rather than the size of the rounding.

This collection has an essay about a reduction that changes the order of a sum and therefore changes the result. That essay’s differences are in the last few digits. Here the same mechanism produces differences in the first.

Relative error of three summation algorithms in binary32A log–log plot of relative error against the number of terms for naive, pairwise and compensated summation, each measured against the exactly rounded sum.10¹10²10³10⁴10⁵10⁶10⁻⁹10⁻⁷10⁻⁵10⁻³10⁻¹number of terms addedrelative error against the exact sumin orderin a treecompensatedbinary32 · terms are 1/icompensated: 3·10⁻⁸
Fig. 16 The mechanism at its ordinary scale: the order of a sum changing its value, by an amount the arithmetic decides.

Two seeds are the cheapest thing available

If nothing else is done, one habit costs almost nothing and catches this.

Run it twice. Perturb the pencil by a different tiny random amount and compare the spectra. On a regular pencil, however ill-conditioned, the two runs agree to whatever the conditioning allows and the disagreement is bounded by it. On a singular pencil they disagree by the size of the answers.

That is this site’s two-routes habit applied to a case where the second route is the same route with a different seed — which is the weakest form of it, and here it is enough, because the failure being detected is precisely a failure of the answer to be a function of the input.

It also has the property that makes a check worth having: it produces a number rather than a verdict. The spread across seeds, divided by the size of the eigenvalues, is a continuous measure that is around 10⁻¹⁶·κ on a regular pencil and around one here.

160 single-probe estimates of one 12×12 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 37.239 — the estimator has no bias at all — and the dashed bands are the standard deviation each distribution is known in closed form to have: 14.18 for the ±1 probe and 21.23 for the normal one. Measured over these draws they come out 14.48 and 20.31.0275481108135322.849442.698862.548282.3975102.247drawzᵀAz from one probenormal±1two routes to one spreadthe trace37±1 spread, predicted14±1 spread, measured14normal ÷ ±11.5no bias in either cloudand one of them is narrower for a reason
Fig. 17 The habit in the field where it is standard, from the randomised essays: the seed changes the answer, so the figure draws a band and the spread is the finding.
Which of the two least-squares methods is more accurate, as a fixed amount of noise moves from b into AOne curve: the ordinary least-squares error over the total least-squares error, at each share of the noise placed in the matrix, as the median over 40 seeds. It runs from 0.28 when all the noise is in b — where ordinary least squares is the more accurate — to 2.47 when all of it is in A. The line at one is where the two methods are equally right.00.250.50.75110⁻¹110¹share of the noise placed in the matrixleast-squares error ÷ total least-squares errorequally accuratetotal leastsquares aheadordinary leastsquares aheadthe model, not the methodadvantage, all noise in b0.28advantage, all noise in A2.5seeds at each share40the same total noise at every pointand only where it sits changes
Fig. 18 And a second: how much of an answer is the data and how much is the draw.

What a library actually does about this

It is worth being fair to the software, because the situation is better than “nobody checks” suggests and worse than the check being present would imply.

LAPACK’s generalised eigenvalue routines return, alongside each eigenvalue, a pair of numbers that bound its conditioning — a reciprocal condition number for the eigenvalue and one for the associated eigenvector, available on request. On a singular pencil those come back tiny, because the eigenvalue is infinitely sensitive and the estimate says so as loudly as a finite number can.

So the diagnosis is available from the standard tool, at the cost of asking for it.

What it is not is the default. The routine most callers reach for returns the eigenvalues and nothing else; the conditioning estimates are behind a different entry point with more arguments; and the documentation that explains what a tiny one means is not the documentation somebody reads while getting a model to run.

That is a recognisable shape and this collection has met it twice already. A condition estimate that is computed and not printed; a residual that is available and not asked for; a backward error that every library could report and none does by default. The information exists, costs almost nothing, and sits one function call away from the person who needs it.

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. 19 The estimator that is one call away in every library, and the matrices on which it under-reports — which are the ones where the number would have been worth having.
What Hager's estimator returns, as a share of the truth, against the size of the matrix built to defeat itThe estimate over the true 1-norm against n, both axes logarithmic. It falls from 0.526 at n = 4 to 0.0658 at n = 32, along the line 1/t where t is the construction's own multiplier — which is bounded only by the size of the matrix, so the ratio has no floor. A line at one marks a correct estimate.10¹10⁻²10⁻¹1size of the matrixestimate ÷ true 1-norma correct estimatewhat it returns1 / tno floorratio at n = 40.53ratio at n = 320.066products, either size5the estimate is always a lower boundwhich is the direction that flatters the matrix
Fig. 20 And how far under-reporting can be pushed on purpose, which bounds how much a default estimate is worth.

The refusal

The assertion is fed a request for the eigenvalues of a singular pencil.

That is the failure this essay is about, in the one place where it can be caught: the routine has the exact characteristic polynomial in front of it, has found that it has no coefficients, and is being asked to return a finite set of roots anyway. Refusing is the difference between a routine that answers a question and one that answers a different question quietly.

The refusal is worth having even though the exact polynomial is rarely available, because it settles what the routine’s contract is. A function that returns eigenvalues of a pencil is claiming the pencil has some, and a function that cannot make that claim should say so rather than return the roots of whatever it was handed.

Eigenvalues of a singular 6×6 pencil after a perturbation of 10^-5, over 8 seedsA and B share a null vector, so det(A − λB) is identically zero and every λ is an eigenvalue: the question has no answer. Perturbing by 10^-5 makes the pencil regular, and a solver then returns 6 eigenvalues with residuals no larger than 1.36·10⁻⁶ — small by any standard a caller would apply, and 4.5·10⁷ times smaller than the disagreement between seeds. Each row here is one seed. No two rows agree, the values are spread over 61.5, and nothing in the output distinguishes this from a well-posed problem. The only thing that says so is the exact characteristic polynomial, which has no coefficients at all.-40-27-14-11225380123456789computed eigenvalueseedevery mark has a residual below 10⁻⁸a small residual, and no answercoefficients of det(A − λB)0worst residual over all seeds1.4·10⁻⁶spread of the answers62seeds drawn8the residual is small at every markand none of the marks means anything
Fig. 21 At 10⁻⁵ the same scatter with slightly larger residuals, which is the only thing about the picture that the perturbation size changes.
Eigenvalues of a singular 6×6 pencil after a perturbation of 10^-11, over 8 seedsA and B share a null vector, so det(A − λB) is identically zero and every λ is an eigenvalue: the question has no answer. Perturbing by 10^-11 makes the pencil regular, and a solver then returns 6 eigenvalues with residuals no larger than 1.89·10⁻¹² — small by any standard a caller would apply, and 2.2·10¹⁶ times smaller than the disagreement between seeds. Each row here is one seed. No two rows agree, the values are spread over 40800, and nothing in the output distinguishes this from a well-posed problem. The only thing that says so is the exact characteristic polynomial, which has no coefficients at all.-40-27-14-11225380123456789computed eigenvalueseedevery mark has a residual below 10⁻⁸a small residual, and no answercoefficients of det(A − λB)0worst residual over all seeds1.9·10⁻¹²spread of the answers4.1·10⁴seeds drawn8the residual is small at every markand none of the marks means anything
Fig. 22 And at 10⁻¹¹, where the residuals are at the roundoff and the answers still do not agree.
Eigenvalues of a singular 6×6 pencil after a perturbation of 10^-8, over 6 seedsA and B share a null vector, so det(A − λB) is identically zero and every λ is an eigenvalue: the question has no answer. Perturbing by 10^-8 makes the pencil regular, and a solver then returns 6 eigenvalues with residuals no larger than 1.8·10⁻⁹ — small by any standard a caller would apply, and 3.7·10¹⁰ times smaller than the disagreement between seeds. Each row here is one seed. No two rows agree, the values are spread over 66.4, and nothing in the output distinguishes this from a well-posed problem. The only thing that says so is the exact characteristic polynomial, which has no coefficients at all.-40-27-14-112253801234567computed eigenvalueseedevery mark has a residual below 10⁻⁸a small residual, and no answercoefficients of det(A − λB)0worst residual over all seeds1.8·10⁻⁹spread of the answers66seeds drawn6the residual is small at every markand none of the marks means anything
Fig. 23 Six seeds rather than eight, which changes how much of the disagreement is visible and none of the disagreement.

What is next

Three essays about a pencil, and the thread running through them has been what a zero in the second matrix means: a constraint, an eigenvalue at infinity, or the whole question falling over. The next pair changes subject and keeps the theme. They are about a claim this site has been making in a source comment since its first month — that one particular algorithm computes small singular values to high relative accuracy — which turns out to be true, and to be true of a class of matrices rather than of the algorithm.

Where this goes next

λ = 0λ = ∞λ = ∞-0.37690.5654.3836.4282 eigenvalues here, and it is one placecounted exactly, in rationalsfinite eigenvalues4at infinity2degree of det(A − λB)4worst residual, either kind5.5·10⁻¹⁴an eigenvalue is a ratioand a ratio has a direction, not a size An eigenvalue with no value If the second matrix of a pencil is singular then some of the eigenvalues are infinite, and that is not a degeneracy — it is the algebraic constraints of the model, one per constraint. What survives is a pair of numbers rather than one, and on the line those pairs live on, infinity is an ordinary point with an ordinary residual. 110²10⁴10⁶10⁸10¹⁰10¹²10¹⁴10⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹condition number κ(A)relative errorforward errorbackward errorpredicted: κ · u8×8, 20 seeds per κ; dashed is the worstthe problem worsens, not the method A small residual is not a small error Substituting the answer back and finding that it fits is the most natural check there is, and it verifies the wrong thing. A residual of 10⁻¹⁷ is entirely compatible with an answer whose second digit is wrong. 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 Where the answer stops being in the data The Picard condition finds the index where a noisy right-hand side stops carrying signal, from the data alone, with no knowledge of the answer. It lands at 32 where the truncation that actually minimises the error is 28 — and at 45 where the best is 38. It overshoots at every noise level from 1% to 0.001%, and it overshoots for a reason.

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.

Backward errorCounterexampleDeterminantExact arithmeticGeneralised eigenvalue problemIll posed problemMatrix pencilResidualSingular pencil