# Simple taylor serie expension that is rather slow

Hey,

I have a function `g`, which definition is the code below, from which I need to extract Taylor coefficients. I am using `TaylorSeries.jl` to do that but that is rather long… Here is the code :

``````using TaylorSeries, BenchmarkTools

function f(t,p,α,s)
x = 1.0 / (1.0 - s*(t+1)/(t-1))
rez = zero(x)
for i in eachindex(p)
rez += p[i] * x^α[i]
end
return rez * sqrt(2) / (1-t)
end
coefficients(m,pars...) = taylor_expand(t -> f(t,pars...),order=m).coeffs

#Random Parameters :
N = 100
m = 20
p,α,s = (rand(N), rand(N), rand())
p ./= sum(p)

# Test :
@btime coefficients(m,p,α,s);
``````

Which outputs on my machine:

`````` 75.100 μs (433 allocations: 92.14 KiB)
``````

Is there a way to make it faster and/or non-allocative ? Since I am using this inside an optimization routine, it is particularly sensitive…

If you want something to play with that isn’t quite ready for prime time, there is a new Taylor-mode automatic differentiation library being built by one of the Julia Lab students.

But it’s missing a lot of overloads right now, so YMMV and you might need to wait for the Symbolic codegen integration. But it’s already quite usable, so maybe worth a try. No docs at this stage of course.

TaylorSeries.jl just isn’t the right design for performance here.

Excellent! If this is much better, than it will improve AsymptoticNumericalMethod.jl which was already quite fast to me 1 Like

Thanks for the link, i took a look and it seems like a lot it still missing. Cool prpoject though.

But year this is the path i’m trying to take, passing `Symbolic.jl` variables for the parameters to precompile the taylor expansion. Without luck sofar, but i am still testing stuffs.

EDIT : The ideal stuff for me would be that the folliwng code works :

``````using TaylorSeries, Symbolics

function fi(t,α,s)
x = (t-1) / (t*(1-s) - (1+s))
return x^α * sqrt(2) / (1-t)
end

@variables t_sym α_sym s_sym

fi(Taylor1(Num,5),α_sym,s_sym)
``````

but right now it outputs :

``````ERROR: MethodError: *(::Taylor1{Num}, ::Num) is ambiguous. Candidates:
*(a::Number, b::Num) in Symbolics at C:\Users\lrnv\.julia\packages\SymbolicUtils\qulQp\src\methods.jl:71
*(b::Taylor1{S}, a::T) where {T<:Union{Real, Complex}, S<:Union{Real, Complex}} in TaylorSeries at C:\Users\lrnv\.julia\packages\TaylorSeries\YHybm\src\arithmetic.jl:231
*(b::Taylor1{T}, a::T) where T<:Number in TaylorSeries at C:\Users\lrnv\.julia\packages\TaylorSeries\YHybm\src\arithmetic.jl:239
Possible fix, define
*(::Taylor1{T}, ::T) where T<:Num
Stacktrace:
 fi(t::Taylor1{Num}, α::Num, s::Num)
@ Main .\Untitled-2:4
 top-level scope
@ Untitled-2:10
``````

@shashi was working with Songchen this week and built a prototype of this to generate optimal Taylor expansion code from DiffRules/ChainRules expressions (that can be traced, etc. etc.). So there’s an example code of it working, but I don’t think it’s in the repo yet. When that’s in, it should in theory be “as powerful as ForwardDiff”, with better (linear) scaling to higher derivatives. That said, it’s brand spanking new so of course, use with caution and without docs If you’re trying to do perturbation theory things, check out the new example in the Symbolics.jl docs:

https://symbolics.juliasymbolics.org/dev/examples/perturbation/

No i am not, I just need a fast function that output the taylors coefs of this one given the parameters But thanks i’ll take a look anyway.

Edit : w.r.t the ambiguity that i encounter, how can i tell the compiler to use the third option, as it clearly seems to be the one i want ?

You need to define that method. So define `*(::Taylor1{T}, ::T) where T<:Num = ...`. You’ll need to do that for a lot of things though.

Ok I found a solution using a branch of `TaylorSeries.jl` that had almost everything I needed. Take a look at this code :

``````# ] add TaylorSeries#lb/Symbolics
using Symbolics, TaylorSeries, BenchmarkTools

############### Parameters
N = 200
m = 10
p,α,s = (rand(N), rand(N), rand())
p ./= sum(p)

################ Old version
function f(t,p,α,s)
x = (t-1) / (t*(1-s) - (1+s))
rez = zero(x)
for i in eachindex(p)
rez += p[i] * x^α[i]
end
return rez * sqrt(2) / (1-t)
end
const tayl = Taylor1(m)
coefficients(args...) = f(tayl,args...).coeffs

##############"" New version
@variables t₀ αₛ sₛ;
tₛ = t₀ + Taylor1(m);
tₛ = tₛ - t₀;
function fi(t,α,s)
x = (t-1) / (t -t*s - 1-s)
return x^α / (1-t)
end

