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

Two matrices and one problem

Ax = λBx is what a finite element model, a structural vibration and a constrained optimisation actually produce, and it is not the one-matrix problem with a change of variables. Everybody is told not to form B⁻¹A because it is not symmetric. That is true, the departure from symmetry is about one, and it is not what decides the accuracy.
28 min read 23 figures Whose fault is itMeasured, not assertedIdentical algebra, different arithmeticTwo routes to a number

Worth reading first: Symmetry is worth more than precision · A factorisation with nothing to pivot for · The condition number is an amplifier.

Every eigenvalue essay in this collection so far has had one matrix in it. That is not what a computation usually has.

Discretise a vibrating structure by finite elements and two matrices come out: a stiffness matrix and a mass matrix, and the natural frequencies are the λ with Kx = λMx. Linearise a constrained optimisation and a saddle-point pencil appears. Write down a small-signal model of a circuit, or a descriptor state-space system, or a stability problem in fluid dynamics, and the same shape turns up. The pencil A − λB is the object, not the matrix.

The first thing anybody does with it is get rid of the second matrix, and there are two ways.

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. 1 Both reductions, against κ(B), on a pencil whose eigenvalues are known exactly by construction. The dashed line is u·κ(B).

The pencil, and why the answer is known

The construction matters here more than usual, because the whole essay is a measurement of error and an error needs a truth.

Take any nonsingular L and set

B = L Lᵀ , A = L Λ Lᵀ

for a diagonal Λ. Then B is symmetric positive definite by construction, A is symmetric, and the generalised eigenvalues are exactly Λ — because L⁻¹AL⁻ᵀ is Λ, with nothing to compute. The eigenvectors are the columns of L⁻ᵀ. And κ(B) is κ(L)², which is the knob: choosing L with a prescribed condition number moves κ(B) across fifteen decades with the eigenvalues fixed.

So every error below is measured against a number that was decided before anything ran, in the way this site’s Hilbert essays measure against a closed-form inverse rather than against a better float computation.

The exact solution of a 13×13 Hilbert system beside the computed oneTwo columns of numbers: the exact answer, which is the integers one to thirteen, and the answer double-precision elimination returns, with the number of correct digits beside each.H13 x = b, b formed exactly so that x = (1, 2, …, 13)12345678910111213exact1.00002.00003.00133.97865.18864.993010.46720.048921.2657-2.575319.21478.906013.5113computed6.6 correct digits4.8 correct digits3.4 correct digits2.3 correct digits1.4 correct digit0.8 correct digitno correct digitsno correct digitsno correct digitsno correct digitsno correct digits0.6 correct digit1.4 correct digitbackward error2.2·10⁻¹⁷κ = 1.7·10¹⁸The algorithm solved a neighbouring problem perfectly. That problem's answer is this one.right-hand side built in BigInt rationalsthe truth is known
Fig. 2 The habit in its original form: a solution known in closed form, so an error is a measurement and not a disagreement between two computations that share a weakness.

The two reductions

Form B⁻¹A. Its eigenvalues are the generalised ones — B⁻¹Ax = λx is Ax = λBx rearranged — and it is a single matrix, so any eigensolver takes it. It is also not symmetric. A and B symmetric does not make B⁻¹A symmetric; that requires A and B to commute, which they do not.

Or factor B = LLᵀ and form C = L⁻¹AL⁻ᵀ. Two triangular solves, which cost what forming B⁻¹A costs, and C is symmetric — in exact arithmetic, exactly. Its eigenvalues are the same λ and its eigenvectors are L⁻¹ times the pencil’s.

The advice every textbook gives is to use the second, and the reason every textbook gives is the symmetry. The measurement below says the advice is right and the reason is not the one that matters.

The symmetry claim is true, and enormous

‖C − Cᵀ‖/‖C‖ for the two matrices, at every conditioning in the sweep:

  • B⁻¹A: between 0.82 and 1.08. The matrix departs from symmetry by about its own size. It is not nearly symmetric; it is a general matrix.
  • L⁻¹AL⁻ᵀ: 3.4·10⁻¹⁶ at the easy end — the unit roundoff — rising to 1.8·10⁻⁴ at the hard end.

