ForwardDiff & mul!()

Is it possible to make ForwardDiff works with mul!() inside function?

Or how to optimize allocations when I try to get gradient/hessian of likelihood function?

i can find mul! in base, where it is?

this is in LinearAlgebra

LinearAlgebra.mul!

You can allocate the result with element type for the appropriate Dual, but I don’t think it is worth the hassle. Just use non-modifying versions of everything.

I function i have part:

function reml(yv, Zv, p, Xv, θvec, β)
...
    for i = 1:n
        θ2m += Xv[i]'*iV*Xv[i]
        r    = yv[i]-Xv[i]*β
        θ3  += r'*iV*r
    end
 ...
end

where Xv - array of matreces, i try to make temp array before cycle and modyfy like this:

     temp = zeros(p,p)
     for i = 1:n
        mul!(temp, Xv[i]', iV)
        θ2m += temp*Xv[i]
        r    = yv[i]-Xv[i]*β
        θ3  += r'*iV*r
    end

and in general case it reduse allocations, bur ForwardDiff get me error:

ERROR: MethodError: no method matching Float64(::ForwardDiff.Dual{ForwardDiff.Tag{typeof(remlf2),Float64},ForwardDiff.Dual{ForwardDiff.Tag{typeof(remlf2),Float64},Float64,5},5})
Closest candidates are:
  Float64(::Real, ::RoundingMode) where T<:AbstractFloat at rounding.jl:194
  Float64(::T<:Number) where T<:Number at boot.jl:741
  Float64(::Int8) at float.jl:60
  ...
Stacktrace:
 [1] convert(::Type{Float64}, ::ForwardDiff.Dual{ForwardDiff.Tag{typeof(remlf2),Float64},ForwardDiff.Dual{ForwardDiff.Tag{typeof(remlf2),Float64},Float64,5},5}) at .\number.jl:7
 [2] setindex! at .\array.jl:769 [inlined]
 [3] _generic_matmatmul!(::Array{Float64,2}, ::Char, ::Char, ::Array{Float64,2}, ::Array{ForwardDiff.Dual{ForwardDiff.Tag{typeof(remlf2),Float64},ForwardDiff.Dual{ForwardDiff.Tag{typeof(remlf2),Float64},Float64,5},5},2}) at C:\cygwin\home\Administrator\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.1\LinearAlgebra\src\matmul.jl:727
 [4] generic_matmatmul! at C:\cygwin\home\Administrator\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.1\LinearAlgebra\src\matmul.jl:581 [inlined]
 [5] mul!(::Array{Float64,2}, ::Adjoint{Float64,Array{Float64,2}}, ::Array{ForwardDiff.Dual{ForwardDiff.Tag{typeof(remlf2),Float64},ForwardDiff.Dual{ForwardDiff.Tag{typeof(remlf2),Float64},Float64,5},5},2}) at C:\cygwin\home\Administrator\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.1\LinearAlgebra\src\matmul.jl:293
 [6] reml2(::Array{Array{Float64,1},1}, ::Array{Array{Float64,2},1}, ::Int64, ::Array{Array{Float64,2},1}, ::Array{ForwardDiff.Dual{ForwardDiff.Tag{typeof(remlf2),Float64},ForwardDiff.Dual{ForwardDiff.Tag{typeof(remlf2),Float64},Float64,5},5},1}, ::Array{Float64,1}) at G:\Cloud\Julia\Mixed\mixbe.jl:376
 [7] remlf2(::Array{ForwardDiff.Dual{ForwardDiff.Tag{typeof(remlf2),Float64},ForwardDiff.Dual{ForwardDiff.Tag{typeof(remlf2),Float64},Float64,5},5},1}) at G:\Cloud\Julia\Mixed\mixbe.jl:386
 [8] vector_mode_gradient(::typeof(remlf2), ::Array{ForwardDiff.Dual{ForwardDiff.Tag{typeof(remlf2),Float64},Float64,5},1}, ::ForwardDiff.GradientConfig{ForwardDiff.Tag{typeof(remlf2),Float64},ForwardDiff.Dual{ForwardDiff.Tag{typeof(remlf2),Float64},Float64,5},5,Array{ForwardDiff.Dual{ForwardDiff.Tag{typeof(remlf2),Float64},ForwardDiff.Dual{ForwardDiff.Tag{typeof(remlf2),Float64},Float64,5},5},1}}) at C:\...\.julia\packages\ForwardDiff\N0wMF\src\apiutils.jl:37
 [9] gradient at C:\Users\vsarn\.julia\packages\ForwardDiff\N0wMF\src\gradient.jl:17 [inlined]
 [10] #94 at C:\Users\vsarn\.julia\packages\ForwardDiff\N0wMF\src\hessian.jl:16 [inlined]
 [11] vector_mode_dual_eval(::getfield(ForwardDiff, Symbol("##94#95")){typeof(remlf2),ForwardDiff.HessianConfig{ForwardDiff.Tag{typeof(remlf2),Float64},Float64,5,Array{ForwardDiff.Dual{ForwardDiff.Tag{typeof(remlf2),Float64},ForwardDiff.Dual{ForwardDiff.Tag{typeof(remlf2),Float64},Float64,5},5},1},Array{ForwardDiff.Dual{ForwardDiff.Tag{typeof(remlf2),Float64},Float64,5},1}}}, ::Array{Float64,1}, ::ForwardDiff.JacobianConfig{ForwardDiff.Tag{typeof(remlf2),Float64},Float64,5,Array{ForwardDiff.Dual{ForwardDiff.Tag{typeof(remlf2),Float64},Float64,5},1}}) at C:\...\.julia\packages\ForwardDiff\N0wMF\src\apiutils.jl:37
 [12] vector_mode_jacobian(::getfield(ForwardDiff, Symbol("##94#95")){typeof(remlf2),ForwardDiff.HessianConfig{ForwardDiff.Tag{typeof(remlf2),Float64},Float64,5,Array{ForwardDiff.Dual{ForwardDiff.Tag{typeof(remlf2),Float64},ForwardDiff.Dual{ForwardDiff.Tag{typeof(remlf2),Float64},Float64,5},5},1},Array{ForwardDiff.Dual{ForwardDiff.Tag{typeof(remlf2),Float64},Float64,5},1}}}, ::Array{Float64,1}, ::ForwardDiff.JacobianConfig{ForwardDiff.Tag{typeof(remlf2),Float64},Float64,5,Array{ForwardDiff.Dual{ForwardDiff.Tag{typeof(remlf2),Float64},Float64,5},1}}) at C:\...\.julia\packages\ForwardDiff\N0wMF\src\jacobian.jl:140
 [13] jacobian(::Function, ::Array{Float64,1}, ::ForwardDiff.JacobianConfig{ForwardDiff.Tag{typeof(remlf2),Float64},Float64,5,Array{ForwardDiff.Dual{ForwardDiff.Tag{typeof(remlf2),Float64},Float64,5},1}}, ::Val{false}) at C:\...\.julia\packages\ForwardDiff\N0wMF\src\jacobian.jl:17
 [14] hessian(::Function, ::Array{Float64,1}, ::ForwardDiff.HessianConfig{ForwardDiff.Tag{typeof(remlf2),Float64},Float64,5,Array{ForwardDiff.Dual{ForwardDiff.Tag{typeof(remlf2),Float64},ForwardDiff.Dual{ForwardDiff.Tag{typeof(remlf2),Float64},Float64,5},5},1},Array{ForwardDiff.Dual{ForwardDiff.Tag{typeof(remlf2),Float64},Float64,5},1}}, ::Val{true}) at C:\...\.julia\packages\ForwardDiff\N0wMF\src\hessian.jl:17
 [15] hessian(::Function, ::Array{Float64,1}, ::ForwardDiff.HessianConfig{ForwardDiff.Tag{typeof(remlf2),Float64},Float64,5,Array{ForwardDiff.Dual{ForwardDiff.Tag{typeof(remlf2),Float64},ForwardDiff.Dual{ForwardDiff.Tag{typeof(remlf2),Float64},Float64,5},5},1},Array{ForwardDiff.Dual{ForwardDiff.Tag{typeof(remlf2),Float64},Float64,5},1}}) at C:\...\.julia\packages\ForwardDiff\N0wMF\src\hessian.jl:15 (repeats 2 times)
 [16] top-level scope at none:0

ForwardDiff works by operating on dual numbers:

http://www.juliadiff.org/ForwardDiff.jl/latest/dev/how_it_works.html

As I said, you should either consider using non-modifying versions of these functions, or (if you are feeling adventurous) pre-allocate buffers of the correct type. It is better to do this for an inner loop only, where the outer part will know the type (Dual has a type tag).

To help with the latter, it would be great to have a full MWE.

1 Like
temp = zeros(promote_type(eltype(first(Xv)),eltype(iV)),p,p)
1 Like