The eigenvalue problem that is not linear

Real or complex by the sign of an error

A hyperbolic quadratic has only real eigenvalues, and an unstructured solver that knows nothing of the class returned a complex pair for one only at the boundary. The prediction was that the first complex pair would appear at a margin between 10⁻¹⁶ and 10⁻¹⁴ above the boundary, where the colliding pair's separation shrinks to the solver's own perturbation of it, and never at 10⁻¹³ or above. Both hold. Below 10⁻¹⁵ the solver returns a complex pair on almost exactly half of a hundred rotated chains — the same at the boundary itself — because its error in the pair's squared separation has the same size at every margin and either sign, and that one distribution, read at the boundary, predicts the whole curve. The structured route that guarantees a real answer loses every digit first.

Worth reading first: Every eigenvalue real, and a test that says so · A matrix that depends on its own eigenvalue.

A guarantee paid for in digits compared two ways of computing a hyperbolic quadratic’s spectrum near the edge of its class. A damped chain of masses (λ2M+λC+K)x=0(\lambda^2 M + \lambda C + K)x = 0 with damping C=βKC = \beta K is hyperbolic — every eigenvalue real, as a property certified by one Cholesky factorisation — for every β\beta above a critical value β∗\beta^*, and at β=β∗(1+ε)\beta = \beta^*(1 + \varepsilon) its two slowest eigenvalues sit 0.98ε0.98\sqrt\varepsilon apart, about to collide. A definite linearisation makes the real spectrum a property of the algorithm: it reduces the problem to a symmetric one, whose eigenvalues cannot be complex. It pays in digits, because its Cholesky divides by a Crawford number that vanishes like 0.228ε0.228\varepsilon. The unstructured companion route — linearise, invert, take a real Schur form — keeps its digits and knows nothing of the class, so in principle it could return the colliding pair as a complex pair. On twenty chains rotated by random orthogonal congruences, it did so only at the boundary itself, ε=0\varepsilon = 0.

The essay asked where, between 10−1210^{-12} and zero, the first complex pair appears. “The companion route moved the colliding pair by about 10−810^{-8} relative at ε=10−12\varepsilon = 10^{-12}, while the pair was 10−610^{-6} apart. The prediction is that on the twenty rotated chains the first complex pair appears at a margin between 10−1610^{-16} and 10−1410^{-14}, where the separation 0.98ε0.98\sqrt\varepsilon falls to the size of the companion route’s perturbation of it, and that no run returns one at 10−1310^{-13} or above.”

Both clauses hold, on a hundred chains rather than twenty. What happens below the margin where the first pair appears is the interesting part: the companion route does not start returning complex pairs. It starts tossing a coin.

A pair that can come back either way

The colliding pair belongs to the chain’s slowest mode, with stiffness k=2−2cos⁡(π/9)k = 2 - 2\cos(\pi/9) for eight masses. Its two eigenvalues are the roots of λ2+βkλ+k=0\lambda^2 + \beta k \lambda + k = 0, a centre −βk/2-\beta k/2 and a half-separation s=12β2k2−4ks = \tfrac12\sqrt{\beta^2k^2 - 4k}, real above the boundary. The pair comes back complex exactly when the computed squared half-separation is negative. So the measurement that matters is not the pair’s error but the error in s2s^2: call it η\eta, the computed squared half-separation minus the exact one, with a complex pair’s squared half-separation counted as minus its imaginary part squared. A run returns a complex pair exactly when η<−s2\eta < -s^2.

At fifteen margins, from 10−1210^{-12} down to the boundary, the chain is rotated by a hundred random orthogonal congruences — the same spectrum, a different matrix each time, so that the rounding is different each time — and the companion route’s colliding pair is located against the closed-form spectrum. The smallest margins are a few units in the last place of β\beta: below about 1.5⋅10−161.5 \cdot 10^{-16} the margin rounds away entirely, and β∗(1+2⋅10−16)\beta^*(1 + 2\cdot 10^{-16}) is stored as the same double as β∗\beta^* itself.

Where the first pair appears

