Include offset in struct definition

[KeithWM]

Hi all, thanks for the great repository!

I've been working on a package to compute the Taylor expansion of the inverse of a function y=f(x) based on the Taylor expansion of f as described by Wolfram I might want to add this here. But that's not my point here.

To me it seems that given a Taylor series expansion is an approximation to a function f(x) around some point x_0, to evaluate the approximation you ideally want to call the Taylor1 object with just the argument x, without having to think about/remember what x_0 was used in defining the expansion.

It might be an idea to implement an OffsetTaylor{N} object that houses a Taylor{N} object along with the used offsets? Curious what you guys think about this. Happy to help implementing things too!

[lbenet]

Hi @KeithWM, thanks for the opening this issue! Let me begin with a "yes", if you think the package can be improved in some way please do open issues and PRs!

To me it seems that given a Taylor series expansion is an approximation to a function f(x) around some point x_0, to evaluate the approximation you ideally want to call the Taylor1 object with just the argument x, without having to think about/remember what x_0 was used in defining the expansion.

I agree partially with what you write, so let me rewrite it in my terms. Yes, the object Taylor1 contains the Taylor series of $f(x)$ around $x_0$ (which is currently not saved at all!), but it is written in terms of the variable $h=x-x_0$ ($x$ is not specified!), rather than in terms of $x$; what I mean is contained in the first equation of the docs here. Then, yes, to evaluate the series expansion, we do not us the value of $x$ (nor $x_0$) separately, but use as $h=x-x_0$.

If I understand correctly your idea/proposal with OffsetTaylor{N} is to have $x_0$ within the struct, and then exploit it for the evaluation, which can be evaluated in terms of $x$ alone ("absolute values"), instead of $h$ ("relative values", as we do it now), but also allowing the latter. I do think this is indeed useful!

This type of idea, perhaps in a less elegantly way, is implemented/exploited in TaylorModels.jl. The subtlety there is, because evaluate uses "relative values" we must remember to subtract the expansion point x_0. Your proposal would simplify this, and make it more clear and even more general.

Finally, I'm not sure if a new type would be needed, or simply add x_0 to the structs already defined (@dpsanders proposed this in the past). My point is, maybe we can still use Taylor1 (which would include x_0) and add methods that distinguish whether the evaluation uses "absolute values" $x$ (absevaluate, or a better name) or "relative values" $h$ (current evaluate).

[KeithWM]

Hi there. Thanks for the swift and extensive answer.

I've now implemented in my own package rather compactly by using the fact that a Taylor1 object can be evaluated with another Taylor1 object as the argument. This means that to translate a Taylor expansion y = Taylor1([...], n) around x_0, which has as its "local" argument h = x - x_0, I can define a second object z = Tayloy1([-x_0, 1], n) and then use y(z)(x) to evaluate my y (defined around x_0) at an x. Of course here the object z needs to "remember" that we are about x_0, but at least it's a nice way to do the mathematic.

It's probably quite sensible to add x_0 as a field with default value of 0. But I think if you have it defined, it should be used. I'm quite convinced there's a good Julia way to create a parametric type that switches on the existence of x_0, but I'm not quite sure what it is. That is, I would suggest not using different function names, but dispatching on whether the Taylor1 (or Taylor{N}) object is defined with an Offset or not.

[lbenet]

I see your point. Aside from the discussion of having x0 inside the struct --I think we agree it may be useful--, and other questions on the implementation, let us have a working example. I'll try to use your notation above, though I think your -x0 should be x0 and y=f(z); do correct me if necessary.

julia> x = Taylor1(5) # independent variable "around 0"; at t=0 --> x = 0
 1.0 t + 𝒪(t⁶)

julia> x0 = 0.25
0.25

julia> z = x0 + x # indep variable around x0 (at t=0 --> z = x0)
 0.25 + 1.0 t + 𝒪(t⁶)

julia> y = sin(z) # TS around x0 for `sin`
 0.24740395925452294 + 0.9689124217106447 t - 0.12370197962726147 t² - 0.1614854036184408 t³ + 0.010308498302271788 t⁴ + 0.00807427018092204 t⁵ + 𝒪(t⁶)

# Your `y(z)` above corresponds to `y` here, I believe
julia> y(-x0+x) # translation of TS to be around 0
 7.207371574557975e-8 + 0.9999983177913969 t + 1.615523415198633e-5 t² - 0.16674748305763631 t³ + 0.00021566057611923858 t⁴ + 0.00807427018092204 t⁵ + 𝒪(t⁶)

