Merge branch 'master' into u/termination_condition

Author: utkarsh530

.github/workflows/CI.yml

@@ -14,11 +14,13 @@ jobs:
   test:
     runs-on: ubuntu-latest
     strategy:
+      fail-fast: false
       matrix:
         group:
           - All
         version:
           - '1'
+          - '~1.10.0-0'
     steps:
       - uses: actions/checkout@v4
       - uses: julia-actions/setup-julia@v1

Project.toml

@@ -1,7 +1,7 @@
 name = "NonlinearSolve"
 uuid = "8913a72c-1f9b-4ce2-8d82-65094dcecaec"
 authors = ["SciML"]
-version = "2.3.0"
+version = "2.5.0"
 
 [deps]
 ADTypes = "47edcb42-4c32-4615-8424-f2b9edc5f35b"
@@ -44,13 +44,13 @@ FiniteDiff = "2"
 ForwardDiff = "0.10.3"
 LeastSquaresOptim = "0.8"
 LineSearches = "7"
-LinearSolve = "2"
+LinearSolve = "2.12"
 NonlinearProblemLibrary = "0.1"
 PrecompileTools = "1"
 RecursiveArrayTools = "2"
 Reexport = "0.2, 1"
 SciMLBase = "2.4"
-SimpleNonlinearSolve = "0.1.22"
+SimpleNonlinearSolve = "0.1.23"
 SparseDiffTools = "2.6"
 StaticArraysCore = "1.4"
 UnPack = "1.0"
@@ -65,6 +65,7 @@ ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 LeastSquaresOptim = "0fc2ff8b-aaa3-5acd-a817-1944a5e08891"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 LinearSolve = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae"
+NaNMath = "77ba4419-2d1f-58cd-9bb1-8ffee604a2e3"
 NonlinearProblemLibrary = "b7050fa9-e91f-4b37-bcee-a89a063da141"
 Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
@@ -77,4 +78,4 @@ Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f"
 DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e"
 
 [targets]
-test = ["Enzyme", "BenchmarkTools", "SafeTestsets", "Pkg", "Test", "ForwardDiff", "StaticArrays", "Symbolics", "LinearSolve", "Random", "LinearAlgebra", "Zygote", "SparseDiffTools", "NonlinearProblemLibrary", "LeastSquaresOptim", "FastLevenbergMarquardt", "DiffEqBase"]
+test = ["Enzyme", "BenchmarkTools", "SafeTestsets", "Pkg", "Test", "ForwardDiff", "StaticArrays", "Symbolics", "LinearSolve", "Random", "LinearAlgebra", "Zygote", "SparseDiffTools", "NonlinearProblemLibrary", "LeastSquaresOptim", "FastLevenbergMarquardt", "NaNMath", "DiffEqBase"]

docs/pages.jl

