From 468a2048ffa0202640f6e985f26b71840e6e2fc3 Mon Sep 17 00:00:00 2001 From: Simon Carlson Date: Sat, 31 Dec 2022 06:24:29 +0100 Subject: [PATCH 1/7] inital commit, started to implement Broyden --- src/SimpleNonlinearSolve.jl | 6 ++++-- src/broyden.jl | 43 +++++++++++++++++++++++++++++++++++++ test/basictests.jl | 4 +++- 3 files changed, 50 insertions(+), 3 deletions(-) create mode 100644 src/broyden.jl diff --git a/src/SimpleNonlinearSolve.jl b/src/SimpleNonlinearSolve.jl index bd47881..213ba23 100644 --- a/src/SimpleNonlinearSolve.jl +++ b/src/SimpleNonlinearSolve.jl @@ -4,6 +4,7 @@ using Reexport using FiniteDiff, ForwardDiff using ForwardDiff: Dual using StaticArraysCore +using LinearAlgebra # TODO check if it is ok to add this import ArrayInterfaceCore @reexport using SciMLBase @@ -18,12 +19,13 @@ include("bisection.jl") include("falsi.jl") include("raphson.jl") include("ad.jl") +include("broyden.jl") import SnoopPrecompile SnoopPrecompile.@precompile_all_calls begin for T in (Float32, Float64) prob_no_brack = NonlinearProblem{false}((u, p) -> u .* u .- p, T(0.1), T(2)) - for alg in (SimpleNewtonRaphson,) + for alg in (SimpleNewtonRaphson, Broyden) solve(prob_no_brack, alg(), tol = T(1e-2)) end @@ -44,6 +46,6 @@ SnoopPrecompile.@precompile_all_calls begin for T in (Float32, Float64) end end # DiffEq styled algorithms -export Bisection, Falsi, SimpleNewtonRaphson +export Bisection, Broyden, Falsi, SimpleNewtonRaphson end # module diff --git a/src/broyden.jl b/src/broyden.jl new file mode 100644 index 0000000..9f9388b --- /dev/null +++ b/src/broyden.jl @@ -0,0 +1,43 @@ +# TODO add docstrings + +# TODO check what the supertype should be +# TODO check if this should be defined as in raphson.jl +struct Broyden <: AbstractSimpleNonlinearSolveAlgorithm end + +function SciMLBase.solve(prob::NonlinearProblem, + alg::Broyden, args...; abstol = nothing, + reltol = nothing, + maxiters = 1000, kwargs...) + f = Base.Fix2(prob.f, prob.p) + xₙ = float(prob.u0) + T = typeof(xₙ) + J⁻¹ = Matrix{T}(I, length(xₙ), length(xₙ)) + + if SciMLBase.isinplace(prob) + error("Broyden currently only supports out-of-place nonlinear problems") + end + + atol = abstol !== nothing ? abstol : + real(oneunit(eltype(T))) * (eps(real(one(eltype(T)))))^(4 // 5) + rtol = reltol !== nothing ? reltol : eps(real(one(eltype(T))))^(4 // 5) + + for _ in 1:maxiters + # TODO check if nameing with heavy use of subscrips is ok + fₙ = f(xₙ) + xₙ₊₁ = xₙ + J⁻¹ * fₙ + Δxₙ = xₙ₊₁ - xₙ + Δfₙ = f(xₙ₊₁) - fₙ + J⁻¹ .+= ((Δxₙ .- J⁻¹ * Δfₙ) ./ (Δxₙ' * J⁻¹ * Δfₙ)) * (Δxₙ' * J⁻¹) + + iszero(fₙ) && + return SciMLBase.build_solution(prob, alg, xₙ₊₁, fₙ; retcode = ReturnCode.Success) + + if isapprox(xₙ₊₁, xₙ, atol = atol, rtol = rtol) + return SciMLBase.build_solution(prob, alg, xₙ₊₁, fₙ; retcode = ReturnCode.Success) + end + xₙ = xₙ₊₁ + end + + @show xₙ, fₙ + return SciMLBase.build_solution(prob, alg, xₙ, fₙ; retcode = ReturnCode.MaxIters) +end diff --git a/test/basictests.jl b/test/basictests.jl index 0c41946..bc5a390 100644 --- a/test/basictests.jl +++ b/test/basictests.jl @@ -77,7 +77,7 @@ for alg in [Bisection(), Falsi()] global g, p g = function (p) probN = IntervalNonlinearProblem{false}(f, tspan, p) - sol = solve(probN, Bisection()) + sol = solve(probN, Bisection()) # TODO check if "alg" should replace "Bisection()" return [sol.left] end @@ -102,6 +102,7 @@ probN = NonlinearProblem(f, u0) @test solve(probN, SimpleNewtonRaphson(); immutable = false).u[end] ≈ sqrt(2.0) @test solve(probN, SimpleNewtonRaphson(; autodiff = false)).u[end] ≈ sqrt(2.0) @test solve(probN, SimpleNewtonRaphson(; autodiff = false)).u[end] ≈ sqrt(2.0) +# TODO check why the 2 lines above are identical for u0 in [1.0, [1, 1.0]] local f, probN, sol @@ -109,6 +110,7 @@ for u0 in [1.0, [1, 1.0]] probN = NonlinearProblem(f, u0) sol = sqrt(2) * u0 + # TODO check why the two lines below are identical @test solve(probN, SimpleNewtonRaphson()).u ≈ sol @test solve(probN, SimpleNewtonRaphson()).u ≈ sol @test solve(probN, SimpleNewtonRaphson(; autodiff = false)).u ≈ sol From be9e2d2d56cd8274a83634dbca59300746aea688 Mon Sep 17 00:00:00 2001 From: Simon Carlson Date: Mon, 2 Jan 2023 20:20:26 +0100 Subject: [PATCH 2/7] Implemented Broyden and tests --- Project.toml | 1 + src/broyden.jl | 45 ++++++++++++++++++---------- test/basictests.jl | 73 +++++++++++++++++++++++++++++----------------- 3 files changed, 77 insertions(+), 42 deletions(-) diff --git a/Project.toml b/Project.toml index fa4f7f0..3f4c0d7 100644 --- a/Project.toml +++ b/Project.toml @@ -7,6 +7,7 @@ version = "0.1.4" ArrayInterfaceCore = "30b0a656-2188-435a-8636-2ec0e6a096e2" FiniteDiff = "6a86dc24-6348-571c-b903-95158fe2bd41" ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" +LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" Reexport = "189a3867-3050-52da-a836-e630ba90ab69" SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462" SnoopPrecompile = "66db9d55-30c0-4569-8b51-7e840670fc0c" diff --git a/src/broyden.jl b/src/broyden.jl index 9f9388b..7664cf6 100644 --- a/src/broyden.jl +++ b/src/broyden.jl @@ -5,13 +5,19 @@ struct Broyden <: AbstractSimpleNonlinearSolveAlgorithm end function SciMLBase.solve(prob::NonlinearProblem, - alg::Broyden, args...; abstol = nothing, - reltol = nothing, - maxiters = 1000, kwargs...) + alg::Broyden, args...; abstol = nothing, + reltol = nothing, + maxiters = 1000, kwargs...) + f = Base.Fix2(prob.f, prob.p) - xₙ = float(prob.u0) - T = typeof(xₙ) - J⁻¹ = Matrix{T}(I, length(xₙ), length(xₙ)) + x = float(prob.u0) + fₙ = f(x) + T = eltype(x) + if length(x) > 1 + J⁻¹ = Matrix{T}(I, length(x), length(x)) + else + J⁻¹ = 1 + end if SciMLBase.isinplace(prob) error("Broyden currently only supports out-of-place nonlinear problems") @@ -21,23 +27,30 @@ function SciMLBase.solve(prob::NonlinearProblem, real(oneunit(eltype(T))) * (eps(real(one(eltype(T)))))^(4 // 5) rtol = reltol !== nothing ? reltol : eps(real(one(eltype(T))))^(4 // 5) - for _ in 1:maxiters - # TODO check if nameing with heavy use of subscrips is ok + xₙ = x + xₙ₋₁ = x + fₙ₋₁ = fₙ + for n in 1:maxiters + + xₙ = xₙ₋₁ - J⁻¹ * fₙ₋₁ fₙ = f(xₙ) - xₙ₊₁ = xₙ + J⁻¹ * fₙ - Δxₙ = xₙ₊₁ - xₙ - Δfₙ = f(xₙ₊₁) - fₙ - J⁻¹ .+= ((Δxₙ .- J⁻¹ * Δfₙ) ./ (Δxₙ' * J⁻¹ * Δfₙ)) * (Δxₙ' * J⁻¹) + Δxₙ = xₙ - xₙ₋₁ + Δfₙ = fₙ - fₙ₋₁ + J⁻¹ += ((Δxₙ - J⁻¹ * Δfₙ) ./ (Δxₙ' * J⁻¹ * Δfₙ)) * (Δxₙ' * J⁻¹) iszero(fₙ) && - return SciMLBase.build_solution(prob, alg, xₙ₊₁, fₙ; retcode = ReturnCode.Success) + return SciMLBase.build_solution(prob, alg, xₙ, fₙ; + retcode = ReturnCode.Success) - if isapprox(xₙ₊₁, xₙ, atol = atol, rtol = rtol) - return SciMLBase.build_solution(prob, alg, xₙ₊₁, fₙ; retcode = ReturnCode.Success) + if isapprox(xₙ, xₙ₋₁, atol = atol, rtol = rtol) + return SciMLBase.build_solution(prob, alg, xₙ, fₙ; + retcode = ReturnCode.Success) end - xₙ = xₙ₊₁ + xₙ₋₁ = xₙ + fₙ₋₁ = fₙ end @show xₙ, fₙ + return SciMLBase.build_solution(prob, alg, xₙ, fₙ; retcode = ReturnCode.MaxIters) end diff --git a/test/basictests.jl b/test/basictests.jl index bc5a390..5518fb4 100644 --- a/test/basictests.jl +++ b/test/basictests.jl @@ -3,6 +3,7 @@ using StaticArrays using BenchmarkTools using Test +# SimpleNewtonRaphson function benchmark_scalar(f, u0) probN = NonlinearProblem{false}(f, u0) sol = (solve(probN, SimpleNewtonRaphson())) @@ -23,38 +24,53 @@ sol = benchmark_scalar(sf, csu0) @test (@ballocated benchmark_scalar(sf, csu0)) == 0 +# Broyden +function benchmark_scalar(f, u0) + probN = NonlinearProblem{false}(f, u0) + sol = (solve(probN, Broyden())) +end + +sol = benchmark_scalar(sf, csu0) +@test sol.retcode === ReturnCode.Success +@test sol.u * sol.u - 2 < 1e-9 +@test (@ballocated benchmark_scalar(sf, csu0)) == 0 + # AD Tests using ForwardDiff # Immutable f, u0 = (u, p) -> u .* u .- p, @SVector[1.0, 1.0] -g = function (p) - probN = NonlinearProblem{false}(f, csu0, p) - sol = solve(probN, SimpleNewtonRaphson(), tol = 1e-9) - return sol.u[end] -end +for alg in [SimpleNewtonRaphson(), Broyden()] + g = function (p) + probN = NonlinearProblem{false}(f, csu0, p) + sol = solve(probN, alg, tol = 1e-9) + return sol.u[end] + end -for p in 1.0:0.1:100.0 - @test g(p) ≈ sqrt(p) - @test ForwardDiff.derivative(g, p) ≈ 1 / (2 * sqrt(p)) + for p in 1.1:0.1:100.0 + @test g(p) ≈ sqrt(p) + @test ForwardDiff.derivative(g, p) ≈ 1 / (2 * sqrt(p)) + end end # Scalar f, u0 = (u, p) -> u * u - p, 1.0 +for alg in [SimpleNewtonRaphson(), Broyden()] -# SimpleNewtonRaphson -g = function (p) - probN = NonlinearProblem{false}(f, oftype(p, u0), p) - sol = solve(probN, SimpleNewtonRaphson()) - return sol.u -end - -@test ForwardDiff.derivative(g, 1.0) ≈ 0.5 + g = function (p) + probN = NonlinearProblem{false}(f, oftype(p, u0), p) + sol = solve(probN, alg) + return sol.u + end -for p in 1.1:0.1:100.0 - @test g(p) ≈ sqrt(p) + p = 1.1 @test ForwardDiff.derivative(g, p) ≈ 1 / (2 * sqrt(p)) + + for p in 1.1:0.1:100.0 + @test g(p) ≈ sqrt(p) + @test ForwardDiff.derivative(g, p) ≈ 1 / (2 * sqrt(p)) + end end tspan = (1.0, 20.0) @@ -77,7 +93,7 @@ for alg in [Bisection(), Falsi()] global g, p g = function (p) probN = IntervalNonlinearProblem{false}(f, tspan, p) - sol = solve(probN, Bisection()) # TODO check if "alg" should replace "Bisection()" + sol = solve(probN, alg) return [sol.left] end @@ -85,16 +101,18 @@ for alg in [Bisection(), Falsi()] @test ForwardDiff.jacobian(g, p) ≈ ForwardDiff.jacobian(t, p) end -gnewton = function (p) - probN = NonlinearProblem{false}(f, 0.5, p) - sol = solve(probN, SimpleNewtonRaphson()) - return [sol.u] +for alg in [SimpleNewtonRaphson(), Broyden()] + global g, p + g = function (p) + probN = NonlinearProblem{false}(f, 0.5, p) + sol = solve(probN, alg) + return [sol.u] + end + @test g(p) ≈ [sqrt(p[2] / p[1])] + @test ForwardDiff.jacobian(g, p) ≈ ForwardDiff.jacobian(t, p) end -@test gnewton(p) ≈ [sqrt(p[2] / p[1])] -@test ForwardDiff.jacobian(gnewton, p) ≈ ForwardDiff.jacobian(t, p) # Error Checks - f, u0 = (u, p) -> u .* u .- 2.0, @SVector[1.0, 1.0] probN = NonlinearProblem(f, u0) @@ -103,6 +121,8 @@ probN = NonlinearProblem(f, u0) @test solve(probN, SimpleNewtonRaphson(; autodiff = false)).u[end] ≈ sqrt(2.0) @test solve(probN, SimpleNewtonRaphson(; autodiff = false)).u[end] ≈ sqrt(2.0) # TODO check why the 2 lines above are identical +@test solve(probN, Broyden()).u[end] ≈ sqrt(2.0) +@test solve(probN, Broyden(); immutable = false).u[end] ≈ sqrt(2.0) for u0 in [1.0, [1, 1.0]] local f, probN, sol @@ -114,6 +134,7 @@ for u0 in [1.0, [1, 1.0]] @test solve(probN, SimpleNewtonRaphson()).u ≈ sol @test solve(probN, SimpleNewtonRaphson()).u ≈ sol @test solve(probN, SimpleNewtonRaphson(; autodiff = false)).u ≈ sol + @test solve(probN, Broyden()).u ≈ sol end # Bisection Tests From 7fff5402b35b3e598e7a19564791c8ef173e3104 Mon Sep 17 00:00:00 2001 From: Simon Carlson Date: Mon, 2 Jan 2023 22:24:29 +0100 Subject: [PATCH 3/7] Cleanup --- src/SimpleNonlinearSolve.jl | 2 +- src/broyden.jl | 21 ++++++++++++--------- test/basictests.jl | 8 ++------ 3 files changed, 15 insertions(+), 16 deletions(-) diff --git a/src/SimpleNonlinearSolve.jl b/src/SimpleNonlinearSolve.jl index 213ba23..0b2fd0b 100644 --- a/src/SimpleNonlinearSolve.jl +++ b/src/SimpleNonlinearSolve.jl @@ -4,7 +4,7 @@ using Reexport using FiniteDiff, ForwardDiff using ForwardDiff: Dual using StaticArraysCore -using LinearAlgebra # TODO check if it is ok to add this +using LinearAlgebra import ArrayInterfaceCore @reexport using SciMLBase diff --git a/src/broyden.jl b/src/broyden.jl index 7664cf6..35d4adb 100644 --- a/src/broyden.jl +++ b/src/broyden.jl @@ -1,14 +1,18 @@ -# TODO add docstrings +""" +```julia +Broyden() +``` + +A low-overhead implementation of Broyden. This method is non-allocating on scalar +and static array problems. +""" -# TODO check what the supertype should be -# TODO check if this should be defined as in raphson.jl struct Broyden <: AbstractSimpleNonlinearSolveAlgorithm end function SciMLBase.solve(prob::NonlinearProblem, - alg::Broyden, args...; abstol = nothing, - reltol = nothing, - maxiters = 1000, kwargs...) - + alg::Broyden, args...; abstol = nothing, + reltol = nothing, + maxiters = 1000, kwargs...) f = Base.Fix2(prob.f, prob.p) x = float(prob.u0) fₙ = f(x) @@ -30,8 +34,7 @@ function SciMLBase.solve(prob::NonlinearProblem, xₙ = x xₙ₋₁ = x fₙ₋₁ = fₙ - for n in 1:maxiters - + for _ in 1:maxiters xₙ = xₙ₋₁ - J⁻¹ * fₙ₋₁ fₙ = f(xₙ) Δxₙ = xₙ - xₙ₋₁ diff --git a/test/basictests.jl b/test/basictests.jl index 5518fb4..6a5f69b 100644 --- a/test/basictests.jl +++ b/test/basictests.jl @@ -57,16 +57,12 @@ end # Scalar f, u0 = (u, p) -> u * u - p, 1.0 for alg in [SimpleNewtonRaphson(), Broyden()] - g = function (p) probN = NonlinearProblem{false}(f, oftype(p, u0), p) sol = solve(probN, alg) return sol.u end - p = 1.1 - @test ForwardDiff.derivative(g, p) ≈ 1 / (2 * sqrt(p)) - for p in 1.1:0.1:100.0 @test g(p) ≈ sqrt(p) @test ForwardDiff.derivative(g, p) ≈ 1 / (2 * sqrt(p)) @@ -120,7 +116,7 @@ probN = NonlinearProblem(f, u0) @test solve(probN, SimpleNewtonRaphson(); immutable = false).u[end] ≈ sqrt(2.0) @test solve(probN, SimpleNewtonRaphson(; autodiff = false)).u[end] ≈ sqrt(2.0) @test solve(probN, SimpleNewtonRaphson(; autodiff = false)).u[end] ≈ sqrt(2.0) -# TODO check why the 2 lines above are identical + @test solve(probN, Broyden()).u[end] ≈ sqrt(2.0) @test solve(probN, Broyden(); immutable = false).u[end] ≈ sqrt(2.0) @@ -130,10 +126,10 @@ for u0 in [1.0, [1, 1.0]] probN = NonlinearProblem(f, u0) sol = sqrt(2) * u0 - # TODO check why the two lines below are identical @test solve(probN, SimpleNewtonRaphson()).u ≈ sol @test solve(probN, SimpleNewtonRaphson()).u ≈ sol @test solve(probN, SimpleNewtonRaphson(; autodiff = false)).u ≈ sol + @test solve(probN, Broyden()).u ≈ sol end From 297c170bbd91ab79544fbb8eb2077c3e8158f8a1 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Mon, 2 Jan 2023 17:28:03 -0500 Subject: [PATCH 4/7] Update src/broyden.jl --- src/broyden.jl | 1 - 1 file changed, 1 deletion(-) diff --git a/src/broyden.jl b/src/broyden.jl index 35d4adb..242ec67 100644 --- a/src/broyden.jl +++ b/src/broyden.jl @@ -6,7 +6,6 @@ Broyden() A low-overhead implementation of Broyden. This method is non-allocating on scalar and static array problems. """ - struct Broyden <: AbstractSimpleNonlinearSolveAlgorithm end function SciMLBase.solve(prob::NonlinearProblem, From 27e2c6befeb362e903df9c4804db921e61b7ce4a Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Mon, 2 Jan 2023 17:34:06 -0500 Subject: [PATCH 5/7] Update src/broyden.jl --- src/broyden.jl | 2 -- 1 file changed, 2 deletions(-) diff --git a/src/broyden.jl b/src/broyden.jl index 242ec67..ddeba8d 100644 --- a/src/broyden.jl +++ b/src/broyden.jl @@ -52,7 +52,5 @@ function SciMLBase.solve(prob::NonlinearProblem, fₙ₋₁ = fₙ end - @show xₙ, fₙ - return SciMLBase.build_solution(prob, alg, xₙ, fₙ; retcode = ReturnCode.MaxIters) end From 4a04a3cee5c1e09c64c3ff8899e7fd2d2ef7a14d Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Mon, 2 Jan 2023 17:43:43 -0500 Subject: [PATCH 6/7] Update src/broyden.jl --- src/broyden.jl | 6 +----- 1 file changed, 1 insertion(+), 5 deletions(-) diff --git a/src/broyden.jl b/src/broyden.jl index ddeba8d..ab4d313 100644 --- a/src/broyden.jl +++ b/src/broyden.jl @@ -16,11 +16,7 @@ function SciMLBase.solve(prob::NonlinearProblem, x = float(prob.u0) fₙ = f(x) T = eltype(x) - if length(x) > 1 - J⁻¹ = Matrix{T}(I, length(x), length(x)) - else - J⁻¹ = 1 - end + J⁻¹ = ArrayInterfaceCore.zeromatrix(x) if SciMLBase.isinplace(prob) error("Broyden currently only supports out-of-place nonlinear problems") From aca04adacdd830e24e245c43bcf573701b685677 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Mon, 2 Jan 2023 17:44:26 -0500 Subject: [PATCH 7/7] Update src/broyden.jl --- src/broyden.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/broyden.jl b/src/broyden.jl index ab4d313..49f5395 100644 --- a/src/broyden.jl +++ b/src/broyden.jl @@ -16,7 +16,7 @@ function SciMLBase.solve(prob::NonlinearProblem, x = float(prob.u0) fₙ = f(x) T = eltype(x) - J⁻¹ = ArrayInterfaceCore.zeromatrix(x) + J⁻¹ = ArrayInterfaceCore.zeromatrix(x) + I if SciMLBase.isinplace(prob) error("Broyden currently only supports out-of-place nonlinear problems")