The vector was what was wanted

Ask why anybody wants a matrix exponential and the answer is always the same. A linear system of ordinary differential equations, ẋ = Ax, has the solution x(t) = eᴬᵗx(0). A continuous-time Markov chain’s distribution at time t is $e^{Qt}p_0$. An exponential integrator advances a stiff problem by applying $e^{hA}$ to a residual. A Kalman filter propagates a state across an interval by $e^{F\Delta t}x$ — and implements it as a solve for the same reason.

In every one of those the object wanted is a vector. The n×n exponential appears only because it is how the vector is written down.

The method in four lines #

Build an orthonormal basis of the Krylov space by the Arnoldi process:

V₁ = b/‖b‖;   for k = 1…m:  w = A Vₖ,  orthogonalise against V₁…Vₖ,  normalise

which produces V with orthonormal columns and an (m+1)×m Hessenberg matrix H with AVₘ = Vₘ Hₘ + hₘ₊₁,ₘ vₘ₊₁eₘᵀ. Then

eᵗᴬb  ≈  ‖b‖ · Vₘ · e^(tHₘ) · e₁

The exponential that is actually computed is m×m, with m in the tens, and it is computed by the scaling-and-squaring of the previous essay. The n×n exponential is never formed.

What A had to supply is m products with a vector. Nothing else.

What it costs and what it returns #

On a 100×100 matrix whose exponential is known in closed form:

and forming the dense exponential and multiplying returns 3.0·10⁻¹⁶. The Krylov route reaches the dense route’s accuracy at m = 20, and it does so at

2mn²  =  0.40 megaflops     against     2n³ ≈ 2.00 megaflops

A fifth of the arithmetic, on a dense matrix. That is the least favourable case for the comparison, because a dense matrix–vector product is the most expensive kind there is.

The convergence is superlinear, which is why m stays small #

The error falls by four orders of magnitude for every five steps at the start, then by six, then by three — it is not a fixed rate. That is the characteristic behaviour of a Krylov method applied to a function of an operator, and it is worth saying why, because it is what makes the step count insensitive to the dimension.

The Krylov approximation is exactly the best polynomial approximation, in a certain weighted sense, to e^z over the spectrum of A, of degree m − 1. Polynomial approximation of an entire function on a bounded set converges superlinearly: the error is bounded by something like (ρ/m)^m, which falls faster than any geometric rate once m passes the norm.

So the step count needed is governed by ‖A‖t and by how tightly the spectrum clusters, and not by n. Doubling the dimension of a discretisation while keeping the operator the same does not change how many steps are needed; it quadruples the cost of forming the exponential.

That is the claim, and it is the kind that is easy to state and cheap to check.

The crossing step is the same integer at every size. At n = 40, 80, 120, 160 and 200 the Krylov route overtakes the dense one at m = 18 — not approximately, not with a trend, the same step — while the dense route’s own accuracy sits at 2.9·10⁻¹⁶ to 3.0·10⁻¹⁶ throughout, so nothing about the comparison is drifting underneath it.

What grows is the size of the saving, not the work needed to get it. The ratio of the two costs runs 2.2, 4.4, 6.7 and 11.1 across those four sizes — climbing as n² against a step count that does not move at all. That is the practical content of “insensitive to the dimension”: the method does not merely survive refinement, it improves under it, and the eighteen steps that were barely worth taking at n = 40 are an eleven-fold saving at n = 200 and a hundred-fold one at the sizes this field’s problems actually reach.

The stopping test, which is the one genuine difficulty #

A Krylov method has to decide when to stop, and here the usual device is unavailable. A linear solver can compute a residual — b − Ax̂ — and it costs one product. There is no residual for a matrix function: eᴬᵗb − y is not something a computation can evaluate without knowing the answer.

What is used instead is the generalised residual, an estimate built from the last subdiagonal entry of the Hessenberg matrix:

error  ≈  ‖b‖ · hₘ₊₁,ₘ · |eₘᵀ e^(tHₘ) e₁|

which is the first term of an exact series expansion of the error, costs nothing beyond the small exponential already being computed, and is asymptotically correct — it becomes tight as the method converges. It is also an estimate rather than a bound, which is a real distinction and the site has been careful about it before: an estimate that can be fooled is about a quantity with the same character, and the same one-sidedness applies. A small generalised residual is evidence; a large one is not proof of anything.

The honest practice is to run a few steps past where the estimate says to stop, which costs a product each and removes the question. On the measurement above the estimate at m = 15 is 3·10⁻¹³ against a true error of 4·10⁻¹³, so it is tight where it matters and it earned no special handling.

Where it becomes not a saving but the only option #

The comparison above is on a dense matrix, where the dense route at least exists. It usually does not.

A discretised diffusion operator on a 100×100×100 grid is a million by a million. Its exponential is a 10¹² matrix — eight terabytes, dense, since the exponential of a sparse matrix is not sparse. A matrix–vector product with it is seven million operations, because the operator is a stencil.

