# Simple taylor serie expension that is rather slow

**URL:** <https://discourse.julialang.org/t/simple-taylor-serie-expension-that-is-rather-slow/91677>\
**Category:** New to Julia\
**Created:** [December 15, 2022, 6:45am UTC](https://discourse.julialang.org/t/simple-taylor-serie-expension-that-is-rather-slow/91677 "2022-12-15T06:45:54Z")\
**Posts on this page:** 13\
**Page:** 1

<div class="post-metadata">

**Author:** ![lrnv](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lrnv/32/19373_2.png) [@lrnv](https://discourse.julialang.org/u/lrnv)\
**Post date:** [December 15, 2022, 6:45am UTC](https://discourse.julialang.org/t/simple-taylor-serie-expension-that-is-rather-slow/91677/1 "2022-12-15T06:45:54Z")

</div>

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 :

```julia
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:

```julia
 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…

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [December 15, 2022, 12:27pm UTC](https://discourse.julialang.org/t/simple-taylor-serie-expension-that-is-rather-slow/91677/2 "2022-12-15T12:27:35Z")

</div>

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.

> **[GitHub - JuliaDiff/TaylorDiff.jl: Taylor-mode automatic differentiation for...](https://github.com/JuliaDiff/TaylorDiff.jl)**
>
> Taylor-mode automatic differentiation for higher-order derivatives - GitHub - JuliaDiff/TaylorDiff.jl: Taylor-mode automatic differentiation for higher-order derivatives

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.

---

<div class="post-metadata">

**Author:** ![rveltz](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rveltz/32/2707_2.png) [@rveltz](https://discourse.julialang.org/u/rveltz)\
**Post date:** [December 15, 2022, 12:41pm UTC](https://discourse.julialang.org/t/simple-taylor-serie-expension-that-is-rather-slow/91677/3 "2022-12-15T12:41:56Z")

</div>

Excellent! If this is much better, than it will improve [AsymptoticNumericalMethod.jl](https://github.com/bifurcationkit/AsymptoticNumericalMethod.jl) which was already quite fast to me 😅

---

<div class="post-metadata">

**Author:** ![lrnv](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lrnv/32/19373_2.png) [@lrnv](https://discourse.julialang.org/u/lrnv)\
**Post date:** [December 15, 2022, 12:43pm UTC](https://discourse.julialang.org/t/simple-taylor-serie-expension-that-is-rather-slow/91677/4 "2022-12-15T12:43:32Z")

</div>

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 :

```julia
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 :

```julia
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:
 [1] fi(t::Taylor1{Num}, α::Num, s::Num)
   @ Main .\Untitled-2:4
 [2] top-level scope
   @ Untitled-2:10

```

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [December 15, 2022, 12:46pm UTC](https://discourse.julialang.org/t/simple-taylor-serie-expension-that-is-rather-slow/91677/5 "2022-12-15T12:46:55Z")

</div>

> [@lrnv](#):
>
> 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.

@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 😅

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [December 15, 2022, 12:48pm UTC](https://discourse.julialang.org/t/simple-taylor-serie-expension-that-is-rather-slow/91677/6 "2022-12-15T12:48:08Z")

</div>

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/](https://symbolics.juliasymbolics.org/dev/examples/perturbation/)

---

<div class="post-metadata">

**Author:** ![lrnv](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lrnv/32/19373_2.png) [@lrnv](https://discourse.julialang.org/u/lrnv)\
**Post date:** [December 15, 2022, 12:50pm UTC](https://discourse.julialang.org/t/simple-taylor-serie-expension-that-is-rather-slow/91677/7 "2022-12-15T12:50:02Z")

</div>

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 ?

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [December 15, 2022, 12:57pm UTC](https://discourse.julialang.org/t/simple-taylor-serie-expension-that-is-rather-slow/91677/8 "2022-12-15T12:57:12Z")

</div>

> [@lrnv](#):
>
> 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.

---

<div class="post-metadata">

**Author:** ![lrnv](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lrnv/32/19373_2.png) [@lrnv](https://discourse.julialang.org/u/lrnv)\
**Post date:** [December 15, 2022, 3:55pm UTC](https://discourse.julialang.org/t/simple-taylor-serie-expension-that-is-rather-slow/91677/9 "2022-12-15T15:55:23Z")

</div>

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

```julia
# ] 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[1]);
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 ?

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [December 15, 2022, 4:51pm UTC](https://discourse.julialang.org/t/simple-taylor-serie-expension-that-is-rather-slow/91677/10 "2022-12-15T16:51:08Z")

</div>

You want to generate it slightly differently and `simplify`

---

<div class="post-metadata">

**Author:** ![tansongchen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tansongchen/32/17310_2.png) [@tansongchen](https://discourse.julialang.org/u/tansongchen)\
**Post date:** [December 16, 2022, 4:16pm UTC](https://discourse.julialang.org/t/simple-taylor-serie-expension-that-is-rather-slow/91677/11 "2022-12-16T16:16:19Z")

</div>

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:

```julia
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):

```bash
$ 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.

---

<div class="post-metadata">

**Author:** ![lrnv](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lrnv/32/19373_2.png) [@lrnv](https://discourse.julialang.org/u/lrnv)\
**Post date:** [December 16, 2022, 4:33pm UTC](https://discourse.julialang.org/t/simple-taylor-serie-expension-that-is-rather-slow/91677/12 "2022-12-16T16:33:13Z")

</div>

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

```julia
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
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).

---

<div class="post-metadata">

**Author:** ![tansongchen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tansongchen/32/17310_2.png) [@tansongchen](https://discourse.julialang.org/u/tansongchen)\
**Post date:** [December 16, 2022, 4:52pm UTC](https://discourse.julialang.org/t/simple-taylor-serie-expension-that-is-rather-slow/91677/13 "2022-12-16T16:52:43Z")

</div>

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

> <https://github.com/JuliaDiff/TaylorDiff.jl/commit/0aca25e275fd84c87b52f9f3b3cd7e058d64d03a>