That second number is worth pausing on. In exact arithmetic the reduced matrix is symmetric exactly. In floating point it is symmetric to about u·κ(B), and at κ(B) = 10¹⁶ that is 10⁻⁴, not 10⁻¹⁶.

No library ever notices, because a symmetric eigensolver reads one triangle and takes the other on trust. The asymmetry is real, it is computable, it is four orders of magnitude at the far end of this axis, and every routine that would be affected by it has been written not to look.

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. 3 The site’s habit applied to the same kind of claim elsewhere: orthogonality printed as a number rather than asserted as an adjective, so that “is orthogonal” becomes a measurement that can be wrong.

And it does not decide the accuracy

Here is the sweep, as relative error in the worst eigenvalue against κ(B):

  • The fitted slope of the B⁻¹A route is 0.92.
  • The fitted slope of the Cholesky route is 0.98.
  • The ratio between the two errors stays between 0.17 and 2.3 across fifteen decades — and it goes both ways. At some conditionings the unsymmetric route is the more accurate one.

Both are u·κ(B), and the dashed reference line in the hero is exactly that.

So the reduction that throws the symmetry away and the reduction that keeps it lose digits at the same rate. What sets the rate is κ(B), and κ(B) is a property of the pencil that no reduction of it can change.

That is the finding, and it has the shape this collection keeps finding: a piece of advice that is correct, universally given, and about a different quantity from the one it is usually justified by.

Why the conditioning belongs to the pencil

The reason is worth stating, because it turns the measurement from a surprise into an expectation.

An eigenvalue of the pencil is a root of det(A − λB). Perturb A and B by relative ε and the root moves by an amount governed by how close the pencil is to being singular — to having no eigenvalues at all — and that closeness is a property of the pair. The quantity that measures it is the Crawford number, the smallest value of √((xᵀAx)² + (xᵀBx)²) over unit x, and on the family here it is small exactly when κ(B) is large: B’s smallest eigenvalue is what makes both terms small on the same direction.

There is no reduction that escapes a property of the pair. There is only a reduction that does not add to it, and both of these are that.

How much a perturbation of the right-hand side is amplified, κ = 10⁶The cumulative distribution of the amplification factor over two hundred random perturbation directions, with the condition number marked as the upper limit.110¹10²10³10⁴10⁵10⁶10⁷00.250.50.751amplification of the input perturbationfraction of directions at or belowκ = 10·10⁵worst found 7.6·10⁵6×6, 200 directionsmedian reaches 0.29 of κ
Fig. 4 The relationship this rests on, from the error field: a condition number is an amplifier, and what it amplifies is whatever the algorithm’s own error happened to be.
Backward and forward error against the condition numberA log–log plot over twelve decades of condition number. The backward error is a flat line at ten to the minus sixteen; the forward error rises in proportion to the condition number.110²10⁴10⁶10⁸10¹⁰10¹²10¹⁴10⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹condition number κ(A)relative errorforward errorbackward errorpredicted: κ · u40×40, 20 seeds per κ; dashed is the worstthe problem worsens, not the method
Fig. 5 And the identity underneath the whole site: forward error is condition number times backward error, and only their product is what anybody sees.

The two errors cross over

The ratio between the two routes’ errors is not merely bounded; it changes sign of advantage along the axis. At κ(B) = 10² the Cholesky route is more than twice as accurate. At 10¹² it is three times worse. At 10¹⁶ it is worse again by a factor of six.

There is nothing to read into any individual crossing — these are single matrices at each stop and the noise in a worst-eigenvalue error is a factor of a few — and that is precisely the point of reporting the range rather than a mean. A difference that changes direction as the axis moves is a difference that is not there.

The comparison that is there is between either route and the dashed line. Both track u·κ(B) over fifteen decades, which is fifteen decades of agreement with a prediction, and neither departs from it in any systematic way.

