ForwardDiff.jl: ERROR: LoadError: MethodError: no method matching Float64(::ForwardDiff.Dual{ForwardDiff.Tag{Zygote.var"#104#105"{var"#f#1"}, Float64}, Float64, 1})

I had a lot of problems with ForwardDiff. I used Zygote to do automatic differential calculations before, so I didn’t know much about ForwardDiff. My original intention was to use Zygote to achieve the gradient of gradients, but there were a lot of problems with twice automatic differentiation in reverse mode. So I used Zygote.hessian to achieve the mixing forward + reverse mode, and then many problems occurred.

I have solved some problems by myself, and some problems have been solved by useful help. However, after solving one problem, a new problem appears. This is my latest problem:

using GenericLinearAlgebra
using LinearAlgebra
using Zygote
using OMEinsum 
using Random
using ForwardDiff, ChainRulesCore
Base.convert(::Type{ForwardDiff.Dual{T, V, N}}, ::ChainRulesCore.AbstractZero) where {T, V, N} = ForwardDiff.Dual{T, V, N}(zero(V))
function ff(x,C,E)
    D = size(C,1)
    W = exp.(- x * [1 -1; -1 1])
    U, S,_ = GenericLinearAlgebra.svd(W)
    W = U * Diagonal(sqrt.(S))
    @ein T[i, j, k, l] := W[a, i] * W[a, j] * W[a, k] * W[a, l] 
    for i in 1:100
        @ein CP[c, l , d, j] := ((C[a, b] * E[a, d, k]) * E[b, c, i]) * T[i, j, k, l]
        CPm = reshape(CP, 2*D, 2*D)
        U, S,_= GenericLinearAlgebra.svd(CPm)
        Z = reshape(U[:, 1:D], D, 2,D)
        @ein C[m, n] :=Z[d, j, n] * (CP[c, l, d, j] * Z[c, l, m])
        @ein E[c,d,j] := ((Z[a,k,c] * E[a,b,i]) * T[i,j,k,l] )* Z[b,l,d]
        C += C'
        E += permutedims( E,(2,1,3))
        C /= maximum(abs.(C))
        E /= maximum(abs.(E))   
    end
    @ein CP[c, l , d, j] := ((C[a, b] * E[a, d, k]) * E[b, c, i]) * T[i, j, k, l]
    @ein CECE[a,d,g,f] := ((C[a,b] * E[b,c,d]) * C[c,e]) * E[e,f,g]
    @ein CECEC[a,d,g,h] := CECE[a,d,g,f] * C[f,h]
    @ein Z[] := CP[c, l , d, j] * CECEC[d,j,l,c]
    return Z[]
end

function test()
    x0=0.5
    D=20
    Random.seed!(123)
    C=rand(D,D)
    Random.seed!(1234)
    E=rand(D,D,2)
    C += C'
    E += permutedims( E,(2,1,3))
    C /= maximum(abs.(C))
    E /= maximum(abs.(E))  
    f(x) = ff(x,C,E)
    @show f(x0)
    @show gradient(f,x0)
    @show hessian(f,x0)
end

test()



the output is

f(x0) = 4.849582676424554
gradient(f, x0) = (21.480979163911602,)
ERROR: LoadError: MethodError: no method matching Float64(::ForwardDiff.Dual{ForwardDiff.Tag{Zygote.var"#104#105"{var"#f#1"}, Float64}, Float64, 1})
Closest candidates are:
  (::Type{T})(::Real, ::RoundingMode) where T<:AbstractFloat at E:\AppData\Julia-1.7.3\share\julia\base\rounding.jl:200
  (::Type{T})(::T) where T<:Number at E:\AppData\Julia-1.7.3\share\julia\base\boot.jl:770
  (::Type{T})(::AbstractChar) where T<:Union{AbstractChar, Number} at E:\AppData\Julia-1.7.3\share\julia\base\char.jl:50
  ...

I’m sorry that my example might be a little complicated, but I don’t know how to simplify again and get the same error.

I don’t know how to solve this problem, and I don’t know what new problems will arise after solving this problem. ForwardDiff makes me feel bad

I got a simpler example:

using Zygote
using OMEinsum 
function ff(x,C,E)
    W = [1 exp(x);exp(-x) -1]
    @ein T[i, j, k, l] := W[a, i] * W[a, j] * W[a, k] * W[a, l] 
    @ein CP[c, l , d, j] := ((C[a, b] * E[a, d, k]) * E[b, c, i]) * T[i, j, k, l]
    @ein Z[] := CP[c, l , c, j] * T[d,j,l,d]
    return Z[]
end

function test()

    C=rand(20,20)
    E=rand(20,20,2)
    f(x) = ff(x,C,E)
    x0=0.5
    @show f(x0)
    @show gradient(f,x0)
    @show hessian(f,x0)
end

test()


then receive the same error:

f(x0) = 370889.75334902847
gradient(f, x0) = (2.092207069892614e6,)
ERROR: LoadError: MethodError: no method matching Float64(::ForwardDiff.Dual{ForwardDiff.Tag{Zygote.var"#104#105"{var"#f#1"{Array{Float64, 3}, Matrix{Float64}}}, Float64}, Float64, 1})
Closest candidates are:
  (::Type{T})(::Real, ::RoundingMode) where T<:AbstractFloat at E:\AppData\Julia-1.7.3\share\julia\base\rounding.jl:200
  (::Type{T})(::T) where T<:Number at E:\AppData\Julia-1.7.3\share\julia\base\boot.jl:770
  (::Type{T})(::AbstractChar) where T<:Union{AbstractChar, Number} at E:\AppData\Julia-1.7.3\share\julia\base\char.jl:50
  ...

The problem occurs on this line , and I don’t know why. How do I fix this?

This looks like an incompatibility between OMEinsum and ForwardDiff. The problem is that CP is allocated as an array of Float64, it needs to be allocated as an array of Dual (usually by using similar and one of the input variables). I am not familiar with OMEinsum and didn’t immediately see a way to control allocations, but perhaps you could allocate CP manually before calling the summation?

I’ll get no error once I write it in this form.

CP = ein"ab,adk,bci,ijkl -> cldj"(C,E,E,T)

The results of these two forms are the same, but the implementation structure seems to be different, I can’t understand the specific

This will also result in an error:

CP = ein"(ab,adk),bci,ijkl -> cldj"(C,E,E,T)

I think I made a mistake. It’s not a question of the two forms of OMEinsum. It’s the fact that I put in parentheses, then create a cache which is not compatible with autodifferentiaion via ForwardDiff.jl