Jacobian of a (multivariate) function

New to Julia
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 Likes
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
7 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)
5 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
5 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
1 Like
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
1 Like
10

Hi, thanks for your solutions. Related to this question: why would the following break?

using ForwardDiff
import ForwardDiff.jacobian

function f((x))
    F = zeros(3)
    F[1] = x[1]^2+x[3]
    F[2] = x[1]+x[2]
    F[3] = x[2]^2 + x[3]^2
    return F
end

x0 = [1,2,3]
J0 = jacobian(f, x0)

Here is the error message:
ERROR: LoadError: MethodError: no method matching Float64(::ForwardDiff.Dual{ForwardDiff.Tag{typeof(f),Int64},Int64,3})

It is easy to verify the solution is

\begin{bmatrix} 2 & 0 & 1 \\ 1 & 1 & 0 \\ 0 & 4 & 6 \end{bmatrix}

This is given correctly by your method, by instead writing f as f((x)) = [x[1]^2+x[3], x[1]+x[2], x[2]^2 + x[3]^2]

I’m writing it the first way since my real function is much more sophisticated and would be impossible to write as an anonymous function. I need to iterate over a list of things and do some calculations, and at the end of each loop write F[i] = loop i result

Hope you could help with this, thanks!

11

  1. You only need function f(x)

  2. Defining F as zeros(3) restricts its entries to be Float64s, whereas to use ForwardDiff
    they need to be Duals.

Change it to F = similar(x) or F = zero.(x)

to get something of the correct type:

function f(x)
    F = zero.(x)

    F[1] = x[1]^2 + x[3]
    F[2] = x[1] + x[2]
    F[3] = x[2]^2 + x[3]^2

    return F
end
6 Likes