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

The vector was what was wanted

Nobody who computes a matrix exponential wants the matrix. They want e^{At}b — one vector, the state of a system at a later time. Twenty matrix–vector products get it to sixteen digits on a hundred-by-hundred problem, without ever forming a hundred-by-hundred exponential, and the exponential that does get computed is twenty by twenty.
22 min read 18 figures Where the algebra stops being the arithmeticMeasured, not assertedTwo routes to a numberExact ground truthAssertions that reject

Worth reading first: A function of a matrix is not a function of its entries · The rate the condition number predicts · The spectrum that predicts nothing.

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^{At}x(0). A continuous-time Markov chain’s distribution at time t is e^{Qt}p₀. 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Δt}x.

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.

Relative error of e^Ab by Krylov against the number of matrix–vector products, n = 100The exponential of a 100×100 matrix is never formed: m products with A build a Krylov basis, the exponential of the m×m Hessenberg matrix is taken, and one combination of the basis vectors is the answer. The error falls superlinearly — 0.026, 8.4·10⁻⁴, 1.5·10⁻⁵, 1.6·10⁻⁷, 1.1·10⁻⁹ at m = 4, 6, 8, 10, 12 — and crosses the accuracy of the full dense exponential, drawn as the flat line at 3·10⁻¹⁶, at m = 18. That is 0.36 against 2.00 megaflops, on a dense matrix; on a sparse one the ratio is far larger.0481216202410⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹matrix–vector products, mrelative error in e^Abforming e^A: 3·10⁻¹⁶crosses at m = 18the vector, not the matrixsteps to the dense answer18dimension100Krylov megaflops0.36dense megaflops2the exponential that is computedis 18×18
Fig. 1 The relative error of e^Ab computed by Krylov, against the number of matrix–vector products, with the accuracy of the full dense exponential drawn as the flat line. Drag the size of the matrix: the step count barely moves and the cost of the dense route grows as n³.

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_m = V_m H_m + h_{m+1,m} v_{m+1}e_mᵀ. Then

e^{tA}b  ≈  ‖b‖ · V_m · e^{tH_m} · 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.

The residual basis of conjugate gradients, at κ = 106A semi-logarithmic plot against iteration count showing the loss of orthogonality among the residual vectors rising while the relative residual falls.05101520253010⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹iteration‖RᵀR − I‖ and ‖r‖/‖b‖step n‖RᵀR − I‖residual30×30, run for exactly n stepsexact arithmetic would end here
Fig. 2 The basis the method builds, and the property it depends on. A Krylov basis that has stopped being orthonormal is not a basis, which is why the Arnoldi step here reorthogonalises once.

What it costs and what it returns

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

m relative error
5 5.1·10⁻³
10 1.6·10⁻⁷
15 4.0·10⁻¹³
20 4.3·10⁻¹⁶

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.

Relative error of e^Ab by Krylov against the number of matrix–vector products, n = 160The exponential of a 160×160 matrix is never formed: m products with A build a Krylov basis, the exponential of the m×m Hessenberg matrix is taken, and one combination of the basis vectors is the answer. The error falls superlinearly — 0.027, 8.7·10⁻⁴, 1.5·10⁻⁵, 1.7·10⁻⁷, 1.2·10⁻⁹ at m = 4, 6, 8, 10, 12 — and crosses the accuracy of the full dense exponential, drawn as the flat line at 2.9·10⁻¹⁶, at m = 18. That is 0.92 against 8.19 megaflops, on a dense matrix; on a sparse one the ratio is far larger.0481216202410⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹matrix–vector products, mrelative error in e^Abforming e^A: 2.9·10⁻¹⁶crosses at m = 18the vector, not the matrixsteps to the dense answer18dimension160Krylov megaflops0.92dense megaflops8.2the exponential that is computedis 18×18
Fig. 3 The same measurement at n = 160. The step count is unchanged and the dense route’s cost has grown by a factor of four, which is the whole economics in two numbers.

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.

