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

A function of a matrix is not a function of its entries

Everybody learns that f(A) means diagonalise, apply f to the eigenvalues, undiagonalise. That is a definition, not a method. On a matrix seven picometres from a defective one — with exact eigenvalues and eigenvectors from a closed form — the definition returns an answer wrong by sixty-five orders of magnitude, and a method that never mentions an eigenvalue returns the right one.
21 min read 15 figures Where the algebra stops being the arithmeticWhose fault is itTwo routes to a numberExact ground truthAssertions that reject

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

The first thing anybody is told about a function of a matrix is that it is not the elementwise version. exp(A) is not the matrix of exp(aᵢⱼ); the exponential of the identity is not the matrix of e’s. Whatever f(A) is, it is not f applied to the entries.

The second thing is what it is, and it comes in the form of a recipe:

A = VΛV⁻¹   ⟹   f(A) = V f(Λ) V⁻¹

Diagonalise, apply f to each eigenvalue, undiagonalise. It is correct wherever A is diagonalisable, which is almost everywhere, and it looks like a method because it is written as three steps.

It is a definition. This essay is about the distance between those two words, and the distance is sixty-five orders of magnitude.

Two routes to e^A as a 8×8 bidiagonal matrix walks up to a defective oneA_δ is upper bidiagonal with diagonal λ, λ+δ, …, so its eigenvalues are exact and its eigenvectors have a closed form. As δ falls the eigendecomposition route's relative error rises to 2.92·10⁶⁵ while scaling and squaring's falls to 4.09·10⁻¹², which is the distance ‖A_δ − A₀‖ = 7·10⁻¹² to the matrix whose exponential is known in closed form. The two agree at the right of the axis, where the matrix really is far from the limit; the divergence is κ(V), which reaches 3.3·10⁸².10⁻¹²10⁻¹⁰10⁻⁸10⁻⁶10⁻⁴10⁻²10⁻¹⁴10⁻¹10¹²10²⁵10³⁸10⁵¹10⁶⁴δ, the gap between consecutive eigenvaluesrelative error in e^Aan answer with no correct digitsV f(Λ) V⁻¹scaling and squaring‖A_δ − A₀‖exact eigenvalues throughoutκ(V) at the smallest δ3.3·10⁸²eigen route2.9·10⁶⁵scaling and squaring4.1·10⁻¹²distance to the limit7·10⁻¹²the eigenvalues are the diagonaland they are exact at every stop
Fig. 1 Two routes to e^A as a bidiagonal matrix walks up to a defective one. The eigenvalues are exact, the eigenvectors come from a closed form, and one route’s answer leaves the page while the other’s tracks the distance to the limit. Drag the entry above the diagonal.

The family, which is built to leave no excuse

The obvious objection to any experiment of this kind is that the eigensolver was bad. So the experiment is arranged so that no eigensolver runs.

A_δ is upper bidiagonal: diagonal λ, λ+δ, λ+2δ, …, and μ above it. Two consequences:

Its eigenvalues are its diagonal, exactly. A triangular matrix’s eigenvalues are on its diagonal — that is not a computation, it is the definition of the characteristic polynomial of a triangular matrix — and the entries are written down.

Its eigenvectors have a closed form. Back-substitution gives, for the eigenvalue λ + kδ,

vᵢ  =  μ^(k−i) / ((k−i)! δ^(k−i)),   i ≤ k

which is a formula, not an iteration.

So the route is handed exact eigenvalues and exact eigenvectors, and it still has to invert V. At δ = 10⁻¹² and n = 8, κ(V) is 3.3·10⁸².

And at δ = 0 the matrix is λI + μS with S the shift, which is not diagonalisable at all — one eigenvalue of multiplicity n, one eigenvector, no V — and whose exponential is known in closed form:

e^{λI + μS}  =  e^λ (I + μS + μ²S²/2! + … + μ^(n−1)S^(n−1)/(n−1)!)

a finite sum, because S is nilpotent. Every entry is e^λ μᵏ/k!.

The family therefore has a theorem at its limit and a route to it that has none.

