Numerical Linear Algebra
where the algebra stops being the arithmetic
Structure, and the solver that cannot see it

An equation whose unknown is a matrix

AX + XB = C is linear in X, so it has a coefficient matrix, and writing it down is the obvious thing to do. At n = 100 that matrix has a hundred million entries for a problem with ten thousand unknowns, and the algorithm everybody uses instead never forms it. Its conditioning is not the eigenvalue gap either, which is the number a reader is invited to consult.
22 min read 17 figures Where the algebra stops being the arithmeticWhose fault is itMeasured, not assertedTwo routes to a numberAssertions that reject

Worth reading first: The eigenvalues that are not there · The condition number is an amplifier.

Most of the equations in this collection have a vector for an unknown. Some have a matrix, and the best known of them is Sylvester’s:

AX + XB  =  C

with A, B and C given and X wanted. It appears wherever two operators have to be reconciled: the transformation that block-diagonalises a Schur form, the sensitivity of an invariant subspace, the Lyapunov equation of a stable linear system when B = Aᵀ, the observability and controllability Gramians of control theory, the coupling term in a domain decomposition.

The map X ↦ AX + XB is linear in X. So it has a matrix, and writing it down is the obvious first move.

sep(A, B) and the smallest |λᵢ(A) + μⱼ(B)| as the departure from normality grows, n = 6A and B are upper triangular with diagonals 1, 2, …, 6 and μ above them, so every eigenvalue sum is at least 2 at every point on the axis and the flat line is exact rather than nearly flat. sep — the smallest ‖AX + XB‖_F over ‖X‖_F = 1, and the reciprocal of the amplification a perturbation of C receives — starts equal to the gap at μ = 0, where the matrices are normal, and falls to 9.47·10⁻⁴ at μ = 8. The number a reader is invited to consult is the one that does not move.01234567810⁻⁵10⁻⁴10⁻³10⁻²10⁻¹110¹μ, the entry above the diagonalsep, and the eigenvalue gapmin |λᵢ + μⱼ|sep(A, B)the spectra never moveeigenvalue gap, throughout2sep at μ = 02sep at μ = 89.5·10⁻⁴amplification there1056solvability is the eigenvaluesand conditioning is not
Fig. 1 The conditioning of a Sylvester equation against the number a reader is invited to consult. The spectra do not move at any point on this axis; the flat line is exact. Drag the size.

The object the equation names

Stacking the columns of X into a vector — vec(X) — turns the equation into an ordinary linear system:

( I ⊗ A  +  Bᵀ ⊗ I ) vec(X)  =  vec(C)

where ⊗ is the Kronecker product. That is exact, it is what the equation means, and it is a system of n² equations in n² unknowns.

n unknowns in X entries in the coefficient matrix flop ratio
4 16 256 1×
20 400 160,000 89×
50 2,500 6,250,000 1,389×
100 10,000 100,000,000 11,111×

At n = 100 the object is 0.8 gigabytes of doubles for a problem whose answer is ten thousand numbers, and eliminating it costs (2/3)n⁶ = 6.7·10¹¹ operations.

The two 6×6 matrices of a Sylvester equation, and the 36×36 matrix it meansA and B are 6×6 with 11 nonzeros each. The coefficient matrix of the map X ↦ AX + XB is 36×36 with 96 nonzeros, drawn at the same scale so the ratio is visible rather than quoted. Bartels and Stewart's algorithm never forms it: two Schur reductions and a back-substitution give the same X to 1.7·10⁻¹⁶, and that X satisfies AX + XB = C to 2·10⁻¹⁶.A6×6B6×6I ⊗ A + Bᵀ ⊗ I36×36one equation, two objectsentries in A and B72entries in the coefficient matrix1296two routes, relative gap1.7·10⁻¹⁶‖AX + XB − C‖/‖C‖2·10⁻¹⁶the small squares are the problemand the large one is the notation
Fig. 2 The two 6×6 matrices of a Sylvester equation, and the 36×36 matrix it means, drawn at the same scale so the ratio is a picture rather than a number.
The coefficient matrix of AX + XB = C against the answer, to n = 100At n = 100 the unknown X has 10000 entries and the coefficient matrix of the linear map has 10⁸ — 0.80 gigabytes of doubles. Eliminating it costs 6.67·10¹¹ operations against the 6·10⁷ Bartels and Stewart's algorithm needs, a ratio of 11111. The Kronecker form is what the equation means and it is not a method.10¹10²10¹10³10⁵10⁷10⁹10¹¹nentries, and operationsKronecker flopsits entriesBartels–Stewartentries in Xat n = 100unknowns10⁴coefficient entries10⁸gigabytes of doubles0.8flop ratio1.1·10⁴the equation is linear in Xand nobody writes down its matrix
Fig. 3 The two sizes against n. Every curve here is a power of n, so the only thing that changes along the axis is how far apart they have got.

