Fast higher order differentiation

[Young Chun]

I have an optimization problem that involves series of high order derivatives of f(x). 

So, to get a decent value, I need to calculate higher order derivative at least to 40~50th order. 

Initially, I tried to use ForwardDiff package with automatic differentiation but kept facing the following error message

<div id="inner-editor"><br class="Apple-interchange-newline">using DualNumbers using ForwardDiff f(x::Vector) = sum(sin, x) + prod(tan, x) * sum(sqrt, x); x = rand(5) g = ForwardDiff.derivative(f);</div>

using DualNumbers

using ForwardDiff

f(x::Vector) = sum(sin, x) + prod(tan, x) * sum(sqrt, x);

x = rand(5)

g = ForwardDiff.derivative(f);

LoadError: MethodError: `derivative` has no method matching derivative(::Function)
Closest candidates are:
  derivative(::Any, !Matched::Any)
while loading In[28], in expression starting on line 4

So, I moved on to D(f) instead from Root package and it works quite well for a small order (like 10th order derivative)

<div id="inner-editor"><br class="Apple-interchange-newline">using DualNumbers using ForwardDiff using Roots using ForwardDiff using JuMP # Need to say it whenever we use JuMP using Gadfly function hder(f, n::Real) # subroutine to calculate higher order derivative temp = f; if n==0 return temp else for i=1:1:n temp = D(temp) end return temp end end</div>

using DualNumbers

using ForwardDiff

using Roots

using ForwardDiff

​

using JuMP  

using Gadfly

​

function hder(f, n::Real) # subroutine to calculate n-th order derivative 

    temp = f;

    if n==0

        return temp 

    else

        for i=1:1:n

            temp = D(temp)

        end

        return temp

    end

end

Out[7]:

hder (generic function with 1 method)

f(x) = (cos(x)) * exp(-1/5 * x)
g0 = hder(f, 0)
g1 = hder(f, 1)
g2 = hder(f, 2)
g3 = hder(f, 3)
g4 = hder(f, 4)
g5 = hder(f, 5)
g6 = hder(f, 6)
g7 = hder(f, 7)
g8 = hder(f, 8)

plot([g0, g1, g2, g3, g4, g5, g6, g7, g8], 0, 5pi)

However, if I push the limit and ask to calculate, like 30th order derivative, the program never ends and keep calculating for hours.

Is there a better way to do this type of task? either by using ForwardDiff package or by modifying my function?

I just learned Julia last week, so please understand if my question sounds stupid. 

[Jeffrey Sarnoff]

 I do not have your answer, your  question is a good one.

[David P. Sanders]

El domingo, 17 de julio de 2016, 17:13:38 (UTC+2), Young Chun escribió:

I have an optimization problem that involves series of high order derivatives of f(x). 

So, to get a decent value, I need to calculate higher order derivative at least to 40~50th order. 

I suggest you use TaylorSeries.jl.

Note that ForwardDiff from version 0.2 removed ForwardDiff.derivative(f); 

this must now be written as the anonymous function x -> ForwardDiff.derivative(f, x).

[David P. Sanders]

Here is how to use TaylorSeries.jl. In the next version, the syntax taylor1_variable will change to just Taylor1:

julia> t = taylor1_variable(40)

 1.0 t + 𝒪(t⁴¹)

julia> f(x) = sin(exp(x))

f (generic function with 1 method)

julia> f(t)

 0.8414709848078965 + 0.5403023058681398 t - 0.15058433946987837 t² - 0.4207354924039483 t³ - 0.3229307265911699 t⁴ - 0.13861923299172246 t⁵ - 0.016963650188307994 t⁶ + 0.027151311649628966 t⁷ + 0.028082409487798127 t⁸ + 0.016534675299332297 t⁹ + 0.006744265666861196 t¹⁰ + 0.0016199249536140723 t¹¹ - 0.00020257707395350685 t¹² - 0.0005026787959383626 t¹³ - 0.00034829060136926197 t¹⁴ - 0.00016780164241401782 t¹⁵ - 6.0803736936094484e-5 t¹⁶ - 1.4645171902851875e-5 t¹⁷ - 1.1346823483603049e-8 t¹⁸ + 2.515198118591517e-6 t¹⁹ + 1.817031683619697e-6 t²⁰ + 8.873827219817266e-7 t²¹ + 3.383697539084375e-7 t²² + 9.90873953752995e-8 t²³ + 1.7331465080481204e-8 t²⁴ - 2.8377431603022023e-9 t²⁵ - 4.478125328289587e-9 t²⁶ - 2.633876550034189e-9 t²⁷ - 1.1444925030062194e-9 t²⁸ - 4.013808037099447e-10 t²⁹ - 1.1102527569232767e-10 t³⁰ - 1.95623844080032e-11 t³¹ + 1.8017173073914597e-12 t³² + 3.715400951669623e-12 t³³ + 2.1778749843678503e-12 t³⁴ + 9.347841028230405e-13 t³⁵ + 3.279923193009391e-13 t³⁶ + 9.423306246867712e-14 t³⁷ + 1.997073344234444e-14 t³⁸ + 1.3114196619133462e-15 t³⁹ - 1.5394356646368331e-15 t⁴⁰ + 𝒪(t⁴¹)

julia> get_coeff(f(t), 40)

-1.5394356646368331e-15    

Note that this coefficient is the 40th derivative divided by 40!.  (The nth Taylor coefficient is f^{(n)} / n!. )

[Young Chun]

Thank you for the comment, David!

I changed the source code as below and it is really fast indeed. 

function hder(f, n, a)
    t = taylor1_variable(n);
    g(x) = f(x+t);
    hod = get_coeff(g(a), n) * gamma(n+1);
    return hod
end

f30(x) = hder(f, 30, x)

f40(x) = hder(f, 40, x)
f60(x) = hder(f, 60, x)

Now it takes about 0.000155 seconds to calculate 40th or 60th order derivatives. 

p.s. I learned Julia from your tutorial video at Scipy 2014 :)

       Thanks for all of the help.  

[Eric Forgy]

Hi Young Chun,

Welcome to Julia! Like Jeffrey, I think your question was a good one and the outcome is even better.

I don't know how others feel, but as an "oldie", it is super encouraging to me that:

1. People can even ask questions like this
2. Julia provides a highly satisfactory solution

When I was a kid, we had to walk uphill to in the snow to school both ways and we could not ask such questions. Even if we could ask, the answer would be less than satisfactory. The possibilities the young scientists today are provided is tremendous and I look forward to seeing what you and others come up with.

Out of curiosity, is anyone using Julia for fast polynomial approximations such as fast-multipole methods in gravitation/electromagnetics? This kind of fast auto-differentiation seems perfectly suited for these large-scale problems.

Cheers!

Eric

[Kristoffer Carlsson]

Take a look at ReverseDiffSource.jl.