Three perturbations of 10⁻¹⁴, three eigenvectors, one planeA circle with three radii at widely different angles, one for each perturbation, and a table of their residuals showing all three satisfy the eigenvalue equation.every radius is an eigenvectorseed 11: 2.0°seed 23: 60.5°seed 41: 18.2°none of them is wrongspread of the returned angles1spread of the returned planes1.3·10⁻¹⁵worst ‖Ax − λx‖5·10⁻¹⁵an exactly repeated eigenvalueno correct answer to fail to return
Fig. 2 What happens to an eigenvector as two eigenvalues merge. The direction is a perfectly good vector at every δ and it has no limit, which is the whole of the difficulty stated geometrically.

What each route returns

δ ‖A_δ − A₀‖ eigen route scaling and squaring
10⁻¹ 0.70 0.545 0.545
10⁻² 0.070 0.042 0.042
10⁻³ 0.0070 4.4·10² 0.0041
10⁻⁴ 7.0·10⁻⁴ 7.2·10⁸ 4.1·10⁻⁴
10⁻⁶ 7.0·10⁻⁶ 1.6·10²³ 4.1·10⁻⁶
10⁻⁸ 7.0·10⁻⁸ 5.5·10³⁷ 4.1·10⁻⁸
10⁻¹² 7.0·10⁻¹² 2.9·10⁶⁵ 4.1·10⁻¹²

Two things to read out of it.

The right-hand column equals the second column. Scaling and squaring returns the exponential of the matrix it was handed, to full precision; the residual error is the distance from A_δ to A₀, which is not an error at all — it is the two matrices being different matrices.

The first two rows are the control. At δ = 0.1 and δ = 0.01 the two routes agree, because the matrix really is that far from the limit and κ(V) is only 10⁵ and 10¹². The separation is not an artefact of the comparison; it appears exactly where κ(V) does.

Two routes to e^A as a 8×8 bidiagonal matrix walks up to a defective oneA_δ is upper bidiagonal with diagonal λ, λ+δ, …, so its eigenvalues are exact and its eigenvectors have a closed form. As δ falls the eigendecomposition route's relative error rises to 4.57·10⁶⁸ while scaling and squaring's falls to 3.76·10⁻¹², which is the distance ‖A_δ − A₀‖ = 7·10⁻¹² to the matrix whose exponential is known in closed form. The two agree at the right of the axis, where the matrix really is far from the limit; the divergence is κ(V), which reaches 5.4·10⁸⁶.10⁻¹²10⁻¹⁰10⁻⁸10⁻⁶10⁻⁴10⁻²10⁻¹⁴110¹⁴10²⁸10⁴²10⁵⁶δ, the gap between consecutive eigenvaluesrelative error in e^Aan answer with no correct digitsV f(Λ) V⁻¹scaling and squaring‖A_δ − A₀‖exact eigenvalues throughoutκ(V) at the smallest δ5.4·10⁸⁶eigen route4.6·10⁶⁸scaling and squaring3.8·10⁻¹²distance to the limit7·10⁻¹²the eigenvalues are the diagonaland they are exact at every stop
Fig. 3 The same family with a larger entry above the diagonal, which makes κ(V) grow faster and the separation open sooner. Nothing about the eigenvalues changes.

The one-sentence explanation

f(A) is a continuous function of A. V is not.

The exponential of a matrix is an entire function of its entries — it is a convergent power series in them — so a perturbation of 7·10⁻¹² in A moves e^A by about 7·10⁻¹². That is a theorem and the right-hand column of the table is it, measured.

The eigenvector matrix has no such property. As two eigenvalues approach each other their eigenvectors approach each other, V approaches singularity, and V⁻¹ blows up. At a defective matrix V does not exist at all. So a formula written in V is a formula in a quantity that has no limit where the answer does.

That is the sentence, and it applies to every f, not just the exponential. The square root, the logarithm, the sign, the sine — all of them have the same definition, all of them are continuous in A where they are defined, and all of them have the same defect if computed through V.

