The eigenvalue problem that is not linear

Six routes to one spectrum

Three linearisations of one quadratic, each reduced to a standard eigenvalue problem two ways. All six have exactly the same eigenvalues in exact arithmetic. On a well-scaled problem they differ by noise; on a badly scaled one by a factor of forty; and two of the six are the same matrix.

Worth reading first: A matrix that depends on its own eigenvalue · The exact answer to a nearby problem.

There is a family of essays on this site with one shape. Two procedures are the same algebra written in a different order; a derivation cannot tell them apart; and a computer can. Two Gram–Schmidts is the original — classical and modified, the same projections in a different sequence, returning 10⁻¹⁵ and 1 for the same measurement. The road that squares the problem is the same shape for a least-squares fit. One sequence and two recurrences is it for a Krylov method.

This is that essay for a quadratic eigenvalue problem, and it has an extra turn in it: there are six procedures rather than two, and on the problem a reader would naturally test them on, none of them can be told from any other.

Six ways of computing one spectrum, at a change of units of 106Three linearisations — the first companion form, the second, and the symmetric member of the DL(Q) family — each reduced to a standard eigenvalue problem in two ways: by inverting the leading coefficient, whose eigenvalues are λ, and by inverting the trailing one, whose eigenvalues are 1/λ. All six have exactly the eigenvalues of the quadratic in exact arithmetic. Measured against the closed form at γ = 106, the best is first/leading at 5.473·10⁻⁶ and the worst second/trailing at 2.234·10⁻⁴, a spread of 40.81. Under the leading reduction the symmetric form and the first companion form are THE SAME MATRIX — −A₁⁻¹A₀ is [[−M⁻¹C, −M⁻¹K], [I, 0]] for both — so the distinction between them is one the inversion discards, and it survives only under the trailing one.first · leading5.47·10⁻⁶first · trailing4.62·10⁻⁵second · leading4.27·10⁻⁵second · trailing2.23·10⁻⁴symmetric · leading5.47·10⁻⁶symmetric · trailing4.05·10⁻⁵all six are the same algebrabest route5.5·10⁻⁶worst route2.2·10⁻⁴spread across the six41condition of the linearisation8.3·10¹²the spectra agreeand the arithmetic does not
Fig. 1 Six ways of computing one spectrum, on a problem whose units are badly chosen. The slider takes the same six back to units that are not.

The six

Three linearisations. The first companion form C₁, which every textbook writes; a second companion form which puts the identity where C₁ puts K; and a symmetric member of the DL(Q) family, λ[[M,0],[0,−K]] + [[C,K],[K,0]], which is symmetric whenever M, C and K are and is therefore the one a structural engineer reaches for.

Two reductions. A linearisation is a pencil λA₁ + A₀, and a Schur factorisation wants a single matrix, so one of the two coefficients has to be inverted:

  • leading — form −A₁⁻¹A₀, whose eigenvalues are λ. The default.
  • trailing — form −A₀⁻¹A₁, whose eigenvalues are 1/λ, and invert them at the end. The route a code takes when the eigenvalues near zero are the ones wanted, because it sends them to the large end of the spectrum where a Schur factorisation resolves them best.

Six combinations, all of them exact in exact arithmetic, all of them ordinary.

Two of the six are the same matrix

The first thing the measurement shows is not about accuracy at all.

Under the leading reduction, the symmetric linearisation and the first companion form reduce to the same matrix. Write it out: for the symmetric form A₁ = [[M,0],[0,−K]] and A₀ = [[C,K],[K,0]], so

−A₁⁻¹A₀ = −[[M⁻¹, 0], [0, −K⁻¹]] [[C, K], [K, 0]] = [[−M⁻¹C, −M⁻¹K], [I, 0]],

which is exactly what C₁ gives. The two pencils are different objects — one is symmetric and the other is not — and the inversion throws the difference away. Whatever a symmetric linearisation is worth, it is worth it only to a method that keeps the pencil, which means QZ rather than QR, and the price of that is a factor of two or three in work.

Measured, the two routes agree to the last digit at every stop, which is the confirmation that this is an identity rather than a coincidence. The distinction between them survives only under the trailing reduction, where the matrices being inverted are genuinely different: A₀ is symmetric indefinite for the symmetric form and unsymmetric for C₁.

What the comparison anybody would run reports

Take the six, run them on the damped chain in sensible units, compare each computed spectrum with the closed form. The forward errors are

4.4·10⁻¹⁵ … 7.5·10⁻¹⁴,

a spread of seventeen at the fourteenth digit. That is noise. Not “small enough to ignore” — noise, in the sense that a different seed for the eigenvector iteration or a different compiler would reorder them. There is nothing in it to read.

Now the same six at a change of units of 10⁶:

route forward error
first companion, leading 5.5·10⁻⁶
symmetric, leading 5.5·10⁻⁶
symmetric, trailing 4.1·10⁻⁵
first companion, trailing 4.6·10⁻⁵
second companion, leading 4.3·10⁻⁵
second companion, trailing 2.2·10⁻⁴