Relative eigenvalue error of two reductions of Ax = λBx, against κ(B), n = 6The 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.07 — about its own size — and reducing through the Cholesky factor gives one that is nearer symmetric by a factor of at least 2.7·10⁶. That difference is real and it does not appear in the accuracy: the two error curves have fitted slopes of 1.04 and 1.01 against κ(B) and stay within a factor of 4.9 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⁻¹A1slope, via Cholesky1worst ratio between them4.9asymmetry of B⁻¹A1.1the symmetry claim is trueand it is not about the accuracy
Fig. 6 At n = 6 the same two slopes and the same dashed line. The constants move a little with the size; the relationship does not.
Relative eigenvalue error of two reductions of Ax = λBx, against κ(B), n = 20The 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.29 — about its own size — and reducing through the Cholesky factor gives one that is nearer symmetric by a factor of at least 643. That difference is real and it does not appear in the accuracy: the two error curves have fitted slopes of 0.92 and 0.96 against κ(B) and stay within a factor of 19.9 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.96worst ratio between them20asymmetry of B⁻¹A1.3the symmetry claim is trueand it is not about the accuracy
Fig. 7 And at twenty. What grows with n is the asymmetry of B⁻¹A — more room for A and B to fail to commute in — while the accuracy it is supposed to predict stays where it is.

Then what does the symmetric route buy

Three things, and they are worth having. None of them is the eigenvalues.

Half the work. A symmetric eigensolver is about half the operations of a general one and produces real eigenvalues without a complex arithmetic path. Given equal accuracy, that settles it.

Real answers, guaranteed. The eigenvalues of a symmetric-definite pencil are real. Handed B⁻¹A, a general eigensolver has no way to know that, so a pair that should be a double real eigenvalue can come back as a complex conjugate pair with a small imaginary part — and the imaginary part is rounding rather than physics. On a vibration problem that is a mode with a spurious damping term in it.

The spectrum in the plane, with the pair at 1 ± 0.10iA complex plane with the real axis marked. Open circles show the eigenvalues the matrix was built from and filled dots show the ones the algorithm returned; two of them sit symmetrically above and below the real axis.-5-3-1135-202real partimaginary partbuilt incomputedthe real line — where a real shift lives‖A − ZTZᵀ‖/‖A‖2.6·10⁻¹⁵‖ZᵀZ − I‖4.3·10⁻¹⁵worst eigenvalue error7.1·10⁻¹⁵departure from normality10⁻¹⁸5×5 real matrix, 1 conjugate pairthe answer is not on the axis
Fig. 8 Where a real matrix keeps a complex pair, from the essay on the real Schur form: a 2×2 block that cannot be split without leaving the reals.
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. 9 And the form a general solver reaches instead, which carries those pairs as blocks.

And B-orthogonal eigenvectors. This is the measurable one and it is the answer to “what is the symmetry actually worth”.

The eigenvectors of a symmetric-definite pencil are orthogonal in the inner product B defines: XᵀBX = I. That is the generalised version of the number this site is named after. Measured, the Cholesky route gives ‖XᵀBX − I‖ that is smaller than the inverse route’s at every conditioning in the sweep, by an average factor of 2.96 and never by more than ten.

A factor of three, not an order — and both of them degrade as u·κ(B) too, so even the orthogonality is governed by the same number.

Loss of orthogonality against condition number, in binary32A log–log plot of the norm of Q-transpose-Q minus the identity against condition number. Classical Gram–Schmidt rises steeply, modified Gram–Schmidt rises gently, and Householder is flat.110²10⁴10⁶10⁸10¹⁰10¹²10⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹condition number κ(A)‖QᵀQ − I‖classicalmodifiedHouseholderκ²uκu8×8, eight seeds per κ, binary32all three reconstruct A
Fig. 10 The one-matrix version of that measurement, which this collection opened with: orthogonality as a number rather than an adjective.
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 5 and 6: 5·10⁻⁵columns 4 and 6: 7.1·10⁻⁷columns 4 and 5: 10⁻⁸modifiedr[i][j] = qᵢ · v ↑ what is LEFT of itv = v − r[i][j] · qᵢthe three worst |qᵢ · qⱼ|:columns 1 and 6: 1.8·10⁻¹⁰columns 2 and 6: 7.3·10⁻¹¹columns 3 and 6: 7.5·10⁻¹²The two R factors agree to 1.2·10⁻¹¹ relative. The two Q factors do not.the 6×6 Hilbert matrixone word, eight orders
Fig. 11 And the pair of algorithms that made the point: identical algebra, different arithmetic, and the difference visible only if the number is printed.