Two Krylov methods against products with A, at a kernel shift of 1Two error curves against the number of products with A, on a logarithmic vertical axis. The Arnoldi method reaches 0.1532 after 4 products and is 9.18 by the end of the run. The bidiagonal method reaches 0.1367 after 42 and degrades far more slowly.16111621263136414651566110⁻¹110¹products with Arelative errorArnoldi's best: 4Arnoldibidiagonalwhat a step buysArnoldi's best0.15products to reach it4bidiagonal's best0.14products to reach it42a tenth of the work to the same answerand no time at all spent there
Fig. 4 The general shape of a Krylov method’s convergence, from the iterative field. The turn is what makes the step count a property of the operator rather than of the discretisation.
Conjugate gradients at κ = 103, against the bound κ permitsA semi-logarithmic plot of the relative A-norm error against iteration count. The measured curve falls below a smooth dashed curve showing the classical condition-number bound.0408012016020024010⁻¹⁶10⁻¹³10⁻¹⁰10⁻⁷10⁻⁴10⁻¹iteration‖e‖_A / ‖e₀‖_Ameasuredκ bound85 steps40×40, spectrum spread evenly in logbound permits 448
Fig. 5 And its best-known instance. The same argument — a Krylov method is a polynomial in A, and the polynomial’s degree is what is being chosen — governs both.

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^{At}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_{m+1,m} · |eₘᵀ e^{tH_m} 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.

How close Hager's estimate is to the true κ₁, over 200 seeded 8×8 matricesFive bars. The estimate is exactly the true condition number on 81% of the sample and inside ten per cent on 87%; the worst underestimate in the whole sample returns 38% of the truth. The last bar is the matrix built to defeat it, at 7.7%, well below anything the random sample reached.each bar is a percentage — of the sample, or of the true condition numberexactly right80.5%inside 10%87.0%inside a factor of 287.0%worst in the sample, ×10037.7%the constructed matrix, ×1007.7%usually exactexact share0.81worst of the sample0.38the constructed matrix0.077a routine that is right most of the timeand never wrong in the safe direction
Fig. 6 The general shape of the problem: an estimate whose errors are one-sided, on a quantity nobody can afford to compute. Every cheap certificate on this site has that property.
Adding 0.1 to 256, a thousand times, at 8 significand bitsThree lines against the number of additions: the exact total rising steadily, a staircase following it, and a flat line at the starting value.02505007501000250275300325350375additionsrunning totalround to nearest: nothing arrivesexactstochasticnearesta thousand additionshalf an ulp at 2561moves, round to nearest0moves, stochastic46relative error, nearest0.28relative error, stochastic0.0228 significand bits, unbounded exponenta flat line is not a small error
Fig. 7 And the failure a stopping test has to survive — a method that stops improving without saying so. A generalised residual reports on the subspace rather than on the answer, which is the gap the extra steps close.

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.

Fill growth under natural: the factor rises as n^1.49A log-log plot of nonzero count against matrix dimension. The matrix's own count is a straight line of slope one; the factor's is steeper; a dense factor is steeper still.10²10².³10²10³10⁴n (dimension)nonzerosdense factorfactor, n^1.49matrix, n^1.05grid Laplacians from 5×5 to 12×12fitted, not quoted
Fig. 8 The same phenomenon in the field that named it. A factor of a sparse matrix is not sparse; an exponential of one is worse, since it fills to the transitive closure immediately.
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. 9 And what the operator’s spectrum does to the step count, which is the quantity that replaced the dimension. A clustered spectrum is a low polynomial degree, and a spread one is a higher degree.

What is given up

Two things, and both are worth stating plainly.

The answer is for one b. Forming e^{At} 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 usually still favours Krylov, because m products per vector is 2mn² and the dense route’s 2n³ is only amortised after n/m of them, but the arithmetic has to be done rather than assumed.

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.