The share of a hundred rotated eight-mass chains on which the companion route returns a complex pair, against the margin above the hyperbolic class's boundary, beside the share predicted from the errors the computation makes at the boundary itselfε 10⁻¹²: 0.00 measured, 0.00 predicted; ε 3·10⁻¹³: 0.00 measured, 0.00 predicted; ε 10⁻¹³: 0.00 measured, 0.00 predicted; ε 5·10⁻¹⁴: 0.00 measured, 0.00 predicted; ε 3·10⁻¹⁴: 0.00 measured, 0.00 predicted; ε 2·10⁻¹⁴: 0.00 measured, 0.00 predicted; ε 10⁻¹⁴: 0.03 measured, 0.04 predicted; ε 7·10⁻¹⁵: 0.11 measured, 0.11 predicted; ε 5·10⁻¹⁵: 0.18 measured, 0.21 predicted; ε 3·10⁻¹⁵: 0.24 measured, 0.33 predicted; ε 2·10⁻¹⁵: 0.37 measured, 0.35 predicted; ε 10⁻¹⁵: 0.45 measured, 0.41 predicted; ε 5·10⁻¹⁶: 0.46 measured, 0.47 predicted; ε 2·10⁻¹⁶: 0.46 measured, 0.47 predicted; ε 0: 0.47 measured, 0.47 predicted.00.10.20.30.40.50.6margin ε above the class's boundaryshare of runs with a complex pair10⁻¹⁶10⁻¹⁵10⁻¹⁴10⁻¹³10⁻¹²0none predicted to the rightmeasured, 100 runsdashed: predicted fromthe boundary's errorsmargin drawn on a logarithmic axis; the boundary at the far lefta coin, below the arithmetic's noise
Fig. 1 The share of a hundred rotated chains on which the companion route returns a complex pair, against the margin above the boundary, with the boundary itself at the far left. The dashed curve is the share predicted from the errors the computation makes at the boundary alone.

No run of a hundred returns a complex pair at any margin from 10−1210^{-12} down to 2⋅10−142 \cdot 10^{-14}. The first appear at 10−1410^{-14}, on three runs. Then the share climbs: eleven at 7⋅10−157 \cdot 10^{-15}, eighteen at 5⋅10−155 \cdot 10^{-15}, twenty-four at 3⋅10−153\cdot 10^{-15}, thirty-seven at 2⋅10−152 \cdot 10^{-15}, forty-five at 10−1510^{-15}. Below that it stops climbing. At 5⋅10−165 \cdot 10^{-16}, at 2⋅10−162 \cdot 10^{-16} and at the boundary itself it is forty-six, forty-six and forty-seven.

The prediction placed the first pair between 10−1610^{-16} and 10−1410^{-14}; it came at the top of that range, and nothing came at 10−1310^{-13} or above. The predicted mechanism — the separation falling to the size of the perturbation — is the right one, and the next figure shows it at work. What the prediction did not say, and what the curve does, is that a margin below the perturbation does not make complex pairs certain. It makes them a coin toss.

The error does not care about the margin