The pattern this belongs to

It is worth putting this beside the two other places on this site where a well-known instruction turned out to be about a different quantity than its usual justification.

The normal equations are not to be used, and the reason given is that forming AᵀA squares the condition number. That reason is correct and the sharper statement is that below a computable value of ε the Cholesky factorisation of AᵀA fails outright, so the method does not degrade, it stops.

An inverse is not to be formed, and the reason given is accuracy. Measured, A⁻¹b and A\b differ by twelve orders of magnitude in backward error and by a factor of fifty in forward error — so the advice is right, and the size of the effect is nothing like what the usual justification suggests.

And here: do not form B⁻¹A, because it is not symmetric. The asymmetry is real and about one; the accuracy is the same either way; and the reason to prefer the other route is work, guaranteed real answers, and eigenvectors that come out B-orthogonal.

Three pieces of advice, all correct, all justified by a mechanism that is not the one doing the work. The common shape is that the stated reason is the one visible without running anything, and the real one needs a measurement.

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. 12 The first of the three, where the failure is a stop rather than a slide.
‖AX̂ − I‖ and ‖X̂A − I‖ for a computed Hilbert inverse, to n = 10Two residuals of the same computed inverse, each divided by ‖A‖‖X̂‖ so both are dimensionless. The exact inverse satisfies both to zero. The computed one satisfies ‖AX̂ − I‖ at 5.1·10⁻¹⁸ — the side its columns were solved along — and ‖X̂A − I‖ at 4.2·10⁻¹⁶, 82 times larger, at n = 10. Each column of X̂ is the exact solution of a slightly perturbed system, but a different perturbation for each column, so there is no single nearby matrix whose inverse X̂ is.4681010⁻²⁰10⁻¹⁹10⁻¹⁸10⁻¹⁷10⁻¹⁶10⁻¹⁵10⁻¹⁴10⁻¹³nresidual / (‖A‖ ‖X̂‖)‖X̂A − I‖‖AX̂ − I‖the two sides, comparedn = 4, ratio1.7n = 6, ratio6.4n = 8, ratio8.3n = 10, ratio82the exact inverse satisfies bothand the computed one satisfies the side it was computed along
Fig. 13 And the second, which this collection priced rather than repeated.

What the sweep does not do

Two honest limitations, both worth naming because they bound what the finding claims.

Scaling does not rescue this family. κ(B) can sometimes be reduced by a diagonal equilibration — this collection has a whole essay on the condition number being a choice of units — and here it cannot: scaling the rows and columns of B to unit diagonal changes κ(B) by less than twenty per cent at every stop. The ill-conditioning of this construction is not a units problem. On a real mass matrix, whose entries carry physical units that a modeller chose, it very often is, and equilibrating first is the cheapest thing available.

Two condition numbers of one 8×8 system, as its rows are put into different unitsFour curves against the spread of the row units, in decades. κ_∞ of the scaled matrix rises from 9.83 to 1.9·10¹⁰ while the componentwise condition number stays at 6.98 throughout — the same system, the same solution, and one of the two numbers is a fact about the units. Hilbert's two numbers are drawn flat beside them at 3.4·10¹⁰ and 1.2·10¹⁰: a matrix whose sensitivity no scaling repairs.0246810110²10⁴10⁶10⁸10¹⁰10¹²spread of the row units (decades)condition numberκ_∞(DA)cond(DA)Hilbert κ_∞Hilbert condone system, two numbersκ_∞ at no spread9.8κ_∞ at 10 decades1.9·10¹⁰cond, either end7Hilbert, equilibrated1.3·10¹⁰the solution is the same at every spreadand one of these curves knows it
Fig. 14 What the units are worth when they are the problem: the same matrix, the same answer, and a condition number that moves by decades.
Two perturbation bounds and the error that was measured, on a 8×8 matrix spread over 8 decades of unitsThree curves against the size of an entrywise relative perturbation. The normwise bound κ_∞·ε is a valid bound and sits 4·10⁷ times above the componentwise one cond(A, x)·ε, which is also a bound and is nearly attained by the worst of forty random perturbations at each size.10⁻¹⁴10⁻¹³10⁻¹²10⁻¹¹10⁻¹⁰10⁻⁹10⁻¹⁵10⁻¹²10⁻⁹10⁻⁶10⁻³1relative size of the entrywise perturbationrelative forward errorκ_∞ · εcond(A,x) · εmeasuredboth bounds holdκ_∞(A)1.9·10⁸cond(A, x)4.8ratio of the bounds4·10⁷both curves above the data are boundsand only one of them is a measurement
Fig. 15 And the condition number that does not move under scaling at all, which is the one to reach for when asking whether an ill-conditioning is real.

