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?
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.
temp = zeros(promote_type(eltype(first(Xv)),eltype(iV)),p,p)