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

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.
21 min read 16 figures Where the algebra stops being the arithmeticIdentical algebra, different arithmeticMeasured, not assertedExact ground truthAssertions that reject

Worth reading first: A function of a matrix is not a function of its entries · Cancellation takes the answer, not a digit · The direction the error leans.

The exponential of a matrix has a series, it converges for every matrix, and the convergence is absolute and fast. Every term is a matrix product and a division by an integer. There is no factorisation, no pivot, no eigenvalue, nothing to choose. It is the most obvious algorithm in the subject.

e^A  =  I + A + A²/2! + A³/3! + …

On Moler and Van Loan’s example — a 2×2 with eigenvalues −1 and −17, entries of size 64 — this converges in 77 terms and returns an answer whose relative error is 5.2·10⁻⁹.

Nine digits from double precision, on a two-by-two, from a series every term of which was computed correctly.

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. 1 The error of scaling and squaring against the number of squarings, measured against a closed-form exponential. Both ends of the axis cost. Drag the entry above the diagonal, which sets ‖A‖.

Where the digits went

The terms of that series rise before they fall. ‖A‖ is 64, so ‖A²/2!‖ is around 2,000, ‖A³/3!‖ around 40,000, and the peak — at k = 26 — is 1.4·10⁷.

The answer they sum to has norm 2.6.

largest term / sum  =  5.4 × 10⁶

Seven digits of cancellation, and cancellation is the arithmetic field’s oldest subject: the digits that cancel were never wrong, and what survives is the rounding of the inputs. A double carries sixteen digits; seven are consumed by the cancellation before truncation is even a question; nine are left, which is precisely what the measurement returns.

‖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. 2 The term norms against k, with the answer drawn as a flat line beneath them. The vertical distance between the peak and the line is the number of digits the method discards.
Relative error of two algebraically identical expressions for (1 − cos x)/x², in binary32A log–log plot of relative error against x. The expression written as it reads loses accuracy as x falls and is eventually wrong in every digit; the rearranged form stays at rounding level.10⁻⁸10⁻⁶10⁻⁴10⁻²110⁻⁹10⁻⁶10⁻³1xrelative error of the computed value(1 − cos x)/x², as written2 sin²(x/2)/x²no digits left at allbinary32 throughoutone function, two spellings · zero below 3.5·10⁻⁴
Fig. 3 The same phenomenon in its simplest form. A subtraction of two nearly equal numbers is exact; what is lost is that the numbers were rounded before it.

The boundary, which is between a half and one

The interesting thing about the series is that it is not always bad. Halve the matrix and the picture changes completely.

matrix largest term / sum series error scaling and squaring
A/4 15 1.9·10⁻¹⁵ 1.7·10⁻¹⁴
A/2 955 1.6·10⁻¹³ 3.9·10⁻¹⁴
A 5.4·10⁶ 5.2·10⁻⁹ 8.5·10⁻¹⁴
1.5A 3.6·10¹⁰ 1.6·10⁻⁵ 3.1·10⁻¹³
2A 2.4·10¹⁴ 0.107 1.7·10⁻¹³

At a quarter of the matrix the series is the better of the two methods. At twice it the answer has no correct digits at all.

Nothing about the matrix has become harder across that table. Its eigenvalues are real, negative, well separated and simply scaled; its condition number is unchanged by a scalar multiple. What changed is a property of the method — its intermediates — and a method whose accuracy depends on a scalar multiple of its input is a method with a defect, because a scalar multiple is the one transformation every user applies without thinking.

‖Aᵏ/k!‖ against ‖e^A‖ for Moler and Van Loan's 2×2, scaled by 0.25The terms of the Taylor series rise to 81.6 at k = 4 and fall back; the sum they produce has norm 5.37, drawn as the flat line. The ratio is 15.2, so that many digits are cancelled away before truncation is a question, and the computed exponential has a relative error of 1.87·10⁻¹⁵ after 36 terms. Scaling and squaring on the same matrix returns 1.65·10⁻¹⁴.0510152025303510⁻⁶10⁻⁴10⁻²110²k‖Aᵏ/k!‖‖e^A‖ = 5.4largest term 82what the series throws awaylargest term82‖e^A‖5.4digits cancelled away15error after the sum1.9·10⁻¹⁵every term is computed correctlyand the sum has lost seven digits
Fig. 4 A quarter of the matrix, where the largest term is fifteen times the answer and the series is a perfectly good method.
‖Aᵏ/k!‖ against ‖e^A‖ for Moler and Van Loan's 2×2, scaled by 2The terms of the Taylor series rise to 2.39·10¹⁴ at k = 33 and fall back; the sum they produce has norm 1.01, drawn as the flat line. The ratio is 2.36·10¹⁴, so that many digits are cancelled away before truncation is a question, and the computed exponential has a relative error of 0.107 after 127 terms. Scaling and squaring on the same matrix returns 1.72·10⁻¹³.016324864809611210⁻⁶10⁻²10²10⁶10¹⁰10¹⁴k‖Aᵏ/k!‖‖e^A‖ = 1largest term 2.4·10¹⁴what the series throws awaylargest term2.4·10¹⁴‖e^A‖1digits cancelled away2.4·10¹⁴error after the sum0.11every term is computed correctlyand the sum has lost seven digits
Fig. 5 And twice it, where the largest term is 2.4·10¹⁴ times the answer and the result has no correct digits.

