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

A spectrum that comes in reciprocal pairs

A palindromic quadratic reads the same backwards, so λ is an eigenvalue exactly when 1/λ is. A general solver discards that, computes the large half of the spectrum perfectly and the small half to seven digits — and the small half is a division away from being perfect too.
37 min read 32 figures Where the algebra stops being the arithmeticMeasured, not assertedIdentical algebra, different arithmetic

Worth reading first: A matrix that depends on its own eigenvalue · A nearby problem of the wrong kind · A factorisation with nothing to pivot for · The condition number is an amplifier.

A palindromic quadratic reads the same backwards:

Q(λ) = λ²A + λB + Aᵀ,      B = Bᵀ.

Multiply Q(1/λ) by λ² and what comes out is Q(λ)ᵀ, so det Q(1/λ)λ^2n = det Q(λ), and λ is an eigenvalue exactly when 1/λ is. The 2n eigenvalues arrive in n pairs whose products are one.

Problems of this shape come from periodic structures — a rail under a moving train, a waveguide, a surface acoustic wave filter — where λ is the factor by which a solution repeats from one period to the next and 1/λ is the same wave going the other way. The symmetry is not decoration; it is the physics.

It is also the third kind of structural constraint this field has met, and the three are worth keeping apart. Hyperbolicity constrains where the eigenvalues are — all real, in two groups. A singular leading coefficient constrains how many there are, as a degree rather than a rank decision. This one constrains how they pair, and it is the only one of the three that a computed spectrum can fail to have while every eigenvalue in it is individually reasonable.

The 12 eigenvalues of a palindromic quadratic, in 6 pairs whose products are oneλ²A + λB + Aᵀ with A and B symmetric reads the same backwards, so λ is an eigenvalue exactly when 1/λ is. The axis is log₁₀|λ|, the marks are the computed eigenvalues and each arc joins a pair; the vertical line is |λ| = 1, about which the whole picture is symmetric. The spectrum spans 9.56 decades at b = 10⁴. Measured, the general solver's spectrum departs from being closed under λ ↦ 1/λ by 3.27·10⁻⁹, and the structured route's — one symmetric generalised eigenvalue problem of size 6 rather than an unsymmetric one of size 12 — by 2.21·10⁻¹⁶, which is the rounding level and is by construction.-5-3.20611-1.412210.3816782.175573.969460log₁₀ |λ||λ| = 1λλ′ = 1, in the algebrapairs6decades of spectrum9.6pairing error, general3.3·10⁻⁹pairing error, structured2.2·10⁻¹⁶a symmetry the solver never knew aboutand the half of the answer it decides
Fig. 1 The eigenvalues on a log-modulus axis, each joined to its partner across |λ| = 1. The picture is symmetric about the vertical line by algebra; whether the computed numbers are is a measurement.

What a general solver returns

Linearise, run the real Schur factorisation, and 2n numbers come back with no relationship between them. The solver was never told about the symmetry and has no way to have inferred it.

The family here is built so that the truth is a closed form. A is the model stiffness matrix and B = b(A + I), so both are polynomials in the same tridiagonal matrix, they are simultaneously diagonalised by the sine basis, and the quadratic falls apart into n scalar palindromic quadratics aλ² + b(a + 1)λ + a = 0. The two roots of each multiply to a/a = one, exactly, and the large one is available without cancellation while the small one is its reciprocal. The parameter b moves the spread: the spectrum spans about 2k decades at b = 10ᵏ.

Measured at n = 6, with b from 10 to 10¹⁰:

b large eigenvalues small eigenvalues departure from λλ′ = 1
10 2.9·10⁻¹⁵ 1.1·10⁻¹⁴ 1.0·10⁻¹⁴
10⁴ 4.4·10⁻¹⁵ 3.3·10⁻⁹ 3.3·10⁻⁹
10⁶ 3.2·10⁻¹⁵ 6.6·10⁻⁸ 6.6·10⁻⁸
10¹⁰ 4.9·10⁻¹⁵ 1.0·10⁻⁷ 1.0·10⁻⁷

The large half is at the rounding level at every stop. The small half loses seven orders. And the third column is the same number as the second — the departure from the pairing is the error in the small eigenvalues, because each small eigenvalue is the reciprocal of a large one that is right.

None of that is surprising on its own. A spectrum spanning twenty decades asks a lot of a similarity transformation, and a small eigenvalue coming back with a poor relative error is the subject of small compared to what and accurate is not a property of a method. It is the setup for the sentence that is surprising.

The small half is a division away

Take the large eigenvalues — the ones that are correct to 5·10⁻¹⁵ — and invert them. What comes back is the small half, correct to 5·10⁻¹⁵.

b small half, directly small half, as 1/λ
10⁴ 3.3·10⁻⁹ 4.4·10⁻¹⁵
10⁶ 6.6·10⁻⁸ 3.2·10⁻¹⁵
10¹⁰ 1.0·10⁻⁷ 4.9·10⁻¹⁵

Seven orders of magnitude, for a division, out of a symmetry the solver discarded.

What a reciprocal pairing is worth: three ways to the small half of a palindromic spectrumA palindromic quadratic of size 6 whose spectrum spreads by two decades for every decade of b. The general solver returns the large eigenvalues to 5.52·10⁻¹⁵ at every stop and the small ones to 1.041·10⁻⁷ at the far end. Taking the reciprocals of the large ones instead returns the small ones to 4.89·10⁻¹⁵ — seven orders, for a division, out of a symmetry the solver discarded. The structure-preserving route, which enforces the pairing exactly, is the worst of the three on the eigenvalues themselves: 6.538·10⁻⁶, because each pair comes from one computed number and both members carry its error. Preserving a symmetry and using one are different acts with different prices.1357910⁻¹⁶10⁻¹³10⁻¹⁰10⁻⁷10⁻⁴log₁₀ b, the couplingrelative errorsmall half, directstructure-preservingsmall half, by 1/λone divisionlarge half, direct5.5·10⁻¹⁵small half, direct10⁻⁷small half, by 1/λ4.9·10⁻¹⁵structure-preserving6.5·10⁻⁶the structure is not decorationit is where half the accuracy is
Fig. 2 Three ways to the small half of the spectrum. The lowest curve is one division applied to numbers the solver already returned.

