Can this library be used to calculate a bunch of taylor coefficients in parallel?

[LeanderK]

I have a family of functions, that is really just the same function with different parameters (so it can be vectorized on the GPU), and I want to calculate higher-order taylor-coefficents of each multivariate function in parallel. Could this be achieved using this library?

I don't know julia enough to really judge this on my own, so I would appreciate the hint. Thanks!

In case this is confusing, I have a family of functions f'(x) = (f_1(x), ...., f_n(x)) that performs the same computation over different constants, I am interested in computing all multivariate taylor-coefficients (up to a certain order) in parallel. I work with these functions on the gpu so ideally a gpu-based solution would be nice but that's optional for the prototype!

[lbenet]

Thanks for openning this issue.

First, I will answer the expansion in Taylor of f_1(x, y). This is indeed possible, with some constrains related to the number of "independent variables" and order of the polynomial (since the actual number of polynomial coefficients grows pretty fast).

As an example (I'm starting Julia as julia -t auto, which in my machine means 8 threads):

julia> using TaylorSeries

julia> x, y = set_variables("δx δy")  # Set the global parameters for your expansion
2-element Vector{TaylorN{Float64}}:
  1.0 δx + 𝒪(‖x‖⁷)
  1.0 δy + 𝒪(‖x‖⁷)

julia> f_1(x, y) = cos(2*x) - sin(3*x*y) # some example for `f_1(x)`
f_1 (generic function with 1 method)

julia> f_1(0.5+x, 0.2+y) # expansion around 0.5 for `x` and around 0.2 for `y`
 0.24478209920680016 - 2.2561438630911566 δx - 1.433004733688409 δy - 1.0274109745372384 δx² - 2.6000412813816123 δx δy + 0.332460232494007 δy² + 1.1563534266857172 δx³ + 0.789877224054325 δx² δy + 1.9746930601358124 δx δy² + 0.5373767751331534 δy³ + 0.3586057281294553 δx⁴ + 0.49992361296811505 δx³ δy + 3.8494066087662437 δx² δy² + 3.1245225810507185 δx δy³ - 0.06233629359262632 δy⁴ - 0.22501132066039242 δx⁵ - 0.0396544078813421 δx⁴ δy + 2.3013460254338645 δx³ δy² + 5.753365063584661 δx² δy³ - 0.6196001231459701 δx δy⁴ - 0.060454887202479754 δy⁵ - 0.048007721923331886 δx⁶ - 0.015189205482960003 δx⁵ δy - 0.39234059337856575 δx⁴ δy² + 2.5277747311278875 δx³ δy³ - 2.4521287086160353 δx² δy⁴ - 0.5933283391781249 δx δy⁵ + 0.004675222019446973 δy⁶ + 𝒪(‖x‖⁷)

Is this the kind of things you are after?

Regarding paralelization, I do not know too much on it, and even less on GPUs, but perhaps the following helps

julia> x0 = -5:1/128:5 ; y0 = -2:1/128:2;

julia> @time begin # "serial" run, after running it twice
       res = Array{TaylorN{Float64}}(undef, length(x0), length(y0))
       for ny0 in eachindex(y0)
           for nx0 in eachindex(x0)
               res[nx0,ny0] = f_1(x0[nx0]+x, y0[ny0]+y)
           end
       end
       end
  8.148394 seconds (330.68 M allocations: 15.356 GiB, 22.35% gc time)

julia> @time begin
       res = Array{TaylorN{Float64}}(undef, length(x0), length(y0))
       @threads for ny0 in eachindex(y0)
           @threads for nx0 in eachindex(x0)
               res[nx0,ny0] = f_1(x0[nx0]+x, y0[ny0]+y)
           end
       end
       end
  2.491567 seconds (330.24 M allocations: 15.348 GiB, 59.57% gc time, 18.09% compilation time)

I hope all this is useful.