Fix return codes for nonlinear least squares

Fixes #459. The crux of the issue is that `f(x) = residual` only applies in the NonlinearProblem and SteadyStateProblem cases. When `f(x)` is a nonlinear least squares problem, finding a local minima is a solution, not a failure of the algorithm. Thus this reclassifies Stalled in NLLSQ to StalledSuccess, which makes it a successful return.

Algorithms which require the NonlinearLeastSquares solution to have `||resid|| < tol` thus need to be careful with the return handling, as is done in the PR that introduces this return code SciML/SciMLBase.jl#1016. However, that's a fairly odd case because it's feasibility checking, while the normal use case for NLLSQ is for optimization, and in an optimization case there's no reason to believe you should always have a solution close to zero.

Author: ChrisRackauckas

lib/NonlinearSolveBase/Project.toml

@@ -1,7 +1,7 @@
 name = "NonlinearSolveBase"
 uuid = "be0214bd-f91f-a760-ac4e-3421ce2b2da0"
 authors = ["Avik Pal <avikpal@mit.edu> and contributors"]
-version = "1.7.0"
+version = "1.8.0"
 
 [deps]
 ADTypes = "47edcb42-4c32-4615-8424-f2b9edc5f35b"
@@ -68,7 +68,7 @@ MaybeInplace = "0.1.4"
 Preferences = "1.4"
 Printf = "1.10"
 RecursiveArrayTools = "3"
-SciMLBase = "2.69"
+SciMLBase = "2.89"
 SciMLJacobianOperators = "0.1.1"
 SciMLOperators = "0.3.13, 0.4"
 SparseArrays = "1.10"

lib/NonlinearSolveBase/src/polyalg.jl

@@ -160,7 +160,10 @@ end
                     InternalAPI.step!($(cache_syms[i]), args...; kwargs...)
                     $(cache_syms[i]).nsteps += 1
                     if !NonlinearSolveBase.not_terminated($(cache_syms[i]))
-                        if SciMLBase.successful_retcode($(cache_syms[i]).retcode)
+                        # If a NonlinearLeastSquaresProblem StalledSuccess, try the next
+                        # solver to see if you get a lower residual
+                        if SciMLBase.successful_retcode($(cache_syms[i]).retcode) && 
+                            $(cache_syms[i]).retcode != ReturnCode.StalledSuccess
                             cache.best = $(i)
                             cache.force_stop = true
                             cache.retcode = $(cache_syms[i]).retcode

lib/NonlinearSolveBase/src/termination_conditions.jl

@@ -21,6 +21,7 @@ const AbsNormModes = Union{
     step_norm_trace
     max_stalled_steps
     u_diff_cache::uType
+    leastsq::Bool
 end
 
 get_abstol(cache::NonlinearTerminationModeCache) = cache.abstol
@@ -36,7 +37,7 @@ function update_u!!(cache::NonlinearTerminationModeCache, u)
 end
 
 function CommonSolve.init(
-        ::AbstractNonlinearProblem, mode::AbstractNonlinearTerminationMode, du, u,
+        prob::AbstractNonlinearProblem, mode::AbstractNonlinearTerminationMode, du, u,
         saved_value_prototype...; abstol = nothing, reltol = nothing, kwargs...
 )
     T = promote_type(eltype(du), eltype(u))
@@ -80,10 +81,12 @@ function CommonSolve.init(
 
     length(saved_value_prototype) == 0 && (saved_value_prototype = nothing)
 
+    leastsq = typeof(prob) <: NonlinearLeastSquaresProblem
+
     return NonlinearTerminationModeCache(
         u_unaliased, ReturnCode.Default, abstol, reltol, best_value, mode,
         initial_objective, objectives_trace, 0, saved_value_prototype,
-        u0_norm, step_norm_trace, max_stalled_steps, u_diff_cache
+        u0_norm, step_norm_trace, max_stalled_steps, u_diff_cache, leastsq
     )
 end
 
@@ -146,6 +149,7 @@ end
 function (cache::NonlinearTerminationModeCache)(
         mode::AbstractSafeNonlinearTerminationMode, du, u, uprev, abstol, reltol, args...
 )
+
     if mode isa AbsNormSafeTerminationMode || mode isa AbsNormSafeBestTerminationMode
         objective = Utils.apply_norm(mode.internalnorm, du)
         criteria = abstol
@@ -177,7 +181,7 @@ function (cache::NonlinearTerminationModeCache)(
     end
 
     # Main Termination Criteria
-    if objective ≤ criteria
+    if !cache.leastsq && objective ≤ criteria
         cache.retcode = ReturnCode.Success
         return true
     end
@@ -195,7 +199,13 @@ function (cache::NonlinearTerminationModeCache)(
             min_obj, max_obj = extrema(cache.objectives_trace)
         end
         if min_obj < mode.min_max_factor * max_obj
-            cache.retcode = ReturnCode.Stalled
+            if cache.leastsq
+                # If least squares, found a local minima thus success
+                cache.retcode = ReturnCode.StalledSuccess
+            else
+                # Not a success if f(x)>0 and residual too high
+                cache.retcode = ReturnCode.Stalled
+            end
             return true
         end
     end
@@ -209,7 +219,7 @@ function (cache::NonlinearTerminationModeCache)(
         end
         du_norm = L2_NORM(cache.u_diff_cache)
         cache.step_norm_trace[mod1(cache.nsteps, length(cache.step_norm_trace))] = du_norm
-        if cache.nsteps > mode.max_stalled_steps
+        if cache.nsteps > mode.max_stalled_steps || iszero(du_norm)
             max_step_norm = maximum(cache.step_norm_trace)
             if mode isa AbsNormSafeTerminationMode ||
                mode isa AbsNormSafeBestTerminationMode
@@ -218,7 +228,11 @@ function (cache::NonlinearTerminationModeCache)(
                 stalled_step = max_step_norm ≤ reltol * (max_step_norm + cache.u0_norm)
             end
             if stalled_step
-                cache.retcode = ReturnCode.Stalled
+                if cache.leastsq
+                    cache.retcode = ReturnCode.StalledSuccess
+                else
+                    cache.retcode = ReturnCode.Stalled
+                end
                 return true
             end
         end

test/core_tests.jl

@@ -427,3 +427,13 @@ end
 
     @test !(solve(nlprob, NewtonRaphson()).alg.autodiff isa AutoPolyesterForwardDiff)
 end
+
+@testitem "NonlinearLeastSquares ReturnCode" tags=[:core] begin
+    f(u,p) = [1.0]
+    nlf = NonlinearFunction(f; resid_prototype=zeros(1))
+    prob = NonlinearLeastSquaresProblem(nlf, [1.0])
+    sol = solve(prob)
+    @test SciMLBase.successful_retcode(sol)
+    @test sol.retcode == ReturnCode.StalledSuccess
+    @test sol.stats.nf == 3
+end