From adde8c0309b83bb4a0986d1f97c45847732e80eb Mon Sep 17 00:00:00 2001 From: Avik Pal Date: Mon, 11 Dec 2023 19:56:19 -0500 Subject: [PATCH 1/4] Use approximate sparsity detection by default --- Project.toml | 4 +++- docs/src/tutorials/large_systems.md | 32 +++++++++++++++++++++++------ ext/NonlinearSolveSymbolicsExt.jl | 7 +++++++ src/jacobian.jl | 10 ++++++++- 4 files changed, 45 insertions(+), 8 deletions(-) create mode 100644 ext/NonlinearSolveSymbolicsExt.jl diff --git a/Project.toml b/Project.toml index d334fd8ba..c7b6fa42a 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "NonlinearSolve" uuid = "8913a72c-1f9b-4ce2-8d82-65094dcecaec" authors = ["SciML"] -version = "3.0.1" +version = "3.0.2" [deps] ADTypes = "47edcb42-4c32-4615-8424-f2b9edc5f35b" @@ -35,6 +35,7 @@ FastLevenbergMarquardt = "7a0df574-e128-4d35-8cbd-3d84502bf7ce" LeastSquaresOptim = "0fc2ff8b-aaa3-5acd-a817-1944a5e08891" MINPACK = "4854310b-de5a-5eb6-a2a5-c1dee2bd17f9" NLsolve = "2774e3e8-f4cf-5e23-947b-6d7e65073b56" +Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7" Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f" [extensions] @@ -43,6 +44,7 @@ NonlinearSolveFastLevenbergMarquardtExt = "FastLevenbergMarquardt" NonlinearSolveLeastSquaresOptimExt = "LeastSquaresOptim" NonlinearSolveMINPACKExt = "MINPACK" NonlinearSolveNLsolveExt = "NLsolve" +NonlinearSolveSymbolicsExt = "Symbolics" NonlinearSolveZygoteExt = "Zygote" [compat] diff --git a/docs/src/tutorials/large_systems.md b/docs/src/tutorials/large_systems.md index 5748cb351..c209195a2 100644 --- a/docs/src/tutorials/large_systems.md +++ b/docs/src/tutorials/large_systems.md @@ -60,12 +60,14 @@ broadcast). Use `dx=dy=1/32`. The resulting `NonlinearProblem` definition is: ```@example ill_conditioned_nlprob -using NonlinearSolve, LinearAlgebra, SparseArrays, LinearSolve, Symbolics +using NonlinearSolve, LinearAlgebra, SparseArrays, LinearSolve const N = 32 const xyd_brusselator = range(0, stop = 1, length = N) + brusselator_f(x, y) = (((x - 0.3)^2 + (y - 0.6)^2) <= 0.1^2) * 5.0 limit(a, N) = a == N + 1 ? 1 : a == 0 ? N : a + function brusselator_2d_loop(du, u, p) A, B, alpha, dx = p alpha = alpha / dx^2 @@ -96,6 +98,7 @@ function init_brusselator_2d(xyd) end u end + u0 = init_brusselator_2d(xyd_brusselator) prob_brusselator_2d = NonlinearProblem(brusselator_2d_loop, u0, p) ``` @@ -120,6 +123,26 @@ include: - BandedMatrix ([BandedMatrices.jl](https://github.com/JuliaLinearAlgebra/BandedMatrices.jl)) - BlockBandedMatrix ([BlockBandedMatrices.jl](https://github.com/JuliaLinearAlgebra/BlockBandedMatrices.jl)) +## Approximate Sparsity Detection & Sparse Jacobians + +In the next section, we will discuss how to declare a sparse Jacobian and how to use +[Symbolics.jl](https://github.com/JuliaSymbolics/Symbolics.jl), to compute exact sparse +jacobians. For this to work it is important to not load `Symbolics.jl`. This is triggered +if you pass in a sparse autodiff type such as `AutoSparseForwardDiff()`. + +```@example ill_conditioned_nlprob +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())); +@btime solve(prob_brusselator_2d, + NewtonRaphson(; autodiff = AutoSparseForwardDiff(), + linsolve = KrylovJL_GMRES())); +``` + ## Declaring a Sparse Jacobian with Automatic Sparsity Detection Jacobian sparsity is declared by the `jac_prototype` argument in the `NonlinearFunction`. @@ -156,7 +179,6 @@ prob_brusselator_2d_sparse = NonlinearProblem(ff, u0, p) Now let's see how the version with sparsity compares to the version without: ```@example ill_conditioned_nlprob -using BenchmarkTools # for @btime @btime solve(prob_brusselator_2d, NewtonRaphson()); @btime solve(prob_brusselator_2d_sparse, NewtonRaphson()); @btime solve(prob_brusselator_2d_sparse, NewtonRaphson(linsolve = KLUFactorization())); @@ -258,10 +280,8 @@ function algebraicmultigrid2(W, du, u, p, t, newW, Plprev, Prprev, solverdata) if newW === nothing || newW A = convert(AbstractMatrix, W) Pl = AlgebraicMultigrid.aspreconditioner(AlgebraicMultigrid.ruge_stuben(A, - presmoother = AlgebraicMultigrid.Jacobi(rand(size(A, - 1))), - postsmoother = AlgebraicMultigrid.Jacobi(rand(size(A, - 1))))) + presmoother = AlgebraicMultigrid.Jacobi(rand(size(A, 1))), + postsmoother = AlgebraicMultigrid.Jacobi(rand(size(A, 1))))) else Pl = Plprev end diff --git a/ext/NonlinearSolveSymbolicsExt.jl b/ext/NonlinearSolveSymbolicsExt.jl new file mode 100644 index 000000000..95f1016e4 --- /dev/null +++ b/ext/NonlinearSolveSymbolicsExt.jl @@ -0,0 +1,7 @@ +module NonlinearSolveSymbolicsExt + +import NonlinearSolve, Symbolics + +NonlinearSolve.is_extension_loaded(::Val{:Symbolics}) = true + +end diff --git a/src/jacobian.jl b/src/jacobian.jl index a1870ffd1..bdd5adfd3 100644 --- a/src/jacobian.jl +++ b/src/jacobian.jl @@ -12,7 +12,15 @@ sparsity_detection_alg(_, _) = NoSparsityDetection() function sparsity_detection_alg(f, ad::AbstractSparseADType) if f.sparsity === nothing if f.jac_prototype === nothing - return SymbolicsSparsityDetection() + if is_extension_loaded(Val(:Symbolics)) + return SymbolicsSparsityDetection() + else + @warn "Symbolics.jl is not loaded and sparse AD mode $(ad) is being used. \ + Using approximate sparsity detection using ForwardDiff. This can \ + potentially fail or generate incorrect sparsity pattern for \ + complicated problems. Use with caution." maxlog=1 + return ApproximateJacobianSparsity() + end else jac_prototype = f.jac_prototype end From 1af69bac1f25ca37070250ab320a2edbc54ba734 Mon Sep 17 00:00:00 2001 From: Avik Pal Date: Tue, 12 Dec 2023 00:01:29 -0500 Subject: [PATCH 2/4] Update docs build --- docs/make.jl | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/docs/make.jl b/docs/make.jl index 959463b09..6a7923c38 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -10,11 +10,10 @@ include("pages.jl") makedocs(; sitename = "NonlinearSolve.jl", authors = "Chris Rackauckas", - modules = [NonlinearSolve, SciMLBase, DiffEqBase, SimpleNonlinearSolve, Sundials, - SteadyStateDiffEq], + modules = [NonlinearSolve, SimpleNonlinearSolve, SteadyStateDiffEq], clean = true, doctest = false, linkcheck = true, linkcheck_ignore = ["https://twitter.com/ChrisRackauckas/status/1544743542094020615"], - warnonly = [:missing_docs, :cross_references], + warnonly = [:cross_references], format = Documenter.HTML(assets = ["assets/favicon.ico"], canonical = "https://docs.sciml.ai/NonlinearSolve/stable/"), pages) From b5692448e6a7e2b6490fd2a0ed5a75ea1363d93e Mon Sep 17 00:00:00 2001 From: Avik Pal Date: Tue, 12 Dec 2023 00:50:34 -0500 Subject: [PATCH 3/4] Fix the docs to check for exports only --- docs/make.jl | 5 +++-- docs/src/api/nonlinearsolve.md | 3 +++ docs/src/api/simplenonlinearsolve.md | 1 + docs/src/tutorials/large_systems.md | 7 +++---- src/lbroyden.jl | 2 +- 5 files changed, 11 insertions(+), 7 deletions(-) diff --git a/docs/make.jl b/docs/make.jl index 6a7923c38..e096f58d7 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -10,10 +10,11 @@ include("pages.jl") makedocs(; sitename = "NonlinearSolve.jl", authors = "Chris Rackauckas", - modules = [NonlinearSolve, SimpleNonlinearSolve, SteadyStateDiffEq], + modules = [NonlinearSolve, SimpleNonlinearSolve, SteadyStateDiffEq, Sundials, + DiffEqBase, SciMLBase], clean = true, doctest = false, linkcheck = true, linkcheck_ignore = ["https://twitter.com/ChrisRackauckas/status/1544743542094020615"], - warnonly = [:cross_references], + warnonly = [:cross_references], checkdocs = :export, format = Documenter.HTML(assets = ["assets/favicon.ico"], canonical = "https://docs.sciml.ai/NonlinearSolve/stable/"), pages) diff --git a/docs/src/api/nonlinearsolve.md b/docs/src/api/nonlinearsolve.md index ad0fd3e14..cefda9ad7 100644 --- a/docs/src/api/nonlinearsolve.md +++ b/docs/src/api/nonlinearsolve.md @@ -10,6 +10,7 @@ PseudoTransient DFSane Broyden Klement +LimitedMemoryBroyden ``` ## Nonlinear Least Squares Solvers @@ -50,4 +51,6 @@ RadiusUpdateSchemes.Hei RadiusUpdateSchemes.Yuan RadiusUpdateSchemes.Bastin RadiusUpdateSchemes.Fan +RadiusUpdateSchemes.NLsolve +RadiusUpdateSchemes.NocedalWright ``` diff --git a/docs/src/api/simplenonlinearsolve.md b/docs/src/api/simplenonlinearsolve.md index 72f738140..f10fb78d6 100644 --- a/docs/src/api/simplenonlinearsolve.md +++ b/docs/src/api/simplenonlinearsolve.md @@ -11,6 +11,7 @@ i.e. `IntervalNonlinearProblem`. ```@docs ITP +Alefeld Bisection Falsi Ridder diff --git a/docs/src/tutorials/large_systems.md b/docs/src/tutorials/large_systems.md index c209195a2..903310c21 100644 --- a/docs/src/tutorials/large_systems.md +++ b/docs/src/tutorials/large_systems.md @@ -136,11 +136,10 @@ 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(), linsolve = KLUFactorization())); @btime solve(prob_brusselator_2d, - NewtonRaphson(; autodiff = AutoSparseForwardDiff(), - linsolve = KrylovJL_GMRES())); + NewtonRaphson(; autodiff = AutoSparseForwardDiff(), linsolve = KrylovJL_GMRES())); +nothing # hide ``` ## Declaring a Sparse Jacobian with Automatic Sparsity Detection diff --git a/src/lbroyden.jl b/src/lbroyden.jl index f31be0d46..811e3400d 100644 --- a/src/lbroyden.jl +++ b/src/lbroyden.jl @@ -14,7 +14,7 @@ An implementation of `LimitedMemoryBroyden` with resetting and line search. - `linesearch`: the line search algorithm to use. Defaults to [`LineSearch()`](@ref), which means that no line search is performed. Algorithms from `LineSearches.jl` can be used here directly, and they will be converted to the correct `LineSearch`. It is - recommended to use [`LiFukushimaLineSearchCache`](@ref) -- a derivative free linesearch + recommended to use [`LiFukushimaLineSearch`](@ref) -- a derivative free linesearch specifically designed for Broyden's method. """ @concrete struct LimitedMemoryBroyden{threshold} <: AbstractNewtonAlgorithm{false, Nothing} From 540ef5b4a85b3e9cb9cebadc5ea57526970f0366 Mon Sep 17 00:00:00 2001 From: Avik Pal Date: Tue, 12 Dec 2023 12:21:13 -0500 Subject: [PATCH 4/4] Handle Approximate Sparsity detection more gracefully --- docs/pages.jl | 1 + docs/src/basics/SparsityDetection.md | 77 ++++++++++++++++++++++++++++ docs/src/tutorials/large_systems.md | 30 +++++++++-- src/jacobian.jl | 20 ++++++-- src/utils.jl | 1 + 5 files changed, 120 insertions(+), 9 deletions(-) create mode 100644 docs/src/basics/SparsityDetection.md diff --git a/docs/pages.jl b/docs/pages.jl index 538d0b098..a8107bea2 100644 --- a/docs/pages.jl +++ b/docs/pages.jl @@ -14,6 +14,7 @@ pages = ["index.md", "basics/NonlinearSolution.md", "basics/TerminationCondition.md", "basics/Logging.md", + "basics/SparsityDetection.md", "basics/FAQ.md"], "Solver Summaries and Recommendations" => Any["solvers/NonlinearSystemSolvers.md", "solvers/BracketingSolvers.md", diff --git a/docs/src/basics/SparsityDetection.md b/docs/src/basics/SparsityDetection.md new file mode 100644 index 000000000..000c9905f --- /dev/null +++ b/docs/src/basics/SparsityDetection.md @@ -0,0 +1,77 @@ +# [(Semi-)Automatic Sparsity Detection](@id sparsity-detection) + +This section describes how to enable Sparsity Detection. For a detailed tutorial on how +to use this in an actual problem, see +[this tutorial on Efficiently Solving Large Sparse Ill-Conditioned Nonlinear Systems](@ref large_systems). + +Notation wise we are trying to solve for `x` such that `nlfunc(x) = 0`. + +## Case I: Sparse Jacobian Prototype is Provided + +Let's say you have a Sparse Jacobian Prototype `jac_prototype`, in this case you can +create your problem as follows: + +```julia +prob = NonlinearProblem(NonlinearFunction(nlfunc; sparsity = jac_prototype), x0) +# OR +prob = NonlinearProblem(NonlinearFunction(nlfunc; jac_prototype = jac_prototype), x0) +``` + +NonlinearSolve will automatically perform matrix coloring and use sparse differentiation. + +Now you can help the solver further by providing the color vector. This can be done by + +```julia +prob = NonlinearProblem(NonlinearFunction(nlfunc; sparsity = jac_prototype, + colorvec = colorvec), x0) +# OR +prob = NonlinearProblem(NonlinearFunction(nlfunc; jac_prototype = jac_prototype, + colorvec = colorvec), x0) +``` + +!!! note + + One thing to be careful about in this case is that `colorvec` is dependent on the + autodiff backend used. Forward Mode and Finite Differencing will assume that the + colorvec is the column colorvec, while Reverse Mode will assume that the colorvec is the + row colorvec. + +## Case II: Sparsity Detection algorithm is provided + +If you don't have a Sparse Jacobian Prototype, but you know the which sparsity detection +algorithm you want to use, then you can create your problem as follows: + +```julia +prob = NonlinearProblem(NonlinearFunction(nlfunc; sparsity = SymbolicsSparsityDetection()), + x0) # Remember to have Symbolics.jl loaded +# OR +prob = NonlinearProblem(NonlinearFunction(nlfunc; sparsity = ApproximateJacobianSparsity()), + x0) +``` + +These Detection Algorithms are from [SparseDiffTools.jl](https://github.com/JuliaDiff/SparseDiffTools.jl), +refer to the documentation there for more details. + +## Case III: Sparse AD Type is being Used + +If you constructed a Nonlinear Solver with a sparse AD type, for example + +```julia +NewtonRaphson(; autodiff = AutoSparseForwardDiff()) +# OR +TrustRegion(; autodiff = AutoSparseZygote()) +``` + +then NonlinearSolve will automatically perform matrix coloring and use sparse +differentiation if none of `sparsity` or `jac_prototype` is provided. If `Symbolics.jl` is +loaded, we default to using `SymbolicsSparsityDetection()`, else we default to using +`ApproximateJacobianSparsity()`. + +`Case I/II` take precedence for sparsity detection and we perform sparse AD based on those +options if those are provided. + +!!! warning + + If you provide a non-sparse AD, and provide a `sparsity` or `jac_prototype` then + we will use dense AD. This is because, if you provide a specific AD type, we assume + that you know what you are doing and want to override the default choice of `nothing`. diff --git a/docs/src/tutorials/large_systems.md b/docs/src/tutorials/large_systems.md index 903310c21..1eab0d88b 100644 --- a/docs/src/tutorials/large_systems.md +++ b/docs/src/tutorials/large_systems.md @@ -60,7 +60,7 @@ broadcast). Use `dx=dy=1/32`. The resulting `NonlinearProblem` definition is: ```@example ill_conditioned_nlprob -using NonlinearSolve, LinearAlgebra, SparseArrays, LinearSolve +using NonlinearSolve, LinearAlgebra, SparseArrays, LinearSolve, SparseDiffTools const N = 32 const xyd_brusselator = range(0, stop = 1, length = N) @@ -100,7 +100,8 @@ function init_brusselator_2d(xyd) end u0 = init_brusselator_2d(xyd_brusselator) -prob_brusselator_2d = NonlinearProblem(brusselator_2d_loop, u0, p) +prob_brusselator_2d = NonlinearProblem(brusselator_2d_loop, u0, p; abstol = 1e-10, + reltol = 1e-10) ``` ## Choosing Jacobian Types @@ -127,8 +128,10 @@ include: In the next section, we will discuss how to declare a sparse Jacobian and how to use [Symbolics.jl](https://github.com/JuliaSymbolics/Symbolics.jl), to compute exact sparse -jacobians. For this to work it is important to not load `Symbolics.jl`. This is triggered -if you pass in a sparse autodiff type such as `AutoSparseForwardDiff()`. +jacobians. This is triggered if you pass in a sparse autodiff type such as +`AutoSparseForwardDiff()`. If `Symbolics.jl` is loaded, then the default changes to +Symbolic Sparsity Detection. See the manual entry on +[Sparsity Detection](@ref sparsity-detection) for more details on the default. ```@example ill_conditioned_nlprob using BenchmarkTools # for @btime @@ -223,6 +226,7 @@ used in the solution of the ODE. An example of this with using ```@example ill_conditioned_nlprob using IncompleteLU + function incompletelu(W, du, u, p, t, newW, Plprev, Prprev, solverdata) if newW === nothing || newW Pl = ilu(W, τ = 50.0) @@ -257,6 +261,7 @@ which is more automatic. The setup is very similar to before: ```@example ill_conditioned_nlprob using AlgebraicMultigrid + function algebraicmultigrid(W, du, u, p, t, newW, Plprev, Prprev, solverdata) if newW === nothing || newW Pl = aspreconditioner(ruge_stuben(convert(AbstractMatrix, W))) @@ -293,5 +298,22 @@ end nothing # hide ``` +## Let's compare the Sparsity Detection Methods + +We benchmarked the solvers before with approximate and exact sparsity detection. However, +for the exact sparsity detection case, we left out the time it takes to perform exact +sparsity detection. Let's compare the two by setting the sparsity detection algorithms. + +```@example ill_conditioned_nlprob +prob_brusselator_2d_exact = NonlinearProblem(NonlinearFunction(brusselator_2d_loop; + sparsity = SymbolicsSparsityDetection()), u0, p; abstol = 1e-10, reltol = 1e-10) +prob_brusselator_2d_approx = NonlinearProblem(NonlinearFunction(brusselator_2d_loop; + sparsity = ApproximateJacobianSparsity()), u0, p; abstol = 1e-10, reltol = 1e-10) + +@btime solve(prob_brusselator_2d_exact, NewtonRaphson()); +@btime solve(prob_brusselator_2d_approx, NewtonRaphson()); +nothing # hide +``` + For more information on the preconditioner interface, see the [linear solver documentation](https://docs.sciml.ai/LinearSolve/stable/basics/Preconditioners/). diff --git a/src/jacobian.jl b/src/jacobian.jl index bdd5adfd3..4bb9d5571 100644 --- a/src/jacobian.jl +++ b/src/jacobian.jl @@ -15,17 +15,27 @@ function sparsity_detection_alg(f, ad::AbstractSparseADType) if is_extension_loaded(Val(:Symbolics)) return SymbolicsSparsityDetection() else - @warn "Symbolics.jl is not loaded and sparse AD mode $(ad) is being used. \ - Using approximate sparsity detection using ForwardDiff. This can \ - potentially fail or generate incorrect sparsity pattern for \ - complicated problems. Use with caution." maxlog=1 return ApproximateJacobianSparsity() end else jac_prototype = f.jac_prototype end - else + elseif f.sparsity isa SparseDiffTools.AbstractSparsityDetection + if f.jac_prototype === nothing + return f.sparsity + else + jac_prototype = f.jac_prototype + end + elseif f.sparsity isa AbstractMatrix jac_prototype = f.sparsity + elseif f.jac_prototype isa AbstractMatrix + jac_prototype = f.jac_prototype + else + error("`sparsity::typeof($(typeof(f.sparsity)))` & \ + `jac_prototype::typeof($(typeof(f.jac_prototype)))` is not supported. \ + Use `sparsity::AbstractMatrix` or `sparsity::AbstractSparsityDetection` or \ + set to `nothing`. `jac_prototype` can be set to `nothing` or an \ + `AbstractMatrix`.") end if SciMLBase.has_colorvec(f) diff --git a/src/utils.jl b/src/utils.jl index 075a1857a..45baf7aba 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -384,6 +384,7 @@ LazyArrays.applied_axes(::typeof(__zero), x) = axes(x) end return pinv(A) end +@inline __safe_inv(A::SparseMatrixCSC) = __safe_inv(Matrix(A)) LazyArrays.applied_eltype(::typeof(__safe_inv), x) = eltype(x) LazyArrays.applied_ndims(::typeof(__safe_inv), x) = ndims(x)