A spread of 41, at the fifth digit, in a stable order. Now there is something to read: the second companion form is worse than the first by about eight, the trailing reduction is worse than the leading by about eight, and the two effects compound.

And the whole table is wrong. Every one of the six has lost ten digits, because the units are badly chosen and the previous essay fixes that with two lines. The best route on this problem is six orders worse than the worst route on the scaled one.

Which is the ordering, and it is the finding

That is the practical content of the essay, and it is worth stating flatly because it is the kind of thing that gets the wrong answer from intuition. There are two knobs. One of them is worth a factor of forty and the other is worth ten orders of magnitude. Fix the units first.

It is also worth noticing what the comparison of spectra could and could not have reported. On the well-scaled problem it reports noise, so a reader who tested there would conclude the choice is free. On the badly scaled problem it reports a real ordering — and a reader who tested only there would conclude that the linearisation is the thing to get right, which is exactly the wrong lesson from the same data. The measurement that separates them is the one that varies the units and holds everything else fixed, which is the shape this site keeps arriving at: a comparison needs a control, and the control here is a change of variable that is exact.

Why the second companion form is the worse one

There is a reason, and it is visible in the norms rather than in the algebra.

C₁ places an identity block where the second companion form places K. So ‖A₀‖ for C₁ is governed by max(‖C‖, ‖K‖, 1), and for the second form by max(‖C‖, ‖K‖) with a K in a second place. When the coefficients are balanced the two are nearly identical; when ‖K‖ is far from one, the second form concentrates more of the imbalance into the matrix being factorised. Measured at γ = 10⁶: ‖A‖ is 2.83 for the first form and 6.78 for the second under the leading reduction, and 8.7·10¹² against 2.8·10¹² under the trailing one. The condition numbers of the two matrices are 8.3·10¹² and 3.9·10¹², within a factor of two of each other — so the difference in the answers is not explained by the conditioning of the matrices being factorised, which is the first thing anybody would check.

What explains it is the same thing that explains the whole field: the map from a perturbation of the linearised matrix back to a perturbation of the three coefficients differs between the forms, and it is worse where the identity blocks are outnumbered.

The trailing reduction, and when it is the right one

The measurement above makes the trailing reduction look like a mistake, and on this problem it is. It is not one in general, and the reason is the one it was invented for.

A Schur factorisation resolves the large end of a spectrum first and best: the smallest eigenvalues of a matrix whose spectrum spans twenty decades come back with a relative error governed by the largest of them. If the eigenvalues wanted are the small ones — the low-frequency modes of a structure, which is nearly always what is wanted — then inverting the problem so that they become the large ones is the correct move, and it is the same idea as the shift-and-invert step in restarting is a filter.

What the measurement here says is narrower and still useful: on a problem whose spectrum is not widely spread, the trailing reduction costs a factor of eight and buys nothing. The reciprocal family where it buys a great deal is a spectrum that comes in reciprocal pairs, where the small half of the spectrum is exactly the reciprocal of the large half and taking it that way is worth seven orders.

What a linearisation has to be, and why it is not just any pencil

The word is doing precise work and it is worth unpacking once, here rather than in the field’s first essay, because this is where six different ones are being compared.

A pencil L(λ) = λA₁ + A₀ of size 2n is a linearisation of Q when there exist matrix polynomials E(λ) and F(λ) with constant nonzero determinants — unimodular, in the same sense the descriptor pencil essay uses the word — such that

E(λ) L(λ) F(λ) = [ Q(λ)  0 ]
                 [  0    I ].

The equality is between polynomial matrices, not between numbers, and that is what makes it stronger than “same eigenvalues”. It forces the elementary divisors to match, so a double eigenvalue of Q is a double eigenvalue of L, a defective one is defective to the same degree, and the Jordan structure carries across intact.

Why insist? Because without it, a linearisation could split a repeated eigenvalue into two simple ones, and a reader looking at two nearby computed numbers would have no way of telling a genuine pair from an artefact of the reformulation. All three used here are linearisations in the strict sense. So when the slider in a matrix that depends on its own eigenvalue walks two eigenvalues into a collision, the collision is a collision in every one of the six routes — and the half of the digits it costs, measured in every eigenvalue real, is a property of the problem rather than of any route to it.

Six routes and one eigenvector

Nothing above has mentioned eigenvectors, and the six differ there too, in a way that is easy to miss.

C₁’s eigenvector is [λx; x], so x can be read from either half. The second companion form’s is [x; Kx/λ], so the top half is the one wanted and the bottom carries a K in it. The symmetric form’s is [x; x/λ]. Under the trailing reduction every one of those is recovered from a vector belonging to 1/λ rather than λ, which exchanges the two halves’ roles again.

Six routes, and each of them offers two candidate eigenvectors with different backward errors, differing by a factor of |λ| — which on the badly scaled problem is 10⁻⁸. A library that takes a fixed half is leaving that factor on the table at every eigenvalue in the wrong group. The measurements in this field always take the better half, deliberately, so that the findings survive the most favourable reading of the method; a code that does not have a closed form to check against can take both and keep the one whose η is smaller, which costs two matrix–vector products.

