Jacobian of a (multivariate) function

#1

I am trying to replicate this using the ForwardDiff module

33%20PM

But this is what I got

julia> using ForwardDiff

julia> f(x,y)=[x^2+y^3-1,x^4 - y^4 + x*y]
f (generic function with 1 method)

julia> a = [1.0,1.0]
2-element Array{Float64,1}:
 1.0
 1.0

julia> ForwardDiff.jacobian(f, a)
ERROR: MethodError: no method matching f(::Array{ForwardDiff.Dual{ForwardDiff.Tag{typeof(f),Float64},Float64,2},1})
Closest candidates are:
  f(::Any, ::Any) at REPL[2]:1
#2

ForwardDiff.jacobian only works with functions that accept an array as an argument, this is pretty easy to fix:

julia> using ForwardDiff

julia> f(x,y)=[x^2+y^3-1,x^4 - y^4 + x*y]
f (generic function with 1 method)

julia> a = [1.0,1.0]
2-element Array{Float64,1}:
 1.0
 1.0

julia> ForwardDiff.jacobian(x ->f(x[1],x[2]), a)
2×2 Array{Float64,2}:
 2.0   3.0
 5.0  -3.0
5 Likes
#3

Note by the way that you can make your code ~7x faster (on my machine) using StaticArrays:

using ForwardDiff
using StaticArrays

f(x,y) = @SVector [x^2+y^3-1,x^4 - y^4 + x*y]
a = @SVector [1.0,1.0]
ForwardDiff.jacobian(x ->f(x[1],x[2]), a)
4 Likes
#4

Apparently this works as well, just use the THREE DOTS

ForwardDiff.jacobian(x -> f(x...), a)

God bless those who invented the THREE DOTS

1 Like
#5

God bless those who invented destructuring

julia> using ForwardDiff
[ Info: Recompiling stale cache file /home/pkofod/.julia/compiled/v1.0/ForwardDiff/k0ETY.ji for ForwardDiff [f6369f11-7733-5829-9624-2563aa707210]

julia> f((x,y))=[x^2+y^3-1,x^4 - y^4 + x*y]
f (generic function with 1 method)

julia> a = [1.0,1.0]
2-element Array{Float64,1}:
 1.0
 1.0

julia> ForwardDiff.jacobian(x ->f(x), a)
2×2 Array{Float64,2}:
 2.0   3.0
 5.0  -3.0
4 Likes
#6

Why not just use this? “ForwardDiff.jacobian(f,a)”

julia> using ForwardDiff

julia> f((x,y)) = [x^2+y^3-1,x^4 - y^4 + x*y]
f (generic function with 1 method)

julia> a = [1.0,1.0]
2-element Array{Float64,1}:
 1.0
 1.0

julia> ForwardDiff.jacobian(f,a)
2×2 Array{Float64,2}:
 2.0   3.0
 5.0  -3.0
#7

That was a mistake, yes!

#8

Related question. What am I doing wrong here? It seems like the docs and everything says it should support SVector but it isn’t managing the conversion.

function BicycleDynamics(plant::Bicycle, x::SVector{5, T}, u::SVector{2, T}) where {T <: Real}
    # of the traction wheel
    turning_radius = plant.wheelbase / tan(x[4])
    translational_velocity = x[5]
    angular_velocity = x[5] / turning_radius
    return SVector(cos(x[3]), sin(x[3]), angular_velocity, u[1], u[2])
end
bike = MakeBike()
ad_adapter(x::SVector{7, T}) where {T <: Real} = BicycleDynamics(bike, x[1:5], x[6:7])
ad_point = SVector{7, Float64}(1.0, 1.0, 0.0, 1.0, 1.0, 1.0, 1.0)
ForwardDiff.jacobian(ad_adapter, ad_point)

Produces the error:

ERROR: LoadError: MethodError: no method matching BicycleDynamics(::Bicycle, ::Array{ForwardDiff.Dual{ForwardDiff.Tag{typeof(ad_adapter),Float64},Float64,7},1}, ::Array{ForwardDiff.Dual{ForwardDiff.Tag{typeof(ad_adapter),Float64},Float64,7},1})

EDIT:
Figured it out. Problem was slicing an SVector with constant indices doesn’t result in another SVector. So I need to convert manually

#9

To expand on this, if what someone wants a callable function, j, that computes the Jacobian then this should work:

julia> f(x,y)=[x^2+y^3-1,x^4-y^4+x*y]
julia> j(x)=ForwardDiff.jacobian(x->f(x[1],x[2]), x)
julia> j([0.5;0.5])
2x2 Array{Float64,2}:
 1.0  0.75
 1.0  0.0