The condition number of one eigenvalue, at an off-diagonal entry of 100A rising curve against the size of an off-diagonal entry, on logarithmic axes, with a flat line at one below it and a vertical marker at the current value.110¹10²10³110¹10²10³off-diagonal entry ccondition number of the eigenvalue√(1 + c²)decoupled: 1measuredthree routes, one number‖A − ZTZᵀ‖/‖A‖1.7·10⁻¹⁵closed form100computed 1/|yᵀx|100worst measured movement46four eigenvalues, two conditioning numbersthe symmetric case has one, and it is 1
Fig. 4 The conditioning of an eigenvalue, which is 1/|yᵀx| for its left and right eigenvectors. That number diverges at a defective matrix, and it is the same divergence the route above inherits.
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. 5 And the gap that governs it. The eigenvector’s sensitivity is set by the distance to the next eigenvalue, so a small gap is not merely awkward — it is the quantity in a denominator.

What a method looks like instead

The method that gets the right answer never mentions an eigenvalue.

Scaling and squaring. Choose s so that ‖A/2ˢ‖ is around a half. Compute a rational (Padé) approximant to e^{A/2ˢ}, which is a polynomial in A divided by another polynomial in A — evaluated as one linear solve. Then square s times:

e^A  =  (e^{A/2ˢ})^(2ˢ)

Every operation is a matrix product or a linear solve. No spectrum is computed, no eigenvector is formed, and nothing in it can notice whether A is diagonalisable.

That is not a coincidence, and it is the design principle: a computation whose only operations are polynomials in A and solves with A automatically respects the continuity that f(A) has and V does not. The site’s whole family of Krylov methods has the same property for the same reason.

Relative error of scaling and squaring against the number of squarings, μ = 4Each point computes e^{A/2ˢ} by a Padé approximant and squares it s times, against a closed-form exponential. At s = 0 the approximant is being used far outside its range and the error is 6·10⁻¹²; the best is 6.88·10⁻¹⁶ at s = 2; at s = 15 it has risen again to 3.16·10⁻¹². The right-hand rise is the squaring phase amplifying its own rounding, which is why a library chooses s from ‖A‖ rather than taking as many as it can afford.0246810121410⁻¹⁷10⁻¹⁶10⁻¹⁵10⁻¹⁴10⁻¹³10⁻¹²10⁻¹¹10⁻¹⁰10⁻⁹10⁻⁸squarings, srelative error against the closed formbest: s = 2the approximant is out of rangethe squarings amplifyboth sides costerror at s = 06·10⁻¹²best, at s = 26.9·10⁻¹⁶error at s = 153.2·10⁻¹²‖A‖ / 2^s at the best1.1the squarings are exact in the algebraand each one doubles the rounding
Fig. 6 The number of squarings, which is the one parameter the method has, and the U-shaped curve it traces against a closed-form answer. The next essay is about both sides of that U.
‖Aᵏ/k!‖ against ‖e^A‖ for Moler and Van Loan's 2×2, scaled by 1The terms of the Taylor series rise to 1.4·10⁷ at k = 16 and fall back; the sum they produce has norm 2.58, drawn as the flat line. The ratio is 5.42·10⁶, so that many digits are cancelled away before truncation is a question, and the computed exponential has a relative error of 5.2·10⁻⁹ after 77 terms. Scaling and squaring on the same matrix returns 8.5·10⁻¹⁴.01020304050607010⁻⁶10⁻³110³10⁶10⁹k‖Aᵏ/k!‖‖e^A‖ = 2.6largest term 1.4·10⁷what the series throws awaylargest term1.4·10⁷‖e^A‖2.6digits cancelled away5.4·10⁶error after the sum5.2·10⁻⁹every term is computed correctlyand the sum has lost seven digits
Fig. 7 And the method that also mentions no eigenvalue and does not work: the Taylor series, whose largest term is 5.4 million times its own sum. Being polynomial in A is necessary and not sufficient.

The Jordan form is not the fix

The natural repair, on being shown the above, is to reach for the Jordan form: if the trouble is that A is nearly defective, use the decomposition that handles defective matrices.

It is worse, and the reason is instructive. The Jordan form is a discontinuous function of the matrix: A₀ has one block of size 8, and A_δ for any δ ≠ 0 has eight blocks of size 1. An arbitrarily small perturbation changes the answer completely, and a computation cannot distinguish δ = 0 from δ = 10⁻²⁰.

