Finding Jacobian using automatic differentiation for a vector function

General Usage
,
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