Vectorization of ForwardDiff.gradient function

New to Julia
1

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])
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?

2

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}:
 0.783327
 0.29552 

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

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:
 [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
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?

4

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.

5

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

6

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?

7

Perfect, thanks.

1 Like
8

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