search

SciML / NonlinearSolve.jl

High-performance and differentiation-enabled nonlinear solvers (Newton methods), bracketed rootfinding (bisection, Falsi), with sparsity and Newton-Krylov support.
https://docs.sciml.ai/NonlinearSolve/stable/
MIT License
235 stars 41 forks source link

LinearAlgebra.SingularException(2) #153

Closed paulflang closed 1 year ago

paulflang commented 1 year ago

We are recently observing a LinearAlgebra.SingularException(2) with the following MWE:

using ModelingToolkit
using SciMLBase: SteadyStateProblem, NonlinearProblem, solve
using DifferentialEquations

iv = only(@variables(t))
sts = @variables s1(t) = 2.0 s1s2(t) = 2.0 s2(t) = 2.0
ps = @parameters k1 = 1.0 c1 = 2.0 [bounds=(0,2), tunable=true] Δt = 2.5
eqs = [
    (Differential(t))(s1) ~ -0.25 * c1 * k1 * s1 * s2
    (Differential(t))(s1s2) ~ 0.25 * c1 * k1 * s1 * s2
    (Differential(t))(s2) ~ -0.25 * c1 * k1 * s1 * s2
]
model = ODESystem(eqs, iv, sts, ps; name = :reactionsystem)

ss_prob = SteadyStateProblem(model, [], [])
prob = NonlinearProblem(ss_prob)
solve(prob)

We think this is due to an update in NonlinearSolve. (cc @SebastianM-C )

SebastianM-C commented 1 year ago

@ChrisRackauckas I tried to look into what causes the issue with the following:

using ModelingToolkit
using SciMLBase
using SteadyStateDiffEq
using NonlinearSolve
using LinearAlgebra

iv = only(@variables(t))
sts = @variables s1(t) = 2.0 s1s2(t) = 2.0 s2(t) = 2.0
ps = @parameters k1 = 1.0 c1 = 2.0 [bounds = (0, 2), tunable = true] Δt = 2.5
eqs = [
       (Differential(t))(s1) ~ -0.25 * c1 * k1 * s1 * s2
       (Differential(t))(s1s2) ~ 0.25 * c1 * k1 * s1 * s2
       (Differential(t))(s2) ~ -0.25 * c1 * k1 * s1 * s2
]
model = ODESystem(eqs, iv, sts, ps; name=:reactionsystem)

ss_prob = SteadyStateProblem(model, [], [])
prob = NonlinearProblem(ss_prob)

cache = init(prob, NewtonRaphson());
det(cache.J)

And it looks like the NewtonRaphson cache is initialized with a singular cache.J. Is this correct? Also, I see that the determinant is sometimes 0 sometimes NaN, not sure if this is relevant, but I assume it's because we initialize the cache with an undef matrix and what's in there is more or less random.

The error that comes from the above mentioned comment is

