# How can I use NonlinearSolve.jl with matrix functions?

I want to use NonlinearSolve.jl with matrix functions, something like:

using NonlinearSolve
using LinearAlgebra

dim = 2
function f(u, p)
return exp(u) - Matrix{Float64}(I(dim))
end

u0 = rand(dim, dim) * 0.01
p = nothing
prob = NonlinearProblem(f, u0, p)
sol = solve(prob)
display(sol)
err = sol - zeros(dim, dim)
display(err)


This, however, casts an error:

ERROR: MethodError: no method matching exp!(::Matrix{ForwardDiff.Dual{ForwardDiff.Tag{NonlinearSolve.NonlinearSolveTag, Float64}, Float64, 4}})

Closest candidates are:
exp!(::StridedMatrix{T}) where T<:Union{Float32, Float64, ComplexF64, ComplexF32}
@ LinearAlgebra ~/.julia/juliaup/julia-1.10.2+0.x64.apple.darwin14/share/julia/stdlib/v1.10/LinearAlgebra/src/dense.jl:649

Stacktrace:
[1] exp(A::Matrix{ForwardDiff.Dual{ForwardDiff.Tag{NonlinearSolve.NonlinearSolveTag, Float64}, Float64, 4}})
@ LinearAlgebra ~/.julia/juliaup/julia-1.10.2+0.x64.apple.darwin14/share/julia/stdlib/v1.10/LinearAlgebra/src/dense.jl:594
[2] f(u::Matrix{ForwardDiff.Dual{ForwardDiff.Tag{NonlinearSolve.NonlinearSolveTag, Float64}, Float64, 4}}, p::Nothing)
@ Main ./Untitled-1:6
[3] (::NonlinearFunction{β¦})(::Matrix{β¦}, ::Vararg{β¦})
@ SciMLBase ~/.julia/packages/SciMLBase/QSc1r/src/scimlfunctions.jl:2197
[4] (::SciMLBase.JacobianWrapper{false, NonlinearFunction{β¦}, Nothing})(u::Matrix{ForwardDiff.Dual{β¦}})
@ SciMLBase ~/.julia/packages/SciMLBase/QSc1r/src/function_wrappers.jl:97
[5] vector_mode_dual_eval!(f::SciMLBase.JacobianWrapper{β¦}, cfg::ForwardDiff.JacobianConfig{β¦}, x::Matrix{β¦})
@ ForwardDiff ~/.julia/packages/ForwardDiff/PcZ48/src/apiutils.jl:24
[6] vector_mode_jacobian!(result::Matrix{β¦}, f::SciMLBase.JacobianWrapper{β¦}, x::Matrix{β¦}, cfg::ForwardDiff.JacobianConfig{β¦})
@ ForwardDiff ~/.julia/packages/ForwardDiff/PcZ48/src/jacobian.jl:145
[7] jacobian!(result::Matrix{β¦}, f::Function, x::Matrix{β¦}, cfg::ForwardDiff.JacobianConfig{β¦}, ::Val{β¦})
@ ForwardDiff ~/.julia/packages/ForwardDiff/PcZ48/src/jacobian.jl:58
[8] jacobian!(result::Matrix{β¦}, f::Function, x::Matrix{β¦}, cfg::ForwardDiff.JacobianConfig{β¦})
@ ForwardDiff ~/.julia/packages/ForwardDiff/PcZ48/src/jacobian.jl:56
[9] sparse_jacobian!(J::Matrix{β¦}, ::AutoForwardDiff{β¦}, cache::SparseDiffTools.ForwardDiffJacobianCache{β¦}, f::SciMLBase.JacobianWrapper{β¦}, x::Matrix{β¦})
@ SparseDiffTools ~/.julia/packages/SparseDiffTools/qxnHN/src/highlevel/forward_mode.jl:60
[10] (::NonlinearSolve.JacobianCache{β¦})(J::Matrix{β¦}, u::Matrix{β¦}, p::Nothing)
@ NonlinearSolve ~/.julia/packages/NonlinearSolve/5yLII/src/internal/jacobian.jl:136
[11] JacobianCache
@ ~/.julia/packages/NonlinearSolve/5yLII/src/internal/jacobian.jl:106 [inlined]
[12] __step!(cache::NonlinearSolve.GeneralizedFirstOrderAlgorithmCache{β¦}; recompute_jacobian::Nothing, kwargs::@Kwargs{})
@ NonlinearSolve ~/.julia/packages/NonlinearSolve/5yLII/src/core/generalized_first_order.jl:207
[13] __step!
@ ~/.julia/packages/NonlinearSolve/5yLII/src/core/generalized_first_order.jl:203 [inlined]
[14] #step!#210
@ ~/.julia/packages/NonlinearSolve/5yLII/src/core/generic.jl:55 [inlined]
[15] step!
@ ~/.julia/packages/NonlinearSolve/5yLII/src/core/generic.jl:50 [inlined]
[16] solve!(cache::NonlinearSolve.GeneralizedFirstOrderAlgorithmCache{β¦})
@ NonlinearSolve ~/.julia/packages/NonlinearSolve/5yLII/src/core/generic.jl:13
[17] __solve(::NonlinearProblem{β¦}, ::GeneralizedFirstOrderAlgorithm{β¦}; kwargs::@Kwargs{β¦})
@ NonlinearSolve ~/.julia/packages/NonlinearSolve/5yLII/src/core/generic.jl:4
[18] __solve
@ ~/.julia/packages/NonlinearSolve/5yLII/src/core/generic.jl:1 [inlined]
[19] macro expansion
@ ~/.julia/packages/NonlinearSolve/5yLII/src/default.jl:278 [inlined]
[20] __solve(::NonlinearProblem{β¦}, ::NonlinearSolvePolyAlgorithm{β¦}; alias_u0::Bool, verbose::Bool, kwargs::@Kwargs{β¦})
@ NonlinearSolve ~/.julia/packages/NonlinearSolve/5yLII/src/default.jl:248
[21] __solve
@ ~/.julia/packages/NonlinearSolve/5yLII/src/default.jl:248 [inlined]
[22] #__solve#329
@ ~/.julia/packages/NonlinearSolve/5yLII/src/default.jl:511 [inlined]
[23] __solve
@ ~/.julia/packages/NonlinearSolve/5yLII/src/default.jl:508 [inlined]
[24] #__solve#61
@ ~/.julia/packages/DiffEqBase/nc5nE/src/solve.jl:1371 [inlined]
[25] __solve
@ ~/.julia/packages/DiffEqBase/nc5nE/src/solve.jl:1364 [inlined]
[26] #solve_call#34
@ ~/.julia/packages/DiffEqBase/nc5nE/src/solve.jl:609 [inlined]
[27] solve_call
@ ~/.julia/packages/DiffEqBase/nc5nE/src/solve.jl:567 [inlined]
[28] #solve_up#42
@ ~/.julia/packages/DiffEqBase/nc5nE/src/solve.jl:1050 [inlined]
[29] solve_up
@ ~/.julia/packages/DiffEqBase/nc5nE/src/solve.jl:1044 [inlined]
[30] #solve#41
@ ~/.julia/packages/DiffEqBase/nc5nE/src/solve.jl:1038 [inlined]
[31] solve(::NonlinearProblem{β¦})
@ DiffEqBase ~/.julia/packages/DiffEqBase/nc5nE/src/solve.jl:1028
[32] top-level scope
@ Untitled-1:12
Some type information was truncated. Use show(err) to see complete types.