What is done instead

Bartels and Stewart’s algorithm, from 1972, and it is three steps.

Reduce. Compute the real Schur forms A = QRQᵀ and B = ZSZᵀ. Two orthogonal similarities, O(n³) each, and orthogonal means nothing is amplified.

Substitute. The equation becomes RY + YS = QᵀCZ with R and S triangular, and then the (i, j) entry reads

(rᵢᵢ + sⱼⱼ) yᵢⱼ  =  cᵢⱼ − Σ_{k>i} rᵢₖ yₖⱼ  −  Σ_{k<j} yᵢₖ sₖⱼ

so Y comes out one column at a time, each column bottom to top, with one division per entry. That is n³ operations and no coefficient matrix exists at any point.

Transform back. X = QYZᵀ, two products.

Total about 60n³ against (2/3)n⁶: at n = 20 a factor of 89, at n = 100 a factor of eleven thousand. And the memory is n² rather than n⁴, which is the difference between running and not running.

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. 4 The reduction the algorithm rests on, and the reason it is safe: every step is an orthogonal similarity, so the residual is bounded by construction.
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. 5 And the first half of it. The whole method is the observation that a triangular version of the equation is a back-substitution.

Two routes, and they agree

At n = 6 both routes can be run. The Kronecker solve eliminates a 36×36 system; Bartels–Stewart performs two Schur reductions and a back-substitution. They share no arithmetic beyond the matrices they start from.

The two answers differ by 6.6·10⁻¹⁷ relative, and the residual ‖AX + XB − C‖/‖C‖ is 1.1·10⁻¹⁶ for both. Neither route can confirm itself; the agreement is evidence, and the residual is the check neither could have done alone.

The two 3×3 matrices of a Sylvester equation, and the 9×9 matrix it meansA and B are 3×3 with 5 nonzeros each. The coefficient matrix of the map X ↦ AX + XB is 9×9 with 21 nonzeros, drawn at the same scale so the ratio is visible rather than quoted. Bartels and Stewart's algorithm never forms it: two Schur reductions and a back-substitution give the same X to 5.5·10⁻¹⁷, and that X satisfies AX + XB = C to 8.8·10⁻¹⁷.A3×3B3×3I ⊗ A + Bᵀ ⊗ I9×9one equation, two objectsentries in A and B18entries in the coefficient matrix81two routes, relative gap5.5·10⁻¹⁷‖AX + XB − C‖/‖C‖8.8·10⁻¹⁷the small squares are the problemand the large one is the notation
Fig. 6 At n = 3, where the coefficient matrix is a nine-by-nine and forming it is merely wasteful rather than impossible. The structure is visible: copies of A down the diagonal, entries of B spread across the blocks.

What the Kronecker form is still for

The n²×n² matrix is not a mistake, and it is worth saying what it is good for before dismissing it as a method.

It is the definition. Every statement about the equation’s solvability, uniqueness and conditioning is a statement about that matrix, and the ones below are derived from it here rather than quoted.

It is how the conditioning is measured. sep is its smallest singular value, and this site computes it by forming the thing — which is the same licence the inverse has for the componentwise condition number: form the object when its entries or its singular values are the answer, never as a step towards something else.

It is the right object at very small n. At n = 3 the coefficient matrix is a nine-by-nine and the Kronecker route is a perfectly reasonable way to solve a Sylvester equation, and easier to write. The ratio at n = 4 is 1.

So the argument is not that the Kronecker form should never be built. It is the same argument the whole phase makes: the object a question names is a fine object and usually the wrong intermediate, and the size at which it stops being a reasonable intermediate is around n = 8.

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. 7 The other case on this site where an object is formed because its entries are the answer. The rule is the same and so is the exception.

