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

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

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