Merge pull request #41 from avik-pal/ap/lbroyden

Add limited memory broyden implementation

Author: ChrisRackauckas

Project.toml

@@ -1,7 +1,7 @@
 name = "SimpleNonlinearSolve"
 uuid = "727e6d20-b764-4bd8-a329-72de5adea6c7"
 authors = ["SciML"]
-version = "0.1.10"
+version = "0.1.11"
 
 [deps]
 ArrayInterfaceCore = "30b0a656-2188-435a-8636-2ec0e6a096e2"

src/SimpleNonlinearSolve.jl

@@ -20,6 +20,7 @@ include("bisection.jl")
 include("falsi.jl")
 include("raphson.jl")
 include("broyden.jl")
+include("lbroyden.jl")
 include("klement.jl")
 include("trustRegion.jl")
 include("ridder.jl")
@@ -52,7 +53,7 @@ SnoopPrecompile.@precompile_all_calls begin for T in (Float32, Float64)
 end end
 
 # DiffEq styled algorithms
-export Bisection, Brent, Broyden, SimpleDFSane, Falsi, Klement, Ridder, SimpleNewtonRaphson,
-       SimpleTrustRegion
+export Bisection, Brent, Broyden, LBroyden, SimpleDFSane, Falsi, Klement,
+       Ridder, SimpleNewtonRaphson, SimpleTrustRegion
 
 end # module

src/lbroyden.jl

