Merge pull request #257 from avik-pal/ap/more_alg

Porting over some low cost Algorithms

Author: ChrisRackauckas

.github/workflows/CI.yml

@@ -16,7 +16,8 @@ jobs:
     strategy:
       matrix:
         group:
-          - All
+          - Core
+          - 23TestProblems
         version:
           - '1'
     steps:

docs/pages.jl

@@ -2,8 +2,7 @@
 
 pages = ["index.md",
     "Getting Started with Nonlinear Rootfinding in Julia" => "tutorials/getting_started.md",
-    "Tutorials" => Any[
-        "Code Optimization for Small Nonlinear Systems" => "tutorials/code_optimization.md",
+    "Tutorials" => Any["Code Optimization for Small Nonlinear Systems" => "tutorials/code_optimization.md",
         "Handling Large Ill-Conditioned and Sparse Systems" => "tutorials/large_systems.md",
         "Symbolic System Definition and Acceleration via ModelingToolkit" => "tutorials/modelingtoolkit.md",
         "tutorials/small_compile.md",
@@ -18,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

@@ -8,6 +8,9 @@ These are the native solvers of NonlinearSolve.jl.
 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/NonlinearSystemSolvers.md

@@ -65,6 +65,12 @@ features, but have a bit of overhead on very small problems.
   - `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.
 
 ### SimpleNonlinearSolve.jl
 

docs/src/tutorials/code_optimization.md

@@ -34,7 +34,7 @@ Take for example a prototypical small nonlinear solver code in its out-of-place
 ```@example small_opt
 using NonlinearSolve
 
-function f(u, p) 
+function f(u, p)
     u .* u .- p
 end
 u0 = [1.0, 1.0]
@@ -54,7 +54,7 @@ using BenchmarkTools
 Note that this way of writing the function is a shorthand for:
 
 ```@example small_opt
-function f(u, p) 
+function f(u, p)
     [u[1] * u[1] - p, u[2] * u[2] - p]
 end
 ```
@@ -69,25 +69,25 @@ In order to avoid this issue, we can use a non-allocating "in-place" approach. W
 this looks like:
 
 ```@example small_opt
-function f(du, u, p) 
+function f(du, u, p)
     du[1] = u[1] * u[1] - p
     du[2] = u[2] * u[2] - p
     nothing
 end
 
 prob = NonlinearProblem(f, u0, p)
-@btime  sol = solve(prob, NewtonRaphson())
+@btime sol = solve(prob, NewtonRaphson())
 ```
 
 Notice how much faster this already runs! We can make this code even simpler by using
 the `.=` in-place broadcasting.
 
 ```@example small_opt
-function f(du, u, p) 
+function f(du, u, p)
     du .= u .* u .- p
 end
 
-@btime  sol = solve(prob, NewtonRaphson())
+@btime sol = solve(prob, NewtonRaphson())
 ```
 
 ## Further Optimizations for Small Nonlinear Solves with Static Arrays and SimpleNonlinearSolve
@@ -140,7 +140,7 @@ want to use the out-of-place allocating form, but this time we want to output
 a static array. Doing it with broadcasting looks like:
 
 ```@example small_opt
-function f_SA(u, p) 
+function f_SA(u, p)
     u .* u .- p
 end
 u0 = SA[1.0, 1.0]
@@ -153,7 +153,7 @@ Note that only change here is that `u0` is made into a StaticArray! If we needed
 for a more complex nonlinear case, then we'd simply do the following:
 
 ```@example small_opt
-function f_SA(u, p) 
+function f_SA(u, p)
     SA[u[1] * u[1] - p, u[2] * u[2] - p]
 end
 
@@ -170,4 +170,4 @@ which are designed for these small-scale static examples. Let's now use `SimpleN
 @btime solve(prob, SimpleNewtonRaphson())
 ```
 
-And there we go, around 100ns from our starting point of almost 6μs!
+And there we go, around 100ns from our starting point of almost 6μs!

docs/src/tutorials/large_systems.md

@@ -1,8 +1,8 @@
 # [Efficiently Solving Large Sparse Ill-Conditioned Nonlinear Systems in Julia](@id large_systems)
 
-This tutorial is for getting into the extra features of using NonlinearSolve.jl. Solving ill-conditioned nonlinear systems 
-requires specializing the linear solver on properties of the Jacobian in order to cut down on the ``\mathcal{O}(n^3)`` 
-linear solve and the ``\mathcal{O}(n^2)`` back-solves. This tutorial is designed to explain the advanced usage of 
+This tutorial is for getting into the extra features of using NonlinearSolve.jl. Solving ill-conditioned nonlinear systems
+requires specializing the linear solver on properties of the Jacobian in order to cut down on the ``\mathcal{O}(n^3)``
+linear solve and the ``\mathcal{O}(n^2)`` back-solves. This tutorial is designed to explain the advanced usage of
 NonlinearSolve.jl by solving the steady state stiff Brusselator partial differential equation (BRUSS) using NonlinearSolve.jl.
 
 ## Definition of the Brusselator Equation
@@ -182,8 +182,8 @@ nothing # hide
 
 Notice that this acceleration does not require the definition of a sparsity
 pattern, and can thus be an easier way to scale for large problems. For more
-information on linear solver choices, see the 
-[linear solver documentation](https://docs.sciml.ai/DiffEqDocs/stable/features/linear_nonlinear/#linear_nonlinear). 
+information on linear solver choices, see the
+[linear solver documentation](https://docs.sciml.ai/DiffEqDocs/stable/features/linear_nonlinear/#linear_nonlinear).
 `linsolve` choices are any valid [LinearSolve.jl](https://linearsolve.sciml.ai/dev/) solver.
 
 !!! note

docs/src/tutorials/modelingtoolkit.md

@@ -120,4 +120,4 @@ sol[u5]
 
 If you're interested in building models in a component or block based form, such as seen in systems like Simulink or Modelica,
 take a deeper look at [ModelingToolkit.jl's documentation](https://docs.sciml.ai/ModelingToolkit/stable/) which goes into
-detail on such features.
+detail on such features.

docs/src/tutorials/small_compile.md

@@ -19,18 +19,18 @@ sol = solve(prob, SimpleNewtonRaphson())
 
 However, there are a few downsides to SimpleNonlinearSolve's `SimpleX` style algorithms to note:
 
-1. SimpleNonlinearSolve.jl's methods are not hooked into the LinearSolve.jl system, and thus do not have
-   the ability to specify linear solvers, use sparse matrices, preconditioners, and all of the other features
-   which are required to scale for very large systems of equations.
-2. SimpleNonlinearSolve.jl's methods have less robust error handling and termination conditions, and thus
-   these methods are missing some flexibility and give worse hints for debugging.
-3. SimpleNonlinearSolve.jl's methods are focused on out-of-place support. There is some in-place support,
-   but it's designed for simple smaller systems and so some optimizations are missing.
+ 1. SimpleNonlinearSolve.jl's methods are not hooked into the LinearSolve.jl system, and thus do not have
+    the ability to specify linear solvers, use sparse matrices, preconditioners, and all of the other features
+    which are required to scale for very large systems of equations.
+ 2. SimpleNonlinearSolve.jl's methods have less robust error handling and termination conditions, and thus
+    these methods are missing some flexibility and give worse hints for debugging.
+ 3. SimpleNonlinearSolve.jl's methods are focused on out-of-place support. There is some in-place support,
+    but it's designed for simple smaller systems and so some optimizations are missing.
 
 However, the major upsides of SimpleNonlinearSolve.jl are:
 
-1. The methods are optimized and non-allocating on StaticArrays
-2. The methods are minimal in compilation
+ 1. The methods are optimized and non-allocating on StaticArrays
+ 2. The methods are minimal in compilation
 
 As such, you can use the code as shown above to have very low startup with good methods, but for more scaling and debuggability
 we recommend the full NonlinearSolve.jl. But that said,
@@ -51,4 +51,4 @@ is not only sufficient but optimal.
 
 Julia has tools for building small binaries via static compilation with [StaticCompiler.jl](https://github.com/tshort/StaticCompiler.jl).
 However, these tools are currently limited to type-stable non-allocating functions. That said, SimpleNonlinearSolve.jl's solvers are
-precisely the subset of NonlinearSolve.jl which are compatible with static compilation.
+precisely the subset of NonlinearSolve.jl which are compatible with static compilation.

src/NonlinearSolve.jl

@@ -9,7 +9,7 @@ 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
@@ -26,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
 
@@ -50,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);
@@ -69,18 +74,22 @@ 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))
@@ -97,10 +106,11 @@ end
 
 export RadiusUpdateSchemes
 
-export NewtonRaphson, TrustRegion, LevenbergMarquardt, DFSane, GaussNewton, PseudoTransient
+export NewtonRaphson, TrustRegion, LevenbergMarquardt, DFSane, GaussNewton, PseudoTransient,
+    GeneralBroyden, GeneralKlement, LimitedMemoryBroyden
 export LeastSquaresOptimJL, FastLevenbergMarquardtJL
 export RobustMultiNewton, FastShortcutNonlinearPolyalg
 
-export LineSearch
+export LineSearch, LiFukushimaLineSearch
 
 end # module