Finding Jacobian using automatic differentiation for a vector function

General Usage
,
1

Hi, I want to use automatic differentiation to calculate the Jacobian of vector functions. However, I am running into errors that I don’t know how to solve. Here is a minimal working example:

using ForwardDiff
using ReverseDiff

lbs = [1.0, 2.0]
ubs = [2.0, 3.0]

θ = [rand(2); rand(2); 0.1]

n_fluxes = 2

function hfunc(θ)
    h = zeros(n_fluxes * 2 + 1)
    h[end] = θ[end]
    for (i, ii) in zip(1:n_fluxes, n_fluxes .+ (1:n_fluxes))
        if θ[i] < 0.5
            h[i] = -lbs[i]
            h[ii] = ubs[i]
        end
    end
    return h
end

h = hfunc(θ) # works

j = ForwardDiff.jacobian(hfunc, θ)
j = ReverseDiff.jacobian(hfunc, θ)

When I run ForwardDiff on this function, I get the following error:

ERROR: MethodError: no method matching Float64(::ForwardDiff.Dual{ForwardDiff.Tag{typeof(hfunc), Float64}, Float64, 5})
Closest candidates are:
  (::Type{T})(::Real, ::RoundingMode) where T<:AbstractFloat at rounding.jl:200
  (::Type{T})(::T) where T<:Number at boot.jl:760
  (::Type{T})(::SymbolicUtils.Symbolic) where T<:Union{AbstractFloat, Integer, Complex{var"#s134"} where var"#s134"<:Integer, Complex{var"#s135"} where var"#s135"<:AbstractFloat} at C:\Users\St. Elmo\.julia\packages\Symbolics\LfLYY\src\Symbolics.jl:129
  ...
Stacktrace:
 [1] convert(#unused#::Type{Float64}, x::ForwardDiff.Dual{ForwardDiff.Tag{typeof(hfunc), Float64}, Float64, 5})
   @ Base .\number.jl:7
 [2] setindex!(A::Vector{Float64}, x::ForwardDiff.Dual{ForwardDiff.Tag{typeof(hfunc), Float64}, Float64, 5}, i1::Int64)
   @ Base .\array.jl:843
 [3] hfunc(θ::Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(hfunc), Float64}, Float64, 5}})
   @ Main .\REPL[76]:3
 [4] vector_mode_dual_eval!
   @ C:\Users\St. Elmo\.julia\packages\ForwardDiff\CkdHU\src\apiutils.jl:37 [inlined]
 [5] vector_mode_jacobian(f::typeof(hfunc), x::Vector{Float64}, cfg::ForwardDiff.JacobianConfig{ForwardDiff.Tag{typeof(hfunc), Float64}, Float64, 5, Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(hfunc), Float64}, Float64, 5}}})
   @ ForwardDiff C:\Users\St. Elmo\.julia\packages\ForwardDiff\CkdHU\src\jacobian.jl:148
 [6] jacobian(f::Function, x::Vector{Float64}, cfg::ForwardDiff.JacobianConfig{ForwardDiff.Tag{typeof(hfunc), Float64}, Float64, 5, Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(hfunc), Float64}, Float64, 5}}}, ::Val{true})
   @ ForwardDiff C:\Users\St. Elmo\.julia\packages\ForwardDiff\CkdHU\src\jacobian.jl:21
 [7] jacobian(f::Function, x::Vector{Float64}, cfg::ForwardDiff.JacobianConfig{ForwardDiff.Tag{typeof(hfunc), Float64}, Float64, 5, Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(hfunc), Float64}, Float64, 5}}}) (repeats 2 times)
   @ ForwardDiff C:\Users\St. Elmo\.julia\packages\ForwardDiff\CkdHU\src\jacobian.jl:19
 [8] top-level scope
   @ REPL[78]:1

When I run ReverseDiff on this function, I get the following error:

ERROR: ArgumentError: Converting an instance of ReverseDiff.TrackedReal{Float64, Float64, ReverseDiff.TrackedArray{Float64, Float64, 1, Vector{Float64}, Vector{Float64}}} to Float64 is not defined. Please use `ReverseDiff.value` instead.
 [1] convert(#unused#::Type{Float64}, t::ReverseDiff.TrackedReal{Float64, Float64, ReverseDiff.TrackedArray{Float64, Float64, 1, Vector{Float64}, Vector{Float64}}})
   @ ReverseDiff C:\Users\St. Elmo\.julia\packages\ReverseDiff\DkqZV\src\tracked.jl:261
 [2] setindex!(A::Vector{Float64}, x::ReverseDiff.TrackedReal{Float64, Float64, ReverseDiff.TrackedArray{Float64, Float64, 1, Vector{Float64}, Vector{Float64}}}, i1::Int64)
   @ Base .\array.jl:843
 [3] hfunc(θ::ReverseDiff.TrackedArray{Float64, Float64, 1, Vector{Float64}, Vector{Float64}})
   @ Main .\REPL[76]:3
 [4] ReverseDiff.JacobianTape(f::typeof(hfunc), input::Vector{Float64}, cfg::ReverseDiff.JacobianConfig{ReverseDiff.TrackedArray{Float64, Float64, 1, Vector{Float64}, Vector{Float64}}, Nothing})
   @ ReverseDiff C:\Users\St. Elmo\.julia\packages\ReverseDiff\DkqZV\src\api\tape.jl:229
 [5] jacobian(f::Function, input::Vector{Float64}, cfg::ReverseDiff.JacobianConfig{ReverseDiff.TrackedArray{Float64, Float64, 1, Vector{Float64}, Vector{Float64}}, Nothing}) (repeats 2 times)
   @ ReverseDiff C:\Users\St. Elmo\.julia\packages\ReverseDiff\DkqZV\src\api\jacobians.jl:23
 [6] top-level scope
   @ REPL[79]:1

I cannot use Zygote for this because the function mutates arrays internally, and I am not sure how to avoid doing this, since the h matrix basically gets constructed from θ and some external variables that are constants (lbs and ubs here).

Does anybody know how to fix these kinds of issues?

2

For ForwardDiff there may be other things going on as well, but it is important to get things fully generic. Here are a few hints on this 5. Introduction to Types and Generic Programming — Quantitative Economics with Julia

In particular, see the zeros

1 Like
3

The problem is that this always allocates a Float64 array, regardless of the type of θ. ForwardDiff works by using “dual numbers” instead of regular numbers, but such numbers can’t be stored in h.

In general, it is a good practice in Julia to make your functions use the types of their arguments. For example, if θ is single precision (Float64) or arbitrary precision (BigFloat), then you want to use the same type. This strategy will automatically work with dual numbers as well.

Also, it is good practice to avoid global variables! (In Julia, globals kill performance, but in general it leads to difficult-to-maintain programs.)

Something like:

function hfunc(θ, lbs, ubs)
    isodd(length(θ)) || throw(ArgumentError("expecting odd-length array"))
    n_fluxes = length(θ) ÷ 2
    length(lbs) == length(ubs) == n_fluxes || throw(DimensionMismatch())
    h = zero(θ) # array of 0's of the same size and type as θ
    h[end] = θ[end]
    for i = 1:n_fluxes
        if θ[i] < 0.5
            h[i] = -lbs[i]
            h[i + n_fluxes] = ubs[i]
        end
    end
    return h
end

and then do

j = ForwardDiff.jacobian(x -> hfunc(x, lbs, ubs), θ)

to pass the extra arguments through to hfunc.

3 Likes
4

Thank you both very much :smiley: