A guarantee paid for in digits

Every eigenvalue real, and a test that says so introduced the class of quadratic eigenvalue problems whose spectrum is real as a property: Q is hyperbolic when there is a μ at which Q(μ) is negative definite, and one Cholesky that completes is the proof. It measured the class’s boundary on a damped chain, where at critical damping the arithmetic loses half its digits with nothing ill conditioned in the matrices, and a class a longer chain takes away found that the damping the class requires grows without bound with the chain.

Neither essay asked what computed the spectrum. It was the site’s general quadratic solver: a companion linearisation, its leading coefficient inverted, a nonsymmetric real Schur form. That route knows nothing about the class. A nonsymmetric eigensolver returns complex pairs whenever rounding pushes two close real eigenvalues onto a pair of complex conjugates, and near the boundary of the class two real eigenvalues are as close as they get. The class’s own theory offers a route that cannot do this, and the expectation a reader brings — the one a spectrum that comes in reciprocal pairs confirmed for another structure — is that a method that keeps the structure is the better method.

Near this boundary it is not. The guarantee is real, and it is paid for in exactly the digits the boundary threatens.

The route on which the spectrum has to be real #

A quadratic $Q(\lambda) = \lambda^2 M + \lambda C + K$, the object a matrix that depends on its own eigenvalue introduced, has linearisations of many kinds, and when its coefficients are symmetric it has symmetric ones: pencils $\lambda X + Y$ of twice the size, with $X$ and $Y$ symmetric, whose eigenvalues are those of $Q$. One family of them, built from an ansatz vector $[v_1, v_2]$, is $X = \begin{bmatrix} v_1 M & v_2 M \\ v_2 M & v_2 C - v_1 K \end{bmatrix}$ and $Y = \begin{bmatrix} v_1 C - v_2 M & v_1 K \\ v_1 K & v_2 K \end{bmatrix}$. When $Q$ is hyperbolic and $\mu$ is a point where $Q(\mu)$ is negative definite, the ansatz $[1, -\mu]$ makes the pencil definite: some combination $\cos t \cdot X + \sin t \cdot Y$ is positive definite.

A definite pencil is solved as a symmetric problem. Take the Cholesky factor $L$ of the positive definite combination, form the symmetric matrix $L^{-1}(-\sin t \cdot X + \cos t \cdot Y)L^{-\mathsf T}$, compute its eigenvalues with a symmetric solver — Jacobi here — and map each back to an eigenvalue of the pencil. The eigenvalues of a symmetric matrix are real, so every eigenvalue the route returns is real, whatever the rounding did. The same μ that certified the class builds the pencil, so the certificate is used twice: once to prove the spectrum is real and once to compute it in a way that keeps it so.

The cost of the route is in one number. The Cholesky divides by the positive definite combination’s least eigenvalue, and the best the pencil allows is its Crawford number: the largest least eigenvalue over all angles $t$. A perturbation of the pencil’s entries of size $\delta$ moves its eigenvalues by up to about $\delta$ over the Crawford number, so a backward-stable computation with errors of size $u$ in the entries delivers eigenvalues accurate to about $u$ over it.

Each route divides by something that vanishes #

The measurement uses the earlier essays’ chain: $n$ unit masses, a stiffness $K$ that is the second-difference matrix, and damping proportional to it, $C = \beta K$. Its spectrum has a closed form, one quadratic $\lambda^2 + \beta\kappa\lambda + \kappa$ per stiffness eigenvalue $\kappa$, and the critical damping is $\beta^* = 1/\sin(\pi/2(n+1))$ — 5.758770483143636 at eight masses. At $\beta = \beta^*(1 + \varepsilon)$ the chain is hyperbolic by a margin $\varepsilon$, and the two roots of the softest mode’s quadratic are about to collide.

The figure at the top of the page is the result at eight masses. At $\varepsilon = 10^{-2}$ both routes are accurate to fourteen digits: the companion route’s worst relative error is $8.9 \cdot 10^{-15}$ and the definite route’s $2.4 \cdot 10^{-14}$. They separate as the margin closes. At $10^{-6}$, $5.3 \cdot 10^{-13}$ against $1.5 \cdot 10^{-10}$; at $10^{-12}$, $2.7 \cdot 10^{-9}$ against $2.0 \cdot 10^{-3}$. Across the last six decades of margin the companion route’s error rises by 3.7 decades and the definite route’s by 7.1.

The two rates are the two denominators. The colliding pair, the quantity the companion route has to resolve, separates as $0.9823\sqrt{\varepsilon}$ at every margin measured, the square-root law the first essay found. A perturbation of size $u$ to a nonsymmetric matrix moves a pair of eigenvalues that close by about $u$ over their separation, so the companion route’s error should rise half a decade for each decade of margin. The Crawford number, the quantity the definite route divides by, is $0.2279\varepsilon$ to four digits from $10^{-4}$ to $10^{-12}$: linear in the margin, so the definite route’s error should rise a full decade per decade. The measured rises are 3.7 and 7.1 over six decades, both a little steeper than their denominators alone, and in the order the denominators set.