And the family has real eigenvalues by construction. A general pencil does not. Two matrices with no definiteness between them can have complex eigenvalues, arbitrarily sensitive ones, and the algorithm libraries actually run for that case is the QZ algorithm — a generalised Schur form reached by unitary transformations of both matrices at once, which never forms B⁻¹A or any factor of B. What this essay measures is the two reductions that are available when the pencil is symmetric-definite, which is the case a great many models produce and the case where the shortcuts exist.

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. 16 The one-matrix version of what QZ generalises: a reduction that changes no eigenvalue and makes every subsequent step cheaper.
Driving the subdiagonal to zero, with λ₄/λ₃ = 0.50A semi-logarithmic plot of the magnitude of the subdiagonal entry against iteration count for three shift strategies. The unshifted curve is a straight line; the two shifted curves plunge to the bottom of the plot within a few steps.0714212835424910⁻¹⁶10⁻¹³10⁻¹⁰10⁻⁷10⁻⁴10⁻¹iteration|subdiagonal entry|no shiftRayleighWilkinsontwo routes to one raterate, from the spectrum0.5rate, measured0.5iterations, none / Wilkinson6.7symmetric 4×4, spectrum 8, 4, 2, 1the dashed line is the prediction
Fig. 17 And the iteration it feeds, whose generalised form works on the pair rather than on either matrix.

What a modeller can do about κ(B)

Since the conditioning belongs to the pencil and not to the reduction, the only useful lever is the pencil, and there are two.

Shift the problem. The eigenvalues of (A − σB, B) are the eigenvalues of (A, B) minus σ, and the sensitivity of an individual eigenvalue depends on where it sits relative to the rest. A shift near the eigenvalues of interest, followed by an inverse iteration, is the standard machinery of every structural eigensolver, and what it buys is not a better κ(B) but a better-separated part of the spectrum. This collection has measured the same trade on the one-matrix problem: a gap decides an eigenvector, and a shift is how a gap is manufactured.

What a perturbation of 10⁻⁶ moves, against the gap it is applied acrossThree curves against the gap between two eigenvalues, both axes logarithmic. Two are flat; the third rises steeply as the gap closes and then levels off.10⁻⁹10⁻⁷10⁻⁵10⁻³10⁻⁸10⁻⁵10⁻²gap between the two eigenvalueshow far it movedthe eigenvectorsthe eigenvaluestheir plane‖E‖ / gapone perturbation, three answerseigenvalue shift, spread over the sweep1plane angle, spread over the sweep1eigenvector angle, spread1.6·10⁵the dashed line is Davis–Kahan's ‖E‖/gaptwo of the three never noticed
Fig. 18 The gap doing that work in the one-matrix case: the eigenvalue barely moves and the eigenvector moves a great deal, and the gap is what decides which.
The computed angle against the true one, two formulationsTwo curves of computed angle against true angle on logarithmic axes. One follows the diagonal all the way down; the other leaves it and flattens at a fixed level.10⁻¹⁴10⁻¹²10⁻¹⁰10⁻⁸10⁻⁶10⁻⁴10⁻²10⁻¹⁵10⁻¹³10⁻¹¹10⁻⁹10⁻⁷10⁻⁵10⁻³the true anglethe angle computed√(2u) = 1.49·10⁻⁸arcsine of ‖(I − QQᵀ)Q₂‖arccosine of σ(Q₁ᵀQ₂)two routes, one of which has a floorfloor of the arccosine route1.5·10⁻⁸√(2u)1.5·10⁻⁸worst overstatement1.5·10⁶angles returned as exactly zero3a plane tilted by a known anglethe flat part is the instrument, not the data
Fig. 19 And the floor underneath it, which is what a shift cannot get past.

