Merge pull request #14 from CCsimon123/main

Add Broyden solver

lib/SimpleNonlinearSolve/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"

lib/SimpleNonlinearSolve/src/SimpleNonlinearSolve.jl

@@ -4,6 +4,7 @@ using Reexport
 using FiniteDiff, ForwardDiff
 using ForwardDiff: Dual
 using StaticArraysCore
+using LinearAlgebra
 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

lib/SimpleNonlinearSolve/src/broyden.jl

@@ -0,0 +1,52 @@
+"""
+```julia
+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,
+                         alg::Broyden, args...; abstol = nothing,
+                         reltol = nothing,
+                         maxiters = 1000, kwargs...)
+    f = Base.Fix2(prob.f, prob.p)
+    x = float(prob.u0)
+    fₙ = f(x)
+    T = eltype(x)
+    J⁻¹ = ArrayInterfaceCore.zeromatrix(x) + I
+
+    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)
+
+    xₙ = x
+    xₙ₋₁ = x
+    fₙ₋₁ = fₙ
+    for _ in 1:maxiters
+        xₙ = xₙ₋₁ - J⁻¹ * fₙ₋₁
+        fₙ = f(xₙ)
+        Δ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)
+
+        if isapprox(xₙ, xₙ₋₁, atol = atol, rtol = rtol)
+            return SciMLBase.build_solution(prob, alg, xₙ, fₙ;
+                                            retcode = ReturnCode.Success)
+        end
+        xₙ₋₁ = xₙ
+        fₙ₋₁ = fₙ
+    end
+
+    return SciMLBase.build_solution(prob, alg, xₙ, fₙ; retcode = ReturnCode.MaxIters)
+end

lib/SimpleNonlinearSolve/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,49 @@ 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()]
+    g = function (p)
+        probN = NonlinearProblem{false}(f, oftype(p, u0), p)
+        sol = solve(probN, alg)
+        return sol.u
+    end
 
-# 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
-
-for p in 1.1: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
 
 tspan = (1.0, 20.0)
@@ -77,24 +89,26 @@ for alg in [Bisection(), Falsi()]
     global g, p
     g = function (p)
         probN = IntervalNonlinearProblem{false}(f, tspan, p)
-        sol = solve(probN, Bisection())
+        sol = solve(probN, alg)
         return [sol.left]
     end
 
     @test g(p) ≈ [sqrt(p[2] / p[1])]
     @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 +117,9 @@ 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)
 
+@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
     f = (u, p) -> u .* u .- 2.0
@@ -112,6 +129,8 @@ 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