Jacobians

The supertype Jacobian comprises different ways of taking Jacobians:

• JacobianFunction,
• JacobianAutodiff,
• JacobianFiniteDifferences.

We first start by showing JacobianAutodiff:

# the input and output dimensions of this function are the same
F(y::AbstractArray, x::AbstractArray) = y .= tanh.(x)
dim = 3
x = rand(dim)
jac = JacobianAutodiff{eltype(x)}(F, dim)

JacobianAutodiff{Float64, SimpleSolvers.var"#F!#32"{typeof(Main.F)}, ForwardDiff.JacobianConfig{ForwardDiff.Tag{SimpleSolvers.var"#F!#32"{typeof(Main.F)}, Float64}, Float64, 3, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{SimpleSolvers.var"#F!#32"{typeof(Main.F)}, Float64}, Float64, 3}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{SimpleSolvers.var"#F!#32"{typeof(Main.F)}, Float64}, Float64, 3}}}}, Vector{Float64}}(SimpleSolvers.var"#F!#32"{typeof(Main.F)}(Main.F), ForwardDiff.JacobianConfig{ForwardDiff.Tag{SimpleSolvers.var"#F!#32"{typeof(Main.F)}, Float64}, Float64, 3, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{SimpleSolvers.var"#F!#32"{typeof(Main.F)}, Float64}, Float64, 3}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{SimpleSolvers.var"#F!#32"{typeof(Main.F)}, Float64}, Float64, 3}}}}((Partials(1.0, 0.0, 0.0), Partials(0.0, 1.0, 0.0), Partials(0.0, 0.0, 1.0)), (ForwardDiff.Dual{ForwardDiff.Tag{SimpleSolvers.var"#F!#32"{typeof(Main.F)}, Float64}, Float64, 3}[Dual{ForwardDiff.Tag{SimpleSolvers.var"#F!#32"{typeof(Main.F)}, Float64}}(6.95175326059346e-310,6.9517532514738e-310,6.9517532514718e-310,0.0), Dual{ForwardDiff.Tag{SimpleSolvers.var"#F!#32"{typeof(Main.F)}, Float64}}(0.0,0.0,0.0,0.0), Dual{ForwardDiff.Tag{SimpleSolvers.var"#F!#32"{typeof(Main.F)}, Float64}}(0.0,0.0,0.0,0.0)], ForwardDiff.Dual{ForwardDiff.Tag{SimpleSolvers.var"#F!#32"{typeof(Main.F)}, Float64}, Float64, 3}[Dual{ForwardDiff.Tag{SimpleSolvers.var"#F!#32"{typeof(Main.F)}, Float64}}(6.95171387608845e-310,6.95175211301296e-310,6.9517457437574e-310,0.0), Dual{ForwardDiff.Tag{SimpleSolvers.var"#F!#32"{typeof(Main.F)}, Float64}}(0.0,0.0,0.0,0.0), Dual{ForwardDiff.Tag{SimpleSolvers.var"#F!#32"{typeof(Main.F)}, Float64}}(0.0,0.0,0.0,0.0)])), [0.0, 0.0, 0.0])

Instead of calling JacobianAutodiff(f, x) we can equivalently do:

jac = Jacobian{eltype(x)}(F, dim; mode = :autodiff)

JacobianAutodiff{Float64, SimpleSolvers.var"#F!#32"{typeof(Main.F)}, ForwardDiff.JacobianConfig{ForwardDiff.Tag{SimpleSolvers.var"#F!#32"{typeof(Main.F)}, Float64}, Float64, 3, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{SimpleSolvers.var"#F!#32"{typeof(Main.F)}, Float64}, Float64, 3}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{SimpleSolvers.var"#F!#32"{typeof(Main.F)}, Float64}, Float64, 3}}}}, Vector{Float64}}(SimpleSolvers.var"#F!#32"{typeof(Main.F)}(Main.F), ForwardDiff.JacobianConfig{ForwardDiff.Tag{SimpleSolvers.var"#F!#32"{typeof(Main.F)}, Float64}, Float64, 3, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{SimpleSolvers.var"#F!#32"{typeof(Main.F)}, Float64}, Float64, 3}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{SimpleSolvers.var"#F!#32"{typeof(Main.F)}, Float64}, Float64, 3}}}}((Partials(1.0, 0.0, 0.0), Partials(0.0, 1.0, 0.0), Partials(0.0, 0.0, 1.0)), (ForwardDiff.Dual{ForwardDiff.Tag{SimpleSolvers.var"#F!#32"{typeof(Main.F)}, Float64}, Float64, 3}[Dual{ForwardDiff.Tag{SimpleSolvers.var"#F!#32"{typeof(Main.F)}, Float64}}(5.0e-324,6.95171515951764e-310,5.0e-324,6.9517151595192e-310), Dual{ForwardDiff.Tag{SimpleSolvers.var"#F!#32"{typeof(Main.F)}, Float64}}(5.0e-324,6.9517151595208e-310,5.0e-324,6.9517151595224e-310), Dual{ForwardDiff.Tag{SimpleSolvers.var"#F!#32"{typeof(Main.F)}, Float64}}(5.0e-324,6.95171515952396e-310,5.0e-324,6.95171515952554e-310)], ForwardDiff.Dual{ForwardDiff.Tag{SimpleSolvers.var"#F!#32"{typeof(Main.F)}, Float64}, Float64, 3}[Dual{ForwardDiff.Tag{SimpleSolvers.var"#F!#32"{typeof(Main.F)}, Float64}}(5.0e-324,6.951715157124e-310,5.0e-324,6.95171515712557e-310), Dual{ForwardDiff.Tag{SimpleSolvers.var"#F!#32"{typeof(Main.F)}, Float64}}(5.0e-324,6.95171515712715e-310,5.0e-324,6.95171515712873e-310), Dual{ForwardDiff.Tag{SimpleSolvers.var"#F!#32"{typeof(Main.F)}, Float64}}(5.0e-324,6.9517151571303e-310,5.0e-324,6.9517151571319e-310)])), [0.0, 0.0, 0.0])

When calling an instance of Jacobian we can use the function [compute_jacobian!]:

j = zeros(dim, dim)
compute_jacobian!(j, x, jac)

3×3 Matrix{Float64}:
 0.770907  0.0       0.0
 0.0       0.721643  0.0
 0.0       0.0       0.493302

This is equivalent to calling:

jac(j, x)

3×3 Matrix{Float64}:
 0.770907  0.0       0.0
 0.0       0.721643  0.0
 0.0       0.0       0.493302