The colliding pair's exact squared half-separation and the median size of the companion route's error in it, against the margin, on a hundred rotated chainsε 10⁻¹²: squared half-separation 2.4·10⁻¹³, error 8.8·10⁻¹⁶; ε 3·10⁻¹³: squared half-separation 7.2·10⁻¹⁴, error 9.8·10⁻¹⁶; ε 10⁻¹³: squared half-separation 2.4·10⁻¹⁴, error 8.5·10⁻¹⁶; ε 5·10⁻¹⁴: squared half-separation 1.2·10⁻¹⁴, error 9.1·10⁻¹⁶; ε 3·10⁻¹⁴: squared half-separation 7.1·10⁻¹⁵, error 9.9·10⁻¹⁶; ε 2·10⁻¹⁴: squared half-separation 4.8·10⁻¹⁵, error 7.3·10⁻¹⁶; ε 10⁻¹⁴: squared half-separation 2.3·10⁻¹⁵, error 7.3·10⁻¹⁶; ε 7·10⁻¹⁵: squared half-separation 1.6·10⁻¹⁵, error 9.6·10⁻¹⁶; ε 5·10⁻¹⁵: squared half-separation 1.1·10⁻¹⁵, error 9·10⁻¹⁶; ε 3·10⁻¹⁵: squared half-separation 6.8·10⁻¹⁶, error 7.4·10⁻¹⁶; ε 2·10⁻¹⁵: squared half-separation 4.2·10⁻¹⁶, error 1.1·10⁻¹⁵; ε 10⁻¹⁵: squared half-separation 1.8·10⁻¹⁶, error 9.9·10⁻¹⁶; ε 5·10⁻¹⁶: squared half-separation 2.8·10⁻¹⁷, error 8.6·10⁻¹⁶. The error's size does not depend on the margin; the separation falls in proportion to it, and they cross between 10⁻¹⁵ and 3·10⁻¹⁵.10⁻¹⁷10⁻¹⁶10⁻¹⁵10⁻¹⁴10⁻¹³10⁻¹²margin εsize10⁻¹⁶10⁻¹⁵10⁻¹⁴10⁻¹³10⁻¹²exact half-separation²error in it, median sizethe crossing is the margin the arithmetic can seea separation the size of the noise
Fig. 2 The colliding pair’s exact squared half-separation and the median size of the companion route’s error in it, against the margin, on a hundred rotated chains. The error’s size is the same at every margin; the separation falls in proportion to the margin; they cross between 10⁻¹⁵ and 3·10⁻¹⁵.

The error η\eta has a median size between 7⋅10−167 \cdot 10^{-16} and 10−1510^{-15} at every one of the fifteen margins, with no trend. It is a property of the arithmetic acting on this matrix, not of how close the matrix is to the boundary. The exact s2s^2 is 0.24ε0.24\varepsilon — 2.4⋅10−132.4 \cdot 10^{-13} at a margin of 10−1210^{-12}, 2.3⋅10−152.3 \cdot 10^{-15} at 10−1410^{-14}, zero at the boundary — and the two curves cross between 10−1510^{-15} and 3⋅10−153 \cdot 10^{-15}. Above the crossing the separation is far larger than any error and the pair stays real; below it the error dominates.

The square root explains the earlier essay’s observation. An error of 10−1510^{-15} in a squared half-separation is an error of about 3⋅10−83 \cdot 10^{-8} in the half-separation itself, when the half-separation is small — the u\sqrt{u} sensitivity of a nearly double eigenvalue that one eigenvalue and two steps measured as a ring of radius u\sqrt u around a defective eigenvalue. At ε=10−12\varepsilon = 10^{-12} the half-separation is 5⋅10−75 \cdot 10^{-7} and an error of 3⋅10−83 \cdot 10^{-8} is the “10−810^{-8}” the earlier essay saw. The pair cannot turn complex until its exact separation is no larger than that, which is ε\varepsilon of a few times 10−1510^{-15}.

Why a coin

The error the companion route makes in the colliding pair's squared half-separation, on a hundred rotated chains at margin 3·10⁻¹⁵, with the threshold below which the pair comes back complexThe exact squared half-separation is 6.8·10⁻¹⁶. 24 of 100 errors fall below its negative, and those runs return a complex pair. The errors are spread from about -4.1·10⁻¹⁵ to 2.1·10⁻⁴.−10⁻¹⁶10⁻¹⁶−10⁻¹⁵10⁻¹⁵−10⁻¹⁴10⁻¹⁴−10⁻¹³10⁻¹³0−s²: complex to the leftη, the error in the pair's squared half-separation, on a signed logarithmic axisred dots: a complex pairthe errors stay; the threshold walks in
Fig. 3 The error in the colliding pair’s squared half-separation on a hundred rotated chains, on a signed logarithmic axis, at the margin on the dial; red where the pair came back complex. The shaded region is below −s2-s^2, where the pair turns complex. The cloud of errors does not move with the margin; the threshold does.

The dial shows why the curve flattens. The hundred errors form two clouds, one negative and one positive, each spread over about a decade around 10−1510^{-15}, with about as many in each. At a margin of 10−1210^{-12} the threshold −s2-s^2 sits at −2.4⋅10−13-2.4 \cdot 10^{-13}, two decades to the left of the negative cloud, and no error reaches it. At 3⋅10−153\cdot 10^{-15} it has walked into the negative cloud and twenty-four errors lie beyond it. At the boundary it is at zero, and every negative error turns the pair complex: forty-seven of a hundred. Once the threshold has passed through the negative cloud there is nothing more for it to take, and the share stays at the share of errors that happen to be negative.

