Merge pull request #105 from avik-pal/ap/minor_fixes

Add compat entries and change default norm

Author: ChrisRackauckas

Project.toml

@@ -18,11 +18,14 @@ SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
 StaticArraysCore = "1e83bf80-4336-4d27-bf5d-d5a4f845583c"
 
 [compat]
+ADTypes = "0.2"
 ArrayInterface = "7"
+ConcreteStructs = "0.2"
 DiffEqBase = "6.126"
 FiniteDiff = "2"
 ForwardDiff = "0.10.3"
 LinearAlgebra = "1.9"
+MaybeInplace = "0.1"
 PrecompileTools = "1"
 Reexport = "1"
 SciMLBase = "2.7"

README.md

@@ -40,7 +40,7 @@ sol = solve(probB, ITP())
 
 For more details on the bracketing methods, refer to the [Tutorials](https://docs.sciml.ai/NonlinearSolve/stable/tutorials/nonlinear/#Using-Bracketing-Methods) and detailed [APIs](https://docs.sciml.ai/NonlinearSolve/stable/api/simplenonlinearsolve/#Solver-API)
 
-## Breaking Changes in v2
+## Breaking Changes in v1.0.0
 
   - Batched solvers have been removed in favor of `BatchedArrays.jl`. Stay tuned for detailed
     tutorials on how to use `BatchedArrays.jl` with `NonlinearSolve` & `SimpleNonlinearSolve`

src/nlsolve/broyden.jl

@@ -7,9 +7,9 @@ and static array problems.
 struct SimpleBroyden <: AbstractSimpleNonlinearSolveAlgorithm end
 
 function SciMLBase.__solve(prob::NonlinearProblem, alg::SimpleBroyden, args...;
-        abstol = nothing, reltol = nothing, maxiters = 1000,
+        abstol = nothing, reltol = nothing, maxiters = 1000, alias_u0 = false,
         termination_condition = nothing, kwargs...)
-    @bb x = copy(float(prob.u0))
+    x = __maybe_unaliased(prob.u0, alias_u0)
     fx = _get_fx(prob, x)
 
     @bb xo = copy(x)

src/nlsolve/dfsane.jl

@@ -51,9 +51,9 @@ Computation, 75, 1429-1448.](https://www.researchgate.net/publication/220576479_
 end
 
 function SciMLBase.__solve(prob::NonlinearProblem, alg::SimpleDFSane, args...;
-        abstol = nothing, reltol = nothing, maxiters = 1000,
+        abstol = nothing, reltol = nothing, maxiters = 1000, alias_u0 = false,
         termination_condition = nothing, kwargs...)
-    x = float(copy(prob.u0))
+    x = __maybe_unaliased(prob.u0, alias_u0)
     fx = _get_fx(prob, x)
     T = eltype(x)
 
@@ -76,11 +76,10 @@ function SciMLBase.__solve(prob::NonlinearProblem, alg::SimpleDFSane, args...;
     history_f_k = fill(fx_norm, M)
 
     # Generate the cache
+    @bb x_cache = similar(x)
     @bb d = copy(x)
     @bb xo = copy(x)
-    @bb x_cache = copy(x)
     @bb δx = copy(x)
-    @bb fxo = copy(fx)
     @bb δf = copy(fx)
 
     k = 0
@@ -91,50 +90,52 @@ function SciMLBase.__solve(prob::NonlinearProblem, alg::SimpleDFSane, args...;
         # Line search direction
         @bb @. d = -σ_k * fx
 
-        η = η_strategy(f_1, k, x, fx)
+        η = η_strategy(f_1, k + 1, x, fx)
         f_bar = maximum(history_f_k)
         α_p = α_1
         α_m = α_1
 
-        @bb @. x += α_p * d
+        @bb @. x_cache = x + α_p * d
 
-        fx = __eval_f(prob, fx, x)
+        fx = __eval_f(prob, fx, x_cache)
         fx_norm_new = NONLINEARSOLVE_DEFAULT_NORM(fx)^nexp
 
         while k < maxiters
-            fx_norm_new ≤ (f_bar + η - γ * α_p^2 * fx_norm) && break
+            Bool(fx_norm_new ≤ (f_bar + η - γ * α_p^2 * fx_norm)) && break
 
-            α_p = α_p^2 * fx_norm / (fx_norm_new + (T(2) * α_p - T(1)) * fx_norm)
-            @bb @. x -= α_m * d
+            α_tp = α_p^2 * fx_norm / (fx_norm_new + (T(2) * α_p - T(1)) * fx_norm)
+            @bb @. x_cache = x - α_m * d
 
-            fx = __eval_f(prob, fx, x)
+            fx = __eval_f(prob, fx, x_cache)
             fx_norm_new = NONLINEARSOLVE_DEFAULT_NORM(fx)^nexp
 
-            fx_norm_new ≤ (f_bar + η - γ * α_m^2 * fx_norm) && break
+            Bool(fx_norm_new ≤ (f_bar + η - γ * α_m^2 * fx_norm)) && break
 
             α_tm = α_m^2 * fx_norm / (fx_norm_new + (T(2) * α_m - T(1)) * fx_norm)
-            α_p = clamp(α_p, τ_min * α_p, τ_max * α_p)
+            α_p = clamp(α_tp, τ_min * α_p, τ_max * α_p)
             α_m = clamp(α_tm, τ_min * α_m, τ_max * α_m)
-            @bb @. x += α_p * d
+            @bb @. x_cache = x + α_p * d
 
-            fx = __eval_f(prob, fx, x)
+            fx = __eval_f(prob, fx, x_cache)
             fx_norm_new = NONLINEARSOLVE_DEFAULT_NORM(fx)^nexp
 
             k += 1
         end
 
+        @bb copyto!(x, x_cache)
+
         tc_sol = check_termination(tc_cache, fx, x, xo, prob, alg)
         tc_sol !== nothing && return tc_sol
 
         # Update spectral parameter
         @bb @. δx = x - xo
-        @bb @. δf = fx - fxo
+        @bb @. δf = fx - δf
 
         σ_k = dot(δx, δx) / dot(δx, δf)
 
         # Take step
         @bb copyto!(xo, x)
-        @bb copyto!(fxo, fx)
+        @bb copyto!(δf, fx)
         fx_norm = fx_norm_new
 
         # Store function value

src/nlsolve/halley.jl

@@ -14,18 +14,22 @@ A low-overhead implementation of Halley's Method.
 
   - `autodiff`: determines the backend used for the Hessian. Defaults to
     `AutoForwardDiff()`. Valid choices are `AutoForwardDiff()` or `AutoFiniteDiff()`.
+
+!!! warning
+
+    Inplace Problems are currently not supported by this method.
 """
 @kwdef @concrete struct SimpleHalley <: AbstractNewtonAlgorithm
     autodiff = AutoForwardDiff()
 end
 
 function SciMLBase.__solve(prob::NonlinearProblem, alg::SimpleHalley, args...;
-        abstol = nothing, reltol = nothing, maxiters = 1000,
+        abstol = nothing, reltol = nothing, maxiters = 1000, alias_u0 = false,
         termination_condition = nothing, kwargs...)
     isinplace(prob) &&
         error("SimpleHalley currently only supports out-of-place nonlinear problems")
 
-    x = copy(float(prob.u0))
+    x = __maybe_unaliased(prob.u0, alias_u0)
     fx = _get_fx(prob, x)
     T = eltype(x)
 

src/nlsolve/klement.jl

@@ -7,9 +7,9 @@ method is non-allocating on scalar and static array problems.
 struct SimpleKlement <: AbstractSimpleNonlinearSolveAlgorithm end
 
 function SciMLBase.__solve(prob::NonlinearProblem, alg::SimpleKlement, args...;
-        abstol = nothing, reltol = nothing, maxiters = 1000,
+        abstol = nothing, reltol = nothing, maxiters = 1000, alias_u0 = false,
         termination_condition = nothing, kwargs...)
-    x = float(prob.u0)
+    x = __maybe_unaliased(prob.u0, alias_u0)
     T = eltype(x)
     fx = _get_fx(prob, x)
 
@@ -21,13 +21,12 @@ function SciMLBase.__solve(prob::NonlinearProblem, alg::SimpleKlement, args...;
     @bb δx = copy(x)
     @bb fprev = copy(fx)
     @bb xo = copy(x)
-    @bb δf = copy(fx)
     @bb d = copy(x)
 
     J = __init_identity_jacobian(fx, x)
-    @bb J_cache = copy(J)
-    @bb δx² = copy(x)
-    @bb J_cache2 = copy(J)
+    @bb J_cache = similar(J)
+    @bb δx² = similar(x)
+    @bb J_cache2 = similar(J)
     @bb F = copy(J)
 
     for _ in 1:maxiters
@@ -45,10 +44,11 @@ function SciMLBase.__solve(prob::NonlinearProblem, alg::SimpleKlement, args...;
             # Singularity test
             if !issuccess(F_)
                 J = __init_identity_jacobian!!(J)
+                @bb copyto!(F, J)
                 if setindex_trait(J) === CanSetindex()
-                    lu!(J; check = false)
+                    F_ = lu!(F; check = false)
                 else
-                    J = lu(J; check = false)
+                    F_ = lu(F; check = false)
                 end
             end
         end
@@ -67,23 +67,18 @@ function SciMLBase.__solve(prob::NonlinearProblem, alg::SimpleKlement, args...;
         tc_sol !== nothing && return tc_sol
 
         @bb δx .*= -1
-        @bb @. δf = fx - fprev
-
-        # Prevent division by 0
+        @bb J_cache .= transpose(J) .^ 2
         @bb @. δx² = δx^2
-        @bb @. J_cache = J^2
-        @bb d = transpose(J_cache) × vec(δx²)
-        @bb @. d = max(d, singular_tol)
-
+        @bb d = J_cache × vec(δx²)
         @bb δx² = J × vec(δx)
-        @bb @. δf = (δf - δx²) / d
-
-        _vδf, _vδx = _vec(δf), _vec(δx)
-        @bb J_cache = _vδf × transpose(_vδx)
+        @bb @. fprev = (fx - fprev - δx²) / ifelse(iszero(d), singular_tol, d)
+        @bb J_cache = vec(fprev) × transpose(_vec(δx))
         @bb @. J_cache *= J
         @bb J_cache2 = J_cache × J
-
         @bb @. J += J_cache2
+
+        @bb copyto!(fprev, fx)
+        @bb copyto!(xo, x)
     end
 
     return build_solution(prob, alg, x, fx; retcode = ReturnCode.MaxIters)

src/nlsolve/lbroyden.jl

@@ -22,9 +22,9 @@ function SimpleLimitedMemoryBroyden(; threshold::Union{Val, Int} = Val(27))
 end
 
 @views function SciMLBase.__solve(prob::NonlinearProblem, alg::SimpleLimitedMemoryBroyden,
-        args...; abstol = nothing, reltol = nothing, maxiters = 1000,
+        args...; abstol = nothing, reltol = nothing, maxiters = 1000, alias_u0 = false,
         termination_condition = nothing, kwargs...)
-    @bb x = copy(float(prob.u0))
+    x = __maybe_unaliased(prob.u0, alias_u0)
     threshold = __get_threshold(alg)
     η = min(SciMLBase._unwrap_val(threshold), maxiters)
 

src/nlsolve/raphson.jl

@@ -24,8 +24,8 @@ const SimpleGaussNewton = SimpleNewtonRaphson
 
 function SciMLBase.__solve(prob::Union{NonlinearProblem, NonlinearLeastSquaresProblem},
         alg::SimpleNewtonRaphson, args...; abstol = nothing, reltol = nothing,
-        maxiters = 1000, termination_condition = nothing, kwargs...)
-    @bb x = copy(float(prob.u0))
+        maxiters = 1000, termination_condition = nothing, alias_u0 = false, kwargs...)
+    x = __maybe_unaliased(prob.u0, alias_u0)
     fx = _get_fx(prob, x)
     @bb xo = copy(x)
     J, jac_cache = jacobian_cache(alg.autodiff, prob.f, fx, x, prob.p)
@@ -37,9 +37,7 @@ function SciMLBase.__solve(prob::Union{NonlinearProblem, NonlinearLeastSquaresPr
         fx, dfx = value_and_jacobian(alg.autodiff, prob.f, fx, x, prob.p, jac_cache; J)
 
         if i == 1
-            if iszero(fx)
-                return build_solution(prob, alg, x, fx; retcode = ReturnCode.Success)
-            end
+            iszero(fx) && build_solution(prob, alg, x, fx; retcode = ReturnCode.Success)
         else
             # Termination Checks
             tc_sol = check_termination(tc_cache, fx, x, xo, prob, alg)

src/nlsolve/trustRegion.jl

@@ -49,9 +49,9 @@ scalar and static array problems.
 end
 
 function SciMLBase.__solve(prob::NonlinearProblem, alg::SimpleTrustRegion, args...;
-        abstol = nothing, reltol = nothing, maxiters = 1000,
+        abstol = nothing, reltol = nothing, maxiters = 1000, alias_u0 = false,
         termination_condition = nothing, kwargs...)
-    @bb x = copy(float(prob.u0))
+    x = __maybe_unaliased(prob.u0, alias_u0)
     T = eltype(real(x))
     Δₘₐₓ = T(alg.max_trust_radius)
     Δ = T(alg.initial_trust_radius)
@@ -85,7 +85,6 @@ function SciMLBase.__solve(prob::NonlinearProblem, alg::SimpleTrustRegion, args.
     @bb Hδ = copy(x)
     dogleg_cache = (; δsd, δN_δsd, δN)
 
-    F = fx
     for k in 1:maxiters
         # Solve the trust region subproblem.
         δ = dogleg_method!!(dogleg_cache, ∇f, fx, g, Δ)
@@ -97,19 +96,19 @@ function SciMLBase.__solve(prob::NonlinearProblem, alg::SimpleTrustRegion, args.
 
         # Compute the ratio of the actual to predicted reduction.
         @bb Hδ = H × vec(δ)
-        r = (fₖ₊₁ - fₖ) / (dot(δ', g) + dot(δ', Hδ) / T(2))
+        r = (fₖ₊₁ - fₖ) / (dot(δ, g) + dot(δ, Hδ) / T(2))
 
         # Update the trust region radius.
-        if r < η₂
+        if r ≥ η₂
+            shrink_counter = 0
+        else
             Δ = t₁ * Δ
             shrink_counter += 1
             shrink_counter > max_shrink_times && return build_solution(prob, alg, x, fx;
                 retcode = ReturnCode.ConvergenceFailure)
-        else
-            shrink_counter = 0
         end
 
-        if r > η₁
+        if r ≥ η₁
             # Termination Checks
             tc_sol = check_termination(tc_cache, fx, x, xo, prob, alg)
             tc_sol !== nothing && return tc_sol
@@ -144,6 +143,7 @@ function dogleg_method!!(cache, J, f, g, Δ)
     @bb δsd .= g
     @bb @. δsd *= -1
     norm_δsd = norm(δsd)
+
     if (norm_δsd ≥ Δ)
         @bb @. δsd *= Δ / norm_δsd
         return δsd

src/utils.jl

@@ -286,14 +286,14 @@ function check_termination(tc_cache, fx, x, xo, prob, alg)
 end
 function check_termination(tc_cache, fx, x, xo, prob, alg,
         ::AbstractNonlinearTerminationMode)
-    if tc_cache(fx, x, xo)
+    if Bool(tc_cache(fx, x, xo))
         return build_solution(prob, alg, x, fx; retcode = ReturnCode.Success)
     end
     return nothing
 end
 function check_termination(tc_cache, fx, x, xo, prob, alg,
         ::AbstractSafeNonlinearTerminationMode)
-    if tc_cache(fx, x, xo)
+    if Bool(tc_cache(fx, x, xo))
         if tc_cache.retcode == NonlinearSafeTerminationReturnCode.Success
             retcode = ReturnCode.Success
         elseif tc_cache.retcode == NonlinearSafeTerminationReturnCode.PatienceTermination
@@ -309,7 +309,7 @@ function check_termination(tc_cache, fx, x, xo, prob, alg,
 end
 function check_termination(tc_cache, fx, x, xo, prob, alg,
         ::AbstractSafeBestNonlinearTerminationMode)
-    if tc_cache(fx, x, xo)
+    if Bool(tc_cache(fx, x, xo))
         if tc_cache.retcode == NonlinearSafeTerminationReturnCode.Success
             retcode = ReturnCode.Success
         elseif tc_cache.retcode == NonlinearSafeTerminationReturnCode.PatienceTermination
@@ -335,3 +335,11 @@ end
 
 @inline __eval_f(prob, fx, x) = isinplace(prob) ? (prob.f(fx, x, prob.p); fx) :
                                 prob.f(x, prob.p)
+
+# Unalias
+@inline __maybe_unaliased(x::Union{Number, SArray}, ::Bool) = x
+@inline function __maybe_unaliased(x::AbstractArray, alias::Bool)
+    # Spend time coping iff we will mutate the array
+    (alias || !ArrayInterface.can_setindex(typeof(x))) && return x
+    return deepcopy(x)
+end

test/23_test_problems.jl

@@ -61,7 +61,7 @@ end
     alg_ops = (SimpleDFSane(),)
 
     broken_tests = Dict(alg => Int[] for alg in alg_ops)
-    broken_tests[alg_ops[1]] = [1, 2, 3, 4, 5, 6, 7, 9, 11, 12, 13, 15, 16, 17, 21, 22]
+    broken_tests[alg_ops[1]] = [1, 2, 3, 4, 5, 6, 11, 21]
 
     test_on_library(problems, dicts, alg_ops, broken_tests)
 end
@@ -82,7 +82,7 @@ end
     alg_ops = (SimpleKlement(),)
 
     broken_tests = Dict(alg => Int[] for alg in alg_ops)
-    broken_tests[alg_ops[1]] = [1, 2, 4, 5, 6, 7, 9, 10, 11, 12, 13, 19, 21, 22]
+    broken_tests[alg_ops[1]] = [1, 2, 4, 5, 6, 7, 11, 13, 22]
 
     test_on_library(problems, dicts, alg_ops, broken_tests)
 end

test/basictests.jl

@@ -198,11 +198,11 @@ end
             res = benchmark_nlsolve_oop(quadratic_f, @SVector[1.0, 1.0], p)
 
             if any(x -> isnan(x) || x <= 1e-5 || x >= 1e5, res)
-                @test_broken all(res .≈ sqrt(p))
+                @test_broken all(abs.(res) .≈ sqrt(p))
                 @test_broken abs.(ForwardDiff.derivative(p -> benchmark_nlsolve_oop(quadratic_f,
                     @SVector[1.0, 1.0], p).u[end], p)) ≈ 1 / (2 * sqrt(p))
             else
-                @test all(res .≈ sqrt(p))
+                @test all(abs.(res) .≈ sqrt(p))
                 @test isapprox(abs.(ForwardDiff.derivative(p -> benchmark_nlsolve_oop(quadratic_f,
                         @SVector[1.0, 1.0], p).u[end], p)), 1 / (2 * sqrt(p)))
             end
@@ -213,12 +213,12 @@ end
         for p in 1.0:0.1:100.0
             res = benchmark_nlsolve_oop(quadratic_f, 1.0, p)
 
-            if any(x -> isnan(x) || x <= 1e-5 || x >= 1e5, res)
-                @test_broken all(res .≈ sqrt(p))
+            if any(x -> isnan(x), res)
+                @test_broken abs(res.u) ≈ sqrt(p)
                 @test_broken abs.(ForwardDiff.derivative(p -> benchmark_nlsolve_oop(quadratic_f,
                         1.0, p).u, p)) ≈ 1 / (2 * sqrt(p))
             else
-                @test all(res .≈ sqrt(p))
+                @test abs(res.u) ≈ sqrt(p)
                 @test isapprox(abs.(ForwardDiff.derivative(p -> benchmark_nlsolve_oop(quadratic_f,
                             1.0, p).u, p)), 1 / (2 * sqrt(p)))
             end