The mechanism is worth being explicit about, because it is the whole essay. The half of the spectrum a Schur factorisation gets right is the large half: relative accuracy at the small end of a widely spread spectrum is not something a similarity transformation preserves, and this site has said so twice before. The reciprocal symmetry says the small half is determined by the large half exactly. So the accurate route to the small eigenvalues does not go through the small eigenvalues at all.

That is the same move as computing the small root of a scalar quadratic from the product of the roots rather than from the quadratic formula — the cancellation that cancellation takes the answer, not a digit is about, avoided by using an identity instead of a subtraction. It is a first-year trick, and it recovers seven orders on a matrix problem where nobody thinks to look for it.

And preserving the structure costs the half that was right

The obvious alternative is a structure-preserving algorithm: one whose arithmetic cannot produce a spectrum that is not reciprocal. There is a clean one here.

Substitute λ = (1 + z)/(1 − z) and multiply by (1 − z)². With A symmetric the z term vanishes and what is left is

z²(2A − B) + (2A + B) = 0,

an ordinary symmetric generalised eigenvalue problem of size n, not 2n. Its eigenvalues come in ± pairs, and the Cayley map sends a ± pair to a reciprocal pair exactly. So the pairing is enforced by construction, at half the size, with a symmetric-definite pencil that this site already knows how to solve accurately.

Measured, its pairing error is 2.3·10⁻¹⁶ at every spread — flat, exact, as promised.

Its eigenvalues are worse. At b = 10¹⁰ the Cayley route’s large eigenvalues have a relative error of 6.5·10⁻⁶, against the general route’s 4.9·10⁻¹⁵. Nine orders, in the wrong direction.

The mechanism is the symmetry itself. Each pair comes from one computed z, so both members carry that number’s error — measured, the two halves of the Cayley spectrum have the same relative error at every spread, to within a factor of two. The general route confined the inaccuracy to the small eigenvalue, where it was harmless because the small eigenvalue was recoverable another way. Enforcing the symmetry spreads it onto the large one, where it is not recoverable at all.

So the best of the three routes is neither algorithm: take the general solver’s large half and divide. Exactly reciprocal, accurate to 5·10⁻¹⁵ in both halves, from the unstructured method and one division.

Preserving a structure and using one are different acts

That is the sentence to carry, and it is not the one the phrase “structure-preserving algorithm” suggests.

A structure can be imposed on the arithmetic, which forces the answer to have the property and constrains what the arithmetic may do to get there. Or it can be used as an identity after the fact, which leaves the arithmetic alone and applies the property where it helps. The first is what the Cayley route does and it costs nine orders here; the second is a division and it gains seven.

There are structures where the first is the right choice — where the property being preserved is the answer, rather than a relation among parts of it. A rotation that must stay a rotation is one: a reflection cannot stop being one is about exactly that, and there the constraint is free because the operation is orthogonal anyway. What separates the cases is whether preserving the property forces two quantities of different accuracy to share a number.

Why the small half is the half that is wanted

There is a temptation to shrug at the small eigenvalues: they are tiny, so who cares whether their relative error is 10⁻⁷ or 10⁻¹⁵. The application answers that, and the answer is why the essay exists.

For a periodic structure λ is a propagation factor — the complex number by which a wave is multiplied from one period to the next. |λ| < 1 is a wave decaying to the right, |λ| > 1 one decaying to the left, and the two members of a reciprocal pair are the same physical wave viewed in the two directions. Which of them is wanted depends on which end of the structure the boundary condition sits at, and both are wanted when the structure has two ends.

So the small half is not a numerical afterthought; it is half of the answer, and half of it is the half a general solver returns badly. In the worked family the small eigenvalues run down to 10⁻¹³, and their relative accuracy decides whether a computed mode decays at the right rate over a hundred periods. That is the same distinction between absolute and relative smallness that deciding that a zero has arrived turns on, arriving here as a question about half of a spectrum rather than about a single entry. A relative error of 10⁻⁷ compounded a hundred times is not a rounding.

What makes this family checkable

The closed form here is doing more work than usual and it is worth saying how it was arranged, since a reader may want to reproduce it.

A and B are chosen to be polynomials in the same tridiagonal matrix — A is the model stiffness and B = b(A + I) — which makes them commute and therefore simultaneously diagonalisable by the sine basis. In that basis the quadratic is n uncoupled scalar palindromic quadratics, each with roots available in closed form, and the whole spectrum is known exactly at every b and every n.

The reference is then computed the stable way, which is itself a demonstration of the essay’s point: the large root comes from the quadratic formula with the sign that avoids cancellation, and the small root comes from the product of the roots — that is, as 1/λ_big — rather than from the formula’s other branch. Computing the reference the naive way would have produced a “truth” with the same seven-order error as the thing being measured, and the measurement would have reported agreement.

That is a trap worth naming, since the site has walked into its relative once: a reference computed by the same route as the thing being checked is not a reference. Here it is avoided by construction, and the construction is the identity the essay is about.

What the Cayley route is still good for

The measurement is unkind to the structure-preserving route, and it should not be read as saying the idea is worthless.

