My function works on it's own but not in NonLinearSolve

General Usage
,
1

I have a set of equations that I would like to solve using the non linear equation solver NonLinearSolve.jl. I have writen a function for this which runs fine on it’s own but when I try to solve it I get an error message.

The equation I would like to solve is the following:

cd("C:/Users/...correct directory etc")

using NonlinearSolve
using LinearAlgebra
using JLD2

mutable struct parameters

    n::Int64
    matrix::Matrix{Float64}
    tensor::Array{Float64}
    norm::Vector{Float64}

end

function equation!(y, z, p::parameters)

    @inbounds for i = 1:p.n
        for j = 1:i
            acc = zero(eltype(p.matrix))
            for k = 1:p.n
                acctmp = zero(eltype(p.matrix))
                @simd for l = 1:p.n
                    acctmp += p.tensor[l, k, j, i] * z[p.n+l]
                end
                acc += acctmp * z[k]
            end
            p.matrix[i, p.n+j] = p.matrix[j, p.n+i] = p.matrix[p.n+i, j] = p.matrix[p.n+j, i] = acc
        end
    end

    p.matrix[p.n+1:2p.n, 1:p.n] .*= p.norm / (z[1] * z[p.n+1])
    p.matrix[1:p.n, p.n+1:2p.n] .*= p.norm / (z[1] * z[p.n+1])

    y[1:2p.n] = (p.matrix - I * z[2p.n+1]) * z[1:2p.n]

end

load_object(string("example u.jld2"))
load_object(string("example p.jld2"))

prob = NonlinearProblem(equation!, u, p)

sol = solve(prob)

Now, I tried to upload the example u and p files to this post but that does not seem to be possible. I hope perhaps someone can just see the problem without needing those files ore perhaps knows how to upload them. :sweat_smile:

When I run this I get the following error message:

ERROR: MethodError: no method matching Float64(::ForwardDiff.Dual{ForwardDiff.Tag{NonlinearSolve.NonlinearSolveTag, Float64}, Float64, 11})

Closest candidates are:
  (::Type{T})(::Real, ::RoundingMode) where T<:AbstractFloat
   @ Base rounding.jl:207
  (::Type{T})(::T) where T<:Number
   @ Core boot.jl:792
  (::Type{T})(::AbstractChar) where T<:Union{AbstractChar, Number}
   @ Base char.jl:50
  ...

With the following stacktrace:

