The series that has to be squared back

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⁻⁹, measured against an answer that is known.

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

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.

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.

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 — the same objection the determinant cannot answer.

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.

“At every distance” is one axis of that figure. The other is μ, the size of the matrix being exponentiated, and turning it is what separates the method from its most common alternative.

That is the pattern for the rest of the slider, and both halves of it are worth stating precisely because the second is the one a reader will disbelieve.

Scaling and squaring does not notice μ, and improves slightly across it. Its error reads 4.2·10⁻¹², 4.2·10⁻¹², 4.1·10⁻¹², 4.0·10⁻¹² and 3.8·10⁻¹² at μ = 0.25, 0.5, 1, 2 and 4 — a ten per cent move over a sixteen-fold range of the matrix. The distance to the defective limit is 7·10⁻¹² at every one of those stops, so the method is sitting at about six tenths of the distance the problem itself allows, and the remaining tenth of a per cent of variation is not about μ.

The eigenvector route notices μ by seven orders of magnitude: 2.5·10⁶¹, 3·10⁶³, 2.9·10⁶⁵, 1.8·10⁶⁷, 4.6·10⁶⁸. And the reason is printed beside it. κ(V), the conditioning of the eigenvector matrix it inverts, reads 2·10⁷⁸, 2.6·10⁸⁰, 3.3·10⁸², 4.2·10⁸⁴ and 5.4·10⁸⁶ — eight orders across the same sweep, tracking the error term for term. Dividing one by the other gives 1.3·10⁻¹⁷, 1.2·10⁻¹⁷, 8.8·10⁻¹⁸, 4.3·10⁻¹⁸ and 8.5·10⁻¹⁹: a ratio that moves by fifteen while its two terms move by a hundred million.

So the failure is not the exponential and not the series. It is the diagonalisation, priced at roughly κ(V) times the unit roundoff, on a family whose whole point is that V becomes singular in the limit. Scaling and squaring wins here by never forming V, which is the same reason it wins on every defective matrix and is a stronger statement than “it is more accurate”.

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.

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.

Where the rule puts s, and where the floor is #

The table above finds the floor by sweeping. A library cannot sweep — it has one matrix and no reference answer — so it chooses s from ‖A‖ by a rule, and the standard rule is scale until ‖A/2ˢ‖ ≤ 0.5. It is worth asking where that lands relative to the floor, because the answer is not on it.

At ‖A‖ = 4.2 the floor is at s = 2, where the scaled norm is 1.05, and the rule takes s = 4. At ‖A‖ = 16.2 the floor is at s = 2, scaled norm 4.05, and the rule takes s = 6. At ‖A‖ = 41.0 the floor is at s = 3, scaled norm 5.12, and the rule takes s = 7.

Two to four extra squarings, every time, and each one doubles the rounding the Padé step made. Priced against the floor the rule costs a factor of 5.7, 15.7 and 76.2 in relative error, and the penalty grows with the norm because a larger matrix buys more unnecessary halvings before it reaches 0.5.

The optimum is not at a scaled norm of 0.5. It is between 1 and 5 — above the rule’s target in every case measured, and increasingly so.

Why the rule is not foolish, and what it is bounding #

It would be easy to read that as an error, and it is not. The rule is chosen so that the Padé approximant’s own truncation error is provably below the unit roundoff at the scaled norm, which is a genuine bound and a good reason to want a small argument. What it bounds is one of the two costs this essay measures, and it says nothing whatever about the other.

That is the general shape and it is worth carrying past this method. A parameter with a cost at both ends cannot be set by bounding one of them. Bounding the left-hand cost and taking the smallest parameter that satisfies the bound puts the answer as far up the right-hand arm as the bound allows — which is exactly what the numbers above show, and it is the same mistake as taking as many squarings as one can afford, made carefully.

The non-normality of this family is a real part of it and not the whole. Repeating the sweep on a normal matrix with the same spectral spread gives penalties of 1.1, 10.5 and 11.3 — so over-scaling costs about an order of magnitude even when ‖A‖ is a fair description of the matrix, and non-normality multiplies that by up to a further factor of seven. The modern repair, which replaces ‖A‖ by ‖Aᵏ‖^(1/k) for a few k, is aimed at that second factor. The measurement says the second factor is the smaller half.

Which leaves the honest summary of what a library is doing: it is choosing a parameter it cannot optimise, using a bound on the only cost it can compute in advance, and paying an order of magnitude for the cost it cannot. That is a defensible engineering decision and it is not the same statement as the rule is right, and this collection’s standing habit — print the quantity that decided something beside the thing it decided — would have a library report s and the scaled norm it settled on, so that a caller who does have a reference answer can see how far up the arm it was left.

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.

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.

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 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.

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.

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