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

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.

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 — the standing this site’s known answers have. 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.

What each route returns #

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.

The sixty-five orders are one number #

The table’s right-hand column is explained — it is the distance between two matrices — and the middle column is left as “wrong by tens of orders of magnitude”. It is a quantity, and it is one that was computable before anything ran.

Put κ(V) beside it. At δ = 10⁻¹ it is 3.63·10⁵, and at each further decade it rises by a factor of 10⁷: 3.29·10¹², 3.29·10¹⁹, 3.29·10²⁶, and on to 3.29·10⁸² at δ = 10⁻¹². Divide the eigen route’s error by κ(V)·u and the answer is 0.12, 0.020, 0.045, 0.15 and 0.080 at δ from 10⁻³ down to 10⁻¹².

The error is κ(V)·u, to within a factor of eight, across sixty-three orders of magnitude. The sixty-five orders in the table are not a mystery and not a defect in the implementation. They are the condition number of a quantity nobody asked for, multiplied by the unit roundoff, exactly as a condition number is supposed to behave.

That the exponent is 10⁷ per decade is the second half of it. κ(V) grows like δ^−(n−1), because the eigenvector for λ + kδ has entries μ^(k−i)/((k−i)!δ^(k−i)) and the largest of them carries δ^−(n−1). Seven, at n = 8, to three digits in the table.

The formula has μ in it to the same power, and the slider is where that gets checked. If the largest eigenvector entry carries (μ/δ)^(n−1), then multiplying μ by sixteen should multiply κ(V) by 16⁷ = 2.68·10⁸ and leave everything else alone.

Across μ = 0.25, 0.5, 1, 2 and 4 the condition number of the eigenvector matrix reads 2·10⁷⁸, 2.6·10⁸⁰, 3.3·10⁸², 4.2·10⁸⁴ and 5.4·10⁸⁶ — a factor of 2.7·10⁸ end to end against the 2.68·10⁸ the exponent predicts, and a factor of 128 ≈ 2⁷ per doubling. The formula for the eigenvectors is being read off the drawing rather than trusted.

And the other column does not move. Scaling and squaring returns 4.2·10⁻¹², 4.2·10⁻¹², 4.1·10⁻¹², 4.0·10⁻¹² and 3.8·10⁻¹² across the same five — a spread of ten per cent, against a factor of 2.7·10⁸ in κ(V) and a factor of 1.8·10⁷ in the other route’s error. The distance to the limit is 7·10⁻¹² throughout, because μ is not in it.

The identification with κ(V)·u survives the sweep and loosens a little: the eigen error divided by κ(V)·u reads 0.11, 0.11, 0.080, 0.039 and 0.0077 across the five values of μ. Over eight orders of magnitude of κ(V) that is a drift of fifteen — still one number describing the whole family, and no longer the factor of eight the δ sweep alone reported.

Which makes the crossover exact #

It is worth noting what this does not license, because a clean identification invites one wrong use. κ(V)·u is the eigen route’s error on this family, where V is written down in closed form and no eigensolver runs. A real code computes V, and a computed eigenvector matrix near a defective limit is not merely ill-conditioned — its columns are the output of an algorithm that has itself lost the directions it was trying to separate. So κ(V)·u is a lower bound on what a real eigen route would do, not an estimate of it, and the sixty-five orders here are the optimistic case.

That is also why the family was built this way. Handing the route exact eigenvalues and exact eigenvectors removes every excuse, and the number that comes back is what the recipe costs when everything upstream of it is perfect.

The first two rows are the control, and the section above says they agree “because κ(V) is only 10⁵ and 10¹²”. That can be made precise, because two quantities are competing and both are now known.

The eigen route carries an error of about κ(V)·u ≈ δ^−(n−1)·u. The comparison against the limit carries a genuine distance of about δ. One rises as δ falls and the other falls, so they cross where

δ⁽ⁿ⁻¹⁾·δ = u, that is δ = u^(1/n).

At n = 8 that is (1.1·10⁻¹⁶)^(1/8) = 1.0·10⁻² — exactly the row where the table’s two columns stop agreeing. Nothing was fitted; the crossover is the unit roundoff raised to one over the size of the matrix.

And the practical reading is the reverse of what the framing invites. A defective limit sounds like a rare corner that only a pathologically close pair of eigenvalues reaches. u^(1/n) grows with n: it is 3·10⁻⁴ at n = 4, 10⁻² at n = 8, and 10⁻¹ at n = 16. On a sixteen-by-sixteen matrix the eigen route is already losing digits when the eigenvalues are a tenth apart, which is not close at all.

So the failure is not about being near a defective matrix in any absolute sense. It is about being within u^(1/n) of one, and that region grows towards the whole space as the matrix grows — which is why the recipe survives in textbooks, where n is three, and fails in codes, where it is not.

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 — which is the same singularity an infinite eigenvalue sits on, reached from the other side. 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 — including the one this collection computes by squaring.

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.

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.

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ᴬᵗ, 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.

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.

The exponential is where the stakes are #

Of all the matrix functions, eᴬᵗ 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ᴬᵗ 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.

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.

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.