So the companion route’s answer below the noise is not “complex”; it is “whichever way the rounding went”. The computation perturbs the matrix by a backward error of the size of the unit roundoff in a direction set by the rotation and the order of operations, and the pair’s discriminant moves by an amount whose sign depends on that direction. Nothing in the algorithm prefers one sign. The direction the error leans found accumulated errors leaning one way only when the rounding mode or the data made them; here neither does, and the errors split almost evenly.

How large a margin is, in a damping coefficient

It helps to translate the margins back into the problem’s own terms. The critical damping for this chain is β∗=5.76\beta^* = 5.76 in the chain’s units. A margin of 10−1410^{-14} is a damping of 5.76+5.8⋅10−145.76 + 5.8 \cdot 10^{-14}; a margin of 10−1510^{-15} is a damping that differs from critical in the fifteenth significant figure. No physical damping is known to fifteen figures, and no measured one to six. The region where the companion route tosses its coin is not a region any data could place a problem in with confidence. It is the region where the problem’s own definition — this chain, this damping — is finer than the number that stores it.

And it is narrow on a logarithmic scale as well. The coin appears between 10−1410^{-14} and 10−1510^{-15}, one decade, and everything above it is certainly real. The earlier essay’s caution about the unstructured route was that it could return a complex pair for a hyperbolic quadratic; the measurement says it does so only inside a decade of margins two orders below anything a caller could meaningfully ask for, and there it does so for a reason that has nothing to do with the route being unstructured. A condition number for one eigenvalue measured how one eigenvalue’s sensitivity can be orders larger than its neighbours’; a nearly double eigenvalue is the extreme case, with a sensitivity that is a square root and therefore unbounded relative to the size of the perturbation. Any route that computes it from rounded data inherits that, structured or not; the definite route’s structure only decides which property it keeps when the digits run out.

Rounding that takes both signs

The even split is the part of the result most worth understanding, because it is what makes the curve predictable. A backward-stable computation returns the exact eigenvalues of a nearby matrix, and the nearby matrix differs from the true one by perturbations of relative size about uu in directions set by the order of the floating-point operations. For the discriminant of the colliding pair, every such perturbation moves it by an amount linear in the perturbation, and the rotation of the chain — a different random orthogonal congruence each run — scrambles which directions the operations favour. Averaged over rotations, the discriminant’s error has no preferred sign, and the two clouds of the dial figure are the result: forty-seven negative errors in a hundred at the boundary, close to the fifty a fair coin would give.

A coin flip that fixes the average used this property on purpose: stochastic rounding makes each rounding error’s sign a fair coin, so that accumulated errors cancel instead of leaning. Here the fairness arrives unbidden, from the rotations, and it has the opposite effect on a caller. An answer that depends on the sign of one rounding error is as likely to be one thing as the other. The structured alternatives the last sections describe are ways to take the coin away — the definite route by forcing the sign, at the cost of digits; a symmetric solver, perhaps, by making the perturbation’s direction systematic.

One distribution predicts the curve

If the error’s distribution is the same at every margin, then the share of complex pairs at margin ε\varepsilon is just the share of errors below −s(ε)2-s(\varepsilon)^2, and the errors measured at any one margin should predict the share at all the others. That is a sharp test of the account above, and the dashed curve in the first figure is it: the hundred errors measured at the boundary, used to predict the share at each of the other fourteen margins.

The prediction is within 0.09 of the measured share at every margin: 0.04 against 0.03 at 10−1410^{-14}, 0.21 against 0.18 at 5⋅10−155 \cdot 10^{-15}, 0.33 against 0.24 at 3⋅10−153 \cdot 10^{-15} — the largest miss — and 0.41 against 0.45 at 10−1510^{-15}. The whole transition, from no complex pairs to a coin, follows from one number per run — the error in the pair’s squared separation — measured where the pair is exactly double. A model with no parameter fitted to the curve reproduces the curve.