My attempts were:

• cast to float
function f(u, p)
return exp(Matrix{Float64}(u)) - Matrix{Float64}(I(dim))
end

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

Closest candidates are:
(::Type{T})(::Real, ::RoundingMode) where T<:AbstractFloat
@ Base rounding.jl:207
(::Type{T})(::T) where T<:Number
@ Core boot.jl:792
Float64(::IrrationalConstants.FourinvΟ)
@ IrrationalConstants ~/.julia/packages/IrrationalConstants/vp5v4/src/macro.jl:112

• use value
function f(u, p)
return exp(u.value) - Matrix{Float64}(I(dim))
end

ERROR: type Array has no field value


How can I fix this?

1 Like

Matrix exponential is not supported by ForwardDiff. I would suggest using ExponentialUtilities.jl and it should be fine.

1 Like

Thanks!

Is there a more elegant way to do this?

function f(u, p)
uu = copy(u)
return exponential!(uu) - Matrix{Float64}(I(dim))
end


exponential!(copy(u)) - I?

Well, yeah, sure. I mean, there is no not-in-place (non-public?) function already?

Not yet. It would be a nice thing to add.

It turned out that only exponential! works with nonlinearsolve.jl, while expv doesnβt. Bug or feature?

solve(prob)
ERROR: MethodError: no method matching exponential!(::Matrix{ForwardDiff.Dual{ForwardDiff.Tag{β¦}, Float64, 3}}, ::ExpMethodHigham2005Base)