At that size the two columns of the comparison are not two costs. One of them is a cost and the other is an impossibility, and the Krylov route’s twenty or fifty products are the same twenty or fifty they were at n = 100.

The exponential of a sparse matrix being dense is worth pausing on, because it is the same fact the sparsity field is built around: eliminating a variable couples everything it touched, and exponentiating couples everything reachable. e^A’s (i, j) entry is nonzero whenever there is a path from i to j of any length, which on a connected graph is everywhere.

The degree grows like the square root of the time #

The caveat below — that the accuracy is not uniform in t, because the polynomial degree needed grows with ‖A‖t — is true and leaves out the exponent, and the exponent is the whole of what a code does about it.

The smallest m reaching a relative error of 10⁻¹², on the same 100×100 family:

at ‖A‖t = 10.0 it is 14. At 19.9, 19. At 39.8, 25. At 79.7, 36. At 159.4, 53. At 318.8, 81.

A thirty-two-fold increase in ‖A‖t costs 5.8 times the degree, against 5.7 for the square root. So m ≈ 4.4·√(‖A‖t), and the growth is sublinear by a whole power.

That reverses the obvious response to a long interval. The instinct a time integrator trains is to split: take k steps of t/k, because each is easier. Here each is easier by less than the factor of k being paid for it — k·√(‖A‖t/k) = √k·√(‖A‖t) — so k substeps cost √k times more products than one whole step. Four substeps cost twice what one step costs, for the same answer.

There are still reasons to substep, and it is worth being exact about which survive. Intermediate values, if they are wanted, have to be computed and one long step does not produce them. And the small exponential is m × m, so a large m eventually meets the squaring count’s own best value and the conditioning of the Hessenberg matrix. Neither of those is accuracy, which is what substepping is usually done for, and for accuracy it is a straight loss.

One more consequence of the square root, because it decides how the method scales into the regime it exists for. A discretised diffusion operator’s norm grows like h⁻² as the grid is refined, so ‖A‖t at a fixed physical time grows like n^(2/d) in n unknowns. Under a linear degree law that would make the Krylov degree grow with the mesh and destroy the claim two sections above — that the step count is insensitive to n. Under a square-root law it grows like n^(1/d), which on a three-dimensional grid is the cube root of the unknown count: eight times as many unknowns costs twice the degree.

So the two claims are consistent, and the reason they are is the exponent. “The step count does not depend on n” is true of the operator held fixed and false of a refinement sequence, and the square root is what keeps the falsehood mild — the same shape as a multigrid rate that is constant in n while its cost per cycle is not.

And the number of right-hand sides favours the dense route #

The paragraph above defers an arithmetic, and doing it changes the conclusion.

Forming eᴬᵗ costs 2n³ and then 2n² a vector. The Krylov route costs 2mn² a vector, every vector, with nothing amortised — each right-hand side generates its own Krylov space. They cross at M = n/(m − 1), which at n = 100 and m = 20 is 5.3.

So at one right-hand side the Krylov route costs 0.40 megaflops against 2.02, and at six the dense route has taken the lead — 2.14 against 2.40 — and it never gives it back. On a dense matrix, “many right-hand sides” is six, and the many-right-hand-sides regime is precisely where the argument for Krylov stops working rather than where it is strongest.

What does favour Krylov is sparsity, and it is a different argument that happens to point the same way at M = 1. A product with a seven-point stencil is 14n rather than 2n², so the Krylov cost per vector falls by a factor of n while the dense route stops existing at all: at n = 10⁶ the crossover is not at 50,000 right-hand sides, it is nowhere, because there is no second route to cross.

Keeping the two apart matters because they have opposite shapes. The sparsity argument has no crossover and gets stronger with n. The right-hand-side argument has a crossover at six and gets weaker with M. A code choosing between the routes on a dense matrix should count its right-hand sides; one on a sparse operator should not bother, and neither should treat the other’s reasoning as support.

What is given up #

Two things, and both are worth stating plainly.

The answer is for one b. Forming eᴬᵗ gives an object that can be applied to any vector; the Krylov route answers one question. Where many right-hand sides are wanted the comparison changes, and the next-but-one section does the arithmetic rather than assuming it — because the answer is not the one this paragraph used to assert.

The accuracy is not uniform in t. The polynomial degree needed grows with ‖A‖t, so a long time interval needs more steps, and past some t the sensible thing is to break the interval and apply the method twice. That is not a defect — it is the same time-stepping decision every integrator makes — but it means the method has a parameter where the dense route does not.

What is not given up is accuracy at a given step count, which is the thing one would expect a shortcut to cost. At m = 20 the two routes agree to the last bit.

The operator does not have to be a matrix #

The Arnoldi process asks A for exactly one thing: the product Av. It never indexes an entry, never compares two of them, never asks how many are zero.