Solvability is the eigenvalues, and conditioning is not

The equation has a unique solution exactly when no eigenvalue of A is the negative of an eigenvalue of B. That is easy to see from the triangular form — the divisor at entry (i, j) is rᵢᵢ + sⱼⱼ — and it is a clean statement about the spectra.

So the natural guess is that the difficulty is measured by how close min|λᵢ + μⱼ| comes to zero.

It is not, and the construction that says so is deliberately blunt. A and B are both upper triangular with diagonals 1, 2, …, n and an entry μ above the diagonal. The eigenvalues are the diagonals and do not move as μ changes, so every eigenvalue sum is at least 2, exactly, at every point on the axis.

μ min |λᵢ + μⱼ| sep(A, B)
0 2 2.000
1 2 1.654
2 2 0.978
4 2 0.131
8 2 9.47·10⁻⁴

The quantity that governs the answer falls by a factor of 2,100 while the quantity a reader is invited to consult does not move at all.

sep(A, B) and the smallest |λᵢ(A) + μⱼ(B)| as the departure from normality grows, n = 8A and B are upper triangular with diagonals 1, 2, …, 8 and μ above them, so every eigenvalue sum is at least 2 at every point on the axis and the flat line is exact rather than nearly flat. sep — the smallest ‖AX + XB‖_F over ‖X‖_F = 1, and the reciprocal of the amplification a perturbation of C receives — starts equal to the gap at μ = 0, where the matrices are normal, and falls to 4.02·10⁻⁴ at μ = 8. The number a reader is invited to consult is the one that does not move.01234567810⁻⁵10⁻⁴10⁻³10⁻²10⁻¹110¹μ, the entry above the diagonalsep, and the eigenvalue gapmin |λᵢ + μⱼ|sep(A, B)the spectra never moveeigenvalue gap, throughout2sep at μ = 02sep at μ = 84·10⁻⁴amplification there2488solvability is the eigenvaluesand conditioning is not
Fig. 8 The same sweep at n = 8, where sep falls by a factor of 5,000. The collapse gets faster with the size because the superdiagonal has more places to act.

What sep is

sep(A, B) is defined as

sep(A, B)  =  min over ‖X‖_F = 1 of ‖AX + XB‖_F

which is the smallest singular value of the Kronecker matrix — the ordinary condition number of an ordinary linear map, arrived at without ever mentioning eigenvalues.

At μ = 0 the matrices are normal and sep is the eigenvalue gap, exactly: 2.000 against 2. That coincidence is why the guess is natural, and it is the same coincidence that makes the spectrum a sufficient description of a normal matrix and an inadequate one otherwise — the eigenvalues that are not there is the general statement and this is a consequence of it.

The amplification is attained rather than bounded. Perturbing C along the direction that achieves the minimum moves X by exactly 1/sep times as much: 0.5 at μ = 0 and 1,056 at μ = 8, measured, with a random perturbation reaching between 4 and 35 per cent of the worst case.

σ_min(zI − A) over the complex plane, for a bidiagonal 6×6 matrix with every eigenvalue at 0.8 and 2 above the diagonalA square of the complex plane shaded by how small σ_min(zI − A) is. Every eigenvalue is at 0.8, marked by the crosshair; the 10⁻³ level reaches out to 1.33, well outside the unit circle drawn through the picture. A perturbation of the matrix at the level of double-precision rounding can put an eigenvalue anywhere in that region.darker is smaller: 10⁻¹⁰, 10⁻⁶, 10⁻³, 10⁻¹one spectrum, two matricesspectral radius0.8reach of the 10⁻³ level1.3eigenvalues, all at0.8the circle is |z| = 1and every eigenvalue is well inside it
Fig. 9 What a non-normal matrix’s spectrum leaves out, drawn directly. sep is the same information in the form a matrix equation needs it.
Where 24 perturbations of size 10⁻⁸ put the eigenvalues of two 6×6 matrices with the same spectrumTwo clouds of eigenvalues in the complex plane. Both matrices have every eigenvalue at 0.8 exactly; both were perturbed by the same 24 random matrices of norm 10⁻⁸. The normal matrix's eigenvalues stay within 8.86·10⁻⁹ of where they were — the size of the perturbation — and the bidiagonal's spread out to 0.0701, a factor of 7.9·10⁶ further on the same data.100.1750.35real partimaginary part0.8bidiagonalnormalone spectrum, two matricesthe perturbation10⁻⁸normal, furthest moved8.9·10⁻⁹bidiagonal, furthest0.07ratio7.9·10⁶the two matrices have identical eigenvaluesand one of them says so under perturbation
Fig. 10 And the same statement made by perturbing rather than by computing a resolvent. The eigenvalues of the pair are exactly where they were and the equation built from them is not.