@@ -17,7 +17,8 @@ pages = ["index.md",
     "Solver Summaries and Recommendations" => Any["solvers/NonlinearSystemSolvers.md",
         "solvers/BracketingSolvers.md",
         "solvers/SteadyStateSolvers.md",
-        "solvers/NonlinearLeastSquaresSolvers.md"],
+        "solvers/NonlinearLeastSquaresSolvers.md",
+        "solvers/LineSearch.md"],
     "Detailed Solver APIs" => Any["api/nonlinearsolve.md",
         "api/simplenonlinearsolve.md",
         "api/minpack.md",

docs/src/api/nonlinearsolve.md

@@ -7,6 +7,10 @@ These are the native solvers of NonlinearSolve.jl.
 ```@docs
 NewtonRaphson
 TrustRegion
+PseudoTransient
+DFSane
+GeneralBroyden
+GeneralKlement
 ```
 
 ## Polyalgorithms

docs/src/solvers/LineSearch.md

@@ -0,0 +1,14 @@
+# [Line Search](@id linesearch)
+
+A convenience wrapper over `LineSearches.jl` and some native Line Search methods, powered
+internally with fast automatic differentiation.
+
+```@docs
+LineSearch
+```
+
+## Native Line Search Methods
+
+```@docs
+LiFukushimaLineSearch
+```

docs/src/solvers/NonlinearLeastSquaresSolvers.md

@@ -19,10 +19,8 @@ Solves the nonlinear least squares problem defined by `prob` using the algorithm
     handling of sparse matrices via colored automatic differentiation and preconditioned
     linear solvers. Designed for large-scale and numerically-difficult nonlinear least squares
     problems.
-  - `SimpleNewtonRaphson()`: Newton Raphson implementation that uses Linear Least Squares
-    solution at every step to compute the descent direction. **WARNING**: This method is not
-    a robust solver for nonlinear least squares problems. The computed delta step might not
-    be the correct descent direction!
+  - `SimpleNewtonRaphson()`: Simple Gauss Newton Implementation with `QRFactorization` to
+    solve a linear least squares problem at each step!
 
 ## Example usage
 

docs/src/solvers/NonlinearSystemSolvers.md

@@ -55,13 +55,26 @@ features, but have a bit of overhead on very small problems.
     improvements suggested in the [paper](https://arxiv.org/abs/1201.5885) "Improvements to
     the Levenberg-Marquardt algorithm for nonlinear least-squares minimization". Designed for
     large-scale and numerically-difficult nonlinear systems.
+  - `PseudoTransient()`: A pseudo-transient method which mixes the stability of Euler-type
+    stepping with the convergence speed of a Newton method. Good for highly unstable
+    systems.
   - `RobustMultiNewton()`: A polyalgorithm that mixes highly robust methods (line searches and
     trust regions) in order to be as robust as possible for difficult problems. If this method
     fails to converge, then one can be pretty certain that most (all?) other choices would
     likely fail.
-  - `FastShortcutNonlinearPolyalg`: The default method. A polyalgorithm that mixes fast methods
+  - `FastShortcutNonlinearPolyalg()`: The default method. A polyalgorithm that mixes fast methods
     with fallbacks to robust methods to allow for solving easy problems quickly without sacrificing
     robustnes on the hard problems.
+  - `GeneralBroyden()`: Generalization of Broyden's Quasi-Newton Method with Line Search and
+    Automatic Jacobian Resetting. This is a fast method but unstable when the condition number of
+    the Jacobian matrix is sufficiently large.
+  - `GeneralKlement()`: Generalization of Klement's Quasi-Newton Method with Line Search and
+    Automatic Jacobian Resetting. This is a fast method but unstable when the condition number of
+    the Jacobian matrix is sufficiently large.
+  - `LimitedMemoryBroyden()`: An advanced version of `LBroyden` which uses a limited memory
+    Broyden method. This is a fast method but unstable when the condition number of
+    the Jacobian matrix is sufficiently large. It is recommended to use `GeneralBroyden` or
+    `GeneralKlement` instead unless the memory usage is a concern.
 
 ### SimpleNonlinearSolve.jl
 

src/NonlinearSolve.jl

@@ -5,10 +5,11 @@ if isdefined(Base, :Experimental) && isdefined(Base.Experimental, Symbol("@max_m
 end
 
 using DiffEqBase, LinearAlgebra, LinearSolve, SparseArrays, SparseDiffTools
+import ArrayInterface: restructure
 import ForwardDiff
 
 import ADTypes: AbstractFiniteDifferencesMode
-import ArrayInterface: undefmatrix, matrix_colors, parameterless_type, ismutable
+import ArrayInterface: undefmatrix, matrix_colors, parameterless_type, ismutable, issingular
 import ConcreteStructs: @concrete
 import EnumX: @enumx
 import ForwardDiff: Dual
@@ -25,6 +26,8 @@ import UnPack: @unpack
 const AbstractSparseADType = Union{ADTypes.AbstractSparseFiniteDifferences,
     ADTypes.AbstractSparseForwardMode, ADTypes.AbstractSparseReverseMode}
 
+abstract type AbstractNonlinearSolveLineSearchAlgorithm end
+
 abstract type AbstractNonlinearSolveAlgorithm <: AbstractNonlinearAlgorithm end
 abstract type AbstractNewtonAlgorithm{CJ, AD} <: AbstractNonlinearSolveAlgorithm end
 
@@ -49,10 +52,13 @@ function SciMLBase.solve!(cache::AbstractNonlinearSolveCache)
         cache.stats.nsteps += 1
     end
 
-    if cache.stats.nsteps == cache.maxiters
-        cache.retcode = ReturnCode.MaxIters
-    else
-        cache.retcode = ReturnCode.Success
+    # The solver might have set a different `retcode`
+    if cache.retcode == ReturnCode.Default
+        if cache.stats.nsteps == cache.maxiters
+            cache.retcode = ReturnCode.MaxIters
+        else
+            cache.retcode = ReturnCode.Success
+        end
     end
 
     return SciMLBase.build_solution(cache.prob, cache.alg, cache.u, get_fu(cache);
@@ -67,18 +73,23 @@ include("trustRegion.jl")
 include("levenberg.jl")
 include("gaussnewton.jl")
 include("dfsane.jl")
+include("pseudotransient.jl")
+include("broyden.jl")
+include("klement.jl")
+include("lbroyden.jl")
 include("jacobian.jl")
 include("ad.jl")
 include("default.jl")
 
 import PrecompileTools
 
-@static if VERSION >= v"1.10"
+@static if VERSION ≥ v"1.10-"
     PrecompileTools.@compile_workload begin
         for T in (Float32, Float64)
             prob = NonlinearProblem{false}((u, p) -> u .* u .- p, T(0.1), T(2))
 
-            precompile_algs = (NewtonRaphson(), TrustRegion(), LevenbergMarquardt())
+            precompile_algs = (NewtonRaphson(), TrustRegion(), LevenbergMarquardt(),
+                PseudoTransient(), GeneralBroyden(), GeneralKlement(), nothing)
 
             for alg in precompile_algs
                 solve(prob, alg, abstol = T(1e-2))
@@ -95,10 +106,11 @@ end
 
 export RadiusUpdateSchemes
 
-export NewtonRaphson, TrustRegion, LevenbergMarquardt, DFSane, GaussNewton
+export NewtonRaphson, TrustRegion, LevenbergMarquardt, DFSane, GaussNewton, PseudoTransient,
+    GeneralBroyden, GeneralKlement, LimitedMemoryBroyden
 export LeastSquaresOptimJL, FastLevenbergMarquardtJL
 export RobustMultiNewton, FastShortcutNonlinearPolyalg
 
-export LineSearch
+export LineSearch, LiFukushimaLineSearch
 
 end # module

src/broyden.jl

@@ -0,0 +1,168 @@
+# Sadly `Broyden` is taken up by SimpleNonlinearSolve.jl
+"""
+    GeneralBroyden(; max_resets = 3, linesearch = LineSearch(), reset_tolerance = nothing)
+
+An implementation of `Broyden` with reseting and line search.
+
+## Arguments
+
+  - `max_resets`: the maximum number of resets to perform. Defaults to `3`.
+  - `reset_tolerance`: the tolerance for the reset check. Defaults to
+    `sqrt(eps(eltype(u)))`.
+  - `linesearch`: the line search algorithm to use. Defaults to [`LineSearch()`](@ref),
+    which means that no line search is performed. Algorithms from `LineSearches.jl` can be
+    used here directly, and they will be converted to the correct `LineSearch`. It is
+    recommended to use [LiFukushimaLineSearchCache](@ref) -- a derivative free linesearch
+    specifically designed for Broyden's method.
+"""
+@concrete struct GeneralBroyden <: AbstractNewtonAlgorithm{false, Nothing}
+    max_resets::Int
+    reset_tolerance
+    linesearch
+end
+
+function GeneralBroyden(; max_resets = 3, linesearch = LineSearch(),
+    reset_tolerance = nothing)
+    linesearch = linesearch isa LineSearch ? linesearch : LineSearch(; method = linesearch)
+    return GeneralBroyden(max_resets, reset_tolerance, linesearch)
+end
+
+@concrete mutable struct GeneralBroydenCache{iip} <: AbstractNonlinearSolveCache{iip}
+    f
+    alg
+    u
+    du
+    fu
+    fu2
+    dfu
+    p
+    J⁻¹
+    J⁻¹₂
+    J⁻¹df
+    force_stop::Bool
+    resets::Int
+    max_resets::Int
+    maxiters::Int
+    internalnorm
+    retcode::ReturnCode.T
+    abstol
+    reset_tolerance
+    reset_check
+    prob
+    stats::NLStats
+    lscache
+end
+
+get_fu(cache::GeneralBroydenCache) = cache.fu
+
+function SciMLBase.__init(prob::NonlinearProblem{uType, iip}, alg::GeneralBroyden, args...;
+    alias_u0 = false, maxiters = 1000, abstol = 1e-6, internalnorm = DEFAULT_NORM,
+    kwargs...) where {uType, iip}
+    @unpack f, u0, p = prob
+    u = alias_u0 ? u0 : deepcopy(u0)
+    fu = evaluate_f(prob, u)
+    J⁻¹ = __init_identity_jacobian(u, fu)
+    reset_tolerance = alg.reset_tolerance === nothing ? sqrt(eps(eltype(u))) :
+                      alg.reset_tolerance
+    reset_check = x -> abs(x) ≤ reset_tolerance
+    return GeneralBroydenCache{iip}(f, alg, u, _mutable_zero(u), fu, zero(fu),
+        zero(fu), p, J⁻¹, zero(_reshape(fu, 1, :)), _mutable_zero(u), false, 0,
+        alg.max_resets, maxiters, internalnorm, ReturnCode.Default, abstol, reset_tolerance,
+        reset_check, prob, NLStats(1, 0, 0, 0, 0),
+        init_linesearch_cache(alg.linesearch, f, u, p, fu, Val(iip)))
+end
+
+function perform_step!(cache::GeneralBroydenCache{true})
+    @unpack f, p, du, fu, fu2, dfu, u, J⁻¹, J⁻¹df, J⁻¹₂ = cache
+    T = eltype(u)
+
+    mul!(_vec(du), J⁻¹, -_vec(fu))
+    α = perform_linesearch!(cache.lscache, u, du)
+    _axpy!(α, du, u)
+    f(fu2, u, p)
+
+    cache.internalnorm(fu2) < cache.abstol && (cache.force_stop = true)
+    cache.stats.nf += 1
+
+    cache.force_stop && return nothing
+
+    # Update the inverse jacobian
+    dfu .= fu2 .- fu
+
+    if all(cache.reset_check, du) || all(cache.reset_check, dfu)
+        if cache.resets ≥ cache.max_resets
+            cache.retcode = ReturnCode.Unstable
+            cache.force_stop = true
+            return nothing
+        end
+        fill!(J⁻¹, 0)
+        J⁻¹[diagind(J⁻¹)] .= T(1)
+        cache.resets += 1
+    else
+        mul!(_vec(J⁻¹df), J⁻¹, _vec(dfu))
+        mul!(J⁻¹₂, _vec(du)', J⁻¹)
+        denom = dot(du, J⁻¹df)
+        du .= (du .- J⁻¹df) ./ ifelse(iszero(denom), T(1e-5), denom)
+        mul!(J⁻¹, _vec(du), J⁻¹₂, 1, 1)
+    end
+    fu .= fu2
+
+    return nothing
+end
+
+function perform_step!(cache::GeneralBroydenCache{false})
+    @unpack f, p = cache
+    T = eltype(cache.u)
+
+    cache.du = _restructure(cache.du, cache.J⁻¹ * -_vec(cache.fu))
+    α = perform_linesearch!(cache.lscache, cache.u, cache.du)
+    cache.u = cache.u .+ α * cache.du
+    cache.fu2 = f(cache.u, p)
+
+    cache.internalnorm(cache.fu2) < cache.abstol && (cache.force_stop = true)
+    cache.stats.nf += 1
+
+    cache.force_stop && return nothing
+
+    # Update the inverse jacobian
+    cache.dfu = cache.fu2 .- cache.fu
+    if all(cache.reset_check, cache.du) || all(cache.reset_check, cache.dfu)
+        if cache.resets ≥ cache.max_resets
+            cache.retcode = ReturnCode.Unstable
+            cache.force_stop = true
+            return nothing
+        end
+        cache.J⁻¹ = __init_identity_jacobian(cache.u, cache.fu)
+        cache.resets += 1
+    else
+        cache.J⁻¹df = _restructure(cache.J⁻¹df, cache.J⁻¹ * _vec(cache.dfu))
+        cache.J⁻¹₂ = _vec(cache.du)' * cache.J⁻¹
+        denom = dot(cache.du, cache.J⁻¹df)
+        cache.du = (cache.du .- cache.J⁻¹df) ./ ifelse(iszero(denom), T(1e-5), denom)
+        cache.J⁻¹ = cache.J⁻¹ .+ _vec(cache.du) * cache.J⁻¹₂
+    end
+    cache.fu = cache.fu2
+
+    return nothing
+end
+
+function SciMLBase.reinit!(cache::GeneralBroydenCache{iip}, u0 = cache.u; p = cache.p,
+    abstol = cache.abstol, maxiters = cache.maxiters) where {iip}
+    cache.p = p
+    if iip
+        recursivecopy!(cache.u, u0)
+        cache.f(cache.fu, cache.u, p)
+    else
+        # don't have alias_u0 but cache.u is never mutated for OOP problems so it doesn't matter
+        cache.u = u0
+        cache.fu = cache.f(cache.u, p)
+    end
+    cache.abstol = abstol
+    cache.maxiters = maxiters
+    cache.stats.nf = 1
+    cache.stats.nsteps = 1
+    cache.resets = 0
+    cache.force_stop = false
+    cache.retcode = ReturnCode.Default
+    return cache
+end

src/default.jl

@@ -159,10 +159,8 @@ end
         ]
     else
         [
-            # FIXME: Broyden and Klement are type unstable
-            #        (upstream SimpleNonlinearSolve.jl issue)
-            !iip ? :(Klement()) : nothing, # Klement not yet implemented for IIP
-            !iip ? :(Broyden()) : nothing, # Broyden not yet implemented for IIP
+            :(GeneralKlement()),
+            :(GeneralBroyden()),
             :(NewtonRaphson(; linsolve, precs, adkwargs...)),
             :(NewtonRaphson(; linsolve, precs, linesearch = BackTracking(), adkwargs...)),
             :(TrustRegion(; linsolve, precs, adkwargs...)),