Two matrices and one problem

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.

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 an answer that is known measures against a closed-form inverse rather than against a better float computation.

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 same shape of blindness a factorisation with nothing to pivot for finds in a definiteness test, where the routine’s only report about a property is whether it managed to finish. 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.

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. It is the same shape as two Gram–Schmidts read backwards: there two procedures that look identical behave differently, and here two that look different behave the same.

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.

Which is computable rather than merely nameable — Stewart’s characterisation makes it maxθ λmin(A cos θ + B sin θ), a scan over one angle — so it can be measured instead of invoked. Measured, at n = 8, with ‖B‖ = 1 at every stop by construction:

So the claim holds and holds more exactly than it was stated. The last column converges to 10, which is the largest eigenvalue of the pencil — this family’s Λ runs from 1 to 10 — so c is λmin(B) times the top of the spectrum, and with ‖B‖ fixed at one that is 10/κ(B). The fitted slope of c against κ(B) is −1.00. The reason offered for the finding is not an analogy; it is the same number.

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 — which is the exact answer to a nearby problem’s separation applied to a pencil rather than to a matrix: the algorithm’s contribution and the problem’s are different numbers, and only the second is at issue here.

What the number does not do #

It does not predict the error, and neither does the dashed line, and saying so is the difference between fifteen decades of agreement and fifteen decades of the right slope.

Both are bounds and both hold at every stop. The measured error sits between 7 and 600 times below u·κ(B), and between 6 and 114 times below the Crawford version — non-monotonically, which is the noise the section above warns about. The Crawford bound is uniformly the tighter of the two, by the same factor of about 6.6 throughout, which is what one expects of two bounds differing by a normalisation rather than by a mechanism.

So the honest form of the essay’s central comparison is: both reductions track u·κ(B) in slope over fifteen decades and neither comes within an order of it in level. That is enough to establish what the essay claims — that the two routes lose digits at the same rate, and that the rate belongs to the pencil — and it is not enough to use either bound to predict what a particular pencil will cost. That is the ordinary condition of a bound on this site and it is worth stating rather than hoping for better: the condition number is an amplifier measures the same slack for a linear system, and a small residual is not a small error is what happens to somebody who reads a bound as a forecast.

One narrower consequence of the Crawford measurement is worth extracting, because it changes what a practitioner should compute. κ(B) is available from a factorisation the reduction is performing anyway; the Crawford number costs a scan over one angle with a symmetric eigendecomposition at each step, which on this family is more work than the eigenproblem itself. Since the two bounds differ by a fixed factor of 6.6 here and both are loose by one to two orders, there is nothing to buy by computing the harder one — on this family. What would change that is a pencil where the two come apart, and the construction says exactly when: c is λmin(B) times the top of the spectrum only because A and B share their eigenvector structure through L. A pencil whose A is nearly singular in a direction where B is not has a small Crawford number and a modest κ(B), and there the harder number is the one that knows something. Nothing in this essay’s family can be that pencil, which is a limit of the family and not of the argument.

That limit is worth naming because it is the boundary between this essay and the next: a pencil whose first matrix carries the degeneracy is where an eigenvalue with no value begins, and the scaling that makes the two matrices comparable in the first place is the units the matrix is measured in’s question asked of a pair.

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.

Then what does the symmetric route buy #

The size the question is asked at turns out to matter, and the sweep says so before the section’s own measurement does.

The advantage is real at every size and it collapses as the size grows. Across n = 8, 10, 12, 16 and 20 the symmetric form’s asymmetry is smaller by at least 2.1·10⁶, 2.5·10⁴, 6,078, 1,331 and 643 times — three and a half orders of magnitude lost over a factor of 2.5 in n, while the accuracy slopes stay between 0.89 and 0.99 throughout and the asymmetry of B⁻¹A itself barely moves (1.33, 1.16, 1.08, 1.29, 1.29).

The advantage decays with n, and decays more slowly as it goes: a factor of 4.6 lost between twelve and sixteen, 2.1 between sixteen and twenty. It is heading somewhere above one rather than towards it.

So what the symmetry buys is a quantity that decays in n, and the decay is in the numerator rather than the denominator. B⁻¹A is about as unsymmetric at twenty as at eight; it is the symmetric form that is losing its exactness, because L⁻¹AL⁻ᵀ is formed by two triangular solves whose rounding accumulates with the size. The advantage is therefore largest exactly where it is least needed — on small pencils a general solver would handle anyway — and it is shrinking towards the sizes where the choice is expensive. Extrapolating the measured decay puts the two forms within a factor of ten of each other somewhere in the low hundreds, which is a size this field’s problems reach routinely.

Three things follow, and 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.

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.

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.

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.

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.

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.

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.

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.

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.