Vectorization of ForwardDiff.gradient function


I’ve got a scalar valued function of two variables, and I would like to be able to vectorize the evaluation of both the function and its gradient, as calculated using ForwardDiff. Writing the function as:

function f(r)
return sin(r[1] + r[2])

I can calculate the gradient as

gradf = (r)->ForwardDiff.gradient(f,r)

and these work fine for single input vectors. However, while f.(rvals) works fine on an array of inputs, gradf.(rvals) fails. What is the recommended way to accomplish this?


Can you post some stacktrace? What is not working? It works fine for me:

julia> using ForwardDiff

julia> f(r) = sin(r[1] + r[2]);

julia> gradf = (r)->ForwardDiff.gradient(f,r);

julia> rvals = [[0.4, 0.5], [0.2, 0.1]];

julia> f.(rvals)
2-element Array{Float64,1}:

julia> gradf.(rvals)
2-element Array{Array{Float64,1},1}:
 [0.62161, 0.62161]  
 [0.955336, 0.955336]


So strange things happen for me with this. First, suppose I apply the vectorized funciton, but on a single point:

julia> r1 = [0., 2.]
2-element Array{Float64,1}:
julia> f(r1)
julia> f.(r1)
ERROR: BoundsError
 [1] getindex at ./number.jl:38 [inlined]
 [2] (::##1#2)(::Float64) at ./REPL[2]:1
 [3] macro expansion at ./broadcast.jl:153 [inlined]
 [4] macro expansion at ./simdloop.jl:73 [inlined]
 [5] macro expansion at ./broadcast.jl:147 [inlined]
 [6] _broadcast!(::##1#2, ::Array{Float64,1}, ::Tuple{Tuple{Bool}}, ::Tuple{Tuple{Int64}}, ::Array{Float64,1}, ::Tuple{}, ::Type{Val{0}}, ::CartesianRange{CartesianIndex{1}}) at ./broadcast.jl:139
 [7] broadcast_t(::Function, ::Type{T} where T, ::Tuple{Base.OneTo{Int64}}, ::CartesianRange{CartesianIndex{1}}, ::Array{Float64,1}) at ./broadcast.jl:268
 [8] broadcast_c at ./broadcast.jl:314 [inlined]
 [9] broadcast(::Function, ::Array{Float64,1}) at ./broadcast.jl:434

But I think I see what my issue was. I had been generating an Nx2 Array{Float64,2} for my input as

julia> x = linspace(0,1,5);

julia> y=x

julia> r = hcat(x,y)
5×2 Array{Float64,2}:
 0.0   0.0 
 0.25  0.25
 0.5   0.5 
 0.75  0.75
 1.0   1.0 

instead of an N array of points. How do I convert between the data types?


First, suppose I apply the vectorized funciton, but on a single point

Basically, just don’t do this. f.(r1) syntactically means "apply f elementwise to r1. If r1 is a vector, then you’re asking Julia to apply f to each element of the vector, which isn’t what you actually want. So use f(r1) when you have one element or f.(rvals) when you have a vector of elements.

How do I convert between the data types?

You can just construct it directly as [x, y], or if you already have the N array of points, you can do [r[:, i] for i in 1:size(r, 2)], for example.


Ok, one last question. [x,y] constructs

julia> [x,y]
2-element Array{StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}},1}:

so I end up with an array of length 2 with each element of length 5. Is there a simple way to get the opposite (which is what I want).



julia> a = [[x[i], y[i]] for i in 1:length(x)]
5-element Array{Array{Float64,1},1}:
 [0.0, 0.0]  
 [0.25, 0.25]
 [0.5, 0.5]  
 [0.75, 0.75]
 [1.0, 1.0]  

julia> f.(a)
5-element Array{Float64,1}:

give you what you want?


Perfect, thanks.


You could also reinterpret a 2d array as a vector of static arrays.