The Schur form is the one that works — every matrix has one, it is computed by orthogonal similarities so nothing amplifies, and f of a triangular matrix can be got by a recurrence. That recurrence has its own difficulty when two diagonal entries are close, which is handled by grouping them into blocks, and the whole business is the standard funm algorithm. It is a real method and it is not the definition.

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. 8 The form that always exists and is reached by orthogonal transformations. Every serious algorithm for a general matrix function goes through this and not through V.
Householder reduction to Hessenberg form, on a symmetric 7×7Two matrices side by side. The first is full; the second has zeros everywhere below the first subdiagonal, and being symmetric is tridiagonal.2.7-1.10.0332.4-4.21.12.6-1.15.21.70.27-0.66-0.36-2.50.0331.7-0.16-1.9-1.4-1.1-0.682.40.27-1.93.31-1.2-2.6-4.2-0.66-1.41-4.1-1.8-0.0911.1-0.36-1.1-1.2-1.81.90.652.6-2.5-0.68-2.6-0.0910.653.6A, symmetric→2.75.7000005.7-1.63.9000003.93.34.6000004.60.72-2.500000-2.51.61.9000001.93.71.5000001.52H = QᵀAQ, tridiagonal‖A − QHQᵀ‖/‖A‖10⁻¹⁵below the subdiagonal0worst eigenvalue movement1.4·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. 9 And the reduction beneath it, whose residual is bounded because every step is a reflection. That boundedness is the property V lacks.

Where this bites in practice

A nearly defective matrix sounds like a laboratory construction. It is the ordinary case in three places.

Repeated roots from symmetry. A structure with a symmetry has degenerate modes, and a small asymmetry — a manufacturing tolerance, a mesh that is not quite symmetric — splits them by a tiny amount. The eigenvectors are then determined by the perturbation rather than by the structure.

Multiple time constants in a linear system. ẋ = Ax with two decay rates that happen to be close. The transition matrix is e^{At}, and computing it through the eigendecomposition puts the difference of two nearly equal rates in a denominator.

Companion and controllability matrices. A companion matrix built from a polynomial with a near-double root is nearly defective by construction, and control theory forms them constantly.

In all three the eigenvalues themselves may be perfectly well determined — a symmetric matrix’s eigenvalues are always well conditioned — while the eigenvectors are not, and it is the vectors the formula needs.

Separating two eigenvalues, against how close they areIterations against the gap between the two largest eigenvalues, on a logarithmic gap axis. The single-vector method needs 17, 18, 20, 22 steps as the gap closes through four decades. The block of two needs 11 at every gap — and spends 22 products with A doing it, which is no less arithmetic. What it saves is synchronisations.10⁻⁴10⁻³10⁻²10⁻¹0510152025gap between the two eigenvaluesiterationsone vectorblock iterationsblock productsiterations, not arithmeticsingle-vector steps at 0.117single-vector steps at 0.000122block iterations, every gap11the same products with Aand half the synchronisations
Fig. 10 Which is an earlier phase’s result, and the reason this essay can lean on it: an eigenvalue and its eigenvector have different condition numbers, and only one of them is governed by the gap.
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. 11 And the non-normal case, where a rounding-sized perturbation moves the eigenvalues a long way. A matrix function computed through them inherits every bit of that movement.

What the definition is still for

None of this makes the eigendecomposition useless, and it is worth being precise about what it is still exactly right for, because the list is longer than the argument above suggests.

Reasoning. f(A)'s eigenvalues are f of A’s eigenvalues; f(A) commutes with A; ‖f(A)‖ ≥ max|f(λ)|. All of those come straight from the definition and none of them is a computation.

Symmetric and normal matrices. Here V is orthogonal, κ(V) = 1, and the definition is a method — the best one there is. The whole objection above is about κ(V), and for a symmetric matrix that number is one by construction. A symmetric matrix function should be computed by diagonalising, always.

Deciding what a function of a matrix even means for a defective one. The Jordan definition — f applied to a block involves f, f′, f″, … — is how the object is defined where there is no basis of eigenvectors, and the closed form this essay measures against is exactly that definition evaluated by hand.