Eight digits either way

The colliding pair as computed by the companion route on a hundred rotated chains, as offsets from the exact pair's centre in units of 10⁻⁷, at a margin of 10⁻¹⁴ and at the boundaryAt 10⁻¹⁴ the exact pair is 9.7·10⁻⁸ apart and 3 of 100 runs return it complex; at the boundary it is a double eigenvalue and 47 of 100 do. Every computed pair, real or complex, sits within about 10⁻⁷ of the exact centre — the eight digits the route keeps.-1-0.500.51-1-0.500.51real offset of the computed pair from its exact centre, units of 10⁻⁷imaginary part, units of 10⁻⁷margin 10⁻¹⁴at the boundaryon the axis: a real pair; off it: complexreal or complex by the sign of an error
Fig. 4 The colliding pair as computed by the companion route on a hundred rotated chains, as offsets from the exact pair’s centre in units of 10⁻⁷, at a margin of 10⁻¹⁴ and at the boundary. Real pairs lie on the horizontal axis, complex pairs on the vertical one; all lie within about 10⁻⁷ of the centre.

Real or complex, the computed pair is accurate. Drawn as offsets from the exact pair’s centre, every computed pair at both margins sits within about 10−710^{-7} of it — eight correct digits, the u\sqrt u the double eigenvalue allows. At 10−1410^{-14}, where the exact pair is 10−710^{-7} apart, most runs return two real eigenvalues straddling the centre at about the right distance and three return a complex pair with imaginary parts of 3⋅10−83 \cdot 10^{-8}. At the boundary, where the exact pair is a double eigenvalue, the runs scatter along both axes, real splits and complex splits of the same size.

A complex pair with an imaginary part of 3⋅10−83\cdot 10^{-8}, for a quadratic whose eigenvalues are certified real, is wrong in kind and right in size. A caller who rounds the result to seven digits sees two equal real eigenvalues either way. A caller who tests whether the imaginary part is exactly zero — the way a program decides whether a mode oscillates — gets a different answer on about half of the equivalent matrices.

The guarantee loses first

Both routes' worst relative error in the eight-mass chain's spectrum against the margin: the definite linearisation, which returns a real spectrum whenever its certificate accepts, and the companion routeε 10⁻¹²: definite 0.002, companion 2.7·10⁻⁹; ε 3·10⁻¹³: definite 0.0022, companion 5.1·10⁻⁹; ε 10⁻¹³: definite 0.015, companion 4.9·10⁻⁹; ε 5·10⁻¹⁴: definite 0.03, companion 7.3·10⁻¹⁰; ε 3·10⁻¹⁴: definite 0.0075, companion 5.6·10⁻⁹; ε 2·10⁻¹⁴: definite 0.0049, companion 1.4·10⁻⁹; ε 10⁻¹⁴: definite 0.0045, companion 3.4·10⁻⁸; ε 7·10⁻¹⁵: definite 0.096, companion 1.7·10⁻⁸; ε 5·10⁻¹⁵: definite 0.23, companion 5.8·10⁻⁸; ε 3·10⁻¹⁵: definite 0.16, companion 2.6·10⁻⁸; ε 2·10⁻¹⁵: definite 2.4, companion 1.9·10⁻⁸; ε 10⁻¹⁵: definite 0.71, companion 7.5·10⁻⁸; ε 5·10⁻¹⁶: definite 0.096, companion 9.9·10⁻⁸; ε 2·10⁻¹⁶: definite refuses, companion 4.6·10⁻⁸; ε 0: definite refuses, companion 3.7·10⁻⁸.10⁻¹⁰10⁻⁸10⁻⁶10⁻⁴10⁻²1margin εworst relative error in the spectrum10⁻¹⁶10⁻¹⁵10⁻¹⁴10⁻¹³10⁻¹²0××definite routecompanion route×: the certificate refusesthe definite route stays real and loses every digit firsttwo routes, two ways to fail
Fig. 5 Both routes’ worst relative error in the unrotated chain’s spectrum against the margin: the definite linearisation, which returns a real spectrum whenever its certificate accepts, and the companion route. Crosses mark the margins where the certificate refuses.