Or change what B is. A mass matrix assembled with a lumped rather than a consistent scheme is diagonal, so κ(B) is a ratio of element masses and the whole problem above evaporates — at the cost of a discretisation that is a little less accurate. That is a modelling decision made long before anybody runs an eigensolver, it is made for other reasons entirely, and it is the single largest lever on the number this essay says decides everything.

Which is the last thing worth saying about it. The quantity that governs the accuracy of a generalised eigenvalue computation is fixed before the computation starts, by somebody who was thinking about elements rather than about conditioning, and no amount of care inside the solver moves it.

The refusal

The assertion is fed a pencil whose second matrix is negative definite.

It is what a sign error produces — a mass matrix assembled with the wrong orientation, a Hessian whose second block was formed as −H, a saddle-point system handed over with its blocks in the other order. The Cholesky factor the reduction is defined in terms of does not exist, and the failure is quiet: an unguarded factorisation takes the square root of a negative pivot, produces NaN, and every subsequent entry of the “factor” is NaN. The eigensolver then returns a spectrum of NaNs, which sorts without complaint and prints as a column of the letters N, a, N.

Refusing at the pivot is the difference between a message naming the matrix and a page of nothing.

Elimination with and without pivoting, ε = 10⁻¹⁴The same two-by-two system solved twice. With a row swap the answer is exact; without one the upper triangular factor contains an entry of order one over epsilon and the second component of the answer is wrong.[ ε 1 ; 1 1 ] x = [ 1 ; 2 ], exact answer (1.000000, 1.000000)with partial pivoting1101U after elimination1.0000001.000000computed xbackward error 0forward error 0without10⁻¹⁴10-10·10¹³U after elimination0.9992011.000000computed xbackward error 2·10⁻⁴forward error 5.7·10⁻⁴no error is raisedgrowth 10·10¹³
Fig. 20 The general form of that failure, from the elimination field: an algorithm that continues past the point where it should have stopped and returns an answer of the right shape.

One more measurement, and it is the uncomfortable one

The asymmetry of the reduced matrix rises with κ(B) — 3·10⁻¹⁶ at the easy end, 10⁻⁴ at the hard one — and the eigensolver it is handed to reads one triangle.

So at the hard end the solver is solving a symmetric problem that is not the problem it was given. It symmetrises by omission: whatever is in the upper triangle is the answer, and the lower triangle, which differs from it in the fourth digit, is discarded without being looked at.

Two things about that are worth separating.

It is not a bug. Symmetrising is the right thing to do, because the exact reduced matrix is symmetric and the asymmetry is entirely rounding. Averaging the two triangles, which is what this site’s own routine does before handing the matrix on, is a slightly better choice than taking one of them, and the difference between the two choices is smaller than the asymmetry itself.

And it is a measurement nobody makes. The size of the discarded difference is a free diagnostic: it is u·κ(B) with the constant included, computable in n² operations from a matrix that is already formed, and it says how much of the accuracy has already gone before the eigensolver starts. A routine that printed it would be telling the caller the one number this whole essay says decides the outcome, at a cost of one pass over the matrix.

That is the site’s rule applied to a place it had not been applied: a quantity that is computed and thrown away is a residual nobody printed.