It halves the size. The symmetric generalised problem is n × n where the linearisation is 2n × 2n, which is a factor of eight in the dense factorisation and a factor of two in memory — and it is symmetric, so the factorisation is the cheaper and better-behaved one this site has spent a whole field on. On a problem where the spectrum is not spread over twenty decades, that is a real win with no accuracy cost, because the shared-error mechanism only bites when the two members of a pair have very different magnitudes.

And it gives the pairing exactly, which is worth something when the pairing is what is being used downstream — a code that classifies modes by whether |λ| is above or below one wants that decision to be consistent, and a pairing error of 10⁻⁷ near |λ| = 1 makes it a coin toss.

The honest summary is a trade rather than a verdict: impose the structure when the spectrum is tight and the pairing is the output; use it as an identity when the spectrum is spread and the eigenvalues are. Which is a sentence that could only be written after measuring both.

The backward error does not see any of this

One more measurement, and it is the negative one that makes the finding sharp.

The structured backward error asks: what is the smallest perturbation of the coefficients that keeps the palindromic form and for which a computed pair is exact? It could in principle be much larger than the unstructured one, in which case ignoring the structure would be a stability question rather than a symmetry question.

It is not. Constructing an explicit symmetric perturbation ΔB with ΔBx = −r/λ — the minimum-norm symmetric matrix mapping x to a given vector, checked to be symmetric and checked to do the job — gives an upper bound on the structured quantity, and it comes out at most 1.41 times the unstructured backward error, with a median of 1.17, over twelve eigenpairs. Every pair is backward stable to begin with, at 10⁻¹⁵ or better.

So nothing in a backward-error analysis would have flagged the general solver. Its answers are exact for a palindromic problem a rounding away. What they are not is reciprocal, and no backward error measures that — which is the same gap a nearby problem of the wrong kind opens from the other side: there the nearby problem is not of the right kind, and here it is, and the computed spectrum still fails to have a property the exact one has.

What to print beside the answer

The operational summary is two lines of output rather than a change of algorithm.

Print the reciprocity. For each computed λ, the distance from 1/λ to the nearest other computed eigenvalue, relative to |1/λ|. It needs no reference answer, no closed form and no knowledge of the spectrum — it is a property of the returned set, computable in O(n²) comparisons, and on this family it is exactly the error in the small half. A solver that prints a residual and not this is printing the quantity that stays flat.

And produce the small half by division. If the reciprocity is large, the large eigenvalues are the trustworthy ones and their reciprocals are the small ones. That is a decision a caller can make from the printed number, and it costs n divisions.

Neither of those is a structure-preserving algorithm and neither requires one. They are what it looks like to use a symmetry rather than to impose it, and the measurement says the difference is seven orders of magnitude.

At other settings