So scale it down, and put it back

The repair writes itself once that table is on the page. e^A obeys

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

for any s, exactly, because A commutes with itself. Choose s so that ‖A/2ˢ‖ is small — around a half — compute the exponential of the small matrix where nothing cancels, and square the result s times.

Squaring is the right way to get back rather than multiplying 2ˢ times, because it takes s products instead of 2ˢ. At ‖A‖ = 64 and a target of 0.5, s = 7, which is seven matrix products.

The small exponential is computed by a Padé approximant rather than a Taylor sum: a ratio of two polynomials in A, chosen so that its series agrees with the exponential’s to as high an order as the degrees allow, evaluated as one linear solve. For the same number of matrix products it is accurate over a wider range, which is exactly what a method wants when it has already decided to shrink its input.

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. 6 And the whole method’s answer against the closed form, on the family from the previous essay. Scaling and squaring returns the exponential of the matrix it was handed, at every distance from the defective limit.

The cost is on both sides

Here is the part that is not obvious and is the reason this essay has a slider.

Squaring is exact in the algebra: (e^{A/2ˢ})^(2ˢ) is e^A for every s, so more squarings can only be better — they take the approximant into a range where it is more accurate. In floating point that is false.

Each squaring is a matrix product, and a matrix product’s rounding is proportional to the norms of its inputs. Squaring a matrix whose entries carry a relative error δ gives one carrying about 2δ, and s squarings give 2ˢδ. So the error the Padé step made is multiplied by roughly 2ˢ on the way back.

s ‖A‖/2ˢ relative error
0 4.2 6.0·10⁻¹²
1 2.1 1.6·10⁻¹⁵
2 1.05 6.9·10⁻¹⁶
5 0.13 4.9·10⁻¹⁵
10 0.004 3.0·10⁻¹⁴
15 1.3·10⁻⁴ 3.2·10⁻¹²

Four orders of magnitude of penalty on the left and a factor of 4,600 on the right, with a floor in between. A parameter whose worst values are at both ends is a different kind of parameter from one that merely needs to be large enough, and it is why libraries choose s from ‖A‖ by a rule rather than taking as many squarings as they can afford.

Relative error of scaling and squaring against the number of squarings, μ = 12Each 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 5.55·10⁻¹¹; the best is 9.08·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 = 05.6·10⁻¹¹best, at s = 29.1·10⁻¹⁶error at s = 153.2·10⁻¹²‖A‖ / 2^s at the best3the squarings are exact in the algebraand each one doubles the rounding
Fig. 7 The same curve for a matrix four times the norm, where the left-hand wall has moved right and the right-hand rise has not. The rule that follows is: take the smallest s that works.
Rounding 1.42 to 4 significand bits, five waysA number line of the representable values in one octave, with the exact value marked between two of them and a row for each rounding mode showing which neighbour it returns.11.131.251.381.51.631.751.882exact 1.42nearest1.375toward zero1.375toward +∞1.5toward −∞1.375stochastic36% upstochastic, measuredposition between neighbours0.36share of draws rounded up0.369 representable values between 1 and 2four functions and one distribution
Fig. 8 The mechanism on the right-hand side, in its simplest form: an operation performed on a rounded input returns an answer carrying that rounding, and doing it repeatedly compounds.

Why Padé rather than a longer Taylor sum

Once the matrix has been scaled down, the cancellation objection has gone: at ‖A/2ˢ‖ ≈ 0.5 the series’ terms fall monotonically from the first and nothing cancels at all. So why not sum the series there, and skip the linear solve?

Because the two are not competing on accuracy, they are competing on how many matrix products buy how much range, and that decides s.