# This errors as the power is not defined.
# fi(tₛ,αₛ,sₛ)

# Let's define it, according to the definition in TaylorSeries.jl
import Base.^
^(a::Taylor1{T}, r::Num) where {T<:Number} = exp(r * log(a)) # quick dirty workaround.

# Let's now extract coefficients from the symbolic Taylor series, at t₀ = 0
symbolic_coefs = fi(tₛ,αₛ,sₛ).coeffs;
# Let's make a few more substitutions :
symbolic_coefs = substitute.(symbolic_coefs,(Dict(log(-1/(-1-sₛ)) => -log(1-sₛ)),));
symbolic_coefs = substitute.(symbolic_coefs,(Dict(exp(-αₛ*log(1 - sₛ)) => (1/(1-sₛ))^αₛ),));
#Compile
coefs_expr= build_function(symbolic_coefs,αₛ,sₛ);
taylor_fi = eval(coefs_expr);
coefficients2(p,α,s) = sqrt(2) * sum(pj .* taylor_fi(αj,s) for (pj,αj) in zip(p,α))

coefficients(p,α,s)
coefficients2(p,α,s)

@btime coefficients(p,α,s); # 52.500 μs (822 allocations: 114.52 KiB)
@btime coefficients2(p,α,s); # 101.800 μs (1601 allocations: 106.27 KiB)... Epic fail.
``````

However, there is still a lot of redundancy is the symbolic expression for the vector `symbolic_coefs`. Do you know how I could produce a function that exploits this redundancy, or is it just not possible yet ? I did a little by hand through the `substitute` calls, but this is not enough, clearly.

For example, the expression of `symbolic_coefs` contains a lot of `1 - sₛ`, `1 + sₛ` and fraction of those two values all over. How could I tell `Symbolics.jl` that, when compiling the function, it should probably pre-compute these two values, their inverses and their ratio to drastically reduce the number of divisions ? Is it even possible, or am I asking too much ?

1 Like

You want to generate it slightly differently and `simplify`

I’m Songchen working on TaylorDiff.jl which is designed to be performance-focused. Not sure if I’m late for this discussion, but for the function you provided TaylorDiff.jl is already working, super-fast and non-allocating:

``````using TaylorDiff, BenchmarkTools

function f(t, p, α, s)
x = 1.0 / (1.0 - s*(t+1)/(t-1))
rez = zero(x)
for i in eachindex(p)
rez += p[i] * x ^ α[i]
end
return rez * sqrt(2) / (1-t)
end

# Random Parameters
N = 100
m = 20
p, α, s = (rand(N), rand(N), rand())
p ./= sum(p)

# Test
t = TaylorScalar{Float64, 20}(0.)
@btime f(t, p, α, s);
``````

On my computer (MBP 2018):

``````\$ julia --project test.jl
26.565 μs (1 allocation: 176 bytes)
``````

Of course this is in early alpha stage, but tests are always passed before pushing to main so it’s somewhat usable. I’d be happy to get connected with use cases and probably add them to benchmarks of this project.

2 Likes

Indeed, this is fast, but it also seems to be wrong (?) :

``````using TaylorSeries, TaylorDiff, BenchmarkTools

function f(t,p,α,s)
x = (t-1) / (t*(1-s) - (1+s))
rez = zero(x)
for i in eachindex(p)
rez += p[i] * x^α[i]
end
return rez * sqrt(2) / (1-t)
end
coefficients(m,args...) = f(Taylor1(eltype(p),m),args...).coeffs
coefficients2(m,args...) = f(TaylorScalar{Float64, m+1}(0.),args...).value
# Random Parameters
N = 100
m = 20
p, α, s = (rand(N), rand(N), rand())
p ./= sum(p)

#Check values:
rez = [coefficients(m,p,α,s) [coefficients2(m,p,α,s)...]]
``````

outputs :

``````julia> rez = [coefficients(m,p,α,s) [coefficients2(m,p,α,s)...]]
21×2 Matrix{Float64}:
1.31636        1.31636
1.14296        1.14296
1.00866        1.00037
0.903913       0.865792
0.821562       0.703275
0.756259       0.429575
0.703992      -0.206211
0.661741      -2.10222
0.627229      -8.92891
0.598733     -37.5301
0.574942    -173.838
0.554859    -901.26
0.537718   -5196.52
0.522932  -32992.5
0.510046      -2.28587e5
0.498706      -1.71539e6
0.488635      -1.38559e7
0.479615      -1.19833e8
0.471476      -1.10462e9
0.464077      -1.081e10
0.457311      -1.11921e11

julia>
``````

Is there a difference in meaning between the two ? Maybe i am issing a simple transformation, like multiplication by a factorial ? (does not look like it).

Okay, I spotted an indexing error at primitive `/`, just fixed. And yes, you need a factorial to translate the results, since the `n + 1` th term in TaylorDiff is the `n` th derivative, not the coefficients

1 Like