Simple taylor serie expension that is rather slow

New to Julia
1

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…

2

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.

3

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

1 Like
4

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

@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 :sweat_smile:

6

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/

7

No i am not, I just need a fast function that output the taylors coefs of this one given the parameters :slight_smile: 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 ?

8

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.

9

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[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 ?

1 Like
10

You want to generate it slightly differently and simplify

11

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
12

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

13

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