So A can be a subroutine. A finite element code that assembles the product element by element, a spectral method that transforms and multiplies, a Jacobian that exists only as a difference quotient — all of them can be exponentiated, and none of them can be factorised, scaled, pivoted, or written down.

That is the property the last essay of this phase is entirely about, and it is worth noticing here that it arrives as a consequence rather than as a design goal. The method was chosen because forming an n×n exponential is expensive. What came with it is that the matrix need not exist.

Three functions that are the same method #

The construction has nothing in it specific to the exponential. Replace $e^{tH_m}$ by $f(H_m)$ and the same four lines compute f(A)b for any f the small matrix can be handed to, which in practice means any f with a scaling-and-squaring or a Schur–Parlett implementation.

Three that are actually computed this way:

The φ-functions of exponential integrators. φ₁(z) = (e^z − 1)/z and its relatives, which appear in every exponential Runge–Kutta method. They are computed on the m×m matrix, where their own cancellation at small z is a manageable problem, and never on the n×n one.

The square root and the sign. Used in lattice quantum chromodynamics for the overlap Dirac operator, where the matrix is of order 10⁷ and the sign function of it has to be applied to a vector several thousand times.

The inverse. A⁻¹b is f(A)b with f(z) = 1/z, and the Arnoldi construction with that f is GMRES, which the site already has. That is worth noticing rather than passing over: the whole iterative field is the special case of this essay where the function is a reciprocal, and every argument here about step counts, spectra and stopping tests is the argument that field makes.

Two routes, and the second one is a theorem #

The measurement above compares Krylov against scaling and squaring, which are both float methods, so neither can confirm the other. The matrix used is λI + μS with S the shift, whose exponential is a finite sum with entries in closed form — $e^{\lambda}\mu^k/k!$ on the kth superdiagonal — so both routes are measured against a theorem rather than against each other.

That is the site’s standing habit and it does real work here. Both float routes could have been wrong in the same direction, as they were in the determinant sweep, where two independent computations agreed with each other and not with the exact answer. Here they agree with the exact answer, and the agreement is evidence.

What the method cannot do, and the two standard repairs #

Two limitations are worth stating, because both have standard answers and both answers change the method’s character.

The basis grows. m steps store m vectors of length n, and the Arnoldi orthogonalisation costs O(m²n) — so a problem needing m = 200 stores 200 full vectors and spends most of its time orthogonalising. The repair is restarting: run k steps, form the approximation, and use the error as the new starting vector. That works for the exponential in a way it does not for an eigenvalue problem, because $e^{(t_1+t_2)A}b = e^{t_2 A}(e^{t_1 A}b)$ — the operator composes, so a restart is a genuine time-stepping rather than an approximation thrown away. What it costs is that each restart’s Krylov space starts again from a single vector.

And the convergence is slow for a spread spectrum. The polynomial degree needed grows with the spread of ‖A‖t, which for a stiff problem is exactly what is large. The repair is shift-and-invert: run the Krylov method on (I − γA)⁻¹ instead of on A, which maps the far-left eigenvalues into a small cluster near zero and converges in a handful of steps. It works extremely well, and it needs a solve per step — so it takes the method straight back out of the matrix-free world it was chosen for. That is a real trade rather than a refinement: a stiff problem can have fast convergence or no matrix, and choosing is the modelling decision.

Both repairs are the ones the iterative field already knows under other names, which is the point this phase keeps arriving at from different directions. A matrix function’s Krylov method is a Krylov method, its restart is a restart, and its preconditioner is a preconditioner; the only thing that is new is which polynomial the space is being asked for.

A third limitation has no repair and is worth being honest about: there is no a priori bound on m. The step count depends on the spectrum through a polynomial approximation problem nobody solves in advance, so a code cannot size its storage before running. Every implementation therefore grows its basis and restarts when it hits a limit somebody chose, and the limit is a memory budget rather than a numerical decision. That is the same situation GMRES has been in since 1986, and for the same reason.

What is worth carrying #

Nobody wants the matrix. eᴬᵗ appears in the notation because that is how the solution operator is written; what a computation consumes is eᴬᵗb, and it is a vector.

The step count is a property of the operator and not of the dimension. Twenty products at n = 100 and twenty at n = 160, while the dense route’s cost went up by a factor of four. Refining a discretisation makes the argument better, which is the opposite of how a shortcut usually behaves.

The stopping test is an estimate rather than a residual, because a matrix function has no residual to compute. It is asymptotically tight and it is one-sided, so a few steps past where it stops is the honest practice.

And it costs nothing in accuracy. At m = 20 the Krylov answer and the dense answer agree to the last bit, both against a closed form. What is given up is generality — one vector, one time — and not digits.

The phase’s remaining essays leave matrix functions. The next one asks a question about orthogonality that turns out to have the same shape: the nearest orthogonal matrix.