The 12 eigenvalues of a palindromic quadratic, in 6 pairs whose products are oneλ²A + λB + Aᵀ with A and B symmetric reads the same backwards, so λ is an eigenvalue exactly when 1/λ is. The axis is log₁₀|λ|, the marks are the computed eigenvalues and each arc joins a pair; the vertical line is |λ| = 1, about which the whole picture is symmetric. The spectrum spans 3.56 decades at b = 10. Measured, the general solver's spectrum departs from being closed under λ ↦ 1/λ by 1.02·10⁻¹⁴, and the structured route's — one symmetric generalised eigenvalue problem of size 6 rather than an unsymmetric one of size 12 — by 1.74·10⁻¹⁶, which is the rounding level and is by construction.-2-1.20615-0.4122940.3815591.175411.969260log₁₀ |λ||λ| = 1λλ′ = 1, in the algebrapairs6decades of spectrum3.6pairing error, general10⁻¹⁴pairing error, structured1.7·10⁻¹⁶a symmetry the solver never knew aboutand the half of the answer it decides
Fig. 3 A weakly coupled problem, whose spectrum spans two decades and whose pairing survives at 10⁻¹⁴.
The 12 eigenvalues of a palindromic quadratic, in 6 pairs whose products are oneλ²A + λB + Aᵀ with A and B symmetric reads the same backwards, so λ is an eigenvalue exactly when 1/λ is. The axis is log₁₀|λ|, the marks are the computed eigenvalues and each arc joins a pair; the vertical line is |λ| = 1, about which the whole picture is symmetric. The spectrum spans 17.6 decades at b = 10⁸. Measured, the general solver's spectrum departs from being closed under λ ↦ 1/λ by 9.08·10⁻⁸, and the structured route's — one symmetric generalised eigenvalue problem of size 6 rather than an unsymmetric one of size 12 — by 2.16·10⁻¹⁶, which is the rounding level and is by construction.-9-5.87277-2.745550.3816783.50896.636130log₁₀ |λ||λ| = 1λλ′ = 1, in the algebrapairs6decades of spectrum18pairing error, general9.1·10⁻⁸pairing error, structured2.2·10⁻¹⁶a symmetry the solver never knew aboutand the half of the answer it decides
Fig. 4 And a strongly coupled one, at sixteen decades.
The 12 eigenvalues of a palindromic quadratic, in 6 pairs whose products are oneλ²A + λB + Aᵀ with A and B symmetric reads the same backwards, so λ is an eigenvalue exactly when 1/λ is. The axis is log₁₀|λ|, the marks are the computed eigenvalues and each arc joins a pair; the vertical line is |λ| = 1, about which the whole picture is symmetric. The spectrum spans 21.6 decades at b = 10¹⁰. Measured, the general solver's spectrum departs from being closed under λ ↦ 1/λ by 1.04·10⁻⁷, and the structured route's — one symmetric generalised eigenvalue problem of size 6 rather than an unsymmetric one of size 12 — by 2.32·10⁻¹⁶, which is the rounding level and is by construction.-11-7.20611-3.412210.3816784.175577.969460log₁₀ |λ||λ| = 1λλ′ = 1, in the algebrapairs6decades of spectrum22pairing error, general10⁻⁷pairing error, structured2.3·10⁻¹⁶a symmetry the solver never knew aboutand the half of the answer it decides
Fig. 5 The far end of the slider, where the pairing error is 10⁻⁷.
What a reciprocal pairing is worth: three ways to the small half of a palindromic spectrumA palindromic quadratic of size 3 whose spectrum spreads by two decades for every decade of b. The general solver returns the large eigenvalues to 2.54·10⁻¹⁵ at every stop and the small ones to 3.242·10⁻¹³ at the far end. Taking the reciprocals of the large ones instead returns the small ones to 2.62·10⁻¹⁵ — seven orders, for a division, out of a symmetry the solver discarded. The structure-preserving route, which enforces the pairing exactly, is the worst of the three on the eigenvalues themselves: 5.871·10⁻⁷, because each pair comes from one computed number and both members carry its error. Preserving a symmetry and using one are different acts with different prices.1357910⁻¹⁶10⁻¹³10⁻¹⁰10⁻⁷10⁻⁴log₁₀ b, the couplingrelative errorsmall half, directstructure-preservingsmall half, by 1/λone divisionlarge half, direct2.5·10⁻¹⁵small half, direct3.2·10⁻¹³small half, by 1/λ2.6·10⁻¹⁵structure-preserving5.9·10⁻⁷the structure is not decorationit is where half the accuracy is
Fig. 6 Three pairs rather than six: a narrower spectrum, and the same ordering of the three routes.
What a reciprocal pairing is worth: three ways to the small half of a palindromic spectrumA palindromic quadratic of size 10 whose spectrum spreads by two decades for every decade of b. The general solver returns the large eigenvalues to 7.31·10⁻¹⁵ at every stop and the small ones to 7.925 at the far end. Taking the reciprocals of the large ones instead returns the small ones to 4.56·10⁻¹⁵ — seven orders, for a division, out of a symmetry the solver discarded. The structure-preserving route, which enforces the pairing exactly, is the worst of the three on the eigenvalues themselves: 2.468·10⁻⁵, because each pair comes from one computed number and both members carry its error. Preserving a symmetry and using one are different acts with different prices.1357910⁻¹⁶10⁻¹³10⁻¹⁰10⁻⁷10⁻⁴log₁₀ b, the couplingrelative errorsmall half, directstructure-preservingsmall half, by 1/λone divisionlarge half, direct7.3·10⁻¹⁵small half, direct7.9small half, by 1/λ4.6·10⁻¹⁵structure-preserving2.5·10⁻⁵the structure is not decorationit is where half the accuracy is
Fig. 7 And ten.
The 16 eigenvalues of an 8 × 8 quadratic eigenvalue problem, computed and in closed formλ²M + λC + K for a chain of 8 masses with C = 0.3M + 0.1K. The crosses are the closed form — one scalar quadratic per eigenvalue of K, whose roots are known exactly — and the discs are the eigenvalues a real Schur factorisation returns from the 16 × 16 first companion linearisation. There are 16 of them for a matrix with 8 rows, of which 16 are complex and arrive in conjugate pairs, so the eigenvectors cannot be independent: 9 vectors in 8 dimensions never are. The worst disagreement between the two routes is 2.12·10⁻¹⁵, and the routes share nothing but the three coefficient matrices.00.0939569-2-1.24498-0.4899610.2650581.020081.7751real partimaginary parttwo routes, one spectrumeigenvalues16rows8complex16against the closed form2.1·10⁻¹⁵n rows and 2n eigenvaluesso the eigenvectors are not a basis
Fig. 8 A quadratic with no such symmetry, from the field’s first essay.
The 16 real eigenvalues of an overdamped chain, in two groups of 8, with the certificate between themA chain of 8 masses with C = 8K, which is above the critical damping 5.75877 = 1/sin(π/2(n+1)), so the quadratic is hyperbolic and every one of its 16 eigenvalues is real. The vertical line is the μ at which Q(μ) is negative definite — a Cholesky that completes, and the whole proof that the spectrum is real. It also lies in the gap: exactly 8 eigenvalues above it and 8 below, with a gap of 0.6698 between the two groups. One factorisation answers both questions, which is the second thing a boolean would not have carried.-33-28.2505-23.5009-18.7514-14.0018-9.25229-4.502750eigenvalueQ(μ) ≺ 0one Cholesky, two answersabove the certificate8below it8the gap0.67critical β for this n5.8the spectrum is real by classnot by outcome
Fig. 9 And one whose structure constrains where the eigenvalues are rather than how they pair.
Relative error of each singular value of a bidiagonal matrix graded over 30 decades, by four routesThe matrix is 8×8, bidiagonal, and every entry is a power of two times a small integer, so the matrix the float routines are handed and the matrix the exact routine reasons about are the same matrix. The reference is a Sturm bisection in BigInt rationals, which never rounds. σ_min is 2.098·10⁻³⁰. One-sided Jacobi holds every value to 4.41·10⁻¹⁶ relative and the zero-shift bidiagonal sweep to 2.21·10⁻¹⁶; the eigenvalues of BᵀB are wrong by 1 — not a rounding error but the whole number — and the damage is entirely at the small end.1234567810⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹singular value, largest firstrelative errorone-sided Jacobizero-shift QRshifted QReigenvalues of BᵀBagainst a rational bisectionσ_min, exactly2.1·10⁻³⁰worst, one-sided Jacobi4.4·10⁻¹⁶worst, zero-shift QR2.2·10⁻¹⁶worst, eigenvalues of BᵀB1a relative error is a ratioand the denominator is the answer
Fig. 10 Relative accuracy in the small eigenvalues of a symmetric matrix, which is the same difficulty one degree down.
Relative error of two algebraically identical expressions for (1 − cos x)/x², in binary64A log–log plot of relative error against x. The expression written as it reads loses accuracy as x falls and is eventually wrong in every digit; the rearranged form stays at rounding level.10⁻¹²10⁻¹⁰10⁻⁸10⁻⁶10⁻⁴10⁻²110⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹xrelative error of the computed value(1 − cos x)/x², as written2 sin²(x/2)/x²no digits left at allbinary64 throughoutone function, two spellings · zero below 1.5·10⁻⁸
Fig. 11 The cancellation the division avoids, in the essay that first measured it.
The same three schemes on a problem with no layer in itWorst nodal error against the grid size, both axes logarithmic, for a manufactured smooth solution on the identical operator at ε = 0.005. Adding no diffusion gives 0.0016, 4·10⁻⁴, 10·10⁻⁵, falling by four at each refinement. The tuned diffusion gives 0.067, 0.022, 0.0061 — 42 times worse at the coarsest grid, and falling more slowly.10²10⁻⁴10⁻³10⁻²10⁻¹grid points nworst nodal errorupwindtunedcentralthe same tuning, another problemtuned ÷ central at n = 3142tuned ÷ central at n = 12761central's error at the finest grid10·10⁻⁵exact on the problem it was derived fromand harmful on the one beside it
Fig. 12 A backward error whose nearby problem is not of the kind that was posed.
The two perturbations that make one computed solution exact, on a 12×12 Toeplitz systemA Kac–Murdock–Szegő matrix at ρ = 0.95, solved by Levinson's recursion. Every backward-error claim on this site says the computed answer solves a nearby problem exactly, and the nearby problem is the top matrix: the smallest perturbation of any kind, 2.18·10⁻¹⁷ relative, rank one, and constant along 0.98 of the way to none of its diagonals. The bottom matrix is the smallest perturbation that is itself a symmetric Toeplitz matrix — the same kind of object the problem was posed with — and it is 3.72·10⁻¹², larger by a factor of 1.71·10⁵. Both explain the same computed answer exactly. Only one of them is a problem anybody could have posed.the smallest perturbation of any kind — 2.18·10⁻¹⁷the smallest Toeplitz one — 3.72·10⁻¹²both exact for the same x̂smallest of any kind2.2·10⁻¹⁷smallest Toeplitz one3.7·10⁻¹²the price of the constraint1.7·10⁵diagonal defect, unconstrained0.98an exact answer to a nearby problemof a kind nobody posed
Fig. 13 The structured backward error against the unstructured one.
Relative eigenvalue error of two reductions of Ax = λBx, against κ(B), n = 12The pencil is built as B = LLᵀ and A = LΛLᵀ, so its generalised eigenvalues are exactly Λ and the error is a measurement rather than a comparison. Forming B⁻¹A gives a matrix whose departure from symmetry reaches 1.08 — about its own size — and reducing through the Cholesky factor gives one that is nearer symmetric by a factor of at least 6078. That difference is real and it does not appear in the accuracy: the two error curves have fitted slopes of 0.92 and 0.98 against κ(B) and stay within a factor of 2.3 of each other over fifteen decades. The conditioning belongs to the pencil, and no reduction of it escapes.10¹10⁴10⁷10¹⁰10¹³10¹⁶10⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹κ(B)relative error, and asymmetryvia B⁻¹Avia Choleskyasymmetry of B⁻¹Au · κ(B)against a spectrum known exactlyslope, via B⁻¹A0.92slope, via Cholesky0.98worst ratio between them2.3asymmetry of B⁻¹A1.1the symmetry claim is trueand it is not about the accuracy
Fig. 14 Two routes to a symmetric-definite pencil, which is what the Cayley substitution produces.
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. 15 A different repair for a different imbalance.
One quadratic eigenvalue problem in 9 systems of units: what the solver reports and what the answer is worthAn overdamped chain of 8 masses, with λ replaced by γμ so that the coefficients become (γ²M, γC, K). That substitution is exact in both directions and divides the spectrum by γ exactly, so the closed form is still available and every error here is measured against it. The backward error of the eigenpair for the LINEARISED MATRIX — the residual a solver's own error analysis is about — is 6.76·10⁻¹⁶ at γ = 1 and 7.59·10⁻¹³ at γ = 108 — it moves by a factor of 1928 while the other two move by 1.68·10¹⁰. The backward error for the QUADRATIC, which is what the person who posed the problem is entitled to, grows by 8.76·10¹⁰ across the same sweep, and the forward error follows it: 7.49·10⁻¹⁴ to 0.00126. Nothing went wrong with the solver at any stop.0246810⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹log₁₀ γ, the change of unitsrelative errorforward errorη, the quadraticη, the linearisationagainst a closed formη(linearisation), worst7.6·10⁻¹³η(quadratic), worst1.2·10⁻⁴forward error, worst0.0013coefficient spread4.2·10¹⁵the solver is right at every stopabout a problem nobody asked
Fig. 16 And the measurement that made the field’s units its first subject.
A double root, approached: the pair separates like √ε and the accuracy fails like √uA chain of 8 masses at β = β*(1 + ε), where β* is the critical damping and the smallest mode's two eigenvalues coincide at ε = 0. The separation of that pair is 0.9657√ε at every one of twelve decades — a spread of 1.02 in the constant — and the computed spectrum's worst relative error against the closed form rises as the pair closes, reaching 3.732·10⁻⁸ at the boundary itself. That is √u times a small constant: half the digits, on a problem where κ(K) is 32.16 at every stop and the coefficients are integers. Nothing here is ill conditioned in any sense this site has used before; what is unbounded is the derivative of the map from a coefficient to a double root.-15-13-11-9-7-5-3-110⁻¹⁶10⁻¹³10⁻¹⁰10⁻⁷10⁻⁴10⁻¹log₁₀ ε, distance past the critical dampingseparation, and relative errorat ε = 0: 3.73·10⁻⁸separationcomputed errornothing is ill conditionedthe √ε constant0.97spread of it, twelve decades1error at the boundary3.7·10⁻⁸κ(K), unchanged throughout32half the digitsand no condition number to blame
Fig. 17 A loss that no identity recovers.
A prediction from three numbers, exact at one end of the spectrum and wrong by n² at the otherThe tropical roots of the max-plus quadratic with coefficients ‖M‖, ‖C‖ and ‖K‖, divided by the extreme moduli they are supposed to estimate, for chains of 4 to 32 masses with C = 6M + 0.5K. The large root is within 5.7 per cent of the largest modulus at every size. The small root is out by 5.36, 17.1, 60.8 and 229.6 — a factor growing like the square of the size. The reason is structural: a norm is a maximum, the small end of this spectrum is set by the SMALLEST eigenvalue of K, which is 4sin²(π/2(n+1)), and no norm of K contains that number. A quantity built out of maxima is exact where the answer is a maximum and silent where it is a minimum.11.31.610⁻¹110¹10²10³log₁₀ ntropical root ÷ actual modulusexactsmallest moduluslargest modulusthree norms, two endslarge root, n = 40.96large root, n = 320.94small root, n = 45.4small root, n = 32230a maximum predicts a maximumand says nothing about a minimum
Fig. 18 Where a prediction from three norms is exact at the top of a spectrum and useless at the bottom.
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| = 2, computed as (1/2πi)∮tr(T(z)⁻¹T′(z))dz by the trapezoidal rule. The true count is 2, from the closed form. The error runs 0.00282, 0.0572, 0.0151, 5.25·10⁻⁴, 3.6·10⁻⁷, 1.35·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 count2at 4 points2at 128 points2finest error1.4·10⁻¹³the integral is an integerand a rounding hides how far it was
Fig. 19 And a spectrum with no ends at all.
The 12 eigenvalues of a palindromic quadratic, in 6 pairs whose products are oneλ²A + λB + Aᵀ with A and B symmetric reads the same backwards, so λ is an eigenvalue exactly when 1/λ is. The axis is log₁₀|λ|, the marks are the computed eigenvalues and each arc joins a pair; the vertical line is |λ| = 1, about which the whole picture is symmetric. The spectrum spans 5.56 decades at b = 100. Measured, the general solver's spectrum departs from being closed under λ ↦ 1/λ by 5.86·10⁻¹³, and the structured route's — one symmetric generalised eigenvalue problem of size 6 rather than an unsymmetric one of size 12 — by 1.88·10⁻¹⁶, which is the rounding level and is by construction.-3-1.87277-0.7455490.3816761.50892.636130log₁₀ |λ||λ| = 1λλ′ = 1, in the algebrapairs6decades of spectrum5.6pairing error, general5.9·10⁻¹³pairing error, structured1.9·10⁻¹⁶a symmetry the solver never knew aboutand the half of the answer it decides
Fig. 20 Four decades of spectrum, where the pairing error is still 6·10⁻¹³.
The 12 eigenvalues of a palindromic quadratic, in 6 pairs whose products are oneλ²A + λB + Aᵀ with A and B symmetric reads the same backwards, so λ is an eigenvalue exactly when 1/λ is. The axis is log₁₀|λ|, the marks are the computed eigenvalues and each arc joins a pair; the vertical line is |λ| = 1, about which the whole picture is symmetric. The spectrum spans 13.6 decades at b = 10⁶. Measured, the general solver's spectrum departs from being closed under λ ↦ 1/λ by 6.56·10⁻⁸, and the structured route's — one symmetric generalised eigenvalue problem of size 6 rather than an unsymmetric one of size 12 — by 1.78·10⁻¹⁶, which is the rounding level and is by construction.-7-4.53944-2.078880.3816782.842245.30280log₁₀ |λ||λ| = 1λλ′ = 1, in the algebrapairs6decades of spectrum14pairing error, general6.6·10⁻⁸pairing error, structured1.8·10⁻¹⁶a symmetry the solver never knew aboutand the half of the answer it decides
Fig. 21 Twelve decades, where it has reached 7·10⁻⁸.
What a reciprocal pairing is worth: three ways to the small half of a palindromic spectrumA palindromic quadratic of size 4 whose spectrum spreads by two decades for every decade of b. The general solver returns the large eigenvalues to 3.46·10⁻¹⁵ at every stop and the small ones to 1.626·10⁻¹¹ at the far end. Taking the reciprocals of the large ones instead returns the small ones to 1.56·10⁻¹⁵ — seven orders, for a division, out of a symmetry the solver discarded. The structure-preserving route, which enforces the pairing exactly, is the worst of the three on the eigenvalues themselves: 3.473·10⁻⁶, because each pair comes from one computed number and both members carry its error. Preserving a symmetry and using one are different acts with different prices.1357910⁻¹⁶10⁻¹³10⁻¹⁰10⁻⁷10⁻⁴log₁₀ b, the couplingrelative errorsmall half, directstructure-preservingsmall half, by 1/λone divisionlarge half, direct3.5·10⁻¹⁵small half, direct1.6·10⁻¹¹small half, by 1/λ1.6·10⁻¹⁵structure-preserving3.5·10⁻⁶the structure is not decorationit is where half the accuracy is
Fig. 22 Four pairs, and the same three curves in the same order.
What a reciprocal pairing is worth: three ways to the small half of a palindromic spectrumA palindromic quadratic of size 8 whose spectrum spreads by two decades for every decade of b. The general solver returns the large eigenvalues to 4.94·10⁻¹⁵ at every stop and the small ones to 4.533·10⁻⁴ at the far end. Taking the reciprocals of the large ones instead returns the small ones to 2.82·10⁻¹⁵ — seven orders, for a division, out of a symmetry the solver discarded. The structure-preserving route, which enforces the pairing exactly, is the worst of the three on the eigenvalues themselves: 9.185·10⁻⁶, because each pair comes from one computed number and both members carry its error. Preserving a symmetry and using one are different acts with different prices.1357910⁻¹⁶10⁻¹³10⁻¹⁰10⁻⁷10⁻⁴log₁₀ b, the couplingrelative errorsmall half, directstructure-preservingsmall half, by 1/λone divisionlarge half, direct4.9·10⁻¹⁵small half, direct4.5·10⁻⁴small half, by 1/λ2.8·10⁻¹⁵structure-preserving9.2·10⁻⁶the structure is not decorationit is where half the accuracy is
Fig. 23 And eight.
The 32 eigenvalues of a circulant, two waysA circulant matrix of size 32 has its eigenvalues in closed form: they are the discrete Fourier transform of its first column. Plotted against an eigensolver's answer for the same matrix, the two curves lie on top of each other to 4.4·10⁻¹⁵ relative. The eigenvectors are the same for every circulant of this size and are known before any entry is looked at.04812162024283210⁻¹1index kλthe transformthe eigensolverC = F* Λ F is a factorisationworst relative disagreement4.4·10⁻¹⁵‖Cx − b‖/‖b‖ from the transform solve9.6·10⁻¹⁶imaginary part of a real spectrum1.9·10⁻¹⁶n = 32, and the whole matrix is 32 numberseigenvectors known in advance
Fig. 24 A structure whose spectrum is known before anything runs, from the field next door.
The preconditioned spectrum on four grids (ρ = 0.9)Every eigenvalue of C⁻¹A, drawn as a point above the grid it belongs to, on a logarithmic vertical axis. The shaded band is within half a unit of one. The number of eigenvalues inside it goes 9, 11, 13, 17 while the number of unknowns goes 16, 36, 64, 100 — so the share clustered falls from 56% to 17%.357911110¹grid side meigenvalue of C⁻¹Awithin ½ of onethe cluster grows like the sidem = 4: inside of 169m = 6: inside of 3611m = 8: inside of 6413m = 10: inside of 10017the cluster grows with the sideand the spectrum with the area
Fig. 25 And one whose spectrum clusters, which is a third thing a structure can do.
Relative error of each singular value of a bidiagonal matrix graded over 2 decades, by four routesThe matrix is 8×8, bidiagonal, and every entry is a power of two times a small integer, so the matrix the float routines are handed and the matrix the exact routine reasons about are the same matrix. The reference is a Sturm bisection in BigInt rationals, which never rounds. σ_min is 0.00478. One-sided Jacobi holds every value to 1.2·10⁻¹⁵ relative and the zero-shift bidiagonal sweep to 1.09·10⁻¹⁵; the eigenvalues of BᵀB are wrong by 3.63·10⁻¹⁴ — not a rounding error but the whole number — and the damage is entirely at the small end.1234567810⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹singular value, largest firstrelative errorone-sided Jacobizero-shift QRshifted QReigenvalues of BᵀBagainst a rational bisectionσ_min, exactly0.0048worst, one-sided Jacobi1.2·10⁻¹⁵worst, zero-shift QR1.1·10⁻¹⁵worst, eigenvalues of BᵀB3.6·10⁻¹⁴a relative error is a ratioand the denominator is the answer
Fig. 26 Relative accuracy across two decades of spectrum rather than twenty.
The worst relative error over the whole spectrum, against how many decades the matrix is graded over, n = 8Each 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 2.1 decades — 3.63·10⁻¹⁴ — and by 11 decades it is at 14.5, 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 grading3.6·10⁻¹⁴and at the widest1.9·10⁷Jacobi, worst over the sweep1.5·10⁻¹⁵zero shift, worst1.1·10⁻¹⁵the definition is not a methodand squaring buries what it squares
Fig. 27 A graded matrix, where relative accuracy in the small entries is available and a general method does not get it.
The 32 eigenvalues of an 16 × 16 quadratic eigenvalue problem, computed and in closed formλ²M + λC + K for a chain of 16 masses with C = 0.3M + 0.1K. The crosses are the closed form — one scalar quadratic per eigenvalue of K, whose roots are known exactly — and the discs are the eigenvalues a real Schur factorisation returns from the 32 × 32 first companion linearisation. There are 32 of them for a matrix with 16 rows, of which 32 are complex and arrive in conjugate pairs, so the eigenvectors cannot be independent: 17 vectors in 16 dimensions never are. The worst disagreement between the two routes is 1.51·10⁻¹⁴, and the routes share nothing but the three coefficient matrices.00.095005-2-1.2366-0.4732080.2901891.053581.81698real partimaginary parttwo routes, one spectrumeigenvalues32rows16complex32against the closed form1.5·10⁻¹⁴n rows and 2n eigenvaluesso the eigenvectors are not a basis
Fig. 28 A quadratic with no pairing at all, for comparison.
A double root, approached: the pair separates like √ε and the accuracy fails like √uA chain of 8 masses at β = β*(1 + ε), where β* is the critical damping and the smallest mode's two eigenvalues coincide at ε = 0. The separation of that pair is 0.9657√ε at every one of twelve decades — a spread of 1.02 in the constant — and the computed spectrum's worst relative error against the closed form rises as the pair closes, reaching 3.732·10⁻⁸ at the boundary itself. That is √u times a small constant: half the digits, on a problem where κ(K) is 32.16 at every stop and the coefficients are integers. Nothing here is ill conditioned in any sense this site has used before; what is unbounded is the derivative of the map from a coefficient to a double root.-15-13-11-9-7-5-3-110⁻¹⁶10⁻¹³10⁻¹⁰10⁻⁷10⁻⁴10⁻¹log₁₀ ε, distance past the critical dampingseparation, and relative errorat ε = 0: 3.73·10⁻⁸separationcomputed errornothing is ill conditionedthe √ε constant0.97spread of it, twelve decades1error at the boundary3.7·10⁻⁸κ(K), unchanged throughout32half the digitsand no condition number to blame
Fig. 29 A structure that constrains where eigenvalues are, and what it costs at its boundary.
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| = 3, computed as (1/2πi)∮tr(T(z)⁻¹T′(z))dz by the trapezoidal rule. The true count is 3, from the closed form. The error runs 0.151, 0.532, 0.0874, 0.0263, 8.44·10⁻⁴, 7.21·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 1 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 count3at 4 points3at 128 points3finest error7.2·10⁻⁷the integral is an integerand a rounding hides how far it was
Fig. 30 And a problem where the spectrum has no ends to pair.
The 12 eigenvalues of a palindromic quadratic, in 6 pairs whose products are oneλ²A + λB + Aᵀ with A and B symmetric reads the same backwards, so λ is an eigenvalue exactly when 1/λ is. The axis is log₁₀|λ|, the marks are the computed eigenvalues and each arc joins a pair; the vertical line is |λ| = 1, about which the whole picture is symmetric. The spectrum spans 3.56 decades at b = 10. Measured, the general solver's spectrum departs from being closed under λ ↦ 1/λ by 1.02·10⁻¹⁴, and the structured route's — one symmetric generalised eigenvalue problem of size 6 rather than an unsymmetric one of size 12 — by 1.74·10⁻¹⁶, which is the rounding level and is by construction.-2-1.20615-0.4122940.3815591.175411.969260log₁₀ |λ||λ| = 1λλ′ = 1, in the algebrapairs6decades of spectrum3.6pairing error, general10⁻¹⁴pairing error, structured1.7·10⁻¹⁶a symmetry the solver never knew aboutand the half of the answer it decides
Fig. 31 The narrowest spread, where nothing has gone wrong yet.
What a reciprocal pairing is worth: three ways to the small half of a palindromic spectrumA palindromic quadratic of size 10 whose spectrum spreads by two decades for every decade of b. The general solver returns the large eigenvalues to 7.31·10⁻¹⁵ at every stop and the small ones to 7.925 at the far end. Taking the reciprocals of the large ones instead returns the small ones to 4.56·10⁻¹⁵ — seven orders, for a division, out of a symmetry the solver discarded. The structure-preserving route, which enforces the pairing exactly, is the worst of the three on the eigenvalues themselves: 2.468·10⁻⁵, because each pair comes from one computed number and both members carry its error. Preserving a symmetry and using one are different acts with different prices.1357910⁻¹⁶10⁻¹³10⁻¹⁰10⁻⁷10⁻⁴log₁₀ b, the couplingrelative errorsmall half, directstructure-preservingsmall half, by 1/λone divisionlarge half, direct7.3·10⁻¹⁵small half, direct7.9small half, by 1/λ4.6·10⁻¹⁵structure-preserving2.5·10⁻⁵the structure is not decorationit is where half the accuracy is
Fig. 32 And the widest problem the closed form is read at.