Why this matters outside the equation

sep is not a curiosity of one equation. It is the condition number of an invariant subspace, and that is where most people meet it without meeting its name.

If a Schur form is partitioned so that the leading block carries the eigenvalues of interest, the sensitivity of the corresponding invariant subspace to a perturbation of the matrix is 1/sep(T₁₁, −T₂₂) — which is exactly this quantity, for the two diagonal blocks. LAPACK’s xTRSEN returns it, and it is the number that says whether the subspace a computation has isolated is a real feature of the matrix or an artefact.

The site has drawn that subspace before and measured its stability directly. This essay is what the number in that measurement is.

The plane, and the two vectors inside it, at a gap of 10⁻⁶A circle representing the invariant plane, with two perpendicular radii for the eigenvectors and two more for the perturbed ones, rotated within the same circle.the invariant planesolid: beforedashed: aftersame perturbation, two questionsthe vectors turned, radians0.029the plane turned, radians7.6·10⁻⁸what left the plane5.6·10⁻⁸drawn in the unperturbed plane's own basisa radius is not determined; the circle is
Fig. 11 The invariant subspace whose conditioning sep governs. The plane can be well determined while its basis vectors are not, and sep is the quantity that decides the plane.
The width of a block of 4, at every stepThree step functions against the step count. A block of 4 independent columns keeps its width at every step. A block with one column copied from another loses that column before the first product with A and keeps 3 for the rest of the run. A block whose columns all lie in the span of two eigenvectors starts at 2 and ends at 0 after one step, having found the whole subspace.01234567891001234stepcolumns in the blockindependent columnsone column repeatedinside an invariant planewhen a block narrowsrandom: width at the end4repeated: width at the end3invariant: vectors found2the mechanism worksand nothing a random start does sets it off
Fig. 12 And the partition of a Schur form that produces the two blocks. Where the two spectra are close in the sep sense, the split is not meaningful however cleanly the eigenvalues separate.

Lyapunov, and the low-rank fact that makes it tractable

The symmetric special case, AX + XAᵀ = −C with A stable, is the Lyapunov equation, and it is the one that is actually solved at scale — in model reduction, in control, in the Gramians that decide which states of a system matter.

At scale it is solved neither by the Kronecker form nor by Bartels–Stewart, because n is 10⁵ and even n³ is too much. What is used instead rests on a fact about the answer: when C is low rank — which it is, because it is usually BBᵀ for a thin input matrix B — the singular values of X decay exponentially, so X can be represented as a low-rank factor and never formed as an n×n matrix.

That is this phase’s sentence again, one level up. The equation names an n×n unknown; the object worth computing is a thin factor of it; and the reason the substitution works is a theorem about the answer rather than a trick.

Error of the best rank-k approximation to a 12×12 matrixApproximation error against k on a logarithmic axis, with the measured error and the next singular value drawn as separate curves lying exactly on top of one another.123456789101110⁻¹⁰10⁻⁸10⁻⁶10⁻⁴10⁻²1rank k of the approximation‖A − Aₖ‖measured, 2-normσₖ₊₁, from theorymeasured, Frobeniusthe first two agreeto 2·10⁻⁸worst |‖A−Aₖ‖₂ − σₖ₊₁| / σₖ₊₁2·10⁻⁸worst Frobenius discrepancy2·10⁻⁸κ = 10⁹; 30 random rank-3 matrices, none closerthe error is σₖ₊₁
Fig. 13 The decay that makes it possible. A matrix whose singular values fall away is one that a thin factor represents, and the Lyapunov solution of a low-rank right-hand side is such a matrix.
What a rank-10 approximation can achieve, by spectrumA semi-logarithmic plot of singular value against index for three spectra — geometric decay, algebraic decay, and flat — with the rank-ten approximation error marked on each.1112131415110⁻³10⁻²10⁻¹1index jσⱼσ11geometric, 0.85ʲalgebraic, j^−1.00flatmeasured, and equal to σ₁₁rank-10 error, geometric0.2rank-10 error, algebraic0.091rank-10 error, flat160×60, spectrum chosen rather than the entriesthe matrix decides, not the method
Fig. 14 And the general shape of that decay, which is what every low-rank method on this site rests on.

