# 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 + r)
end
``````

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 + r);

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

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

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}:
0.0
2.0
julia> f(r1)
0.9092974268256817
julia> f.(r1)
ERROR: BoundsError
Stacktrace:
 getindex at ./number.jl:38 [inlined]
 (::##1#2)(::Float64) at ./REPL:1
 macro expansion at ./broadcast.jl:153 [inlined]
 macro expansion at ./simdloop.jl:73 [inlined]
 macro expansion at ./broadcast.jl:147 [inlined]
 _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
 broadcast_t(::Function, ::Type{T} where T, ::Tuple{Base.OneTo{Int64}}, ::CartesianRange{CartesianIndex{1}}, ::Array{Float64,1}) at ./broadcast.jl:268
``````

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
0.0:0.25:1.0

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}:
0.0:0.25:1.0
0.0:0.25:1.0
``````

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).

Does

``````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}:
0.0
0.479426
0.841471
0.997495
0.909297
``````

give you what you want?

Perfect, thanks.

1 Like

You could also reinterpret a 2d array as a vector of static arrays.
http://juliaarrays.github.io/StaticArrays.jl/latest/pages/api.html#Arrays-of-static-arrays-1