Autodifferentiation of a recursive function

Specific Domains Optimization (Mathematical)
, , ,
21

Ok, so using Memoize.jl it’s not too bad until you get to the clincher… If you’re ok with a finite difference approximation, then it’s doable… let’s see what happens.

2 Likes
22

Using the hard stuff here, already;)

23

But, OTOH derivatives are somehow recursive too…

That’s what the finite differences approach leverages to get O(N) :wink:

1 Like
24

Yes, it seems to be the way. I’ll will implement it and share the results here.

25

No need, start with this:

using Memoize

@memoize function wanddw(n,t,p)

    if n < 0
        error("wanddw: n must be an integer > 0")
    end
    if n == 0
        if t < -1
            return (complex(0.0),complex(0.0))
        elseif t <= 1
            return ((exp(p+p*t) - 1.0)/p, (t+1)*exp(p+p*t)/p - (exp(p+p*t)-1)/p^2)
        else
            return (2*exp(p*t)*sinh(p)/p, 2*sinh(p)*t*exp(p*t)/p - 2*sinh(p)*exp(p*t)/p^2 +2*cosh(p)*exp(p*t)/p)
        end
    elseif n == 1
        t < -1 && return (complex(0.),complex(0.))
        t <= 1 && return (-(1.0+exp(p+p*t)*(p-1)+p*t)/p^2, (-(p-1.0)*(t+1)*exp(p*t+p)-exp(p*t+p)-t)/p^2 - 2*((-(p-1)*exp(p*t+p)-p*t-1)/p^3))
        return (-2.0*exp(p*t)*(p*cosh(p)-sinh(p))/p^2, -(2.0*(p*cosh(p)-sinh(p))*t*exp(p*t)/p^2)-2.0*sinh(p)*exp(p*t)/p - 4.0*(p*cosh(p)-sinh(p))*exp(p*t)/p^3)
    elseif n >= 2
        (wprev,dwprev) = wanddw(n-1,t,p)
        (wprev2,dwprev2) = wanddw(n-2,t,p)
        (wprev2d, dwprev2d) =  wanddw(n-2,t,p+complex(0.0125,.0125))
        return (2*t*wprev - 2*dwprev-wprev2, 2*t*dwprev - (dwprev2d-dwprev2)/complex(0.0125,0.0125))
    end
    
end

The finite difference is hard coded, so you’ll want to fix that, but it returns IMMEDIATELY.

I used Maxima to get the symbolic derivatives for the base cases (because I’m lazy).

2 Likes
26

Thanks! I like the idea of returning both the value and the derivative in the same function for extra memoization efficiency!

I wonder if I can get this to work in terms of accuracy.

Thank you very much again.

27

It was fun! I think you just add the finite difference step size as a parameter and then use that in the code. Let us know how it works out, show some graphs or something.

I did that, and then tried it like this:


julia> wanddw(4,.1,complex(-.1),.0005)
(0.2959821757039954 - 0.0015914894716217987im, 0.4360000703131307 + 0.002182010051928633im)

julia> wanddw(4,.1,complex(-.1),.00005)
(0.2959821757039954 - 0.0015914894716217987im, 0.43592222276367454 + 0.0021033784685235563im)

julia> wanddw(4,.1,complex(-.1),.000005)
(0.2959821757039954 - 0.0015914894716217987im, 0.43591441315285184 + 0.002095691848207612im)

so at least as the step size gets smaller, it isn’t doing anything much

here’s the code with the dp:

using Memoize

@memoize function wanddw(n,t,p,dp)

    if n < 0
        error("wanddw: n must be an integer > 0")
    end
    if n == 0
        if t < -1
            return (complex(0.0),complex(0.0))
        elseif t <= 1
            return ((exp(p+p*t) - 1.0)/p, (t+1)*exp(p+p*t)/p - (exp(p+p*t)-1)/p^2)
        else
            return (2*exp(p*t)*sinh(p)/p, 2*sinh(p)*t*exp(p*t)/p - 2*sinh(p)*exp(p*t)/p^2 +2*cosh(p)*exp(p*t)/p)
        end
    elseif n == 1
        t < -1 && return (complex(0.),complex(0.))
        t <= 1 && return (-(1.0+exp(p+p*t)*(p-1)+p*t)/p^2, (-(p-1.0)*(t+1)*exp(p*t+p)-exp(p*t+p)-t)/p^2 - 2*((-(p-1)*exp(p*t+p)-p*t-1)/p^3))
        return (-2.0*exp(p*t)*(p*cosh(p)-sinh(p))/p^2, -(2.0*(p*cosh(p)-sinh(p))*t*exp(p*t)/p^2)-2.0*sinh(p)*exp(p*t)/p - 4.0*(p*cosh(p)-sinh(p))*exp(p*t)/p^3)
    elseif n >= 2
        (wprev,dwprev) = wanddw(n-1,t,p)
        (wprev2,dwprev2) = wanddw(n-2,t,p)
        (wprev2d, dwprev2d) =  wanddw(n-2,t,p+complex(dp,dp))
        return (2*t*wprev - 2*dwprev-wprev2, 2*t*dwprev - (dwprev2d-dwprev2)/complex(dp,dp))
    end
    
end

Even order 100 returns immediately, thanks to Memoization.jl !

julia> wanddw(100,.1,complex(-.1),.000005)
(6.329687336395586e86 - 2.2156952437460406e86im, -3.069266938528541e90 + 4.369343530297239e90im)

28

For what it’s worth, the above code is not numerically stable (and contains an error). I found another way to calculate the convolution analytically and use the recursive function to define the partial derivative with respect to p. Thanks again for the help!

1 Like