ERROR: SingularException(2)
Stacktrace:
  [1] checknonsingular
    @ C:\Users\sebastian\.julia\juliaup\julia-1.9.0-rc2+0.x64.w64.mingw32\share\julia\stdlib\v1.9\LinearAlgebra\src\factorization.jl:19 [inlined]
  [2] generic_lufact!(A::Matrix{Float64}, pivot::RowMaximum; check::Bool)
    @ LinearAlgebra C:\Users\sebastian\.julia\juliaup\julia-1.9.0-rc2+0.x64.w64.mingw32\share\julia\stdlib\v1.9\LinearAlgebra\src\lu.jl:190
  [3] generic_lufact!
    @ C:\Users\sebastian\.julia\juliaup\julia-1.9.0-rc2+0.x64.w64.mingw32\share\julia\stdlib\v1.9\LinearAlgebra\src\lu.jl:133 [inlined]
  [4] do_factorization
    @ C:\Users\sebastian\.julia\packages\LinearSolve\LD2dF\src\factorization.jl:65 [inlined]
  [5] #solve#7
    @ C:\Users\sebastian\.julia\packages\LinearSolve\LD2dF\src\factorization.jl:11 [inlined]
  [6] solve
    @ C:\Users\sebastian\.julia\packages\LinearSolve\LD2dF\src\factorization.jl:9 [inlined]
  [7] #solve#6
    @ C:\Users\sebastian\.julia\packages\LinearSolve\LD2dF\src\common.jl:224 [inlined]
  [8] solve
    @ C:\Users\sebastian\.julia\packages\LinearSolve\LD2dF\src\common.jl:223 [inlined]
  [9] #dolinsolve#4
    @ C:\Users\sebastian\.julia\packages\NonlinearSolve\dQWiD\src\utils.jl:95 [inlined]
 [10] dolinsolve
    @ C:\Users\sebastian\.julia\packages\NonlinearSolve\dQWiD\src\utils.jl:69 [inlined]
 [11] perform_step!(cache::NonlinearSolve.NewtonRaphsonCache{true, NonlinearFunction{true, SciMLBase.AutoSpecialize, SciMLBase.var"#41#50"{ODEFunction{true, SciMLBase.AutoSpecialize, ModelingToolkit.var"#f#520"{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x3ef33a87, 0xc20c61d7, 0x66b4e26f, 0xb00ac2e7, 0x36aa7f8b)}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0xd9485a69, 0x1115f17e, 0x9689970c, 0xaeb9e622, 0xb605bfef)}}, UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Symbol, Vector{Symbol}, ModelingToolkit.var"#generated_observed#525"{ODESystem, Dict{Any, Any}}, Nothing, ODESystem}}, UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Vector{Symbol}, ModelingToolkit.var"#generated_observed#525"{ODESystem, Dict{Any, Any}}, Nothing, Nothing}, NewtonRaphson{0, true, Val{:forward}, Nothing, typeof(NonlinearSolve.DEFAULT_PRECS), true, nothing}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, typeof(NonlinearSolve.DEFAULT_NORM), Float64, NonlinearProblem{Vector{Float64}, true, Vector{Float64}, NonlinearFunction{true, SciMLBase.AutoSpecialize, SciMLBase.var"#41#50"{ODEFunction{true, SciMLBase.AutoSpecialize, ModelingToolkit.var"#f#520"{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x3ef33a87, 0xc20c61d7, 0x66b4e26f, 0xb00ac2e7, 0x36aa7f8b)}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0xd9485a69, 0x1115f17e, 0x9689970c, 0xaeb9e622, 0xb605bfef)}}, UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Symbol, Vector{Symbol}, ModelingToolkit.var"#generated_observed#525"{ODESystem, Dict{Any, Any}}, Nothing, ODESystem}}, UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Vector{Symbol}, ModelingToolkit.var"#generated_observed#525"{ODESystem, Dict{Any, Any}}, Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardNonlinearProblem}, NonlinearSolve.JacobianWrapper{NonlinearFunction{true, SciMLBase.AutoSpecialize, SciMLBase.var"#41#50"{ODEFunction{true, SciMLBase.AutoSpecialize, ModelingToolkit.var"#f#520"{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x3ef33a87, 0xc20c61d7, 0x66b4e26f, 0xb00ac2e7, 0x36aa7f8b)}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0xd9485a69, 0x1115f17e, 0x9689970c, 0xaeb9e622, 0xb605bfef)}}, UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Symbol, Vector{Symbol}, ModelingToolkit.var"#generated_observed#525"{ODESystem, Dict{Any, Any}}, Nothing, ODESystem}}, UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Vector{Symbol}, ModelingToolkit.var"#generated_observed#525"{ODESystem, Dict{Any, Any}}, Nothing, Nothing}, Vector{Float64}}, LinearSolve.LinearCache{Matrix{Float64}, Vector{Float64}, Vector{Float64}, SciMLBase.NullParameters, LinearSolve.GenericLUFactorization{RowMaximum}, LU{Float64, Matrix{Float64}, Vector{Int64}}, LinearSolve.InvPreconditioner{Diagonal{Float64, Vector{Float64}}}, Diagonal{Float64, Vector{Float64}}, Float64, true, LinearSolve.OperatorCondition.IllConditioned}, Matrix{Float64}, SparseDiffTools.ForwardColorJacCache{Vector{ForwardDiff.Dual{ForwardDiff.Tag{NonlinearSolve.NonlinearSolveTag, Float64}, Float64, 3}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{NonlinearSolve.NonlinearSolveTag, Float64}, Float64, 3}}, Vector{Float64}, Vector{Vector{Tuple{Float64, Float64, Float64}}}, UnitRange{Int64}, Nothing}})
    @ NonlinearSolve C:\Users\sebastian\.julia\packages\NonlinearSolve\dQWiD\src\raphson.jl:163
 [12] solve!(cache::NonlinearSolve.NewtonRaphsonCache{true, NonlinearFunction{true, SciMLBase.AutoSpecialize, SciMLBase.var"#41#50"{ODEFunction{true, SciMLBase.AutoSpecialize, ModelingToolkit.var"#f#520"{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x3ef33a87, 0xc20c61d7, 0x66b4e26f, 0xb00ac2e7, 0x36aa7f8b)}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0xd9485a69, 0x1115f17e, 0x9689970c, 0xaeb9e622, 0xb605bfef)}}, UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Symbol, Vector{Symbol}, ModelingToolkit.var"#generated_observed#525"{ODESystem, Dict{Any, Any}}, Nothing, ODESystem}}, UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Vector{Symbol}, ModelingToolkit.var"#generated_observed#525"{ODESystem, Dict{Any, Any}}, Nothing, Nothing}, NewtonRaphson{0, true, Val{:forward}, Nothing, typeof(NonlinearSolve.DEFAULT_PRECS), true, nothing}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, typeof(NonlinearSolve.DEFAULT_NORM), Float64, NonlinearProblem{Vector{Float64}, true, Vector{Float64}, NonlinearFunction{true, SciMLBase.AutoSpecialize, SciMLBase.var"#41#50"{ODEFunction{true, SciMLBase.AutoSpecialize, ModelingToolkit.var"#f#520"{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x3ef33a87, 0xc20c61d7, 0x66b4e26f, 0xb00ac2e7, 0x36aa7f8b)}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0xd9485a69, 0x1115f17e, 0x9689970c, 0xaeb9e622, 0xb605bfef)}}, UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Symbol, Vector{Symbol}, ModelingToolkit.var"#generated_observed#525"{ODESystem, Dict{Any, Any}}, Nothing, ODESystem}}, UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Vector{Symbol}, ModelingToolkit.var"#generated_observed#525"{ODESystem, Dict{Any, Any}}, Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardNonlinearProblem}, NonlinearSolve.JacobianWrapper{NonlinearFunction{true, SciMLBase.AutoSpecialize, SciMLBase.var"#41#50"{ODEFunction{true, SciMLBase.AutoSpecialize, ModelingToolkit.var"#f#520"{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x3ef33a87, 0xc20c61d7, 0x66b4e26f, 0xb00ac2e7, 0x36aa7f8b)}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0xd9485a69, 0x1115f17e, 0x9689970c, 0xaeb9e622, 0xb605bfef)}}, UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Symbol, Vector{Symbol}, ModelingToolkit.var"#generated_observed#525"{ODESystem, Dict{Any, Any}}, Nothing, ODESystem}}, UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Vector{Symbol}, ModelingToolkit.var"#generated_observed#525"{ODESystem, Dict{Any, Any}}, Nothing, Nothing}, Vector{Float64}}, LinearSolve.LinearCache{Matrix{Float64}, Vector{Float64}, Vector{Float64}, SciMLBase.NullParameters, LinearSolve.GenericLUFactorization{RowMaximum}, LU{Float64, Matrix{Float64}, Vector{Int64}}, LinearSolve.InvPreconditioner{Diagonal{Float64, Vector{Float64}}}, Diagonal{Float64, Vector{Float64}}, Float64, true, LinearSolve.OperatorCondition.IllConditioned}, Matrix{Float64}, SparseDiffTools.ForwardColorJacCache{Vector{ForwardDiff.Dual{ForwardDiff.Tag{NonlinearSolve.NonlinearSolveTag, Float64}, Float64, 3}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{NonlinearSolve.NonlinearSolveTag, Float64}, Float64, 3}}, Vector{Float64}, Vector{Vector{Tuple{Float64, Float64, Float64}}}, UnitRange{Int64}, Nothing}})
    @ NonlinearSolve C:\Users\sebastian\.julia\packages\NonlinearSolve\dQWiD\src\raphson.jl:188