Stacktrace:
  [1] convert(#unused#::Type{Float64}, x::ForwardDiff.Dual{ForwardDiff.Tag{NonlinearSolve.NonlinearSolveTag, Float64}, Float64, 11})
    @ Base .\number.jl:7
  [2] setindex!
    @ .\array.jl:971 [inlined]
  [3] equation!(y::Vector{ForwardDiff.Dual{ForwardDiff.Tag{NonlinearSolve.NonlinearSolveTag, Float64}, Float64, 11}}, z::Vector{ForwardDiff.Dual{ForwardDiff.Tag{NonlinearSolve.NonlinearSolveTag, Float64}, Float64, 11}}, p::parameters)
    @ Main c:\Users\laura\Documents\Universiteit\Doctoraat\Impurities in Fermi Fluids\Julia\example code.jl:29
  [4] NonlinearFunction
    @ C:\Users\laura\.julia\packages\SciMLBase\2HZ5m\src\scimlfunctions.jl:2356 [inlined]
  [5] JacobianWrapper
    @ C:\Users\laura\.julia\packages\SciMLBase\2HZ5m\src\function_wrappers.jl:99 [inlined]
  [6] chunk_mode_jacobian!(result::Matrix{Float64}, f!::SciMLBase.JacobianWrapper{true, NonlinearFunction{true, SciMLBase.FullSpecialize, typeof(equation!), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED_NO_TIME), Nothing, SymbolicIndexingInterface.SymbolCache{Nothing, Nothing, Nothing}, Nothing}, parameters}, y::Vector{Float64}, x::Vector{Float64}, cfg::ForwardDiff.JacobianConfig{ForwardDiff.Tag{NonlinearSolve.NonlinearSolveTag, Float64}, Float64, 11, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{NonlinearSolve.NonlinearSolveTag, Float64}, Float64, 11}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{NonlinearSolve.NonlinearSolveTag, Float64}, Float64, 11}}}})
    @ ForwardDiff C:\Users\laura\.julia\packages\ForwardDiff\PcZ48\src\jacobian.jl:183
  [7] jacobian!
    @ C:\Users\laura\.julia\packages\ForwardDiff\PcZ48\src\jacobian.jl:80 [inlined]
  [8] jacobian!
    @ C:\Users\laura\.julia\packages\ForwardDiff\PcZ48\src\jacobian.jl:76 [inlined]
  [9] sparse_jacobian!
    @ C:\Users\laura\.julia\packages\SparseDiffTools\CsonI\src\highlevel\forward_mode.jl:70 [inlined]
 [10] #jacobian!!#130
    @ C:\Users\laura\.julia\packages\NonlinearSolve\KlGj2\src\jacobian.jl:62 [inlined]
 [11] jacobian!!
    @ C:\Users\laura\.julia\packages\NonlinearSolve\KlGj2\src\jacobian.jl:53 [inlined]
 [12] perform_step!(cache::NonlinearSolve.BroydenCache{true, :true_jacobian, :good_broyden, NonlinearFunction{true, SciMLBase.FullSpecialize, typeof(equation!), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED_NO_TIME), Nothing, SymbolicIndexingInterface.SymbolCache{Nothing, Nothing, Nothing}, Nothing}, Broyden{:true_jacobian, :good_broyden, true, AutoForwardDiff{nothing, ForwardDiff.Tag{NonlinearSolve.NonlinearSolveTag, Float64}}, Nothing, LineSearch{Nothing, Nothing, Bool}, Nothing}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, parameters, SciMLBase.JacobianWrapper{true, NonlinearFunction{true, SciMLBase.FullSpecialize, typeof(equation!), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED_NO_TIME), Nothing, SymbolicIndexingInterface.SymbolCache{Nothing, Nothing, Nothing}, Nothing}, parameters}, Matrix{Float64}, Nothing, Vector{Float64}, Float64, Nothing, typeof(DiffEqBase.NONLINEARSOLVE_DEFAULT_NORM), Float64, Float64, Float64, NonlinearSolve.var"#118#119"{Float64}, SparseDiffTools.ForwardDiffJacobianCache{SparseDiffTools.NoMatrixColoring, ForwardDiff.JacobianConfig{ForwardDiff.Tag{NonlinearSolve.NonlinearSolveTag, Float64}, Float64, 11, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{NonlinearSolve.NonlinearSolveTag, Float64}, Float64, 11}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{NonlinearSolve.NonlinearSolveTag, Float64}, Float64, 11}}}}, Nothing, Vector{Float64}, Vector{Float64}}, NonlinearProblem{Vector{Float64}, true, parameters, NonlinearFunction{true, SciMLBase.FullSpecialize, typeof(equation!), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED_NO_TIME), Nothing, SymbolicIndexingInterface.SymbolCache{Nothing, Nothing, Nothing}, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardNonlinearProblem}, NonlinearSolve.NoLineSearchCache{Float64}, DiffEqBase.NonlinearTerminationModeCache{Vector{Float64}, Float64, AbsSafeBestTerminationMode{Nothing, Int64, Float64}, Float64, Vector{Float64}, Nothing}, NonlinearSolve.NonlinearSolveTrace{false, false, TraceMinimal, Nothing}})
    @ NonlinearSolve C:\Users\laura\.julia\packages\NonlinearSolve\KlGj2\src\broyden.jl:171
 [13] solve!(cache::NonlinearSolve.BroydenCache{true, :true_jacobian, :good_broyden, NonlinearFunction{true, SciMLBase.FullSpecialize, typeof(equation!), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED_NO_TIME), Nothing, SymbolicIndexingInterface.SymbolCache{Nothing, Nothing, Nothing}, Nothing}, Broyden{:true_jacobian, :good_broyden, true, AutoForwardDiff{nothing, ForwardDiff.Tag{NonlinearSolve.NonlinearSolveTag, Float64}}, Nothing, LineSearch{Nothing, Nothing, Bool}, Nothing}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, parameters, SciMLBase.JacobianWrapper{true, NonlinearFunction{true, SciMLBase.FullSpecialize, typeof(equation!), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED_NO_TIME), Nothing, SymbolicIndexingInterface.SymbolCache{Nothing, Nothing, Nothing}, Nothing}, parameters}, Matrix{Float64}, Nothing, Vector{Float64}, Float64, Nothing, typeof(DiffEqBase.NONLINEARSOLVE_DEFAULT_NORM), Float64, Float64, Float64, NonlinearSolve.var"#118#119"{Float64}, SparseDiffTools.ForwardDiffJacobianCache{SparseDiffTools.NoMatrixColoring, ForwardDiff.JacobianConfig{ForwardDiff.Tag{NonlinearSolve.NonlinearSolveTag, Float64}, Float64, 11, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{NonlinearSolve.NonlinearSolveTag, Float64}, Float64, 11}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{NonlinearSolve.NonlinearSolveTag, Float64}, Float64, 11}}}}, Nothing, Vector{Float64}, Vector{Float64}}, NonlinearProblem{Vector{Float64}, true, parameters, NonlinearFunction{true, SciMLBase.FullSpecialize, typeof(equation!), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED_NO_TIME), Nothing, SymbolicIndexingInterface.SymbolCache{Nothing, Nothing, Nothing}, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardNonlinearProblem}, NonlinearSolve.NoLineSearchCache{Float64}, DiffEqBase.NonlinearTerminationModeCache{Vector{Float64}, Float64, AbsSafeBestTerminationMode{Nothing, Int64, Float64}, Float64, Vector{Float64}, Nothing}, NonlinearSolve.NonlinearSolveTrace{false, false, TraceMinimal, Nothing}})
    @ NonlinearSolve C:\Users\laura\.julia\packages\NonlinearSolve\KlGj2\src\NonlinearSolve.jl:149
 [14] __solve(::NonlinearProblem{Vector{Float64}, true, parameters, NonlinearFunction{true, SciMLBase.FullSpecialize, typeof(equation!), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED_NO_TIME), Nothing, SymbolicIndexingInterface.SymbolCache{Nothing, Nothing, Nothing}, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardNonlinearProblem}, ::Broyden{:true_jacobian, :good_broyden, true, Nothing, Nothing, LineSearch{Nothing, Nothing, Bool}, Nothing}; kwargs::Base.Pairs{Symbol, Bool, Tuple{Symbol, Symbol}, NamedTuple{(:default_set, :second_time), Tuple{Bool, Bool}}})
    @ NonlinearSolve C:\Users\laura\.julia\packages\NonlinearSolve\KlGj2\src\NonlinearSolve.jl:135
 [15] __solve
    @ C:\Users\laura\.julia\packages\NonlinearSolve\KlGj2\src\NonlinearSolve.jl:132 [inlined]
 [16] macro expansion
    @ C:\Users\laura\.julia\packages\NonlinearSolve\KlGj2\src\default.jl:123 [inlined]
 [17] __solve(::NonlinearProblem{Vector{Float64}, true, parameters, NonlinearFunction{true, SciMLBase.FullSpecialize, typeof(equation!), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED_NO_TIME), Nothing, SymbolicIndexingInterface.SymbolCache{Nothing, Nothing, Nothing}, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardNonlinearProblem}, ::NonlinearSolvePolyAlgorithm{:NLS, 7, Tuple{Broyden{:identity, :good_broyden, false, Nothing, Nothing, LineSearch{Nothing, Nothing, Bool}, Nothing}, Broyden{:true_jacobian, :good_broyden, true, Nothing, Nothing, LineSearch{Nothing, Nothing, Bool}, Nothing}, Klement{:identity, false, Nothing, Nothing, typeof(NonlinearSolve.DEFAULT_PRECS), LineSearch{Nothing, Nothing, Bool}, Bool}, NewtonRaphson{nothing, Nothing, Nothing, typeof(NonlinearSolve.DEFAULT_PRECS), LineSearch{Nothing, Nothing, Bool}}, NewtonRaphson{nothing, Nothing, Nothing, typeof(NonlinearSolve.DEFAULT_PRECS), LineSearch{BackTracking{Float64, Int64}, Nothing, Bool}}, TrustRegion{nothing, Nothing, Rational{Int64}, Nothing, typeof(NonlinearSolve.DEFAULT_PRECS), Rational{Int64}, Nothing}, TrustRegion{nothing, Nothing, Rational{Int64}, Nothing, typeof(NonlinearSolve.DEFAULT_PRECS), Rational{Int64}, Nothing}}}; kwargs::Base.Pairs{Symbol, Bool, Tuple{Symbol, Symbol}, NamedTuple{(:default_set, :second_time), Tuple{Bool, Bool}}})
    @ NonlinearSolve C:\Users\laura\.julia\packages\NonlinearSolve\KlGj2\src\default.jl:115
 [18] __solve
    @ C:\Users\laura\.julia\packages\NonlinearSolve\KlGj2\src\default.jl:115 [inlined]
 [19] #__solve#196
    @ C:\Users\laura\.julia\packages\NonlinearSolve\KlGj2\src\default.jl:386 [inlined]
 [20] __solve
    @ C:\Users\laura\.julia\packages\NonlinearSolve\KlGj2\src\default.jl:383 [inlined]
 [21] #__solve#61
    @ C:\Users\laura\.julia\packages\DiffEqBase\eTCPy\src\solve.jl:1370 [inlined]
 [22] __solve
    @ C:\Users\laura\.julia\packages\DiffEqBase\eTCPy\src\solve.jl:1363 [inlined]
 [23] #solve_call#34
    @ C:\Users\laura\.julia\packages\DiffEqBase\eTCPy\src\solve.jl:608 [inlined]
 [24] solve_call
    @ C:\Users\laura\.julia\packages\DiffEqBase\eTCPy\src\solve.jl:566 [inlined]
 [25] #solve_up#42
    @ C:\Users\laura\.julia\packages\DiffEqBase\eTCPy\src\solve.jl:1049 [inlined]
 [26] solve_up
    @ C:\Users\laura\.julia\packages\DiffEqBase\eTCPy\src\solve.jl:1043 [inlined]
 [27] solve(::NonlinearProblem{Vector{Float64}, true, parameters, NonlinearFunction{true, SciMLBase.FullSpecialize, typeof(equation!), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED_NO_TIME), Nothing, SymbolicIndexingInterface.SymbolCache{Nothing, Nothing, Nothing}, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardNonlinearProblem}; sensealg::Nothing, u0::Nothing, p::Nothing, wrap::Val{true}, kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ DiffEqBase C:\Users\laura\.julia\packages\DiffEqBase\eTCPy\src\solve.jl:1037
 [28] solve(::NonlinearProblem{Vector{Float64}, true, parameters, NonlinearFunction{true, SciMLBase.FullSpecialize, typeof(equation!), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED_NO_TIME), Nothing, SymbolicIndexingInterface.SymbolCache{Nothing, Nothing, Nothing}, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardNonlinearProblem})
    @ DiffEqBase C:\Users\laura\.julia\packages\DiffEqBase\eTCPy\src\solve.jl:1027
 [29] top-level scope
    @ c:\Users\laura\Documents\Universiteit\Doctoraat\Impurities in Fermi Fluids\Julia\example code.jl:45

In my limited understanding of what is happening here, the non linear solver tries to calculate the jacobian of my function by propagating some object of the type ::ForewardDiff through my function which does not work for some reason.

Does anyone know if there is a way to change my function in such a way that it is possible for the solver to propagate this ::ForewardDiff object through it or should I attempt to make a seperate jacobian function?

Or perhaps the problem is something else all together?

Best regards,
Laura

2

It’s the same issue as:

You’ll want to make use of PreallocationTools.jl if you’re precaching values that are then being autodiff’d:

3

In this case your code seems fully generic with respect to element types, the only thing you need to do is allow the same genericity in your struct fields:

mutable struct parameters{T1,T2,T3}
    n::Int64
    matrix::Matrix{T1}
    tensor::Array{T2}
    norm::Vector{T3}
end

You probably also want to specify the order of the tensor to make the field type-inferrable (at the moment Julia doesn’t know how many dimensions the Array has).

4

That doesn’t solve it because the p is constructed outside of the function, and so it needs to be generic to the multiple cache sizes. matrix, tensor, and norm need to be PreallocationTools.DualCache objects.

5

Making the the element types generic and specifying the array dimensions does indeed not seem to solve this. I’ll have a look at the Autodifferentiation and Dual Numbers.

6

The concrete problem is the acc and acctmp variables, which you set to zero(eltype(p.matrix)), i.e. a zero of type Float64. In the loop you create the acctmp by multiplying with z[...]. Now, when differentiating, both y and z contains Numbers of type ForwardDiff.Dual (which can be seen near the top of the stack dump). Now, a Dual can be multiplied by Float64 to yield a new Dual. This makes acctmp a Dual, and then acc becomes a Dual. Then you try to put it into p.matrix[...], which contains Float64. This fails, that’s the error message, because a Dual can’t be converted to a Float64.

To work with autodiff, your function can’t assume the inputs are Float64, it must return a ForwardDiff.Dual when the inputs are ForwardDiff.Dual.

7

Ah, I think I understand what you mean.

So if for example, my function made a new temporary matrix every itteration instead of storing those values in the p.matrix then it would all work. The reason I made it store the temporary values in p.matrix was to save on allocations if that makes sense.

Or alternatively I turn p.matrix into a Dual numer. But won’t that cause problems when I just try to evaluate the function on it’s own?

8

Yes, you’re right. Allocating would solve it, and making p.matrix a Dual will create problems with ordinary evaluations.

This is what Chris suggested a solution to above. The PreallocationTools package (which I haven’t used) apparently can help with this problem.

9

It’s a little hard to tell exactly what your code does without any example data, but here is a (rather nonsensical) example of how to use PreallocationTools.jl with NonlinearSolve.jl:

using NonlinearSolve
using LinearAlgebra
using PreallocationTools

struct Parameters{T}
    matrix::T
end

function equation!(y, z, p::Parameters)

    # Cached Matrix 
    cmatrix = PreallocationTools.get_tmp(p.matrix, z)

    for i ∈ axes(cmatrix, 1), j ∈ axes(cmatrix, 2)
        cmatrix[i,j] = exp( z[i] + 2*z[j] )
    end

    for i ∈ eachindex(y)
        y[i] = 2.5 - sum( cmatrix[i,j]*z[j] for j ∈ eachindex(z) )
    end
    
    return nothing
end

u = zeros(3)
m = zeros(3,3)
cmatrix = PreallocationTools.DiffCache(m)
p = Parameters(cmatrix)

prob = NonlinearProblem(equation!, u, p)

sol = solve(prob)
10

This ended up almost working for me except before converting p.matrix using get_tmp I first had to turn it into a DiffCache object.

Using my own example it would look like the code below. I will probably still fiffle around with it a bit because I have the feeling I could further optimize this by passing the p.matrix as a DiffCache from the beginning.

using NonlinearSolve
using LinearAlgebra
using JLD2
using PreallocationTools

mutable struct parameters{T}

    n::Int64
    norm::Vector{Float64}
    matrix::Matrix{T}
    tensor::Array{Float64,4}

end

function equation!(y, z, p::parameters)

    cache_matrix = PreallocationTools.get_tmp(DiffCache(p.matrix), z)

    @inbounds for i = 1:p.n
        for j = 1:i
            acc = zero(eltype(p.matrix))
            for k = 1:p.n
                acctmp = zero(eltype(p.matrix))
                @simd for l = 1:p.n
                    acctmp += p.tensor[l, k, j, i] * z[p.n+l]
                end
                acc += acctmp * z[k]
            end
            cache_matrix[i, p.n+j] = cache_matrix[j, p.n+i] = cache_matrix[p.n+i, j] = cache_matrix[p.n+j, i] = acc
        end
    end

    cache_matrix[p.n+1:2p.n, 1:p.n] .*= p.norm / (z[1] * z[p.n+1])
    cache_matrix[1:p.n, p.n+1:2p.n] .*= p.norm / (z[1] * z[p.n+1])

    y[1:2p.n] = (cache_matrix - I * z[2p.n+1]) * z[1:2p.n]

end

u = load_object(string("example u.jld2"))
p = load_object(string("example p.jld2"))

prob = NonlinearProblem(equation!, u, p)

sol = solve(prob)
11

Right, you won’t get much benefit out of PreallocationTools if you’re recreating the DiffCache at every iteration. The following might work for you:

using Accessors
# Functions for the problem
# [...]
u = load_object(string("example u.jld2"))
p_in = load_object(string("example p.jld2"))
p = @set p_in.matrix = PreallocationTools.DiffCache(p_in.matrix)
prob = NonlinearProblem(equation!, u, p)
sol = solve(prob)

and changing the line getting the cache to:

cache_matrix = PreallocationTools.get_tmp(p.matrix, z)