A Taylor sum of degree m is accurate to about ‖A‖^(m+1)/(m+1)!, and costs m − 1 products naively or about 2√m with a clever evaluation scheme. A diagonal Padé approximant of the same total degree costs about the same and is accurate to about ‖A‖^(2p+1) times a much smaller constant, because it matches twice as many terms of the series for the same number of coefficients. In practice that is worth two to three extra squarings avoided at a fixed accuracy — and every squaring avoided is a factor of two of rounding not amplified.

So the Padé step is not there because the series is inaccurate at small norms. It is there because it lets s be smaller, and s is the parameter whose right-hand cost the previous section measured. The two halves of the method are coupled: the better the approximant, the fewer squarings, and the squarings are where the rounding is.

The proved bound against the error it bounds, at 24 significand bitsTwo curves against the Hilbert size on a logarithmic vertical axis: the error the computation actually made, which is known because the exact answer is the integers by construction, and the bound proved for it in interval arithmetic. The bound is above the error at every size, by 9.949 to 23.12 times. Past n = 5 the method refuses to return a bound at all.345678910111210⁻⁵10⁻³10⁻¹Hilbert size nerrorrefusedrefusedrefusedrefusedrefusedthe proved boundthe actual errorbound ÷ errorslack at n = 313slack at n = 423slack at n = 59.9the error is known, not estimatedthe exact answer is the integers 1 … n
Fig. 9 The general question a truncation bound answers and the one it does not: how far a series is from its limit, against how far the computed sum is from the series. Choosing an approximant is choosing the first; choosing s is choosing the second.
Three filters at their own optima, on one problemFilter factors against the singular-value index. Truncation at K = 21 is a step from one to zero; Tikhonov at λ = 0.0215 is a smooth descent; the 20-step conjugate gradient filter rises above one and changes direction 5 times. The three reach errors within a few per cent of each other while differing by up to 0.99 in individual factors.081624324048566400.250.50.7511.25index kfilter factortruncation K = 21TikhonovCGLS, 20 stepsone floor, three weightingstruncation vs Tikhonov0.77truncation vs CGLS0.99reversals in the CGLS filter5the same error, to a few per centby three different weightings
Fig. 10 And the general shape of choosing a rational approximation over a polynomial one, which recurs across this site wherever a function of an operator is being applied.

Two orders that are the same algebra

The Padé approximant is evaluated as a solve rather than as an inverse, per the phase’s standing rule, and there is a second choice inside it that matters as much.

r(A) = D(A)⁻¹N(A) can be computed as a solve with D against the columns of N, or as N(A)D(A)⁻¹, or by a partial-fraction expansion. All three are the same rational function of A. They are not the same computation, and the difference shows up on matrices with a large norm — which are exactly the ones this method has been designed to avoid producing.

That is a small instance of the site’s oldest thread, and it is worth noticing that it recurs inside the repair for a different instance of the same thread. Choosing to scale and square removed one cancellation and introduced a choice of association; the choice of association is a second place where algebraically identical routes differ.

Relative error of three summation algorithms in binary32A log–log plot of relative error against the number of terms for naive, pairwise and compensated summation, each measured against the exactly rounded sum.10¹10²10³10⁴10⁵10⁶10⁻⁹10⁻⁷10⁻⁵10⁻³10⁻¹number of terms addedrelative error against the exact sumin orderin a treecompensatedbinary32 · terms are 1/icompensated: 3·10⁻⁸
Fig. 11 The oldest instance on the site: one sum, several orders, several answers. Every reassociation in this essay is a version of that.
Forward error of a 8×8 Hilbert solve at eight precisionsA bar for each significand width from 12 to 53 bits showing the relative error in the computed solution, with the condition number times the unit roundoff marked as a prediction.κ = 1.5·10¹⁰ · the exact answer is (1, 2, …, 8)12 bits4.116 bits1.620 bits1024 bits0.530 bits0.5636 bits0.005943 bits5.5·10⁻⁶53 bits3.2·10⁻⁸dashed: κ · unit roundoffone matrix, eight arithmeticsmeasured against a known answer
Fig. 12 And what a careful order buys, which is the same amount as it buys here — a few digits, for free, from a decision nobody sees in the formula. The alternative is to spend a format on the same digits.

What the method costs, counted

The arithmetic is worth stating, because the whole method is an exchange of one kind of work for another and the exchange is favourable by a wide margin.