What is the same across all six

It is worth listing what does not vary, because the list is longer than the one that does and it is what makes the comparison fair.

Every route computes the same 2n eigenvalues in exact arithmetic. Every route is backward stable for its own matrix, at the rounding level, at every stop — the Schur factorisation does not know which linearisation it was handed and its error analysis does not care. Every route inherits the same problem-level conditioning: κ for the quadratic is 4.98 whichever route is used, because it is a property of Q and the routes are all routes to Q.

And every route is destroyed by the same thing. At γ = 10⁸ all six have lost every digit; at γ = 1 all six are correct to fourteen. There is no route that is robust to badly chosen units, which is worth knowing precisely because a reader offered six options will look for the one that solves the problem. None of them does. Two lines of scaling do.

The one comparison that is not available

A reader might reasonably ask for the obvious experiment: perturb the coefficients, run all six, and see which one’s answer moves least. That would be a direct measurement of which linearisation inherits the least sensitivity, and it is not what is drawn here.

The reason is that it measures the wrong thing. A perturbation of (M, C, K) is a perturbation of the problem, and every route inherits it identically, because they all have exactly the eigenvalues of the perturbed quadratic. What differs between the routes is the perturbation each one manufactures out of the rounding in its own arithmetic, and the only way to see that is to compare against an answer that is exactly known. Which is why every number in this field is measured against a closed form and none of them against a better computation — a habit this site has had since an answer that is known and that this field needs more than any before it.

The seventh route, which is not on the list

Everything above forms a 2n × 2n object and factorises it. There is a route that does not, and it belongs in the list even though it is not measured here: build a subspace directly from the coefficients, project the quadratic onto it, and solve a small quadratic instead of a large matrix. That is a Krylov space for a problem that is not linear, and its arithmetic never forms a linearisation at all.

It is the right route when n is large enough that a dense 2n × 2n factorisation is out of reach, which is the ordinary case in the application this field comes from — a finite-element model of a damped structure has n in the hundreds of thousands and wants twenty eigenvalues near zero. What it inherits from this essay is the ordering: the subspace method is handed the same coefficients and the same units, and scaling them is as free and as necessary there as it is here.

What to take from six numbers that agree

The finding a reader should carry is not the ranking. It is that the experiment that would naturally be run to establish the ranking cannot see it.

Six spectra, on a problem in sensible units, agreeing to fourteen digits, is a measurement with no information in it — and it is the measurement anybody would make, because a well-scaled test problem is what a careful person constructs. The thing that has to vary for the difference to appear is the one nobody thinks of as a variable at all, since it is a change of units and changes of units are supposed to be free.

That is the same trap the site’s oldest essay is about, one level up. A residual can be small for a computation that has lost every digit; here a comparison of two computations can be empty for a choice that costs a factor of forty. In both cases what is missing is an answer known independently of the arithmetic being tested, and in both cases the fix is the same: construct the problem so that the truth is a closed form, and measure against it.

At other settings