So the split is clean: normal matrices, diagonalise; everything else, stay polynomial in A. That is one of the sharper practical rules the site has, and it falls out of a single number.

How far a perturbation of size ε moves an eigenvalue, 10×10A log–log plot of eigenvalue movement against perturbation size. The symmetric case lies on a line of slope one; the non-symmetric case lies on a line of slope one eighth, and at a perturbation of ten to the minus sixteen it has already moved by a hundredth.10⁻¹⁶10⁻¹³10⁻¹⁰10⁻⁷10⁻⁴10⁻¹10⁻¹⁶10⁻¹³10⁻¹⁰10⁻⁷10⁻⁴10⁻¹size of the perturbation ‖δA‖how far the eigenvalues moveJordan block, ε^(1/10)symmetric, ≤ ‖δA‖rounding error alone moves it to 10⁻²six seeds per symmetric point; Jordan is closed formsymmetry beats precision
Fig. 12 The site’s earlier measurement of what symmetry buys: a symmetric matrix gives up its eigenvalues without complaint, and it gives up its eigenvectors in an orthogonal basis. That is the other half of the rule above, and it is one number — κ(V) = 1.

The exponential is where the stakes are

Of all the matrix functions, e^{At} is the one that is actually computed, because it is the solution operator of ẋ = Ax. Every linear system of ordinary differential equations, every Markov chain in continuous time, every exponential integrator, every state transition in a Kalman filter over an interval — all of them are asking for e^{At} or for its action on a vector.

And the standard failure has a name. Moler and Van Loan’s Nineteen Dubious Ways to Compute the Exponential of a Matrix is fifty years old, was reissued twenty-five years later with the count unchanged, and its conclusion was that no method is reliable for all classes of matrices. Two of the nineteen have survived as defaults: scaling and squaring with a Padé approximant, for the matrix; and Krylov, for the action on a vector.

Both of them are polynomials in A. That is not a coincidence either.

‖e^{tA}‖₂ for a 10×10 matrix with every eigenvalue at −1 and 4 above the diagonalThe norm rises to 3.68·10⁴ at t = 9.00 — which is (n − 1)/|λ| = 9.0 — before turning over and decaying to 1.8·10⁻¹⁸. The spectral prediction e^{λt}, drawn beneath it, is 1.2·10⁻⁴ at the peak: wrong by a factor of 3·10⁸. Every value is from the closed form, so the hump is a property of the matrix and not of a method.0102030405060708010⁻²⁰10⁻¹⁶10⁻¹²10⁻⁸10⁻⁴110⁴t‖e^{tA}‖₂‖e^{tA}‖ = 1‖e^{tA}‖e^{λt}every eigenvalue is −1peak norm3.7·10⁴at t =9e^{λt} at the peak1.2·10⁻⁴‖e^{tA}‖ at t = 801.8·10⁻¹⁸the spectrum says it decaysand it does, eventually
Fig. 13 And the reason e^{At} is worth computing carefully rather than estimating from the spectrum: a matrix with every eigenvalue at −1 whose exponential rises by four orders of magnitude before it decays.
‖Aᵏ‖ for a 6×6 matrix whose spectral radius is 0.8, with 2 above the diagonalTwo curves against the power. The norm of Aᵏ rises to 19800 at step 24 before turning over and decaying to 2.5·10⁻⁵ by step 160; ρᵏ = 0.8ᵏ, drawn beside it, falls from the start. The two horizontal lines are the Kreiss bracket: the peak is at least K = 6757 and at most e·n·K = 1.1·10⁵, both computed from the resolvent norms outside the unit circle and not from the powers at all.027548110813510⁻⁶10⁻⁴10⁻²110²10⁴10⁶power‖Aᵏ‖Kreiss constant 6760e · n · K‖Aᵏ‖ρᵏtwo routes to one peakspectral radius0.8peak of ‖Aᵏ‖2·10⁴Kreiss constant6757e · n · K1.1·10⁵everything here decays in the endand one of these curves says how much first
Fig. 14 The discrete-time twin, from an earlier phase. The powers of a matrix with spectral radius 0.8 grow by a factor of 19,800 first, and the spectrum says nothing about it.

The other functions, and where each one stops existing