A Padé approximant of order p costs about p/2 matrix products to build the two polynomials — Paterson and Stockmeyer’s scheme, which evaluates a degree-p polynomial in about 2√p products — plus one n×n solve for the division. The squarings cost s products. So the whole thing is roughly

(2√p + s + 1) × 2n³   flops

At p = 13 and s = 7 that is about fourteen matrix products, which is 14 × 2n³. For comparison, one LU factorisation is (2/3)n³ and one eigendecomposition is around 25n³. The exponential is therefore about the cost of a Schur reduction, which is the right order of magnitude for an operation of this kind and is emphatically not free.

The s in that count is the part a caller controls, and it grows as log₂‖A‖. Doubling a matrix costs one more squaring, so a badly scaled problem is a linear cost in the number of decades and not an exponential one — which is the second thing the squaring construction buys, after accuracy.

The same arithmetic, twice, with M = 144 words of fast memoryThree counters against the matrix size on logarithmic axes. The operation count is one curve, because both orderings perform exactly the same operations — 338,120 of them at n = 80. The words moved between fast and slow memory are two curves a factor of 3.24 apart. The factorisations they produce are identical to the last bit.10²10³10⁴10⁵matrix size ncountoperations, bothwords, unblockedwords, blocked (b = 6)the answer does not move‖PA − LU‖/‖A‖, unblocked2.8·10⁻¹⁶‖PA − LU‖/‖A‖, blocked2.8·10⁻¹⁶difference between them0the dashed curve is both orderings' operation countthe solid pair is what they cost
Fig. 13 Counting operations, which is what this exchange is denominated in. The relevant comparison is against a factorisation, and the exponential costs a small multiple of one.

The coefficients are a trap in themselves

The diagonal Padé coefficients have a closed form:

cⱼ  =  (2p − j)! p! / [ (2p)! j! (p − j)! ]

and the standard order is p = 13. (26)! is 4·10²⁶, which is not representable as an integer in a double and rounds when written as a floating-point literal.

Computing the coefficients from the formula as written is therefore a defect that only appears at the order everybody actually uses. The recurrence

c₀ = 1,   cⱼ = cⱼ₋₁ (p − j + 1) / ((2p − j + 1) j)

has no factorial in it, keeps every intermediate near one, and is what this site’s implementation does.

It is a small thing and it belongs here, because it is the essay’s own subject appearing one level down: a formula that is exact in the algebra, written out directly, producing a number the format cannot hold.

Where each norm works, for a vector of 8Four horizontal rows, one per format. Each carries a pale bar for the format's whole range, a bar above it for the scaled norm and a shorter bar below it for the naive one.-47-37-27-17-7313233301234log₁₀ of the vector's normfp1611 bitsbfloat168 bitstf3211 bitsbinary3224 bitspale: the format's range · blue: √(Σ(xᵢ/m)²)·m · red: √(Σxᵢ²)fp16 and tf32 have the same eleven significand bitsand their bars do not overlap
Fig. 14 The same failure in the site’s oldest example of it: a quantity that is perfectly representable, computed by a route whose intermediates are not.
Refinement thresholds by format, against a problem at κ = 10⁴A horizontal bar chart of the condition number at which iterative refinement stops working, one bar per floating-point format, with a marked line at the condition number of the problem.bfloat16 · 8 bits1/u = 256fp16 · 11 bits1/u = 2048tf32 · 11 bits1/u = 2048fp32 · 24 bits1/u = 1.7·10⁷fp64 · 53 bits1/u = 9·10¹⁵κ·u = 39.063 — past the thresholdκ·u = 4.883 — past the thresholdκ·u = 4.883 — past the thresholdκ·u = 6·10⁻⁴ — refinement recoversκ·u = 1.1·10⁻¹² — refinement recoversthe problem is at κ = 10⁴; a format works when κ·u < 1thresholds are 1/u, a property of the arithmetic2 of 5 formats clear this κ
Fig. 15 And what the format can hold. (26)! is not in it, and a coefficient formula that forms it is asking for a number that does not exist.

How a library chooses, which is a backward-error argument

The parameters this essay has been sweeping by hand are chosen by a rule in every serious implementation, and the rule is worth stating because it is not the one the pictures suggest.

The obvious rule is forward: choose s and p so that the approximant’s truncation error is below the unit roundoff on the scaled matrix. That is what the classical analysis does and it over-scales — it takes more squarings than necessary, and the previous section measured what extra squarings cost.