@@ -0,0 +1,89 @@
+"""
+    LBroyden(threshold::Int = 27)
+
+A limited memory implementation of Broyden. This method applies the L-BFGS scheme to
+Broyden's method.
+"""
+Base.@kwdef struct LBroyden <: AbstractSimpleNonlinearSolveAlgorithm
+    threshold::Int = 27
+end
+
+@views function SciMLBase.__solve(prob::NonlinearProblem,
+                                  alg::LBroyden, args...; abstol = nothing,
+                                  reltol = nothing,
+                                  maxiters = 1000, kwargs...)
+    threshold = min(maxiters, alg.threshold)
+    x = float(prob.u0)
+
+    if x isa Number
+        restore_scalar = true
+        x = [x]
+        f = u -> prob.f(u[], prob.p)
+    else
+        f = Base.Fix2(prob.f, prob.p)
+        restore_scalar = false
+    end
+
+    fₙ = f(x)
+    T = eltype(x)
+
+    if SciMLBase.isinplace(prob)
+        error("LBroyden currently only supports out-of-place nonlinear problems")
+    end
+
+    U = fill!(similar(x, (threshold, length(x))), zero(T))
+    Vᵀ = fill!(similar(x, (length(x), threshold)), zero(T))
+
+    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ₙ
+    update = fₙ
+    for i in 1:maxiters
+        xₙ = xₙ₋₁ .+ update
+        fₙ = f(xₙ)
+        Δxₙ = xₙ .- xₙ₋₁
+        Δfₙ = fₙ .- fₙ₋₁
+
+        if iszero(fₙ)
+            xₙ = restore_scalar ? xₙ[] : xₙ
+            return SciMLBase.build_solution(prob, alg, xₙ, fₙ; retcode = ReturnCode.Success)
+        end
+
+        if isapprox(xₙ, xₙ₋₁; atol, rtol)
+            xₙ = restore_scalar ? xₙ[] : xₙ
+            return SciMLBase.build_solution(prob, alg, xₙ, fₙ; retcode = ReturnCode.Success)
+        end
+
+        _U = U[1:min(threshold, i), :]
+        _Vᵀ = Vᵀ[:, 1:min(threshold, i)]
+
+        vᵀ = _rmatvec(_U, _Vᵀ, Δxₙ)
+        mvec = _matvec(_U, _Vᵀ, Δfₙ)
+        Δxₙ = (Δxₙ .- mvec) ./ (sum(vᵀ .* Δfₙ) .+ convert(T, 1e-5))
+
+        Vᵀ[:, mod1(i, threshold)] .= vᵀ
+        U[mod1(i, threshold), :] .= Δxₙ
+
+        update = -_matvec(U[1:min(threshold, i + 1), :], Vᵀ[:, 1:min(threshold, i + 1)], fₙ)
+
+        xₙ₋₁ = xₙ
+        fₙ₋₁ = fₙ
+    end
+
+    xₙ = restore_scalar ? xₙ[] : xₙ
+    return SciMLBase.build_solution(prob, alg, xₙ, fₙ; retcode = ReturnCode.MaxIters)
+end
+
+function _rmatvec(U::AbstractMatrix, Vᵀ::AbstractMatrix,
+                  x::Union{<:AbstractVector, <:Number})
+    return -x .+ dropdims(sum(U .* sum(Vᵀ .* x; dims = 1)'; dims = 1); dims = 1)
+end
+
+function _matvec(U::AbstractMatrix, Vᵀ::AbstractMatrix,
+                 x::Union{<:AbstractVector, <:Number})
+    return -x .+ dropdims(sum(sum(x .* U'; dims = 1) .* Vᵀ; dims = 2); dims = 2)
+end

test/basictests.jl

@@ -78,7 +78,7 @@ using ForwardDiff
 # Immutable
 f, u0 = (u, p) -> u .* u .- p, @SVector[1.0, 1.0]
 
-for alg in (SimpleNewtonRaphson(), Broyden(), Klement(), SimpleTrustRegion(),
+for alg in (SimpleNewtonRaphson(), Broyden(), LBroyden(), Klement(), SimpleTrustRegion(),
             SimpleDFSane())
     g = function (p)
         probN = NonlinearProblem{false}(f, csu0, p)
@@ -94,7 +94,7 @@ end
 
 # Scalar
 f, u0 = (u, p) -> u * u - p, 1.0
-for alg in (SimpleNewtonRaphson(), Broyden(), Klement(), SimpleTrustRegion(),
+for alg in (SimpleNewtonRaphson(), Broyden(), LBroyden(), Klement(), SimpleTrustRegion(),
             SimpleDFSane())
     g = function (p)
         probN = NonlinearProblem{false}(f, oftype(p, u0), p)
@@ -160,13 +160,13 @@ for alg in [Bisection(), Falsi(), Ridder(), Brent()]
     @test ForwardDiff.jacobian(g, p) ≈ ForwardDiff.jacobian(t, p)
 end
 
-for alg in (SimpleNewtonRaphson(), Broyden(), Klement(), SimpleTrustRegion(),
+for alg in (SimpleNewtonRaphson(), Broyden(), LBroyden(), Klement(), SimpleTrustRegion(),
             SimpleDFSane())
     global g, p
     g = function (p)
         probN = NonlinearProblem{false}(f, 0.5, p)
         sol = solve(probN, alg)
-        return [sol.u]
+        return [abs(sol.u)]
     end
     @test g(p) ≈ [sqrt(p[2] / p[1])]
     @test ForwardDiff.jacobian(g, p) ≈ ForwardDiff.jacobian(t, p)
@@ -181,6 +181,7 @@ probN = NonlinearProblem(f, u0)
 @test solve(probN, SimpleTrustRegion()).u[end] ≈ sqrt(2.0)
 @test solve(probN, SimpleTrustRegion(; autodiff = false)).u[end] ≈ sqrt(2.0)
 @test solve(probN, Broyden()).u[end] ≈ sqrt(2.0)
+@test solve(probN, LBroyden()).u[end] ≈ sqrt(2.0)
 @test solve(probN, Klement()).u[end] ≈ sqrt(2.0)
 @test solve(probN, SimpleDFSane()).u[end] ≈ sqrt(2.0)
 
@@ -199,6 +200,7 @@ for u0 in [1.0, [1, 1.0]]
     @test solve(probN, SimpleTrustRegion(; autodiff = false)).u ≈ sol
 
     @test solve(probN, Broyden()).u ≈ sol
+    @test solve(probN, LBroyden()).u ≈ sol
     @test solve(probN, Klement()).u ≈ sol
     @test solve(probN, SimpleDFSane()).u ≈ sol
 end