12 restarts keeping 4 of 8, on a 40×40 matrixThe residual bound of the worst wanted eigenvalue and its true error, against the number of products with A. The bound falls from 1.21 to 1.08·10⁻¹³ across 12 cycles and 140 products, and the true error reaches 7.11·10⁻¹⁵. The basis is 8 vectors at every cycle and never grows.0183654729010812610⁻¹⁶10⁻¹³10⁻¹⁰10⁻⁷10⁻⁴10⁻¹products with Asizeresidual boundtrue errorbounded memorybasis vectors kept8products with A140worst error in the k wanted7.1·10⁻¹⁵the bound is free and the error is notand the basis never grows
Fig. 10 The standard remedy when a Krylov space has to be kept small — restarting — which applies here in the same form and for the same reason: a basis that grows without bound is the method’s real cost.
What each reorthogonalisation costs and what it leaves, n = 40Two quantities for three variants on one logarithmic axis: the number of projections performed, and the departure from orthogonality left behind. Doing nothing costs nothing and leaves a basis that is not one. Full reorthogonalisation costs 820 projections and leaves 1.9·10⁻¹⁵. Selective costs 137 — 6.0 times fewer — and leaves 10·10⁻⁸, which is √u and is enough.010⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹count, and ‖QᵀQ − I‖√u — semi-orthogonalitynoneselectivefullprojections‖QᵀQ − I‖two counters, three methodsghosts, none4ghosts, selective0selective is cheaper by6the cheap repair leaves a basis that is not orthogonaland that is the claim, not a shortfall
Fig. 11 And what keeping the basis orthogonal costs, which is the other half of it. Both are the price of a subspace method and neither depends on n.

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.

Error of a difference quotient for J(x)v against ε, on the Bratu problem at n = 64The Jacobian of this problem is a formula, so the error of each quotient is measured against a derivative that is exact rather than against a better quotient. The forward difference falls with a slope of 1.01 — first order — reaches 1.28·10⁻¹⁰ at ε = 10^-6, and rises again with a slope of -1.00 as cancellation takes over. The central difference falls with a slope of 2.00 and bottoms at 1.11·10⁻¹² for two residual evaluations instead of one. Neither gets near the unit roundoff at any ε.10⁻¹⁷10⁻¹⁵10⁻¹³10⁻¹¹10⁻⁹10⁻⁷10⁻⁵10⁻³10⁻¹10⁻¹⁴10⁻¹²10⁻¹⁰10⁻⁸10⁻⁶10⁻⁴10⁻²1εrelative error in J(x)vforwardcentralcancellationtruncationagainst a derivative that is exactforward floor1.3·10⁻¹⁰central floor1.1·10⁻¹²truncation slope, forward1truncation slope, central2no ε reaches the roundoffand the analytic derivative is free of the choice
Fig. 12 The commonest such operator: a Jacobian that is a difference quotient, accurate to ten digits rather than sixteen, and available where no matrix is.
Newton on the Bratu problem: analytic Jacobian, differenced Jacobian at ε = 10^-6, and matrix-freeThree residual sequences, all starting from zero. The analytic Jacobian gives 8, 0.052, 3·10⁻⁶, 5.2·10⁻¹³; a Jacobian differenced at ε = 10^-6, whose entries are correct to about 1.3·10⁻¹⁰, gives 8, 0.052, 3·10⁻⁶, 4.6·10⁻¹³. They stop at the same residual. The matrix-free run, where GMRES sees only a closure and no entry exists anywhere, lands 8.1·10⁻¹⁶ from the analytic answer in 438 products.0123456710⁻¹⁴10⁻¹²10⁻¹⁰10⁻⁸10⁻⁶10⁻⁴10⁻²1Newton step‖F(x)‖analytic Jdifferenced Jmatrix-freeone fixed point, three derivativesanalytic floor4.2·10⁻¹³differenced floor3.9·10⁻¹³matrix-free floor3.8·10⁻¹³products used, matrix-free438the derivative chooses the stepand the residual decides the answer
Fig. 13 And what that costs a nonlinear solve, which is the subject of the essay this one hands on to.

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.

GMRES on the Laplacian and on the cyclic shift, both 12×12A semi-logarithmic plot of relative residual against iteration. One curve falls steadily; the other is flat at one for every step until the last, where it drops to zero.02468101210⁻¹⁴10⁻¹²10⁻¹⁰10⁻⁸10⁻⁶10⁻⁴10⁻²1iteration‖r‖ / ‖b‖Laplaciancyclic shiftno progress at allevery eigenvalue of the shift is on the unit circleand it predicts nothing
Fig. 14 The special case the site already had. A Krylov solve is a Krylov matrix function with f(z) = 1/z, and the spectrum decides the polynomial degree in exactly the same way.
Two filters on one sum, λ = 0.01The weight each term of the solution is given, against its index. Truncation is a step: one for the first 26 terms and zero after. Tikhonov is σ²/(σ² + λ²), which falls smoothly through the same place. The unregularised solution is the constant one, which is why it divides noise by a σ of 1.7·10⁻¹³.081624324048566400.250.50.751index kfilter factor fₖno regularisation: fₖ = 1truncationTikhonovthe same sum, three weightsTikhonov, relative error0.11truncation, relative error0.11no filter at all5.5·10⁸both filters are one expression with a different weightfₖ = 1 is the catastrophe
Fig. 15 And the same object viewed as a filter on the spectrum, which is what any polynomial in A is. The choice of f changes which part of the spectrum is emphasised and nothing about the machinery.

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^λμᵏ/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.

