Numerical precision of coefficients

[MikeInnes]

In a simple case we get poor precision on the 177th term, and it zeros out on the 178th.

julia> x = 5 + Taylor1(200)
 5.0 + 1.0 t + 𝒪(t²⁰¹)

julia> derivative(sin(x), 177)
 0.3461310333647123 + 𝒪(t²⁰¹)

julia> derivative(sin(x), 178)
 0.0 + 𝒪(t²⁰¹)

Getting accurate and extremely fast 176th derivatives is hardly unimpressive, but I figured I'd open this both so that anyone else who hits this isn't confused, and to speculate about whether TaylorSeries can be made to work to arbitrary order.

Presumably the core issue is that the factorial terms 1/k! get too small for Float64 to handle. Would it be possible/performant to store the gradients without these factors, and dynamically include them when multiplying series together?

[lbenet]

Thanks for reporting!

Indeed, the normalized coefficients of the Taylor series are simply too small to be handled by Float64; note that we compute the coefficients using recurrence relations, which include the 1/k!s without computing them explicitly/separately. Using your example:

julia> x = 5 + Taylor1(200)
 5.0 + 1.0 t + 𝒪(t²⁰¹)

julia> sin_x = sin(x);

julia> sin_x[177]  # same as getcoeff(sin_x, 177)
1.0e-323

julia> sin_x[178]
0.0

Note that the 177th (normalized) coefficient of the Taylor series is already a subnormal Float64, so I guess there may be precision problems already before.

Is it of any help to use BigFloats?

julia> xb = 5 + Taylor1(BigFloat, 200)
 5.0 + 1.0 t + 𝒪(t²⁰¹)

julia> sin_xb = sin(xb);

julia> sin_xb[177]
8.097958712298242025294651104096755715810462343573419929878033256108316811354375e-324

julia> sin_xb[178]
1.537936570049971023666298556096633985711543026211028342942105542375668955243689e-325

Let me think it over if there is a simple (recursive) way to compute the derivatives directly.

[dpsanders]

I think what Mike is asking is if we could have a version where the coefficients stored are n! times the current coefficients. (I had already asked this same question.)

This should be easy to implement (just don't divide by k when calculating the coefficients), but it's not so clear what API to use. Basically this would be a different type.

[dpsanders]

Something like the exponential generating function of the derivatives (?)

[MikeInnes]

Right, exactly. It could be the same type if you modified the display and getindex functions to apply the 1/n! term, so that there's no user-visible change to how TaylorSeries behaves. The potential downside is that some taylor series' coefficients might now overflow the other way, but I don't know if that's more or less common than underflow.

[dpsanders]

Hmm, I think this would be much harder that it seems at a first glance, since the product of two series would no longer be the simple Cauchy product.

[lbenet]

In terms of recurrence relations, things are not too far from what we already have. As for products of series, the Cauchy product will not be so simple, indeed, but it is doable.

As Mike points out, overflow may be a concern, which is obviously controlled by the 1/k! factors. With regards to this point, simply note that for Float64, factorial(21) already throws an OverflowError, so I fear that it may/will constrain things much more.