From 95695357bc52f43156b95603d1f6c194d3de3ea2 Mon Sep 17 00:00:00 2001 From: Avik Pal Date: Mon, 2 Oct 2023 20:42:06 -0400 Subject: [PATCH 1/8] Fix bug in Levenberg --- Project.toml | 2 +- src/levenberg.jl | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/Project.toml b/Project.toml index c7decafdb..2d9474b20 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "NonlinearSolve" uuid = "8913a72c-1f9b-4ce2-8d82-65094dcecaec" authors = ["SciML"] -version = "2.0.1" +version = "2.0.2" [deps] ADTypes = "47edcb42-4c32-4615-8424-f2b9edc5f35b" diff --git a/src/levenberg.jl b/src/levenberg.jl index 7264c127f..30665178f 100644 --- a/src/levenberg.jl +++ b/src/levenberg.jl @@ -216,7 +216,7 @@ function perform_step!(cache::LevenbergMarquardtCache{true}) # The following lines do: cache.a = -J \ cache.fu_tmp mul!(cache.du, J, v) @. cache.fu_tmp = (2 / h) * ((cache.fu_tmp - fu1) / h - cache.du) - linres = dolinsolve(alg.precs, linsolve; A = J, b = _vec(cache.fu_tmp), + linres = dolinsolve(alg.precs, linsolve; A = cache.mat_tmp, b = _vec(cache.fu_tmp), linu = _vec(cache.du), p = p, reltol = cache.abstol) cache.linsolve = linres.cache @. cache.a = -cache.du @@ -225,7 +225,7 @@ function perform_step!(cache::LevenbergMarquardtCache{true}) # Require acceptable steps to satisfy the following condition. norm_v = norm(v) - if (2 * norm(cache.a) / norm_v) < α_geodesic + if 1 + log2(norm(cache.a)) - log2(norm_v) ≤ log2(α_geodesic) @. cache.δ = v + cache.a / 2 @unpack δ, loss_old, norm_v_old, v_old, b_uphill = cache f(cache.fu_tmp, u .+ δ, p) From 0b4b5b9262907bf184183efc52fa9ab32e00f3f5 Mon Sep 17 00:00:00 2001 From: Avik Pal Date: Tue, 3 Oct 2023 10:14:53 -0400 Subject: [PATCH 2/8] Use `@fastmath` --- src/levenberg.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/levenberg.jl b/src/levenberg.jl index 30665178f..c71c9016c 100644 --- a/src/levenberg.jl +++ b/src/levenberg.jl @@ -225,7 +225,7 @@ function perform_step!(cache::LevenbergMarquardtCache{true}) # Require acceptable steps to satisfy the following condition. norm_v = norm(v) - if 1 + log2(norm(cache.a)) - log2(norm_v) ≤ log2(α_geodesic) + if @fastmath((1 + log2(norm(cache.a)) - log2(norm_v))≤log2(α_geodesic)) @. cache.δ = v + cache.a / 2 @unpack δ, loss_old, norm_v_old, v_old, b_uphill = cache f(cache.fu_tmp, u .+ δ, p) From 4bdf983246a3c4cb6d76bda2ce409044b9f25a04 Mon Sep 17 00:00:00 2001 From: Avik Pal Date: Tue, 3 Oct 2023 10:40:23 -0400 Subject: [PATCH 3/8] Patch the OOP version as well --- src/levenberg.jl | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/src/levenberg.jl b/src/levenberg.jl index c71c9016c..d4eb23a3f 100644 --- a/src/levenberg.jl +++ b/src/levenberg.jl @@ -274,18 +274,19 @@ function perform_step!(cache::LevenbergMarquardtCache{false}) end @unpack u, p, λ, JᵀJ, DᵀD, J = cache + cache.mat_tmp = JᵀJ + λ * DᵀD # Usual Levenberg-Marquardt step ("velocity"). - cache.v = -(JᵀJ + λ * DᵀD) \ (J' * fu1) + cache.v = -cache.mat_tmp \ (J' * fu1) @unpack v, h, α_geodesic = cache # Geodesic acceleration (step_size = v + a / 2). - cache.a = -J \ ((2 / h) .* ((f(u .+ h .* v, p) .- fu1) ./ h .- J * v)) + cache.a = -cache.mat_tmp \ ((2 / h) .* ((f(u .+ h .* v, p) .- fu1) ./ h .- J * v)) cache.stats.nsolve += 1 cache.stats.nfactors += 1 # Require acceptable steps to satisfy the following condition. norm_v = norm(v) - if (2 * norm(cache.a) / norm_v) < α_geodesic + if @fastmath((1 + log2(norm(cache.a)) - log2(norm_v))≤log2(α_geodesic)) cache.δ = v .+ cache.a ./ 2 @unpack δ, loss_old, norm_v_old, v_old, b_uphill = cache fu_new = f(u .+ δ, p) From c17ccbf352d752da411ee376a440f1bebff88731 Mon Sep 17 00:00:00 2001 From: Avik Pal Date: Thu, 5 Oct 2023 19:38:08 -0400 Subject: [PATCH 4/8] Revert log2 change --- src/levenberg.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/levenberg.jl b/src/levenberg.jl index d4eb23a3f..2a21f9936 100644 --- a/src/levenberg.jl +++ b/src/levenberg.jl @@ -225,7 +225,7 @@ function perform_step!(cache::LevenbergMarquardtCache{true}) # Require acceptable steps to satisfy the following condition. norm_v = norm(v) - if @fastmath((1 + log2(norm(cache.a)) - log2(norm_v))≤log2(α_geodesic)) + if 2 * norm(cache.a) ≤ α_geodesic * norm_v @. cache.δ = v + cache.a / 2 @unpack δ, loss_old, norm_v_old, v_old, b_uphill = cache f(cache.fu_tmp, u .+ δ, p) @@ -286,7 +286,7 @@ function perform_step!(cache::LevenbergMarquardtCache{false}) # Require acceptable steps to satisfy the following condition. norm_v = norm(v) - if @fastmath((1 + log2(norm(cache.a)) - log2(norm_v))≤log2(α_geodesic)) + if 2 * norm(cache.a) ≤ α_geodesic * norm_v cache.δ = v .+ cache.a ./ 2 @unpack δ, loss_old, norm_v_old, v_old, b_uphill = cache fu_new = f(u .+ δ, p) From 9d98446416b9c60f8c784245e8658b7836004645 Mon Sep 17 00:00:00 2001 From: Avik Pal Date: Thu, 5 Oct 2023 19:59:16 -0400 Subject: [PATCH 5/8] Add 23 test problems test Co-authored-by: Flemming Holtorf --- Project.toml | 4 ++- test/23_test_problems.jl | 62 ++++++++++++++++++++++++++++++++++++++++ test/runtests.jl | 2 ++ 3 files changed, 67 insertions(+), 1 deletion(-) create mode 100644 test/23_test_problems.jl diff --git a/Project.toml b/Project.toml index 2d9474b20..e521cc75b 100644 --- a/Project.toml +++ b/Project.toml @@ -35,6 +35,7 @@ FiniteDiff = "2" ForwardDiff = "0.10.3" LineSearches = "7" LinearSolve = "2" +NonlinearProblemLibrary = "0.1" PrecompileTools = "1" RecursiveArrayTools = "2" Reexport = "0.2, 1" @@ -52,6 +53,7 @@ Enzyme = "7da242da-08ed-463a-9acd-ee780be4f1d9" ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" LinearSolve = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae" +NonlinearProblemLibrary = "b7050fa9-e91f-4b37-bcee-a89a063da141" Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f" @@ -62,4 +64,4 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f" [targets] -test = ["Enzyme", "BenchmarkTools", "SafeTestsets", "Pkg", "Test", "ForwardDiff", "StaticArrays", "Symbolics", "LinearSolve", "Random", "LinearAlgebra", "Zygote", "SparseDiffTools"] +test = ["Enzyme", "BenchmarkTools", "SafeTestsets", "Pkg", "Test", "ForwardDiff", "StaticArrays", "Symbolics", "LinearSolve", "Random", "LinearAlgebra", "Zygote", "SparseDiffTools", "NonlinearProblemLibrary"] diff --git a/test/23_test_problems.jl b/test/23_test_problems.jl new file mode 100644 index 000000000..f13182119 --- /dev/null +++ b/test/23_test_problems.jl @@ -0,0 +1,62 @@ +using NonlinearSolve, LinearAlgebra, NonlinearProblemLibrary, Test + +problems = NonlinearProblemLibrary.problems +dicts = NonlinearProblemLibrary.dicts + +function test_on_library(problems, dicts, alg_ops, broken_tests, ϵ = 1e-5) + for (idx, (problem, dict)) in enumerate(zip(problems, dicts)) + x = dict["start"] + res = similar(x) + nlprob = NonlinearProblem(problem, x) + @testset "$(dict["title"])" begin + for alg in alg_ops + sol = solve(nlprob, alg, abstol = 1e-15, reltol = 1e-15) + problem(res, sol.u, nothing) + broken = idx in broken_tests[alg] ? true : false + @test norm(res)≤ϵ broken=broken + end + end + end +end + +# NewtonRaphson +@testset "NewtonRaphson test problem library" begin + alg_ops = (NewtonRaphson(),) + + # dictionary with indices of test problems where method does not converge to small residual + broken_tests = Dict(alg => Int[] for alg in alg_ops) + broken_tests[alg_ops[1]] = [1, 6] + + test_on_library(problems, dicts, alg_ops, broken_tests) +end + +@testset "TrustRegion test problem library" begin + alg_ops = (TrustRegion(; radius_update_scheme = RadiusUpdateSchemes.Simple), + TrustRegion(; radius_update_scheme = RadiusUpdateSchemes.Fan), + TrustRegion(; radius_update_scheme = RadiusUpdateSchemes.Hei), + TrustRegion(; radius_update_scheme = RadiusUpdateSchemes.Yuan), + TrustRegion(; radius_update_scheme = RadiusUpdateSchemes.Bastin), + TrustRegion(; radius_update_scheme = RadiusUpdateSchemes.NLsolve)) + + # dictionary with indices of test problems where method does not converge to small residual + broken_tests = Dict(alg => Int[] for alg in alg_ops) + broken_tests[alg_ops[1]] = [6, 11, 21] + broken_tests[alg_ops[2]] = [6, 11, 21] + broken_tests[alg_ops[3]] = [1, 6, 11, 12, 15, 16, 21] + broken_tests[alg_ops[4]] = [1, 6, 8, 11, 15, 16, 21, 22] + broken_tests[alg_ops[5]] = [6, 21] + broken_tests[alg_ops[6]] = [6, 21] + + test_on_library(problems, dicts, alg_ops, broken_tests) +end + +@testset "TrustRegion test problem library" begin + alg_ops = (LevenbergMarquardt(), LevenbergMarquardt(; α_geodesic = 0.5)) + + # dictionary with indices of test problems where method does not converge to small residual + broken_tests = Dict(alg => Int[] for alg in alg_ops) + broken_tests[alg_ops[1]] = [3, 6, 11, 17, 21] + broken_tests[alg_ops[2]] = [3, 6, 11, 21] + + test_on_library(problems, dicts, alg_ops, broken_tests) +end diff --git a/test/runtests.jl b/test/runtests.jl index a84fc3cb1..58ef647f0 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -15,6 +15,8 @@ end if GROUP == "All" || GROUP == "Core" @time @safetestset "Basic Tests + Some AD" include("basictests.jl") @time @safetestset "Sparsity Tests" include("sparse.jl") + + @time @safetestset "23 Test Problems" include("23_test_problems.jl") end if GROUP == "GPU" From a1877cd83b73eb9911c1f7edf8f7768106c62c49 Mon Sep 17 00:00:00 2001 From: Avik Pal Date: Thu, 5 Oct 2023 20:31:03 -0400 Subject: [PATCH 6/8] Fix Jacobian Construction for concrete_jac --- Project.toml | 2 +- src/jacobian.jl | 11 ++++++++--- 2 files changed, 9 insertions(+), 4 deletions(-) diff --git a/Project.toml b/Project.toml index e521cc75b..a07d1b735 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "NonlinearSolve" uuid = "8913a72c-1f9b-4ce2-8d82-65094dcecaec" authors = ["SciML"] -version = "2.0.2" +version = "2.0.3" [deps] ADTypes = "47edcb42-4c32-4615-8424-f2b9edc5f35b" diff --git a/src/jacobian.jl b/src/jacobian.jl index 01c030462..82f2ef2bb 100644 --- a/src/jacobian.jl +++ b/src/jacobian.jl @@ -1,4 +1,4 @@ -@concrete struct JacobianWrapper{iip} +@concrete struct JacobianWrapper{iip} <: Function f p end @@ -73,9 +73,9 @@ function jacobian_caches(alg::AbstractNonlinearSolveAlgorithm, f, u, p, ::Val{ii jac_cache = nothing end - J = if !linsolve_needs_jac + J = if !(linsolve_needs_jac || alg_wants_jac) # We don't need to construct the Jacobian - JacVec(uf, u; autodiff = alg.ad) + JacVec(uf, u; autodiff = __get_nonsparse_ad(alg.ad)) else if has_analytic_jac f.jac_prototype === nothing ? undefmatrix(u) : f.jac_prototype @@ -98,6 +98,11 @@ function jacobian_caches(alg::AbstractNonlinearSolveAlgorithm, f, u, p, ::Val{ii return uf, linsolve, J, fu, jac_cache, du end +__get_nonsparse_ad(::AutoSparseForwardDiff) = AutoForwardDiff() +__get_nonsparse_ad(::AutoSparseFiniteDiff) = AutoFiniteDiff() +__get_nonsparse_ad(::AutoSparseZygote) = AutoZygote() +__get_nonsparse_ad(ad) = ad + ## Special Handling for Scalars function jacobian_caches(alg::AbstractNonlinearSolveAlgorithm, f, u::Number, p, ::Val{false}; kwargs...) From cb105fcbbbd3f172d7d6ce2620dc751b55d34cee Mon Sep 17 00:00:00 2001 From: Avik Pal Date: Thu, 5 Oct 2023 21:00:39 -0400 Subject: [PATCH 7/8] Fix Jacobian Construction --- src/NonlinearSolve.jl | 2 +- src/jacobian.jl | 26 +++++++++++++++++++++++++- 2 files changed, 26 insertions(+), 2 deletions(-) diff --git a/src/NonlinearSolve.jl b/src/NonlinearSolve.jl index 90839dea5..168658bfe 100644 --- a/src/NonlinearSolve.jl +++ b/src/NonlinearSolve.jl @@ -4,7 +4,7 @@ if isdefined(Base, :Experimental) && isdefined(Base.Experimental, Symbol("@max_m @eval Base.Experimental.@max_methods 1 end -using DiffEqBase, LinearAlgebra, LinearSolve, SparseDiffTools +using DiffEqBase, LinearAlgebra, LinearSolve, SparseArrays, SparseDiffTools import ForwardDiff import ADTypes: AbstractFiniteDifferencesMode diff --git a/src/jacobian.jl b/src/jacobian.jl index 82f2ef2bb..d9327e701 100644 --- a/src/jacobian.jl +++ b/src/jacobian.jl @@ -80,7 +80,11 @@ function jacobian_caches(alg::AbstractNonlinearSolveAlgorithm, f, u, p, ::Val{ii if has_analytic_jac f.jac_prototype === nothing ? undefmatrix(u) : f.jac_prototype else - f.jac_prototype === nothing ? init_jacobian(jac_cache) : f.jac_prototype + if f.jac_prototype === nothing + __safe_init_jacobian(jac_cache) + else + f.jac_prototype + end end end @@ -98,6 +102,26 @@ function jacobian_caches(alg::AbstractNonlinearSolveAlgorithm, f, u, p, ::Val{ii return uf, linsolve, J, fu, jac_cache, du end +@generated function __getfield(c::T, ::Val{S}) where {T, S} + hasfield(T, S) && return :(c.$(S)) + return :(nothing) +end + +function __safe_init_jacobian(c::SparseDiffTools.AbstractMaybeSparseJacobianCache) + T = promote_type(eltype(c.fx), eltype(c.x)) + return __safe_init_jacobian(__getfield(c, Val(:jac_prototype)), T, c.fx, c.x) +end +function __safe_init_jacobian(::Nothing, ::Type{T}, fx, x) where {T} + return similar(fx, T, length(fx), length(x)) +end +function __safe_init_jacobian(J::SparseMatrixCSC, ::Type{T}, fx, x) where {T} + @assert size(J, 1) == length(fx) && size(J, 2) == length(x) + return T.(J) +end +function __safe_init_jacobian(J, ::Type{T}, fx, x) where {T} + return similar(fx, T, length(fx), length(x)) # This is not safe for sparse jacobians +end + __get_nonsparse_ad(::AutoSparseForwardDiff) = AutoForwardDiff() __get_nonsparse_ad(::AutoSparseFiniteDiff) = AutoFiniteDiff() __get_nonsparse_ad(::AutoSparseZygote) = AutoZygote() From b5b929824b2bc922ecc828357969ffffbe0284fa Mon Sep 17 00:00:00 2001 From: Avik Pal Date: Thu, 5 Oct 2023 21:22:01 -0400 Subject: [PATCH 8/8] upstream fix to SparseDiffTools.jl --- src/jacobian.jl | 26 +------------------------- 1 file changed, 1 insertion(+), 25 deletions(-) diff --git a/src/jacobian.jl b/src/jacobian.jl index d9327e701..82f2ef2bb 100644 --- a/src/jacobian.jl +++ b/src/jacobian.jl @@ -80,11 +80,7 @@ function jacobian_caches(alg::AbstractNonlinearSolveAlgorithm, f, u, p, ::Val{ii if has_analytic_jac f.jac_prototype === nothing ? undefmatrix(u) : f.jac_prototype else - if f.jac_prototype === nothing - __safe_init_jacobian(jac_cache) - else - f.jac_prototype - end + f.jac_prototype === nothing ? init_jacobian(jac_cache) : f.jac_prototype end end @@ -102,26 +98,6 @@ function jacobian_caches(alg::AbstractNonlinearSolveAlgorithm, f, u, p, ::Val{ii return uf, linsolve, J, fu, jac_cache, du end -@generated function __getfield(c::T, ::Val{S}) where {T, S} - hasfield(T, S) && return :(c.$(S)) - return :(nothing) -end - -function __safe_init_jacobian(c::SparseDiffTools.AbstractMaybeSparseJacobianCache) - T = promote_type(eltype(c.fx), eltype(c.x)) - return __safe_init_jacobian(__getfield(c, Val(:jac_prototype)), T, c.fx, c.x) -end -function __safe_init_jacobian(::Nothing, ::Type{T}, fx, x) where {T} - return similar(fx, T, length(fx), length(x)) -end -function __safe_init_jacobian(J::SparseMatrixCSC, ::Type{T}, fx, x) where {T} - @assert size(J, 1) == length(fx) && size(J, 2) == length(x) - return T.(J) -end -function __safe_init_jacobian(J, ::Type{T}, fx, x) where {T} - return similar(fx, T, length(fx), length(x)) # This is not safe for sparse jacobians -end - __get_nonsparse_ad(::AutoSparseForwardDiff) = AutoForwardDiff() __get_nonsparse_ad(::AutoSparseFiniteDiff) = AutoFiniteDiff() __get_nonsparse_ad(::AutoSparseZygote) = AutoZygote()