‖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. 16 The closed form in question, used to draw the transient. The same formula supplies the yardstick for every measurement in these three essays.
The exact solution of a 13×13 Hilbert system beside the computed oneTwo columns of numbers: the exact answer, which is the integers one to thirteen, and the answer double-precision elimination returns, with the number of correct digits beside each.H13 x = b, b formed exactly so that x = (1, 2, …, 13)12345678910111213exact1.00002.00003.00133.97865.18864.993010.46720.048921.2657-2.575319.21478.906013.5113computed6.6 correct digits4.8 correct digits3.4 correct digits2.3 correct digits1.4 correct digit0.8 correct digitno correct digitsno correct digitsno correct digitsno correct digitsno correct digits0.6 correct digit1.4 correct digitbackward error2.2·10⁻¹⁷κ = 1.7·10¹⁸The algorithm solved a neighbouring problem perfectly. That problem's answer is this one.right-hand side built in BigInt rationalsthe truth is known
Fig. 17 And the site’s other exact answer, from the field where the habit started. A method measured against a better method is measured against nothing.

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₁+t₂)A}b = e^{t₂A}(e^{t₁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.

What 3 shifts do to each direction of the starting vectorAmplification of each eigen-component against its eigenvalue, on a logarithmic vertical axis. The points are measured — the component of the vector along each eigenvector, before and after — and the curve is |Π(λ − θ)| with the roots at the 3 discarded Ritz values. They agree to 3.9·10⁻¹¹. The 3 wanted directions are amplified by at least 1 and everything else by at most 0.802.1357910⁻⁵10⁻⁴10⁻³10⁻²10⁻¹1eigenvalue λamplificationkeptdiscardeda filter on the starting vectormeasured against the polynomial3.9·10⁻¹¹worst kept direction1best discarded direction0.8the roots are the discarded Ritz valuesand the vertical lines are where they sit
Fig. 18 The restart, from the field that named it. For a linear solve a restart discards information; for an exponential it is a time step, because the operator composes.

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^{At} appears in the notation because that is how the solution operator is written; what a computation consumes is e^{At}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.

Where this goes next

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 The series that has to be squared back The Taylor series for the matrix exponential is not wrong — every term is computed correctly — and on Moler and Van Loan's two-by-two its largest term is 5.4 million times the answer it sums to. The method that replaces it scales the matrix down and squares the result back, and both halves of that sentence cost: too few squarings and the approximant is out of range, too many and each one doubles the rounding. 05101520253010⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹iteration‖RᵀR − I‖ and ‖r‖/‖b‖step n‖RᵀR − I‖residual30×30, run for exactly n stepsexact arithmetic would end here An orthogonalisation nobody calls one Conjugate gradients are derived as a minimisation and behave as an orthogonalisation, which is why the finite-termination property in every textbook is not a property the method has in floating point. 10⁻¹⁷10⁻¹⁵10⁻¹³10⁻¹¹10⁻⁹10⁻⁷10⁻⁵10⁻³10⁻¹10⁻¹⁴10⁻¹²10⁻¹⁰10⁻⁸10⁻⁶10⁻⁴10⁻²1εrelative error in J(x)vforwardcentralcancellationtruncationagainst a derivative that is exactforward floor1.3·10⁻¹⁰central floor1.1·10⁻¹²truncation slope, forward1truncation slope, central2no ε reaches the roundoffand the analytic derivative is free of the choice An operator with no entries At the sizes where linear algebra is expensive the matrix does not exist. What exists is a subroutine that returns Av. Every Krylov method survives that unchanged; every algorithm that reads an entry disappears. And the derivative such a code computes is accurate to ten digits instead of sixteen, which turns out to cost nothing at all.

What links here

Computed from the collection, not written here: the essays that point at this one.

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.

ArnoldiExact ground truthFlop countKrylov subspaceMatrix exponentialMatrix-freeMatrix functionSparsity