The exponential is the easy case in one respect that is worth being explicit about: it is entire, so e^A exists for every matrix and there is nothing to choose. The functions that follow it are not, and what they lose is worth separating from what this essay has been about.

The square root. Every matrix with no eigenvalue on the negative real axis has a principal square root — the one whose eigenvalues have positive real part — and it is unique. But a matrix with n distinct eigenvalues has 2ⁿ square roots, one per choice of sign, and if any eigenvalue is repeated there are infinitely many. So “the square root of A” is a convention, and a computation has to say which one it is producing. The identity matrix’s square roots include every reflection.

The logarithm. Same structure, worse: log A exists when A has no eigenvalue on the closed negative real axis, the principal branch is the one whose eigenvalues have imaginary part in (−π, π], and a real matrix with negative real eigenvalues has no real logarithm at all. That last case is not exotic — a rotation by π has one.

The sign function. sign(A) is defined for any matrix with no eigenvalue on the imaginary axis, and it is not continuous there: a matrix one rounding error either side of it returns a completely different answer. It is also the one whose definition is most obviously a computation rather than a formula, since sign(A) = A(A²)^(−1/2) is a composition of two other functions.

All three inherit this essay’s objection in full, and all three inherit something extra: where the exponential’s difficulty is the conditioning of a similarity, theirs is also the existence of the answer. A method has to decline rather than return a plausible matrix, which is the same discipline a proof that declines has everywhere on this site — and the interval field’s whole subject.

And the easy case is worth naming last, because it is most of what anybody computes. If A is symmetric, V is orthogonal, κ(V) = 1, and the definition is the best method there is. Diagonalise, apply f to the eigenvalues, undiagonalise; the answer is accurate to the conditioning of f on the spectrum and nothing else. Every objection in this essay is an objection about non-normal matrices, and the split is that sharp.

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. 15 The quantity that decides which side of the split a matrix is on. Where the eigenvectors are orthogonal there is nothing to lose in the similarity, and the definition and the method coincide.

One consequence of that split is worth stating on its own, because it decides what a practitioner should do rather than what they should believe. Test whether the matrix is normal before choosing a method. ‖AᵀA − AAᵀ‖/‖A‖² costs two matrix products and answers it; if the answer is at rounding, the eigendecomposition is both the best method and the most informative one, and every objection above is void. If it is not, the number that matters is κ(V), which cannot be computed without the eigendecomposition it is a warning about — so the honest move is not to compute it at all and use a method that does not need it.

What is worth carrying

A definition and a method are different objects, and the definition of f(A) is written in a quantity that is not continuous in A. V has no limit where f(A) has one, so a formula in V is not a computation of f.

Handing the route exact eigenvalues and exact eigenvectors does not save it. The experiment above does exactly that — a triangular matrix’s spectrum is its diagonal and its eigenvectors are a closed form — and the answer is still wrong by sixty-five orders of magnitude at δ = 10⁻¹². The defect is the inversion of V, and nothing upstream of it can help.

A method that is a polynomial in A inherits the continuity that f(A) has. Scaling and squaring, the Schur–Parlett recurrence, and Krylov all have that property; the eigendecomposition and the Jordan form do not.

The next essay is about the method that does work, and about the fact that it has a parameter with a worst value on both sides: the series that has to be squared back.

Where this goes next

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. 027548110813510⁻⁶10⁻⁴10⁻²110²10⁴10⁶power‖Aᵏ‖Kreiss constant 6760e · n · K‖Aᵏ‖ρᵏtwo routes to one peakspectral radius0.8peak of ‖Aᵏ‖2·10⁴Kreiss constant6757e · n · K1.1·10⁵everything here decays in the endand one of these curves says how much first A spectral radius that grows first ρ(A) below one guarantees that the powers of A go to zero and says nothing about what they do on the way. Here they rise by a factor of twenty thousand before turning over, and the peak is bracketed above and below by a constant computed from the resolvent norms outside the unit circle — two routes to one number, one through the plane and one through the powers. 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.

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.

Condition numberDefective matrixEigenvaluesEigenvectorsExact ground truthJordan formMatrix exponentialMatrix functionNon-normality