The first essay found that “the window of μ that works closes like the square root of the margin” while “the depth of the well” is a quarter of the margin, and that the certificate’s search is harder on the depth than on the width. The definite pencil inherits the depth. A μ always exists inside the class, which is why the route works at every margin, but the pencil it builds is only as definite as $Q(\mu)$ is negative, and that is linear in the margin.

Each route meets its own bound #

The denominators predict the rates; they also predict the errors, and the check is to multiply each error by its denominator and divide by the unit roundoff $u = 2^{-53}$. If a route’s error is its denominator’s bound and nothing more, the product is a small number that does not move with the margin.

For the definite route the product is between 0 and 4 at every margin on every length: at eight masses it reads 0, 1, 0, 1, 2 and 4 from $\varepsilon = 10^{-2}$ down to $10^{-12}$. The route delivers exactly what a backward-stable reduction of a pencil with that Crawford number can deliver, and not a digit less. Nothing in the Cholesky, the Jacobi sweeps or the map back from the reduced eigenvalues adds error of its own. For the companion route the product, error times pair separation over $u$, is between 0.7 and 36 across the same eighteen cells, larger and noisier because a nonsymmetric Schur form’s perturbation of a close pair depends on the angle between the pair’s left and right eigenvectors as well as on their distance — the quantity two condition numbers of one matrix separates from a norm — but bounded, with no trend in the margin.

So neither route is doing badly. Each is as good as its own analysis says it can be, and the comparison is between the analyses: a backward error of $u$ divided by a quantity linear in the margin, or divided by a quantity that is its square root. Near the boundary the square root is the larger denominator by a factor of $1/\sqrt{\varepsilon}$, which is a million at $\varepsilon = 10^{-12}$, and the two errors differ by roughly that.

No angle escapes the margin #

The Crawford number is the best of a family, and it is natural to suspect that the family was searched badly: perhaps some combination of $X$ and $Y$ other than the one found keeps a healthy least eigenvalue.

None does. Divided by the margin, the least eigenvalue of every combination traces one curve at $\varepsilon = 10^{-2}$, $10^{-4}$ and $10^{-6}$: 0.2153 at $t = 0$, a peak of 0.2279 at $t = 0.334$, 0.0901 at $t = 1.5$. The three curves cannot be told apart in the figure. So the entire definite range of the pencil shrinks in proportion to the margin, and the search was not the problem: the best angle is only six per cent better than the simplest, $t = 0$, which is $X$ itself. The pencil is definite at every margin and barely definite at all of them, and the shortfall is a property of the chain’s distance from its boundary, not of the linearisation’s construction.

The other choice in the construction is μ, and there the certificate already does the best that can be done. Any μ in the window where $Q(\mu)$ is negative definite builds a definite pencil, and the window is the interval between the two roots that collide. The certificate’s search, which minimises the largest eigenvalue of $Q(\mu)$, lands at the window’s midpoint to three digits at both $\varepsilon = 10^{-2}$ and $10^{-6}$. Measured across the window at those margins, the Crawford number over $\varepsilon$ is 0.043 a twentieth of the way in from either end, 0.171 a quarter of the way in, and 0.228 at the middle, the same at both margins. A μ chosen anywhere but the middle costs up to a factor of five more, and the middle is still linear in the margin. The certificate’s μ is the right one, and the right one is not enough.

Both pairs real, one far off #

At $\varepsilon = 10^{-12}$ the colliding pair sits near $-0.347296$, and the closed form puts its two members $9.82 \cdot 10^{-7}$ apart. The companion route returns them $9.8 \cdot 10^{-7}$ apart, both real, at the right place to the resolution of this axis and beyond: its worst error over the whole spectrum is $2.7 \cdot 10^{-9}$. The definite route returns them $6.94 \cdot 10^{-4}$ apart, both real, as it must. One of them is seven hundred times further from its partner than the closed form says. The route kept its promise exactly. Every eigenvalue is real; one of them is real and wrong in the fourth digit.

That is the shape of the trade at this boundary. The structured route guards a qualitative property and spends accuracy to do it. The unstructured route guards nothing and, measured, happens to keep the qualitative property as well as more of the accuracy.

Where the danger it guards against actually is #

The case for the structured route rests on the companion route’s ability to return a complex pair, so the next measurement looks for one. The proportional chain is a weak test, because its matrices are all diagonal in one basis; so each chain is rotated by a random orthogonal congruence $Q^{\mathsf T}(\cdot)Q$, which leaves the spectrum and the class untouched and gives the solver a dense problem with no structure to find.

