Jacobian calculation with FastDifferentiation.jl

Hi there,

I wanted to test this package to compute the Jacobian of a relatively lengthy function.

However, I encounter an issue because the result is apparently not correct if I compare with other methods. Yet, the code works on “simpler” functions.

Am I missing something?

Here is my code:

import FastDifferentiation as fad
import ForwardDiff as fod
import Symbolics as sym

fad.@variables x1 x2 x3 x4 u1 u2 u3 u4 u5 u6

### Quite simple function
#f(x1,x2,x3,x4,u1,u2,u3,u4,u5,u6) = x1*x2/cos(x3)+sin(x4)*tan(u1)*u2^2

### More complex function
f(x1,x2,x3,x4,u1,u2,u3,u4,u5,u6) = (((0.0588 / (0.010425239999999999 - (0.008230118399999998 * (cos(x2) ^ 2)))) * ((((((-0.04 * (((810.8 + (1621.6 * u1)) * (((-0.034910000000000004 * u1) - 0.09076000000000001) - (-0.04 * x1))) - ((54.1 + (108.1 * u1)) * (-0.04 * x3)))) + (-0.04 * (((810.8 + (1621.6 * u2)) * (((0.034910000000000004 * u2) - -0.027930000000000003) - (-0.04 * x1))) - ((54.1 + (108.1 * u2)) * (-0.04 * x3))))) + (-0.027999999999999997 * (((810.8 + (1621.6 * u5)) * (((-0.05498 * u5) - 0.14294) - ((-0.027999999999999997 * x1) + (-0.035 * x2)))) - ((54.1 + (108.1 * u5)) * ((-0.027999999999999997 * x3) + (-0.035 * x4)))))) + (-0.027999999999999997 * (((810.8 + (1621.6 * u6)) * (((0.05498 * u6) - -0.01343) - ((-0.027999999999999997 * x1) + (-0.035 * x2)))) - ((54.1 + (108.1 * u6)) * ((-0.027999999999999997 * x3) + (-0.035 * x4)))))) + (((-0.35 * sin((x1 + x2))) - (0.3 * sin(x1))) * ((((0.3 * x3) * cos(x1)) + ((0.35 * (x3 + x4)) * cos((x1 + x2)))) * (1.5 * 0)))) - (((-0.18144 * (x4 * sin(x2))) * x3) + (((0.09072 * sin(x2)) * -(x4)) * x4)))) + (((-0.0588 - (0.09072 * cos(x2))) / (0.010425239999999999 - (0.008230118399999998 * (cos(x2) ^ 2)))) * ((((((-0.025 * (((810.8 + (1621.6 * u3)) * (((-0.02182 * u3) - 0.05672) - (-0.025 * x2))) - ((54.1 + (108.1 * u3)) * (-0.025 * x4)))) + (-0.025 * (((810.8 + (1621.6 * u4)) * (((0.02182 * u4) - 0.00436) - (-0.025 * x2))) - ((54.1 + (108.1 * u4)) * (-0.025 * x4))))) + (-0.035 * (((810.8 + (1621.6 * u5)) * (((-0.05498 * u5) - 0.14294) - ((-0.027999999999999997 * x1) + (-0.035 * x2)))) - ((54.1 + (108.1 * u5)) * ((-0.027999999999999997 * x3) + (-0.035 * x4)))))) + (-0.035 * (((810.8 + (1621.6 * u6)) * (((0.05498 * u6) - -0.01343) - ((-0.027999999999999997 * x1) + (-0.035 * x2)))) - ((54.1 + (108.1 * u6)) * ((-0.027999999999999997 * x3) + (-0.035 * x4)))))) + ((-0.35 * sin((x1 + x2))) * ((((0.3 * x3) * cos(x1)) + ((0.35 * (x3 + x4)) * cos((x1 + x2)))) * (1.5 * 0)))) - (((0.09072 * sin(x2)) * x3) * x3))))

f_expr = f(x1,x2,x3,x4,u1,u2,u3,u4,u5,u6)

x = [x1,x2,x3,x4]
u = [u1,u2,u3,u4,u5,u6]

∇f = fad.jacobian([f_expr],x)

∇f_exe! = fad.make_function(∇f, x, u, in_place=true)

x0 = ones(4)
u0 = ones(6)

∇f0 = zeros(1,4)
∇f_exe!(∇f0, vcat(x0,u0))

println("with FastDifferentiation:")
@show ∇f0

∇f20 = fod.jacobian( x -> [f(x...,u0...)], x0)

println("with ForwardDiff:")
@show ∇f20

sym.@variables x1 x2 x3 x4 u1 u2 u3 u4 u5 u6
f_expr = f(x1,x2,x3,x4,u1,u2,u3,u4,u5,u6)
x = [x1,x2,x3,x4]
u = [u1,u2,u3,u4,u5,u6]

∇f = sym.jacobian([f_expr],x)

∇f_exe = sym.build_function(∇f, x, u, expression=Val{false}, nanmath=false)[1]

∇f30 = ∇f_exe(x0,u0)