Where this goes next

the smallest perturbation of any kind — 1.89·10⁻¹⁷the smallest Toeplitz one — 2.65·10⁻¹³both exact for the same x̂smallest of any kind1.9·10⁻¹⁷smallest Toeplitz one2.7·10⁻¹³the price of the constraint1.4·10⁴diagonal defect, unconstrained0.97an exact answer to a nearby problemof a kind nobody posed A nearby problem of the wrong kind A good algorithm returns the exact answer to a nearby problem. A hundred and eighteen essays have measured the distance and not one has asked what the nearby problem looks like. On a Toeplitz system it is a rank-one matrix that is constant along none of its diagonals — and the smallest one that is Toeplitz is two and a half million times larger. -33-28.2505-23.5009-18.7514-14.0018-9.25229-4.502750eigenvalueQ(μ) ≺ 0one Cholesky, two answersabove the certificate8below it8the gap0.67critical β for this n5.8the spectrum is real by classnot by outcome Every eigenvalue real, and a test that says so A quadratic eigenvalue problem has no reason to have real eigenvalues. One class does, as a property rather than an outcome, and the proof is a Cholesky that completes. The boundary of the class has a closed form, and at the boundary the arithmetic loses half its digits with nothing ill conditioned anywhere. 1234567810⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹singular value, largest firstrelative errorone-sided Jacobizero-shift QRshifted QReigenvalues of BᵀBagainst a rational bisectionσ_min, exactly2.1·10⁻³⁰worst, one-sided Jacobi4.4·10⁻¹⁶worst, zero-shift QR2.2·10⁻¹⁶worst, eigenvalues of BᵀB1a relative error is a ratioand the denominator is the answer Small compared to what This site's own singular value routine has carried a sentence since the month it was written — that one-sided Jacobi computes the small singular values to high relative accuracy and the standard method does not. It has never been measured here, because measuring it needs a σ that is known rather than computed. A bidiagonal matrix and a Sturm count in exact rationals supply one.

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.

Named objects

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

Cayley transformGeneralised eigenvalue problemPalindromic quadraticQuadratic eigenvalue problemReciprocal pairRelative accuracyStructure preservingStructured backward error