I’m interested in using AD, and so have been something like 90% successful
using TaylorSeries.jl to generate the inverse function to almost any general arbitrary function.
For example TaylorSeries.jl can invert asin from sin , sin from asin etc. - but not,
strangely enough, log/ln for exp - might you have an insight as to why this MWE snippet below works for all except the exp function ?
using Micrologging
using TaylorSeries
# @warn "MSC NOTE TIP TaylorSeries.inverse — Function. https://www.juliadiff.org/TaylorSeries.jl/v0.5/api.html "
# MSC NOTE TIP inverse(f)
# Return the Taylor expansion of f^−1(t)
# , of order N = f.order, **for f::Taylor1 polynomial** **if the first coefficient of f is zero.** **Otherwise, an ArgumentError is thrown.**
# The algorithm implements Lagrange inversion at t=0
# if f(0)=0
# f−1(t)=∑n=1Ntnn!dn−1dzn−1(zf(z))n∣∣∣z=0.
sin(asin(0.9))
# 0.9
asin(sin(0.9))
# 0.9
t = Taylor1(15)
@debug "t = Taylor1(15) ; 1.0 t + 𝒪(t¹⁶)"
p = sin(t)
TaylorSeries.evaluate(inverse(p), TaylorSeries.evaluate(p, 0.9))
# 0.8999969178405351
tBig = Taylor1(BigFloat, 50) # Independent variable with BigFloats, order 50
# julia> # With order 50 precision
tBig = Taylor1(BigFloat, 50) # Independent variable with BigFloats, order 50
# 1.0 t + 𝒪(t⁵¹)
# julia> # TaylorSeries.evaluate(inverse(p), TaylorSeries.evaluate(p, 0.9))
TaylorSeries.evaluate(inverse(exp(tBig)), TaylorSeries.evaluate(exp(tBig), 0.9))
# DEBUG ERROR:
# ERROR: ArgumentError: **Evaluation of Taylor1 series at 0 is non-zero.**
# For high accuracy, revert a Taylor1 series with first coefficient 0 and re-expand about f(0).
At first glance it seems like inverting exp with Taylor Series would be the easiest - not the hardest, So hopefully you/someone on this thread might have some hints about that TY
The error message is giving you the clue: the value of exp(t) at t=0 is 1.
When you reflect this in the line y=x to find the inverse function, that function (log) cannot be expanded in a power series around t=0.
Also, please prefer to start a new thread for questions that are not really related to the original one.
As noted by the error message, to invert a series using Lagrange inversion formula you need it to have a zero constant term. The following I think is what you are looking for:
julia> using TaylorSeries
julia> t = Taylor1(5) # independent variable
1.0 t + 𝒪(t⁶)
julia> p = exp(t) # Taylor series of `exp` around zero
1.0 + 1.0 t + 0.5 t² + 0.16666666666666666 t³ + 0.041666666666666664 t⁴ + 0.008333333333333333 t⁵ + 𝒪(t⁶)
julia> inverse(p-1) # t -> t-1, so the new expansion is around 1
1.0 t - 0.5 t² + 0.3333333333333333 t³ - 0.25 t⁴ + 0.2 t⁵ + 𝒪(t⁶)
julia> log(1+t) # Taylor series of `log` around 1
1.0 t - 0.5 t² + 0.3333333333333333 t³ - 0.25 t⁴ + 0.2 t⁵ + 𝒪(t⁶)
Thank you for your help.
Here’s a summary of what we have figured out so far:
# Tip howto Use Taylor Series Polynomial to automatically calculate the Inverse function of an arbitrary function.
# This Example Code uses Taylor Series Polynomial inverse(exp) to calculate the natural log
# julia> CLI aka REPL
shift_taylor(a) = a + Taylor1(typeof(a),127)
p = exp(shift_taylor(0.0)) # Taylor series of `exp` around zero
exp(TaylorSeries.evaluate(inverse(p - 1.0), 1.0))
# 2.00785841499612
exp(TaylorSeries.evaluate(inverse(p - 1.0), 0.9))
# 1.9000000109002069
exp(TaylorSeries.evaluate(inverse(p - 1.0), 0.5))
# 1.5
exp(TaylorSeries.evaluate(inverse(p - 1.0), 0.0))
# 1.0
Next issue to tackle and some observations that may lead us to
potentially useful techniques to improve accuracy near endpoint extrema=2.0 .
Issue Example Code:
shift_taylor(a) = a + Taylor1(typeof(a),127)
p = exp(shift_taylor(0.0)) # Taylor series of `exp` around zero
exp(TaylorSeries.evaluate(inverse(p - 1.0), 1.0))
# 2.00785841499612
Issue Description: Despite Taylor1 Polynomial @very high order=127,
there is a reduction of accuracy/loss of significant figures near discontinuity endpoint=2.0
IOW : Taylor Series Polynomials are Subject to Runges Phenomenon / Phenomena aka Ringing at the endpoint(s) discontinuitie(s)
SOLN TIPs Visualize the issue and possible mitigations solutions
Details
Runge’s function is given by … When this function is sampled at equally spaced points in the range , as the number of sample points increase the error difference the interpolating polynomial and the function grow without bound.
When Cheybshev points are used the error diminishes as the number of points increases.
SOLN TIPs Operate strictly using continuous Fourier FFT in Frequency space.
Taylor Series ( as Truncated Polynomials ) for Trig functions such as sin(x) cos(x) are Subject to Gibbs Phenomena aka Ringing at the endpoint(s) discontinuitie(s)
In the polynomial interpolation setting, the Gibbs phenomenon can be >> mitigated using the S-Gibbs algorithm [14]. <<
A Python implementation of this procedure can be found here.
The Gibbs phenomenon is an overshoot (or “ringing”) of Fourier series and other eigenfunction series occurring at simple discontinuities. It can be >> reduced with the Lanczos sigma factor. <<
SOLN TIPs Possible Additional resources
Book: Numerical Methods that Work (Spectrum) 1st Edition classic originally published in 1970
by Forman S. Acton (Author)
Part II presents basic tools:
=============
Moving Code Compilation vs REPL CLI v1.04 issue here for now, probably a separate Post later .
Issue Example Code Description: When you execute ALL Three Lines AT ONCE in the ( v1.04 specific bug ?) REPL the answer is different than if you pause before running TaylorSeries.evaluate( , which points to a “race condition” which is a common programming issue described in Race Condition Issue Description:** after the example code here:
shift_taylor(a) = a + Taylor1(typeof(a),127)
p = exp(shift_taylor(0.0)) # Taylor series of `exp` around zero
exp(TaylorSeries.evaluate(inverse(p - 1.0), 1.0))
julia> exp(TaylorSeries.evaluate(inverse(p - 1.0), 1.0))
2.9131380935020
julia> exp(TaylorSeries.evaluate(inverse(p - 1.0), 1.0))
2.00785841499612
Race Condition Issue Description: One definition of race condition describes it as “anomalous behavior due to unexpected critical dependence on the relative timing of events.” Race Condition (Software)
There is a bunch of issues reported here. I’ll try to clarify some of them, while I get to understand the others. Thanks a lot for the reference. As for the race condition you report, I can’t check it now; please do open an issue in TaylorSeries.jl to follow this up.
What I understand from the code with exp and its inverse (log), and the fact that you are using a quite high order, is that you are interested in checking how the series behaves as you approach the radius of convergence. The following shows that the series does converge, though slowly:
for ord = 100:100:500
t = Taylor1(ord) # independent variable of order `ord`
p = exp(t)
q = inverse(p-1)
l2 = q(1.0) # equivalent to `evaluate(q, 1.0)`
@show(ord, l2-log(2.0))
println()
end
ord = 100
l2 - log(2.0) = -0.004975001249750366
ord = 200
l2 - log(2.0) = -0.0024937500781213595
ord = 300
l2 - log(2.0) = -0.0016638889043206762
ord = 400
l2 - log(2.0) = -0.0012484375048821272
ord = 500
l2 - log(2.0) = -0.0009990000019985956
As it can be seen, the difference of evaluating the Taylor polynomial to log(2) gets smaller as the order is increased.
I should clarify that evaluate uses Horner’s method, so there is no interpolation for computing the evaluation of the polynomial itself, nor in the computation of the coefficients of the series. Certainly, there is a loss of accuracy, which is due to the fact that the coefficients have finite precision, and that you are evaluating the series at the radius of convergence.
Thank You Luis @lbenet I believe it is critical to be able to inverse functions so it’s important not to “lose information”, so suggest consider that in light of “… should lower the order. Yet, to change this may have other consequences due to #226.” here @@ https://github.com/JuliaDiff/TaylorSeries.jl/issues/230#issuecomment-808425320 … and so Luis @lbenet whatever changes you make (if any), I request you regression test against the presently working use case above where “… the TaylorSeries.jl code works with exp and its inverse (log)” And of course I would be happy to write that code into a PR if that helps.
Dear Marc, I agree with you that inverse functions are important and that we don’t want to loose any information, but I also agree (with others) that mathematical correctness should be preserved, and that at the moment we are not being fully consecuent. I’ll keep in mind your post, and before merging TaylorSeries.jl/pull/248 we will do some thorough testing, including the case used above, which I think it will not be affected by the PR.
I’ll also would like to invite you to participate in the discussion of that issue and associated PR, so we are all happy with the outcome. Finally, due to lack of time I haven’t been able to complete that PR, but I am truly willing to do.
At first glance it seems like inverting exp with Taylor Series would be the easiest - not the hardest, So hopefully you/someone on this thread might have some hints about that TY 