The structured route makes the opposite trade, and it reaches its failure sooner. The definite linearisation certifies the quadratic and returns a real spectrum at every margin down to 5⋅10−165\cdot 10^{-16}, and its error grows as its Crawford number shrinks: 2⋅10−32 \cdot 10^{-3} at a margin of 10−1210^{-12}, 4.5⋅10−34.5 \cdot 10^{-3} at 10−1410^{-14}, 0.160.16 at 3⋅10−153 \cdot 10^{-15}, 2.42.4 at 2⋅10−152\cdot 10^{-15} — an eigenvalue more than twice wrong — and 0.710.71 at 10−1510^{-15}. At 2⋅10−162 \cdot 10^{-16} and at the boundary the Cholesky certificate refuses: as stored, β\beta is β∗\beta^* and the quadratic is not hyperbolic. The companion route’s error over the same range is between 10−910^{-9} and 10−710^{-7}.

So below a margin of about 10−1410^{-14} neither route delivers both properties a caller might want. The definite route delivers a real spectrum with no correct digit; the companion route delivers eight correct digits and the wrong kind of number half the time. Every eigenvalue real, and a test that says so introduced the certificate as the class’s defining computation, and a class a longer chain takes away found the class shrinking as chains grow; this is where, for one chain, the class stops being distinguishable from its boundary in double precision.

What the coin means for a caller

The useful reading is not that the companion route is unreliable near the boundary. It is that near the boundary the question “real or complex” has no answer in double precision, and the companion route reports that honestly by answering at random. A margin of 10−1510^{-15} is a damping coefficient correct to fifteen digits, and no measured damping is. Any real problem near this boundary is a problem whose class is undetermined by its data, and an algorithm that returned a confident “real” there — as the definite route does, with 70 per cent errors — would be reporting a property of its own construction rather than of the problem. Six routes to one spectrum found the routes agreeing to rounding far from the boundary; here they disagree in kind, and the disagreement is the information.

What a caller can do with that is test the gap rather than the type. A pair whose computed separation, real or imaginary, is below about u\sqrt u times the eigenvalues’ size is a pair the arithmetic cannot classify, and it should be reported as a double eigenvalue of uncertain type — which is what the number that moves when the problem does would call a condition number of infinity for the question being asked.

What a hundred chains do not show

Every chain has eight masses and proportional damping, so its colliding pair is one mode’s and its centre is known in closed form. With non-proportional damping the colliding pair mixes modes, and the error in its separation would be expected to be larger and its sign distribution possibly lopsided; that is not measured. The rotations are random orthogonal congruences, which change the rounding without changing the problem; a structured solver that preserved symmetry throughout might make the error’s sign systematic, which is the second question below. And the near-even split of signs is measured, not derived: forty-five to forty-seven of a hundred is consistent with a fair coin, and why the computation’s backward error should move the discriminant symmetrically is argued from the absence of a preferred direction, not proved.

Still open: a sign the structure fixes, and a test on the gap

A symmetric solver’s sign. The companion route destroys the quadratic’s symmetry in its first step. A solver that keeps MM, CC and KK symmetric throughout — a symmetric linearisation reduced by congruences only — perturbs only symmetric data, and a symmetric perturbation of an overdamped pair’s discriminant may have a preferred sign. The prediction with a sign is that such a solver returns a complex pair on fewer than ten of a hundred rotated chains at the boundary, against the companion route’s forty-seven, while keeping its errors within a factor of ten of the companion route’s at every margin.

Testing the gap instead of the type. A classifier that reports “double, type undetermined” whenever the computed half-separation is below cuc\sqrt u times the pair’s centre would be right on every run where it answers at all. The prediction is that with c=2c = 2 it answers “real” on every run at margins of 3⋅10−143 \cdot 10^{-14} and above, “undetermined” on every run at 10−1510^{-15} and below, and never answers “complex” for a certified hyperbolic quadratic.

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 errorCompanion formDefective matrixDouble rootHyperbolic quadraticLinearisationOverdampingQuadratic eigenvalue problem