ChrisRackauckas commented 1 year ago

That test there makes no sense. Yes, the Jacobian before it's evaluated doesn't have values in it. But when it does get values in it, it's:

J = [-1.0 -0.0 -1.0; 1.0 0.0 1.0; -1.0 -0.0 -1.0]

which is singular. The Jacobian is singular though: and that's clearly matching the analytical solution of what the starting Jacobian is. I'm not sure what you're looking for though as the solution to those equations is clearly not unique because s1s2 is not used in any of the equations. It's a nonlinear solve for 3 equations in 2 variables which is why its Jacobian is singular.

This problem is well-defined as a steady state of an ODE though:

using DifferentialEquations, ModelingToolkit, SteadyStateDiffEq

iv = only(@variables(t))
sts = @variables s1(t) = 2.0 s1s2(t) = 2.0 s2(t) = 2.0
ps = @parameters k1 = 1.0 c1 = 2.0 [bounds=(0,2), tunable=true] Δt = 2.5
eqs = [
    (Differential(t))(s1) ~ -0.25 * c1 * k1 * s1 * s2
    (Differential(t))(s1s2) ~ 0.25 * c1 * k1 * s1 * s2
    (Differential(t))(s2) ~ -0.25 * c1 * k1 * s1 * s2
]
model = ODESystem(eqs, iv, sts, ps; name = :reactionsystem)

ss_prob = SteadyStateProblem(model, [], [])
solve(ss_prob, DifferentialEquations.SteadyStateDiffEq.DynamicSS(Tsit5()))

That works fine.

SebastianM-C commented 1 year ago

Ah, I see, thanks for clearing this up. In this case would it make sense for DifferentialEquations.jl to choose DynamicSS instead of NewtonRaphson if we have a singular Jacobian or should I just do that check on my own?