Higham’s 2005 algorithm asks a backward question instead: for which ‖A‖ is the Padé approximant of order p the exact exponential of a matrix within u of A? That has a computable answer, one θ_p per order, and the rule becomes: pick the smallest s for which ‖A‖/2ˢ ≤ θ_p. For p = 13 the threshold is θ₁₃ ≈ 5.37, which is an order of magnitude larger than the “around a half” the naive reading suggests, and every factor of two saved there is a squaring not performed.

That is this site’s own spine deciding a library’s default: the question is not how close the approximant is to the exponential, it is which matrix’s exponential the approximant is. The two give different answers because the second is allowed to charge the error to the input.

The second refinement is the same argument applied to a norm. ‖A‖ is a crude bound on what the approximant sees — for a matrix whose powers decay it is far too pessimistic, and the scaling is then chosen for a matrix that is not the one in hand. Replacing ‖A‖ with estimates of ‖A^k‖^(1/k) for a few k gives a smaller number for exactly the non-normal matrices where a norm over-states the growth, and saves one or two squarings on them. The estimates are computed by the one-norm estimator this site already has an essay about — the one whose errors are one-sided in the flattering direction — which is safe here for a reason worth noticing: an under-estimate of ‖A^k‖ leads to under-scaling, which the θ_p test then catches, so the one-sidedness fails towards the conservative branch.

How close Hager's estimate is to the true κ₁, over 200 seeded 12×12 matricesFive bars. The estimate is exactly the true condition number on 84% of the sample and inside ten per cent on 90%; the worst underestimate in the whole sample returns 36% 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 right83.5%inside 10%89.5%inside a factor of 289.5%worst in the sample, ×10035.7%the constructed matrix, ×1007.7%usually exactexact share0.83worst of the sample0.36the constructed matrix0.077a routine that is right most of the timeand never wrong in the safe direction
Fig. 16 The estimator in question, whose one-sidedness is a liability where it is used as a certificate and an asset here, because the branch it errs towards is the one that gets checked.

There is a third refinement in the same family and it is the one that reads oddly: the algorithm also declines to be too clever. Higham’s later analysis found that the earlier version’s aggressive scaling could lose accuracy on matrices with widely spread eigenvalues, because the squaring phase amplifies the rounding of a matrix whose entries span many orders of magnitude, and the fix was a test that reverts to a more conservative s when the cheaper estimate cannot be trusted. A parameter with a worst value at both ends does not have a rule that is monotone in anything, which is what makes it worth a picture rather than a formula.

What is worth carrying

A convergent series is not a method. The exponential series converges for every matrix and converges quickly, and on an ordinary 2×2 it discards seven digits before truncation is a question. Convergence is a statement about the limit and says nothing about the intermediates.

The failure depends on a scalar multiple of the input, which is the sharpest available evidence that it belongs to the method. Nothing about a matrix becomes harder when it is doubled, and the series’ answer goes from thirteen correct digits to none.

And a closed form for the coefficients is a trap at the order everybody uses. (26)! is 4·10²⁶ and is not an integer any double holds; the recurrence keeps every intermediate near one and computes the same numbers.

And the repair has a cost on both sides. Too few squarings leaves the approximant outside its range; too many multiplies its rounding by 2ˢ. The parameter has a floor, not a direction, and the rule that follows is to take the smallest s that works rather than the largest that is affordable.

The next essay stops computing the matrix altogether. What a differential equation solver needs is e^{At}b, one vector, and there is a method that gets it in twenty matrix–vector products without ever forming an n×n exponential: the vector was what was wanted.

Where this goes next

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 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. 10⁻¹²10⁻¹⁰10⁻⁸10⁻⁶10⁻⁴10⁻²110⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹xrelative error of the computed value(1 − cos x)/x², as written2 sin²(x/2)/x²no digits left at allbinary64 throughoutone function, two spellings · zero below 1.5·10⁻⁸ Cancellation takes the answer, not a digit Subtracting two nearly equal numbers is exact. That is what makes it dangerous — the subtraction introduces no error at all, it exposes error the operands were already carrying, and the exposure can consume every significant figure at once. 10¹10²10³10⁴10⁵10⁶10⁻⁹10⁻⁷10⁻⁵10⁻³10⁻¹number of terms addedrelative error against the exact sumin orderin a treecompensatedbinary32 · terms are 1/icompensated: 3·10⁻⁸ The order they are added in Addition is associative in the algebra and is not associative in the arithmetic. The same million numbers, added in a different order, give answers that differ in the third significant figure — and the fix is not a wider float, it is a different order.

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.

CancellationExact ground truthFlop countMatrix exponentialMatrix functionNon-normalityPadé approximantRounding