Trying to understand the math of the power function

[tansongchen]

Hi there,

I am reading the documentation of this package which mentioned several recurrence formula for elementary function. However, I couldn't fully understand how they are derived.

https://juliadiff.org/TaylorSeries.jl/latest/background/#Elementary-functions-of-polynomials

For example, could you demonstrate in more detail how you derived

from differential equations theory?

[lrnv]

This is simply Faa di bruno formula no ? Edit: It would work but does not correspond to the differential equation derivation proposed in the docs.

[tansongchen]

From my understanding Faa di Bruno's formula simply decompose the n-th derivative of $f(g(x))$ to combinations of various order derivative of $f$ and $g$. However, this is too complicated to implement generally, and to simplify one need specific recurrence relations for a specific $f$

[tansongchen]

There is a reference (mentioned in the README of this repo) also talked about differential equations, but is equally obscure that I didn't understand either :)

[lbenet]

Hi there,

I am reading the documentation of this package which mentioned several recurrence formula for elementary function. However, I couldn't fully understand how they are derived.

https://juliadiff.org/TaylorSeries.jl/latest/background/#Elementary-functions-of-polynomials

Hi!

Thanks for asking! Yes, actually it is not difficult. (While I was writing, you were much faster and quoted the reference; example 2.5 is truly useful!)

The first thing, is to propose an initial value problem. If you take the derivative of $p(x) = (f(x))^\alpha$ wrt $x$ (the independent variable), you can arrive to the ODE $p'(x) = \alpha p(x) f'(x)/f(x)$, with initial condition $p(x_0) = (f(x_0))^\alpha$.

You write $f(x) = \sum_k f_{[k]}(x_0) (x-x_0)^k$ and $p(x) = \sum_k p_{[k]}(x_0) (x-x_0)^k$, the Taylor series of $f(x)$ and $p(x)$ around $x_0$; here, the coefficients $f_{[k]}(x_0)$ are assumed to be known, while the coefficients $p_{[k]}(x_0)$ are to be determined. Note that $p_{[0]}(x_0) = p(x_0) = (f(x_0))^\alpha$ is the initial condition. You now plug in the expansions in the differential equation, and exploit the product and division rules of polynomials. If you solve for the $k$-th coefficient of $p$, you obtain the desired formula.

[lbenet]

I hope this minimalistic explanation helps!

[lbenet]

Incidentally, much easier is to derieve the formula for $\exp(f(x))$...

[lrnv]

Also relevant, a derivation of $exp(f(x_1,...,x_n))$ can be found there https://arxiv.org/pdf/1911.11722

[tansongchen]

Thank you so much for this explanation! I finally managed to get the derivation.

My method was: write the ODE as $fp'=\alpha pf'$, insert the expansions, and relate the coefficient of $x^{n-1}$ on both sides, assuming that the expansion on $f$ and $p$ is up to order $x^n$. This is a bit weird, because originally I thought this should be like relating the coefficient of $x^n$ on both sides...

But anyway, I really appreciate that.

[tansongchen]

Another question is the handling of ^(a, b::Real) with zeroth order coefficient to be 0. At https://github.com/JuliaDiff/TaylorSeries.jl/blob/master/src/power.jl#L271 , it seems that some extra handling is needed to make this mathematically correct. So there are two questions:

1. Could you explain what is the high-level idea to handle this, and how to reason about the corresponding part of pow!
2. Line 271 requires that r * l0 to be an integer, which implies I cannot do that:

julia> a = Taylor1(Float64, 1)
 1.0 t + 𝒪(t²)

julia> a^1.5
ERROR: DomainError with  1.0 t + 𝒪(t²):
The 0-th order Taylor1 coefficient must be non-zero
to raise the Taylor1 polynomial to a non-integer exponent.

But mathematically $\mathrm d(x^{3/2})=(3/2) x^{1/2}\mathrm dx$, so this is defined at 0. Also

julia> ForwardDiff.derivative(x -> x^1.5, 0.)
0.0

So could this be made more general?

[lbenet]

1. Could you explain what is the high-level idea to handle this, and how to reason about the corresponding part of `pow!`

2. Line 271 requires that `r * l0` to be an integer, which implies I cannot do that:

Thanks again for asking! Perhaps you are right and we may need to handle these cases more carefully. Yet, I think we do it right (and consistently wrt other programs doing the same, i.e. gp/pari). The crutial point is to note that the power function with non-integer power has issues at 0: it is discontinuous there (and its derivatives too!), and thus we cannot apply Taylor's theorem.

Let me use your example, $f(x)=x^{1.5}$, then $f'(x) = 1.5 x^{0.5}$, $f''(x)=0.75*x^{-0.5}$, etc. So, if we would like to construct the Taylor series around $x_0=0$, we have $f(x_0)=0=f'(x_0)$, as you correctly note, but $f''(x_0)$ is undefined, and therefore we have no Taylor expansion. Yet, if you ensure continuity (and also of its derivatives), which requires for this $f(x)$ to be non-zero, then you can use Taylor's theorem and:

julia> using TaylorSeries

julia> t = Taylor1(5)
 1.0 t + 𝒪(t⁶)

julia> (0.1+t)^1.5 # expansion of f(x) around 0.1
 0.0316227766016838 + 0.47434164902525694 t + 1.1858541225631423 t² - 1.9764235376052368 t³ + 7.4115882660196375 t⁴ - 37.05794133009819 t⁵ + 𝒪(t⁶)

Very naively, and perhaps too cryptically, this is encoded in the line you quoted. Going into the details, l0 is the first nonzero coefficient of the Taylor expansion of the argument of the pow function, i.e., of $f(x)$, and hence an Int, and r::Real is the power. Then, r*l0 is integer if l0==0 (so the constant term is non-zero, and you can apply Taylor's theorem), or if the constant-term is zero, then r must be an integer.

Since you are playing with this function, please please, have a look into #369. I'm trying to modify how the pow function is handled, in particular to reduce allocations.