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

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.

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

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:

1 Like

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).

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.

1 Like

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.

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 `Number`s 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`.

1 Like

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?

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.

1 Like

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)
``````
1 Like

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)
``````

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)
``````
1 Like