The count needs care, and the care is a finding of its own. A real Schur form’s 2 × 2 blocks hold complex pairs, and the routine that reads them off decides by a tolerance which subdiagonals count as zero; at its default of $10^{-12}$ relative to the matrix, a genuinely complex pair whose block has a tiny subdiagonal would be read as two real numbers and hidden. So every block is read at tolerance zero, where a block is reported complex only if its discriminant is negative.

On twenty rotated chains at each of $\varepsilon = 10^{-4}$, $10^{-8}$ and $10^{-12}$ the companion route returns no complex value on any run, read either way. At $\varepsilon = 0$, the boundary itself, where the colliding pair is an exact double root, it returns a complex pair on five of twenty, with an imaginary part up to $5.3 \cdot 10^{-8}$, and on six at the default tolerance. The danger the guarantee guards against is real, and it is confined to a margin about as wide as the rounding: at $\varepsilon = 10^{-12}$ the pair is already $10^{-6}$ apart, which is a hundred times more than the companion route’s perturbation of it, and it stays real.

Every length, the same order #

The comparison does not depend on the chain. At $\varepsilon = 10^{-12}$ the companion route’s error is $9.0 \cdot 10^{-10}$, $2.7 \cdot 10^{-9}$ and $2.5 \cdot 10^{-9}$ on four, eight and twelve masses, and the definite route’s $5.6 \cdot 10^{-4}$, $2.0 \cdot 10^{-3}$ and $1.1 \cdot 10^{-3}$. The constants move with the length — the Crawford number is $0.6498\varepsilon$, $0.2279\varepsilon$ and $0.1130\varepsilon$, the pair separation $1.75\sqrt{\varepsilon}$, $0.982\sqrt{\varepsilon}$ and $0.682\sqrt{\varepsilon}$ — and the powers do not. A class a longer chain takes away found the critical damping growing with the length; at a fixed relative margin from it, both denominators shrink with the length and neither changes its law.

Why structure did not help here #

A spectrum that comes in reciprocal pairs is the case where a structured method is decisively better: a general solver computes the large half of a palindromic spectrum to full accuracy and the small half to seven digits, and a solver that pairs them recovers the rest. The structure there is a symmetry of the spectrum, and enforcing it removes an error the general solver makes. The structure here is a sign: every eigenvalue real. The general solver does not make the error the structure forbids, except within a rounding of the boundary, and the method that enforces it does so through a factorisation whose pivot is the margin itself.

Six routes to one spectrum found that linearisations equal in exact arithmetic differ in floating point by how they are reduced, and this is a seventh route with the same lesson, the one a backward-stable answer to a problem nobody asked put most sharply: a pencil can be exactly right, provably structured and backward stable in every step, and still divide by a number the problem has made small. The symmetric eigensolver at its heart is not the weak point; Jacobi computes the reduced matrix’s eigenvalues to full accuracy. The weak point is the reduction to it, a Cholesky of a combination whose least eigenvalue is $0.23\varepsilon$.

A code choosing between the routes has both denominators in hand before it computes anything, because the certificate’s search reports them. The first essay found that the window of μ that works is exactly as wide as the colliding pair’s separation, and that the well’s depth is a quarter of the margin; the Crawford number follows the depth. So the search that proves the spectrum real also says which route will compute it more accurately, and on every chain and margin measured here the answer is the one that does not enforce the proof.

What three chains do not show #

One family: proportional damping on a uniform chain, where the closed form makes every error exact. With damping not proportional to the stiffness there is no closed form, and the comparison would need a reference computed in higher precision. One structured algorithm, the Cholesky reduction of the best combination; there are algorithms for definite pencils that avoid forming that Cholesky, working with the indefinite symmetric structure directly, and they are designed for exactly this weakness. One unstructured algorithm, the same real Schur form used throughout, with its convergence tolerance at $10^{-14}$. The complex-pair survey is twenty rotations at four margins; a margin between $10^{-12}$ and zero, where the pair’s separation falls through the companion route’s perturbation of it, would locate the first complex pair more finely than this does.

Still open: a reduction that does not divide by the margin, and the margin where the pair turns complex #

A definite route without the Cholesky. Methods for definite pencils that work with the indefinite matrix $\cos t \cdot X + \sin t \cdot Y$ through an indefinite factorisation, or with hyperbolic rotations that preserve its inertia, do not divide by its least eigenvalue. The prediction with a sign is that a hyperbolic Jacobi method on this pencil, at $\varepsilon = 10^{-12}$ on the eight-mass chain, returns the spectrum with a worst relative error no larger than ten times the companion route’s — so that the guarantee and the accuracy stop being a trade.

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