println("with Symbolics:")
@show ∇f30

nothing

When I execute this, I get:

with FastDifferentiation:
∇f0 = [-66.94400802991359 62.90102233528086 -0.36119179526925516 7.973221105027532]
with ForwardDiff:
∇f20 = [-20.931923273862566 102.94501544413166 1.775014645202905 7.973221105027534]
with Symbolics:
∇f30 = [-20.93192327386258 102.94501544413171 1.7750146452029059 7.973221105027532]

using:
FastDifferentiation v0.3.17
ForwardDiff v0.10.36
Symbolics v6.3.0

Thanks for any feedback!

Can you try it with DifferentiationInterface.jl?

Sure. I simplified a bit (this is likely not a MWE but it’s simple enough I believe).

using DifferentiationInterface
import ForwardDiff, FastDifferentiation, Symbolics

f(x) = [(((0.0588 / (0.010425239999999999 - (0.008230118399999998 * (cos(x[2]) ^ 2)))) * ((((((-0.04 * (((810.8 + (1621.6 * 1)) * (((-0.034910000000000004 * 1) - 0.09076000000000001) - (-0.04 * x[1]))) - ((54.1 + (108.1 * 1)) * (-0.04 * x[3])))) + (-0.04 * (((810.8 + (1621.6 * 1)) * (((0.034910000000000004 * 1) - -0.027930000000000003) - (-0.04 * x[1]))) - ((54.1 + (108.1 * 1)) * (-0.04 * x[3]))))) + (-0.027999999999999997 * (((810.8 + (1621.6 * 1)) * (((-0.05498 * 1) - 0.14294) - ((-0.027999999999999997 * x[1]) + (-0.035 * x[2])))) - ((54.1 + (108.1 * 1)) * ((-0.027999999999999997 * x[3]) + (-0.035 * x[4])))))) + (-0.027999999999999997 * (((810.8 + (1621.6 * 1)) * (((0.05498 * 1) - -0.01343) - ((-0.027999999999999997 * x[1]) + (-0.035 * x[2])))) - ((54.1 + (108.1 * 1)) * ((-0.027999999999999997 * x[3]) + (-0.035 * x[4])))))) + (((-0.35 * sin((x[1] + x[2]))) - (0.3 * sin(x[1]))) * ((((0.3 * x[3]) * cos(x[1])) + ((0.35 * (x[3] + x[4])) * cos((x[1] + x[2])))) * (1.5 * 0)))) - (((-0.18144 * (x[4] * sin(x[2]))) * x[3]) + (((0.09072 * sin(x[2])) * -(x[4])) * x[4])))) + (((-0.0588 - (0.09072 * cos(x[2]))) / (0.010425239999999999 - (0.008230118399999998 * (cos(x[2]) ^ 2)))) * ((((((-0.025 * (((810.8 + (1621.6 * 1)) * (((-0.02182 * 1) - 0.05672) - (-0.025 * x[2]))) - ((54.1 + (108.1 * 1)) * (-0.025 * x[4])))) + (-0.025 * (((810.8 + (1621.6 * 1)) * (((0.02182 * 1) - 0.00436) - (-0.025 * x[2]))) - ((54.1 + (108.1 * 1)) * (-0.025 * x[4]))))) + (-0.035 * (((810.8 + (1621.6 * 1)) * (((-0.05498 * 1) - 0.14294) - ((-0.027999999999999997 * x[1]) + (-0.035 * x[2])))) - ((54.1 + (108.1 * 1)) * ((-0.027999999999999997 * x[3]) + (-0.035 * x[4])))))) + (-0.035 * (((810.8 + (1621.6 * 1)) * (((0.05498 * 1) - -0.01343) - ((-0.027999999999999997 * x[1]) + (-0.035 * x[2])))) - ((54.1 + (108.1 * 1)) * ((-0.027999999999999997 * x[3]) + (-0.035 * x[4])))))) + ((-0.35 * sin((x[1] + x[2]))) * ((((0.3 * x[3]) * cos(x[1])) + ((0.35 * (x[3] + x[4])) * cos((x[1] + x[2])))) * (1.5 * 0)))) - (((0.09072 * sin(x[2])) * x[3]) * x[3]))))]

x = [1.0, 2.0, 3.0, 4.0]

backends = [AutoSymbolics(), AutoForwardDiff(), AutoFastDifferentiation()]

for backend in backends
    println("With " * string(backend))
    ∇fx = jacobian(f, backend, x)
    @show ∇fx
end

nothing

Result is:

With AutoSymbolics()
∇fx = [-64.62263875764114 -141.61479381476425 1.1598116002053955 6.871693628806222]
With AutoForwardDiff()
∇fx = [-64.62263875764113 -141.61479381476425 1.1598116002053955 6.871693628806222]
With AutoFastDifferentiation()
∇fx = [-53.473435731408145 -127.678290031973 1.1598116002053955 6.871693628806222]

Well, it’s perhaps related to this known issue:
https://github.com/brianguenter/FastDifferentiation.jl/issues/65