Six ways of computing one spectrum, at a change of units of 100Three linearisations — the first companion form, the second, and the symmetric member of the DL(Q) family — each reduced to a standard eigenvalue problem in two ways: by inverting the leading coefficient, whose eigenvalues are λ, and by inverting the trailing one, whose eigenvalues are 1/λ. All six have exactly the eigenvalues of the quadratic in exact arithmetic. Measured against the closed form at γ = 100, the best is first/trailing at 4.383·10⁻¹⁵ and the worst first/leading at 7.489·10⁻¹⁴, a spread of 17.08. Under the leading reduction the symmetric form and the first companion form are THE SAME MATRIX — −A₁⁻¹A₀ is [[−M⁻¹C, −M⁻¹K], [I, 0]] for both — so the distinction between them is one the inversion discards, and it survives only under the trailing one.first · leading7.49·10⁻¹⁴first · trailing4.38·10⁻¹⁵second · leading2.28·10⁻¹⁴second · trailing8.86·10⁻¹⁵symmetric · leading7.49·10⁻¹⁴symmetric · trailing4.71·10⁻¹⁴all six are the same algebrabest route4.4·10⁻¹⁵worst route7.5·10⁻¹⁴spread across the six17condition of the linearisation452the spectra agreeand the arithmetic does not
Fig. 2 The six routes where the units are sensible: seventeen at the fourteenth digit, which is nothing.
Six ways of computing one spectrum, at a change of units of 102Three linearisations — the first companion form, the second, and the symmetric member of the DL(Q) family — each reduced to a standard eigenvalue problem in two ways: by inverting the leading coefficient, whose eigenvalues are λ, and by inverting the trailing one, whose eigenvalues are 1/λ. All six have exactly the eigenvalues of the quadratic in exact arithmetic. Measured against the closed form at γ = 102, the best is first/trailing at 4.159·10⁻¹³ and the worst first/leading at 5.18·10⁻¹², a spread of 12.46. Under the leading reduction the symmetric form and the first companion form are THE SAME MATRIX — −A₁⁻¹A₀ is [[−M⁻¹C, −M⁻¹K], [I, 0]] for both — so the distinction between them is one the inversion discards, and it survives only under the trailing one.first · leading5.18·10⁻¹²first · trailing4.16·10⁻¹³second · leading1.2·10⁻¹²second · trailing6.2·10⁻¹³symmetric · leading5.18·10⁻¹²symmetric · trailing7.12·10⁻¹³all six are the same algebrabest route4.2·10⁻¹³worst route5.2·10⁻¹²spread across the six12condition of the linearisation8.3·10⁴the spectra agreeand the arithmetic does not
Fig. 3 Two decades of badly chosen units, where the ordering starts to be stable.
Six ways of computing one spectrum, at a change of units of 104Three linearisations — the first companion form, the second, and the symmetric member of the DL(Q) family — each reduced to a standard eigenvalue problem in two ways: by inverting the leading coefficient, whose eigenvalues are λ, and by inverting the trailing one, whose eigenvalues are 1/λ. All six have exactly the eigenvalues of the quadratic in exact arithmetic. Measured against the closed form at γ = 104, the best is first/trailing at 2.638·10⁻⁹ and the worst second/trailing at 1.549·10⁻⁸, a spread of 5.872. Under the leading reduction the symmetric form and the first companion form are THE SAME MATRIX — −A₁⁻¹A₀ is [[−M⁻¹C, −M⁻¹K], [I, 0]] for both — so the distinction between them is one the inversion discards, and it survives only under the trailing one.first · leading8.57·10⁻⁹first · trailing2.64·10⁻⁹second · leading8.61·10⁻⁹second · trailing1.55·10⁻⁸symmetric · leading8.57·10⁻⁹symmetric · trailing5.87·10⁻⁹all six are the same algebrabest route2.6·10⁻⁹worst route1.5·10⁻⁸spread across the six5.9condition of the linearisation8.3·10⁸the spectra agreeand the arithmetic does not
Fig. 4 Four, where the spread is a factor of ten.
Six ways of computing one spectrum, at a change of units of 108Three linearisations — the first companion form, the second, and the symmetric member of the DL(Q) family — each reduced to a standard eigenvalue problem in two ways: by inverting the leading coefficient, whose eigenvalues are λ, and by inverting the trailing one, whose eigenvalues are 1/λ. All six have exactly the eigenvalues of the quadratic in exact arithmetic. Measured against the closed form at γ = 108, the best is first/leading at 0.001258 and the worst second/trailing at 1.545, a spread of 1229. Under the leading reduction the symmetric form and the first companion form are THE SAME MATRIX — −A₁⁻¹A₀ is [[−M⁻¹C, −M⁻¹K], [I, 0]] for both — so the distinction between them is one the inversion discards, and it survives only under the trailing one.first · leading0.00126first · trailing0.491second · leading0.813second · trailing1.55symmetric · leading0.00126symmetric · trailing0.402all six are the same algebrabest route0.0013worst route1.5spread across the six1229condition of the linearisation8.3·10¹⁶the spectra agreeand the arithmetic does not
Fig. 5 And eight, where every route has lost every digit.
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. 6 The repair that makes all six of them agree again.
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. 7 Why they disagree in the first place.
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 The spectrum all six are computing.
The condition number the problem has, and the one the solver's error analysis is written againstThe same eigenvalue of the same overdamped chain of 8 masses, in seven systems of units. Its condition number as an eigenvalue of the QUADRATIC — Tisseur's, with the three coefficient norms in the numerator and yᵀQ′(λ)x in the denominator — is 4.98 at γ = 1 and 4.98 at γ = 10⁶, a spread of 1 over six decades: it cannot move, because a change of units is not a change of problem. Its condition number as an eigenvalue of the LINEARISED MATRIX runs 49.66 to 7.844·10¹¹, a factor of 1.58·10¹⁰. The forward error follows the second one, and the first one is the honest description of the problem — so the substitution has manufactured an ill conditioning that belongs to the algorithm rather than to the question.012345610⁻¹10²10⁵10⁸10¹¹log₁₀ γ, the change of unitscondition numberthe linearisationthe quadraticone problem, two amplifiersκ(quadratic), first5κ(quadratic), last5κ(linearisation), last7.8·10¹¹how far the first moved1the problem is as well conditioned as everand the method is not
Fig. 9 And the condition number none of the six changes.
QᵀQ from classical Gram–Schmidt and from Householder on the 8×8 Hilbert matrixTwo eight-by-eight tables of QᵀQ. The upper one has ones on the diagonal and entries as large as one off it; the lower one is the identity to three decimal places everywhere.A = H8 · κ = 1.5·10¹⁰ · both factorisations reconstruct A to 6·10⁻¹⁷the diagonal is 1 in both — every column is a unit vector either way1.000000000001.000000000001.000000000001.000000000001.00000.002-0.002000001.0000.125-0.13300000.0020.1251.000-1.0000000-0.002-0.133-1.0001.000classical Gram–Schmidt1.000000000001.000000000001.000000000001.000000000001.000000000001.000000000001.000000000001.000Householderclassical ‖QᵀQ − I‖1.4Householder ‖QᵀQ − I‖1.4·10⁻¹⁵largest off-diagonal 1 against 3.1·10⁻¹⁶length is not angle
Fig. 10 Two Gram–Schmidts, the essay this one is shaped after.
Least squares by QR and by the normal equations in binary32Relative error of the computed coefficients against epsilon, on log axes. The QR route is a flat line near the bottom; the normal-equations route climbs and then stops, at the epsilon where the cross-product matrix becomes exactly singular.10⁻⁸10⁻⁷10⁻⁶10⁻⁵10⁻⁴10⁻³10⁻²10⁻¹10⁻⁹10⁻⁷10⁻⁵10⁻³10⁻¹10¹ε in Läuchli's matrix (smaller ε, larger κ)relative error in the coefficientsAᵀA exactly singularnormal equationsQRκ from 1.7·10⁸ to 17the cliff is at √u = 2.4·10⁻⁴
Fig. 11 The same shape for a least-squares fit.
Steps to a relative error of 10⁻⁶, on one 60×12 problem at five condition numbersTwo curves against the condition number on a logarithmic horizontal axis, with four seeds drawn at each. The two recurrences compute the same iterates in exact arithmetic. At κ = 10² they cost 16 and 16 steps; at κ = 10¹⁰ they cost 110 and 209.10²10⁴10⁶10⁸10¹⁰0285684112140168196224condition numberstepsnormal equationsbidiagonalisationone sequence, two costssteps at κ = 10², both16at κ = 10⁶, ratio1.1at κ = 10¹⁰, ratio1.9the same iterates in the algebraand twice the work at κ = 10¹⁰
Fig. 12 And for a Krylov recurrence.
The real Schur form with 2 conjugate pairs: 2 blocks that cannot be splitA square matrix drawn as a grid. Everything below the diagonal is zero except for a small number of two-by-two boxes on the diagonal, which are highlighted.3000000120000-21000000-0.5-1.500001.5-0.5000000-2T = ZᵀAZthe highlighted boxes each hold one conjugate pair, and no real rotation removes themthe form, and that it is one‖A − ZTZᵀ‖/‖A‖1.8·10⁻¹⁵‖ZᵀZ − I‖2.5·10⁻¹⁵worst eigenvalue error2.7·10⁻¹⁵surviving subdiagonal26×6, spectrum chosen before the matrix was builtquasi-triangular is as far as the reals go
Fig. 13 The factorisation all six routes end in.
Householder reduction to Hessenberg form, on a symmetric 6×6Two matrices side by side. The first is full; the second has zeros everywhere below the first subdiagonal, and being symmetric is tridiagonal.-0.16-0.63-0.493.43-0.025-0.632.3-0.69-0.95-1.70.076-0.49-0.694.3-1.6-1.41.93.4-0.95-1.65.40.0141.63-1.7-1.40.0144.60.055-0.0250.0761.91.60.0553.1A, symmetric-0.164.600004.65.92.500002.51-1.20000-1.24.90.3800000.385.80.2400000.242H = QᵀAQ, tridiagonal‖A − QHQᵀ‖/‖A‖1.1·10⁻¹⁵below the subdiagonal0worst eigenvalue movement7.1·10⁻¹⁵a similarity, so the spectrum is untouched — and every later step is O(n²) rather than O(n³)one reduction, then every iteration is cheapthe eigenvalues did not move
Fig. 14 And the reduction that precedes it.
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. 15 Two routes to a generalised eigenvalue problem, which is the same question one degree down.
Along a sequence of exactly rank-two tensors: the step's condition number, and what a fit from a random start achievesEvery A_n on this sequence *is* a rank-two tensor and its two rank-one terms are written down, so nothing here is about existence. The rising curve is the condition number of the r × r system each alternating sweep solves, which on this sequence has a closed form in n whose asymptote is 2n² — the marks are measured and the dashed line is that closed form, agreeing to 1.1·10⁻¹¹. The lower marks are what a three-hundred-sweep fit from a random start returns: 2.1·10⁻⁴ at n = 2 rising to 0.0118 at n = 1024. Started at the answer instead, the same code stays within 9.2·10⁻⁹ of it at every n — not the rounding level, because the drift from an exact start is itself about κ times the unit roundoff, but nine orders below what a random start reaches. That is the control that says the failure is the conditioning and not the implementation. A tensor away from the boundary conditions its step at 6.03.110¹10²10³10⁻¹²10⁻⁸10⁻⁴110⁴ncondition number, and residual reachedmarks above: κ of the step · dashes: its closed formmiddle: a fit from a random startbelow: the same fit started at the answera decomposition that is ill-conditionedκ at n = 10242.1·10⁶its closed form2.1·10⁶cosine of the terms1from a random start0.012from the answer9.2·10⁻⁹the answer existsand cannot be found
Fig. 16 What a decomposition inherits from the matrix it is taken of.
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. 17 The family where the trailing reduction earns its keep.
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. 18 And what taking the reciprocal is worth there.
Where binary32 stops being able to hold the linearisation, and where the scaled problem does notForming the first companion linearisation needs γ²M to be representable. The upper row is whether every entry of (γ²M, γC, K) is finite in binary32, at twelve changes of units; the lower row is the same question after Fan–Lin–Van Dooren scaling. The boundary is √(largest finite value ÷ largest entry of M) = 1.8447·10¹⁹, computed from the format's parameters and nothing measured — and the sweep agrees with it at every stop: the last change of units that survives is 10¹⁹ and the first that does not is 10²⁰. The scaled row never fails, because δγ²M has norm about one whatever γ was. In a narrow format the scaling is not an accuracy device: it is what makes the problem exist.02468101201change of units, by exponent10⁰10¹10²10³10⁴10⁶10¹⁰10¹⁶10¹⁹10²⁰10⁴⁰10¹⁵⁰as writtenafter scalinga range questionlargest finite value3.4·10³⁸predicted boundary γ1.8·10¹⁹last γ that forms10¹⁹stops where scaling fails0not a poor answerno answer at all
Fig. 19 Where the same badly chosen units stop being an accuracy question.
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. 20 The family every route is run on.
Two Krylov methods for one quadratic, priced by the numbers they storeA chain of 40 masses. Arnoldi on the 80 × 80 linearisation stores m vectors of length 80 and returns m Ritz values; the second-order recurrence stores m vectors of length 40 and its projection is taken onto the QUADRATIC, which has 2m eigenvalues for an m × m problem. So at every subspace size the second route holds half as many numbers and returns twice as many approximations — a factor of four in Ritz values per number stored, exactly and by construction. Plotted against storage, its curve is the lower one at every point: at 240 numbers it is 0.01379 from the dominant eigenvalue where the linearised route needs 480 to reach 0.07294. Both converge slowly — twenty-six vectors buy three digits — because the chain's spectrum is an arc with no isolated extreme, and that is the honest half: this is a comparison of prices rather than a fast method.33.310⁻⁴10⁻³10⁻²10⁻¹1log₁₀ numbers storeddistance to the dominant eigenvalueArnoldi, linearisedprojected quadraticper number heldstorage, linearised2080storage, second-order1040Ritz values, linearised26Ritz values, second-order52half the storageand twice the approximations
Fig. 21 And the seventh route, which does not form a linearisation at all.
One quadratic eigenvalue problem in 7 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 1.26·10⁻¹³ at γ = 106 — it moves by a factor of 320 while the other two move by 7.31·10⁷. The backward error for the QUADRATIC, which is what the person who posed the problem is entitled to, grows by 1.86·10⁸ across the same sweep, and the forward error follows it: 7.49·10⁻¹⁴ to 5.47·10⁻⁶. Nothing went wrong with the solver at any stop.012345610⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹log₁₀ γ, the change of unitsrelative errorforward errorη, the quadraticη, the linearisationagainst a closed formη(linearisation), worst1.3·10⁻¹³η(quadratic), worst2.6·10⁻⁷forward error, worst5.5·10⁻⁶coefficient spread4.2·10¹¹the solver is right at every stopabout a problem nobody asked
Fig. 22 The loss all six routes share, over six decades.
The condition number the problem has, and the one the solver's error analysis is written againstThe same eigenvalue of the same overdamped chain of 12 masses, in seven systems of units. Its condition number as an eigenvalue of the QUADRATIC — Tisseur's, with the three coefficient norms in the numerator and yᵀQ′(λ)x in the denominator — is 6.076 at γ = 1 and 6.076 at γ = 10⁶, a spread of 1 over six decades: it cannot move, because a change of units is not a change of problem. Its condition number as an eigenvalue of the LINEARISED MATRIX runs 60.31 to 9.461·10¹¹, a factor of 1.57·10¹⁰. The forward error follows the second one, and the first one is the honest description of the problem — so the substitution has manufactured an ill conditioning that belongs to the algorithm rather than to the question.012345610⁻¹10²10⁵10⁸10¹¹log₁₀ γ, the change of unitscondition numberthe linearisationthe quadraticone problem, two amplifiersκ(quadratic), first6.1κ(quadratic), last6.1κ(linearisation), last9.5·10¹¹how far the first moved1the problem is as well conditioned as everand the method is not
Fig. 23 And the condition number none of them changes.
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 4 masses, before and after Fan–Lin–Van Dooren scaling. The unscaled curve runs 7.01·10⁻¹⁵ to 4.331·10⁻⁴ — every digit gone by the far end — and the scaled one runs 7.67·10⁻¹⁵ to 4.77·10⁻¹⁵, flat to within a factor of 2.25. 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, worst4.3·10⁻⁴scaled, worst7.7·10⁻¹⁵orders recovered11scaled coefficient spread4.6the answer was never the problemthe units were
Fig. 24 The repair that makes the choice of route stop mattering.
Reversing the columns of a 6×6 matrix: what each orthogonalisation returnsBoth grids print the entries of Q₁P − Q₂, where Q₁ orthogonalises A and Q₂ orthogonalises A with its columns reversed. If the answer permuted with the columns, both grids would be zero. The polar factor's is: ‖U₁P − U₂‖_F = 1.53·10⁻¹⁵. QR's is not: ‖Q₁P − Q₂‖_F = 2.353 on matrices whose own Frobenius norm is 2.449, so the two answers are essentially unrelated. Gram–Schmidt and Householder both build Q one column at a time, and the first column is treated differently from the last.polar: U₁P − U₂QR: Q₁P − Q₂0000000000000000000000000000000000000.56-0.15-0.21-0.330.06-0.64-0.430.310.310.120.20.73-0.130.64-0.031-0.0028-0.420.075-0.18-0.43-0.130.21-0.29-0.420.640.630.670.62-0.52-0.420.26-0.520.022-0.0870.2-0.24the same columns, reordered‖U₁P − U₂‖_F1.5·10⁻¹⁵‖Q₁P − Q₂‖_F2.4‖Q‖_F, for scale2.4κ of the matrix10the column space did not moveand one of the two answers did
Fig. 25 Two orders for one orthogonalisation, which is this essay’s shape in the field it came from.
Classical and modified Gram–Schmidt: the same subtraction, in a different orderTwo panels of pseudocode differing in one argument, with the resulting pairwise dot products of the computed columns listed beneath each.for each previous column i, subtract the projection of column j onto qᵢclassicalr[i][j] = qᵢ · a[j] ↑ the ORIGINAL columnv = v − r[i][j] · qᵢthe three worst |qᵢ · qⱼ|:columns 7 and 8: 1columns 6 and 8: 0.13columns 6 and 7: 0.13modifiedr[i][j] = qᵢ · v ↑ what is LEFT of itv = v − r[i][j] · qᵢthe three worst |qᵢ · qⱼ|:columns 1 and 8: 4.4·10⁻⁷columns 2 and 8: 2.7·10⁻⁷columns 3 and 8: 2.4·10⁻⁸The two R factors agree to 1.2·10⁻⁶ relative. The two Q factors do not.the 8×8 Hilbert matrixone word, eight orders
Fig. 26 And the same pair measured as a loss of orthogonality.
The eigenvalues of a 6×6 pencil with 0 algebraic constraints, on the projective lineλ = α/β is a ratio, so an eigenvalue of a pencil is a direction rather than a number and it lives on a line whose two ends are the same point. Drawn as an angle φ = arctan λ, the 6 finite eigenvalues — -0.7663, -0.01384, 0.9593, 4.611, 5.61, 7.6 — sit inside the arc, and the 0 infinite ones sit at its ends, which is one point and not two. Nothing about them is degenerate: each carries a residual ‖βAx − αBx‖ in the same scaling as every other, the largest being 8.9·10⁻¹⁵, and an infinite eigenvalue's residual is ‖Bx‖/‖B‖ — the statement that its eigenvector is a null vector of B. The count is not a rank decision here: it is n minus the degree of det(A − λB), computed in exact rational arithmetic.λ = 0λ = ∞λ = ∞-0.7663-0.013840.95934.6115.617.6no constraints: B is nonsingular and nothing is at the polescounted exactly, in rationalsfinite eigenvalues6at infinity0degree of det(A − λB)6worst residual, either kind8.9·10⁻¹⁵an eigenvalue is a ratioand a ratio has a direction, not a size
Fig. 27 A pencil with no infinite eigenvalues, for comparison with the reduction that inverts the trailing coefficient.
The eigenvalues of a 6×6 pencil with 4 algebraic constraints, on the projective lineλ = α/β is a ratio, so an eigenvalue of a pencil is a direction rather than a number and it lives on a line whose two ends are the same point. Drawn as an angle φ = arctan λ, the 2 finite eigenvalues — 0.4384, 4.562 — sit inside the arc, and the 4 infinite ones sit at its ends, which is one point and not two. Nothing about them is degenerate: each carries a residual ‖βAx − αBx‖ in the same scaling as every other, the largest being 1.55·10⁻¹⁵, and an infinite eigenvalue's residual is ‖Bx‖/‖B‖ — the statement that its eigenvector is a null vector of B. The count is not a rank decision here: it is n minus the degree of det(A − λB), computed in exact rational arithmetic.λ = 0λ = ∞λ = ∞0.43844.5624 eigenvalues here, and it is one placecounted exactly, in rationalsfinite eigenvalues2at infinity4degree of det(A − λB)2worst residual, either kind1.5·10⁻¹⁵an eigenvalue is a ratioand a ratio has a direction, not a size
Fig. 28 And one with two, where inverting the leading coefficient is not available at all.
After reflector 1 of 5: the bulge in column 0A square matrix drawn as a grid with one nonzero subdiagonal. A small triangle of entries below that subdiagonal is highlighted, and moves down the matrix from panel to panel.0.461.90.14-0.0070.2-0.18-1.20.78-1.40.54-0.18-0.0791.40.440.120.870.19-0.077-0.49-0.320.932.30.90.0380000.982.11.500001.52.3after reflector 1the highlighted entries are the bulge — the only thing that is not Hessenbergentries below the subdiagonal3reflectors used1reflectors in a whole step5the shifts are never formed — only their sum and productand both of those are real
Fig. 29 The implicit step inside the factorisation every route ends in.
Where the missing eigenvalues went: an 6 × 6 quadratic with 2 of its masses removedA chain of 6 masses with the first 2 of them set to zero, so M is singular of rank 4. det Q(λ) is a polynomial of degree at most 12; interpolated exactly in BigInt rationals at 13 nodes it has degree 10, so 2 of the 12 eigenvalues are at infinity — the same object a descriptor pencil has, arriving here because a degree of freedom with no inertia is an algebraic constraint. The float route counts the singular values of M judged to be zero and reaches 2, backed by a gap of −∞ between consecutive singular values. One integer, two routes, and only the second of them is a decision.finite eigenvalues (degree of det Q)10at infinity (2n − degree)2at infinity, by the rank of M22n, if M were nonsingular12a degree, not a decisiondegree of det Q10at infinity2by the rank of M2singular-value gap−∞the count is a degreeand the other route is a judgement
Fig. 30 The count that decides whether a leading coefficient can be inverted.
A double root, approached: the pair separates like √ε and the accuracy fails like √uA chain of 12 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.6909√ε at every one of twelve decades — a spread of 1.013 in the constant — and the computed spectrum's worst relative error against the closed form rises as the pair closes, reaching 7.465·10⁻⁸ at the boundary itself. That is √u times a small constant: half the digits, on a problem where κ(K) is 67.83 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: 7.47·10⁻⁸separationcomputed errornothing is ill conditionedthe √ε constant0.69spread of it, twelve decades1error at the boundary7.5·10⁻⁸κ(K), unchanged throughout68half the digitsand no condition number to blame
Fig. 31 A loss every route shares, because it belongs to the problem.
Where a quadratic stops being hyperbolic, located by a Cholesky and by a sineThe overdamping margin min over modes of (βκ)² − 4κ, for a chain of 12 masses, against β. It reaches zero at β* = 1/sin(π/2(n+1)) = 8.29622981056, which is the closed form. Bisecting on a completely different question — does a Cholesky of −Q(μ) complete for some μ — gives 8.29622981056, agreeing to 13 digits. Neither route computes an eigenvalue. The marks below the axis are the largest imaginary part in the computed spectrum, which is zero to the rounding level above β* and not below it, so a third route agrees with the other two about where the boundary is.56.38277.765419.1481110.530811.913510⁻¹110¹10²stiffness damping βoverdamping marginβ* = 8.29623two routes to a boundaryclosed form β*8.3by certificate8.3difference3.5·10⁻¹³bisection steps44a factorisation that completesand a sine, agreeing to twelve digits
Fig. 32 And the boundary at which it arrives.
Six ways of computing one spectrum, at a change of units of 106Three linearisations — the first companion form, the second, and the symmetric member of the DL(Q) family — each reduced to a standard eigenvalue problem in two ways: by inverting the leading coefficient, whose eigenvalues are λ, and by inverting the trailing one, whose eigenvalues are 1/λ. All six have exactly the eigenvalues of the quadratic in exact arithmetic. Measured against the closed form at γ = 106, the best is first/leading at 5.473·10⁻⁶ and the worst second/trailing at 2.234·10⁻⁴, a spread of 40.81. Under the leading reduction the symmetric form and the first companion form are THE SAME MATRIX — −A₁⁻¹A₀ is [[−M⁻¹C, −M⁻¹K], [I, 0]] for both — so the distinction between them is one the inversion discards, and it survives only under the trailing one.first · leading5.47·10⁻⁶first · trailing4.62·10⁻⁵second · leading4.27·10⁻⁵second · trailing2.23·10⁻⁴symmetric · leading5.47·10⁻⁶symmetric · trailing4.05·10⁻⁵all six are the same algebrabest route5.5·10⁻⁶worst route2.2·10⁻⁴spread across the six41condition of the linearisation8.3·10¹²the spectra agreeand the arithmetic does not
Fig. 33 The table’s own stop, drawn again beside the ones on either side of it.
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 + 3K. 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 2 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 8.65·10⁻¹⁵, and the routes share nothing but the three coefficient matrices.-13-10.6516-8.30312-5.95468-3.60624-1.2578-13-8.55186-4.103730.3444124.792559.24069real partimaginary parttwo routes, one spectrumeigenvalues16rows8complex2against the closed form8.7·10⁻¹⁵n rows and 2n eigenvaluesso the eigenvectors are not a basis
Fig. 34 The overdamped spectrum all six routes are compared on.

What links here

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

Reads more easily once this is understood

Essays that name this one as worth reading first.

Shares its objects with

Essays that name at least two of the same things, and that neither author linked.

Named objects

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

Backward errorCompanion formCondition numberLinearisationMatrix polynomialQuadratic eigenvalue problemScalingSchur form