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.