# "exact" result
julia> sin(x)
 1.0 t - 0.16666666666666666 t³ + 0.008333333333333333 t⁵ + 𝒪(t⁶)

# the *same* translation using `update!` on a copy of `y`
julia> yy = y; 

julia> update!(yy, -x0);

julia> yy
 - 0.2474001108545715 + 0.9688693089044778 t + 0.12389603552765488 t² - 0.16191672477067928 t³ - 0.00987717715003331 t⁴ + 0.00807427018092204 t⁵ + 𝒪(t⁶)

Three comments:

• Exploit update!(y, -x0) accomplishes the translation by -x0 to the "new expansion point" 0, since the expansion point of yy is x0
• The result (of this translation) is an approximation, because is based in the Taylor polynomial used which is itself an approximation, and there are corrections from higher order terms. To make this point clear, repeating this but using expansions of order 15, eventually gets you to:

  julia> y - sin(x)
   2.7755575615628914e-17 - 1.1102230246251565e-16 t - ... + 𝒪(t¹⁶)

• How good the resulting translation is, aside from the order of the polynomials, also depends on how far you are translating. This is actually quite subtle, since it is related to the convergence properties of the series. So the trick works, if you stay within the convergence radius of the series (around x0).

[KeithWM]

I'm not entirely sure we're on the same page about the offsets, but maybe I'm misunderstanding the notation.

Instead of sin(x), let's say we're doing a Taylor approximation of y = log(x) around x_0 = 1 by defining an h = x - 1. Then we have the approximation Y(h) = -h + h^2/2 - h^3 - h + O(h^4). This is also more or less the function signature of the current Taylor1 implementation. But if the objcet knows about x_0, we can substitute h = x - x_0 to create a call signature Y2(x) = Y(x_0 + h) [I think I made a mistake here before?]. A user can then use this instead of having to call Y(x - x_0).

This does indeed create a polynomial expansion about x - 0 that is an approximation of log(x), but this approximation will be terribly inaccurate around x = 0 where the log is undefined. The approximation is of course just as accurate as the one centered about x_0 = 1 if you look at the whole range of values for x.

[lbenet]

I agree with you that there is an issue with the notation... Let us then expand log(x) around x0 = 1, i.e. the "small" quantity is h = x - x0 (which below appears as t):

julia> x0 = 1.0
1.0

julia> hT = Taylor1(5) # small quantity
 1.0 t + 𝒪(t⁶)

julia> x = x0 + hT
 1.0 + 1.0 t + 𝒪(t⁶)

julia> Y = log(x)
 1.0 t - 0.5 t² + 0.3333333333333333 t³ - 0.25 t⁴ + 0.2 t⁵ + 𝒪(t⁶)

Note that there is an overall multiplicative -1 with respect to what you wrote as Y(h). Now, let us check the result by evaluating Y at h=0.01, and compare it with log(1.01):

julia> Y(0.01)
0.009950330853333333

julia> log(1.01)
0.009950330853168092

Not the same, but a good approximation. So indeed evaluating Y at h approximates log(x0+h).

Let us obtain (an approximation of) the Taylor expansion of log(x) again, but now around another point z0=1.25, but exploiting the series Y we obtained before; I though that this was your motivation. My point is: beside Y we also require the value of x0, which we currently do not store in Y, so from Y alone, this is not possible. Assuming we know x0, the Taylor expansion we are after is Y(z0+h-x0), which is a translation of h:

julia> z0 = 1.25
1.25

julia> Y(z0+hT-x0) # Same result, updating the definition of Y, as `update!(Y, z0-x0)`
 0.22317708333333333 + 0.80078125 t - 0.3125 t² + 0.20833333333333331 t³ - 1.3877787807814457e-17 t⁴ + 0.2 t⁵ + 𝒪(t⁶)

julia> log(z0+hT) # for comparison
 0.22314355131420976 + 0.8 t - 0.32 t² + 0.1706666666666667 t³ - 0.1024 t⁴ + 0.06553600000000001 t⁵ + 𝒪(t⁶)

In these expressions t is the small quantity, now referred to z0 instead to x0. The coefficients of the series, while not identical, show a good correspondence. Interestingly, from the Taylor expansion at x0, we have obtained an approximation of some higher-order derivatives at z0.