Merge pull request #304 from avik-pal/ap/fix_docs

Fix storing the trace

Author: ChrisRackauckas

Project.toml

@@ -75,6 +75,7 @@ SciMLOperators = "0.3"
 SimpleNonlinearSolve = "0.1.23"
 SparseArrays = "<0.0.1, 1"
 SparseDiffTools = "2.14"
+StableRNGs = "1"
 StaticArrays = "1"
 Symbolics = "5"
 Test = "1"
@@ -98,11 +99,12 @@ NonlinearProblemLibrary = "b7050fa9-e91f-4b37-bcee-a89a063da141"
 Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f"
+StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3"
 SparseDiffTools = "47a9eef4-7e08-11e9-0b38-333d64bd3804"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f"
 
 [targets]
-test = ["Aqua", "Enzyme", "BenchmarkTools", "SafeTestsets", "Pkg", "Test", "ForwardDiff", "StaticArrays", "Symbolics", "LinearSolve", "Random", "LinearAlgebra", "Zygote", "SparseDiffTools", "NonlinearProblemLibrary", "LeastSquaresOptim", "FastLevenbergMarquardt", "NaNMath", "BandedMatrices", "DiffEqBase"]
+test = ["Aqua", "Enzyme", "BenchmarkTools", "SafeTestsets", "Pkg", "Test", "ForwardDiff", "StaticArrays", "Symbolics", "LinearSolve", "Random", "LinearAlgebra", "Zygote", "SparseDiffTools", "NonlinearProblemLibrary", "LeastSquaresOptim", "FastLevenbergMarquardt", "NaNMath", "BandedMatrices", "DiffEqBase", "StableRNGs"]

docs/Project.toml

@@ -30,7 +30,7 @@ NonlinearSolve = "1, 2"
 NonlinearSolveMINPACK = "0.1"
 SciMLBase = "2.4"
 SciMLNLSolve = "0.1"
-SimpleNonlinearSolve = "0.1.5"
+SimpleNonlinearSolve = "0.1.5, 1"
 StaticArrays = "1"
 SteadyStateDiffEq = "1.10, 2"
 Sundials = "4.11"

src/dfsane.jl

@@ -6,9 +6,7 @@
 
 A low-overhead and allocation-free implementation of the df-sane method for solving large-scale nonlinear
 systems of equations. For in depth information about all the parameters and the algorithm,