Closest candidates are:
exponential!(::StridedMatrix{T}, ::ExpMethodHigham2005Base, ::Any) where T<:Union{Float32, Float64, ComplexF64, ComplexF32}
@ ExponentialUtilities ~/.julia/packages/ExponentialUtilities/xLH9y/src/exp_baseexp.jl:23
exponential!(::StridedMatrix{T}, ::ExpMethodHigham2005Base) where T<:Union{Float32, Float64, ComplexF64, ComplexF32}
@ ExponentialUtilities ~/.julia/packages/ExponentialUtilities/xLH9y/src/exp_baseexp.jl:23
exponential!(::Any, ::ExpMethodNative, ::Any)
@ ExponentialUtilities ~/.julia/packages/ExponentialUtilities/xLH9y/src/exp.jl:70
...

Stacktrace:
[1] expv!(w::Vector{β¦}, t::Float64, Ks::KrylovSubspace{β¦}; cache::Nothing, expmethod::ExpMethodHigham2005Base)
@ ExponentialUtilities ~/.julia/packages/ExponentialUtilities/xLH9y/src/krylov_phiv.jl:95
[2] _expv_hb(t::Float64, A::Matrix{β¦}, b::Vector{β¦}; expmethod::ExpMethodHigham2005Base, cache::Nothing, kwargs_arnoldi::@Kwargs{})
@ ExponentialUtilities ~/.julia/packages/ExponentialUtilities/xLH9y/src/krylov_phiv.jl:51
[3] _expv_hb(t::Float64, A::Matrix{ForwardDiff.Dual{ForwardDiff.Tag{β¦}, Float64, 3}}, b::Vector{Float64})
@ ExponentialUtilities ~/.julia/packages/ExponentialUtilities/xLH9y/src/krylov_phiv.jl:45
[4] expv(t::Float64, A::Matrix{ForwardDiff.Dual{β¦}}, b::Vector{Float64}; mode::Symbol, kwargs::@Kwargs{})
@ ExponentialUtilities ~/.julia/packages/ExponentialUtilities/xLH9y/src/krylov_phiv.jl:38
[5] expv(t::Float64, A::Matrix{ForwardDiff.Dual{ForwardDiff.Tag{β¦}, Float64, 3}}, b::Vector{Float64})
@ ExponentialUtilities ~/.julia/packages/ExponentialUtilities/xLH9y/src/krylov_phiv.jl:35
[6] optgoal(v::Vector{ForwardDiff.Dual{ForwardDiff.Tag{β¦}, Float64, 3}}, p::Vector{Float64})
@ Main ./REPL[13]:2
[7] NonlinearFunction
@ ~/.julia/packages/SciMLBase/QSc1r/src/scimlfunctions.jl:2197 [inlined]
[8] JacobianWrapper
@ ~/.julia/packages/SciMLBase/QSc1r/src/function_wrappers.jl:97 [inlined]
[9] vector_mode_dual_eval!
@ ~/.julia/packages/ForwardDiff/PcZ48/src/apiutils.jl:24 [inlined]
[10] vector_mode_jacobian!(result::Matrix{β¦}, f::SciMLBase.JacobianWrapper{β¦}, x::Vector{β¦}, cfg::ForwardDiff.JacobianConfig{β¦})
@ ForwardDiff ~/.julia/packages/ForwardDiff/PcZ48/src/jacobian.jl:145
[11] jacobian!
@ ~/.julia/packages/ForwardDiff/PcZ48/src/jacobian.jl:58 [inlined]
[12] jacobian!
@ ~/.julia/packages/ForwardDiff/PcZ48/src/jacobian.jl:56 [inlined]
[13] sparse_jacobian!
@ ~/.julia/packages/SparseDiffTools/qxnHN/src/highlevel/forward_mode.jl:60 [inlined]
[14] JacobianCache
@ ~/.julia/packages/NonlinearSolve/5yLII/src/internal/jacobian.jl:136 [inlined]
[15] JacobianCache
@ ~/.julia/packages/NonlinearSolve/5yLII/src/internal/jacobian.jl:106 [inlined]
[16] __step!(cache::NonlinearSolve.GeneralizedFirstOrderAlgorithmCache{β¦}; recompute_jacobian::Nothing, kwargs::@Kwargs{})
@ NonlinearSolve ~/.julia/packages/NonlinearSolve/5yLII/src/core/generalized_first_order.jl:207
[17] __step!
@ ~/.julia/packages/NonlinearSolve/5yLII/src/core/generalized_first_order.jl:203 [inlined]
[18] #step!#210
@ ~/.julia/packages/NonlinearSolve/5yLII/src/core/generic.jl:55 [inlined]
[19] step!
@ ~/.julia/packages/NonlinearSolve/5yLII/src/core/generic.jl:50 [inlined]
[20] solve!(cache::NonlinearSolve.GeneralizedFirstOrderAlgorithmCache{β¦})
@ NonlinearSolve ~/.julia/packages/NonlinearSolve/5yLII/src/core/generic.jl:13
[21] __solve(::NonlinearProblem{β¦}, ::GeneralizedFirstOrderAlgorithm{β¦}; kwargs::@Kwargs{β¦})
@ NonlinearSolve ~/.julia/packages/NonlinearSolve/5yLII/src/core/generic.jl:4
[22] __solve
@ ~/.julia/packages/NonlinearSolve/5yLII/src/core/generic.jl:1 [inlined]
[23] macro expansion
@ ~/.julia/packages/NonlinearSolve/5yLII/src/default.jl:278 [inlined]
[24] __solve(::NonlinearProblem{β¦}, ::NonlinearSolvePolyAlgorithm{β¦}; alias_u0::Bool, verbose::Bool, kwargs::@Kwargs{β¦})
@ NonlinearSolve ~/.julia/packages/NonlinearSolve/5yLII/src/default.jl:248
[25] __solve
@ ~/.julia/packages/NonlinearSolve/5yLII/src/default.jl:248 [inlined]
[26] #__solve#329
@ ~/.julia/packages/NonlinearSolve/5yLII/src/default.jl:511 [inlined]
[27] __solve
@ ~/.julia/packages/NonlinearSolve/5yLII/src/default.jl:508 [inlined]
[28] #__solve#61
@ ~/.julia/packages/DiffEqBase/nc5nE/src/solve.jl:1371 [inlined]
[29] __solve
@ ~/.julia/packages/DiffEqBase/nc5nE/src/solve.jl:1364 [inlined]
[30] #solve_call#34
@ ~/.julia/packages/DiffEqBase/nc5nE/src/solve.jl:609 [inlined]
[31] solve_call
@ ~/.julia/packages/DiffEqBase/nc5nE/src/solve.jl:567 [inlined]
[32] #solve_up#42
@ ~/.julia/packages/DiffEqBase/nc5nE/src/solve.jl:1050 [inlined]
[33] solve_up
@ ~/.julia/packages/DiffEqBase/nc5nE/src/solve.jl:1044 [inlined]
[34] #solve#41
@ ~/.julia/packages/DiffEqBase/nc5nE/src/solve.jl:1038 [inlined]
[35] solve(::NonlinearProblem{β¦})
@ DiffEqBase ~/.julia/packages/DiffEqBase/nc5nE/src/solve.jl:1028
[36] top-level scope
@ REPL[17]:1
Some type information was truncated. Use show(err) to see complete types.


Unexpected. Worth an issue. Shouldnβt be too hard to solve, looks like itβs just a missing dispatch or a poorly chosen default.