Implemented Broyden and tests

Author: Deltadahl

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"

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

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,24 +93,26 @@ 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
 
     @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 +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