The two-dimensional case, which is where it becomes unavoidable

Sylvester equations are usually met as an abstraction. They arrive concretely, in enormous numbers, from separable partial differential equations.

A Poisson problem on a rectangle, discretised with the same stencil in each direction, is exactly

A U + U Aᵀ  =  F

where A is the one-dimensional operator and U is the grid of unknowns as a matrix. Vectorising it — which is what every finite-difference code does — turns it into the Kronecker system, and the Kronecker system is what the site’s own two-dimensional model problem is.

That is worth noticing, because it means the site has been solving Sylvester equations since the depth phase without calling them that. A 31×31 grid is 961 unknowns; as a matrix equation it is two 31×31 operators. Every multigrid, every conjugate gradient, every algebraic hierarchy the iterative field built was operating on the Kronecker form of a matrix equation that has a direct solution in n³.

The direct route is not always better — the iterative methods scale to problems that are not separable, and separability is a strong assumption — but where it holds, a fast Poisson solver via the eigendecomposition of the one-dimensional operator is O(n² log n) for n² unknowns, and it is this essay’s algorithm with the Schur reduction replaced by a sine transform.

V-cycle convergence factor against problem size, ε = 0.001Three flat curves of convergence factor against the number of unknowns on a logarithmic axis. At small anisotropy one sits near one and the others near a tenth.10²10².⁵10³10³.⁵00.250.50.751unknownsresidual reduction per cyclepointy-linesemi-yy-line spread, 4× in size1.1·10⁻⁴point at the largest grid0.97y-line at the largest grid0.03715×15, 31×31, 63×63 interior gridsflat in the size, whatever the anisotropy
Fig. 15 The iterative machinery the site built for exactly this problem, which never noticed that its operator was a Kronecker sum.
A 4-level hierarchy on 15 points, and the coarse operator it impliesRows of dots, each row half the length of the one above it, with lines joining every coarse point to the three fine points it interpolates to.15731pointsfinest gridone unknown — the recursion bottoms out in a divisionthe coarse operator, two ways‖RA_hP − A_2h‖/‖A_2h‖10⁻¹⁸unknowns / finest grid1.7cycle cost, in fine sweeps12each coarse point reaches three fine ones½, 1, ½ — and the restriction is its transpose
Fig. 16 And the hierarchy underneath it. Separability is the structure a direct method exploits and an iterative one does not need, which is the trade in one sentence.

What the equation is for

Sylvester and Lyapunov equations are not a curiosity of numerical linear algebra; they are how several whole subjects state their central quantities, and it is worth naming them because the sizes involved are what make the Kronecker form a non-starter rather than merely wasteful.

Controllability and observability. A linear system ẋ = Ax + Bu has a controllability Gramian defined by AP + PAᵀ + BBᵀ = 0, and its observability Gramian by the transposed equation. Those two matrices answer the questions which states can be reached with bounded energy and which states leave a trace in the output, and their eigenvalues are the numbers a control engineer actually looks at. The system is stabilisable when P is positive definite; the directions where it is nearly singular are the ones no input can excite.

Model reduction. Balanced truncation — the standard way to replace a system of order 10⁵ by one of order 20 with a certified error bound — is exactly a simultaneous diagonalisation of those two Gramians. The Hankel singular values it truncates on are the square roots of the eigenvalues of PQ, and the error bound is twice the sum of the ones discarded. Two Lyapunov solves are the whole of the preprocessing, and at n = 10⁵ that is not an equation anybody writes a coefficient matrix for.

Block-diagonalising a Schur form. Given a Schur form partitioned into two diagonal blocks, the similarity that decouples them is the solution of T₁₁X − XT₂₂ = −T₁₂, which is a Sylvester equation and is where sep came from in the first place. That is what LAPACK’s xTRSEN does, and its returned condition estimate is 1/sep.

