Update the FAQ and tutorials

avik-pal · avik-pal · commit 2a6f1a98ec8f · 2024-01-11T13:36:23.000-05:00
diff --git a/docs/Project.toml b/docs/Project.toml
@@ -6,6 +6,7 @@ DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 DocumenterCitations = "daee34ce-89f3-4625-b898-19384cb65244"
 IncompleteLU = "40713840-3770-5561-ab4c-a76e7d0d7895"
+InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
 LinearSolve = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae"
 ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78"
 NonlinearSolve = "8913a72c-1f9b-4ce2-8d82-65094dcecaec"
diff --git a/docs/src/basics/diagnostics_api.md b/docs/src/basics/diagnostics_api.md
@@ -1,3 +1,5 @@
+# [Diagnostics API](@id diagnostics_api)
+
 # Logging the Solve Process
 
 All NonlinearSolve.jl native solvers allow storing and displaying the trace of the nonlinear
diff --git a/docs/src/basics/faq.md b/docs/src/basics/faq.md
@@ -76,7 +76,8 @@ sol = solve(prob_oop, LevenbergMarquardt(; autodiff = AutoFiniteDiff()); maxiter
 ```
 
 This worked but, Finite Differencing is not the recommended approach in any scenario.
-Instead, rewrite the function to use
+
+  2. Rewrite the function to use
 [PreallocationTools.jl](https://github.com/SciML/PreallocationTools.jl) or write it as
 
 ```@example dual_error_faq
@@ -93,4 +94,58 @@ sol = solve(prob_oop, LevenbergMarquardt(); maxiters = 10000, abstol = 1e-8)
 
 ## I thought NonlinearSolve.jl was type-stable and fast. But it isn't, why?
 
+It is hard to say why your code is not fast. Take a look at the
+[Diagnostics API](@ref diagnostics_api) to pin-point the problem. One common issue is that
+there is type instability.
+
+If you are using the defaults for the autodiff and your problem is not a scalar or using
+static arrays, ForwardDiff will create type unstable code. See this simple example:
+
+```@example type_unstable
+using NonlinearSolve, InteractiveUtils
+
+f(u, p) = @. u^2 - p
+
+prob = NonlinearProblem{false}(f, 1.0, 2.0)
+
+@code_warntype solve(prob, NewtonRaphson())
+nothing # hide
+```
+
+Notice that this was type-stable, since it is a scalar problem. Now what happens for static
+arrays
+
+```@example type_unstable
+using StaticArrays
+
+prob = NonlinearProblem{false}(f, @SVector([1.0, 2.0]), 2.0)
+
+@code_warntype solve(prob, NewtonRaphson())
+nothing # hide
+```
+
+Again Type-Stable! Now let's try using a regular array:
+
+```@example type_unstable
+prob = NonlinearProblem(f, [1.0, 2.0], 2.0)
+
+@code_warntype solve(prob, NewtonRaphson())
+nothing # hide
+```
+
+Oh no! This is type unstable. This is because ForwardDiff.jl will chunk the jacobian
+computation and the type of this chunksize can't be statically inferred. To fix this, we
+directly specify the chunksize:
+
+```@example type_unstable
+@code_warntype solve(prob, NewtonRaphson(; autodiff = AutoForwardDiff(; chunksize = 2)))
+nothing # hide
+```
+
+And boom! Type stable again. For selecting the chunksize the method is:
+
+  1. For small inputs `≤ 12` use `chunksize = <length of input>`
+  2. For larger inputs, use `chunksize = 12`
 
+In general, the chunksize should be `≤ length of input`. However, a very large chunksize
+can lead to excessive compilation times and slowdown.
diff --git a/docs/src/basics/nonlinear_problem.md b/docs/src/basics/nonlinear_problem.md
@@ -35,7 +35,7 @@ that `f(u) = 0`, the `NonlinearProblem` does not have a preferred solution, whil
 `SteadyStateProblem` the preferred solution is the `u(∞)` that would arise from solving the
 ODE `u' = f(u,t)`.
 
-!!! warn
+!!! warning
     
     Most solvers for `SteadyStateProblem` do not guarantee the preferred solution and
     instead will solve for some `u` in the set of solutions. The documentation of the
diff --git a/docs/src/native/globalization.md b/docs/src/native/globalization.md
@@ -12,6 +12,7 @@ Pages = ["globalization.md"]
 LiFukushimaLineSearch
 LineSearchesJL
 RobustNonMonotoneLineSearch
+NoLineSearch
 ```
 
 ## Radius Update Schemes for Trust Region
diff --git a/docs/src/refs.bib b/docs/src/refs.bib
@@ -1,3 +1,13 @@
+@article{bastin2010retrospective,
+  title     = {A retrospective trust-region method for unconstrained optimization},
+  author    = {Bastin, Fabian and Malmedy, Vincent and Mouffe, M{\'e}lodie and Toint, Philippe L and Tomanos, Dimitri},
+  journal   = {Mathematical programming},
+  volume    = {123},
+  pages     = {395--418},
+  year      = {2010},
+  publisher = {Springer}
+}
+
 @article{broyden1965class,
   title   = {A class of methods for solving nonlinear simultaneous equations},
   author  = {Broyden, Charles G},
@@ -19,6 +29,37 @@ @article{coffey2003pseudotransient
   publisher = {SIAM}
 }
 
+@article{fan2006convergence,
+  title     = {Convergence rate of the trust region method for nonlinear equations under local error bound condition},
+  author    = {Fan, Jinyan},
+  journal   = {Computational Optimization and Applications},
+  volume    = {34},
+  number    = {2},
+  pages     = {215--227},
+  year      = {2006},
+  publisher = {Springer}
+}
+
+@article{fan2016retrospective,
+  title     = {A retrospective trust region algorithm with trust region converging to zero},
+  author    = {Fan, Jinyan and Pan, Jianyu and Song, Hongyan},
+  journal   = {Journal of Computational Mathematics},
+  volume    = {34},
+  number    = {4},
+  pages     = {421--436},
+  year      = {2016},
+  publisher = {JSTOR}
+}
+
+@article{hei2003self,
+  title     = {A self-adaptive trust region algorithm},
+  author    = {Hei, Long},
+  journal   = {Journal of Computational Mathematics},
+  pages     = {229--236},
+  year      = {2003},
+  publisher = {JSTOR}
+}
+
 @article{kelley1998convergence,
   title     = {Convergence analysis of pseudo-transient continuation},
   author    = {Kelley, Carl Timothy and Keyes, David E},
@@ -78,6 +119,16 @@ @article{yuan2015recent
   publisher = {Springer}
 }
 
+@article{yuan2015recent,
+  title     = {Recent advances in trust region algorithms},
+  author    = {Yuan, Ya-xiang},
+  journal   = {Mathematical Programming},
+  volume    = {151},
+  pages     = {249--281},
+  year      = {2015},
+  publisher = {Springer}
+}
+
 @article{ziani2008autoadaptative,
   title     = {An autoadaptative limited memory Broyden’s method to solve systems of nonlinear equations},
   author    = {Ziani, Mohammed and Guyomarc’h, Fr{\'e}d{\'e}ric},
diff --git a/docs/src/tutorials/large_systems.md b/docs/src/tutorials/large_systems.md
@@ -137,11 +137,14 @@ Symbolic Sparsity Detection. See the manual entry on
 using BenchmarkTools # for @btime
 
 @btime solve(prob_brusselator_2d, NewtonRaphson());
-@btime solve(prob_brusselator_2d, NewtonRaphson(; autodiff = AutoSparseForwardDiff()));
 @btime solve(prob_brusselator_2d,
-    NewtonRaphson(; autodiff = AutoSparseForwardDiff(), linsolve = KLUFactorization()));
+    NewtonRaphson(; autodiff = AutoSparseForwardDiff(; chunksize = 32)));
 @btime solve(prob_brusselator_2d,
-    NewtonRaphson(; autodiff = AutoSparseForwardDiff(), linsolve = KrylovJL_GMRES()));
+    NewtonRaphson(; autodiff = AutoSparseForwardDiff(; chunksize = 32),
+        linsolve = KLUFactorization()));
+@btime solve(prob_brusselator_2d,
+    NewtonRaphson(; autodiff = AutoSparseForwardDiff(; chunksize = 32),
+        linsolve = KrylovJL_GMRES()));
 nothing # hide
 ```
 
@@ -175,7 +178,7 @@ ff = NonlinearFunction(brusselator_2d_loop; sparsity = jac_sparsity)
 Build the `NonlinearProblem`:
 
 ```@example ill_conditioned_nlprob
-prob_brusselator_2d_sparse = NonlinearProblem(ff, u0, p)
+prob_brusselator_2d_sparse = NonlinearProblem(ff, u0, p; abstol = 1e-10, reltol = 1e-10)
 ```
 
 Now let's see how the version with sparsity compares to the version without:
diff --git a/docs/src/tutorials/optimizing_parameterized_ode.md b/docs/src/tutorials/optimizing_parameterized_ode.md
@@ -54,7 +54,7 @@ Now, we can use any NLLS solver to solve this problem.
 
 ```@example parameterized_ode
 res = solve(nlls_prob, LevenbergMarquardt(); maxiters = 1000, show_trace = Val(true),
-    trace_level = TraceAll())
+    trace_level = TraceWithJacobianConditionNumber(25))
 nothing # hide
 ```
 
@@ -66,7 +66,7 @@ We can also use Trust Region methods.
 
 ```@example parameterized_ode
 res = solve(nlls_prob, TrustRegion(); maxiters = 1000, show_trace = Val(true),
-    trace_level = TraceAll())
+    trace_level = TraceWithJacobianConditionNumber(25))
 nothing # hide
 ```
 
diff --git a/src/NonlinearSolve.jl b/src/NonlinearSolve.jl
@@ -10,7 +10,7 @@ import PrecompileTools: @recompile_invalidations, @compile_workload, @setup_work
 @recompile_invalidations begin
     using ADTypes, ConcreteStructs, DiffEqBase, FastBroadcast, FastClosures, LazyArrays,
         LineSearches, LinearAlgebra, LinearSolve, MaybeInplace, Preferences, Printf,
-        SciMLBase, SimpleNonlinearSolve, SparseArrays, SparseDiffTools, TimerOutputs
+        SciMLBase, SimpleNonlinearSolve, SparseArrays, SparseDiffTools
 
     import ArrayInterface: undefmatrix, can_setindex, restructure, fast_scalar_indexing
     import DiffEqBase: AbstractNonlinearTerminationMode,
diff --git a/src/globalization/trust_region.jl b/src/globalization/trust_region.jl
@@ -21,12 +21,7 @@ of specifying a trust region radius.
     iteration ``i``. Reasonable choices for `b_uphill` are `1.0` or `2.0`, with
     `b_uphill = 2.0` allowing higher uphill moves than `b_uphill = 1.0`. When
     `b_uphill = 0.0`, no uphill moves will be accepted. Defaults to `1.0`. See Section 4 of
-    [1].
-
-### References
-
-[1] Transtrum, Mark K., and James P. Sethna. "Improvements to the Levenberg-Marquardt
-algorithm for nonlinear least-squares minimization." arXiv preprint arXiv:1201.5885 (2012).
+    [transtrum2012improvements](@ref).
 """
 @concrete struct LevenbergMarquardtTrustRegion <: AbstractTrustRegionMethod
     β_uphill
@@ -120,6 +115,7 @@ end
 
 const T = AbstractRadiusUpdateScheme
 
+struct __Simple <: AbstractRadiusUpdateScheme end
 """
     RadiusUpdateSchemes.Simple
 
@@ -128,93 +124,70 @@ follows the conventional approach to update the trust region radius, i.e. if the
 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.
 """
-struct __Simple <: AbstractRadiusUpdateScheme end
 const Simple = __Simple()
 
+struct __NLsolve <: AbstractRadiusUpdateScheme end
 """
     RadiusUpdateSchemes.NLsolve
 
 The same updating scheme as in NLsolve's (https://github.com/JuliaNLSolvers/NLsolve.jl)
 trust region dogleg implementation.
 """
-struct __NLsolve <: AbstractRadiusUpdateScheme end
 const NLsolve = __NLsolve()
 
+struct __NocedalWright <: AbstractRadiusUpdateScheme end
 """
     RadiusUpdateSchemes.NocedalWright
 
 Trust region updating scheme as in Nocedal and Wright [see Alg 11.5, page 291].
 """
-struct __NocedalWright <: AbstractRadiusUpdateScheme end
 const NocedalWright = __NocedalWright()
 
+struct __Hei <: AbstractRadiusUpdateScheme end
 """
     RadiusUpdateSchemes.Hei
 
-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
+This scheme is proposed in [hei2003self](@citet). 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.
-
-### References
-
-[1] Hei, Long. "A self-adaptive trust region algorithm." Journal of Computational
-Mathematics (2003): 229-236.
 """
-struct __Hei <: AbstractRadiusUpdateScheme end
 const Hei = __Hei()
 
+struct __Yuan <: AbstractRadiusUpdateScheme end
 """
     RadiusUpdateSchemes.Yuan
 
-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.
-
-### References
-
-[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.
+This scheme is proposed by [yuan2015recent](@citet). 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.
 """
-struct __Yuan <: AbstractRadiusUpdateScheme end
 const Yuan = __Yuan()
 
+struct __Bastin <: AbstractRadiusUpdateScheme end
 """
     RadiusUpdateSchemes.Bastin
 
-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
+This scheme is proposed by [bastin2010retrospective](@citet). 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.
-
-### References
-
-[1] Bastin, Fabian, et al. "A retrospective trust-region method for unconstrained
-optimization." Mathematical programming 123 (2010): 395-418.
 """
-struct __Bastin <: AbstractRadiusUpdateScheme end
 const Bastin = __Bastin()
 
+struct __Fan <: AbstractRadiusUpdateScheme end
 """
     RadiusUpdateSchemes.Fan
 
-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
+This scheme is proposed by [fan2006convergence](@citet). 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.
-
-### References
-
-[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.
 """
-struct __Fan <: AbstractRadiusUpdateScheme end
 const Fan = __Fan()
 
 end
@@ -248,13 +221,11 @@ the value used in the respective paper.
   - `step_threshold`: the threshold for taking a step. In every iteration, the threshold is
     compared with a value `r`, which is the actual reduction in the objective function
     divided by the predicted reduction. If `step_threshold > r` the model is not a good
-    approximation, and the step is rejected. Defaults to `nothing`. For more details, see
-    [2].
+    approximation, and the step is rejected. Defaults to `nothing`.
   - `shrink_threshold`: the threshold for shrinking the trust region radius. In every
     iteration, the threshold is compared with a value `r` which is the actual reduction in
     the objective function divided by the predicted reduction. If `shrink_threshold > r` the
-    trust region radius is shrunk by `shrink_factor`. Defaults to `nothing`. For more
-    details, see [2].
+    trust region radius is shrunk by `shrink_factor`. Defaults to `nothing`.
   - `expand_threshold`: the threshold for expanding the trust region radius. If a step is
     taken, i.e `step_threshold < r` (with `r` defined in `shrink_threshold`), a check is
     also made to see if `expand_threshold < r`. If that is true, the trust region radius is
@@ -263,13 +234,6 @@ the value used in the respective paper.
     `shrink_threshold > r` (with `r` defined in `shrink_threshold`). Defaults to `0.25`.
   - `expand_factor`: the factor to expand the trust region radius with if
     `expand_threshold < r` (with `r` defined in `shrink_threshold`). Defaults to `2.0`.
-
-### References
-
-[1] Yuan, Ya-xiang. "Recent advances in trust region algorithms." Mathematical Programming
-151 (2015): 249-281.
-[2] Rahpeymaii, Farzad. "An efficient line search trust-region for systems of nonlinear
-equations." Mathematical Sciences 14.3 (2020): 257-268.
 """
 @kwdef @concrete struct GenericTrustRegionScheme{
     M <: RadiusUpdateSchemes.AbstractRadiusUpdateScheme}
@@ -358,7 +322,7 @@ end
         u0_norm, fu_norm) where {T}
     method isa RUS.__NLsolve && return T(ifelse(u0_norm > 0, u0_norm, 1))
     (method isa RUS.__Hei || method isa RUS.__Bastin) && return T(1)
-    method isa RUS.__Fan && return T((fu_norm^0.99) // 10)
+    method isa RUS.__Fan && return T((fu_norm^0.99) / 10)
     return T(max_tr / 11)
 end
 
diff --git a/src/internal/linear_solve.jl b/src/internal/linear_solve.jl
diff --git a/src/timer_outputs.jl b/src/timer_outputs.jl
diff --git a/test/misc/jacobian_reuse.jl b/test/misc/jacobian_reuse.jl

Original file line number	Diff line number	Diff line change
`@@ -1,3 +1,5 @@`
	`1`	`+# [Diagnostics API](@id diagnostics_api)`
	`2`	`+`
`1`	`3`	`# Logging the Solve Process`
`2`	`4`
`3`	`5`	`All NonlinearSolve.jl native solvers allow storing and displaying the trace of the nonlinear`