ForwardDiff of in-place functions?

[ChrisRackauckas]

Is there any way to ForwardDiff in-place functions? Most of mine look like this:

g = (du,x,t)-> du[:]=f(x,t))

for some f not in place. Is my best bet just to use f directly and then build an in-place derivative of g from the derivative of f?

[KristofferC]

http://www.juliadiff.org/ForwardDiff.jl/api.html#jacobians-of

Second section. Is that what you want?

[ChrisRackauckas]

That works in some of cases, but not for all. The best fix would be to find some way to pre-allocate du as a vector of the dual numbers that will be coming out of f (in order to use in conjunction with NLsolve). Is that possible?

[jrevels]

That works in some of cases, but not for all.

Can you give an example where the provided ForwardDiff.jacobian! doesn't work? For the example you gave in the OP:

du, x, t = ... # define some values
g = (du,x,t)-> du[:]=f(x,t))

# jacobian w.r.t. to x
out = similar(du, length(du), length(x))
ForwardDiff.jacobian!(out, (du, x) -> g(du, x, t), du, x)

# jacobian w.r.t. to t
out = similar(du, length(du), length(t))
ForwardDiff.jacobian!(out, (du, t) -> g(du, x, t), du, t)

[ChrisRackauckas]

I can't do that inside of NLsolve. A pre-allocated du would work though, but I don't know how to make that happen.

[ChrisRackauckas]

Nevermind, I think I thought of a way to make it work.

[jrevels]

I can't do that inside of NLsolve.

Can you provide a minimal working example that demonstrates the limitation you ran into?

Nevermind, I think I thought of a way to make it work.

If you wouldn't mind sharing, I'd be interested in hearing it. We might be able to use your feedback to make the package better/easier to use (or at least improve the documentation).

[ChrisRackauckas]

There is just a problem with allowing for closures. For example, in your example, what if only a function of (x,t) is exposed via some function? Then you would want to use the closure (x,t)->g(du,x,t) as the function call. But it will fail unless you can pre-allocate du.

[jrevels]

I'm not totally sure I understand without an example, but it kind of sounds like the problem is that you're given a non-in-place function f(x,t) from an external source, and you want to take it's Jacobian. If that's the case, why not just do ForwardDiff.jacobian!(out, x -> f(x, t), x) (or the equivalent w.r.t. to t)?

[ChrisRackauckas]

I can do that by modifying their source code. But that seems like a heavy-handed solution for not being able to pre-allocate a vector of duals...

[jrevels]

I can do that by modifying their source code.

Okay, I think I'm misunderstanding your problem - in the example I'm talking about, no source code modification is needed. You are given the differentiable function f(x,t).

not being able to pre-allocate a vector of duals...

I'm confused as to why you wish to do this when ForwardDiff is already doing this kind of work for you (or at least, it should be). If for some reason you really did want to allocate a vector of Duals, it's as easy as calling Vector{ForwardDiff.Dual{N,T}}(n).

Without a code example so that I can understand the problem, I don't think I'm going to be able to help much here.

[ChrisRackauckas]

I was trying to clean up the examples, but I have to go teach. The code is here:

https://github.com/ChrisRackauckas/DifferentialEquations.jl/blob/master/src/ode/ode_integrators.jl#L618

You can see the Implicit Euler code can be fixed by better use of resid which is the only allocated vector of duals (tried and it works). This doesn't generalize very well, as you can see in the trapezoid method which uses a closure with du1 and du2, but I need those instead allocated as the "right vector of duals" for that to work. In general, for longer functions, you need to do it to cut out allocations anyways.

[KristofferC]

Btw, just pointing out that lines like: https://github.com/ChrisRackauckas/DifferentialEquations.jl/blob/master/src/ode/ode_integrators.jl#L588 will have very bad performance due to being lowered to +(a,b,c,d...) and then having too many arguments and get splatted.

f(a) = a[1] + a[2] + a[3] + a[4] + a[5] + a[6] + a[7] + a[8]
g(a) = a[1] + a[2] + a[3] + a[4]

f is like 20 times slower here because it has too many arguments to +. To solve it you need to add paranthesis to group up the expressions in like 4 by 4. Ugly I know... There is an issue about it somewhere.
I wrote a macro to do it automatically...

[jrevels]

Thanks for the link. I think my confusion came from not realizing that you were working in a situation where you were handling Duals directly.

Since this code is manually injecting dual numbers into a separate differentiable context, I just to make sure you're aware of perturbation confusion. In other words, the dual numbers you manually construct will mix with the dual numbers injected when NLSolve performs autodiff, possibly leading to the wrong result.

I'm still a bit confused by what this code is doing. For example, take this line:

du = map((x)->ForwardDiff.Dual{length(u),typeof(x)}(x,length(u)),zeros(u))

This will never work for length(u) != 1, since you're not passing in the same number of perturbations specified by the type parameter. Even if it did work, it would still fail for length(u) > 10, since the Dual type is limited to 10 simultaneous perturbations for performance reasons (i.e. more than that can cause thrashing). I also don't understand why the dual number's perturbation should be seeded with the length of the input vector...

What is the motivation for injecting dual numbers manually before passing the function to NLSolve? AFAICT, du is never used after the initial differentiation, so the only effect of making du's elements dual numbers is to induce perturbation confusion in NLSolve's autodiff code.

[ChrisRackauckas]

f is like 20 times slower here because it has too many arguments to +. To solve it you need to add paranthesis to group up the expressions in like 4 by 4. Ugly I know... There is an issue about it somewhere.

Thanks for pointing that out. I'll go find the issue. I assumed the compiler would turn those into a bunch of FMA's. Maybe it's best just to loop over the appropriate indices tmp[:]+=b[i]*k[i], that would let the compiler optimize it, right? That's a different story though.

Since this code is manually injecting dual numbers into a separate differentiable context, I just to make sure you're aware of perturbation confusion. In other words, the dual numbers you manually construct will mix with the dual numbers injected when NLSolve performs autodiff, possibly leading to the wrong result.

I'm not making Dual numbers that mean anything... just trying to make du a pre-allocated vector of Duals so that it doesn't make the vector each time. Obviously using du[:] = f(u,t) works, but I just want to do f(du,u,t) for better performance... and all that it requires is for du to be an appropriately typed array.

And obviously the code that doesn't work doesn't work, I said to look at du1 and du2 of the trapezoidal method. The Euler method only needs one temporary around, so it will work without the full closure. The new master has the Euler method fully in-place using the 1 temporary that's given, resid. The trapezoidal rule can't do that, so unless there's a way to properly pre-allocate an array of Duals and push it in via the closure, it's not going to work (without having a temporary).

[KristofferC]

Dual vectors are cached inside ForwardDiff for a given length and element type.

Why are you not just passing du in? You are differentiating with respect to u, right?