And the nonlinear cousin. The algebraic Riccati equation, AᵀX + XA − XBR⁻¹BᵀX + Q = 0, is the same shape with a quadratic term, and it is the equation behind every linear-quadratic regulator and every steady-state Kalman filter. It is solved by an eigenvalue method on a 2n×2n Hamiltonian matrix rather than by anything resembling elimination — which is the same move once more, since the object named in the equation is not the object the algorithm forms.

In every one of those the unknown is a matrix, the coefficient matrix of the linear map is never built, and the reduction that makes the problem tractable is a Schur form. The pattern is worth carrying rather than the equation: where an unknown has two indices, the thing to reduce is the operators acting on each of them separately.

The width of a block of 4, at every stepThree step functions against the step count. A block of 4 independent columns keeps its width at every step. A block with one column copied from another loses that column before the first product with A and keeps 3 for the rest of the run. A block whose columns all lie in the span of two eigenvectors starts at 2 and ends at 0 after one step, having found the whole subspace.01234567891001234stepcolumns in the blockindependent columnsone column repeatedinside an invariant planewhen a block narrowsrandom: width at the end4repeated: width at the end3invariant: vectors found2the mechanism worksand nothing a random start does sets it off
Fig. 17 The partition that produces the two blocks, and the equation whose solution decouples them. The condition number of that decoupling is sep between the blocks, and LAPACK returns it beside the subspace it isolated.

There is one more thing the Kronecker form is genuinely the right tool for, and it is the reason this essay could measure anything at all: it is how the conditioning is defined. sep is the smallest singular value of that matrix and there is no other way to say what it is. Forming a 36×36 object to learn a number about two 6×6 ones is exactly the licence the phase’s other essays give the componentwise condition number and the diagonal of a hat matrix — form it when its singular values or its entries are the answer, never as a step towards one.

What is worth carrying

An equation that is linear in a matrix has a coefficient matrix, and forming it is n² times the problem. At n = 100 it is 0.8 gigabytes for ten thousand unknowns, and the algorithm everybody uses is three orders of magnitude cheaper and never builds it.

Solvability is the eigenvalues and conditioning is not. λᵢ(A) + μⱼ(B) ≠ 0 decides whether there is an answer; sep decides how well determined it is, and the two agree only when the matrices are normal.

And sep is a familiar quantity in unfamiliar clothes. It is the condition number of an invariant subspace, LAPACK returns it, and it is the number that says whether a computed subspace is a feature of the matrix or an accident of the perturbation.

The last essay of this phase removes the matrix altogether: an operator with no entries.

Where this goes next

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 The form a real matrix can reach A real matrix with complex eigenvalues has no real triangular form, and the reason is one line — a real triangular matrix has a real diagonal, and a similarity does not move the spectrum. What it has instead is triangular except for one two-by-two block per conjugate pair, and the count is decided by the matrix rather than by where the iteration stopped. darker is smaller: 10⁻¹⁰, 10⁻⁶, 10⁻³, 10⁻¹one spectrum, two matricesspectral radius0.8reach of the 10⁻³ level1.3eigenvalues, all at0.8the circle is |z| = 1and every eigenvalue is well inside it The eigenvalues that are not there For a normal matrix the resolvent norm is exactly one over the distance to the nearest eigenvalue, so a picture of it carries nothing the spectrum did not. Move one entry above the diagonal and the region a perturbation of 10⁻⁸ can put an eigenvalue into stops being a disc and reaches out past the unit circle, while every eigenvalue stays at 0.8. 024681012141610⁻¹1index kλthe transformthe eigensolverC = F* Λ F is a factorisationworst relative disagreement2.6·10⁻¹⁵‖Cx − b‖/‖b‖ from the transform solve4.7·10⁻¹⁶imaginary part of a real spectrum1.2·10⁻¹⁶n = 16, and the whole matrix is 16 numberseigenvectors known in advance The matrix that is one row A circulant of size 16 is sixteen numbers, has no zero entry anywhere, and hands over its entire spectrum in closed form — the discrete Fourier transform of its first column, exactly. An eigensolver spends a sweep of Jacobi rotations over 256 entries arriving at the same answer, and agrees to 1.2·10⁻¹⁵.

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.

Back-substitutionCondition numberEigenvaluesFlop countKronecker productMatrix equationNon-normalitySchur form