Automatic differentiation of mapping from R^N to R^N?

[Patrick Kofod Mogensen]

So I have a map from a set of strategies of players to the best response of each player given the original strategies (a discrete probability distribution for each player). I want to find the jacobian of this mapping, because I need to calculate the spectral radius of it. I am using finite differences at the moment, but I would like to use automatic differentiation.

Is this possible? I can't seem to find a JuliaDiff package that does it. Should I just do it for each element of my best-response mapping one at a time? So I get N functions from R^N to R ? (really it's from (0,1)^N to (0,1), but...)

Best,
Patrick

[Kristoffer Carlsson]

https://github.com/JuliaDiff/ForwardDiff.jl does this.

I use it all the time to compute jacobians of R^N -> R^N functions (the jacobian of the residual equations in non linear material modeling).

[Patrick Kofod Mogensen]

Thanks for the hint, I will have to look at bit closer at this. I can run my function, say phi(), on a vector, say P, and get the proper best response out. But when I try to use ForwardDiff.Jacobian, I get

ERROR: InexactError()
 in setindex! at array.jl:313
 in copy! at abstractarray.jl:307
 in phi at none:29
 in _calc_jacobian at /home/pkm/.julia/v0.4/ForwardDiff/src/api/jacobian.jl:101
 in jacobian at /home/pkm/.julia/v0.4/ForwardDiff/src/api/jacobian.jl:27

This strikes me a kind of weird, since the original method works. Line 29 assigns a vector into a vector of vectors.

[Kristoffer Carlsson]

You are most likely trying to copy a "gradient number" into a vector that only takes Floats or something.

If you initialize your vector with zeros(10) or something, you should instead initialize it like zeros(eltype(x), 10) where x is your input vector.

[Patrick Kofod Mogensen]

I am off to bed, but I think that might do it. Thank you!

[Patrick Kofod Mogensen]

Of course you were right, that I needed to make the eltypes a bit more flexible, so forwarddiff can do its magic.

With that fixed, it's still not working though. It runs for ~10 seconds, and then it just says:

julia> j = ForwardDiff.jacobian(phi, P)
Killed

I tried to search for "Killed" in the github repository but did not find anything. Any idea what the problem is? Am I trying to AD things that cannot be? (I have BLAS calls in my mapping for example)

[Kristoffer Carlsson]

Do you have explicit blas calls? That won't work. Instead you should just use the more high level functions (either * or the Ax_mul_Bx versions) so that Julias multiple dispatch can do its thing and run the generic versions (julia implementations) of the functions when it is called with "graident numbers".

Easiest would be if you could post a small runnable example of something that models what you actually want to do.

[Patrick Kofod Mogensen]

I know help is much easier to provide with a MWE, but it must be the explicit BLAS-calls, I'm sure. I will re-write the function using Julia functions only. If it doesn't work, I'll have to try to reduce my mapping to a simpler version I can show here. Thanks for taking your time so far.

[Jarrett Revels]

The restriction that one needs to use generic functions (and as such, can't use explicit BLAS calls) isn't clearly stated in ForwardDiff.jl's documentation; I just made an issue to track documenting it. Feel free to open up an issue on ForwardDiff.jl's repository if you run into further errors or unexpected performance bugs.

Best,

Jarrett