How to use AD for function with multiple inputs?

Specific Domains Optimization (Mathematical)
,
1

Over the past several days I have been struggling with AD for a particular function. I have boiled down the problem to this:

using ForwardDiff

function F!(resx, resy, x, y)
    resx[1] = 1 - x[1]
    resx[2] = 10(x[2]-y[1]^2)
    resy[1] = 1 - y[1]
    resy[2] = 10(y[2]-x[1]^2)
    return resx, resy
end

resx = zeros(5)
resy = zeros(5)
x = ones(5)
y = ones(5)

len = length(x)
xcache = zeros(eltype(x), len)
Fcache = zeros(eltype(x), len)
Jcache = zeros(eltype(x), len, len)

function f!(resx, x)
   F!(resx, resy, x, y) 
end

ForwardDiff.jacobian!(Jcache, f!, Fcache, x)

that last line results in this error:

ERROR: LoadError: MethodError: no method matching Float64(::ForwardDiff.Dual{ForwardDiff.Tag{typeof(f!), 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})(::AbstractChar) where T<:Union{AbstractChar, Number} at char.jl:50
  ...
Stacktrace:
 [1] convert(#unused#::Type{Float64}, x::ForwardDiff.Dual{ForwardDiff.Tag{typeof(f!), Float64}, Float64, 5})
   @ Base ./number.jl:7
 [2] setindex!(A::Vector{Float64}, x::ForwardDiff.Dual{ForwardDiff.Tag{typeof(f!), Float64}, Float64, 5}, i1::Int64)
   @ Base ./array.jl:839
 [3] F!(resx::Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(f!), Float64}, Float64, 5}}, resy::Vector{Float64}, x::Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(f!), Float64}, Float64, 5}}, y::Vector{Float64})
   @ Main ~/discourseQ5.jl:9
 [4] f!(resx::Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(f!), Float64}, Float64, 5}}, x::Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(f!), Float64}, Float64, 5}})
   @ Main ~/discourseQ5.jl:24
 [5] vector_mode_dual_eval
   @ ~/.julia/packages/ForwardDiff/QOqCN/src/apiutils.jl:44 [inlined]
 [6] vector_mode_jacobian!(result::Matrix{Float64}, f!::typeof(f!), y::Vector{Float64}, x::Vector{Float64}, cfg::ForwardDiff.JacobianConfig{ForwardDiff.Tag{typeof(f!), Float64}, Float64, 5, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(f!), Float64}, Float64, 5}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(f!), Float64}, Float64, 5}}}})
   @ ForwardDiff ~/.julia/packages/ForwardDiff/QOqCN/src/jacobian.jl:172
 [7] jacobian!(result::Matrix{Float64}, f!::Function, y::Vector{Float64}, x::Vector{Float64}, cfg::ForwardDiff.JacobianConfig{ForwardDiff.Tag{typeof(f!), Float64}, Float64, 5, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(f!), Float64}, Float64, 5}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(f!), Float64}, Float64, 5}}}}, ::Val{true})
   @ ForwardDiff ~/.julia/packages/ForwardDiff/QOqCN/src/jacobian.jl:78
 [8] jacobian!(result::Matrix{Float64}, f!::Function, y::Vector{Float64}, x::Vector{Float64}, cfg::ForwardDiff.JacobianConfig{ForwardDiff.Tag{typeof(f!), Float64}, Float64, 5, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(f!), Float64}, Float64, 5}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(f!), Float64}, Float64, 5}}}}) (repeats 2 times)
   @ ForwardDiff ~/.julia/packages/ForwardDiff/QOqCN/src/jacobian.jl:76
 [9] top-level scope
   @ ~/discourseQ5.jl:27
in expression starting at /Users/amrods/discourseQ5.jl:27

I think that happens because I am “mixing” tracked and non-tracked inputs in resx[2] and resy[2].
Is there a way to make AD work for such a function?

2

Your issue is that resy is a Vector{Float64}, but you are trying to store a ForwardDiff.Dual in it.

You can see this if you print the types of the arguments to F!.

function F!(resx, resy, x, y)
    @show typeof(resx)
    @show typeof(resy)
    @show typeof(x)
    @show typeof(y)
    resx[1] = 1 - x[1]
    resx[2] = 10(x[2]-y[1]^2)
    resy[1] = 1 - y[1]
    resy[2] = 10(y[2]-x[1]^2)
    return resx, resy
end
2 Likes
3

I see, jacobian! converts resx but not resy.

4

You can make caches of the appropriate type using similar, based on the input. Or you can let the function make them, but it’s not really converting the global resx, since the cache needed is larger. I think that if you use JacobianConfig these can be saved & re-used, but they won’t be as written here:

function fx1(x)
    resx = similar(x,5) .= 0
    @show eltype(resx)  # Dual{Tag{...}, ...}
    resy = similar(x,5) .= 0
    F!(resx, resy, x, y)[1]
end

ForwardDiff.jacobian(fx1, x)
ForwardDiff.jacobian!(Jcache, fx1, x)

function fx1!(resx, x)
    @show pointer(resx) # not the same as global resx
    resy = similar(x,5) .= 0
    F!(resx, resy, x, y)[1]
end

pointer(resx)
ForwardDiff.jacobian(fx1!, resx, x)
ForwardDiff.jacobian!(Jcache, fx1!, resx, x)
1 Like
5

Isn’t that pointer a function for the C interface? I find it strange that it works here.