-see the paper: [W LaCruz, JM Martinez, and M Raydan (2006), Spectral Residual Method without
-Gradient Information for Solving Large-Scale Nonlinear Systems of Equations, Mathematics of
-Computation, 75, 1429-1448.](https://www.researchgate.net/publication/220576479_Spectral_Residual_Method_without_Gradient_Information_for_Solving_Large-Scale_Nonlinear_Systems_of_Equations)
+see the paper [1].
 
 ### Keyword Arguments
 
@@ -40,6 +38,12 @@ Computation, 75, 1429-1448.](https://www.researchgate.net/publication/220576479_
     ``fn_1 / n^2``.
   - `max_inner_iterations`: the maximum number of iterations allowed for the inner loop of the
     algorithm. Defaults to `100`.
+
+### References
+
+[1] W LaCruz, JM Martinez, and M Raydan (2006), Spectral Residual Method without Gradient
+Information for Solving Large-Scale Nonlinear Systems of Equations, Mathematics of
+Computation, 75, 1429-1448.
 """
 @kwdef @concrete struct DFSane <: AbstractNonlinearSolveAlgorithm
     σ_min = 1e-10

src/trace.jl

@@ -184,7 +184,7 @@ function __init_trace_history(::Val{show_trace}, trace_level, ::Val{store_trace}
     !store_trace && !show_trace && return nothing
     entry = __trace_entry(trace_level, 0, u, fu, J, δu)
     show_trace && show(entry)
-    store_trace && return [entry]
+    store_trace && return NonlinearSolveTraceEntry[entry]
     return nothing
 end
 
@@ -218,29 +218,6 @@ function update_trace!(trace::NonlinearSolveTrace{ShT, StT}, iter, u, fu, J, δu
     return trace
 end
 
-# Needed for Algorithms which directly use `inv(J)` instead of `J`
-function update_trace_with_invJ!(trace::NonlinearSolveTrace{ShT, StT}, iter, u, fu, J, δu,
-        α = 1; last::Val{L} = Val(false)) where {ShT, StT, L}
-    !StT && !ShT && return nothing
-
-    if L
-        entry = NonlinearSolveTraceEntry(-1, norm(fu, Inf), NaN32, nothing, nothing,
-            nothing, nothing, nothing)
-        show(entry)
-        return trace
-    end
-
-    show_now = ShT && (mod1(iter, trace.trace_level.print_frequency) == 1)
-    store_now = StT && (mod1(iter, trace.trace_level.store_frequency) == 1)
-    if show_now || store_now
-        J_ = trace.trace_level isa TraceMinimal ? J : inv(J)
-        entry = __trace_entry(trace.trace_level, iter, u, fu, J_, δu, α)
-    end
-    store_now && push!(trace.history, entry)
-    show_now && show(entry)
-    return trace
-end
-
 function update_trace!(cache::AbstractNonlinearSolveCache, α = true)
     trace = __getproperty(cache, Val(:trace))
     trace === nothing && return nothing
@@ -252,8 +229,8 @@ function update_trace!(cache::AbstractNonlinearSolveCache, α = true)
             update_trace!(trace, cache.stats.nsteps + 1, get_u(cache), get_fu(cache),
                 nothing, cache.du, α)
         else
-            update_trace_with_invJ!(trace, cache.stats.nsteps + 1, get_u(cache),
-                get_fu(cache), J_inv, cache.du, α)
+            update_trace!(trace, cache.stats.nsteps + 1, get_u(cache), get_fu(cache),
+                ApplyArray(inv, J_inv), cache.du, α)
         end
     else
         update_trace!(trace, cache.stats.nsteps + 1, get_u(cache), get_fu(cache), J,

src/trustRegion.jl

@@ -1,34 +1,36 @@
 """
     RadiusUpdateSchemes
 
-`RadiusUpdateSchemes` is the standard enum interface for different types of radius update schemes
-implemented in the Trust Region method. These schemes specify how the radius of the so-called trust region
-is updated after each iteration of the algorithm. The specific role and caveats associated with each
-scheme are provided below.
+`RadiusUpdateSchemes` is the standard enum interface for different types of radius update
+schemes implemented in the Trust Region method. These schemes specify how the radius of the
+so-called trust region is updated after each iteration of the algorithm. The specific role
+and caveats associated with each scheme are provided below.
 
 ## Using `RadiusUpdateSchemes`
 
-`RadiusUpdateSchemes` uses the standard EnumX interface (https://github.com/fredrikekre/EnumX.jl),
-and hence inherits all properties of being an EnumX, including the type of each constituent enum
-states as `RadiusUpdateSchemes.T`. Simply put the desired scheme as follows:
+`RadiusUpdateSchemes` uses the standard
+[EnumX Interface](https://github.com/fredrikekre/EnumX.jl), and hence inherits all
+properties of being an EnumX, including the type of each constituent enum states as
+`RadiusUpdateSchemes.T`. Simply put the desired scheme as follows:
 `TrustRegion(radius_update_scheme = your desired update scheme)`. For example,
 `sol = solve(prob, alg=TrustRegion(radius_update_scheme = RadiusUpdateSchemes.Hei))`.
 """
 @enumx RadiusUpdateSchemes begin
     """
         RadiusUpdateSchemes.Simple
 
-    The simple or conventional radius update scheme. This scheme is chosen by default
-    and follows the conventional approach to update the trust region radius, i.e. if the
-    trial step is accepted it increases the radius by a fixed factor (bounded by a maximum radius)
+    The simple or conventional radius update scheme. This scheme is chosen by default and
+    follows the conventional approach to update the trust region radius, i.e. if the trial
+    step is accepted it increases the radius by a fixed factor (bounded by a maximum radius)
     and if the trial step is rejected, it shrinks the radius by a fixed factor.
     """
     Simple
 
     """
         RadiusUpdateSchemes.NLsolve
 
-    The same updating scheme as in NLsolve's (https://github.com/JuliaNLSolvers/NLsolve.jl) trust region dogleg implementation.
+    The same updating scheme as in NLsolve's (https://github.com/JuliaNLSolvers/NLsolve.jl)
+    trust region dogleg implementation.
     """
     NLsolve
 
@@ -42,40 +44,57 @@ states as `RadiusUpdateSchemes.T`. Simply put the desired scheme as follows:
     """
         RadiusUpdateSchemes.Hei
 
-    This scheme is proposed by [Hei, L.] (https://www.jstor.org/stable/43693061). The trust region radius
-    depends on the size (norm) of the current step size. The hypothesis is to let the radius converge to zero
-    as the iterations progress, which is more reliable and robust for ill-conditioned as well as degenerate
-    problems.
+    This scheme is proposed by Hei, L. [1]. The trust region radius depends on the size
+    (norm) of the current step size. The hypothesis is to let the radius converge to zero as
+    the iterations progress, which is more reliable and robust for ill-conditioned as well
+    as degenerate problems.
+
+    [1] Hei, Long. "A self-adaptive trust region algorithm." Journal of Computational
+    Mathematics (2003): 229-236.
     """
     Hei
 
     """
         RadiusUpdateSchemes.Yuan
 
-    This scheme is proposed by [Yuan, Y.] (https://www.researchgate.net/publication/249011466_A_new_trust_region_algorithm_with_trust_region_radius_converging_to_zero).
-    Similar to Hei's scheme, the trust region is updated in a way so that it converges to zero, however here,
-    the radius depends on the size (norm) of the current gradient of the objective (merit) function. The hypothesis
-    is that the step size is bounded by the gradient size, so it makes sense to let the radius depend on the gradient.
+    This scheme is proposed by Yuan, Y [1]. Similar to Hei's scheme, the trust region is
+    updated in a way so that it converges to zero, however here, the radius depends on the
+    size (norm) of the current gradient of the objective (merit) function. The hypothesis is
+    that the step size is bounded by the gradient size, so it makes sense to let the radius
+    depend on the gradient.
+
+    [1] Fan, Jinyan, Jianyu Pan, and Hongyan Song. "A retrospective trust region algorithm
+    with trust region converging to zero." Journal of Computational Mathematics 34.4 (2016):
+    421-436.
     """
     Yuan
 
     """
         RadiusUpdateSchemes.Bastin
 
-    This scheme is proposed by [Bastin, et al.] (https://www.researchgate.net/publication/225100660_A_retrospective_trust-region_method_for_unconstrained_optimization).
-    The scheme is called a retrospective update scheme as it uses the model function at the current
-    iteration to compute the ratio of the actual reduction and the predicted reduction in the previous
-    trial step, and use this ratio to update the trust region radius. The hypothesis is to exploit the information
-    made available during the optimization process in order to vary the accuracy of the objective function computation.
+    This scheme is proposed by Bastin, et al. [1]. The scheme is called a retrospective
+    update scheme as it uses the model function at the current iteration to compute the
+    ratio of the actual reduction and the predicted reduction in the previous trial step,
+    and use this ratio to update the trust region radius. The hypothesis is to exploit the
+    information made available during the optimization process in order to vary the accuracy
+    of the objective function computation.
+
+    [1] Bastin, Fabian, et al. "A retrospective trust-region method for unconstrained
+    optimization." Mathematical programming 123 (2010): 395-418.
     """
     Bastin
 
     """
         RadiusUpdateSchemes.Fan
 
-    This scheme is proposed by [Fan, J.] (https://link.springer.com/article/10.1007/s10589-005-3078-8). It is very much similar to
-    Hei's and Yuan's schemes as it lets the trust region radius depend on the current size (norm) of the objective (merit)
-    function itself. These new update schemes are known to improve local convergence.
+    This scheme is proposed by Fan, J. [1]. It is very much similar to Hei's and Yuan's
+    schemes as it lets the trust region radius depend on the current size (norm) of the
+    objective (merit) function itself. These new update schemes are known to improve local
+    convergence.
+
+    [1] Fan, Jinyan. "Convergence rate of the trust region method for nonlinear equations
+    under local error bound condition." Computational Optimization and Applications 34.2
+    (2006): 215-227.
     """
     Fan
 end

test/nonlinear_least_squares.jl

@@ -1,4 +1,5 @@
-using NonlinearSolve, LinearSolve, LinearAlgebra, Test, Random, ForwardDiff, Zygote
+using NonlinearSolve,
+    LinearSolve, LinearAlgebra, Test, StableRNGs, Random, ForwardDiff, Zygote
 import FastLevenbergMarquardt, LeastSquaresOptim
 
 true_function(x, θ) = @. θ[1] * exp(θ[2] * x) * cos(θ[3] * x + θ[4])
@@ -21,7 +22,7 @@ function loss_function(resid, θ, p)
     return resid
 end
 
-θ_init = θ_true .+ randn!(similar(θ_true)) * 0.1
+θ_init = θ_true .+ randn!(StableRNG(0), similar(θ_true)) * 0.1
 prob_oop = NonlinearLeastSquaresProblem{false}(loss_function, θ_init, x)
 prob_iip = NonlinearLeastSquaresProblem(NonlinearFunction(loss_function;
         resid_prototype = zero(y_target)), θ_init, x)

test/runtests.jl

@@ -14,11 +14,11 @@ end
 @time begin
     if GROUP == "All" || GROUP == "Core"
         @time @safetestset "Quality Assurance" include("qa.jl")
+        @time @safetestset "Nonlinear Least Squares" include("nonlinear_least_squares.jl")
         @time @safetestset "Basic Tests + Some AD" include("basictests.jl")
         @time @safetestset "Sparsity Tests" include("sparse.jl")
         @time @safetestset "Polyalgs" include("polyalgs.jl")
         @time @safetestset "Matrix Resizing" include("matrix_resizing.jl")
-        @time @safetestset "Nonlinear Least Squares" include("nonlinear_least_squares.jl")
         @time @safetestset "Infeasible Problems" include("infeasible.jl")
     end