Relative eigenvalue error of two reductions of Ax = λBx, against κ(B), n = 8The 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.33 — about its own size — and reducing through the Cholesky factor gives one that is nearer symmetric by a factor of at least 2.1·10⁶. That difference is real and it does not appear in the accuracy: the two error curves have fitted slopes of 0.89 and 0.90 against κ(B) and stay within a factor of 5.0 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.89slope, via Cholesky0.9worst ratio between them5asymmetry of B⁻¹A1.3the symmetry claim is trueand it is not about the accuracy
Fig. 21 At eight, where the constants differ and the two slopes do not.
The growth factor of a 10×10 elimination, against the condition number of the matrixThree curves against κ. Cholesky's growth factor is exactly 1 at every condition number drawn — the elimination never produces an entry larger than the matrix already had. Partial pivoting on the same matrices reaches the same growth, and reaches it by making up to 5 row interchanges, each of which destroys the symmetry that was the reason to use a symmetric factorisation. The bound the general theory allows, 2^9 = 512, is drawn above them both.10¹10³10⁵10⁷10⁹10¹¹110¹10²10³10⁴condition number of the matrixgrowth factorbound 2^9partial pivotingCholeskyno pivot to gain fromCholesky growth, every κ1Cholesky interchanges0partial pivoting, worst5the bound, 2^9512both eliminations reach the same growthand only one of them had to swap to get there
Fig. 22 And what the factorisation the reduction rests on costs, from the elimination field.
Relative eigenvalue error of two reductions of Ax = λBx, against κ(B), n = 16The 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.29 — about its own size — and reducing through the Cholesky factor gives one that is nearer symmetric by a factor of at least 1331. That difference is real and it does not appear in the accuracy: the two error curves have fitted slopes of 0.93 and 0.99 against κ(B) and stay within a factor of 54.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.93slope, via Cholesky0.99worst ratio between them54asymmetry of B⁻¹A1.3the symmetry claim is trueand it is not about the accuracy
Fig. 23 At sixteen the ratio between the two routes reaches its widest in this family — fifty-four one way and a half the other, which is what a difference that is not there looks like when it is swept.

What is next

The pencil here has a nonsingular B, so it has n finite eigenvalues and the reductions exist. The next essay is about what happens when B is singular — where some of the eigenvalues are infinite, that is not a degeneracy, and the representation that survives is a pair of numbers rather than one.

Where this goes next

10⁻¹⁶10⁻¹³10⁻¹⁰10⁻⁷10⁻⁴10⁻¹10⁻¹⁶10⁻¹³10⁻¹⁰10⁻⁷10⁻⁴10⁻¹size of the perturbation ‖δA‖how far the eigenvalues moveJordan block, ε^(1/8)symmetric, ≤ ‖δA‖rounding error alone moves it to 10⁻²six seeds per symmetric point; Jordan is closed formsymmetry beats precision Symmetry is worth more than precision A symmetric matrix gives up its eigenvalues to full accuracy however ill-conditioned it is. An unsymmetric one can move them by the eighth root of a perturbation, so the rounding involved in merely storing the matrix shifts the spectrum by a hundredth. 10¹10³10⁵10⁷10⁹10¹¹110¹10²10³10⁴condition number of the matrixgrowth factorbound 2^11partial pivotingCholeskyno pivot to gain fromCholesky growth, every κ1Cholesky interchanges0partial pivoting, worst9the bound, 2^112048both eliminations reach the same growthand only one of them had to swap to get there A factorisation with nothing to pivot for Cholesky's growth factor is not bounded by one. It is equal to one, at every size and every condition number, and the two-line reason is why the algorithm needs no pivoting at all — not "usually gets away without it". Its only failure is the square root of a non-positive number, which is exactly the test for definiteness, and in floating point that test moves with the precision. 0246810110²10⁴10⁶10⁸10¹⁰10¹²spread of the row units (decades)condition numberκ_∞(DA)cond(DA)Hilbert κ_∞Hilbert condone system, two numbersκ_∞ at no spread9.8κ_∞ at 10 decades1.9·10¹⁰cond, either end7Hilbert, equilibrated1.3·10¹⁰the solution is the same at every spreadand one of these curves knows it The units the matrix is measured in One linear system, written twice. The rows of the second are the rows of the first in different units, the solution is identical to the last bit, and the condition number has moved by eight orders of magnitude. One of those two numbers is a fact about the problem and the other is a fact about the notation.

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 errorCholeskyCondition numberEigenvalue condition numberEigenvaluesGeneralised eigenvalue problemMatrix pencilOrthogonalitySymmetry