Forward Diff to cast Dual type to Float64

Specific Domains
1

Hello, I’m a complete beginner to Julia. I’m currently using ForwardDiff to compute the Jacobian matrix of a vector-valued function, and that works very well. Now, I want to continue using ForwardDiff to compute ∂JAC/∂τ[j], but it doesn’t work. Here is my code

function f(τ::Vector)
        T = eltype(τ)
        F = Vector{T}()

        alpha = zeros(T, size(τ))
        Xτ = X.(τ, alpha, system.param)
        Yτ = Y.(τ, alpha, system.param)
        RT = T(RT)
 
        eq_cons  = A * Xτ 
        eq_omega = sum(Xτ) 
        eq_equil =  RT .* Yτ
        F = vcat(eq_cons, eq_omega, eq_equil)
        return F
    end
function JAC(τ)
            τ = ForwardDiff.partials.(τ)
            cfg = ForwardDiff.JacobianConfig(f, τ, ForwardDiff.Chunk{length(τ)}())
            out = DiffResults.JacobianResult(τ)
            ForwardDiff.jacobian!(out, f, τ, cfg)
            return DiffResults.jacobian(out)
end
function dJAC_dτj(τ, j::Int)
            τ = ForwardDiff.partials.(τ)
            g(x) = vec(JAC(x))                 
            G = ForwardDiff.jacobian(g, τ)     
            result = reshape(dJAC_dτj(τ,j), n, n)     # nxn is size of matrix JAC(τ). This line meets error.
            return result
end

The error I get is:

no method matching Float64(::ForwardDiff.Dual{ForwardDiff.Tag{Main.name.var"#f#13"{Speciation, Vector{Float64}, Vector{Float64}, Matrix{Real}, Matrix{Real}}, Float64}, Float64, 4})

It seems that Julia is trying to cast a Dual type to Float64. The type Float64 exists, but there’s no method for converting this kind of Dual number to Float64. However, when I try adding this line before the problematic part:

τ = ForwardDiff.value.(τ)

the code runs without error, but ∂JAC/∂τ[j] ends up being all zeros. After checking, I realized that this happens because calling ForwardDiff.value.(τ) erases the derivative information, so the result is computed as if everything is constant.

My question: How can I compute ∂JAC/∂τ[j] given that I already have JAC and the function f as described above? Thank you very much for your time and concerns.

2

Hi @Kate, welcome to the community!

When you post, try to make sure your code examples are self-contained, this will be useful for others to lend a hand. Here we’re missing some imports and definitions, so I used a simpler function.
What you are after is a directional derivative but ForwardDiff.jl doesn’t have a direct API for that. Instead, you can interpret \partial J / \partial \tau_j as the derivative of \alpha \mapsto J(\tau + \alpha e_j) at \alpha = 0, where e_j is a Euclidean basis vector.
You could also compute the full Jacobian of J but this would be wasteful, requiring n times more computation than necessary to get just the j-th partial derivative.

using ForwardDiff

function basis(n, j)
    e = zeros(Bool, n)
    e[j] = 1
    return e
end

f(τ::Vector) = sin.(τ) + cos.(reverse(τ))  # replace with your function
∂f(τ::Vector) = ForwardDiff.jacobian(f, τ)  # returns a matrix

function ∂²f(τ::Vector, j::Int)
    ej = basis(length(τ), j)
    return ForwardDiff.derivative(α -> ∂f(τ + α * ej), 0)  # returns a matrix
end

julia> ∂²f([3.0, 4.0], 1)
2×2 Matrix{Float64}:
 -0.14112   0.0
  0.989992  0.0
1 Like
3

Thank you very much for your time and support. I will be more careful in the future. Wish you a good day.

1 Like